run_mst: Multistage-Adaptive Test (MST) Simulation

In irtQ: Unidimensional Item Response Theory Modeling

Multistage-Adaptive Test (MST) Simulation

Description

Simulates an MST administration for a panel of examinees. At each stage, examinees are routed to a module based on their current ability estimate (or a set of cut scores), and a final ability estimate is computed from the accumulated responses across all administered stages. The function supports all IRT models available in irtQ (1PLM, 2PLM, 3PLM, GRM, GPCM) and all scoring methods available in est_score.

The overall structure - separate route_score and final_score list arguments, and the route_map/module matrix conventions - was inspired by the randomMST() function in the mstR package (Magis et al., 2017).

Usage

run_mst(
  x,
  route_map,
  module,
  theta = NULL,
  response = NULL,
  D = 1,
  ini_mod = NULL,
  route_method = "bmat",
  cut_score = NULL,
  route_score = list(method = "ML", range = c(-5, 5), norm.prior = c(0, 1), nquad = 41L,
    tol = 1e-04, max.iter = 100L, fence.a = 3, fence.b = NULL, intpol = TRUE, range.tcc =
    c(-7, 7), max.it = 500L),
  final_score = list(method = "ML", range = c(-5, 5), norm.prior = c(0, 1), nquad = 41L,
    tol = 1e-04, max.iter = 100L, fence.a = 3, fence.b = NULL, intpol = TRUE, range.tcc =
    c(-7, 7), max.it = 500L),
  se = TRUE,
  missing = NA,
  verbose = TRUE,
  return_full_resp = FALSE
)

Arguments

