Changes in version 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.

                        Changes in version 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.

                        Changes in version 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.