Package 'birtRcpp'

Description

Built-in models of response accuracy and response times (ports of ExtendedRtIrtModeling.jl), named by model and fitted with ecm() (the main estimator) or gibbs(); algorithm() and equations() show the steps and the model, check_ecm() and check_sampler() check them:

dat <- input_data(resp = items, time = times, cov = covs, id = "IDSTUD", data = raw)
fit_ecm <- ecm(dat, model = "cross")
fit <- gibbs(dat, model = "cross", init = fit_ecm, n_iter = 10000, n_chain = 3, seed = 1)
summary(fit)
estimates(fit)
scores(fit)

Details

To start, these are all you need (the other help pages are for diagnostics, simulation, the building blocks of the samplers, the Julia names and deprecated functions):

The tutorial (Traditional Chinese) explains every line: https://github.com/jiewenTsai/birtExamples/blob/main/birtRcpp/tutorial/birtRcpp_intro.md.

The declarative block_model(), vi() and mml() are deprecated (they still work). The accessors have the names of the birt package. Julia names are aliases (julia_aliases). The generics shared with birt (estimates, scores, reliability, fit_indices, equations, convergence, score_test, algorithm) dispatch to the methods of both packages whichever was attached last: when both are loaded, each package's methods are registered in the other's generics.

Author(s)

Maintainer: Jie-Wen Tsai [email protected]

See Also

Useful links:

• https://github.com/jiewenTsai/birtRcpp

• https://github.com/jiewenTsai/birtExamples

• Report bugs at https://github.com/jiewenTsai/birtRcpp/issues

Description

Lists the steps of the algorithm behind a fit: for gibbs() every full conditional distribution in the order of one iteration (with collapse = TRUE the partially collapsed version), for ecm() the E- and M-steps. The prior constants of the fit are substituted into the formulas, and each step names the building block that draws it (blocks), so a step can be run, replaced or checked from R. Works before estimation too: give the data and the model as in ecm() and gibbs().

Usage

algorithm(object, ...)

## S3 method for class 'rtirt_qset'
algorithm(
  object,
  method = NULL,
  collapse = NULL,
  format = c("unicode", "latex"),
  detail = TRUE,
  ...
)

## S3 method for class 'rtirt'
algorithm(
  object,
  method = NULL,
  collapse = NULL,
  format = c("unicode", "latex"),
  detail = TRUE,
  ...
)

## S3 method for class 'rtirt_algorithm'
print(x, ...)

## S3 method for class 'input_data'
algorithm(
  object,
  model = NULL,
  itemtype = "2pl",
  speed_var = c("free", "fixed"),
  quantile = NULL,
  cov = NULL,
  intercept = FALSE,
  method = NULL,
  collapse = NULL,
  format = c("unicode", "latex"),
  detail = TRUE,
  ...
)

Arguments

Value

An rtirt_algorithm object (printed); as.character() gives the lines.

See Also

equations() for the model, blocks for the steps as R functions.

Examples

cond <- set_cond(n_subj = 100, n_item = 5)
dat <- sim_data(cond, sim_para(cond, "cross"), "cross")
algorithm(dat, model = "cross") # Gibbs full conditionals
algorithm(dat, model = "cross", collapse = TRUE) # partially collapsed Gibbs
algorithm(dat, model = "cross", method = "ecm") # ECM steps
algorithm(dat, model = "cross", format = "latex")
algorithm(dat, model = "cross", detail = FALSE) # the steps only

Description

Fits a block_model() by coordinate ascent variational inference with the same steps as gibbs() and ecm(). The three engines use the same full conditionals:

Usage

## S3 method for class 'block_model_vi'
algorithm(object, format = c("unicode", "latex"), ...)

## S3 method for class 'block_model'
vi(
  model,
  steps = NULL,
  latent = NULL,
  nodes = NULL,
  init = NULL,
  cycles = 6,
  max_iter = 1000,
  tol = 1e-09,
  verbose = FALSE,
  quadrature = c("adaptive", "ba81"),
  ...
)

## S3 method for class 'block_model_vi'
print(x, ...)

Arguments

Details

Expectations of products such as $E_q[(a(\theta - b))^2]$ are exact: the expressions are multilinear, and under the factorized q every product is a product of first and second moments. The objective is the evidence lower bound (ELBO); every update raises it. A latent variable integrated exactly (normal and linear, e.g. speed) uses the steps of ecm() when its parts are affine in all latent variables: q(u | grid) from a quadratic per person, and in the updates each person's mean and covariance of the latent variables through 2L points (exact for these quadratic parts; Julier & Uhlmann, 2004).

Known behavior of mean-field VB, also here: posterior means are close to the posterior, but posterior SDs of the item parameters are too small (the correlations between parameters are ignored); use gibbs() for intervals. step_mala() has no CAVI update (the parameter needs a conditionally conjugate full conditional) and step_shift() is skipped.

Value

A list of class block_model_vi: summary (posterior means and SDs of q), latent (posterior means and SDs of the person parameters), q, elbo (trace, one value per sweep), elbo_drop (the largest relative drop between sweeps of one round: 0 when every update is right), iterations, converged, time, model, steps; algorithm(fit) lists the updates.

References

Blei, D. M., Kucukelbir, A., & McAuliffe, J. D. (2017). Variational inference: A review for statisticians. Journal of the American Statistical Association, 112, 859–877.

Bock, R. D., & Aitkin, M. (1981). Marginal maximum likelihood estimation of item parameters: Application of an EM algorithm. Psychometrika, 46(4), 443–459. doi:10.1007/BF02293801

Durante, D., & Rigon, T. (2019). Conditionally conjugate mean-field variational Bayes for logistic models. Statistical Science, 34, 472–485.

Julier, S. J., & Uhlmann, J. K. (2004). Unscented filtering and nonlinear estimation. Proceedings of the IEEE, 92(3), 401–422. doi:10.1109/JPROC.2003.823141

Examples

set.seed(1); N <- 300; K <- 6
th <- rnorm(N); a0 <- runif(K, 0.8, 2); b0 <- seq(-1.5, 1.5, length.out = K)
Y <- matrix(rbinom(N * K, 1, plogis(sweep(outer(th, a0), 2, a0 * b0))), N)
m <- block_model(Y = Y,
  theta[person] ~ normal(0, 1),
  a[item] ~ normal(1, 1, lower = 0),
  b[item] ~ normal(0, 3),
  Y ~ bernoulli_logit(a * (theta - b)))
fit <- vi(m, list(step_gibbs(theta, a, b)))
fit
algorithm(fit)

Description

anova(fit0, fit1, ...) compares ecm() fits of the same responses and response times, as anova() of lm(): the number of parameters, log-likelihood, AIC and BIC of each fit, and the likelihood ratio test of each fit against the one in the row above (⁠Chisq = 2 * (logLik - logLik above)⁠, Df the difference in the number of parameters). The rows are ordered by the number of parameters. The test assumes that each fit is nested in the next, e.g. "null" in "latreg" (the regression coefficients 0) or in "cross" (equal rho), itemtype = "1pl" in "2pl", speed_var = "fixed" in "free"; AIC and BIC also compare fits that are not nested. A fit with more parameters but a lower log-likelihood gives a warning (not nested, or not converged).

Usage

## S3 method for class 'rtirt'
anova(object, ...)

## S3 method for class 'rtirt_anova'
print(x, digits = 3, signif.stars = getOption("show.signif.stars", TRUE), ...)

Arguments

Details

For gibbs() fits the table has the criteria of fit_indices() (DIC, WAIC, LOOIC with its standard error), smaller is better, and no test. Fits with prior = "map" (no likelihood), quantile fits (an ALD working likelihood) and a mix of ECM and Gibbs fits are refused.

Value

A data frame of class rtirt_anova (and anova), printed with a heading that names the fits.

Examples

cond <- set_cond(n_subj = 300, n_item = 8, n_feat = 1)
dat <- sim_data(cond, sim_para(cond, "latreg"), "latreg")
fit0 <- ecm(dat, model = "null")
fit1 <- ecm(dat, model = "latreg")
anova(fit0, fit1)

Description

Deprecated: block_model() and its steps still work (with a warning once per session) but get no new features. Use the built-in models, ecm(data, model = ...) and gibbs(data, model = ...).

Usage

block_model(...)

## S3 method for class 'block_model'
print(x, ...)

## S3 method for class 'block_model'
equations(object, format = c("unicode", "latex"), ...)

Arguments

Details

block_model() declares a model once, the way PyMC does: the data, and one line name ~ distribution(...) per parameter and per observed matrix. Samplers are then stacked as steps that only name the parameters they update (step_gibbs(), step_mala(), step_shift(), run by gibbs()). Each step derives what it needs from the model, so

• the steps always target the same posterior (a prior cannot be forgotten in one step);

• the Pólya–Gamma variables $\omega$ of a logistic part are hidden: they are drawn whenever a Gibbs step needs them and are fresh (drawn once per iteration and again after a MALA step has moved a parameter of that part);

• positive parameters are moved on the log scale by MALA, with the Jacobian added;

• dimensions are declared once and checked.

Data: named arguments. Matrices are persons x items (all the same size); vectors have one value per person (length N) or per item (length K); missing responses or log times are NA.

Parameters: name[person] ~ ..., name[item] ~ ... or name ~ ... (a single value).

Distributions

• normal(mean, sd, lower = -Inf, upper = Inf): mean and sd may use parameters and data, e.g. zeta[person] ~ normal(beta * theta, sqrt(tau2)). For an observed matrix the mean is written with the person and item parameters as they are; they are expanded to the matrix.

• bernoulli_logit(eta): logistic responses, e.g. Y ~ bernoulli_logit(a * (theta - b)).

