Multistage-Adaptive Test (MST) Simulation

Simulates an MST administration for a set 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, with separate route_score and final_score list arguments and the route_map and 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
)

# S3 method for class 'run_mst'
print(x, digits = 3L, ...)

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 item metadata are checked and completed internally before the simulation.

route_map

A binary square matrix defining the MST transition structure. An entry of 1 in row i and 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)). Only the responses to the items of the administered modules are used, so the entries for the other items are ignored and may be NA. A missing response to an administered item (NA or the value given in missing) is excluded from the likelihood for "ML", "WL", "MLF", "MAP", and "EAP", and counted as 0 in the sum score for "EAP.SUM" and "INV.TCC". 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 scaling constant used in IRT models to make the logistic function closely approximate the normal ogive function. A value of 1.702 is commonly used for this purpose. Default is 1.

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. This is appropriate for panels that use spiraling or booklet-based stage-1 assignment, and it is equivalent to a fixed assignment for panels with a single stage-1 module. 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, from $-0.5$ up to but not including $0.5$ to the medium module, and $0.5$ or above to the hardest module (the modules numbered from easiest to hardest). Each element must contain finite values in strictly ascending order and have one fewer value than the number of modules in stage s+1; cut score k separates the k-th and the (k+1)-th module of that stage in the order of the module indices. A NULL element is read as an empty vector for a stage with a single module. The output element cut_score of find_cut can be passed directly. 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; see est_score). For "EAP.SUM" and "INV.TCC", a sum-score-to-theta lookup table is pre-computed once per partial pathway (the modules of stages 1 to s) 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 of the ability estimates for "ML", "WL", "MLF", and "MAP"; also the default fence locations for "MLF". Default: c(-5, 5).
norm.prior: Numeric vector of length 2: mean and SD of the normal prior for "MAP", "EAP", and "EAP.SUM". Default: c(0, 1).
nquad: Integer: number of quadrature points for "EAP" and "EAP.SUM". Default: 41.
tol: A single positive number: convergence tolerance of the Fisher scoring iterations for "ML", "WL", "MLF", and "MAP", and of the bisection search for "INV.TCC". Default: 1e-4.
max.iter: A single positive whole number: maximum number of Fisher scoring iterations for "ML", "WL", "MLF", and "MAP". 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: ability estimates assigned to the lowest and the maximum possible sum scores for "INV.TCC" when intpol = TRUE (see est_score). Default: c(-7, 7).
max.it: Integer: maximum bisection iterations for "INV.TCC". Default: 500L.

As in est_score, the function stops with an error when tol or max.iter is not valid for the chosen method. 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). It takes the same fields, methods, and defaults as route_score. 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 return the standard error of the final ability estimate. The standard errors are computed as in est_score: from the test information for "ML", "WL", "MLF", and "MAP", as the posterior standard deviation for "EAP" and "EAP.SUM", and, for "INV.TCC", as the standard deviation of the inverse TCC estimates over the conditional sum score distribution at the estimate. For "ML", "WL", "MLF", and "MAP", the standard error is set to 99.9999 when the estimate equals a bound of range. The standard error is NA when se = FALSE, when all responses are missing, or, for "INV.TCC", when the sum score has no ability estimate (for example, with intpol = FALSE). Default is TRUE.

missing

A scalar specifying the value that marks a missing response in response. It is used only when response is supplied. 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.

digits

The number of decimal places used for the summary of the ability estimates and for the bias and RMSE. Default is 3.

...

Additional arguments, currently not used.

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 errors of the final ability estimates (see the se argument).
theta.route: An N x n.stage numeric matrix of ability estimates. Column s (s < n.stage) is the cumulative routing estimate after stage s, computed from the responses to all modules administered in stages 1 to s; column n.stage is the final estimate (equal to est.theta). Column names are "stage.1", "stage.2", etc.
path: An N x n.stage integer matrix of administered module indices. Column names are "stage.1", "stage.2", etc.
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.

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 spiraling 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 responses to all modules administered so far 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). The cut scores are tied to the modules of stage s+1 in order of module index: cut score k separates the k-th and the (k+1)-th module of the stage. When only some modules of the stage can be reached from the current module, only the cut scores that separate the reachable modules are used, as in reval_mst, and the module whose interval contains the estimate is administered. An estimate equal to a cut score is assigned to the higher module. When every module of the next stage is reachable, all cut scores are used.

Final scoring: Responses from all administered stages are concatenated, and the final ability estimate is computed using the final_score method. With "EAP.SUM" or "INV.TCC", a missing response is counted as 0 in the sum score, as in est_score, but no warning is issued. The same rule applies to the routing estimates. When all responses up to a non-final stage are missing, the routing estimate is 0 after stage 1 and the previous routing estimate after a later stage. An examinee with no observed response to the administered items receives NA as the final estimate, and one warning reports the number of such examinees.

