estimates(), scores(), fit_indices(), convergence(), compare() and the tables of
summary() return tibbles. The values are unchanged and as.data.frame() gives the old
form; tibbles print 3 significant digits (options(pillar.sigfig = 7) shows more) and the
first 10 rows of long tables (print(x, n = Inf)).ipcr(): individual parameter contribution regression (Arnold et al., 2020) for maximum
likelihood fits. It regresses each person's contribution to every estimate on all
covariates at once, continuous and categorical, from one fit; the slopes are shrunk toward
0 (one-step approximation), so it screens and the moderated model gives effect sizes.
dif_tree() is deprecated in its favour: it still works and warns once per session
(options(birt.deprecation = "always") warns every time).Error in birt(model, d),
Error in hgrm(...), Error in compare(...)), also from inside a function of the user or
lapply(), instead of an internal helper (Error in build_spec(model, data, ...),
Error in parse_model(model), Error in hgrm_syntax(...)); the messages are unchanged.plot() of an RTMB fit without quantile effects that change with p (no V() statement or
shash() residual) draws the person scores (type = "scores") instead of only a message,
and print() points to plot() for the quantile effects only when there are such effects.
plot(fit, type = "information") now works for binary and continuous indicators: the test
information about the first latent variable from all indicators, the accuracy items and the
response times (from cond_reliability()). The RTMB plots have their own help page,
?plot.birt_rtmb, linked from ?plot.birt.estimates() has the same column names for every fit, as as.data.frame(summary(fit)):
est, se (posterior SD for ELGM and JAGS, standard error for maximum likelihood), lower,
upper (and z, pvalue for maximum likelihood; se_adj, lower_adj, upper_adj for ALD
fits). The names of 0.1.0 (sd, q025, q975, ci.lower, ci.upper) still work with $,
[[ and [, with a warning once per session.log sd(speed), atanh cor(ability, speed) instead of lsd2, atr; the prior
of a variance row reads on SD: dt(0, 1, 3) T(0,) (it is a prior on the SD); the ELGM
parameter table ends with a note that there are no p-values and how to get them (maximum
likelihood); convergence() of an ELGM fit labels the inner rounds -2 log posterior
(they include the prior) instead of -2logL.print() of ipcr() shows the slopes with Holm-adjusted p < 0.05, sorted by p (all of
them with print(x, all = TRUE); x$coefficients is unchanged), and says how to read a
slope on the internal scale (log loadings, log SDs, atanh correlation).sim_rtirt() and sim_hgrm() have seed = NULL by default: without a seed they draw from
the session's random numbers (as rnorm() does), so set.seed() and repeated calls in a
simulation give new data. Code that relied on the old default seed = 1 gets the same data
with seed = 1; the examples, tests, README, vignette and tutorial now give it.compare() of ELGM fits returns BF_best as a number (it was text); the printed table is
as before.?birt, the README and the tutorial state the rule for cross-loadings as
the code applies it (a factor split over consecutive =~ lines keeps its measurement
loadings; =~ lines written after a line of the other factor are cross-loadings), list
E(y3) ~ x2 and V(t1) ~ f in the first syntax table, start the examples of ?birt with
the default ELGM fit, and say what progress does for the RTMB engine. The README starts
with the questions birt answers and a "question -> function" table, and says that the
default fit has no p-values and how to get them; the tutorial opens with the same table and
a short path for a first result, installs from GitHub as the README does, and folds the
outer-grid details into an "advanced" box. ?score_test and ?ipcr have the concepts DIF,
measurement invariance and response time models for help.search().6e-17 prints 0.000; numbers of 1e10 or
more, e.g. the bound of a diverging estimate, print as 5.950e+99):
print() and summary() of every fit (ELGM, maximum likelihood, JAGS), estimates(),
scores(), fit_indices(), quantiles(), ordinal_information(), convergence(),
compare(), reliability(), coef(), cond_reliability(), quantile_test(),
quantile_score(), score_test(), ipcr(), ppc() and check_sbc(). Counts (N, npar, df,
nodes, ESS) print without decimals, p-values below 0.001 as <0.001, Bayes factors above
1e5 as e.g. 6.26e+11. Only the display is rounded: the returned tables are tibbles (or data
frames) of class birt_tbl and the named numbers of reliability() and coef() have class
birt_values (see ?birt_print); x$est, arithmetic, dplyr verbs and as.data.frame() use
the full values. The tables print all rows when they have at most 20 (otherwise the first 10;
print(x, n = Inf) for all), and wide tables continue below. The tables inside returned lists
(quantile_information()$best, rt_information()$bands, quantile_test()$test,
ppc()$accuracy, check_sbc()$test, ...) are birt_tbl too.inst/tutorial/birt_tutorial.qmd) starts with a chapter for readers new to R
and IRT: installing R, RStudio and birt; the few R basics needed (<-, functions and their
arguments, data frames, head(), reading an error message); what the data look like; ability,
difficulty and discrimination with an item curve; why response times help; and a first
complete analysis explained line by line. The README points to it.~~ or on the left of ~ that is neither a latent variable nor an indicator
(e.g. abilty =~ ... and then ability ~~ speed) is reported as a likely typo with the
closest defined name, instead of "residual covariances between indicators are not
supported"; a covariate missing from data names the closest column.rtmb_control(method = "ml") (also "ML", "mle"; the same as
"aghq"). birt(..., method = "ml") says to write control = rtmb_control(method = "ml"),
other unknown arguments of birt() are refused with a typo hint, control = list(method = "ml") is passed to rtmb_control(), and control = "aghq" says what to write.ipcr() refuses covariates named like indicators of the model (ipcr(fit, d) with the
whole data frame) and gives the impact message of score_test() (new argument
impact_check): a covariate related to speed shows up as slopes on every response time.ordered =) get a warning with the ordered = call to write. The warnings of the fit name
parameters as the syntax does (log(ability =~ y1) instead of lpos).log_rt: times that look already log-transformed get a warning; times <= 0 say how many
and, when negative, that the data may already hold log times.ability ~= y1, ability =- y1, e(y3) ~ x2, a line
ending in +) names the line of the model, says what is wrong (=~ misspelled, E() /
V() in lower case, a missing term) and lists the operators, instead of the message of
lavaan or of R's parser. A line that lavaan only warns about ("invalid formula") is an error.ability ~ gender) gets the 0/1 recoding to write (dummy variables for more than two
values) and missing covariate values are reported separately, with the rows; a non-numeric
indicator gets the 0/1 or integer recoding; y2 y3 without + asks about the missing +;
a formula as model says to quote the syntax; score_test(fit, "x1") says to give the
values (data$x1).scores() no longer prints in scientific notation.rtirt_syntax() is -1).First release. birt fits item response and response time models written in lavaan-style syntax, with persons j and items i (y_ij is the response of person j to item i).
birt(model, data) fits in 'RTMB' (engine = "rtmb") with
rtmb_control(method = "elgm"): approximate Bayesian inference for an extended latent
Gaussian model (Stringer, Brown & Stafford, 2023). The persons are integrated out by
per-person adaptive Gauss-Hermite quadrature, the other parameters by an outer adaptive
quadrature at the posterior mode, with the priors of dpriors(). It gives posterior means,
SDs and intervals, EAP scores and the log marginal likelihood for compare() (Bayes
factors). convergence() reports the outer grid, and a warning flags a posterior that is
too wide for the grid along a hyperparameter. Checked: with a Gaussian outer layer the log
marginal likelihood equals the Laplace value at the mode (to 1e-8). hyper = "auto" puts
outer quadrature on at most 3 hyperparameters (SHASH shape first) and says which it treats as
Gaussian: with a 4th direction (5^4 = 625 nodes) the outer step took 5 times longer for
posterior means within 0.034 SD. Timing (N = 1000-2000, 20-40 items, two latent variables,
normal or SHASH; 81 outer nodes): ELGM took 1.4-2.0 times the maximum likelihood fit.rtmb_control(method = "aghq"): the same marginal likelihood
maximized, with standard errors, AIC / BIC and the casewise scores used by score_test(),
quantile_score() and dif_tree(). Checked against brute-force integration (to 1e-4) and
against mirt's GPCM and nominal models (-2logL within 1e-5).engine = "jags" (R2jags and rjags in Suggests) samples the person
latents with the parameters (conditional likelihood): spike-and-slab priors, ald()
quantile models, DIC / WAIC / PSIS-LOO (marginal likelihood per draw by adaptive
Gauss-Hermite quadrature; Liu & Pierce, 1994), posterior predictive checks (ppc()).
Starting values come from the data, and each chain starts the scale of each latent variable
at a different place so that R-hat can see slow mixing along the scale.equations() and fit$likelihood say whether a fit uses the conditional
or the marginal likelihood (Merkle, Furr & Rabe-Hesketh, 2019).birt(): one or two latent variables, binary, ordinal and continuous (e.g. log response
time) indicators, latent regression, factor correlation, labels for equality constraints,
prior() modifiers and fixed values. rtirt_syntax() writes the usual response time IRT
models (van der Linden, 2007; cross-loadings, path or correlation structure, covariates).dpriors(): a hierarchical lognormal prior for the measurement loadings of
each latent variable (log lambda ~ N(mu_f, 1), mu_f ~ N(0, 10)), and priors for thresholds,
steps, moderation effects, SDs and the correlation; the prior of every parameter is shown in
summary().V(y) ~ z, all continuous indicators are
family = list(.continuous = ). Ambiguous or ignored specifications are errors (a label
with prior() on one term, labels shared across parameter types, a label times a number).log_rt = TRUE for raw seconds, id = for person identifiers, 1/2 coding of
binary items, and a warning when continuous indicators look like raw times.ordered = TRUE (or names) makes integer indicators with 3-10 categories ordinal, combined
freely with binary and continuous indicators, by every engine. itemtype = chooses the
graded response model (default), the generalized partial credit model (Muraki, 1992) or the
two-parameter partial credit model (Yu, 1991).E(y) ~ z
(threshold or intercept shift), E(y) ~ z:f (log loading), V(y) ~ z (log residual SD),
E(f) ~ z and V(f) ~ z, with priors, labels, anchors (0*z) and spike-and-slab
(prior("ssp"), JAGS). Moderation of the partial credit models is experimental (not yet
checked by simulation-based calibration).hgrm() writes the moderated graded response model from the options a, b, impact and
anchor and fits it by ELGM (or maximum likelihood for score_test()).fit$dif, fit$impact, fit$precision (score precision by moderator tercile),
ordinal_information() and plot(fit, type = "dif" / "information").shash() residuals (sinh-arcsinh; Jones & Pewsey, 2009) for continuous indicators and
latent variables, with skewness moderation shash(skew = ).quantiles(): quantile effects of the predictors at every level p, with a simultaneous 95%
band over the levels (sup-t; Montiel Olea & Plagborg-Moller, 2019) and their joint
covariance (posterior covariance with ELGM). Checked against numerical derivatives of the
conditional quantiles (to 1e-6).quantile_test(): does an effect change with p? It does exactly when the predictor
moderates the scale or skewness of the target, so the test is the Wald test of those
parameters (Koenker & Bassett, 1982). Checked: 5 of 56 rejections at .05 under the null.quantile_score(): casewise influence functions of the quantile effects (Hampel, 1974),
robust standard errors, case influence, conditional or unconditional on covariates
(Graubard & Korn, 1999), and score_test(quantile_score(fit, p), z) for invariance of an
effect. Checked: robust SEs agree with 100 bootstrap fits (0.056 / 0.049 vs 0.058 / 0.051).quantile_information(): which quantile regression of the response times says the most about
the first latent variable: the information of the tau-th regression (the Godambe information of
its check-loss estimating function, f(q)^2 beta(tau)^2 / (tau (1 - tau)); Godambe, 1960;
Koenker, 2005), speed projected out, and its efficiency against the whole RT distribution.
Checked: beta equals the cross-loading and the item values equal the closed form for normal
residuals; for SHASH beta and f(q) match the quantile function (to 1e-5) and the information
agrees with Monte Carlo of the dichotomized times (within 3%).ald() (JAGS): Bayesian quantile regression with the asymmetric Laplace working likelihood
(Yu & Zhang, 2005; Kozumi & Kobayashi, 2011), measurement or latent level, with intervals
adjusted as in Yang, Wang & He (2016).cond_reliability(): conditional reliability of the first latent variable,
I / (I + 1 / var) (Nicewander, 2018), with speed as a nuisance, the information of the
response-time location, scale and skewness, and delta-method standard errors. Checked
against the 2PL information, a closed form and Monte Carlo. Ordinal items (GRM, GPCM,
TPPCM) and moderated loadings and shifts are included at the moderator values at (a
censored fit: the information under its limits, Escobar & Meeker, 1998). Checked
against ordinal_information() (to 1e-8) and the numerical Fisher information of the
likelihood code (to 1e-6; censored and SHASH items to 1e-5).rt_information() (non-core): where in the response-time distribution the information lies (bands of
RT quantile levels; Zheng & Gastwirth, 2000) and how much is kept when times beyond a limit,
or fast responses, are censored (Efron & Johnstone, 1990). Checked: the bands sum to the RT
information (to 3e-7) and the censored information agrees with Monte Carlo (within 0.1%)
and with the curvature of the censored likelihood of censor = (to 1e-5). Moderated
loadings and shifts and ordinal items are included, as in cond_reliability().ordinal_information(): test information of ordinal indicators at chosen moderator values.score_test(): score-based tests of parameter invariance (Merkle & Zeileis, 2013) from one
maximum likelihood fit, with casewise scores from estfun(): DIF along a numeric, ordinal
or categorical moderator, per item with Holm adjustment, a message when the moderator is
related to the scores (impact), and cluster-level moderators. Valid for the measurement
parameters in a joint accuracy and response-time model.score_test(),
quantile_score() and estfun() (and dif_tree() with an ELGM control) stop with the
instruction to refit with rtmb_control(method = "aghq"). At the posterior mode the scores
sum to minus the gradient of the log prior and MAP versions need centred scores with
covariance-dependent error rates (Debelak, Pawel, Strobl & Merkle, 2022); the Bayesian
route for one effect is compare() of ELGM fits with and without it (Bayes factor).dif_tree(): DIF trees by model-based recursive partitioning (Zeileis, Hothorn & Hornik,
2008; Strobl, Kopf & Zeileis, 2015) for binary and ordinal items, with latent mean
differences modelled in each node (impact = TRUE). Numeric partitioning variables are cut
into nbins = 10 quantile groups and tested as ordinal (Merkle, Fan & Zeileis, 2014):
a 2PL tree (N = 600) with a numeric variable took 26 s instead of 15 minutes for the
exhaustive search (nbins = NULL), and split at 6.01 for a true cut point at 6 (exhaustive:
6.05).p_incl, the learned prior inclusion probability and BF10 in summary(fit)$selection.check_sbc(): simulation-based calibration of the JAGS engine (Talts et al., 2018), with
an R simulator written independently of the generated JAGS code.convergence(): R-hat and ESS (JAGS); optimizer, gradient, quadrature accuracy and the
outer grid (RTMB).summary() returns the tables of a fit as one object ($fit, $items, $se,
$parameters with the prior of each parameter, $moderation, $selection, ...), printed
in lavaan / blavaan sections; as.data.frame() gives the parameter table.equations() writes the model as equations with the priors (Unicode, or LaTeX with
format = "latex"); algorithm() writes how the fit was computed, with the numbers of the
fit (JAGS samplers, AGHQ rounds, ELGM outer quadrature).jags_code() and rtmb_code() export a script that reproduces the fit, from a fit or from
syntax (with ordered =, itemtype =). rtmb_code() writes ordinal items (GRM, GPCM, TPPCM),
E() moderation and censored indicators; checked: the
script reproduces -2logL of the fit (to 1e-5) and its joint density equals birt's (to 1e-10).
scores(),
estimates(), posterior_draws(), compare(), plot().birt(..., censor = list(t1 = c(upper = 60))), with
limits on the scale of the data columns (log-transformed with the times under log_rt). A
value at or beyond a limit contributes log S(c | f) or log F(c | f) instead of the density
(Meeker & Escobar, 1998; Lee & Ying, 2015); normal and SHASH residuals, ELGM and maximum
likelihood (not JAGS). Checked: -logL equals brute-force integration (normal to 1e-13; SHASH:
the likelihood code on a fine grid to 1e-8); in a time-limit simulation (N = 1000, 10 items,
18% censored, 8 replications) the mean bias of the time intercepts, residual SDs and speed SD
is -0.0003, 0.0003 and 0.004, against -0.066, -0.063 and -0.079 when the censored times are
treated as observed.