Nothing
inst/doc/BasicsIdentifiabilityPower.R
In mpt: Multinomial Processing Tree Models
## -----------------------------------------------------------------------------
library(mpt)
## -----------------------------------------------------------------------------
dat <- read.table(header = TRUE, text = "
lag age resp freq
0 young E1 90 # adjacent recall: apple, banana
0 young E2 14 # non-adjacent recall: apple, ..., banana
0 young E3 84 # single recall: only apple or only banana
0 young E4 212 # non-recall: neither apple nor banana
n/a young F1 102 # recall
n/a young F2 298 # non-recall
0 old E1 42
0 old E2 5
0 old E3 63
0 old E4 290
n/a old F1 64
n/a old F2 336
15 young E1 67
15 young E2 18
15 young E3 123
15 young E4 192
15 old E1 30
15 old E2 13
15 old E3 95
15 old E4 262
")
dat$age <- factor(dat$age, levels = c("young", "old"))
## -----------------------------------------------------------------------------
# Pairs Singletons
#
# /r -- E1
# c
# / \1-r -- E4 a -- F1 c: cluster storage
# / / r: cluster retrieval
# \ /u -- E2 \ u: non-clustered storage/retrieval
# \ u 1-a -- F2 a: singleton storage/retrieval
# \ / \1-u -- E3
# 1-c
# \ /u -- E3
# 1-u
# \1-u -- E4
mptspec( # specify model
E.1 = c*r,
E.2 = (1 - c)*u^2, # dot notation:
E.3 = 2 * (1 - c)*u*(1 - u), # tree.response
E.4 = c*(1 - r) + (1 - c)*(1 - u)^2,
F.1 = a,
F.2 = 1 - a,
.restr = list(a = u)
) |>
mpt(data = dat[dat$lag != 15 & dat$age == "young", ]) # fit to data
## -----------------------------------------------------------------------------
m <- mpt(mptspec("SR"),
data = dat[dat$lag != 15 & dat$age == "young", ])
## -----------------------------------------------------------------------------
print(m)
summary(m)
logLik(m)
AIC(m)
BIC(m)
coef(m)
confint(m, logit = FALSE)
## -----------------------------------------------------------------------------
# Target Distractor
#
# r -- Hit
# / } 55 b -- FA 45
# / b -- Hit / r: recognition
# \ / \ b: bias
# 1-r 1-b -- CR 765
# \
# 1-b -- Miss 35
## -----------------------------------------------------------------------------
nll <- function(theta, data, treeid) { # negative log-likelihood
c <- theta[1]
r <- theta[2]
u <- theta[3]
a <- theta[3] # a = u
pcat <- c(
E.1 = c*r,
E.2 = (1 - c)*u^2,
E.3 = 2 * (1 - c)*u*(1 - u),
E.4 = c*(1 - r) + (1 - c)*(1 - u)^2,
F.1 = a,
F.2 = 1 - a
)
-sum(
dmultinom(data[treeid == 1], size = sum(data[treeid == 1]),
prob = pcat[treeid == 1], log = TRUE),
dmultinom(data[treeid == 2], size = sum(data[treeid == 2]),
prob = pcat[treeid == 2], log = TRUE)
)
}
## -----------------------------------------------------------------------------
optim(par = c(.5, .5, .5), # starting values
fn = nll, # function to be minimized
data = dat[dat$lag != 15 &
dat$age == "young", "freq"], # arguments passed to fn
treeid = c(1, 1, 1, 1, 2, 2),
method = "L-BFGS-B", # algorithm
lower = rep(.001, 3), # boundaries
upper = rep(.999, 3),
hessian = TRUE)
## -----------------------------------------------------------------------------
m0 <-
update(m$spec, .restr = list(c = 0.5)) |> # specify restriction(s)
mpt(data = m$y) |> # refit
print()
## -----------------------------------------------------------------------------
anova(m0, m)
## -----------------------------------------------------------------------------
m1 <- mptspec("SR", .replicates = 2) |>
mpt(data = dat[dat$lag != 15, ]) |>
print()
## -----------------------------------------------------------------------------
update(m1$spec, .restr = list(c2 = s * c1)) |>
mpt(data = m1$y)
## -----------------------------------------------------------------------------
m5 <-
update(m1$spec, .restr = list(c2 = s_c * c1, # two shrinkage pars
r2 = s_r * r1)) |>
mpt(data = m1$y)
m6 <-
update(m5$spec, .restr = list(s_c = s_r)) |> # single shrinkage par
mpt(data = m1$y)
anova(m6, m5)
## -----------------------------------------------------------------------------
s1 <- mptspec( # not identifiable: 5 parameters
E.1 = c * r,
E.2 = (1 - c) * u1 * u2,
E.3 = (1 - c) * u1 * (1 - u2) + (1 - c) * u2 * (1 - u1),
E.4 = c * (1 - r) + (1 - c) * (1 - u1) * (1 - u2),
F.1 = a,
F.2 = 1 - a
)
## -----------------------------------------------------------------------------
runif(5) |> s1$par2deriv() |> _$deriv |> print() |> qr() |> _$rank
## -----------------------------------------------------------------------------
s2 <- update(s1, .restr = list(u1 = u, u2 = u)) # make it identifiable
runif(4) |> s2$par2deriv() |> _$deriv |> print() |> qr() |> _$rank
## -----------------------------------------------------------------------------
suppressWarnings(replicate(10,
mpt(s1, data = dat[dat$lag != 15 & dat$age == "young", ],
start = runif(5)) |>
coef()
)) |> t() # not identifiable: parameter trade-offs
replicate(10,
mpt(s2, data = dat[dat$lag != 15 & dat$age == "young", ],
start = runif(4)) |>
coef()
) |> t() # identifiable: unique estimates
## -----------------------------------------------------------------------------
s <- mptspec(
E.1 = c * r,
E.2 = (1 - c) * u^2,
E.3 = 2 * (1 - c) * u * (1 - u),
E.4 = c * (1 - r) + (1 - c) * (1 - u)^2,
F.1 = a,
F.2 = 1 - a
)
s$par # before you use par2prob(), carefully check position of parameters!
s$par2prob(c(c = 0.5, r = 0.5, u = 0.4, a = 0.6)) # evaluate model equations
dataGen <- function(nn, d) {
structure(list( # stub mpt object
treeid = s$treeid,
n = setNames((nn * c(2, 1)/3)[s$treeid], s$treeid), # 2:1 ratio
pcat = s$par2prob(c(c = 0.5, r = 0.5,
u = 0.5 - d/2, a = 0.5 + d/2))
), class = "mpt") |>
simulate()
}
dataGen(960, 0.2)
## ----fig.width = 5------------------------------------------------------------
replicate(20, dataGen(960, 0.2) |> mpt(s, data = _) |> coef()) |>
t() |>
boxplot(horizontal = TRUE)
## -----------------------------------------------------------------------------
testFun <- function(nn, d) {
y <- dataGen(nn, d) # generate data with effect
m1 <- mpt(s, y)
m2 <- mpt(update(s, .restr = list(a = u)), y)
anova(m2, m1)$"Pr(>Chi)"[2] # test H0, return p value
}
## -----------------------------------------------------------------------------
# pwrsim1 <- expand.grid(d = seq(0, 0.5, 0.1), n = 30 * 2^(0:5))
#
# system.time( # about 20 min
# pwrsim1$pval <-
# mapply(function(nn, d) replicate(500, testFun(nn, d)),
# nn = pwrsim1$n, d = pwrsim1$d, SIMPLIFY = FALSE)
# )
# pwrsim1$pwr <- sapply(pwrsim1$pval, function(p) mean(p < .05))
# }
## ----fig.width = 5------------------------------------------------------------
data(pwrsim) # provides pwrsim1 and pwrsim2
if(requireNamespace("lattice")) {
lattice::xyplot(pwr ~ d, pwrsim1, groups = n, type = c("g", "b"),
auto.key = list(space = "right"),
xlab = "Parameter difference a - u",
ylab = "Simulated power")
}
## ----fig.width = 5------------------------------------------------------------
s <- mptspec("SR", .replicates = 2)
dataGen <- function(nn, d) {
structure(list(
treeid = s$treeid,
n = setNames((nn/2 * c(2, 1, 2, 1)/3)[s$treeid], s$treeid),
pcat = s$par2prob(c(c1 = 0.5, r1 = 0.4 + d/2, u1 = 0.3, # young
c2 = 0.5, r2 = 0.4 - d/2, u2 = 0.3)) # old
), class = "mpt") |>
simulate(m)
}
testFun <- function(nn, d) {
y <- dataGen(nn, d)
m1 <- mpt(s, y)
m2 <- mpt(update(s, .restr = list(r1 = r2)), y)
anova(m2, m1)$"Pr(>Chi)"[2]
}
# pwrsim2 <- expand.grid(d = seq(0, 0.4, 0.1), n = 120 * 2^(0:2))
if(requireNamespace("lattice")) {
lattice::xyplot(pwr ~ d, pwrsim2, groups = n, type = c("g", "b"),
auto.key = list(space = "right"),
xlab = "Parameter difference r1 - r2",
ylab = "Simulated power")
}
Try the mpt package in your browser
Any scripts or data that you put into this service are public.
mpt documentation built on Oct. 11, 2024, 5:07 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## -----------------------------------------------------------------------------
library(mpt)
## -----------------------------------------------------------------------------
dat <- read.table(header = TRUE, text = "
lag age resp freq
0 young E1 90 # adjacent recall: apple, banana
0 young E2 14 # non-adjacent recall: apple, ..., banana
0 young E3 84 # single recall: only apple or only banana
0 young E4 212 # non-recall: neither apple nor banana
n/a young F1 102 # recall
n/a young F2 298 # non-recall
0 old E1 42
0 old E2 5
0 old E3 63
0 old E4 290
n/a old F1 64
n/a old F2 336
15 young E1 67
15 young E2 18
15 young E3 123
15 young E4 192
15 old E1 30
15 old E2 13
15 old E3 95
15 old E4 262
")
dat$age <- factor(dat$age, levels = c("young", "old"))
## -----------------------------------------------------------------------------
# Pairs Singletons
#
# /r -- E1
# c
# / \1-r -- E4 a -- F1 c: cluster storage
# / / r: cluster retrieval
# \ /u -- E2 \ u: non-clustered storage/retrieval
# \ u 1-a -- F2 a: singleton storage/retrieval
# \ / \1-u -- E3
# 1-c
# \ /u -- E3
# 1-u
# \1-u -- E4
mptspec( # specify model
E.1 = c*r,
E.2 = (1 - c)*u^2, # dot notation:
E.3 = 2 * (1 - c)*u*(1 - u), # tree.response
E.4 = c*(1 - r) + (1 - c)*(1 - u)^2,
F.1 = a,
F.2 = 1 - a,
.restr = list(a = u)
) |>
mpt(data = dat[dat$lag != 15 & dat$age == "young", ]) # fit to data
## -----------------------------------------------------------------------------
m <- mpt(mptspec("SR"),
data = dat[dat$lag != 15 & dat$age == "young", ])
## -----------------------------------------------------------------------------
print(m)
summary(m)
logLik(m)
AIC(m)
BIC(m)
coef(m)
confint(m, logit = FALSE)
## -----------------------------------------------------------------------------
# Target Distractor
#
# r -- Hit
# / } 55 b -- FA 45
# / b -- Hit / r: recognition
# \ / \ b: bias
# 1-r 1-b -- CR 765
# \
# 1-b -- Miss 35
## -----------------------------------------------------------------------------
nll <- function(theta, data, treeid) { # negative log-likelihood
c <- theta[1]
r <- theta[2]
u <- theta[3]
a <- theta[3] # a = u
pcat <- c(
E.1 = c*r,
E.2 = (1 - c)*u^2,
E.3 = 2 * (1 - c)*u*(1 - u),
E.4 = c*(1 - r) + (1 - c)*(1 - u)^2,
F.1 = a,
F.2 = 1 - a
)
-sum(
dmultinom(data[treeid == 1], size = sum(data[treeid == 1]),
prob = pcat[treeid == 1], log = TRUE),
dmultinom(data[treeid == 2], size = sum(data[treeid == 2]),
prob = pcat[treeid == 2], log = TRUE)
)
}
## -----------------------------------------------------------------------------
optim(par = c(.5, .5, .5), # starting values
fn = nll, # function to be minimized
data = dat[dat$lag != 15 &
dat$age == "young", "freq"], # arguments passed to fn
treeid = c(1, 1, 1, 1, 2, 2),
method = "L-BFGS-B", # algorithm
lower = rep(.001, 3), # boundaries
upper = rep(.999, 3),
hessian = TRUE)
## -----------------------------------------------------------------------------
m0 <-
update(m$spec, .restr = list(c = 0.5)) |> # specify restriction(s)
mpt(data = m$y) |> # refit
print()
## -----------------------------------------------------------------------------
anova(m0, m)
## -----------------------------------------------------------------------------
m1 <- mptspec("SR", .replicates = 2) |>
mpt(data = dat[dat$lag != 15, ]) |>
print()
## -----------------------------------------------------------------------------
update(m1$spec, .restr = list(c2 = s * c1)) |>
mpt(data = m1$y)
## -----------------------------------------------------------------------------
m5 <-
update(m1$spec, .restr = list(c2 = s_c * c1, # two shrinkage pars
r2 = s_r * r1)) |>
mpt(data = m1$y)
m6 <-
update(m5$spec, .restr = list(s_c = s_r)) |> # single shrinkage par
mpt(data = m1$y)
anova(m6, m5)
## -----------------------------------------------------------------------------
s1 <- mptspec( # not identifiable: 5 parameters
E.1 = c * r,
E.2 = (1 - c) * u1 * u2,
E.3 = (1 - c) * u1 * (1 - u2) + (1 - c) * u2 * (1 - u1),
E.4 = c * (1 - r) + (1 - c) * (1 - u1) * (1 - u2),
F.1 = a,
F.2 = 1 - a
)
## -----------------------------------------------------------------------------
runif(5) |> s1$par2deriv() |> _$deriv |> print() |> qr() |> _$rank
## -----------------------------------------------------------------------------
s2 <- update(s1, .restr = list(u1 = u, u2 = u)) # make it identifiable
runif(4) |> s2$par2deriv() |> _$deriv |> print() |> qr() |> _$rank
## -----------------------------------------------------------------------------
suppressWarnings(replicate(10,
mpt(s1, data = dat[dat$lag != 15 & dat$age == "young", ],
start = runif(5)) |>
coef()
)) |> t() # not identifiable: parameter trade-offs
replicate(10,
mpt(s2, data = dat[dat$lag != 15 & dat$age == "young", ],
start = runif(4)) |>
coef()
) |> t() # identifiable: unique estimates
## -----------------------------------------------------------------------------
s <- mptspec(
E.1 = c * r,
E.2 = (1 - c) * u^2,
E.3 = 2 * (1 - c) * u * (1 - u),
E.4 = c * (1 - r) + (1 - c) * (1 - u)^2,
F.1 = a,
F.2 = 1 - a
)
s$par # before you use par2prob(), carefully check position of parameters!
s$par2prob(c(c = 0.5, r = 0.5, u = 0.4, a = 0.6)) # evaluate model equations
dataGen <- function(nn, d) {
structure(list( # stub mpt object
treeid = s$treeid,
n = setNames((nn * c(2, 1)/3)[s$treeid], s$treeid), # 2:1 ratio
pcat = s$par2prob(c(c = 0.5, r = 0.5,
u = 0.5 - d/2, a = 0.5 + d/2))
), class = "mpt") |>
simulate()
}
dataGen(960, 0.2)
## ----fig.width = 5------------------------------------------------------------
replicate(20, dataGen(960, 0.2) |> mpt(s, data = _) |> coef()) |>
t() |>
boxplot(horizontal = TRUE)
## -----------------------------------------------------------------------------
testFun <- function(nn, d) {
y <- dataGen(nn, d) # generate data with effect
m1 <- mpt(s, y)
m2 <- mpt(update(s, .restr = list(a = u)), y)
anova(m2, m1)$"Pr(>Chi)"[2] # test H0, return p value
}
## -----------------------------------------------------------------------------
# pwrsim1 <- expand.grid(d = seq(0, 0.5, 0.1), n = 30 * 2^(0:5))
#
# system.time( # about 20 min
# pwrsim1$pval <-
# mapply(function(nn, d) replicate(500, testFun(nn, d)),
# nn = pwrsim1$n, d = pwrsim1$d, SIMPLIFY = FALSE)
# )
# pwrsim1$pwr <- sapply(pwrsim1$pval, function(p) mean(p < .05))
# }
## ----fig.width = 5------------------------------------------------------------
data(pwrsim) # provides pwrsim1 and pwrsim2
if(requireNamespace("lattice")) {
lattice::xyplot(pwr ~ d, pwrsim1, groups = n, type = c("g", "b"),
auto.key = list(space = "right"),
xlab = "Parameter difference a - u",
ylab = "Simulated power")
}
## ----fig.width = 5------------------------------------------------------------
s <- mptspec("SR", .replicates = 2)
dataGen <- function(nn, d) {
structure(list(
treeid = s$treeid,
n = setNames((nn/2 * c(2, 1, 2, 1)/3)[s$treeid], s$treeid),
pcat = s$par2prob(c(c1 = 0.5, r1 = 0.4 + d/2, u1 = 0.3, # young
c2 = 0.5, r2 = 0.4 - d/2, u2 = 0.3)) # old
), class = "mpt") |>
simulate(m)
}
testFun <- function(nn, d) {
y <- dataGen(nn, d)
m1 <- mpt(s, y)
m2 <- mpt(update(s, .restr = list(r1 = r2)), y)
anova(m2, m1)$"Pr(>Chi)"[2]
}
# pwrsim2 <- expand.grid(d = seq(0, 0.4, 0.1), n = 120 * 2^(0:2))
if(requireNamespace("lattice")) {
lattice::xyplot(pwr ~ d, pwrsim2, groups = n, type = c("g", "b"),
auto.key = list(space = "right"),
xlab = "Parameter difference r1 - r2",
ylab = "Simulated power")
}
Try the mpt package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.