Relation to reval_mst(): With route_method = NULL, a cut_score list, and route_score = list(method = "INV.TCC"), the routing rule is the one that reval_mst evaluates analytically: cut scores are applied to the inverse TCC estimate of the cumulative sum score. With final_score = list(method = "INV.TCC") as well, the simulated bias and conditional standard error of measurement approach the reval_mst() values as the number of simulated examinees grows.

Methods (by generic)

print(run_mst): Print method for an object of class run_mst (given as x). It prints the panel structure, the routing and scoring methods, a summary of the final ability estimates, the bias and RMSE when the true abilities are available, and the module frequencies by stage.

References

Bock, R. D., & Mislevy, R. J. (1982). Adaptive EAP estimation of ability in a microcomputer environment. Applied Psychological Measurement, 6(4), 431-444. 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. 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. 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. 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. doi:10.1177/014662169501900105 .

Warm, T. A. (1989). Weighted likelihood estimation of ability in item response theory. Psychometrika, 54(3), 427-450. doi:10.1007/BF02294627 .

Examples

# \donttest{
## ---------------------------------------------------------
## 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
)
#> [run_mst] Panel: 3 stages, 7 modules (1-3-3 design)
#> [run_mst] Simulating responses for 500 examinees...
#> [run_mst] Scoring 500 examinees...
#>   Processing examinee 1 / 500 ...
#>   Processing examinee 50 / 500 ...
#>   Processing examinee 100 / 500 ...
#>   Processing examinee 150 / 500 ...
#>   Processing examinee 200 / 500 ...
#>   Processing examinee 250 / 500 ...
#>   Processing examinee 300 / 500 ...
#>   Processing examinee 350 / 500 ...
#>   Processing examinee 400 / 500 ...
#>   Processing examinee 450 / 500 ...
#>   Processing examinee 500 / 500 ...
#> [run_mst] Done.

# Print simulation summary
print(result_bmat)
#> 
#> Call:
#> 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) 
#> 
#> MST Simulation Results
#> ======================================== 
#> 
#> Panel structure:
#>   Stages              : 3
#>   Modules per stage   : 1 - 3 - 3
#>   Valid pathways      : 7
#>   Routing method      : bmat (b-matching)
#> 
#> Number of examinees : 500
#> 
#> Ability estimation:
#>   Routing method      : ML
#>   Final method        : ML
#> 
#> Final ability estimates (est.theta):
#>   Mean : -0.011
#>   SD   : 1.077
#>   Min  : -4.000
#>   Max  : 4.000
#> 
#> Estimation accuracy (est.theta - true.theta):
#>   Bias : 0.019
#>   RMSE : 0.319
#> 
#> Module frequency by stage:
#>   Stage 1: Module 1: 500 (100.0%)
#>   Stage 2: Module 2: 177 (35.4%),  Module 3: 106 (21.2%),  Module 4: 217 (43.4%)
#>   Stage 3: Module 5: 185 (37.0%),  Module 6: 105 (21.0%),  Module 7: 210 (42.0%)
#> 

# Final ability estimates
head(result_bmat$est.theta)
#> [1]  1.6882857 -0.8039112  0.6168219  0.5994916  0.4140971  0.2068442

# Module pathway taken by each examinee
head(result_bmat$path)
#>      stage.1 stage.2 stage.3
#> [1,]       1       4       7
#> [2,]       1       2       5
#> [3,]       1       4       7
#> [4,]       1       4       7
#> [5,]       1       3       7
#> [6,]       1       4       6

## ---------------------------------------------------------
## 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
)
#> [run_mst] Panel: 3 stages, 7 modules (1-3-3 design)
#> [run_mst] Simulating responses for 500 examinees...
#> [run_mst] Scoring 500 examinees...
#>   Processing examinee 1 / 500 ...
#>   Processing examinee 50 / 500 ...
#>   Processing examinee 100 / 500 ...
#>   Processing examinee 150 / 500 ...
#>   Processing examinee 200 / 500 ...
#>   Processing examinee 250 / 500 ...
#>   Processing examinee 300 / 500 ...
#>   Processing examinee 350 / 500 ...
#>   Processing examinee 400 / 500 ...
#>   Processing examinee 450 / 500 ...
#>   Processing examinee 500 / 500 ...
#> [run_mst] Done.
print(result_cut)
#> 
#> Call:
#> 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) 
#> 
#> MST Simulation Results
#> ======================================== 
#> 
#> Panel structure:
#>   Stages              : 3
#>   Modules per stage   : 1 - 3 - 3
#>   Valid pathways      : 7
#>   Routing method      : cut-score based
#> 
#> Number of examinees : 500
#> 
#> Ability estimation:
#>   Routing method      : EAP
#>   Final method        : ML
#> 
#> Final ability estimates (est.theta):
#>   Mean : -0.030
#>   SD   : 1.002
#>   Min  : -4.000
#>   Max  : 4.000
#> 
#> Estimation accuracy (est.theta - true.theta):
#>   Bias : 0.000
#>   RMSE : 0.292
#> 
#> Module frequency by stage:
#>   Stage 1: Module 1: 500 (100.0%)
#>   Stage 2: Module 2: 142 (28.4%),  Module 3: 219 (43.8%),  Module 4: 139 (27.8%)
#>   Stage 3: Module 5: 171 (34.2%),  Module 6: 189 (37.8%),  Module 7: 140 (28.0%)
#> 

