BIvariate Response Time Modeling with 'Rcpp'
birtRcpp fits joint models of item responses and response times in R, from their full conditionals.
Name a built-in model (2PL with log-normal response times, latent regressions of ability and
speed, cross-relations of ability with each response time, and their quantile versions), as in
mirt: ecm(dat, model = "cross"). The main estimator is ecm(): maximum marginal likelihood by
ECM, with standard errors, AIC and BIC, usually in about a second. The same conditional steps drive
the Gibbs sampler gibbs(), and every draw runs in C++ (RcppArmadillo). The built-in models port the Julia package ExtendedRtIrtModeling.jl with corrected
full conditionals; the Julia function names are kept as aliases.
中文版（完整說明）：README.zh-TW.md
Start here. The birtRcpp intro tutorial (Traditional Chinese; read it in the browser, no Quarto needed) explains every line and every output; the vignette is the short English version.
| You are | Read |
|---|---|
| new to R or to IRT | section 0 of the intro tutorial: the data, IRT and response times in plain words, installing, R basics, a first analysis in ten lines |
| a student or new to these models | the intro tutorial, then the Quick start below |
| a researcher who wants the methods | the Highlights below and ?ecm, ?gibbs, ?algorithm
|
| coming from ExtendedRtIrtModeling.jl | the Julia names (?julia_aliases) and README.zh-TW.md
|
| coming from older birtRcpp code | the deprecation messages and NEWS.md
|
# install.packages("remotes")
remotes::install_github("jiewenTsai/birtRcpp", build_vignettes = TRUE) # needs a C++ compiler
install.packages("birtRcpp") # from CRAN, once released
A C++ compiler means Rtools on Windows and the Xcode command line tools on macOS.
build_vignettes = TRUE installs vignette("birtRcpp"); without it install_github() skips
the vignette and vignette("birtRcpp") says "not found" (read it on GitHub instead).
library(birtRcpp) reports that vi masks utils::vi: that is the deprecated vi()
(variational Bayes) of birtRcpp, and utils::vi() (an editor) is still there under its full
name.
Three steps: input_data(), ecm() with the model named, and the tables.
library(birtRcpp)
set.seed(1)
cond <- set_cond(n_subj = 500, n_item = 10)
sim <- sim_data(cond, sim_para(cond, "cross"), "cross") # simulated responses and seconds
raw <- data.frame(id = sprintf("P%03d", 1:500), sim$Y, sim$T) # a data frame as it usually comes
names(raw) <- c("id", paste0("q", 1:10), paste0("q", 1:10, "_sec"))
dat <- input_data(resp = paste0("q", 1:10), time = paste0("q", 1:10, "_sec"), id = "id", data = raw)
fit <- ecm(dat, model = "cross") # logit P(y = 1) = a(θ − b); log t = λ − ζ − ρθ + ε
fit
model = "cross" (rtirt_cross)
500 persons x 10 items
maximum marginal likelihood (ECM, 41 Gauss-Hermite nodes): 14 iterations, 0.2 s
logLik -8309.032, npar 51, AIC 16720.064, BIC 16935.009
Reliability: ability 0.746, speed 0.975
Use summary(), coef(), scores(), estimates(), fit_indices(); gibbs(data, model = ..., init = fit) to sample from here.
summary(fit)$items # one row per item (summary(fit)$se has the standard errors)
item parameters (ECM estimates)
item a b lambda sigma2t rho
q1 1.013 -1.113 2.415 0.269 -0.091
q2 0.911 -0.711 3.608 0.305 0.077
q3 0.669 -0.533 3.061 0.240 0.163
q4 0.829 -0.192 2.818 0.293 0.026
q5 1.247 -0.215 3.851 0.299 0.338
...
head(scores(fit), 3) # EAP and posterior SD of ability and speed, with the id
# A tibble: 3 x 5
id ability ability_psd speed speed_psd
P001 -0.126 0.487 0.473 0.171
P002 -0.770 0.494 0.140 0.171
P003 -0.404 0.487 -1.109 0.171
Tables print with 3 decimals (counts without decimals, p-values below 0.001 as <0.001); the
values themselves keep full precision (scores(fit)$ability, as.data.frame(), dplyr).
input_data() takes matrices or column names of a data frame; time is in seconds (the package
takes the log; use log_time = for logged times), cov gives centred covariates. It checks the
coding (0/1), missing values, logged times and uncentred covariates, and says what to do.
ecm() gives estimates, standard errors, logLik(), AIC(), BIC() and likelihood-ratio tests
(anova()).
gibbs(dat, model = "cross", init = fit) starts the chains at that solution, so a short burn-in
is enough.
f0 <- ecm(dat, model = "null") # nested: equal rho for every item
anova(f0, fit) # npar, logLik, AIC, BIC, LR test
#> npar logLik AIC BIC Chisq Df Pr(>Chisq)
#> f0 42 -8458.130 17000.260 17177.274
#> fit 51 -8309.032 16720.064 16935.009 298.196 9 <0.001 ***
g <- gibbs(dat, model = "cross", init = fit, n_iter = 2000, n_chain = 3, seed = 1)
g
model = "cross" (rtirt_cross)
500 persons x 10 items
3 chains x 2000 iterations (burn-in 1000, thin 1); 4.2 s
Reliability: ability 0.738, speed 0.970
Marginal: DIC 16720.571 (pD 51.053), WAIC 16721.320, LOOIC 16722.056
Warning: only 98.0% of 51 parameters have ESS > 400 and R-hat < 1.05 (max R-hat 1.021, min ESS 318); run longer (n_iter) or see convergence()
Use summary(), coef(), scores(), estimates(), convergence(), fit_indices(), plot().
print() warns when a parameter fails the criteria of convergence() (ESS > 400 and R-hat < 1.05),
or when one chain gives no R-hat; this short run is for illustration. In practice start from the
defaults (n_iter = 5000, n_chain = 4) and raise n_iter (e.g. to 10000) until print() gives
no warning. While sampling, gibbs() prints a line when it starts and one per chain with the time
left (control = gibbs_control(verbose = FALSE) turns them off).
estimates(g, pars = "^rho")[1:3, ] # posterior mean, SD, 95% interval, R-hat, ESS
# A tibble: 3 x 8
parameter est sd q025 q975 rhat ess sig
rho[q1] -0.094 0.062 -0.216 0.023 1.000 1717
rho[q2] 0.076 0.062 -0.048 0.193 1.001 1765
rho[q3] 0.164 0.062 0.045 0.286 1.002 1666 *
convergence(g)
#> 98.0% of 51 parameters with ESS > 400 and R-hat < 1.05
plot(g) # item parameters with 95% intervals; plot(g, "trace") for trace plots
| Need | ecm() |
gibbs() |
|---|---|---|
| estimates and SEs | 1-2 s | tens of seconds |
| AIC, BIC, likelihood-ratio test | yes | marginal DIC, WAIC, LOO instead |
| small samples, parameters near a bound | prior = "map" |
priors regularize |
| quantile (asymmetric Laplace) models | no | yes |
| full uncertainty of derived quantities | Wald | posterior draws |
equations() writes the model with the parameter names of the tables; algorithm() lists every
step of the estimator with the prior constants substituted. Both print Unicode and give LaTeX with
format = "latex" (an aligned block for a paper).
equations(g)
Measurement:
logit P(yᵢⱼ = 1) = aᵢ(abilityⱼ − bᵢ) i ∈ {q1, q2, …, q10}
log tᵢⱼ = λᵢ − speedⱼ − ρᵢ·abilityⱼ + εᵢⱼ i ∈ {q1, q2, …, q10}
εᵢⱼ ∼ N(0, σᵢ²)
Person variables:
abilityⱼ ∼ N(0, 1)
speedⱼ ∼ N(0, s) s = var_speed
Priors:
aᵢ ∼ N⁺(1, 1)
bᵢ ∼ N(0, 3²)
λᵢ ∼ N(0, 10²)
σᵢ ∼ half-t(3, 1)
ρᵢ ∼ N(0, 1)
√s ∼ half-t(3, 1)
algorithm(g) # the full conditionals of one Gibbs iteration, in order (excerpt)
# algorithm(g, detail = FALSE): the steps only
Step 1. Cross-relations:
ρᵢ | · ∼ N(mᵢ/pᵢ, 1/pᵢ) normal regression of the log times on theta
pᵢ = 1 + Σⱼ wᵢⱼθⱼ²
mᵢ = Σⱼ wᵢⱼθⱼ(λᵢ − ζⱼ − log tᵢⱼ)
Step 3. Accuracy items:
ωᵢⱼ | · ∼ PG(1, aᵢ(θⱼ − bᵢ)) draw_omega()
bᵢ | ω, a, θ ∼ N(mᵢ/pᵢ, 1/pᵢ) draw_items_pg(): b first
pᵢ = 1/9 + aᵢ² Σⱼ ωᵢⱼ
mᵢ = aᵢ Σⱼ (aᵢωᵢⱼθⱼ − κᵢⱼ)
aᵢ | ω, b, θ ∼ N⁺(mᵢ/pᵢ, 1/pᵢ) then a (truncated normal)
pᵢ = 1 + Σⱼ ωᵢⱼ(θⱼ − bᵢ)²
mᵢ = 1 + Σⱼ κᵢⱼ(θⱼ − bᵢ)
...
Step 9. Cross-relation move:
(ρᵢ, ζⱼ) ← (ρᵢ − δ, ζⱼ + δθⱼ) leaves the likelihood unchanged
δ | · ∼ N(m/p, 1/p)
The note on the right of each step names the building block (draw_*()) that runs it, so a new
sampler can be composed in R with the same functions (?blocks, inst/examples/ex07_blocks.R).
quantile = replaces the mean of the log times by a quantile (asymmetric Laplace working
likelihood). A vector of levels fits one model per level, by gibbs() or by ecm() (maximum
working likelihood with sandwich standard errors).
qs <- gibbs(dat, model = "cross", quantile = c(0.1, 0.5, 0.9), n_iter = 2000, n_chain = 2, seed = 2)
qs
model = "cross" (rtirt_cross) at quantiles 0.1, 0.5, 0.9
parameter q=0.1 q=0.5 q=0.9
rho[q1] -0.179 -0.085 -0.206
rho[q2] 0.100 0.092 0.140
rho[q3] 0.209 0.193 0.287
rho[q4] -0.051 0.044 0.099
rho[q5] 0.352 0.348 0.442
...
Warning (q = 0.1): only 39.2% of 51 parameters have ESS > 400 and R-hat < 1.05 (max R-hat 1.118, min ESS 54); run longer (n_iter) or see convergence()
Warning (q = 0.5): only 88.2% of 51 parameters have ESS > 400 and R-hat < 1.05 (max R-hat 1.057, min ESS 253); run longer (n_iter) or see convergence()
Warning (q = 0.9): only 41.2% of 51 parameters have ESS > 400 and R-hat < 1.05 (max R-hat 1.026, min ESS 110); run longer (n_iter) or see convergence()
The chains here are short for the sake of the example (the defaults are n_iter = 5000,
n_chain = 4); the print warns for each quantile level where some parameter has ESS <= 400
or R-hat >= 1.05 (convergence(qs) lists them).
plot(qs) # rho of each item against the quantile level, with 95% intervals
Deprecated.
block_model(), its steps (step_gibbs(),step_mala(),step_shift()),vi()andmml()still work, with a warning once per session, but get no new features. New models are added to the built-in menu ofecm()andgibbs()instead.
block_model() declares a model in the style of JAGS or PyMC. A list of steps names the blocks;
the same model and steps go to ecm() (conditional maximization), gibbs() (conditional draws)
and vi() (coordinate-ascent variational inference). Polya-Gamma variables, Jacobians,
truncation and dimensions ([person], [item]) are handled for you.
set.seed(11)
N <- 500; K <- 10
a_true <- seq(0.6, 2, length.out = K); b_true <- seq(-2.5, 2.5, length.out = K)
theta_true <- rnorm(N)
Y <- matrix(rbinom(N * K, 1, plogis(sweep(outer(theta_true, a_true), 2, a_true * b_true))), N, K)
m <- block_model(Y = Y,
theta[person] ~ normal(0, 1),
a[item] ~ normal(1, 1, lower = 0),
b[item] ~ normal(0, 3),
Y ~ bernoulli_logit(a * (theta - b)))
m
Data (j = person, N = 500; i = item, K = 10):
logit P(Yⱼᵢ = 1) = aᵢ(θⱼ − bᵢ)
Priors:
θⱼ ∼ N(0, 1)
aᵢ ∼ N⁺(1, 1)
bᵢ ∼ N(0, 9)
steps <- list(step_mala(a, b), step_gibbs(theta, a, b), step_shift(theta = 1, b = 1))
fit_e <- ecm(m, steps) # posterior mode and SEs (prior = "ml" for ML)
fit_g <- gibbs(m, steps, n_iter = 2000, seed = 1) # posterior draws
fit_v <- vi(m, list(step_gibbs(theta, a, b))) # CAVI: one factor q per block
fit_e
ECM for the declared model: converged after 18 iterations (2 rounds of the grid), log posterior -2408.095
parameter est se
a[1] 0.500 0.136
a[2] 0.794 0.165
a[3] 0.939 0.168
...
b[9] 1.849 0.191
b[10] 2.527 0.339
The three engines on the discriminations (true is the simulated value):
parameter true ecm ecm_se gibbs gibbs_sd vi vi_sd
a[1] 0.600 0.500 0.136 0.484 0.125 0.497 0.032
a[2] 0.756 0.794 0.165 0.756 0.185 0.789 0.048
a[3] 0.911 0.939 0.168 0.936 0.174 0.931 0.065
a[4] 1.067 1.049 0.179 1.036 0.177 1.040 0.071
a[5] 1.222 1.038 0.164 1.031 0.168 1.022 0.096
...
The CAVI means sit next to the posterior means; its SDs are too small, as for any mean-field
approximation, and vi() says so. Each engine has a check that it does what it claims:
check_ecm(m, steps) # every step raises the objective, Fisher identity, a maximum
check_sampler(m, steps = steps) # Geweke (2004) joint distribution test; method = "sbc" for SBC
check_vi(m) # the ELBO never decreases and the fit is a maximum
ECM check (10 starting values): all checks passed
check worst limit passed
each step raises Q(. | x) (relative drop) 0.00e+00 1e-08 TRUE
one iteration never lowers the objective (relative drop) 0.00e+00 1e-08 TRUE
Fisher identity: gradient of the objective = gradient of Q (relative difference) 7.16e-09 1e-04 TRUE
ECM stops at a maximum: gain of a general optimizer (relative) 1.13e-15 1e-06 TRUE
ECM stops at a maximum: largest |gradient| 8.94e-06 1e-02 TRUE
The tolerances of these checks (worst, limit) stay in scientific notation; every other
table prints 3 decimals.
A declared response-time model reproduces the built-in one. Speed enters only normal parts
linearly, so ecm() integrates it exactly and uses quadrature for ability alone:
mrt <- block_model(Y = unname(sim$Y), logT = unname(sim$log_t),
theta[person] ~ normal(0, 1),
zeta[person] ~ normal(0, sqrt(v)),
v ~ half_t(3, 1),
a[item] ~ normal(1, 1, lower = 0),
b[item] ~ normal(0, 3),
lambda[item] ~ normal(0, 10),
rho[item] ~ normal(0, 1),
sigma2[item] ~ half_t(3, 1),
Y ~ bernoulli_logit(a * (theta - b)),
logT ~ normal(lambda - zeta - rho * theta, sqrt(sigma2)))
c(builtin = as.numeric(logLik(fit)), declared = ecm(mrt, prior = "ml")$loglik)
#> builtin declared
#> -8309.032 -8309.032
step_gibbs(theta, cut = "logT") gives a cut (the times never feed back into ability), and
score_test(fit_e, add = Y ~ delta[item] * z) tests an added term from the null fit alone.
group <- rbinom(500, 1, 0.5) == 1 # a grouping or continuous moderator (gender, position, total time, ...)
score_test(fit, group, by_item = TRUE) # score-based test of parameter invariance (DIF), per item with Holm
Score-based test of parameter invariance (LM test, factor with 2 levels moderator, block decorrelation)
parameters n statistic p.value p.holm
accuracy 20 25.072 0.199
q1 2 2.747 0.253 1.000
q2 2 0.317 0.853 1.000
...
ci_test(fit) # conditional independence of responses and times (LM test)
#> Conditional independence of responses and response times (van der Linden & Glas, 2010)
#> all items: LM = 13.224, df = 10, p = 0.211
person_fit(g) # response-time person fit (Marianti et al., 2014), from a Gibbs fit
preknowledge_test(f0, compromised = c("q1", "q2")) # item preknowledge (Sinharay, 2020); a fit of "null", "latreg" or "latent"
score_test() also returns the casewise scores through estfun(), so the fit can go straight
into strucchange or partykit.
rtirt_priors(a = c(mean = 1, sd = 0.5), rho = c(sd = 0.3)) # the Prior column of summary() shows them
ecm(dat, model = "cross", prior = "map", priors = rtirt_priors(a = c(1, 0.5))) # posterior mode in the item parameters
gibbs(dat, model = "cross", init = fit, control = gibbs_control(collapse = TRUE)) # partially collapsed Gibbs: same posterior, higher ESS
Priors (gibbs(); item parameters also for ecm(prior = "map")):
a N+(1, 0.5^2)
b N(0, 3^2)
lambda N(0, 10^2)
sigma2t SD ~ half-t(3, 1)
beta N(0, 1)
rho N(0, 0.3^2)
var_speed SD ~ half-t(3, 1)
cov N(0, 1)
The defaults are common data-independent priors and every full conditional stays conjugate.
collapse = TRUE integrates speed and the Polya-Gamma variables out of the slow updates: on the
example above the effective sample size of the discriminations goes from 318-1099 to 939-1447
with the same number of iterations.
model = |
model |
|---|---|
"irt" |
2PL (or 1PL); with covariates, a latent regression of ability |
"null" |
2PL + log-normal response times, ability and speed correlated (van der Linden, 2007) |
"latreg" |
... with latent regressions of ability and speed on covariates |
"latent" |
speed regressed on ability and covariates; quantile = for quantile regression |
"cross" |
a cross-relation of ability with each log response time; quantile =
|
"mixture" |
the mixture model of Wang and Xu (2015); planned, not available yet |
All are fitted by ecm(dat, model = ...) or gibbs(dat, model = ...); the quantile versions by
gibbs() only. The details are arguments: itemtype ("2pl", "1pl"; "grm" planned),
speed_var, quantile, cov; the rarely changed engine settings go in ecm_control() and
gibbs_control(). Tables (estimates(), scores(), the tables of summary()) are tibbles.
The old constructors (mlirt(), rtirt_null(), rtirt_latreg(), rtirt_latent(),
rtirt_cross()) and gibbs(fit, init = "ecm") still work, with a deprecation warning; the
Julia aliases (GibbsRtIrtCross(), ...) are kept as they were. The
accessors summary(), coef(), estimates(), scores(), reliability(), convergence(),
fit_indices() and plot() work on both kinds of fit.
vignette("birtRcpp"): getting started (installed with build_vignettes = TRUE; also
on GitHub).inst/tutorial/birtRcpp_intro.qmd). The second tutorial,
birtRcpp_blocks.qmd (declaring models and composing samplers and EM algorithms), covers
the deprecated block_model().system.file("examples", package = "birtRcpp"): 13 scripts, from the quick start to testlet,
multidimensional, LNIRT, variable-speed and preknowledge models composed from the blocks.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.