NEWS

---

birtRcpp 0.3.0

• ecm() fits the quantile models: ecm(dat, model = "cross", quantile = 0.9) (and model = "latent") maximizes the marginal ALD working likelihood, with ability on the Gauss-Hermite nodes and speed integrated in closed form given ability. The E-step and the item step are those of the normal models; the response-time and structural parameters are updated by nlminb() with the analytic gradient. The standard errors are sandwich standard errors (H^-1 J H^-1, J from the scores of the persons), since the ALD is not the distribution of the data. A vector of levels gives an rtirt_qset of ECM fits (print(), coef(), estimates(), plot(), convergence()). print() and summary() say that the logLik, AIC and BIC are of the working likelihood; anova() refuses quantile fits from ecm() as from gibbs(); prior = "map" is for the normal models only.
• gibbs() of a quantile model starts at an ecm() fit at the same quantile (init = fit), or at its own ECM fit (init = "ecm", one per level for several levels); before, quantile models started only from random values.
• check_ecm() checks the ECM of a quantile model (one level): one iteration never lowers the working log-likelihood, the analytic score is its gradient, and the fit is a maximum. algorithm(fit) writes the ECM of a quantile model.

birtRcpp 0.2.0

• Built-in models are fitted from the data with a model name: ecm(dat, model = "cross"), gibbs(dat, model = "cross", init = fit) (fit an ecm() fit of the same data and model), with model = "irt", "null", "latreg", "latent" or "cross". quantile = belongs to gibbs() (latent and cross models). Rarely changed settings moved to control = ecm_control(...) and control = gibbs_control(...) (nodes, max_iter, tol, accelerate; n_burnin, n_thin, collapse, fit_indices, verbose). equations(dat, model = ...) and algorithm(dat, model = ...) write a model out before it is fitted.
• Deprecated (they still work and warn once per session; options(birtRcpp.deprecation = "always") warns every time): the constructors mlirt(), rtirt_null(), rtirt_latreg(), rtirt_latent(), rtirt_cross(); gibbs(fit, init = "ecm"); block_model() and step_gibbs(), step_mala(), step_shift(); vi(), check_vi(); mml(). The building blocks (draw_*(), irt_part(), normal_part(), ...) and the Julia aliases are not deprecated.
• estimates(), scores(), fit_indices(), compare_para() and the tables of summary() return tibbles (values unchanged; as.data.frame() gives the old form).
• Printed numbers have 3 decimals everywhere, in fixed notation (numbers of 1e10 or more, e.g. the interval bound of a diverging estimate, as 5.950e+99): print() and summary() of ecm() and gibbs() fits (print(fit) showed AIC, BIC, DIC, WAIC and LOOIC with 1 decimal), estimates(), scores(), fit_indices(), compare_para(), precis(), dic(), anova(), convergence(), score_test(), ci_test(), person_fit(), preknowledge_test(), check_sampler(), run_simulation(), coef(), reliability(), marginal_loglik() and the prints of the deprecated declared-model fits. Counts (persons, npar, N, Df, df, ESS, ...) have no decimals; p-values below 0.001 print as <0.001 (not < 2.22e-16; anova() keeps its significance stars); NA is blank. Only the display changes: the objects keep full precision (estimates(fit)$est, as.data.frame(), dplyr verbs). The tables are tibbles (or, for person_fit(), preknowledge_test() and dic(), data frames) of class rtirt_tbl (?rtirt_tbl), printed in the plain layout of summary() with all rows up to 40 (else the first 20; print(x, n = Inf) for all). In R Markdown and Quarto documents these tables print as at the console, also with df-print: kable (which printed them with kable() and 7 digits). Run times and the share of converged parameters keep one decimal; the tolerances of check_ecm() (1e-08) stay in scientific notation, since 3 decimals would show every check as 0.000. The regression tests compare printed text with each number checked against the rounding of both printings, so the deliberate change of format does not hide a changed value.
• The R code is rewritten in a small-example style and formatted with styler; regression tests compare every built-in model, through the old and the new interface, with a snapshot taken before the rewrite.
• Documentation for new users: the README starts with where to begin (a link to the intro tutorial on GitHub, readable without Quarto, and a "you are ... read ..." table), installs with build_vignettes = TRUE and says why vignette("birtRcpp") is otherwise "not found", and explains the message that vi masks utils::vi. ?birtRcpp lists the functions needed to start and appears in the help index. The intro tutorial is written for readers who know the 2PL (with a one-paragraph recap), has a runnable section "from a CSV file to the report" (long format, milliseconds, missing times, dummy-coded and centred covariates, ecm(), anova(), scores merged back), states at the start of appendix B that check_sampler() checks the sampler code and not the convergence of an analysis, and no longer refers to row numbers 41-43 of a tibble. The vignette sketches the same data steps and links the tutorial.
• The intro tutorial starts with a beginner's introduction (section 0, for readers who have barely used R and never fitted an IRT model): what the data look like, IRT and response times in plain words with an item characteristic curve, installing R, RStudio, the C++ tools and birtRcpp, the R basics the tutorial needs (<-, data frames, head(), arguments, reading an error), and a first analysis in ten lines with the four numbers to read first. Readers who know R and IRT are sent to section 1; the 2PL recap moved into section 0, and the glossary is its last part (0.8). The README points to it.
• algorithm(fit, detail = FALSE) lists only the steps of one iteration (Step 1. Cross-relations, ...), without the formulas; the default is unchanged.
• check_sampler(dat, model = "cross") and check_ecm(dat, model = ...) say that they check the sampler or ECM steps of a fit, with the call (check_sampler(ecm(dat, model = "cross", se = FALSE))) and, for check_sampler(), that the convergence of an analysis is checked with convergence(); they failed with "'arg' must be NULL or a character vector". The Julia constructors called with data only (GibbsRtIrtCross(dat)) say that the conditions come first and give the R call (gibbs(dat, model = "cross")), instead of "argument "Data" is missing".
• gibbs() progress (gibbs_control(verbose = TRUE), the default): a line as soon as sampling starts (model = "cross": sampling 4 chains x 5000 iterations, one after the other; chain 1 of 4 ...), one per chain with its time and an estimate of the time left (chain 1 of 4 done (13.8 s); about 41 s left), and one before the marginal log-likelihood for fit_indices() is computed. Before, nothing appeared until the first chain had finished. The draws are unchanged. print() of set_cond() says that its MCMC settings are for run_simulation() and the Julia aliases, not for gibbs(data, model = ...); the README and the tutorial recommend the same chain lengths (start from the defaults, lengthen n_iter until print() gives no convergence warning).
• anova(fit0, fit1, ...) compares fits of the same data, as for lm(): for ecm() fits the number of parameters, log-likelihood, AIC, BIC and the likelihood ratio test of each fit against the row above (rows ordered by the number of parameters; a warning when a larger fit has a lower log-likelihood, i.e. the fits are not nested); for gibbs() fits DIC, WAIC and LOOIC. Fits of other data, a model without response times next to one with, prior = "map" fits, quantile fits and a mix of ECM and Gibbs fits are refused with the reason. The README, the vignette and the intro tutorial use it instead of computing 2 * (logLik(f1) - logLik(f0)), f1$ecm$npar - f0$ecm$npar and pchisq() by hand.
• Printing: precis() returns its table (class rtirt_precis, printed in the plain layout) instead of printing it and returning it invisibly, so p <- precis(fit) prints nothing; summary() of an rtirt_qset likewise returns an object that prints. print(summary(fit), tables = "item") is an error with the closest table name (it printed only the header); estimates(fit, pars = "rho") that selects nothing warns and lists the parameters of the fit. plot(ecm_fit, "trace") says that an ecm() fit has no draws and how to sample from it (was "run gibbs() first"). New plot(fit, "structural") (estimates and 95% intervals of the regression coefficients, cor_ability_speed, b_ability, var_speed) and as.data.frame(fit) (the table of estimates(), for ggplot2); trace plots label the x axis "iteration". The title of the Gibbs $fit table says to compare models by WAIC or LOOIC and what DIC_complete is.
• input_data() messages: responses given as text or factors show the values found and the recoding ((values: "wrong", "correct"); recode them to 0/1, e.g. data[resp] <- lapply(data[resp], function(x) as.integer(x == "correct"))); missing times (and log times and covariates) are counted per person with the code to drop them, as missing responses were; times with a median above 1000 warn that they look like milliseconds. ?input_data explains that the unit of time only shifts the time intensities (lambda by log 1000 = 6.91 for milliseconds). ecm(cov = ) / gibbs(cov = ) given the values of a covariate instead of its name says to name covariates of the data or add them in input_data() (it printed every value).
• gibbs(init = fit) with a fit of other data, another model or other settings says what differs (init is an ecm() fit of model = "cross", but gibbs() was asked for model = "latreg", init was fitted to 200 persons x 10 items, gibbs() has 300 x 10, init was fitted with speed_var = "fixed", ...) and the call to use, instead of "init must be an ecm() fit of the same data, model, speed_var and itemtype". init with a quantile model says to drop it (it spoke of init = "ecm" when a fit was given); init = "ECM" is accepted, and init = TRUE gets a readable error instead of the one of match.arg(). gibbs(fit, model = "cross") and ecm(fit, model = "cross") say that the first argument is the data and how to start from the fit (gibbs(data, model = "cross", init = fit)); they failed with "$ operator is invalid for atomic vectors" or refitted with a warning.
• ecm_control() and gibbs_control() return classed settings (ecm_control, gibbs_control, with a print() method): ecm(control = gibbs_control()) and gibbs(control = ecm_control()) say which function the settings belong to (they failed with "argument is of length zero" or "collapse must be TRUE or FALSE"). A name that is not a setting is an error with the closest setting (gibbs_control(burnin = 50): unknown setting burnin of gibbs_control(); did you mean n_burnin?) or the function it belongs to (gibbs_control(n_iter = 200): an argument of gibbs()). A plain list of settings is still accepted.
• Model names: ecm() and gibbs() take the model name in any case (model = "Cross"); an unknown name gets the closest one (unknown model "crosss"; did you mean "cross"?), a class name its model name ("mlirt": did you mean "irt"?), an item type the hint that it belongs in itemtype =, and the message lists only the models that are available, one line each saying what the model is ("latent" speed regressed on ability and the covariates). sim_para(), sim_data() and run_simulation() take the same names ("irt"; "mlirt" still works, with the same draws). print() and summary() start with the model as it was named, then its class: model = "cross" (rtirt_cross), model = "cross", quantile = 0.9 (rtirt_cross); so do the messages about missing covariates or times. itemtype lists only "2pl" and "1pl" and names the Rasch model; a number as the third argument (gibbs(dat, "cross", 1000)) says that it is itemtype and to name n_iter =.
• gibbs(data, prior = "map") (the prior = of ecm(), matched by R to priors =) says that prior = "ml" or "map" is a setting of ecm() and that gibbs() samples the posterior under rtirt_priors(), instead of only "priors must come from rtirt_priors()". The z column of estimates() of an ecm() fit, and z and P(>|z|) of summary(), are left empty, like sig, for discriminations, time intensities and variances (a Wald test of a = 0 or sigma2t = 0 is no hypothesis).
• Documentation and messages no longer point to the deprecated interface: ?check_ecm (an ecm() fit, and a hand-written EM, instead of rtirt_cross() and block_model()), ?check_sampler (an ecm() or gibbs() fit; a note that it checks the sampler code, not the convergence of an analysis), ?sim_data, ?score_test, ?rtirt_priors, ?person_fit, ?preknowledge_test, and the messages of check_sampler(), check_ecm(), score_test() (impact: ecm(input_data(..., cov = z), model = "latreg")), person_fit() and preknowledge_test(). The intro tutorial checks ecm() fits in appendix B, and its FAQ no longer sends missing data to the deprecated blocks; it covers long-format data and categorical covariates.
• input_data() refuses data in long format (one row per person and item, as exported by PsychoPy or E-Prime), recognized by a single item or by ids that repeat on most rows, with the call to reshape them (the data look like long format (1 item, 300 distinct ids on 3000 rows): ... wide <- tidyr::pivot_wider(data, id_cols = id, names_from = item, values_from = c(correct, rt)) ...); such data became "3000 persons x 1 item" without a message. A few repeated ids give a warning. A non-numeric covariate gets how to dummy-code it (model.matrix(~ grade, data)[, -1], k - 1 columns for k categories, then centre them).
• The sig column of estimates() (and precis()) is filled only for parameters where 0 is a meaningful value (b, beta, rho, cor_ability_speed, b_ability); it is NA for discriminations, variances (sigma2t, var_speed) and time intensities. The intervals of the positive parameters are on the log scale and never contain 0, so every a was starred, even a = 0.003 with an interval from 1e-105 to 6e99.
• print() and summary() of an ecm() fit warn when the estimates suggest too little data for maximum likelihood: |b| > 6, a < 0.1 or a > 10, or a standard error above 5 (or none) for a, b or rho (Warning: extreme estimates or standard errors (b[item02] = -34.1, SE 4105; ...): the sample may be too small for maximum likelihood; try ecm(..., prior = "map") or gibbs(); 40 persons x 10 items). The thresholds are in ?ecm. The header of summary() of an ecm() fit shows the method line again (it read a field that did not exist).
• print() of a gibbs() fit warns with the criteria of convergence() (ESS > 400 and R-hat < 1.05 for every parameter that varies: Warning: only 98.0% of 51 parameters have ESS > 400 and R-hat < 1.05 (max R-hat 1.02, min ESS 318); run longer (n_iter) or see convergence()), not only when R-hat exceeds 1.05; with one chain it says that R-hat cannot be computed. convergence() returns an rtirt_convergence object that its print() method shows (it printed even when assigned, and print() of the result showed the raw list); it leaves out fixed parameters and works on several quantile levels (convergence(fits) for an rtirt_qset, one block per level), and print() of an rtirt_qset names the level of each warning (Warning (q = 0.25): ...).
• ecm(data, model = ...) and gibbs(data, model = ...) refuse arguments they do not have, before anything is fitted (they were warned about after the run, or failed with "formal argument matched by multiple actual arguments"). A setting of ecm_control() / gibbs_control() written directly gets the call to use (nodes is a setting of ecm_control(), not an argument of ecm(): use ecm(dat, model = "cross", control = ecm_control(nodes = 61))), an argument of the other engine says so, and a misspelled name gets the closest argument (unknown argument niter; did you mean n_iter?). The deprecated ecm(model_object, nodes = ) keeps its arguments.
• gibbs(data, model = "null" | "latreg" | "irt", quantile = ) is an error, as in ecm(): only "latent" and "cross" have quantile versions (the quantile was ignored and the fit printed as "(quantile 0.9)"). Also for equations() and algorithm() of data.