## ---------------------------------------------------------
## 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
)
#> [run_mst] Panel: 3 stages, 7 modules (1-3-3 design)
#> [run_mst] Scoring 500 examinees...
#>   Processing examinee 1 / 500 ...
#>   Processing examinee 50 / 500 ...
#>   Processing examinee 100 / 500 ...
#>   Processing examinee 150 / 500 ...
#>   Processing examinee 200 / 500 ...
#>   Processing examinee 250 / 500 ...
#>   Processing examinee 300 / 500 ...
#>   Processing examinee 350 / 500 ...
#>   Processing examinee 400 / 500 ...
#>   Processing examinee 450 / 500 ...
#>   Processing examinee 500 / 500 ...
#> [run_mst] Done.

# 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
)
#> [run_mst] Panel: 3 stages, 7 modules (1-3-3 design)
#> [run_mst] Scoring 500 examinees...
#>   Processing examinee 1 / 500 ...
#>   Processing examinee 50 / 500 ...
#>   Processing examinee 100 / 500 ...
#>   Processing examinee 150 / 500 ...
#>   Processing examinee 200 / 500 ...
#>   Processing examinee 250 / 500 ...
#>   Processing examinee 300 / 500 ...
#>   Processing examinee 350 / 500 ...
#>   Processing examinee 400 / 500 ...
#>   Processing examinee 450 / 500 ...
#>   Processing examinee 500 / 500 ...
#> [run_mst] Done.

# 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))
#> RMSE (bmat): 0.3068
cat(sprintf("RMSE (cut-score): %.4f\n", rmse_B))
#> RMSE (cut-score): 0.3050

## ---------------------------------------------------------
## Example 4: TIF-based cut scores via find_cut()
## find_cut() derives the cut scores from the TIF crossings of adjacent
## modules; cut-score routing with these values avoids the path
## reversals that MFI routing can produce.
## ---------------------------------------------------------
# Derive cut scores from TIF crossings between adjacent modules
cut_result <- find_cut(
  x         = x,
  module    = module,
  route_map = route_map,
  D         = 1.702
)

# Inspect the derived cut scores
print(cut_result)
#> MST TIF-Crossing Cut Score Results
#> ==================================================
#> 
#> stage.2
#> -----------------------------------
#>   Modules (index order)    : 2, 3, 4
#>   Modules (difficulty order): 2, 3, 4
#>   Mean item locations      : -0.87, -0.0574, 0.7858
#>   Selected cut score(s)    : -0.4451, 0.4589
#>   Pair details:
#>     [mod2_vs_mod3]
#>       Proper crossing(s)    : 1  [theta = -0.4451]
#>       Anomalous crossing(s) : 1  [theta = 3.3893]  (excluded)
#>       Selected cut score    : -0.4451
#>     [mod3_vs_mod4]
#>       Proper crossing(s)    : 1  [theta = 0.4589]
#>       Selected cut score    : 0.4589
#> 
#> stage.3
#> -----------------------------------
#>   Modules (index order)    : 5, 6, 7
#>   Modules (difficulty order): 5, 6, 7
#>   Mean item locations      : -0.5586, -0.042, 0.5066
#>   Selected cut score(s)    : -0.4171, 0.4532
#>   Pair details:
#>     [mod5_vs_mod6]
#>       Proper crossing(s)    : 1  [theta = -0.4171]
#>       Anomalous crossing(s) : 1  [theta = 1.617]  (excluded)
#>       Selected cut score    : -0.4171
#>     [mod6_vs_mod7]
#>       Proper crossing(s)    : 1  [theta = 0.4532]
#>       Anomalous crossing(s) : 1  [theta = -1.9874]  (excluded)
#>       Selected cut score    : 0.4532
#> ==================================================
#> Usage: run_mst(..., route_method = NULL, cut_score = <result>$cut_score)