• half_t(df, scale) and inv_gamma(shape, scale): priors of a variance v that is used as sd = sqrt(v); half_t() is the prior of the SD $\sqrt v$ (drawn through the inverse-gamma auxiliary variable of Huang & Wand, 2013).

Expressions may use ⁠+ - * / ^⁠, exp(), log() and sqrt(). A line name <- expression defines a name for later lines, as in rethinking::ulam(), e.g. eta <- a * (theta - b) and then Y ~ bernoulli_logit(eta); definitions are substituted (they are not parameters).

Value

An object of class block_model; equations() prints the model.

See Also

gibbs() and the steps; run_sampler() for samplers written by hand.

Examples

set.seed(1); N <- 300; K <- 6
th <- rnorm(N); a0 <- runif(K, 0.8, 2); b0 <- seq(-1.5, 1.5, length.out = K)
Y <- matrix(rbinom(N * K, 1, plogis(sweep(outer(th, a0), 2, a0 * b0))), N)
m <- block_model(Y = Y,
  theta[person] ~ normal(0, 1),
  a[item] ~ normal(1, 1, lower = 0),
  b[item] ~ normal(0, 3),
  Y ~ bernoulli_logit(a * (theta - b)))
m
fit <- gibbs(m, list(step_mala(theta), step_gibbs(a, b)), n_iter = 600)
summary(fit)[1:4, ]
algorithm(fit)

Description

Steps name the parameters they update; gibbs() runs them in the order given, once per iteration. What a step does is derived from the block_model():

Usage

step_gibbs(..., cut = character())

step_mala(..., step = 0.1)

step_shift(...)

Arguments

Details

• step_gibbs(...): draws each parameter from its full conditional, one after the other. A parameter with a normal prior gets its normal full conditional: the precision and linear term are added up over every part of the model that uses it (normal parts directly; a logistic part through its Pólya–Gamma variables, which are drawn automatically when needed). A variance with a half_t() or inv_gamma() prior gets its conjugate update. The model must be linear in the parameter (otherwise use step_mala()). With cut = "logT" the parts of the data named there do not inform this step (a cut; Plummer, 2015): e.g. step_gibbs(theta, cut = "logT") draws the ability from the responses only, while the response-time parameters still use the ability. The chain then targets the cut distribution, not the posterior (so the Geweke check of check_sampler() does not apply), and ecm() and vi() refuse a cut.

• step_mala(..., step = 0.1): one MALA step (Roberts & Tweedie, 1996) per person or item, for the named parameters jointly (they must have the same dimension); the logistic parts use their exact likelihood (no Pólya–Gamma variables). A parameter with a lower bound is moved on the log scale with the Jacobian; the step sizes adapt during burn-in.

• step_shift(...): moves several parameters together along a direction, e.g. step_shift(theta = 1, b = 1) (theta + c and b + c) or step_shift(rho = -1, zeta = theta), with c drawn from its exact conditional (a generalized Gibbs step; Liu & Sabatti, 2000). The logistic parts must not change along the direction, and the normal parts must be linear in c. It speeds up directions that only the prior identifies.

Value

A step (a list) for gibbs().

Description

Each ⁠draw_*()⁠ function draws one block from its full conditional and is named after the C++ function of the samplers that it calls (irt_part() draws nothing: it returns the accuracy part of the conditional of ability). algorithm() names the function of every step of gibbs().

Usage

draw_gauss(prec, num)

draw_invgamma(shape, scale)

draw_tnorm(mean, sd, lower = -Inf, upper = Inf)

draw_omega(a, b, theta)

draw_probit_z(Y, eta)

draw_pg(eta)

irt_part(Y, omega, a, b)

draw_items_pg(
  Y,
  omega,
  theta,
  a,
  priors = rtirt_priors(),
  one_pl = FALSE,
  b_mean = NULL
)

draw_items_joint(Y, theta, a, b, priors = rtirt_priors(), b_mean = NULL)

draw_lambda(L, m, k1nu, W, mu, sd, lower = 0)

draw_sigma2(
  L,
  m,
  lambda,
  current,
  priors = rtirt_priors(),
  nu = NULL,
  q = NULL
)

draw_var(alpha, beta, current, prior = c(df = 3, scale = 1))

draw_nu(r, s, q)

draw_theta_collapsed(prec, num, W, y, rho, m0, g, tau)

draw_rt_items_collapsed(
  L,
  k1nu,
  W,
  theta,
  mz,
  tau,
  with_rho,
  mu,
  sd,
  rho_sd,
  lambda,
  rho,
  lower = 0
)

Arguments

Details

The full conditionals used by gibbs(), exposed so that new samplers (other blockings, cut or semi-modular schemes, nested chains) can be put together in R while every draw still runs in C++. All of them use R's random number generator, so set.seed() makes a composed sampler reproducible, and a sampler written with them reproduces gibbs() step by step. Notation as in the samplers: persons j = 1..N, items i = 1..K, Y the 0/1 responses, L the log response times, logit P(y = 1) = a (theta - b), log T = lambda - zeta - rho theta + e.

• draw_gauss(prec, num): independent draws from N(num / prec, 1 / prec), elementwise (any conditionally normal scalar: ability, speed, rho_i, a shift).

• draw_invgamma(shape, scale), draw_tnorm(mean, sd, lower, upper): inverse gamma and truncated normal draws, elementwise.

• draw_omega(a, b, theta): Pólya–Gamma variables PG(1, a_i (theta_j - b_i)), an N x K matrix.

• draw_pg(eta): Pólya–Gamma variables PG(1, eta) for any matrix of linear predictors (models whose logit is not a (theta - b), e.g. several abilities or covariates in the item response).

• draw_probit_z(Y, eta): for a probit response model, the latent normal responses of Albert & Chib (1993), N(eta, 1) truncated to the positive (y = 1) or negative (y = 0) half line; missing y (NA) are not truncated. Given them, the probit model is a normal linear model.

• irt_part(Y, omega, a, b): the accuracy part of the conditional of ability given omega, a normal part (precision prec and linear term num, see normal_parts): draw(irt_part(Y, omega, a, b) + normal_prior(0, 1)).

• draw_items_pg(Y, omega, theta, a, priors, one_pl, b_mean): b, then a, given omega (Gibbs); returns list(a, b).

• draw_items_joint(Y, theta, a, b, priors, b_mean): (a, b) jointly with omega integrated out (Metropolis–Hastings with a Fisher-scoring proposal in (a, -a b)); returns list(a, b) with the attribute accepted (one logical per item). It mixes much better than draw_items_pg() for very easy or very difficult items; draw omega anew after it.

The item priors are those of rtirt_priors() (a: truncated normal on a > 0, b: normal); b_mean gives b a prior mean per item (e.g. W %*% delta in an LLTM with item covariates W), in place of the mean of priors$b.

• draw_lambda(L, m, k1nu, W, mu, sd, lower): time intensities given the RT means m (N x K, e.g. -zeta - rho theta), the ALD shifts k1nu (zeros for normal times) and the precisions W (N x K); prior N(mu, sd^2) on (lower, Inf).

• draw_sigma2(L, m, lambda, current, priors, nu = NULL, q = NULL): residual variances (normal) or ALD scales (with the mixing variables nu and the quantile q), from their current values (a half-t prior draws its inverse-gamma auxiliary variable from them).

• draw_var(alpha, beta, current, prior): a variance with likelihood var^-alpha exp(-beta / var) (n residuals u ~ N(0, var): alpha = n / 2, beta = sum(u^2) / 2) and the prior c(df, scale) (half-t on its square root, two inverse-gamma steps) or c(shape, scale) (inverse gamma).

• draw_nu(r, s, q): ALD mixing variables for residuals r and scales s (matrices).

• draw_theta_collapsed(prec, num, W, y, rho, m0, g, tau) and draw_rt_items_collapsed(L, k1nu, W, theta, mz, tau, with_rho, mu, sd, rho_sd, lambda, rho): the collapsed steps of gibbs(collapse = TRUE) (speed integrated out).

Missing responses (NA in Y) are skipped by irt_part(), draw_items_pg() and draw_items_joint(). draw_lambda() and draw_sigma2() need complete log times (they stop on NA); with missing times use calc_normal_part() with draw() and draw_var(). run_sampler() runs a sampler written from these blocks (chains, burn-in, storage, summary) and check_sampler() checks it.

These are low-level: arguments are checked only for their dimensions, and the caller is responsible for a valid sampler (each block must be drawn from its conditional given the current values of everything else, or from a marginal followed by the integrated variable). inst/examples/ex07_blocks.R writes the rtirt_cross() sampler with them; ex08_testlet.R, ex09_mirt_rt.R and ex10_conditional_dependence.R compose models from the literature.

Value

Draws (vectors, matrices or lists, see Details).

Examples

set.seed(1)
draw_gauss(prec = c(1, 4), num = c(0, 2))          # N(0, 1) and N(0.5, 1/4)
theta <- rnorm(200); Y <- matrix(rbinom(600, 1, plogis(theta)), 200)
om <- draw_omega(a = rep(1, 3), b = rep(0, 3), theta = theta)
str(draw_items_pg(Y, om, theta, a = rep(1, 3), priors = rtirt_priors()))

Description

The EM counterpart of the Geweke check of check_sampler(). A Gibbs sampler is right when every step leaves the posterior unchanged; an EM algorithm is right when

Usage

check_ecm(
  model = NULL,
  steps = NULL,
  latent = NULL,
  nodes = NULL,
  prior = c("map", "ml"),
  n_start = 10,
  data = NULL,
  init = NULL,
  step = NULL,
  loglik = NULL,
  positive = character(),
  seed = NULL,
  quadrature = c("adaptive", "ba81")
)

## S3 method for class 'em_check'
print(x, ...)

Arguments

Details

