birtRcpp fits joint
models of response accuracy and response time. The main estimator is
ecm(): maximum marginal likelihood by the ECM algorithm
(Meng & Rubin, 1993), with standard errors, AIC and BIC, usually in
about a second. gibbs() samples the posterior of the same
model, starting from the ECM solution. Both are written in C++. The
package ports the Julia package ExtendedRtIrtModeling.jl, with corrected
full conditionals. Notation: j indexes persons, i items.
The model is named by model = in ecm() and
gibbs():
model |
model |
|---|---|
"irt" |
2PL (or 1PL, itemtype = "1pl") with a latent regression
of ability |
"null" |
2PL + log-normal response times (van der Linden, 2007) |
"latreg" |
… with latent regressions of ability and speed |
"latent" |
speed regressed on ability and covariates; quantile =
(Gibbs) for quantile regression |
"cross" |
cross-relations of ability with each log response time;
quantile = (Gibbs) |
input_data() takes responses, response times (in
seconds) and centred covariates, either as matrices or as column names
of a data frame:
library(birtRcpp)
#>
#> Attaching package: 'birtRcpp'
#> The following object is masked from 'package:utils':
#>
#> vi
set.seed(1)
cond <- set_cond(n_subj = 300, n_item = 8, n_feat = 2)
sim <- sim_data(cond, sim_para(cond, "latent"), "latent")
raw <- data.frame(sim$Y, sim$T, sim$X)
names(raw) <- c(paste0("y", 1:8), paste0("t", 1:8), "x1", "x2")
dat <- input_data(resp = paste0("y", 1:8), time = paste0("t", 1:8), cov = c("x1", "x2"), data = raw)
dat
#> input_data: 300 persons x 8 items, response times, 2 covariates
#> items: y1, y2, y3, y4, y5, y6 ...
#> covariates: x1, x2Data from an experiment program often come in long format (one row
per person and item), with times in milliseconds, text responses and
categorical covariates. input_data() refuses long data and
text responses with the code to fix them; the steps are:
wide <- reshape(long[c("pid", "item", "correct", "rt_ms")], idvar = "pid", timevar = "item",
direction = "wide") # one row per person
wide <- merge(wide, unique(long[c("pid", "group")]), by = "pid") # person-level columns back
items <- paste0("correct.q", 1:8); times <- paste0("rt_ms.q", 1:8) # names made by reshape()
wide[times] <- wide[times] / 1000 # milliseconds -> seconds (only lambda changes)
wide <- wide[stats::complete.cases(wide[c(items, times)]), ] # complete data only
groups <- stats::model.matrix(~ group, wide)[, -1] # k - 1 dummies for k groups
wide[colnames(groups)] <- scale(groups, scale = FALSE) # centred
dat <- input_data(resp = items, time = times, cov = colnames(groups), id = "pid", data = wide)Section 11 of the intro tutorial runs this from a CSV file to a table of the regression coefficients and person scores merged back into the data.
ECM integrates ability over Gauss-Hermite nodes and speed in closed form. Each M step is split into conditional maximizations, one per block of parameters, and each has a closed form: the logistic item block through the Polya-Gamma bound (Polson, Scott & Windle, 2013), the response-time and structural blocks as normal models. The log-likelihood never decreases; SQUAREM accelerates it, and the standard errors come from the observed information.
fit <- ecm(dat, model = "latent")
fit
#> model = "latent" (rtirt_latent)
#> 300 persons x 8 items; covariates: x1, x2
#> maximum marginal likelihood (ECM, 41 Gauss-Hermite nodes): 10 iterations, 0.2 s
#> logLik -3844.078, npar 36, AIC 7760.155, BIC 7893.491
#> Reliability: ability 0.679, speed 0.959
#> Use summary(), coef(), scores(), estimates(), fit_indices(); gibbs(data, model = ..., init = fit) to sample from here.
estimates(fit, pars = "beta|b_ability") # est, se, z, Wald 95% interval
#> # A tibble: 3 x 7
#> parameter est se z q025 q975 sig
#> beta[x1] -0.031 0.024 -1.314 -0.078 0.015
#> beta[x2] -0.735 0.022 -33.409 -0.778 -0.692 *
#> b_ability -0.223 0.029 -7.603 -0.281 -0.166 *
head(scores(fit)) # EAP and posterior SD of ability and speed
#> # A tibble: 6 x 5
#> id ability ability_psd speed speed_psd
#> 1 -0.164 0.542 -0.219 0.174
#> 2 -0.086 0.541 1.140 0.174
#> 3 -0.197 0.542 -1.328 0.174
#> 4 1.186 0.598 -0.236 0.175
#> 5 -0.699 0.558 -1.378 0.174
#> 6 -0.812 0.564 -1.064 0.174
fit_indices(fit) # logLik, npar, N, AIC, BIC (also logLik(fit), AIC(fit), BIC(fit))
#> # A tibble: 1 x 6
#> method logLik npar N AIC BIC
#> ECM -3844.078 36 300 7760.155 7893.491Tables print with 3 decimals (counts without decimals, p-values below
0.001 as <0.001); the values keep full precision:
estimates(fit)$est,
as.data.frame(estimates(fit)).
anova() compares fits of the same data: AIC and BIC, and
for nested models a likelihood-ratio test; model = "null"
is model = "latreg" without covariates (it leaves them
out):
f0 <- ecm(dat, model = "null")
f1 <- ecm(dat, model = "latreg")
anova(f0, f1)
#> Likelihood ratio tests of ecm() fits (each row against the row above; valid for nested fits)
#>
#> f0: model = "null" (rtirt_null)
#> f1: model = "latreg" (rtirt_latreg)
#>
#> npar logLik AIC BIC Chisq Df Pr(>Chisq)
#> f0 34 -4077.213 8222.425 8348.354
#> f1 38 -3842.950 7761.899 7902.643 468.526 4 <0.001 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Other settings: se, and prior = "map" with
priors = rtirt_priors(...) for the posterior mode in the
item parameters (small samples); control = ecm_control()
sets nodes (Gauss-Hermite nodes, default 41),
max_iter, tol and accelerate.
init = fit starts every chain near the ECM solution:
quiet <- gibbs_control(verbose = FALSE)
g <- gibbs(dat, model = "latent", init = fit, n_iter = 1500, n_chain = 2, seed = 1, control = quiet)
estimates(g, pars = "beta|b_ability") # posterior mean, SD, 95% interval, R-hat, ESS
#> # A tibble: 3 x 8
#> parameter est sd q025 q975 rhat ess sig
#> beta[x1] -0.031 0.025 -0.081 0.016 1.006 719
#> beta[x2] -0.734 0.023 -0.779 -0.692 1.004 621 *
#> b_ability -0.224 0.029 -0.283 -0.168 1.001 292 *
convergence(g)
#> 61.1% of 36 parameters with ESS > 400 and R-hat < 1.05
#> parameter rhat ess
#> b[y5] 1.017 398
#> a[y3] 1.015 260
#> lambda[y8] 1.009 738
#> b[y1] 1.007 389
#> beta[x1] 1.006 719Quantile models use the asymmetric Laplace (ALD) working likelihood:
the estimates describe a quantile of the log times
("cross") or of speed ("latent"), not a model
of the whole distribution. ecm() maximizes the working
likelihood and gives sandwich standard errors:
fit_q <- ecm(dat, model = "latent", quantile = 0.25)
estimates(fit_q, pars = "beta|b_ability")
#> # A tibble: 3 x 7
#> parameter est se z q025 q975 sig
#> beta[x1] -0.031 0.024 -1.300 -0.078 0.016
#> beta[x2] -0.736 0.023 -32.053 -0.781 -0.691 *
#> b_ability -0.250 0.033 -7.523 -0.315 -0.185 *gibbs() samples the same model. For several quantile
levels, give a vector; the result holds one model per level:
qs <- gibbs(dat,
model = "latent", quantile = c(0.25, 0.5, 0.75), n_iter = 1000, n_chain = 2,
seed = 1, control = quiet
)
qs
#> model = "latent" (rtirt_latent) at quantiles 0.25, 0.50, 0.75
#> parameter q=0.25 q=0.50 q=0.75
#> beta[x1] -0.028 -0.031 -0.025
#> beta[x2] -0.737 -0.722 -0.727
#> b_ability -0.253 -0.236 -0.273
#> Warning (q = 0.25): only 27.8% of 36 parameters have ESS > 400 and R-hat < 1.05 (max R-hat 1.164, min ESS 30); run longer (n_iter) or see convergence()
#> Warning (q = 0.5): only 22.2% of 36 parameters have ESS > 400 and R-hat < 1.05 (max R-hat 1.060, min ESS 51); run longer (n_iter) or see convergence()
#> Warning (q = 0.75): only 30.6% of 36 parameters have ESS > 400 and R-hat < 1.05 (max R-hat 1.323, min ESS 23); run longer (n_iter) or see convergence()Use longer chains in practice (the defaults are
n_iter = 5000, n_chain = 4) and check
convergence(). The intro
tutorial (in Traditional Chinese; read it in the browser) explains
every line and output; its source is
system.file("tutorial", "birtRcpp_intro.qmd", package = "birtRcpp").
Meng, X.-L., & Rubin, D. B. (1993). Maximum likelihood estimation via the ECM algorithm: A general framework. Biometrika, 80, 267-278.
Polson, N. G., Scott, J. G., & Windle, J. (2013). Bayesian inference for logistic models using Polya-Gamma latent variables. Journal of the American Statistical Association, 108, 1339-1349.
van der Linden, W. J. (2007). A hierarchical framework for modeling speed and accuracy on test items. Psychometrika, 72, 287-308.