# 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
)
#> [run_mst] Panel: 3 stages, 7 modules (1-3-3 design)
#> [run_mst] Simulating responses for 500 examinees...
#> [run_mst] Scoring 500 examinees...
#>   Processing examinee 1 / 500 ...
#>   Processing examinee 50 / 500 ...
#>   Processing examinee 100 / 500 ...
#>   Processing examinee 150 / 500 ...
#>   Processing examinee 200 / 500 ...
#>   Processing examinee 250 / 500 ...
#>   Processing examinee 300 / 500 ...
#>   Processing examinee 350 / 500 ...
#>   Processing examinee 400 / 500 ...
#>   Processing examinee 450 / 500 ...
#>   Processing examinee 500 / 500 ...
#> [run_mst] Done.
print(result_findcut)
#> 
#> Call:
#> 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) 
#> 
#> MST Simulation Results
#> ======================================== 
#> 
#> Panel structure:
#>   Stages              : 3
#>   Modules per stage   : 1 - 3 - 3
#>   Valid pathways      : 7
#>   Routing method      : cut-score based
#> 
#> Number of examinees : 500
#> 
#> Ability estimation:
#>   Routing method      : EAP
#>   Final method        : ML
#> 
#> Final ability estimates (est.theta):
#>   Mean : 0.000
#>   SD   : 1.058
#>   Min  : -4.000
#>   Max  : 4.000
#> 
#> Estimation accuracy (est.theta - true.theta):
#>   Bias : 0.030
#>   RMSE : 0.340
#> 
#> Module frequency by stage:
#>   Stage 1: Module 1: 500 (100.0%)
#>   Stage 2: Module 2: 136 (27.2%),  Module 3: 225 (45.0%),  Module 4: 139 (27.8%)
#>   Stage 3: Module 5: 170 (34.0%),  Module 6: 169 (33.8%),  Module 7: 161 (32.2%)
#> 

## Visualize 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 to theta
# lookup table for every partial pathway (routing) and every complete
# pathway (final scoring). This example uses bmat routing; the design
# that reval_mst() evaluates analytically uses cut-score routing
# (route_method = NULL) with INV.TCC routing and final scores.
# SE is taken from the pre-computed lookup 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
)
#> [run_mst] Panel: 3 stages, 7 modules (1-3-3 design)
#> [run_mst] Simulating responses for 500 examinees...
#> [run_mst] Scoring 500 examinees...
#>   Processing examinee 1 / 500 ...
#>   Processing examinee 50 / 500 ...
#>   Processing examinee 100 / 500 ...
#>   Processing examinee 150 / 500 ...
#>   Processing examinee 200 / 500 ...
#>   Processing examinee 250 / 500 ...
#>   Processing examinee 300 / 500 ...
#>   Processing examinee 350 / 500 ...
#>   Processing examinee 400 / 500 ...
#>   Processing examinee 450 / 500 ...
#>   Processing examinee 500 / 500 ...
#> [run_mst] Done.
print(result_inv)
#> 
#> Call:
#> 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) 
#> 
#> MST Simulation Results
#> ======================================== 
#> 
#> Panel structure:
#>   Stages              : 3
#>   Modules per stage   : 1 - 3 - 3
#>   Valid pathways      : 7
#>   Routing method      : bmat (b-matching)
#> 
#> Number of examinees : 500
#> 
#> Ability estimation:
#>   Routing method      : INV.TCC
#>   Final method        : INV.TCC
#> 
#> Final ability estimates (est.theta):
#>   Mean : -0.021
#>   SD   : 1.267
#>   Min  : -7.000
#>   Max  : 7.000
#> 
#> Estimation accuracy (est.theta - true.theta):
#>   Bias : 0.009
#>   RMSE : 0.599
#> 
#> Module frequency by stage:
#>   Stage 1: Module 1: 500 (100.0%)
#>   Stage 2: Module 2: 187 (37.4%),  Module 3: 84 (16.8%),  Module 4: 229 (45.8%)
#>   Stage 3: Module 5: 199 (39.8%),  Module 6: 94 (18.8%),  Module 7: 207 (41.4%)
#> 

# standard errors from the lookup table
head(result_inv$se.theta)
#> [1] 0.5213622 0.2798838 0.2850462 0.2797603 0.2756490 0.2739662

# Compare RMSE: INV.TCC routing and scoring vs. ML routing and scoring
# (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))
#> RMSE (INV.TCC routing + scoring): 0.5995
cat(sprintf("RMSE (ML    routing + scoring):   %.4f\n", rmse_ml))
#> RMSE (ML    routing + scoring):   0.3188
# }