Getting started with birt

1. Model and fit (ELGM)

library(birt)
d <- sim_rtirt(N = 500, K = 10, seed = 1)        # y1..y10 (0/1), t1..t10 (log RT), x1, x2
items <- paste0("y", 1:10)
times <- paste0("t", 1:10)

model <- rtirt_syntax(items, times)    # van der Linden (2007)
cat(model)                             # ability =~ y's; speed =~ -1*t's; ability ~~ speed

fit <- birt(model, data = d)           # ELGM (default), about 20 s
fit                                    # log marginal likelihood, inner / outer quadrature
summary(fit)                           # posterior means, SDs, 95% intervals, priors

• The first =~ line of a factor gives its measurement loadings; later lines add cross-loadings (e.g. rtirt_syntax(items, times, cross = "free"): ability on each log RT, Bolsinova & Tijmstra, 2018).
• ability has variance 1 (it has free loadings); speed, whose loadings are all fixed, has a free variance; they are uncorrelated unless ability ~~ speed is added.
• Priors: defaults from dpriors() (dp = dpriors(cross = "dnorm(0, 100)") changes a class), or prior("dnorm(0, 100)")*t1 for one parameter (JAGS precision parameterization).

2. What the object holds

estimates(fit)          # lavaan-like table: est, se (posterior SD), lower, upper, prior
scores(fit)             # EAP and posterior SD per person and latent variable
reliability(fit)        # var(EAP) / (var(EAP) + mean PSD^2)
summary(fit)$items      # a, b (IRT parameterization), RT intercepts, loadings, variances
fit_indices(fit)        # log marginal likelihood
convergence(fit)        # optimizer, outer nodes, quadrature check
equations(fit)          # the model as equations (format = "latex" for a paper)
algorithm(fit)          # how the fit was computed, step by step

These tables and numbers print with 3 decimals (p-values below 0.001 as <0.001, counts without decimals); the returned values keep full precision, e.g. estimates(fit)$est.

3. Model comparison

fc <- birt(rtirt_syntax(items, times, cross = "free"), d)      # item cross-loadings
fi <- birt(sub("ability ~~ speed", "", model, fixed = TRUE), d)  # independent
compare(correlated = fit, cross = fc, independent = fi)          # logML, Bayes factors

# maximum likelihood: AIC and BIC
ml  <- birt(model, d, control = rtmb_control(method = "aghq"))
mlc <- birt(rtirt_syntax(items, times, cross = "free"), d, control = rtmb_control(method = "aghq"))
compare(correlated = ml, cross = mlc)

Marginal likelihoods depend on the priors; report them with the priors.

4. Distributional models and quantiles

Instead of normal residuals, shash() gives sinh-arcsinh residuals (skewness, tail weight); V(y) ~ z in the syntax moderates a log scale. quantiles(fit, p) then gives the effects on every conditional quantile from one fit, with simultaneous bands over p; quantile_test() tests whether an effect changes with p.

reg <- rtirt_syntax(items, times, structure = "path", cov = c("x1", "x2"))  # speed ~ ability + x1 + x2
fs <- birt(reg, d, family = list(speed = shash()), control = rtmb_control(k_hyper = 13))

m5 <- paste(rtirt_syntax(items, times, cross = "free"),
            sprintf("V(%s) ~ ability", paste(times, collapse = " + ")), sep = "\n")
f5 <- birt(m5, d)
quantiles(f5, p = c(0.1, 0.5, 0.9))
quantile_test(f5, p = c(0.1, 0.5, 0.9))
cond_reliability(f5)    # reliability of ability given its value: accuracy, RT, both
rt_information(f5)      # where the RT information lies; what a time limit keeps

If ELGM warns that the posterior is too wide for the outer grid along a hyperparameter (more mass on its two outermost nodes than a Gaussian posterior 1.5 times wider would put there, and at least 2.5%), check whether that parameter is weakly determined and increase rtmb_control(k_hyper = ).