1. every step never lowers the objective (from many starting values, not only from the path that EM happened to take);

2. the E-step agrees with the objective (Fisher's identity: the gradient of the log marginal likelihood equals the gradient of the expected complete-data log likelihood $Q(\cdot \mid x)$ at $x$); and

3. ECM stops at a maximum (a general optimizer started at the ECM estimate finds nothing better, and the gradient is zero there).

For an ecm() fit of a built-in model (a gibbs() fit or a model object of the deprecated constructors also works) the three are checked for the EM of ecm() with its own E- and M-step, from n_start starting values jittered around those of ecm(): one whole M-step (which updates all blocks) never lowers the objective, the analytic score (Fisher's identity) equals the numerical gradient of the objective, and the fit is a maximum. For a quantile model (one level, e.g. check_ecm(ecm(dat, model = "cross", quantile = 0.9))) the objective is the ALD working log-likelihood, the starts are jittered around the fit, and prior is not used. For an EM written by hand, give data, init, step and loglik as for run_em(); checks 1 (from n_start jittered starts) and 3 are done. For a block_model() (deprecated) all three are checked, with starting values drawn from the declared priors.

Value

A data frame with one row per check (check, worst, limit, passed), class em_check.

Examples

# the EM of ecm() for the cross-relation model
cond <- set_cond(n_subj = 200, n_item = 5)
dat <- sim_data(cond, sim_para(cond, "cross"), "cross")
fit <- ecm(dat, model = "cross", se = FALSE)
check_ecm(fit, prior = "ml", n_start = 3, seed = 1)
# an EM written by hand (the 2PL of run_em())
set.seed(1)
Y <- matrix(rbinom(300 * 5, 1, plogis(outer(rnorm(300), seq(-1, 1, length.out = 5), "-"))), 300)
step <- function(s, d) {
  e <- estep_irt(d$Y, s$a, s$b)                 # E-step
  s[c("a", "b")] <- mstep_items(e, s$a, s$b)    # M-step
  s$.loglik <- e$loglik
  s
}
loglik <- function(s, d) estep_irt(d$Y, s$a, s$b)$loglik
check_ecm(data = list(Y = Y), init = list(a = rep(1, 5), b = rep(0, 5)), step = step,
          loglik = loglik, positive = "a", n_start = 3, seed = 1)

Description

Tests whether the sampler of gibbs() for a model draws from the correct posterior, with two general methods that need only the prior and the likelihood:

Usage

check_sampler(
  model = NULL,
  method = c("geweke", "sbc"),
  n_persons = 20,
  n_iter = NULL,
  n_rep = 200,
  n_draws = 99,
  collapse = NULL,
  itemtype = NULL,
  priors = rtirt_priors(),
  sampler_priors = priors,
  probs = c(0.1, 0.25, 0.5, 0.75, 0.9),
  n_prior = 1e+05,
  seed = NULL,
  steps = NULL,
  monitor = NULL,
  step = NULL,
  prior = NULL,
  simulate = NULL,
  ...
)

## S3 method for class 'rtirt_check'
print(x, n = 10, ...)

Arguments

Details

• method = "geweke" (Geweke, 2004): the successive-conditional simulator. One sweep of the sampler given the data alternates with new person values and a new data set simulated given the current parameters. (Simulating new data from the old person values is also exact, but the persons are then held by their own data and the chain crawls: a false alarm in practice.) If the sampler is correct, the chain of the parameters has the prior as its distribution. For every parameter and each probability in probs, the share of the chain below the prior quantile is compared with the probability (z-test with the effective sample size of the chain); p-values are Holm-adjusted over all comparisons. Cheap: one sweep per iteration on a small data set.

• method = "sbc", simulation-based calibration (Cook, Gelman & Rubin, 2006; Talts et al., 2018; Modrák et al., 2025): n_rep times, draw the parameters from the prior, simulate data, run the sampler and record the rank of the true value among n_draws thinned posterior draws. The ranks are uniform for a correct sampler; a chi-square test over 10 bins per parameter (Holm-adjusted).

Both use the model only as a template: its type, quantile, settings, number of items and (the first n_persons rows of) its covariates. Data are simulated, so a small n_persons keeps the posterior wide and the check fast. The priors must be proper and must not depend on the data (rtirt_priors_julia() centres lambda on the data and cannot be used).

sampler_priors lets the sampler use other priors than the simulation: a check of the check, which must then report a failure.

Chain length matters for both. The successive-conditional chain must travel over the whole prior, and its parameters move by about one posterior width per iteration, so parameters with wide priors (lambda ~ N(0, 10^2)) need long chains; the print method warns when an effective sample size is small. In SBC the posterior draws of each fit must be close to independent after thinning: with short fits a slowly mixing parameter (such as rho) piles its ranks up at one end, which looks like an error of the sampler. Both effects disappear with longer chains, a real error does not. Samplers written by the user from the building blocks (run_sampler()'s step()) are checked with step, prior, simulate and monitor (Geweke test): prior() returns a state drawn from the prior (all parameters and the person parameters), simulate(s) the data list given the state, and step(s, data) is the sampler. Keep the model small (10-20 persons, a few items) and the priors proper and not too wide, so that the chain mixes. The test finds wrong full conditionals and forgotten or wrong prior terms, which summaries of a fit cannot show; its power is that of a test on a small model (an error that moves the posterior by a fraction of a posterior SD can pass), so passing is evidence, not proof. This checks the code of the sampler, not the convergence of your own analysis: for that use convergence() and plot(fit, "trace").

Value

A data frame of class rtirt_check with one row per comparison (Geweke: parameter, probability, share in the chain, effective sample size, z, p-value, Holm-adjusted p-value; SBC: parameter, chi-square statistic, p-value, Holm-adjusted p-value) and attributes method, passed (all Holm-adjusted p-values above 0.01), draws (Geweke chain) or ranks (SBC).

References

Cook, S. R., Gelman, A., & Rubin, D. B. (2006). Validation of software for Bayesian models using posterior quantiles. Journal of Computational and Graphical Statistics, 15, 675-692.

Geweke, J. (2004). Getting it right: Joint distribution tests of posterior simulators. Journal of the American Statistical Association, 99, 799-804.

Modrák, M., Moon, A. H., Kim, S., Bürkner, P., Huurre, N., Faltejsková, K., Gelman, A., & Vehtari, A. (2025). Simulation-based calibration checking for Bayesian computation: The choice of test quantities shapes sensitivity. Bayesian Analysis, 20(2).

Talts, S., Betancourt, M., Simpson, D., Vehtari, A., & Gelman, A. (2018). Validating Bayesian inference algorithms with simulation-based calibration. arXiv:1804.06788.

Examples

cond <- set_cond(n_subj = 100, n_item = 3)
dat <- sim_data(cond, sim_para(cond, "cross"), "cross")
fit <- ecm(dat, model = "cross", se = FALSE) # the model to check (its fit is not used)
narrow <- rtirt_priors(lambda = c(mean = 0, sd = 1)) # any proper prior; narrow mixes fast
check_sampler(fit, n_iter = 1e5, priors = narrow, seed = 1) # Geweke
check_sampler(fit, method = "sbc", n_rep = 100, priors = narrow, seed = 1)
# a wrong prior in the sampler is detected
wrong <- rtirt_priors(lambda = c(mean = 0, sd = 1), a = c(mean = 1, sd = 0.3))
check_sampler(fit, n_iter = 2e4, priors = narrow, sampler_priors = wrong, seed = 1)

Description

Deprecated with vi().

Usage

check_vi(
  model,
  steps = NULL,
  latent = NULL,
  nodes = NULL,
  n_start = 5,
  sweeps = 20,
  seed = NULL,
  quadrature = c("adaptive", "ba81")
)

Arguments

Details

The CAVI counterpart of check_ecm() and the Geweke check: coordinate ascent is right when

1. the ELBO never decreases from one sweep to the next (from n_start starting values drawn from the declared priors, on a fixed grid); and

2. the fit is a maximum of the ELBO: moving the mean of any q (normal factors) a little up or down lowers it.

Value

A data frame with one row per check (class em_check, printed like check_ecm()).

Description

The Lagrange multiplier (score) test of van der Linden & Glas (2010): under conditional independence the distribution of the log time on item i does not depend on whether the response is correct, given ability and speed. The alternative shifts the log times of correct responses, $\log t_{ij} = \lambda_i + \delta_i y_{ij} - \zeta_j [- \rho_i\theta_j] + e_{ij}$, and the test is of $\delta_i = 0$ from the maximum likelihood fit of the model without the shift (ecm()). The score of $\delta_i$ for person j is $y_{ij} E[\log t_{ij} - \lambda_i + \zeta_j (+ \rho_i\theta_j) | \text{data}_j] / \sigma_i^2$ (Fisher identity). The information is the observed information of the model extended by the shifts (central differences of the analytic score, as for the standard errors of ecm()) with the parameters of the model projected out, so the statistic is chi-square(1) per item, and chi-square(K) for all items together. (The outer product of the casewise scores is simpler but rejects too often: about .065 per item and .10 overall at the 5% level with 500 persons and 10 items.) A positive delta means that correct responses took longer than incorrect ones.

Usage

ci_test(object, ...)

Arguments

Value

A data frame of class rtirt_ci_test per item: item, delta (one-step estimate of the shift), statistic, p.value, p.holm; attribute overall (joint test of all items).

References

van der Linden, W. J., & Glas, C. A. W. (2010). Statistical tests of conditional independence between responses and/or response times on test items. Psychometrika, 75, 120-139.

Examples

cond <- set_cond(n_subj = 500, n_item = 8)
d <- sim_data(cond, sim_para(cond, "cross"), "cross")
ci_test(ecm(d, model = "null", se = FALSE))

Description

As comparePara() in Julia (for simulations with true_para): bias, RMSE and correlation of the posterior means (or of the ECM estimates) for one parameter block.

Usage

compare_para(object, par = "a")

Arguments

Description

General steps that come back in many samplers written from the blocks (blocks); with them the samplers of the examples ex08-ex12 (testlet, multidimensional, conditional dependence, LNIRT, variable speed) are a few lines per step. Persons j are the rows and items i the columns of the matrices.

Usage

calc_normal_part(x, z, w = 1, rest = 0, by = c("item", "person"))

draw_mvn(prec, num, lower = -Inf, upper = Inf, max_tries = 100)

draw_items_mh(
  Y,
  theta,
  a,
  b,
  prior_mean,
  prior_prec,
  link = c("logit", "probit"),
  form = c("difficulty", "intercept")
)

draw_invwishart(df, scale)

draw_mean_mvn(X, Sigma, mean0, cov0)

draw_cor(x1, x2, cor, steps = 5, step = 0.05, method = c("rw", "mala"))

calc_person_fit(L, fitted, sigma2)

Arguments

Details

• calc_normal_part(x, z, w, rest, by): the part of a full conditional that comes from normal (pseudo-)observations $z_{ij} = x_{ij}\beta + \text{rest}_{ij} + e_{ij}$, $e_{ij} \sim N(0, 1/w_{ij})$, for a parameter $\beta$ of each item (by = "item", sums over persons) or of each person (by = "person", sums over items): precision $\sum w x^2$ and linear term $\sum w x (z - \text{rest})$. Add the prior and draw with draw_gauss() or draw_tnorm(). Pólya–Gamma: z = (Y - 1/2) / omega, w = omega; probit: z from draw_probit_z(), w = 1; log times: ⁠z = log t⁠, w = 1 / sigma2. Missing z (NA) are skipped.

• draw_mvn(prec, num, lower, upper): multivariate normal $N(Q^{-1}h, Q^{-1})$ from the precision prec (Q) and linear term num (h); prec can be an array p x p x n and num a p x n matrix for n independent draws (one per item or person), optionally truncated to the box ⁠[lower, upper]⁠ (rejection; a slice that fails is NA).

• draw_items_mh(Y, theta, a, b, prior_mean, prior_prec, link, form): item parameters (a, b) with the latent responses integrated out, Metropolis–Hastings with a Fisher-scoring proposal (the cure for the slow mixing of easy and difficult items under data augmentation); logit or probit, form = "difficulty" for a (theta - b) or "intercept" for a theta - b; prior N((a, b); prior_mean (a row per item or one row), prior_prec^-1). Returns list(a, b) with the attribute accepted (as draw_items_joint()); draw the latent responses anew afterwards.

• draw_invwishart(df, scale): inverse Wishart.

• draw_mean_mvn(X, Sigma, mean0, cov0): the common mean of the rows of X ~ N(mu, Sigma) with prior mu ~ N(mean0, cov0) (hyperparameters of item- or person-level distributions).

• draw_cor(x1, x2, cor, steps, step, method): the correlation of two standard-normal variables (unit variances, the identification of several abilities), steps Metropolis steps on atanh(cor) with a uniform prior: random walk ("rw") or MALA ("mala", see draw_mala()).

• calc_person_fit(L, fitted, sigma2): the RT person-fit statistic of Marianti et al. (2014) $l^t_j = \sum_i (\log t_{ij} - \text{fitted}_{ij})^2 / \sigma_i^2$ per person, with the number of observed times df (missing times are skipped); chi-square with df degrees of freedom.

Value

calc_normal_part(): a normal part (prec, num; see normal_parts); draw_mvn(): a vector or a p x n matrix; draw_items_mh(): a list; draw_invwishart(): a matrix; draw_mean_mvn(): a vector; draw_cor(): a number with attribute accepted; calc_person_fit(): a data frame.

References

Marianti, S., Fox, J.-P., Avetisyan, M., Veldkamp, B. P., & Tijmstra, J. (2014). Testing for aberrant behavior in response time modeling. Journal of Educational and Behavioral Statistics, 39, 426-451.

Examples

set.seed(1)
N <- 300; K <- 5; theta <- rnorm(N); a <- rep(1.2, K); b <- seq(-1, 1, length.out = K)
Y <- matrix(rbinom(N * K, 1, plogis(sweep(outer(theta, a), 2, a * b))), N)
om <- draw_omega(a, b, theta)
# b_i | omega: eta = a (theta - b) enters the PG pseudo-observations with x = -a
g <- calc_normal_part(matrix(-a, N, K, byrow = TRUE), (Y - 0.5) / om, om, rest = outer(theta, a))
draw_gauss(1 / 9 + g$prec, g$num)                       # prior N(0, 3^2)
# (a, b) jointly with the latent responses integrated out
draw_items_mh(Y, theta, a, b, prior_mean = c(1, 0), prior_prec = diag(c(1, 1 / 9)))
draw_mvn(diag(2) * 4, c(1, 2), lower = 0)

Description

The share of item and structural parameters with effective sample size above ess and R-hat below rhat (over independent chains), and the worst parameters. With diagnostics = "rank" (see estimates()) both the bulk and the tail ESS must exceed ess, and the default R-hat threshold is 1.01 (Vehtari et al., 2021) instead of 1.05. Parameters that do not vary (fixed ones) are left out. With one chain R-hat cannot be computed and only the ESS is checked. print() of a gibbs() fit warns with the same criteria (default thresholds) when a parameter fails them. Julia: checkConvergence().

Usage

convergence(object, ...)

## S3 method for class 'rtirt'
convergence(
  object,
  ess = 400,
  rhat = NULL,
  diagnostics = c("classic", "rank"),
  ...
)

## S3 method for class 'rtirt_qset'
convergence(object, ...)

## S3 method for class 'rtirt_convergence'
print(x, ...)

Arguments

Details

For several quantile levels (an rtirt_qset from gibbs(..., quantile = c(.25, .75))) each level is checked; the model of one level is fits$models[[k]].

Value

An object of class rtirt_convergence, printed by its print() method: share (share of parameters that pass), worst (the five parameters with the largest R-hat, or the smallest ESS with one chain), n (parameters checked), n_chain and the thresholds; for an ECM fit converged and iterations; for an rtirt_qset one such object per level (levels) with quantile.

Examples

cond <- set_cond(n_subj = 200, n_item = 6)
dat <- sim_data(cond, sim_para(cond, "null"), "null")
fit <- gibbs(dat, model = "null", n_iter = 1000, n_chain = 2, seed = 1)
check <- convergence(fit)   # printed only when not assigned
check
check$share

Description

As getDic(MCMC) in Julia: $\bar D$ is the posterior mean of the complete-data deviance (after burn-in; the Julia version averaged over the burn-in as well), $\hat D$ the deviance at the posterior means, $p_D = \bar D - \hat D$. The deviance conditions on the person parameters (and, for the ALD models, on the mixing variables), so DIC values are not comparable between mean and quantile models or across quantile levels.

Usage

dic(object, ...)

## S3 method for class 'rtirt'
dic(object, ...)

Arguments

Description

Some directions of a model are identified only by the priors: adding the same constant to the abilities, the difficulties and the intercept of a latent regression leaves the likelihood of a 2PL unchanged. A Gibbs sampler moves along such a direction very slowly. draw_shift() draws the shift $c$ of the move $x_k \leftarrow x_k + c\, d_k$ exactly from its conditional (a generalized Gibbs step; Liu & Sabatti, 2000): only the normal priors of the shifted parameters depend on $c$. Parameters listed in values without a prior are shifted with a prior that moves along (e.g. abilities whose prior mean contains the shifted intercept).

Usage

draw_shift(values, prior, along = list(), lower = list())

Arguments

Details

The caller states that the likelihood is unchanged by the shift. Common cases:

Value

values shifted, with attribute shift.

References

Liu, J. S., & Sabatti, C. (2000). Generalised Gibbs sampler and multigrid Monte Carlo for Bayesian computation. Biometrika, 87, 353-369.

Examples

s <- draw_shift(list(theta = rnorm(5), b = c(-1, 0, 1), g0 = 0.2),
                prior = list(b = c(0, 3), g0 = c(0, 10)))
attr(s, "shift")

Description

Fits a built-in normal model by maximum marginal likelihood: the EM algorithm of Bock and Aitkin (1981) with ability integrated by Gauss-Hermite quadrature and speed in closed form. The item step uses the Polya-Gamma bound of the logistic likelihood (Polson, Scott & Windle, 2013), so every step is closed form and the log-likelihood never decreases; SQUAREM speeds up convergence. Standard errors come from the observed information (the numerical derivative of the analytic score).

Usage

ecm(data, ...)

## S3 method for class 'input_data'
ecm(
  data,
  model = NULL,
  itemtype = "2pl",
  speed_var = c("free", "fixed"),
  quantile = NULL,
  cov = NULL,
  prior = c("ml", "map"),
  priors = NULL,
  se = TRUE,
  control = ecm_control(),
  ...
)

Arguments

Details

Quantile models. With ⁠quantile =⁠ (models "latent" and "cross", as in gibbs()) ecm() maximizes the marginal ALD working likelihood (Yu & Moyeed, 2001): the q-quantile of the log times ("cross") or of speed given ability and the covariates ("latent") is linear in the parameters. Ability is integrated on the nodes and speed in closed form, so the working log-likelihood is smooth. 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 (two iterations per ECM step, which never lowers the working log-likelihood). The fit starts at the ECM fit of the mean model, moved to the q-quantile. The standard errors are sandwich (robust) standard errors, $H^{-1} J H^{-1}$ with $J$ the sum of the outer products of the scores of the persons: the ALD is not the distribution of the data, so the observed information alone would be wrong. As in gibbs(), sigma2t of "cross" is the ALD scale of the log times and var_speed of "latent" the ALD scale of speed. The working logLik, AIC and BIC compare only fits at the same quantile; anova() refuses quantile fits. A vector of levels gives one fit per level (an rtirt_qset).

The parameterization and identification are those of gibbs(): ability has mean ⁠x'beta⁠ and variance 1, speed mean ⁠x'beta⁠ (no intercepts: the covariates are centred and the item parameters carry the location). The fit starts the Gibbs sampler:

fit_ecm <- ecm(dat, model = "cross")
fit_gibbs <- gibbs(dat, model = "cross", init = fit_ecm)

Small samples. print() and summary() warn when the estimates suggest that the data are too few for maximum likelihood: a difficulty with |b| > 6, a discrimination below 0.1 or above 10, or a standard error above 5 (or missing) for a, b or rho. These thresholds flag only clearly broken fits (with ability N(0, 1), |b| > 6 means almost everyone or no one solves the item, and an SE of 5 a 95% interval about 20 wide); prior = "map" or gibbs() regularize such fits.

Value

A fitted model with ecm: estimates (est, se, z, Wald 95% interval on the log scale for positive parameters and the Fisher-z scale for the correlation), logLik, npar, AIC, BIC, iterations, converged, trace (log-likelihood per iteration), vcov, person (EAP and posterior SD) and secs. Accessors: print, summary, coef, estimates(), scores(), reliability(), fit_indices(), logLik, vcov, AIC, BIC.

References

Bock, R. D., & Aitkin, M. (1981). Marginal maximum likelihood estimation of item parameters: Application of an EM algorithm. Psychometrika, 46, 443-459.

Polson, N. G., Scott, J. G., & Windle, J. (2013). Bayesian inference for logistic models using Polya-Gamma latent variables. JASA, 108, 1339-1349.

White, H. (1982). Maximum likelihood estimation of misspecified models. Econometrica, 50, 1-25.

Yu, K., & Moyeed, R. A. (2001). Bayesian quantile regression. Statistics & Probability Letters, 54, 437-447.

Varadhan, R., & Roland, C. (2008). Simple and globally convergent methods for accelerating the convergence of any EM algorithm. Scandinavian Journal of Statistics, 35, 335-353.

Examples

cond <- set_cond(n_subj = 300, n_item = 8)
dat <- sim_data(cond, sim_para(cond, "cross"), "cross")
fit <- ecm(dat, model = "cross")
fit
head(estimates(fit))

fit_q <- ecm(dat, model = "cross", quantile = 0.9) # the 0.9-quantile of the log times
estimates(fit_q, pars = "^rho")

Description

The settings of ecm() and gibbs() that are rarely changed.

Usage

ecm_control(
  max_iter = 1000,
  tol = 1e-10,
  nodes = 41,
  accelerate = TRUE,
  verbose = FALSE,
  ...
)

gibbs_control(
  n_burnin = NULL,
  n_thin = 1,
  collapse = FALSE,
  fit_indices = TRUE,
  verbose = TRUE,
  ...
)

## S3 method for class 'ecm_control'
print(x, ...)

## S3 method for class 'gibbs_control'
print(x, ...)

Arguments

Value

A list of settings of class ecm_control or gibbs_control; ecm() refuses a gibbs_control() and gibbs() an ecm_control(). A plain list with the names of the settings is also accepted.

Examples

fit_settings <- gibbs_control(n_burnin = 500, collapse = TRUE)
fit_settings

Description

Deprecated with block_model(): use the built-in models, ecm(data, model = ...).

Usage

## S3 method for class 'block_model'
ecm(
  data,
  steps = NULL,
  latent = NULL,
  nodes = NULL,
  prior = c("map", "ml"),
  init = NULL,
  cycles = 6,
  max_iter = 1000,
  tol = 1e-10,
  se = TRUE,
  analytic = TRUE,
  quadrature = c("adaptive", "ba81"),
  ...
)

## S3 method for class 'block_model_ecm'
print(x, ...)

## S3 method for class 'block_model_ecm'
algorithm(object, format = c("unicode", "latex"), ...)

Arguments

Details

Fits a block_model() by EM with the same steps as gibbs(). Gibbs and EM share the full conditionals of the model: a step that a sampler uses to draw a parameter is used by EM to move the parameter to the mode of its full conditional averaged over the latent variables.

• E-step: the person parameters (latent, by default every ⁠[person]⁠ parameter) are integrated out by adaptive Gauss–Hermite quadrature: a product grid of nodes points per latent variable, centred at each person's posterior mean and scaled by the posterior SD. The centres are fixed while EM runs (so the objective is one fixed function and EM never lowers it), then recentred, for up to cycles rounds. The Pólya–Gamma variables of a logistic part are replaced by their expectation $E[\omega] = \tanh(\eta/2)/(2\eta)$ at each node (a minorizer: the step still never lowers the objective).

• step_gibbs() becomes a conditional maximization: the normal part of the parameter is summed over the model and over the nodes (weighted by the posterior of each person), and the parameter moves to its mode $m/p$ (clipped to the bounds); a variance moves to the mode of its inverse-gamma conditional (prior = "ml": $\beta/\alpha$).

• step_mala() becomes a few gradient-ascent steps with step halving on the same weighted objective (generalized EM), on the log scale for bounded parameters (no Jacobian: a mode does not need one).

• step_shift() is skipped: the latent variables are integrated out.

The objective is the log marginal likelihood (prior = "ml"; the priors of the latent variables stay, they define the marginal likelihood) or the log marginal posterior (prior = "map", the default: also the priors of the other parameters). Standard errors come from its Hessian (Louis' formula: expected complete-data Hessian plus the covariance of the complete-data score, from the derivatives of the declared expressions); the posterior means and SDs of the latent variables are in fit$latent. EM runs through run_em() (SQUAREM, the checks of run_em()). Use check_ecm() to check the steps.

Value

A block_em (see run_em()) with latent (posterior means and SDs per person), model, steps and cycles; algorithm(fit) lists the E- and M-steps.

Theory of each step

• ECM. EM needs only the conditional expectations of the complete-data sufficient statistics (Dempster, Laird & Rubin, 1977, Sec. 2) and never lowers the objective when every M-step raises Q (their Theorem 1, generalized EM); the M-step may be a sequence of conditional maximizations (Meng & Rubin, 1993).

• Adaptive quadrature. Nodes centred and scaled at each person's posterior (Naylor & Smith, 1982; Liu & Pierce, 1994; Bock & Gibbons, 2021, Sec. 6.4). Schilling & Bock (2005) recentre at every EM iteration; here the centres are kept fixed within a round and recentred between rounds, as Bock & Gibbons (2021, pp. 210–211) suggest to save time and as in the pseudo-adaptive rule of Rizopoulos (2012). With fixed centres the approximated objective is one fixed function, so the monotonicity of Dempster et al. (1977, Theorem 1) applies to it (our argument; no source states it).

• quadrature = "ba81". Common nodes and expected counts per node ("artificial data"): Bock & Aitkin (1981); Bock & Gibbons (2021, Sec. 6.4, Eqs. 6.17–6.22).

• Logistic parts. $E[\omega] = \tanh(\eta/2)/(2\eta)$ (Polson, Scott & Windle, 2013, Sec. 2.2) makes the M-step an exact EM step for the Pólya–Gamma augmented model (Scott & Sun, 2013); it equals the Jaakkola–Jordan bound, a minorize–maximize step that never lowers the likelihood (Durante & Rigon, 2019, Appendix).

• Exact integration (analytic). A latent variable u that enters only normal parts linearly, with variances free of u, has a normal conditional $u \mid \cdot \sim N(h/p, 1/p)$ (completing the square; Bishop, 2006, Sec. 2.3; for speed in the hierarchical RT model, van der Linden, 2007, Eq. 42), and the integral over u is the integrand at the mode times $\sqrt{2\pi/p}$, exact for a quadratic log integrand (Bishop, 2006, Eq. 4.135). The RT-IRT marginal likelihood papers we found integrate speed numerically together with ability (Glas & van der Linden, 2010; Molenaar & Bolsinova, 2017); integrating it exactly is our reduction, valid only under the condition above (checked from the declaration; not for heteroscedastic or skewed speed).

• Exact variables in the M-step. When every part of an exact variable is linear in all latent variables jointly, with variances free of them, these parts are quadratic in the latent variables, and the expectation of a quadratic depends only on the mean and covariance ($E[x'Ax] = \mathrm{tr}(A\Sigma) + \mu'A\mu$, for any distribution). The M-step therefore uses each person's posterior mean and covariance of the latent variables through the 2L symmetric points $m \pm \sqrt{L}$ x columns of the Cholesky factor, weight 1/(2L) each, which have exactly that mean and covariance (Julier & Uhlmann, 2004, Eq. 11). Logistic parts stay on the nodes.

• Standard errors. Louis (1982): observed information = expected complete-data information minus the covariance of the complete-data score. Under prior = "map" the log prior of the fixed parameters does not involve the latent variables, so its Hessian is added to that of the log marginal likelihood (Chalmers, 2018, notes that how priors enter Louis' form needs saying). The squared scores are quartic in an exact variable, so three Gauss–Hermite points per exact variable are used (n points are exact to degree 2n - 1; DLMF Sec. 3.5(v)).

References

Bishop, C. M. (2006). Pattern recognition and machine learning. Springer.

Bock, R. D., & Aitkin, M. (1981). Marginal maximum likelihood estimation of item parameters: Application of an EM algorithm. Psychometrika, 46(4), 443–459. doi:10.1007/BF02293801

Bock, R. D., & Gibbons, R. D. (2021). Item response theory. Wiley. doi:10.1002/9781119716723

Chalmers, R. P. (2018). Numerical approximation of the observed information matrix with Oakes' identity. British Journal of Mathematical and Statistical Psychology, 71(3), 415–436. doi:10.1111/bmsp.12127

Dempster, A. P., Laird, N. M., & Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society B, 39(1), 1–22. doi:10.1111/j.2517-6161.1977.tb01600.x

DLMF. NIST Digital Library of Mathematical Functions, Sec. 3.5(v), Gauss quadrature. https://dlmf.nist.gov/3.5

Durante, D., & Rigon, T. (2019). Conditionally conjugate mean-field variational Bayes for logistic models. Statistical Science, 34(3), 472–485. doi:10.1214/19-STS712

Glas, C. A. W., & van der Linden, W. J. (2010). Marginal likelihood inference for a model for item responses and response times. British Journal of Mathematical and Statistical Psychology, 63(3), 603–626. doi:10.1348/000711009X481360

Julier, S. J., & Uhlmann, J. K. (2004). Unscented filtering and nonlinear estimation. Proceedings of the IEEE, 92(3), 401–422. doi:10.1109/JPROC.2003.823141

Liu, Q., & Pierce, D. A. (1994). A note on Gauss–Hermite quadrature. Biometrika, 81(3), 624–629. doi:10.1093/biomet/81.3.624

Louis, T. A. (1982). Finding the observed information matrix when using the EM algorithm. Journal of the Royal Statistical Society B, 44(2), 226–233. doi:10.1111/j.2517-6161.1982.tb01203.x

Meng, X.-L., & Rubin, D. B. (1993). Maximum likelihood estimation via the ECM algorithm: A general framework. Biometrika, 80(2), 267–278. doi:10.1093/biomet/80.2.267

Molenaar, D., & Bolsinova, M. (2017). A heteroscedastic generalized linear model with a non-normal speed factor for responses and response times. British Journal of Mathematical and Statistical Psychology, 70(2), 297–316. doi:10.1111/bmsp.12087

Naylor, J. C., & Smith, A. F. M. (1982). Applications of a method for the efficient computation of posterior distributions. Applied Statistics, 31(3), 214–225. doi:10.2307/2347995

Polson, N. G., Scott, J. G., & Windle, J. (2013). Bayesian inference for logistic models using Pólya–Gamma latent variables. Journal of the American Statistical Association, 108(504), 1339–1349. doi:10.1080/01621459.2013.829001

Rizopoulos, D. (2012). Fast fitting of joint models for longitudinal and event time data using a pseudo-adaptive Gaussian quadrature rule. Computational Statistics & Data Analysis, 56(3), 491–501. doi:10.1016/j.csda.2011.09.007

Schilling, S., & Bock, R. D. (2005). High-dimensional maximum marginal likelihood item factor analysis by adaptive quadrature. Psychometrika, 70(3), 533–555. doi:10.1007/s11336-003-1141-x

Scott, J. G., & Sun, L. (2013). Expectation-maximization for logistic regression. arXiv:1306.0040.

van der Linden, W. J. (2007). A hierarchical framework for modeling speed and accuracy on test items. Psychometrika, 72(3), 287–308. doi:10.1007/s11336-006-1478-z

Examples

set.seed(1); N <- 300; K <- 6
th <- rnorm(N); a0 <- runif(K, 0.8, 2); b0 <- seq(-1.5, 1.5, length.out = K)
Y <- matrix(rbinom(N * K, 1, plogis(sweep(outer(th, a0), 2, a0 * b0))), N)
m <- block_model(Y = Y,
  theta[person] ~ normal(0, 1),
  a[item] ~ normal(1, 1, lower = 0),
  b[item] ~ normal(0, 3),
  eta <- a * (theta - b),
  Y ~ bernoulli_logit(eta))
fit <- ecm(m, list(step_gibbs(a, b)), prior = "ml")
fit
algorithm(fit)

Description

The steps of ecm() as functions that can be combined into an EM for another model (the counterpart of the sampling blocks, blocks and compose: ⁠draw_*()⁠ draws a block from its full conditional, ⁠update_*()⁠ maximizes, or takes a monotone MM step, in the same block).

Usage

calc_nodes(Q = 41)

calc_loglik_nodes(Y, a, b, theta)

calc_normal_given_nodes(o, w, m, v)

calc_node_weights(logL, logw = 0)

calc_artificial_data(W, Y)

calc_item_sums(Y, a, b, theta, W, score = FALSE)

update_items_logistic(S, T0, T1, a, b, prior_a = NULL, prior_b = NULL)

em_squarem(
  x,
  Fmap,
  max_iter = 1000,
  tol = 1e-10,
  accelerate = TRUE,
  progress = NULL
)

calc_hessian(score, x, h = 1e-04)

Arguments

Details

E-step

• calc_nodes(Q): Gauss–Hermite nodes z and log weights logw of N(0, 1).

• calc_loglik_nodes(Y, a, b, theta): the 2PL log-likelihood of every person at every node, $\sum_i y_{ij}\eta - \log(1 + e^\eta)$, $\eta = a_i(\theta - b_i)$, an N x Q matrix. theta is the vector of nodes when they are common to all persons (two matrix products, Bock–Aitkin), or an N x Q matrix of person-specific nodes (latent regressions; C++).

• calc_normal_given_nodes(o, w, m, v): a normal latent variable u (e.g. speed) with prior N(m, v) at each node (m: N x Q, or a vector of Q), observed through normal pseudo-observations $o_{ij} = u_j + e_{ij}$, $e_{ij} \sim N(0, 1/w_i)$: integrated in closed form. Returns logI (N x Q, add to the log-likelihood), the conditional mean (N x Q) and var of u.

• calc_node_weights(logL, logw): posterior weights W of the nodes (N x Q) and the log-likelihood ll of every person.

• calc_artificial_data(W, Y): Bock–Aitkin's expected counts: n (Q x K, persons at each node answering each item) and r (Q x K, correct answers).

• calc_item_sums(Y, a, b, theta, W, score): weighted item sums over persons and nodes, either the Pólya–Gamma minorizer sums $S_r = \sum W E[\omega|\eta]\theta^r$ (score = FALSE) or the score sums $\sum W (y - P)\theta^r$ (score = TRUE), r = 0, 1, 2; common nodes go through the artificial data.

M-step

• update_items_logistic(S, T0, T1, a, b, prior_a, prior_b): the Pólya–Gamma (Jaakkola–Jordan) minorizer step for 2PL items: a 2 x 2 weighted least-squares problem per item from the sums S of calc_item_sums() and $T_0 = \sum_j (y_{ij} - 1/2)$, $T_1 = \sum_j (y_{ij} - 1/2) E[\theta_j]$; the marginal log-likelihood never decreases. With prior_a = c(mean, sd) and prior_b = c(mean, sd, lower, upper) (or a list whose mean has one value per item) the conditional modes (MAP).

Driver and information

• em_squarem(x, Fmap, ...): SQUAREM (Varadhan & Roland, 2008; SqS3) around a fixed-point map Fmap(x) returning list(x = , obj = ) (one EM step from x and the log-likelihood at x), with a fall-back to plain EM steps when the extrapolation lowers the objective.

• calc_hessian(score, x): the symmetrized numerical Jacobian of an analytic score function (central differences); $(-H)^{-1}$ is the covariance of the estimates.

Value

See Details.

References

Bock, R. D., & Aitkin, M. (1981). Marginal maximum likelihood estimation of item parameters: Application of an EM algorithm. Psychometrika, 46, 443-459.

Varadhan, R., & Roland, C. (2008). Simple and globally convergent methods for accelerating the convergence of any EM algorithm. Scandinavian Journal of Statistics, 35, 335-353.

Examples

set.seed(1)
N <- 500; K <- 10; a <- runif(K, 0.8, 2); b <- rnorm(K)
Y <- matrix(rbinom(N * K, 1, plogis(sweep(outer(rnorm(N), a), 2, a * b))), N)
# an EM for the 2PL from the components
nd <- calc_nodes(41); kap <- Y - 0.5
Fmap <- function(x) {
  a <- exp(x[1:K]); b <- x[K + 1:K]
  nw <- calc_node_weights(calc_loglik_nodes(Y, a, b, nd$z), nd$logw)
  S <- calc_item_sums(Y, a, b, nd$z, nw$W)
  new <- update_items_logistic(S, colSums(kap), colSums(kap * drop(nw$W %*% nd$z)), a, b)
  list(x = c(log(new$a), new$b), obj = sum(nw$ll))
}
fit <- em_squarem(c(rep(0, K), rep(0, K)), Fmap)
round(cbind(a, a_hat = exp(fit$x[1:K]), b, b_hat = fit$x[K + 1:K]), 2)

Description

Writes the model as equations, with the parameter names of summary() and the priors of gibbs(). The default is Unicode text for the console; format = "latex" gives an aligned block for a paper or an appendix (e.g. writeLines(as.character(equations(fit, "latex")), "model.tex")). Works before estimation too: give the data and the model as in ecm() and gibbs().

Usage

equations(object, ...)

## S3 method for class 'rtirt_equations'
as.character(x, ...)

## S3 method for class 'rtirt_equations'
print(x, ...)

## S3 method for class 'rtirt'
equations(object, format = c("unicode", "latex"), priors = TRUE, ...)

## S3 method for class 'rtirt_qset'
equations(object, format = c("unicode", "latex"), priors = TRUE, ...)

## S3 method for class 'input_data'
equations(
  object,
  model = NULL,
  itemtype = "2pl",
  speed_var = c("free", "fixed"),
  quantile = NULL,
  cov = NULL,
  intercept = FALSE,
  format = c("unicode", "latex"),
  priors = TRUE,
  ...
)

Arguments

Value

An rtirt_equations object (printed); as.character() gives the lines.

Examples

cond <- set_cond(n_subj = 100, n_item = 5)
dat <- sim_data(cond, sim_para(cond, "cross"), "cross")
equations(dat, model = "cross")
equations(dat, model = "cross", quantile = 0.25, format = "latex")

Description

estfun() returns, for a model fitted by ecm(), the matrix of casewise scores: one row per person, one column per parameter, the derivative of that person's marginal log-likelihood (ability and speed integrated out). The columns sum to zero at the estimates. The parameters are on the unconstrained scale used by ecm(): log a, b, lambda, log sigma2t, then the structural parameters (for model = "latreg" / "null" with free speed variance: the regression coefficients, the covariance c and the log residual variance log_v of speed given ability). These are the estimating functions used by score_test() and by strucchange::gefp().

Usage

## S3 method for class 'rtirt'
estfun(x, ...)

## S3 method for class 'rtirt'
bread(x, ...)

Arguments

Value

A persons x parameters matrix.

bread(): the inverse of the average observed information (unconstrained scale), as in the 'sandwich' package.

Description

estimates() returns the posterior summary of the item and structural parameters as a data frame: posterior mean (est), sd, 2.5% and 97.5% quantiles (q025, q975), R-hat and effective sample size (after burn-in, over the independent chains); sig = "*" when the 95% interval excludes 0. sig is given only for parameters where 0 is a meaningful value: difficulties b, regression coefficients beta, cross-relations rho, cor_ability_speed and b_ability. It is NA for discriminations, variances (sigma2t, var_speed), which are positive and whose intervals never contain 0, and time intensities lambda. precis() (the Julia name) returns the same table, printed in a plain layout (assigning it prints nothing); as.data.frame() gives it as a plain data frame, e.g. for ggplot2. coef() returns the posterior means as a named vector. For a model fitted by ecm() the table has the maximum likelihood estimates with standard errors (est, se, z, Wald 95% interval q025, q975); z, like sig, is NA where 0 is not a meaningful value, and summary() then shows no z or p.

Usage

estimates(object, ...)

## S3 method for class 'rtirt'
estimates(
  object,
  pars = NULL,
  method = c("auto", "gibbs", "ecm"),
  diagnostics = c("classic", "rank"),
  ...
)

precis(object, ...)

## S3 method for class 'rtirt'
precis(object, pars = NULL, digits = 3, ...)

## S3 method for class 'rtirt_precis'
print(x, digits = attr(x, "digits"), ...)

## S3 method for class 'rtirt'
as.data.frame(x, row.names = NULL, optional = FALSE, ...)

## S3 method for class 'rtirt'
coef(object, ...)

Arguments

Value

estimates(): a tibble of class rtirt_tbl (printed with 3 decimals; x$est and as.data.frame(x) keep the full precision). coef(): a named vector, printed with 3 decimals.

References

Gelman, A., & Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science, 7, 457-472.

Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., & Bürkner, P.-C. (2021). Rank-normalization, folding, and localization: An improved R-hat for assessing convergence of MCMC. Bayesian Analysis, 16, 667-718.

Description

Criteria from the marginal log-likelihood (person parameters integrated out; see marginal_loglik()), which compare models with different latent structures: DIC (with pD = mean deviance - deviance at the posterior means of the item and structural parameters), WAIC (p_waic) and PSIS-LOO (LOOIC, SE_LOOIC; needs the 'loo' package). DIC_complete is the complete-data DIC of dic() (conditional on the person parameters; for comparisons within one model only). Quantile models use the ALD as a working likelihood, so their criteria are not comparable with mean models or across quantile levels. For a model fitted by ecm(): the log-likelihood, AIC and BIC.

Usage

fit_indices(object, ...)

## S3 method for class 'rtirt'
fit_indices(object, ...)

Arguments

Description

Bias, RMSE and correlation

Usage

get_bias(est, true)

get_rmse(est, true)

get_corr(est, true)

Arguments

Description

Deprecated with block_model(): use the built-in models, gibbs(data, model = ...).

Usage

## S3 method for class 'block_model'
gibbs(
  data,
  steps = NULL,
  n_iter = 2000,
  n_burnin = floor(n_iter/2),
  n_chain = 2,
  n_thin = 1,
  seed = 1,
  init = NULL,
  monitor = NULL,
  ...
)

## S3 method for class 'block_model'
algorithm(object, steps = NULL, format = c("unicode", "latex"), ...)

## S3 method for class 'block_model_fit'
algorithm(object, format = c("unicode", "latex"), ...)

Arguments

Details

Runs the steps once per iteration in the order given (see block_steps) with run_sampler() (chains, burn-in, storage, summary). Without steps, every parameter that has a Gibbs step gets one (in the order of the declarations) and the others get MALA.

Value

A block_fit (see run_sampler()) with the model and the steps; algorithm(fit) lists the steps with their full conditionals.

Description

Samples a built-in model from its full conditionals: n_chain independent chains, each from its own random starting values or from an ecm() fit.

Usage

## S3 method for class 'input_data'
gibbs(
  data,
  model = NULL,
  itemtype = "2pl",
  speed_var = c("free", "fixed"),
  quantile = NULL,
  cov = NULL,
  intercept = FALSE,
  priors = NULL,
  n_iter = 5000,
  n_chain = 4,
  seed = NULL,
  init = "random",
  control = gibbs_control(),
  ...
)

gibbs(data, ...)

Arguments

Details

fit_ecm <- ecm(dat, model = "cross")
fit_gibbs <- gibbs(dat, model = "cross", init = fit_ecm)

Value

A fitted model with post: draws (a coda::mcmc.list of the item and structural parameters, all iterations), mean (posterior means after burn-in, an input_para()), person (posterior mean and SD of ability and speed), loglik (complete-data log-likelihood per iteration and chain) and secs. Accessors: summary(), coef(), estimates(), scores(), reliability(), convergence(), fit_indices(), plot().

References

Wang, C., & Xu, G. (2015). A mixture hierarchical model for response times and response accuracy. British Journal of Mathematical and Statistical Psychology, 68, 456-477.

Examples

cond <- set_cond(n_subj = 300, n_item = 8)
dat <- sim_data(cond, sim_para(cond, "cross"), "cross")
fit <- gibbs(dat, model = "cross", n_iter = 1000, n_chain = 2, seed = 1)
summary(fit)
head(scores(fit))

Description

Collects the item responses, response times and covariates of one test (mirrors InputData() of the Julia package). Each block is either a matrix / data frame or, with data, a vector of column names:

Usage

input_data(
  resp,
  time = NULL,
  cov = NULL,
  data = NULL,
  id = NULL,
  log_time = NULL
)

Arguments

Details

dat <- input_data(resp = items, time = paste0(items, "_S"), cov = covs, id = "IDSTUD", data = raw)
dat <- input_data(resp = Y, time = RT, cov = X)

Item names (the column names of resp) label all output. Missing values are not supported: drop the incomplete persons first.

The data must be in wide format: one row per person, one column per item (for responses and for times). Data with one row per person and item (long format, as exported by PsychoPy or E-Prime) are refused, recognized by a single item or by ids that repeat on most rows; reshape them first, e.g. tidyr::pivot_wider(long, id_cols = id, names_from = item, values_from = c(correct, rt)) or stats::reshape(long, idvar = "id", timevar = "item", direction = "wide"). A few duplicated ids give a warning. Covariates must be numeric: dummy-code categorical ones first, e.g. model.matrix(~ grade, data)[, -1] (k - 1 columns for k categories), and centre them.

Responses coded as text or factors ("correct" / "wrong") must be recoded to 0/1 first; the message shows the values found and the recoding. TRUE / FALSE are accepted.

Units of time. Seconds are expected, but the unit only shifts the time intensities: times in milliseconds give every lambda larger by log(1000) = 6.91 (minutes: smaller by log(60) = 4.09), and a, b, rho, the residual variances and the correlations do not change. Report the unit with lambda. The default prior of lambda in gibbs() (N(0, 10^2)) is meant for seconds, so a median time above 1000 gives a warning that suggests milliseconds.

Value

An object of class input_data with Y, kappa = Y - 1/2, T, log_t, X, items and id.

Description

Mirrors InputPara: a list of the model parameters (theta, a, b, zeta, lambda, sigma2t, nu, beta, rho, b_theta, sigma_p). Used for true values in simulations and for posterior means (fit$post$mean).

Usage

input_para(...)

Arguments

Description

The function names of ExtendedRtIrtModeling.jl, mapped to the snake_case functions:

Usage

setCond(
  nSubj = 2000,
  nItem = 15,
  nFeat = 3,
  nIter = 5000,
  nChain = 4,
  nBurnin = NULL,
  nThin = 1,
  nRep = 10,
  qRa = 0.5,
  qRt = 0.5
)

InputData(Y, T = NULL, X = NULL)

InputData4R(Y, kappa = Y - 0.5, logT = NULL, X = NULL)

InputPara(...)

GibbsMlIrt(Cond, Data, truePara = NULL, Para = NULL)

GibbsRtIrtNull(Cond, Data, truePara = NULL, Para = NULL)

GibbsRtIrt(Cond, Data, truePara = NULL, Para = NULL)

GibbsRtIrtLatent(Cond, Data, truePara = NULL, Para = NULL)

GibbsRtIrtLatentQr(Cond, Data, truePara = NULL, Para = NULL)

GibbsRtIrtCross(Cond, Data, truePara = NULL, Para = NULL)

GibbsRtIrtCrossQr(Cond, Data, truePara = NULL, Para = NULL)

sampleGibbs(MCMC, cov2one = NULL, ...)

getDic(MCMC)

checkConvergence(MCMC, ...)

comparePara(MCMC, par = "a")

getBias(est, true)

getRmse(est, true)

getCorr(est, true)

setTrueParaMlIrt(Cond)

setTrueParaRtIrt(Cond, trueCorr = 0.3)

setTrueParaRtIrtLatent(Cond)

setTrueParaRtIrtCross(Cond)

setDataMlIrt(Cond, truePara)

setDataRtIrt(Cond, truePara)

setDataRtIrtNull(Cond, truePara)

setDataRtIrtLatent(Cond, truePara, type = "norm")

setDataRtIrtCross(Cond, truePara, type = "norm")

runSimulation(
  Cond,
  truePara,
  funcData = setDataRtIrt,
  funcGibbs = GibbsRtIrt,
  par = c("a", "b", "theta"),
  ...
)

Arguments

Details

The constructors take ⁠(Cond, Data, truePara, Para)⁠ as in Julia and keep the Julia defaults: the speed variance is fixed at 1 (speed_var = "fixed") in GibbsRtIrtNull, GibbsRtIrt and GibbsRtIrtCross (sampleGibbs(M, cov2one = FALSE) frees it), and the Qr versions use the quantile level Cond$q_rt. setCond() keeps the Julia default dimensions (2000 persons, 15 items, 3 covariates), which must match the data.

Description

A Metropolis–Hastings step whose proposal follows the gradient of the log density (MALA; Roberts & Tweedie, 1996): $x^* = x + \frac{\epsilon^2}{2}\nabla\log\pi(x) + \epsilon z$, accepted with the usual MH probability (the proposal is not symmetric, so its density in both directions enters the ratio). It replaces a random-walk Metropolis step where a block has no conjugate full conditional: it moves in the direction of higher density and needs far fewer rejections, especially in several dimensions (optimal acceptance about 0.574 against 0.234).

Usage

draw_mala(x, logpost, step, joint = FALSE)

mcmc_adapt(step, accepted, iter, target = 0.574)

Arguments

Details

• draw_mala(x, logpost, step, joint): by default one step for many independent targets at once, one per row of x (or per element of a vector x, e.g. one per person); logpost(x) returns list(lp = , grad = ): the log density of every row (a vector) and its gradient (like x); step is one number or one per row. With joint = TRUE, x is one target and lp one number. Returns x with the attribute accepted.

• A parameter on a transformed scale (e.g. ⁠u = log a⁠) needs the log Jacobian in lp (+ u) and in grad (+ 1); variables integrated out by the proposal (e.g. Pólya–Gamma omega) must be drawn anew before the next step uses them.

• mcmc_adapt(step, accepted, iter, target): Robbins–Monro adaptation of the step size during burn-in (Andrieu & Thoms, 2008), $\log\epsilon \leftarrow \log\epsilon + t^{-0.6}(\alpha - \alpha^*)$; stop adapting after burn-in so that the chain is a proper Markov chain.

draw_cor() uses it with method = "mala".

Value

draw_mala(): the new x with attribute accepted; mcmc_adapt(): the new step(s).

References

Roberts, G. O., & Tweedie, R. L. (1996). Exponential convergence of Langevin distributions and their discrete approximations. Bernoulli, 2, 341-363.

Andrieu, C., & Thoms, J. (2008). A tutorial on adaptive MCMC. Statistics and Computing, 18, 343-373.

Examples

# abilities of a 2PL with known items: one MALA step per person, all persons at once
set.seed(1); N <- 200; a <- c(1, 1.5, 2); b <- c(-1, 0, 1)
theta0 <- rnorm(N); Y <- matrix(rbinom(N * 3, 1, plogis(sweep(outer(theta0, a), 2, a * b))), N)
lp <- function(th) { p <- plogis(sweep(outer(th, a), 2, a * b))
  list(lp = rowSums(Y * log(p) + (1 - Y) * log(1 - p)) - th^2 / 2,
       grad = drop((Y - p) %*% a) - th) }
th <- rep(0, N); eps <- rep(1, N)
for (t in 1:2000) {
  th <- draw_mala(th, lp, eps)
  if (t <= 500) eps <- mcmc_adapt(eps, attr(th, "accepted"), t)
}
mean(attr(th, "accepted"))

Description

marginal_loglik() returns the log-likelihood of each person with the person parameters (ability, speed) integrated out, for posterior draws of the item and structural parameters. fit_indices() uses it for DIC, WAIC and PSIS-LOO (with the 'loo' package); these, unlike the complete-data DIC of dic(), compare models with different latent structures (Merkle, Furr & Rabe-Hesketh, 2019). Ability is integrated by Gauss-Hermite quadrature and speed in closed form (normal log times; normal-ALD convolution for rtirt_latent quantile models; piecewise-linear ALD log-likelihood for rtirt_cross). Quantile models use the ALD as a working likelihood, so their criteria are not comparable with mean models or across levels.

Usage

marginal_loglik(object, n_draws = NULL, Q = NULL)

Arguments

Value

A draws x persons matrix (printed with 3 decimals; see rtirt_tbl).

Description

Deprecated with block_model(): use the built-in models, ecm(data, model = ...).

Usage

mml(model, ...)

## S3 method for class 'block_model'
mml(
  model,
  latent = NULL,
  nodes = NULL,
  prior = c("map", "ml"),
  init = NULL,
  cycles = 6,
  max_iter = 1000,
  tol = 1e-09,
  se = TRUE,
  analytic = TRUE,
  quadrature = c("adaptive", "ba81"),
  method = c("newton", "quasi-newton"),
  ...
)

Arguments

Details

Maximizes the objective of ecm.block_model() (the marginal log-likelihood, or the log posterior with prior = "map") directly with stats::nlminb() (within the bounds of the parameters), all parameters together, with the gradient from Fisher's identity and, for method = "newton", the Hessian from Louis' formula (Glas & van der Linden, 2010, Sec. 3.1, use Newton-Raphson or EM). No closed-form step is needed, only the derivatives of the declared expressions, so this helps when a parameter has no closed-form step. It reaches the same maximum as ecm() inside the bounds of the parameters (ECM keeps a bounded parameter 1e-3 inside its bound; a parameter at its bound gets no standard error, with a warning); ECM never lowers the objective and its iterations are cheap, Newton needs few iterations.

Value

As ecm.block_model().

References

Glas, C. A. W., & van der Linden, W. J. (2010). Marginal likelihood inference for a model for item responses and response times. British Journal of Mathematical and Statistical Psychology, 63, 603-626.

Description

Deprecated: name the model in ecm() or gibbs() instead, e.g. rtirt_cross(dat) becomes ecm(dat, model = "cross") or gibbs(dat, model = "cross"). The constructors still work (with a warning once per session) and get no new features.

Usage

mlirt(
  data,
  itemtype = c("2pl", "1pl"),
  intercept = FALSE,
  cond = NULL,
  true_para = NULL,
  para = NULL
)

rtirt_null(
  data,
  speed_var = c("free", "fixed"),
  itemtype = c("2pl", "1pl"),
  cond = NULL,
  true_para = NULL,
  para = NULL
)

rtirt_latreg(
  data,
  speed_var = c("free", "fixed"),
  itemtype = c("2pl", "1pl"),
  intercept = FALSE,
  cond = NULL,
  true_para = NULL,
  para = NULL
)

rtirt_latent(
  data,
  quantile = NULL,
  itemtype = c("2pl", "1pl"),
  intercept = FALSE,
  cond = NULL,
  true_para = NULL,
  para = NULL
)

rtirt_cross(
  data,
  quantile = NULL,
  speed_var = c("free", "fixed"),
  itemtype = c("2pl", "1pl"),
  cond = NULL,
  true_para = NULL,
  para = NULL
)

Arguments

Details

fit <- gibbs(dat, model = "cross", quantile = 0.25, n_iter = 10000, n_chain = 3, seed = 1)

• mlirt: 2PL with a latent regression of ability (no response times).

• rtirt_null: joint RT-IRT (van der Linden, 2007), ability and speed correlated.

• rtirt_latreg: joint RT-IRT with latent regressions of ability and speed on covariates (dissertation Chapter 3).

• rtirt_latent: speed regressed on ability and covariates, mean regression or, with quantile, quantile regression (ALD; Chapter 4).

• rtirt_cross: cross-relations of ability with each log response time, normal or, with quantile, ALD residuals (Chapter 5).

The Julia constructors (GibbsRtIrtCross(Cond, Data) etc.) are aliases with the Julia argument order and defaults (see julia_aliases).

Value

A model object of class ⁠c(<model>, "rtirt")⁠; see gibbs().

Description

Most full conditionals of the samplers are normal. A normal part stores one as its precision $Q$ and its linear term $h$ (prec and num; mean $Q^{-1}h$, variance $Q^{-1}$), so that the contributions of the data and of the prior simply add up:

Usage

normal_part(prec, num)

normal_prior(mean = 0, sd = 1)

reg_part(X, y, w = 1)

## S3 method for class 'normal_part'
e1 + e2

draw(part, lower = -Inf, upper = Inf)

## S3 method for class 'normal_part'
print(x, ...)

Arguments

Details

• irt_part() and calc_normal_part() return the data part of independent scalars (one per person or item);

• reg_part(X, y, w) returns the data part of a vector of regression coefficients $g$ in $y = Xg + e$, $e_j \sim N(0, 1/w_j)$ ($Q = X'WX$, $h = X'Wy$), e.g. a latent regression of the abilities on covariates;

• normal_prior(mean, sd) turns a normal prior into a part ($Q = 1/sd^2$, $h = mean/sd^2$);

• part1 + part2 adds parts (elementwise, or a prior added to the diagonal of a regression part);

• draw(part, lower, upper) draws from it (truncated to ⁠[lower, upper]⁠ when given; a regression part is drawn with draw_mvn());

• print(part) shows the conditional means and SDs.