birtRcpp 0.1.0

• score_test(fit, add = Y ~ delta[item] * z): Lagrange multiplier (score) tests for declared models, from the ecm() (or mml()) fit of the null model only (Rao, 1948; for IRT Glas, 1998, 1999; Glas & van der Linden, 2010). add names an observed part and a term added to its linear predictor (bernoulli_logit()) or mean (normal()); new parameters are written with their dimension (delta[item]) or as one value, and are 0 under the null. Covariates come from the data of block_model() or from data = list(z = ). The score is the expected complete-data score of the E-step (Fisher's identity) and the information is the observed information by Louis' formula, with the null parameters projected out; one joint test and one test per item (Holm, with the one-step estimate). Directions that the null parameters absorb (e.g. a shift of every item with a free latent mean) are found from the casewise scores and left out of the joint test. Checks (2PL, 500 persons, 8 items, 100 replications): LM and LR statistics correlate 0.9995 without DIF and 0.997-0.999 with DIF; size under the null 5.0% (joint, KS p = 0.62) and 4.7% per item (KS p = 0.38), also with impact and a latent mean (4.0%, 5.2%; df 7); power for a uniform DIF of 0.5 / 0.8 logits on one item 57% / 92%. It reproduces ci_test() of the built-in cross model (add = logT ~ delta[item] * Y; statistics within 4e-4) and the Merkle-Zeileis LM test of score_test() for b along a binary moderator (correlation 1.000).
• Louis' information for declared models with person-level parameters that are not integrated out (ecm(latent = "theta") with a fixed speed per person) is analytic: the covariance of the complete-data scores is block diagonal by person and is computed person by person. Standard errors equal the numerical derivative of the score to 3e-8 (relative); the fit of 200 persons x 6 items takes 0.3 s instead of 1.5 s, and its Hessian 0.01 s instead of 1.5 s. Estimates and standard errors of the other fits are unchanged (identical).
• algorithm() of declared Gibbs (gibbs(), algorithm(model, steps = )) and ECM fits is written as formulas, as for vi(): the model p, then each step as its full conditional (Gibbs: aᵢ | · ∼ N⁺(mᵢ/pᵢ, 1/pᵢ) with p and m written out, Pólya-Gamma and half-t auxiliary steps, MALA, shifts, cuts) or, for ECM, the objective, the E-step (quadrature, exact integration, E[omega]) and each conditional mode with its expectations; mml() fits show one maximization. Unicode or format = "latex". Expressions are printed with the parentheses they need (x − (a − b); the residuals of the CAVI variance updates were printed without them).
• vi(quadrature = "ba81"): common nodes with the Bock-Aitkin collapse of the logistic parts, as in ecm(). When a latent variable is integrated exactly and its parts are affine in the latent variables (speed in the RT models), vi() now uses the steps of ecm(): the exact step from a quadratic per person instead of at every node, and each person's mean and covariance through 2L points in the updates. Same fixed point as before (1e-13), the ELBO never decreases (check_vi(quadrature = )). 1000 persons x 20 items (cross RT model): 203 s before, 56 s now, 21 s with quadrature = "ba81" (same 771 sweeps and ELBO; ba81 within 6e-4; machine load 1.7-2.9 for all three).
• First CRAN release. ECM is the main estimator: ecm() (Meng & Rubin, 1993; was em(), which was ECM all along), with check_ecm() and gibbs(init = "ecm"); direct maximization by Newton or quasi-Newton steps is mml() (was em(method = "newton")). Declared models (block_model()) are fitted with the same verbs as the built-in models: ecm(), gibbs() and vi() (variational Bayes); ecm() takes nodes and prior = "ml" | "map" for both; gibbs() takes one set of argument names (n_iter, n_chain, n_burnin, n_thin; Julia's cov2one only in sampleGibbs()). Samplers written by the user are checked with check_sampler(step = ). draw_shift(lower = ) draws a shift under truncated priors exactly (replaces draw_location_shift()). DIF trees moved to birt (birt::dif_tree()).
• The built-in models port the Julia package ExtendedRtIrtModeling.jl with corrected full conditionals.
• Faster ECM for declared models (same maxima): starting values from the data for intercept-like parameters (e.g. lambda at the mean log time, b at the logit of the proportion correct; gibbs(), ecm(), mml(); vi() keeps the prior means so that the ELBO rises from the first sweep), one Polya-Gamma expectation per M-step shared by the parameters of a logistic part, and quadrature cells of negligible weight skipped. On 1000 persons x 20 items (cross RT model, no starting values) ecm() went from 44 s (182 iterations) to 5-7 s, below mml() (Newton, 16-20 s). The default tol of ecm() for declared models is now 1e-10, as for the built-in models. Printing and fit_indices() of ECM fits say "ECM"; prior = "map" reports no AIC / BIC; convergence(), compare_para() and plot(type = "items") work on ECM fits.
• summary() returns the tables of a fit as one S3 object ($fit, $items, $se, $parameters with the prior of each parameter, ...), printed in the layout of tppcm's irt_pars() and blavaan; as.data.frame() gives the parameter table.
• algorithm() lists the computation steps of a fit: every full conditional of the Gibbs sampler in the order of one iteration (also the partially collapsed version), or the E- and M-steps, SQUAREM and standard errors of ecm(), with the prior constants substituted and the building block (draw_*()) of each step; Unicode or format = "latex".
• equations() writes the model as equations (Unicode by default, format = "latex" for an aligned block), with the parameter symbols, the items they apply to and the priors.
• Models: mlirt(), rtirt_null(), rtirt_latreg(), rtirt_latent() and rtirt_cross(), each with quantile versions (quantile =).
• Data input with checks: input_data(resp, time, cov, id, data).
• Sampling with gibbs(); gibbs(init = "ecm") starts the chains at the EM fit.
• gibbs(collapse = TRUE): a partially collapsed sampler (van Dyk & Park, 2008). Speed is integrated out of the ability, RT item and structural regression updates (ability out of the mlirt() regression), and the 2PL item parameters (a, b) are drawn jointly with the Polya-Gamma variables integrated out (Metropolis-Hastings with a Fisher-scoring proposal in (a, -a b), followed by the usual Polya-Gamma steps). Same posterior; on TIMSS the smallest effective sample size grows 1.9-7 times (1.3-5 times per second); no gain for 1PL models or rtirt_latreg(speed_var = "fixed"). The default sampler is unchanged.
• rtirt_priors(): the prior distributions of gibbs() (and of the item parameters for ecm(prior = "map")) can be changed; the Prior column of summary() and equations() show the priors used. The defaults are common data-independent priors (a ~ N+(1, 1), b ~ N(0, 3^2), lambda ~ N(0, 10^2), half-t(3, 1) on the residual and speed SDs, N(0, 1) for regression coefficients, rho and c) and every full conditional stays conjugate: the half-t priors are sampled through the inverse-gamma auxiliary variable of Huang & Wand (2013). rtirt_priors_julia() gives the priors of ExtendedRtIrtModeling.jl (b truncated to [-4, 4], lambda centred at the mean and SD of the log times, IG(0.001, 0.001) variances) and reproduces the earlier draws exactly.
• score_test(): score-based tests of parameter invariance (Merkle & Zeileis, 2013) from one ecm() fit, with casewise scores from estfun() (and bread()) for 'strucchange' and 'partykit': DIF along gender, item position or the total response time (valid for the accuracy parameters in the joint model, with decorrelation within the tested block), per item with Holm adjustment, a message when the moderator is related to ability (impact), and cluster-level moderators.
• Building blocks (?blocks, draw_*()): the full conditionals of the samplers callable from R, so new samplers can be composed in R with every draw in C++; inst/examples/ex07_blocks.R rewrites the rtirt_cross() sampler with them and reproduces gibbs() to rounding. draw_pg() draws Pólya–Gamma variables for any matrix of logits, so models from the literature can be composed: ex08_testlet.R (Bradlow, Wainer & Wang, 1999), ex09_mirt_rt.R (compensatory multidimensional IRT with response times; Man, Harring, Jiao & Zhan, 2019) and ex10_conditional_dependence.R (Bolsinova, De Boeck & Tijmstra, 2017).
• check_sampler(): simulation checks that a sampler draws from the right posterior, for any model and setting of gibbs(): the joint distribution test of Geweke (2004; successive-conditional simulator compared with the prior through prior quantiles, ESS-based z-tests, Holm adjustment) and simulation-based calibration (Cook, Gelman & Rubin, 2006; Talts et al., 2018). sampler_priors checks the check (a wrong prior must be detected). scripts/60_check_samplers.R runs the Geweke test for all samplers.
• person_fit(): the RT person-fit statistic of Marianti et al. (2014), evaluated in every iteration of gibbs() (posterior mean and posterior probability of exceeding the chi-square 95% quantile; normal log times).
• ci_test(): Lagrange multiplier test of conditional independence between responses and response times (van der Linden & Glas, 2010; shift of the log times of correct responses), from an ecm() fit with the observed information of the extended model (correct size; the outer-product version rejected too often).
• Helpers for composing samplers (?compose): calc_normal_part() (precision and linear term of a parameter entering normal pseudo-observations: Pólya–Gamma, probit or log times), draw_mvn() (batched multivariate normal draws in precision form, with truncation), draw_items_mh() (item parameters with the latent responses integrated out, logit or probit), draw_invwishart(), draw_mean_mvn(), draw_cor(), calc_person_fit(). The examples ex09-ex12 are rewritten with them (same results, 1.5-2.5 times faster).
• Notation: persons j, items i (y_ij is the response of person j to item i) in equations(), algorithm(), the documentation and the tutorial, as in birt.
• Default item names are zero-padded (item01, ...).
• preknowledge_test(): statistics of Sinharay (2020) for item preknowledge (RT person fit chi_pf, signed likelihood ratio of speed Lambda_s and of ability L_s on a compromised item set). inst/examples/ex13_preknowledge.R reproduces Table 2 of Sinharay (2020) on the public licensure data, and ex12_variable_speed.R the variable-speed model of Fox & Marianti (2016) on the Amsterdam Chess data.
• estimates() / convergence() with diagnostics = "rank": rank-normalized split R-hat and bulk/tail ESS (Vehtari et al., 2021; package posterior).
• draw_probit_z(): latent normal responses of Albert & Chib (1993) for probit models. inst/examples/ex11_lnirt.R composes the LNIRT model (probit, time discrimination, item-level multivariate normal) and reproduces Fox & Marianti (2017, Table 2) on the public credential data, with the latent responses integrated out for (a, b) (Fisher-scoring Metropolis-Hastings) against the slow mixing of easy items; about 17 times the minimum ESS per second of LNIRT.
• ecm(): maximum marginal likelihood for the normal models (Bock-Aitkin EM with Polya-Gamma item steps, SQUAREM, observed-information standard errors), with logLik(), vcov(), AIC() and BIC(). Without a latent regression all persons share the quadrature nodes and the E- and M-steps collapse to matrix products on the Bock-Aitkin artificial data (the same iterates; 5000 x 60: fit 10.5 -> 1.0 s, standard errors 143 -> 20 s); person-specific nodes use C++ kernels.
• Internal clean-up (results unchanged up to rounding): one formula for the normal integral over speed (ecm(), marginal_loglik(), calc_normal_given_nodes()), one 2PL log-likelihood at the nodes (calc_loglik_nodes(), with softplus: dic() and marginal_loglik() no longer overflow for very large a(theta - b)), Gauss-Hermite nodes only from calc_nodes(), pg_mean() in C++, one Fisher-scoring Metropolis-Hastings step behind draw_items_joint(), draw_items_mh() and gibbs(collapse = TRUE). runSimulation() calls run_simulation() (new argument speed_var). ecm(quadrature = "ba81", analytic = FALSE) warns with the gap to adaptive quadrature and advises analytic = TRUE.
• EM components (?em_components): calc_nodes(), calc_loglik_nodes(), calc_normal_given_nodes(), calc_node_weights(), calc_artificial_data(), calc_item_sums(), update_items_logistic(), em_squarem() and calc_hessian(); ecm() is built from them and they compose new EM algorithms.
• Three estimation engines for declared models, sharing the model and the step list (Gibbs draws from each full conditional, ECM moves to its mode with the latent parts averaged, CAVI sets q to it with expected coefficients):
  • ecm(): ECM. The person parameters are integrated by adaptive Gauss-Hermite quadrature (frozen centres, recentred between rounds; all nodes stacked); a person parameter that enters only normal parts linearly (e.g. speed) is integrated exactly given the others; Polya-Gamma variables are replaced by E[omega]; step_mala() becomes gradient ascent (GEM); standard errors from Louis' formula (expected complete-data Hessian plus the covariance of the complete-data score, from the first and second derivatives of the declared expressions; run_em(hessian =)), equal to the numerical derivative of the Fisher-identity score to 1e-7 and about 10 times faster. quadrature = "ba81": Bock & Aitkin (1981): the same nodes for every person; data parts that depend on the persons only through the grid variables are collapsed to artificial data per node and item (E-step by matrix products, M-step on nodes x items); still ECM. The normal parts of a step are summed in C++ (.ml_part). When every part of an exactly integrated variable is linear in all latent variables with a variance free of them (speed in the cross and hierarchical RT models), the exact step is computed from sums over items once per person (p, and h linear in the grid variables) instead of at every node, and the M-step of those parts uses each person's posterior mean and covariance of the latent variables through 2L points (exact, since the parts are quadratic in them); otherwise the general path is used. The declared cross model reproduces ecm(rtirt_cross()) (log-likelihood equal, estimates within 1e-4, same standard errors) in 0.25 s for 300 x 8 and 2.1 s (quadrature = "ba81": 1.7 s) for 1000 x 20 (15.4 s before, with numerical standard errors).
  • vi(): CAVI with q(person) free-form on the same grid, truncated-normal q for normal priors, inverse-gamma q for variances (half-t through the Huang-Wand auxiliary) and q(omega | theta) = PG(1, c) (Durante & Rigon, 2019); exact mean-field expectations from the multilinear form of the expressions; the ELBO never decreases. Means close to the posterior; SDs of item parameters too small, as for any mean-field approximation.
  • ecm(method = "newton") or "quasi-newton": maximize the same objective directly with nlminb() within the parameter bounds (gradient by Fisher's identity, Hessian by Louis' formula; Glas & van der Linden, 2010); no closed-form steps needed. Same maximum as ECM; ECM is faster on large models (1000 x 20 RT model: 2.2 s vs 12.9 s), Newton is faster when a parameter has no closed-form step (2PL with exp(la): 2.9 s vs 0.4 s).
  • Fixes from an independent review of the blocks layer: a truncated prior whose mean or sd uses a parameter is refused (its normalizing constant would depend on the parameter); observed data given as person vectors work in ecm() and vi() (masks stacked; the C++ kernels check every length); with prior = "ml" the priors of fixed parameters are left out of every part (score, M-step, Hessian) and a parameter that appears only in such priors is refused; parameters at a bound of their prior get no standard error, with a warning; quadrature = "ba81" compares its objective with adaptive nodes at the estimates (fit$adaptive_gap) and warns when the common nodes are too coarse. check_vi() also checks the scales of the variance factors of q.
  • equations() of a block_model() and algorithm() of a vi() fit: Unicode (default) or format = "latex" (an aligned block, as equations() of the built-in models); indices j (person) and i (item) added. The CAVI algorithm is written the mean-field way: the model p as a product of its parts, the factorized q, the update rule log q*(x) = E_-x[log p] + const (Bishop, 2006, Eq. 10.9), and the form of each factor of q with its expected coefficients.
  • check_ecm() and check_vi(): the counterparts of the Geweke check (every step raises the objective from many starts; Fisher identity; the fit is a maximum). check_ecm() also checks hand-written EMs and the EM of ecm() for the built-in normal models (check_ecm(rtirt_cross(data))).
  • step_gibbs(..., cut = "logT"): a cut (Plummer, 2015).
  • Definitions as in rethinking::ulam(): eta <- a * (theta - b).
  • Expressions are compiled to a multilinear form (sum of monomials) and evaluated in C++ without expanding person and item vectors to matrices.
• Declared models for the building blocks, in the style of PyMC: block_model() declares the model once (theta[person] ~ normal(0, 1), Y ~ bernoulli_logit(a * (theta - b)), variances with half_t() or inv_gamma()), and gibbs() runs stacked steps that only name parameters: step_gibbs() (normal full conditionals summed over every part of the model that uses the parameter; conjugate variances), step_mala() (exact likelihood; log scale with the Jacobian for bounded parameters; adapted step sizes) and step_shift() (generalized Gibbs translations, e.g. step_shift(rho = -1, zeta = theta)). The Polya-Gamma variables are hidden and redrawn when a step needs them (once per iteration and after MALA moved a parameter of the logistic part), so stale-omega errors cannot occur; steps that do not cover every parameter, mix dimensions, or ask Gibbs for a nonlinear parameter are refused with a hint. algorithm() and equations() describe the composed sampler and the model; check_sampler(model, steps = ) runs the Geweke test with parameters and data simulated from the declaration. Checks: the declared cross model reproduces gibbs(rtirt_cross()) (41 parameters, max 2.0 Monte Carlo SEs); 58 random step stacks for a 2PL all match a long reference (hand-written steps: 23 of 240 random stacks gave a wrong posterior).
• A workflow layer for the building blocks, after a usability test with a graduate-student reader (samplers and EM written as "a state list + one step function"):
  • run_sampler(data, init, step, monitor): chains, burn-in, thinning, storage, running person means and the summary table (mean, SD, interval, R-hat, ESS); n_iter, n_burnin, n_thin and n_chain as in gibbs().
  • Normal parts: irt_part() and calc_normal_part() return a normal part (precision and linear term); normal_prior(mean, sd) adds a prior with +, and draw(part, lower, upper) draws (truncated when bounds are given). No hand conversion of priors to the precision form.
  • check_sampler(step = , prior = , simulate = , monitor = ): Geweke's joint distribution test for samplers written by the user (rejects a wrong prior).
  • draw_shift(): exact shift along directions identified only by the priors (e.g. the intercept of a latent regression with the difficulties).
  • run_em(data, init, step, loglik, positive): packing, log scale for positive parameters, SQUAREM, convergence and standard errors from the Hessian of loglik(); estep_irt() and mstep_items() are the E- and M-steps of a 2PL (person-specific prior means, missing responses, ML or MAP). A 2PL written this way equals mirt's estimates and standard errors.
  • Readable errors: dimension checks in draw_omega() (argument order), irt_part(), draw_items_pg(), draw_items_joint(), calc_normal_part(); draw_gauss() refuses non-positive precisions; draw_mala() checks logpost() and takes joint = TRUE for one target (a scalar lp with several elements used to be applied elementwise).
  • Missing responses (NA) are skipped by irt_part(), draw_items_pg() and draw_items_joint().
  • reg_part(X, y, w): the normal part of regression coefficients (a precision matrix), draw(reg_part(X, theta) + normal_prior(0, 10)).
  • b_mean in draw_items_pg(), draw_items_joint() and mstep_items(): a prior mean of b per item (LLTM with item covariates). mstep_items() takes the item priors as priors = rtirt_priors(), as the samplers do.
  • print() of a block_fit warns when R-hat > 1.05 or ESS < 100 and names the usual cures; run_em() warns when the objective decreases, when s$.loglik differs from loglik(), when the estimates are not a maximum of loglik() (a wrong M-step converges to the wrong values; also when the Hessian is not negative definite), and when step() changes no parameter (no M-step); SQUAREM stops on non-finite steps. draw_sigma2() and draw_lambda() stop with a message when the log times contain NA (they returned NaN before); use calc_normal_part() with draw() and draw_var() for missing times.
  • draw_items_joint() and draw_items_mh() return list(a, b) with the attribute accepted, as draw_mala() and draw_cor().
  • Help pages are written with roxygen markdown (code, lists, tables and % render correctly).
  • The blocks tutorial: a JAGS-to-blocks table, the Pólya–Gamma idea before the first example, check_sampler(step = ) (and what it cannot detect), identification and draw_shift() (with the shift-invariant quantities to report), an LLTM with b_mean and reg_part(), and a table of common mistakes, and writing an EM with run_em() (with the warnings for a wrong M-step).
• draw_mala() and mcmc_adapt() (?mala): vectorized Metropolis-adjusted Langevin steps with Robbins-Monro step adaptation, for blocks without a conjugate full conditional; draw_cor(method = "mala") (about 2.5-3 times the ESS of the random walk; the default stays "rw").
• marginal_loglik() and marginal DIC / WAIC / PSIS-LOO in fit_indices().
• The generics shared with birt dispatch to both packages whichever was attached last.
• Two tutorials (inst/tutorial/): birtRcpp_intro (the ready-made models, for beginners) and birtRcpp_blocks (samplers and EM from the building blocks), each code chunk explained line by line.
• convergence(): the share of converged parameters was NaN without tail ESS (classic diagnostics); the R-hat threshold is 1.05, as in print().
• Accessors: summary(), estimates(), scores(), reliability(), fit_indices() and convergence(). summary() lists every parameter with its prior (as blavaan). The person variables are named ability and speed in scores, reliabilities and parameter names (beta_ability, b_ability, cor_ability_speed).
• Julia-style aliases.