5. Maximum likelihood tools

score_test(ml, d$x1 > 0, by_item = TRUE)   # measurement invariance from one fit
score_test(ml, rowSums(d[times]))          # along the total log time
qs <- quantile_score(birt(m5, d, control = rtmb_control(method = "aghq")), p = 0.9)
qs                                         # influence functions, robust SEs
score_test(qs, d$x2)
ipcr(ml, d[c("x1", "x2")])                 # which parameters change with which covariates
rtmb_code(ml)                              # a stand-alone RTMB script

ipcr() regresses each person’s contribution to every estimate on the covariates (individual parameter contribution regression; Arnold et al., 2020), so one call screens all parameters against several covariates. dif_tree(model, data, partition) (package partykit) is deprecated.

6. Ordinal indicators moderated by a response time

With ordered = TRUE the indicators are ordinal (graded response model) and person-level moderators enter as moderation statements on conditional moments, as in MNLFA: E(y1 + y2) ~ z shifts the thresholds (uniform DIF), E(y1) ~ z:f moderates the discrimination (nonuniform DIF), E(f) ~ z and V(f) ~ z the latent mean and variance. hgrm() writes this syntax from a few options; both fit by ELGM by default.

dh <- sim_hgrm(N = 600, J = 8, seed = 1)
dh$logT <- as.numeric(scale(dh$logT))
m <- 'SA =~ y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8
      E(y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8) ~ pa*logT:SA
      E(y1 + y2) ~ logT
      E(SA) ~ logT'
f_ord <- birt(m, dh, ordered = TRUE)
summary(f_ord)
ordinal_information(f_ord)

h0 <- hgrm("SA =~ y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8", dh, moderators = "logT",
           control = rtmb_control(method = "aghq"))
score_test(h0, dh$logT, by_item = TRUE)    # DIF screen, item by item

7. The JAGS engine (optional)

fit_ssp <- birt(rtirt_syntax(items, times, cross = "ssp"), d, engine = "jags")
summary(fit_ssp)$selection    # spike-and-slab: inclusion probabilities, Bayes factors
compare(M2b = fit_ssp, M1 = birt(model, d, engine = "jags"))   # LOOIC, ELPD
ppc(fit_ssp)                  # posterior predictive checks incl. conditional dependence
jags_code(fit_ssp)            # the JAGS model
q <- birt(reg, d, family = list(speed = ald(c(0.1, 0.5, 0.9))), engine = "jags")

8. Notes

• Covariates are used as given (standardize them) and may not be missing; indicators may be missing at random.
• The log-likelihood of response times is that of the log times (no Jacobian), the same for log_rt = TRUE and for data that are already logged; do not compare it with models for times in seconds.
• With ALD indicators the JAGS fit indices are slow; use fit_indices = FALSE and add_fit_indices(fit, n_draws = 200) later. Posterior intervals under ALD working likelihoods can be too narrow (Yang, Wang & He, 2016).

References

Arnold, M., Oberski, D. L., Brandmaier, A. M., & Voelkle, M. C. (2020). Identifying heterogeneity in dynamic panel models with individual parameter contribution regression. Structural Equation Modeling, 27, 613-628. https://doi.org/10.1080/10705511.2019.1667240

Bolsinova, M., & Tijmstra, J. (2018). Improving precision of ability estimation: Getting more from response times. British Journal of Mathematical and Statistical Psychology, 71, 13-38.

Stringer, A., Brown, P., & Stafford, J. (2023). Fast, scalable approximations to posterior distributions in extended latent Gaussian models. Journal of Computational and Graphical Statistics, 32, 84-98. https://doi.org/10.1080/10618600.2022.2099403

van der Linden, W. J. (2007). A hierarchical framework for modeling speed and accuracy on test items. Psychometrika, 72, 287-308.

Yang, Y., Wang, H. J., & He, X. (2016). Posterior inference in Bayesian quantile regression with asymmetric Laplace likelihood. International Statistical Review, 84, 327-344.