x	A data frame of item bank metadata following the irtQ format (columns: id, cats, model, par.1, par.2, ...). Use shape_df to create this frame or bring.flexmirt / bring.bilog / bring.parscale / bring.mirt to import from IRT software output. The function calls confirm_df(x) internally to validate and normalise the metadata.
route_map	A binary square matrix defining the MST transition structure. A 1 at row i, column j means that a test taker can be routed from module i to module j. Equivalent to the transMatrix argument in randomMST() from mstR (Magis et al., 2017). See reval_mst for details.
module	A binary matrix mapping items (rows) to modules (columns). module[j, m] = 1 means item j belongs to module m. Equivalent to the modules argument in randomMST() from mstR (Magis et al., 2017). Must have nrow(module) == nrow(x) and ncol(module) == total number of modules in the panel.
theta	A numeric vector of true ability values, one per examinee (length N). Used to simulate item responses when response is NULL, and stored in the result for bias/RMSE computation. Either theta or response (or both) must be provided. If theta is omitted and only response is given, true.theta in the result will be NULL and bias/RMSE cannot be computed.
response	An optional N x J matrix of observed item responses, where N is the number of examinees and J is the total number of items in the item bank (nrow(x)). Responses for items not administered at a given stage may be NA. If NULL and theta is provided, responses are simulated internally via simdat. If both theta and response are provided, response is used as-is and theta is retained for evaluation (bias/RMSE computation) only. Providing a pre-generated response matrix (e.g., from simdat(x, theta, D)) is useful when you want to compare different routing methods or scoring options on the same set of item responses.
D	A numeric scaling constant for the IRT model. Use D = 1 (default) for the logistic metric or D = 1.702 to approximate the normal ogive.
ini_mod	Either NULL (default) or a single positive integer index into the stage-1 module list (panel_info(route_map)$config[[1]]). When NULL, each examinee is independently and randomly assigned to one of the stage-1 modules with equal probability - appropriate for MST panels that use spiralling or booklet-based stage-1 assignment, or for panels with only a single stage-1 module (where random assignment is equivalent to a fixed assignment). When an integer is supplied, all examinees begin with that stage-1 module; use ini_mod = 1 to force all examinees to the first (and typically only) stage-1 module.
route_method	A character string specifying the adaptive routing method used to select the next module at each stage transition. "bmat" B-matching: selects the module whose average item location (difficulty) is closest to the current ability estimate. For dichotomous items the b parameter is used; for polytomous items (GRM, GPCM) the mean of the threshold / step parameters is used. "mfi" Maximum Fisher Information: selects the module that provides the highest test information function (TIF) value at the current ability estimate. NULL Cut-score routing: uses the cut scores supplied in cut_score to route examinees to the next module. Each examinee's current routing score is compared against the cut scores to determine which of the reachable next-stage modules to administer. Default is "bmat".
cut_score	A list of numeric vectors specifying the routing cut scores used when route_method = NULL. The list should have n.stage - 1 elements, one per stage transition. Element s contains the cut scores for routing from stage s to stage s+1. For example, in a 1-3-3 MST, cut_score = list(c(-0.5, 0.5), c(-0.6, 0.6)) routes examinees whose stage-1 score is below -0.5 to the easiest stage-2 module, between -0.5 and 0.5 to the medium module, and above 0.5 to the hardest module. Ignored when route_method is "bmat" or "mfi". Default is NULL.
route_score	A named list specifying the scoring method and options used to obtain intermediate ability estimates for routing decisions (applied after every stage except the last). Fields: method Character. One of "ML" (maximum likelihood), "WL" (weighted likelihood; Warm, 1989), "MLF" (maximum likelihood with fences; Han, 2016), "MAP" (maximum a posteriori), "EAP" (expected a posteriori; Bock & Mislevy, 1982), "EAP.SUM" (EAP summed scoring; Thissen et al., 1995), or "INV.TCC" (inverse test characteristic curve scoring; Lim et al., 2021). For "EAP.SUM" and "INV.TCC", a sum-score-to-theta lookup table is pre-computed once per module before the simulation loop; routing theta is then obtained by a single named-vector lookup, making the approach efficient for large N. Default: "ML". range Numeric vector of length 2: lower and upper bounds for the ability scale. Default: c(-5, 5). norm.prior Numeric vector of length 2: mean and SD of the normal prior for MAP/EAP. Default: c(0, 1). nquad Integer: number of quadrature points for EAP. Default: 41. tol Numeric: convergence tolerance for iterative methods. Default: 1e-4. max.iter Integer: maximum iterations. Default: 100. fence.a Numeric: discrimination parameter of fence items for MLF. Default: 3.0. fence.b Numeric vector of length 2 or NULL: difficulty parameters of the lower and upper fence items. If NULL, range is used. Default: NULL. intpol Logical: enable linear interpolation for "INV.TCC". Default: TRUE. range.tcc Numeric vector of length 2: theta search range for "INV.TCC". Default: c(-7, 7). max.it Integer: maximum bisection iterations for "INV.TCC". Default: 500L. Unspecified fields take their defaults; the full default is list(method = "ML", range = c(-5, 5), norm.prior = c(0, 1), nquad = 41L, tol = 1e-4, max.iter = 100L, fence.a = 3.0, fence.b = NULL, intpol = TRUE, range.tcc = c(-7, 7), max.it = 500L). Example: route_score = list(method = "INV.TCC", range.tcc = c(-5, 5)).
final_score	A named list specifying the scoring method and options for the final ability estimate (applied to all items accumulated across the entire administered pathway). Supports all fields in route_score plus: method Character. Supports all route_score methods, plus "EAP.SUM" (EAP summed scoring; Thissen et al., 1995) and "INV.TCC" (inverse test characteristic curve scoring; Lim et al., 2021). Default: "ML". intpol Logical: enable linear interpolation for "INV.TCC". Default: TRUE. range.tcc Numeric vector of length 2: theta search range for "INV.TCC". Default: c(-7, 7). max.it Integer: maximum bisection iterations for "INV.TCC". Default: 500. The full default is equivalent to list(method = "ML", range = c(-5, 5), norm.prior = c(0, 1), nquad = 41L, tol = 1e-4, max.iter = 100L, fence.a = 3.0, fence.b = NULL, intpol = TRUE, range.tcc = c(-7, 7), max.it = 500L).
se	Logical. Whether to compute and return the standard error of the final ability estimate. SE is computed for all final_score methods: for "ML", "WL", "MLF", "MAP", and "EAP" via the Fisher information / posterior variance; for "EAP.SUM" and "INV.TCC" via the posterior standard deviation stored in the pre-computed lookup table. Returns NA when se = FALSE, when all responses are missing, or when the sum score falls outside the estimable range (e.g. a perfect or zero score beyond range.tcc). Default is TRUE.
missing	A scalar specifying how missing (not-administered) responses are represented in response. Default is NA.
verbose	Logical. If TRUE, prints progress messages to the console during the simulation (approximately every 10\ Default is TRUE.
return_full_resp	Logical. If TRUE, the full.resp element of the returned list is an N x J integer matrix, where N is the number of examinees and J = nrow(x) (all items in the item bank). Cell [i, j] contains examinee i's response to item j if that item was administered, or NA if it was not. Row names are "examinee.1", "examinee.2", etc.; column names are the item IDs from x$id. If FALSE (default), full.resp is NULL. Because this matrix can be large for large panels and many examinees, it is not populated by default.

Details

Panel structure: The function first calls panel_info to identify stages, module membership, and valid pathways from route_map. Module m's items are those where module[, m] == 1.

Response simulation: If response is NULL, item responses are simulated from theta using simdat.

Stage-1 module assignment: When ini_mod = NULL (default), each examinee is independently assigned to a stage-1 module by uniform random sampling from panel_info(route_map)$config[[1]]. For panels with a single stage-1 module (the most common design), this is equivalent to a fixed assignment. For panels with multiple stage-1 modules (e.g., a 2-3-3 or 3-2-3 design), random assignment approximates the spiralling or booklet-based allocation used in operational testing. When ini_mod is an integer, all examinees start at the specified stage-1 module.

Routing (stages 1 to n.stage - 1): After each non-final stage, an intermediate ability estimate is obtained from the current stage's responses using the route_score method. This estimate is used to select the next-stage module:

• "bmat": the reachable module (via route_map) whose mean item location is closest to the routing estimate is selected. For dichotomous items, the location is the b parameter; for GRM/GPCM items, it is the mean of the threshold / step parameters.

• "mfi": the reachable module with the highest test information at the routing estimate is selected.

• NULL: the routing estimate is compared against the cut scores in cut_score[[s]] (for the transition from stage s to stage s+1) to assign a rank, and the rank-th reachable module (ordered by module index) is administered. If the rank exceeds the number of reachable modules, the last module is used.

Final scoring: Responses from all administered stages are concatenated, and the final ability estimate is computed using the final_score method.

Value

An object of class "run_mst", which is a named list containing:

call

The matched function call.

est.theta

A numeric vector of length N containing the final ability estimate for each examinee.

se.theta

A numeric vector of length N containing the standard error of the final ability estimate. All seven final_score methods return SE; for "EAP.SUM" and "INV.TCC", SE is the posterior standard deviation from the pre-computed lookup table. NA when se = FALSE, when all responses are missing, or when the sum score falls outside the estimable range.

theta.route

An N x n.stage numeric matrix of ability estimates. Columns 1 to n.stage - 1 are the intermediate routing estimates; column n.stage is the final estimate (equal to est.theta).

path

An N x n.stage integer matrix of administered module indices.

true.theta

The theta argument (true abilities), or NULL if not provided.

panel

The output of panel_info(route_map).

D

The scaling constant used.

route.method

The route_method argument.

route.score

The fully-merged route_score list (with all defaults applied).

final.score

The fully-merged final_score list (with all defaults applied).

N

Number of examinees.

full.resp

Always present in the returned list. When return_full_resp = TRUE, an N x J integer matrix of all item responses, where N is the number of examinees and J = nrow(x) (total items in the item bank). Row names are "examinee.1", "examinee.2", etc.; column names are the item IDs from x$id. NA for items not administered to a given examinee. When return_full_resp = FALSE (default), this element is NULL.

References

Bock, R. D., & Mislevy, R. J. (1982). Adaptive EAP estimation of ability in a microcomputer environment. Applied Psychological Measurement, 6(4), 431-444. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1177/014662168200600405")}

Han, K. T. (2016). Maximum likelihood score estimation method with fences for short-length tests and computerized adaptive tests. Applied Psychological Measurement, 40(4), 289-301. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1177/0146621616631317")}

Lim, H., Davey, T., & Wells, C. S. (2021). A recursion-based analytical approach to evaluate the performance of MST. Journal of Educational Measurement, 58(2), 154-178. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1111/jedm.12276")}

Magis, D., Yan, D., & von Davier, A. A. (2017). Computerized adaptive and multistage testing with R: Using packages catR and mstR. Springer. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1007/978-3-319-69218-0")}

Thissen, D., Pommerich, M., Billeaud, K., & Williams, V. S. L. (1995). Item response theory for scores on tests including polytomous items with ordered responses. Applied Psychological Measurement, 19(1), 39-49. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1177/014662169501900105")}

Warm, T. A. (1989). Weighted likelihood estimation of ability in item response theory. Psychometrika, 54(3), 427-450. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1007/BF02294627")}

See Also

panel_info, reval_mst, est_score, simdat, simMST

Examples

## ---------------------------------------------------------
## Example 1: bmat routing with ML scoring (default)
## Using the simMST 1-3-3 panel (7 modules, 8 items each)
## ---------------------------------------------------------
# Load item bank metadata
x        <- simMST$item_bank
module   <- simMST$module
route_map <- simMST$route_map

# True abilities for 500 examinees
set.seed(42)
theta_true <- rnorm(500, mean = 0, sd = 1)

# Run MST simulation with bmat routing and ML scoring
result_bmat <- run_mst(
  x          = x,
  route_map  = route_map,
  module     = module,
  theta      = theta_true,
  D          = 1.702,
  route_method = "bmat",
  route_score  = list(method = "ML", range = c(-4, 4)),
  final_score  = list(method = "ML", range = c(-4, 4)),
  se = TRUE
)

# Print simulation summary
print(result_bmat)

# Final ability estimates
head(result_bmat$est.theta)

# Module pathway taken by each examinee
head(result_bmat$path)

## ---------------------------------------------------------
## Example 2: Cut-score routing with EAP routing and ML final scoring
## EAP keeps the routing estimate finite even from a short Stage 1
## module (e.g., a perfect or zero score), while ML avoids the shrinkage
## that EAP's prior would otherwise introduce in the reported final score.
## ---------------------------------------------------------
cut_score <- simMST$cut_score   # list(c(-0.45, 0.46), c(-0.42, 0.45))

result_cut <- run_mst(
  x            = x,
  route_map    = route_map,
  module       = module,
  theta        = theta_true,
  D            = 1.702,
  route_method = NULL,
  cut_score    = cut_score,
  route_score  = list(method = "EAP", norm.prior = c(0, 1), nquad = 41),
  final_score  = list(method = "ML", range = c(-4, 4)),
  se = TRUE
)
print(result_cut)

## ---------------------------------------------------------
## Example 3: Using a pre-generated response matrix
## (useful when comparing multiple routing methods on the
##  same simulated responses, or when using real observed data)
## ---------------------------------------------------------
# Generate responses once from theta using simdat()
set.seed(123)
theta_true2 <- rnorm(500, mean = 0, sd = 1)
resp_matrix <- simdat(x = x, theta = theta_true2, D = 1.702)
# resp_matrix is an N x J matrix (N examinees x J total items in bank)

# Now run run_mst() twice with the SAME responses but different routing methods
# Method A: bmat routing
result_A <- run_mst(
  x            = x,
  route_map    = route_map,
  module       = module,
  theta        = theta_true2,   # kept for bias/RMSE evaluation
  response     = resp_matrix,   # pre-generated N x J response matrix
  D            = 1.702,
  route_method = "bmat",
  route_score  = list(method = "EAP", norm.prior = c(0, 1), nquad = 41),
  final_score  = list(method = "ML", range = c(-4, 4)),
  se = FALSE
)

# Method B: cut-score routing using the same responses
result_B <- run_mst(
  x            = x,
  route_map    = route_map,
  module       = module,
  theta        = theta_true2,
  response     = resp_matrix,   # identical response matrix -> fair comparison
  D            = 1.702,
  route_method = NULL,
  cut_score    = simMST$cut_score,
  route_score  = list(method = "EAP", norm.prior = c(0, 1), nquad = 41),
  final_score  = list(method = "ML", range = c(-4, 4)),
  se = FALSE
)

# Compare RMSE of the two routing strategies
rmse_A <- sqrt(mean((result_A$est.theta - theta_true2)^2))
rmse_B <- sqrt(mean((result_B$est.theta - theta_true2)^2))
cat(sprintf("RMSE (bmat): %.4f\n", rmse_A))
cat(sprintf("RMSE (cut-score): %.4f\n", rmse_B))

## ---------------------------------------------------------
## Example 4: TIF-based cut scores via find_cut()
## find_cut() computes principled cut scores from TIF crossings.
## This prevents anomalous routing under MFI-like heuristics.
## ---------------------------------------------------------
# Derive cut scores from TIF crossings between adjacent modules
cut_result <- find_cut(
  x         = x,
  module    = module,
  route_map = route_map
)

# Inspect the derived cut scores
print(cut_result)

# Use the derived cut scores for routing (route_method = NULL)
result_findcut <- run_mst(
  x            = x,
  route_map    = route_map,
  module       = module,
  theta        = theta_true,
  D            = 1.702,
  route_method = NULL,
  cut_score    = cut_result$cut_score,
  route_score  = list(method = "EAP", norm.prior = c(0, 1), nquad = 41),
  final_score  = list(method = "ML", range = c(-4, 4)),
  se = TRUE
)
print(result_findcut)

## Visualise the TIF-based cut scores used for routing
plot(cut_result)

## ---------------------------------------------------------
## Example 5: INV.TCC routing and final scoring
## (pre-computed lookup tables; efficient for large N)
## ---------------------------------------------------------
# INV.TCC can be used as both the routing and final scoring method.
# Before the examinee loop, run_mst() pre-computes a sum_score -> theta
# lookup table for every module (routing) and for every unique complete
# pathway (final scoring), consistent with the reval_mst() approach.
# SE is now returned correctly for INV.TCC (posterior SD from the table).

result_inv <- run_mst(
  x            = x,
  route_map    = route_map,
  module       = module,
  theta        = theta_true,
  D            = 1.702,
  route_method = "bmat",
  route_score  = list(method = "INV.TCC", range.tcc = c(-7, 7)),
  final_score  = list(method = "INV.TCC", range.tcc = c(-7, 7)),
  se           = TRUE
)
print(result_inv)

# SE is now populated (posterior SD from the pre-computed table)
head(result_inv$se.theta)

# Compare RMSE: INV.TCC routing vs. ML routing (result_bmat from Example 1)
rmse_inv <- sqrt(mean((result_inv$est.theta  - theta_true)^2, na.rm = TRUE))
rmse_ml  <- sqrt(mean((result_bmat$est.theta - theta_true)^2, na.rm = TRUE))
cat(sprintf("RMSE (INV.TCC routing + scoring): %.4f\n", rmse_inv))
cat(sprintf("RMSE (ML    routing + scoring):   %.4f\n", rmse_ml))

irtQ documentation built on July 27, 2026, 9:08 a.m.