diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index 2fd552b..b3872e3 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -36,25 +36,42 @@ jobs: - uses: r-lib/actions/setup-pandoc@d3c5be51b12e724e68f33216ca3c148b66d5f0b6 # v2 - - uses: r-lib/actions/setup-r@6f6e5bc62fba3a704f74e7ad7ef7676c5c6a2590 # v2 + - uses: r-lib/actions/setup-r@d3c5be51b12e724e68f33216ca3c148b66d5f0b6 # v2 with: r-version: ${{ matrix.config.r }} http-user-agent: ${{ matrix.config.http-user-agent }} use-public-rspm: true - - uses: r-lib/actions/setup-r-dependencies@6f6e5bc62fba3a704f74e7ad7ef7676c5c6a2590 # v2 + # Build stringfish before dependency setup because packages built during + # resolution can load it and fail on a prebuilt oneTBB ABI mismatch. + - name: Source-build stringfish before dependency resolution + run: >- + Rscript -e 'install.packages(c("RcppParallel", "stringfish"), type = "source", + repos = "https://cloud.r-project.org")' + + - uses: r-lib/actions/setup-r-dependencies@d3c5be51b12e724e68f33216ca3c148b66d5f0b6 # v2 with: + # Refresh the dependency cache after the macOS qs2/RcppParallel TBB ABI mismatch. + cache-version: '2' extra-packages: | any::rcmdcheck any::testthat needs: check + # Rebuild every installed package that links to RcppParallel against the + # same local oneTBB ABI before loading kaefa or its test dependencies. + - name: Rebuild RcppParallel-linked packages from source + run: >- + Rscript -e 'r <- "https://cloud.r-project.org"; install.packages("RcppParallel", type = "source", repos = r); + linked <- tools::dependsOnPkgs("RcppParallel", dependencies = "LinkingTo"); + if (length(linked)) install.packages(linked, type = "source", repos = r)' + - name: Run Zh formula regression tests run: | Rscript -e 'install.packages(".", repos = NULL, type = "source")' Rscript -e 'library(kaefa); testthat::test_file("tests/testthat/test-zh-misfit-decision-rule.R")' - - uses: r-lib/actions/check-r-package@6f6e5bc62fba3a704f74e7ad7ef7676c5c6a2590 # v2 + - uses: r-lib/actions/check-r-package@d3c5be51b12e724e68f33216ca3c148b66d5f0b6 # v2 with: error-on: '"error"' upload-snapshots: true diff --git a/.github/workflows/test-fast.yaml b/.github/workflows/test-fast.yaml index ea1cf69..32c4c54 100644 --- a/.github/workflows/test-fast.yaml +++ b/.github/workflows/test-fast.yaml @@ -16,15 +16,30 @@ jobs: steps: - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1 - - uses: r-lib/actions/setup-r@6f6e5bc62fba3a704f74e7ad7ef7676c5c6a2590 # v2 + - uses: r-lib/actions/setup-r@d3c5be51b12e724e68f33216ca3c148b66d5f0b6 # v2 with: use-public-rspm: true - - uses: r-lib/actions/setup-r-dependencies@6f6e5bc62fba3a704f74e7ad7ef7676c5c6a2590 # v2 + # Build stringfish before dependency setup because packages built during + # resolution can load it and fail on a prebuilt oneTBB ABI mismatch. + - name: Source-build stringfish before dependency resolution + run: >- + Rscript -e 'install.packages(c("RcppParallel", "stringfish"), type = "source", + repos = "https://cloud.r-project.org")' + + - uses: r-lib/actions/setup-r-dependencies@d3c5be51b12e724e68f33216ca3c148b66d5f0b6 # v2 with: extra-packages: any::testthat needs: check + # Rebuild every installed package that links to RcppParallel against the + # same local oneTBB ABI before loading kaefa or its test dependencies. + - name: Rebuild RcppParallel-linked packages from source + run: >- + Rscript -e 'r <- "https://cloud.r-project.org"; install.packages("RcppParallel", type = "source", repos = r); + linked <- tools::dependsOnPkgs("RcppParallel", dependencies = "LinkingTo"); + if (length(linked)) install.packages(linked, type = "source", repos = r)' + - name: Install kaefa package for fast tests run: R CMD INSTALL . @@ -38,4 +53,8 @@ jobs: reporter = reporter) testthat::test_file("tests/testthat/test-core-api-contract.R", reporter = reporter) + testthat::test_file("tests/testthat/test-aefa-parameter-recovery.R", + reporter = reporter) + testthat::test_file("tests/testthat/test-mixedmirt-parameter-recovery.R", + reporter = reporter) RSCRIPT diff --git a/.github/workflows/test-suite.yaml b/.github/workflows/test-suite.yaml index 439b9b7..0694152 100644 --- a/.github/workflows/test-suite.yaml +++ b/.github/workflows/test-suite.yaml @@ -1,7 +1,8 @@ name: test-suite -# Runs the complete testthat suite (26 files), which the previous CI setup -# skipped: R-CMD-check runs with --no-tests and test-fast runs only 3 files. +# Runs the complete testthat suite, which the previous CI setup skipped: +# R-CMD-check runs with --no-tests and test-fast runs a small productization +# subset. # On push/PR the suite's own CI guards (.skip_expensive_ci_calls, skip_on_ci) # keep expensive AEFA estimations out. The weekly scheduled run and manual # dispatch set RUN_FULL_AEFA_TESTS=1, which unlocks the estimations guarded @@ -32,15 +33,30 @@ jobs: steps: - uses: actions/checkout@3d3c42e5aac5ba805825da76410c181273ba90b1 # v7.0.1 - - uses: r-lib/actions/setup-r@6f6e5bc62fba3a704f74e7ad7ef7676c5c6a2590 # v2 + - uses: r-lib/actions/setup-r@d3c5be51b12e724e68f33216ca3c148b66d5f0b6 # v2 with: use-public-rspm: true - - uses: r-lib/actions/setup-r-dependencies@6f6e5bc62fba3a704f74e7ad7ef7676c5c6a2590 # v2 + # Build stringfish before dependency setup because packages built during + # resolution can load it and fail on a prebuilt oneTBB ABI mismatch. + - name: Source-build stringfish before dependency resolution + run: >- + Rscript -e 'install.packages(c("RcppParallel", "stringfish"), type = "source", + repos = "https://cloud.r-project.org")' + + - uses: r-lib/actions/setup-r-dependencies@d3c5be51b12e724e68f33216ca3c148b66d5f0b6 # v2 with: extra-packages: any::testthat needs: check + # Rebuild every installed package that links to RcppParallel against the + # same local oneTBB ABI before loading kaefa or its test dependencies. + - name: Rebuild RcppParallel-linked packages from source + run: >- + Rscript -e 'r <- "https://cloud.r-project.org"; install.packages("RcppParallel", type = "source", repos = r); + linked <- tools::dependsOnPkgs("RcppParallel", dependencies = "LinkingTo"); + if (length(linked)) install.packages(linked, type = "source", repos = r)' + - name: Install kaefa package run: R CMD INSTALL . diff --git a/AGENTS.md b/AGENTS.md index e83e8e3..71aaee8 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -11,6 +11,7 @@ Use this file as the local operating guide when modifying `kaefa`. - Core orchestration: `R/kaefa.R` - Estimation engine: `R/newEngine.R` - Shared helpers: `R/utils.R` +- Recovery protocol: `R/recovery.R` - Shiny app: `inst/shiny-app/app.R` - Tests: `tests/testthat/` - Local CI workflow definitions: `.github/workflows/` diff --git a/ARCHITECTURE.md b/ARCHITECTURE.md index 6ddd063..d20f1e9 100644 --- a/ARCHITECTURE.md +++ b/ARCHITECTURE.md @@ -1,6 +1,6 @@ # kaefa Architecture -Last updated: 2026-02-14 +Last updated: 2026-08-17 ## Purpose @@ -32,7 +32,11 @@ explicitly requires vendored source integration. - `R/kaefa.R`: public orchestration entry points and exported runtime behavior. - `R/newEngine.R`: candidate-model estimation engine used by the AEFA loop. -- `R/utils.R`: helper routines and shared utilities. +- `R/utils.R`: helper routines and shared utilities. `.mirt()` and + `.mixedmirt()` treat a missing (`NA`) second-order test as + non-convergence when `leniency` is false. +- `R/recovery.R`: internal true-parameter RMSE helpers used by the AEFA + and nested mixedmirt recovery protocols. Not a public API. - `inst/shiny-app/app.R`: bundled Shiny interface logic. - `inst/shiny-app/README.md`: Shiny usage and minimal UI configuration guide. - `tests/testthat/*.R`: functional, regression, and integration tests. @@ -58,6 +62,13 @@ explicitly requires vendored source integration. - PR merge requires review approval and resolved conversations. - Required checks include R-CMD-check matrix and dependency review. +- Buyer-facing recovery evidence is the five-repeat RMSE protocol in + `tests/testthat/test-aefa-parameter-recovery.R` and the nested mixedmirt + contracts in `tests/testthat/test-mixedmirt-parameter-recovery.R`. The + formula, alignment, and coverage-exclusion contracts run in `test-fast`. + Live `aefa()` and five-seed `.mixedmirt` recovery are gated by + `RUN_FULL_AEFA_TESTS=1`. One-seed `.mirt` and `.mixedmirt` recovery fits + run in the full suite. - If code scanning is enabled later, alerts can be tracked via GitHub code scanning APIs. diff --git a/CLAUDE.md b/CLAUDE.md index 9747fed..d3cf59b 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -57,6 +57,10 @@ testthat::test_file("tests/testthat/test-shiny-product-surface.R", reporter = reporter) testthat::test_file("tests/testthat/test-core-api-contract.R", reporter = reporter) +testthat::test_file("tests/testthat/test-aefa-parameter-recovery.R", + reporter = reporter) +testthat::test_file("tests/testthat/test-mixedmirt-parameter-recovery.R", + reporter = reporter) RSCRIPT # R CMD check as CI runs it (note: tests are skipped here) @@ -124,6 +128,8 @@ psychometrics (see AGENTS.md). - `R/newEngine.R` — `engineAEFA()`: the candidate-model estimation engine used by the aefa loop (MCMC/estimation controls such as NCYCLES, BURNIN, SEMCYCLES, rotation choices, etc.). +- `R/recovery.R` — internal true-parameter RMSE helpers for the AEFA and + nested mixedmirt recovery protocols (not exported). - `R/utils.R` — shared internals (`.mirt`/`.mixedmirt` wrappers around mirt, `.covdataClassifieder`, `.covdataFixedEffectComb`, `.exportParmsEME`) and the theta-prior utilities `fitThetaPrior()`, `testThetaPriorCalibration()`, diff --git a/R/recovery.R b/R/recovery.R new file mode 100644 index 0000000..c772c3e --- /dev/null +++ b/R/recovery.R @@ -0,0 +1,618 @@ +# Internal true-parameter recovery helpers for AEFA / IRT Monte Carlo evidence. +# These functions are not exported. Buyer-facing recovery is the test protocol +# in tests/testthat/test-aefa-parameter-recovery.R, +# tests/testthat/test-mixedmirt-parameter-recovery.R, and the provenance notes +# in docs/traceability/. + +#' Root-mean-square error for a recovered parameter vector +#' +#' @param estimated Numeric vector of recovered parameter values. +#' @param truth Numeric vector of known true parameter values. +#' @return A single finite RMSE value. +#' @noRd +.parameterRecoveryRMSE <- function(estimated, truth) { + estimated <- as.numeric(estimated) + truth <- as.numeric(truth) + if (length(estimated) != length(truth)) { + stop("Estimated and true vectors must have the same length.", call. = FALSE) + } + if (!length(estimated)) { + stop("Cannot compute RMSE from empty parameter vectors.", call. = FALSE) + } + if (any(!is.finite(estimated)) || any(!is.finite(truth))) { + stop("RMSE requires finite estimated and true values.", call. = FALSE) + } + sqrt(mean((estimated - truth)^2)) +} + +# R data.frames never keep rownames as NULL: `rownames(x) <- NULL` resets them +# to the sequential defaults "1", "2", .... Treat those as unnamed so callers +# cannot silently align on positional indices. + +#' Item names used for recovery alignment +#' +#' @param x A matrix or data.frame of item parameters. +#' @return Character item names, or `NULL` when names are missing or positional. +#' @noRd +.irtItemNames <- function(x) { + rn <- rownames(x) + if (is.null(rn) || identical(rn, as.character(seq_len(nrow(x))))) { + return(NULL) + } + rn +} + +#' Align estimated and true IRT item tables by identical item names +#' +#' @param estimated Matrix or data.frame of estimated item parameters. +#' @param truth Matrix or data.frame of known true item parameters. +#' @param columns Character vector of required recovery columns. +#' @return A list with aligned `estimated`, `truth`, and `items`. +#' @noRd +.alignIrtItemParameters <- function(estimated, truth, columns = c("a", "b")) { + if (!is.data.frame(estimated) && !is.matrix(estimated)) { + stop("estimated must be a matrix or data.frame of item parameters.", call. = FALSE) + } + if (!is.data.frame(truth) && !is.matrix(truth)) { + stop("truth must be a matrix or data.frame of item parameters.", call. = FALSE) + } + estimated <- as.data.frame(estimated, stringsAsFactors = FALSE) + truth <- as.data.frame(truth, stringsAsFactors = FALSE) + est_names <- .irtItemNames(estimated) + true_names <- .irtItemNames(truth) + if (is.null(est_names) || is.null(true_names)) { + stop( + "Estimated and true parameter tables must have item names as row names.", + call. = FALSE + ) + } + missing_estimated <- setdiff(columns, colnames(estimated)) + missing_truth <- setdiff(columns, colnames(truth)) + if (length(missing_estimated) || length(missing_truth)) { + stop( + "Missing recovery columns: ", + paste(unique(c(missing_estimated, missing_truth)), collapse = ", "), + call. = FALSE + ) + } + missing_estimated_items <- setdiff(true_names, est_names) + missing_truth_items <- setdiff(est_names, true_names) + if (length(missing_estimated_items) || length(missing_truth_items)) { + stop( + "Estimated and true parameter tables must contain the same item names.", + call. = FALSE + ) + } + shared <- est_names + list( + estimated = estimated[shared, columns, drop = FALSE], + truth = truth[shared, columns, drop = FALSE], + items = shared + ) +} + +#' Summarise exactly five complete recovery repeats +#' +#' @param rmse_by_repeat A data.frame with `seed`, `parameter`, and `rmse`. +#' @return A list with `n_repeats`, `seeds`, `per_run`, and `summary`. +#' @noRd +.summariseRecoveryRepeats <- function(rmse_by_repeat) { + if (!is.data.frame(rmse_by_repeat)) { + stop("rmse_by_repeat must be a data.frame.", call. = FALSE) + } + required <- c("seed", "parameter", "rmse") + missing <- setdiff(required, names(rmse_by_repeat)) + if (length(missing)) { + stop( + "rmse_by_repeat must contain columns: ", + paste(required, collapse = ", "), + call. = FALSE + ) + } + seeds <- unique(rmse_by_repeat$seed) + if (length(seeds) != 5L) { + stop("Recovery protocol requires exactly 5 repeats.", call. = FALSE) + } + if (anyNA(rmse_by_repeat$seed) || + anyNA(rmse_by_repeat$parameter) || + anyNA(rmse_by_repeat$rmse)) { + stop("Recovery repeats cannot contain missing values.", call. = FALSE) + } + parameters <- unique(as.character(rmse_by_repeat$parameter)) + repeat_counts <- table( + as.character(rmse_by_repeat$parameter), + rmse_by_repeat$seed + ) + if (!length(repeat_counts) || any(repeat_counts != 1L)) { + stop( + "Each parameter must have exactly one RMSE value for each recovery seed.", + call. = FALSE + ) + } + summary_rows <- lapply(parameters, function(parameter_name) { + values <- rmse_by_repeat$rmse[as.character(rmse_by_repeat$parameter) == parameter_name] + data.frame( + parameter = parameter_name, + n_repeats = length(values), + mean_rmse = mean(values), + sd_rmse = stats::sd(values), + stringsAsFactors = FALSE + ) + }) + list( + n_repeats = 5L, + seeds = seeds, + per_run = rmse_by_repeat[order(rmse_by_repeat$seed, rmse_by_repeat$parameter), ], + summary = do.call(rbind, summary_rows) + ) +} + +#' Extract IRT item parameters from an aefa history or mirt fit +#' +#' `MixedClass` is accepted so the fail-closed contract can be tested: supported +#' `mirt` versions do not implement `IRTpars` / `simplify` for `coef,MixedClass`. +#' That path therefore returns no `$items` table and stops. Use +#' `.extractMixedmirtIrtItems()` for mixedmirt recovery. +#' +#' @param fit An `aefa` history, `SingleGroupClass`, or `MixedClass` object. +#' @return A data.frame of item parameters with item names as row names. +#' @noRd +.extractAefaIrtItems <- function(fit) { + if (inherits(fit, "aefa") || (is.list(fit) && !is.null(fit$estModelTrials))) { + trials <- fit$estModelTrials + if (!length(trials)) { + stop("aefa history has no estimated model trials.", call. = FALSE) + } + fit <- trials[[length(trials)]] + } + if (!methods::is(fit, "SingleGroupClass") && !methods::is(fit, "MixedClass")) { + stop("Recovery extraction requires an aefa history or a mirt model.", call. = FALSE) + } + items <- mirt::coef(fit, IRTpars = TRUE, simplify = TRUE)$items + if (is.null(items)) { + stop("Could not extract IRT item parameters.", call. = FALSE) + } + as.data.frame(items, stringsAsFactors = FALSE) +} + +#' Locate the `par` row in a mixedmirt coefficient block +#' +#' Supported `mirt` `coef,MixedClass` methods always emit a `par` row. A missing +#' `par` row is fail-closed rather than silently using the first row. +#' +#' @param mat A coefficient matrix from `mirt::coef()` on a `MixedClass` fit. +#' @return The string `"par"`. +#' @noRd +.mixedmirtParRow <- function(mat) { + rn <- rownames(mat) + if (is.null(rn) || !("par" %in% rn)) { + stop("mixedmirt coefficient block has no 'par' row.", call. = FALSE) + } + "par" +} + +#' Extract IRT item parameters from a mixedmirt MixedClass fit +#' +#' Supported `mirt` versions do not implement `IRTpars` / `simplify` on +#' `coef,MixedClass-method`. Extraction therefore reads the named item list +#' from `mirt::coef(fit)` and accepts either IRT `a`/`b` or slope-intercept +#' `a1`/`d` (`b = -d / a1` for Rasch). Other column sets fail closed. +#' +#' @param fit A `MixedClass` object returned by `.mixedmirt()`. +#' @return A data.frame with at least columns `a` and `b`. +#' @noRd +.extractMixedmirtIrtItems <- function(fit) { + if (!methods::is(fit, "MixedClass")) { + stop("mixedmirt recovery extraction requires a MixedClass fit.", call. = FALSE) + } + coefs <- mirt::coef(fit) + reserved <- c("GroupPars", "group", "items", "Theta") + item_names <- setdiff(names(coefs), reserved) + if (!length(item_names)) { + stop("Could not extract mixedmirt item parameters.", call. = FALSE) + } + rows <- lapply(item_names, function(item_name) { + block <- as.matrix(coefs[[item_name]]) + par_row <- .mixedmirtParRow(block) + cols <- colnames(block) + if (all(c("a", "b") %in% cols)) { + data.frame( + a = as.numeric(block[par_row, "a"]), + b = as.numeric(block[par_row, "b"]), + row.names = item_name, + stringsAsFactors = FALSE + ) + } else if (all(c("a1", "d") %in% cols)) { + a_hat <- as.numeric(block[par_row, "a1"]) + d_hat <- as.numeric(block[par_row, "d"]) + data.frame( + a = a_hat, + b = -d_hat / a_hat, + row.names = item_name, + stringsAsFactors = FALSE + ) + } else { + stop( + "mixedmirt item '", item_name, "' lacks IRT a/b or a1/d columns.", + call. = FALSE + ) + } + }) + do.call(rbind, rows) +} + +#' Extract mixedmirt item Wald intervals after name alignment +#' +#' @param fit A `MixedClass` object returned by `.mixedmirt()`. +#' @param items Character vector of item names to extract, in that order. +#' @return A data.frame with `item`, `estimate`, `lower`, and `upper` for `b`. +#' @noRd +.extractMixedmirtItemIntervals <- function(fit, items) { + if (!methods::is(fit, "MixedClass")) { + stop("mixedmirt interval extraction requires a MixedClass fit.", call. = FALSE) + } + if (!is.character(items) || !length(items)) { + stop("items must be a non-empty character vector.", call. = FALSE) + } + coefs <- mirt::coef(fit) + rows <- lapply(items, function(item_name) { + block <- coefs[[item_name]] + if (is.null(block)) { + stop("mixedmirt coefficients do not contain item '", item_name, "'.", call. = FALSE) + } + mat <- as.matrix(block) + if (!all(c("CI_2.5", "CI_97.5") %in% rownames(mat))) { + stop( + "mixedmirt item '", item_name, "' has no Wald interval rows.", + call. = FALSE + ) + } + par_row <- .mixedmirtParRow(mat) + if ("b" %in% colnames(mat)) { + estimate <- as.numeric(mat[par_row, "b"]) + lower <- as.numeric(mat["CI_2.5", "b"]) + upper <- as.numeric(mat["CI_97.5", "b"]) + } else if (all(c("a1", "d") %in% colnames(mat))) { + a_hat <- as.numeric(mat[par_row, "a1"]) + estimate <- -as.numeric(mat[par_row, "d"]) / a_hat + # Intervals are on the intercept d; map through b = -d / a1. + lower <- -as.numeric(mat["CI_97.5", "d"]) / a_hat + upper <- -as.numeric(mat["CI_2.5", "d"]) / a_hat + } else { + stop( + "mixedmirt item '", item_name, "' lacks interval columns for b or d.", + call. = FALSE + ) + } + data.frame( + item = item_name, + estimate = estimate, + lower = min(lower, upper), + upper = max(lower, upper), + stringsAsFactors = FALSE + ) + }) + do.call(rbind, rows) +} + +#' Extract the nested group-variance estimand from a mixedmirt fit +#' +#' @param fit A `MixedClass` object returned by `.mixedmirt()`. +#' @param group_name Name of the grouping factor used in `random = ~ 1|G`. +#' @return A list with `estimate`, `se`, `lower`, `upper`, and `name`. +#' @noRd +.extractMixedmirtGroupVariance <- function(fit, group_name = "group") { + if (!methods::is(fit, "MixedClass")) { + stop("Group-variance extraction requires a MixedClass fit.", call. = FALSE) + } + if (!is.character(group_name) || length(group_name) != 1L || !nzchar(group_name)) { + stop("group_name must be a single non-empty string.", call. = FALSE) + } + coefs <- mirt::coef(fit) + block <- coefs[[group_name]] + if (is.null(block)) { + stop( + "mixedmirt coefficients do not contain the '", group_name, "' random effect.", + call. = FALSE + ) + } + mat <- as.matrix(block) + cov_cols <- grep("^COV_", colnames(mat), value = TRUE) + if (!length(cov_cols)) { + stop( + "mixedmirt '", group_name, "' block has no COV_ variance column.", + call. = FALSE + ) + } + if (length(cov_cols) > 1L) { + stop( + "Crossed or multiple random-effect variances are outside this recovery slice.", + call. = FALSE + ) + } + par_row <- .mixedmirtParRow(mat) + estimate <- as.numeric(mat[par_row, cov_cols]) + lower <- if ("CI_2.5" %in% rownames(mat)) as.numeric(mat["CI_2.5", cov_cols]) else NA_real_ + upper <- if ("CI_97.5" %in% rownames(mat)) as.numeric(mat["CI_97.5", cov_cols]) else NA_real_ + se <- if (is.finite(lower) && is.finite(upper)) { + (upper - lower) / (2 * stats::qnorm(0.975)) + } else { + NA_real_ + } + list( + estimate = estimate, + se = se, + lower = lower, + upper = upper, + name = cov_cols + ) +} + +#' Record whether known-true values fall inside supplied intervals +#' +#' This is interval-inclusion evidence, not a Monte Carlo coverage rate. +#' Five seeds cannot support a nominal 95% coverage claim. +#' +#' @param estimate Numeric recovered values. +#' @param truth Numeric known-true values. +#' @param lower Numeric interval lower bounds. +#' @param upper Numeric interval upper bounds. +#' @return A data.frame with estimates, bounds, and a `covers` indicator. +#' @noRd +.recoveryIntervalInclusion <- function(estimate, truth, lower, upper) { + estimate <- as.numeric(estimate) + truth <- as.numeric(truth) + lower <- as.numeric(lower) + upper <- as.numeric(upper) + n <- length(estimate) + if (length(truth) != n || length(lower) != n || length(upper) != n) { + stop("Interval inclusion vectors must have the same length.", call. = FALSE) + } + if (!n) { + stop("Cannot evaluate interval inclusion on empty vectors.", call. = FALSE) + } + if (any(!is.finite(c(estimate, truth, lower, upper)))) { + stop("Interval inclusion requires finite estimates, truth, and bounds.", call. = FALSE) + } + if (any(upper < lower)) { + stop("Interval upper bound must be at least the lower bound.", call. = FALSE) + } + data.frame( + estimate = estimate, + truth = truth, + lower = lower, + upper = upper, + covers = truth >= lower & truth <= upper, + stringsAsFactors = FALSE + ) +} + +#' Refuse a single-level mixedmirt recovery design +#' +#' Collapsing persons across groups and fitting a single-level IRT model would +#' treat clustered observations as independent (the atomistic fallacy). +#' +#' @param covdata Person-level covariate data.frame containing the grouping factor. +#' @param random A random-effect formula or list of formulas. +#' @return Invisibly `TRUE` when the design is nested and non-atomistic. +#' @noRd +.rejectAtomisticMixedmirtRecovery <- function(covdata, random) { + if (is.null(covdata) || !is.data.frame(covdata) || !nrow(covdata)) { + stop( + "mixedmirt recovery refuses a single-level fit; that would be the atomistic fallacy.", + call. = FALSE + ) + } + if (is.null(random)) { + stop( + "mixedmirt recovery requires a nested random-effect formula.", + call. = FALSE + ) + } + if (inherits(random, "list") && !inherits(random, "formula")) { + if (length(random) != 1L) { + stop( + "Crossed random effects are outside this recovery slice.", + call. = FALSE + ) + } + random <- random[[1]] + } + if (!inherits(random, "formula")) { + stop( + "mixedmirt recovery requires a nested random-effect formula.", + call. = FALSE + ) + } + random_text <- paste(deparse(random), collapse = " ") + if (!grepl("|", random_text, fixed = TRUE)) { + stop( + "mixedmirt recovery requires a nested grouping term of the form ~ 1 | G.", + call. = FALSE + ) + } + if (grepl(":", random_text, fixed = TRUE) || + length(gregexpr("|", random_text, fixed = TRUE)[[1]]) > 1L) { + stop( + "Crossed or interaction random effects are outside this recovery slice.", + call. = FALSE + ) + } + invisible(TRUE) +} + +#' Simulate a nested two-level Rasch design with known true parameters +#' +#' Persons are nested in groups. Ability is +#' `\theta_{ig} = u_g + e_{ig}` with `u_g ~ N(0, \tau_{00})` and +#' `e_{ig} ~ N(0, \sigma^2)`. This is the Chalmers (2015) / `mixedmirt` +#' random-groups data-generating process, not a collapsed single-level sample. +#' +#' @param seed Integer random seed for the response draw. +#' @param n_groups Number of groups. Must be at least 30 (Maas & Hox, 2005). +#' @param n_per_group Number of persons per group. Must be at least 5. +#' @param tau00 True group-level intercept variance. +#' @param residual_variance True person-level residual variance. +#' @param true_b Known Rasch difficulties. +#' @param group_name Name of the grouping factor in `covdata`. +#' @return A list with responses, covariates, true item table, and `tau00`. +#' @noRd +.mixedmirtNestedRecoveryDesign <- function(seed, + n_groups = 40L, + n_per_group = 20L, + tau00 = 0.5, + residual_variance = 0.5, + true_b = c(-1.2, -0.6, 0, 0.6, 1.2), + group_name = "group") { + if (length(seed) != 1L || !is.finite(seed)) { + stop("seed must be a single finite value.", call. = FALSE) + } + if (n_groups < 30L) { + stop( + "Nested mixedmirt recovery requires at least 30 groups (Maas & Hox, 2005).", + call. = FALSE + ) + } + if (n_per_group < 5L) { + stop("Nested mixedmirt recovery requires at least 5 persons per group.", call. = FALSE) + } + if (!is.finite(tau00) || tau00 <= 0) { + stop("Group variance tau00 must be a positive finite value.", call. = FALSE) + } + if (!is.finite(residual_variance) || residual_variance <= 0) { + stop("Person residual variance must be a positive finite value.", call. = FALSE) + } + if (!is.numeric(true_b) || !length(true_b) || any(!is.finite(true_b))) { + stop("true_b must be a non-empty finite numeric vector.", call. = FALSE) + } + if (!is.character(group_name) || length(group_name) != 1L || !nzchar(group_name)) { + stop("group_name must be a single non-empty string.", call. = FALSE) + } + + set.seed(seed) + n <- as.integer(n_groups) * as.integer(n_per_group) + group <- factor(rep(paste0("G", seq_len(n_groups)), each = n_per_group)) + u <- stats::rnorm(n_groups, mean = 0, sd = sqrt(tau00)) + theta <- u[as.integer(group)] + stats::rnorm(n, mean = 0, sd = sqrt(residual_variance)) + true_a <- matrix(rep(1, length(true_b)), ncol = 1) + true_d <- matrix(-true_a[, 1] * true_b, ncol = 1) + # mirt >= 1.46 rejects itemtype = "Rasch" in simdata as ambiguous. + # a = 1 with 2PL is the dichotomous Rasch generating model. + response_data <- as.data.frame(mirt::simdata( + a = true_a, + d = true_d, + itemtype = "2PL", + Theta = matrix(theta, ncol = 1) + )) + names(response_data) <- paste0("Item", seq_len(ncol(response_data))) + covdata <- data.frame(group = group, stringsAsFactors = TRUE) + names(covdata)[1] <- group_name + random_formula <- stats::as.formula(paste("~ 1 |", group_name)) + .rejectAtomisticMixedmirtRecovery(covdata, random_formula) + list( + response_data = response_data, + covdata = covdata, + truth_items = data.frame( + a = true_a[, 1], + b = as.numeric(true_b), + row.names = names(response_data), + stringsAsFactors = FALSE + ), + truth_tau00 = as.numeric(tau00), + residual_variance = as.numeric(residual_variance), + n_groups = as.integer(n_groups), + n_per_group = as.integer(n_per_group), + group_name = group_name, + random_formula = random_formula, + estimands = c("b", "tau00") + ) +} + +#' Fit the registered nested mixedmirt recovery model +#' +#' Uses `fixed = ~ 0 + items` and `random = ~ 1 | group` (Chalmers, 2015). +#' Does not fit `lr.random`, crossed lists, multiple-membership weights, or +#' time-indexed membership. +#' +#' @param design Output of `.mixedmirtNestedRecoveryDesign()`. +#' @param NCYCLES MH-RM cycle budget passed to `.mixedmirt()`. +#' @param BURNIN MH-RM burn-in. Must exceed the default `RANDSTART` of 100. +#' @param SEMCYCLES MH-RM SEM cycle budget. +#' @return A `MixedClass` fit, or `NULL` when `.mixedmirt()` rejects the fit. +#' @noRd +.fitMixedmirtNestedRecovery <- function(design, + NCYCLES = 400, + BURNIN = 200, + SEMCYCLES = 100) { + if (!is.list(design) || is.null(design$response_data) || is.null(design$covdata)) { + stop("design must come from .mixedmirtNestedRecoveryDesign().", call. = FALSE) + } + .rejectAtomisticMixedmirtRecovery(design$covdata, design$random_formula) + .mixedmirt( + data = design$response_data, + model = 1, + itemtype = "Rasch", + SE = TRUE, + GenRandomPars = FALSE, + calcNull = FALSE, + covdata = design$covdata, + fixed = ~ 0 + items, + random = design$random_formula, + lr.fixed = ~ 1, + lr.random = NULL, + leniency = FALSE, + NCYCLES = NCYCLES, + BURNIN = BURNIN, + SEMCYCLES = SEMCYCLES + ) +} + +#' Honest recovery-coverage log for buyer-facing estimands +#' +#' @return A data.frame with `surface`, `status`, and `reason`. +#' @noRd +.recoveryCoverageExclusions <- function() { + data.frame( + surface = c( + "unidimensional 2PL via .mirt", + "AEFA greedy search on unidimensional 2PL", + "mixedmirt nested two-level random intercept", + "mixedmirt 2PL lr.random multilevel", + "multiple-membership weights", + "crossed random effects", + "time-flow / longitudinal membership" + ), + status = c( + "covered", + "covered when RUN_FULL_AEFA_TESTS=1", + "covered", + "excluded", + "excluded", + "excluded", + "excluded" + ), + reason = c( + "Known-true 2PL simulation with IRT a/b RMSE.", + "Known-true 2PL simulation through aefa() with five seeds.", + paste( + "Known-true nested Rasch via .mixedmirt(fixed=~0+items, random=~1|group);", + "recovers item b and group variance tau00 without collapsing clusters." + ), + paste( + "Chalmers (2015) supports lr.random for non-Rasch multilevel IRT;", + "no registered 2PL mixedmirt recovery design yet." + ), + paste( + "kaefa-core has no multiple-membership weight matrix;", + "Chung and Beretvas (2012) show that ignoring multiple membership", + "biases variance components." + ), + paste( + "mixedmirt can estimate list(~1|G1, ~1|G2), but no known-true", + "crossed-variance design is registered." + ), + "No time-indexed or time-varying membership design is implemented in kaefa-core." + ), + stringsAsFactors = FALSE + ) +} diff --git a/R/utils.R b/R/utils.R index 0c9ffdf..95fe14c 100644 --- a/R/utils.R +++ b/R/utils.R @@ -197,11 +197,13 @@ } if(exists('mod')){ - if(mod@OptimInfo$converged){ + if(isTRUE(mod@OptimInfo$converged)){ if(leniency){ mod } else { - if(mod@OptimInfo$secondordertest){ + # mirt can leave secondordertest as NA when the Hessian is unavailable; + # `if (NA)` would abort the caller instead of treating the fit as failed. + if(isTRUE(mod@OptimInfo$secondordertest)){ mod } else { NULL @@ -229,11 +231,11 @@ SEMCYCLES = SEMCYCLES, symmetric = symmetric)) if(exists('mod')){ - if(mod@OptimInfo$converged){ + if(isTRUE(mod@OptimInfo$converged)){ if(leniency){ mod } else { - if(mod@OptimInfo$secondordertest){ + if(isTRUE(mod@OptimInfo$secondordertest)){ mod } else { NULL diff --git a/docs/papers/README.md b/docs/papers/README.md index cc2caa8..1e7af46 100644 --- a/docs/papers/README.md +++ b/docs/papers/README.md @@ -119,10 +119,71 @@ cited with its DOI. Open-access / preprint links are noted where available. needed to reconstruct it. Kaefa accepts DIC only when the fitted model supplies a finite DIC value and never relabels AIC as DIC. +## 7. Monte Carlo parameter recovery RMSE + +- **Source:** Harwell, M. R., Stone, C. A., Hsu, T.-C., & Kirisci, L. (1996). + Monte Carlo studies in item response theory. *Applied Psychological + Measurement, 20*(2), 101-125. + DOI: [10.1177/014662169602000201](https://doi.org/10.1177/014662169602000201) +- **Canonical equation.** For recovered parameters \(\hat{\theta}\) and known + true parameters \(\theta\), + + RMSE = sqrt( mean( (hat_theta - theta)^2 ) ) + + Harwell et al. treat RMSE (and related Monte Carlo error summaries) as the + standard way to judge whether an IRT estimator recovers a known generating + model. kaefa uses that definition on IRT `a` and `b` after name alignment, + and on the nested mixedmirt group-variance estimand `tau00`. +- **Usage in kaefa:** internal helpers in `R/recovery.R` and the five-repeat + protocols in `docs/traceability/aefa-parameter-recovery.md` and + `docs/traceability/mixedmirt-parameter-recovery.md`. The formula, + alignment, and five-repeat schema are pinned by + `tests/testthat/test-aefa-parameter-recovery.R` and + `tests/testthat/test-mixedmirt-parameter-recovery.R`. Live `.mirt`, + `aefa()`, and `.mixedmirt` recovery fits live in + `tests/testthat/test-aefa-recovery-fits.R` and + `tests/testthat/test-mixedmirt-recovery-fits.R`. +- **Boundary.** `lr.random` 2PL multilevel, multiple-membership weights, + crossed random-effect recovery, and time-flow designs remain exclusions + until a true-parameter design is registered. The engine remains R/`mirt`; + this protocol does not introduce a Rust or GPU numeric core. + +## 8. Nested mixedmirt random-effect recovery + +- **Sources.** + - Chalmers, R. P. (2015). Extended mixed-effects item response models with + the MH-RM algorithm. *Journal of Educational Measurement, 52*(2), + 200-222. + DOI: [10.1111/jedm.12072](https://doi.org/10.1111/jedm.12072) + - Kamata, A. (2001). Item analysis by the hierarchical generalized linear + model. *Journal of Educational Measurement, 38*(1), 79-93. + DOI: [10.1111/j.1745-3984.2001.tb01117.x](https://doi.org/10.1111/j.1745-3984.2001.tb01117.x) + - Maas, C. J. M., & Hox, J. J. (2005). Sufficient sample sizes for + multilevel modeling. *Methodology, 1*(3), 86-92. + DOI: [10.1027/1614-2241.1.3.86](https://doi.org/10.1027/1614-2241.1.3.86) + - Diez-Roux, A. V. (1998). Bringing context back into epidemiology: + Variables and fallacies in multilevel analysis. *American Journal of + Public Health, 88*(2), 216-222. + DOI: [10.2105/AJPH.88.2.216](https://doi.org/10.2105/AJPH.88.2.216) +- **Canonical design.** Persons nested in groups; Rasch intercepts as + `fixed = ~ 0 + items`; group intercept as `random = ~ 1 | group`. Ability + is `theta_ig = u_g + e_ig` with known `tau00 = Var(u_g)`. A single-level + fit of the same responses is refused (atomistic fallacy). +- **Usage in kaefa:** `.mixedmirtNestedRecoveryDesign()`, + `.fitMixedmirtNestedRecovery()`, and the RMSE / interval helpers in + `R/recovery.R`. Live evidence is + `tests/testthat/test-mixedmirt-recovery-fits.R`. +- **Boundary.** Multiple membership (Chung & Beretvas, 2012), crossed + random effects, `lr.random` 2PL, and longitudinal membership are logged + as exclusions. Copyrighted PDFs are not redistributed; each source is + cited with its DOI. + ## Audit note kaefa does **not** re-implement `P(theta)`, the MML-EM E-/M-step, `S-X2`, `infit`, or `outfit`; those are delegated verbatim to `mirt` and remain subject to `mirt`'s validation. Package-local formulas and decision rules are pinned above: -the `Zh` cutoff, the exact Hurvich-Tsai AICc correction, and the explicit -posterior-information boundary that prevents DIC from being fabricated. +the `Zh` cutoff, the exact Hurvich-Tsai AICc correction, the explicit +posterior-information boundary that prevents DIC from being fabricated, and +the Harwell et al. RMSE recovery definition, and the nested mixedmirt +`b` / `tau00` recovery design. diff --git a/docs/product/kaefa-core-api-contract.md b/docs/product/kaefa-core-api-contract.md index 5d81b81..a9c0eeb 100644 --- a/docs/product/kaefa-core-api-contract.md +++ b/docs/product/kaefa-core-api-contract.md @@ -172,6 +172,10 @@ exist and pass in CI: - theta-prior helpers cover fit success, unsupported distribution, and fallback behavior. - dot-prefixed helpers are either documented or covered by deprecation tests. +- true-parameter recovery reports five-repeat RMSE for unidimensional 2PL + `a`/`b` after item-name alignment, and for nested mixedmirt Rasch `b` plus + group variance `tau00`. Multiple-membership, crossed, `lr.random` 2PL, and + time-flow surfaces remain exclusions until a registered design exists. ## Split Readiness Rule diff --git a/docs/superpowers/plans/2026-07-02-kaefa-2b-krw-sale-readiness.md b/docs/superpowers/plans/2026-07-02-kaefa-2b-krw-sale-readiness.md index d5192fe..4be08d6 100644 --- a/docs/superpowers/plans/2026-07-02-kaefa-2b-krw-sale-readiness.md +++ b/docs/superpowers/plans/2026-07-02-kaefa-2b-krw-sale-readiness.md @@ -251,11 +251,11 @@ Guardrails: steps: - uses: actions/checkout@de0fac2e4500dabe0009e67214ff5f5447ce83dd # v6.0.2 - - uses: r-lib/actions/setup-r@6f6e5bc62fba3a704f74e7ad7ef7676c5c6a2590 # v2 + - uses: r-lib/actions/setup-r@d3c5be51b12e724e68f33216ca3c148b66d5f0b6 # v2 with: use-public-rspm: true - - uses: r-lib/actions/setup-r-dependencies@6f6e5bc62fba3a704f74e7ad7ef7676c5c6a2590 # v2 + - uses: r-lib/actions/setup-r-dependencies@d3c5be51b12e724e68f33216ca3c148b66d5f0b6 # v2 with: extra-packages: any::testthat needs: check diff --git a/docs/traceability/aefa-parameter-recovery.md b/docs/traceability/aefa-parameter-recovery.md new file mode 100644 index 0000000..489310e --- /dev/null +++ b/docs/traceability/aefa-parameter-recovery.md @@ -0,0 +1,70 @@ +# AEFA true-parameter recovery protocol + +## Decision + +Buyer-facing recovery evidence for `kaefa-core` is a five-repeat RMSE protocol +against known-true item parameters. The protocol is implemented by internal +helpers in `R/recovery.R` and pinned by +`tests/testthat/test-aefa-parameter-recovery.R` (formula, alignment, schema, +coverage exclusions) plus `tests/testthat/test-aefa-recovery-fits.R` (live +fits). + +This is not a public API. External behaviour of `aefa()` and `engineAEFA()` is +unchanged. `.mirt()` / `.mixedmirt()` now treat a missing +(`NA`) second-order test as non-convergence when `leniency` is false, instead +of aborting on `if (NA)`. + +## RMSE definition + +For a recovered parameter vector \(\hat{\theta}\) and a true vector \(\theta\), + + RMSE = sqrt( mean( (hat_theta - theta)^2 ) ) + +Items are aligned by name before the difference is taken. Estimated and true +tables must contain the same item names; a partial intersection is rejected. +The required IRT columns for the unidimensional 2PL case are `a` +(discrimination) and `b` (difficulty) from +`mirt::coef(..., IRTpars = TRUE, simplify = TRUE)$items`. + +Exactly five seeds are required, and each parameter must have one RMSE value +for every seed. The summary schema is: per-run `seed`, `parameter`, `rmse`, +plus `mean_rmse` and `sd_rmse` by parameter. + +## Current coverage + +Covered: + +- Unidimensional 2PL recovery through `kaefa::.mirt` (`N = 1500`, + `SE = TRUE`, same cycle budget as the FIIFM stability fit). +- AEFA greedy search on the same 2PL design when `RUN_FULL_AEFA_TESTS=1`. +- Nested two-level Rasch random-intercept recovery through + `kaefa::.mixedmirt` (`fixed = ~ 0 + items`, `random = ~ 1 | group`). + See `docs/traceability/mixedmirt-parameter-recovery.md`. + +Explicitly excluded until a registered design exists: + +- `mixedmirt` 2PL / `lr.random` multilevel recovery. +- Multiple-membership weights (person in more than one group). +- Crossed random effects claimed as recovered crossed variances. +- Time-flow / longitudinal membership. + +The exclusion log is asserted in the recovery contract so a later claim cannot +silently treat those surfaces as covered. + +## Compatibility and rollback + +The helpers are unused by the estimation loop. Removing `R/recovery.R` and the +two test files restores the previous evidence surface. Do not export the +helpers or substitute a different error metric without updating this note and +the formula tests together. + +## References + +Harwell, M. R., Stone, C. A., Hsu, T.-C., & Kirisci, L. (1996). Monte Carlo +studies in item response theory. *Applied Psychological Measurement, 20*(2), +101-125. + + +Chalmers, R. P. (2012). mirt: A multidimensional item response theory package +for the R environment. *Journal of Statistical Software, 48*(6), 1-29. + diff --git a/docs/traceability/mixedmirt-parameter-recovery.md b/docs/traceability/mixedmirt-parameter-recovery.md new file mode 100644 index 0000000..8506006 --- /dev/null +++ b/docs/traceability/mixedmirt-parameter-recovery.md @@ -0,0 +1,169 @@ +# mixedmirt nested random-effect recovery protocol + +## Decision + +The next buyer-visible slice of issue #46, after the unidimensional 2PL +protocol in `docs/traceability/aefa-parameter-recovery.md` (PR #79), is +honest `mixedmirt` recovery of a **nested two-level Rasch random intercept**. +Persons are nested in groups. The fit uses kaefa's `.mixedmirt()` wrapper +with Chalmers' (2015) registered specification + + fixed = ~ 0 + items + random = ~ 1 | group + +and does **not** collapse clusters into a single-level sample. Collapsing +would treat clustered responses as independent and is refused as the +atomistic fallacy (Diez-Roux, 1998). + +Helpers live in `R/recovery.R` and reuse PR #79's RMSE definition, identical +item-name alignment, and complete five-seed summary. They are not a public +API. External behaviour of `aefa()` and `engineAEFA()` is unchanged. + +## Estimands + +Known-true generating process (Chalmers, 2015, random-groups example): + + theta_ig = u_g + e_ig + u_g ~ N(0, tau00) + e_ig ~ N(0, sigma^2) + +with `tau00 = 0.5`, `sigma^2 = 0.5`, 40 groups, 20 persons per group, and +five Rasch items whose difficulties are + + b = (-1.2, -0.6, 0, 0.6, 1.2) + +Recovered estimands: + +- item difficulty `b` after identical-name alignment +- group-level intercept variance `tau00` from the mixedmirt `COV_group` + block + +Rasch slopes are constrained to 1 and are not treated as free estimands. +Person residual variance is part of the generating process but is not a +buyer-facing recovery target in this slice. + +## RMSE and interval evidence + +RMSE is the Harwell, Stone, Hsu, and Kirisci (1996) definition already +pinned by PR #79: + + RMSE = sqrt( mean( (hat_theta - theta)^2 ) ) + +Accuracy bounds, applied to both the one-seed live fit and the five-seed +mean: + +- `RMSE(b) < 0.35` +- `RMSE(tau00) < 0.35` + +The item bound matches the unidimensional 2PL protocol. The variance bound +is the same numerical gate, not a loosened one. Maas and Hox (2005) showed +that with at least 30 groups the sampling error of a level-2 variance is +large relative to a fixed effect but still finite; 40 groups of 20 is above +that floor, so a 0.35 RMSE remains a real recovery claim rather than a +vacuous ceiling. + +Where mixedmirt returns Wald rows (`CI_2.5`, `CI_97.5`), the protocol +records whether the known-true value falls inside the interval. That is +**interval-inclusion evidence**, not a Monte Carlo coverage rate. Five +repeats cannot support a nominal 95% coverage claim (Harwell et al., 1996). +The live tests therefore require: + +- item-`b` inclusion on more than 60% of estimable item intervals +- a finite `tau00` interval on every completed seed + +They do not require the group-variance interval to cover on every seed. + +## Explicit exclusions + +Implemented now: nested two-level random intercept only. + +Honestly excluded until a registered known-true design exists: + +- `lr.random` 2PL / non-Rasch multilevel IRT (supported by mixedmirt; + Chalmers, 2015, notes that non-Rasch multilevel models use `lr.random`) +- multiple-membership weights (a person in more than one group). kaefa-core + has no membership-weight matrix. Chung and Beretvas (2012) show that + ignoring multiple membership biases variance components. +- crossed random effects such as `list(~ 1|school, ~ 1|rater)` or + `~ 1|group + ~ 1|items` claimed as recovered crossed variances +- time-flow / longitudinal / time-varying membership (Singer & Willett, + 2003; te Marvelde et al., 2006) + +The exclusion log is asserted in +`tests/testthat/test-mixedmirt-parameter-recovery.R` so these surfaces +cannot be claimed as covered. The engine remains R/`mirt`. This protocol +does not introduce a Rust or GPU numeric core. + +## MixedClass extraction contract + +Supported `mirt` versions implement `coef,MixedClass-method` without +`IRTpars` or `simplify`. Calling +`mirt::coef(fit, IRTpars = TRUE, simplify = TRUE)$items` on a +`MixedClass` fit therefore returns `NULL`. That is a fail-closed +contract for `.extractAefaIrtItems()`, not a recovery path. + +`.extractMixedmirtIrtItems()` reads the named item list from +`mirt::coef(fit)` and accepts either IRT `a`/`b` or slope-intercept +`a1`/`d` (`b = -d / a1` for Rasch). Missing item names, missing `par` +rows, and item blocks that lack both column sets stop. Interval +extraction likewise requires Wald rows (`CI_2.5`, `CI_97.5`) and `b` +or `d`. Group-variance extraction requires exactly one `COV_` column; +zero or more than one is fail-closed (the latter is the crossed-effects +exclusion). These branches are pinned by deterministic `MixedClass` +fixtures in `tests/testthat/test-mixedmirt-parameter-recovery.R`. + +## Tests + +- Always-on contracts (also in `test-fast`): + `tests/testthat/test-mixedmirt-parameter-recovery.R` +- One-seed live `.mixedmirt` fit in the full suite: + `tests/testthat/test-mixedmirt-recovery-fits.R` +- Five-seed live recovery behind `RUN_FULL_AEFA_TESTS=1`, using the same + complete-seed helper as PR #79 + +## Compatibility and rollback + +The helpers are unused by the estimation loop. Removing the mixedmirt +functions from `R/recovery.R` and the two mixedmirt test files restores the +PR #79 evidence surface. Do not export the helpers, do not treat a +single-level `.mirt` fit as multilevel recovery, and do not substitute a +different error metric without updating this note and the formula tests +together. + +## References + +Chalmers, R. P. (2015). Extended mixed-effects item response models with +the MH-RM algorithm. *Journal of Educational Measurement, 52*(2), 200-222. + + +Chung, H., & Beretvas, S. N. (2012). The impact of ignoring multiple +membership data structures in multilevel models. *British Journal of +Mathematical and Statistical Psychology, 65*(2), 185-200. + + +Diez-Roux, A. V. (1998). Bringing context back into epidemiology: Variables +and fallacies in multilevel analysis. *American Journal of Public Health, +88*(2), 216-222. + + +Harwell, M. R., Stone, C. A., Hsu, T.-C., & Kirisci, L. (1996). Monte Carlo +studies in item response theory. *Applied Psychological Measurement, +20*(2), 101-125. + + +Kamata, A. (2001). Item analysis by the hierarchical generalized linear +model. *Journal of Educational Measurement, 38*(1), 79-93. + + +Maas, C. J. M., & Hox, J. J. (2005). Sufficient sample sizes for multilevel +modeling. *Methodology, 1*(3), 86-92. + + +Singer, J. D., & Willett, J. B. (2003). *Applied longitudinal data +analysis: Modeling change and event occurrence*. Oxford University Press. + +te Marvelde, J. M., Glas, C. A. W., Van Landeghem, G., & Van Damme, J. +(2006). Application of multidimensional item response theory models to +longitudinal data. *Educational and Psychological Measurement, 66*(1), +5-34. + diff --git a/docs/traceability/r-lib-actions-v2.12.1.md b/docs/traceability/r-lib-actions-v2.12.1.md new file mode 100644 index 0000000..2960000 --- /dev/null +++ b/docs/traceability/r-lib-actions-v2.12.1.md @@ -0,0 +1,29 @@ +# r-lib/actions v2.12.1 공급망 추적 + +## 결정 + +세 개의 R 검증 workflow에서 사용하는 setup-pandoc, setup-r, +setup-r-dependencies, check-r-package를 v2.12.1 release commit +d3c5be51b12e724e68f33216ca3c148b66d5f0b6으로 통일한다. 전체 commit SHA 외의 태그·브랜치·짧은 SHA는 +회귀 계약이 거부한다. + +## 호환성 범위 + +공식 NEWS에 따르면 v2.12는 Node.js 24 전환, public RSPM 기본값 조정, +아키텍처별 cache key와 Pandoc 3.8.3을 포함하고, v2.12.1은 setup-r URL parser +경고와 Quarto action을 갱신한다. 현재 workflow의 R matrix, 권한, testthat 실행, +--no-tests 분리와 scheduled full-suite 계약은 변경하지 않는다. + +## 되돌리기 + +runner 또는 package 호환성 회귀가 확인되면 네 action을 함께 마지막 검증 SHA로 +되돌리고 R-CMD-check, fast/full test suite와 중앙 보안 검사를 같은 헤드에서 다시 +수행한다. 일부 action만 되돌리거나 이동 태그로 우회하지 않는다. + +## 참고문헌 + +R-lib. (2026, June 23). *r-lib/actions v2.12.1* [Software release notes]. +[NEWS.md](https://github.com/r-lib/actions/blob/d3c5be51b12e724e68f33216ca3c148b66d5f0b6/NEWS.md) + +R-lib. (2026). *Update NEWS for v1.12.1* [Source code commit]. +[d3c5be51b12e724e68f33216ca3c148b66d5f0b6](https://github.com/r-lib/actions/commit/d3c5be51b12e724e68f33216ca3c148b66d5f0b6) diff --git a/inst/WORDLIST b/inst/WORDLIST index f03f48c..f525206 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -2,6 +2,7 @@ 00273171 005 011 +02023 032408 10 100 @@ -19,6 +20,7 @@ 1500 153 153553 +16142241 1745 175 18 @@ -35,6 +37,7 @@ 2016 2017 21 +2105 222 2364 27 @@ -60,6 +63,7 @@ 51 518 52 +5282490 537 549 56 @@ -77,29 +81,40 @@ 76 82 83 +88 9218 AEFA AIC AICc +Adams AppVeyor BIC BURNIN Behavioral Behavioural Bentler +Beretvas Bifactor +Carlo +Chung DIC DIF +Diez DiscreteClass EM FUN +GPU GenRandomPars +Harwell +Hox +Hsu IRT Jennrich Jiao Jin Kamata Kang +Kirisci LAN LCA Leiman @@ -111,9 +126,12 @@ ML MM MMMM MPI +Maas Mansolf +Marvelde Mels MixedClass +Monte MultipleGroupClass NCYCLES NROW @@ -126,18 +144,26 @@ Psychometrika Q1 QMC RM +RMSE Rasch +Reckase Reise RemoteClusters Robbin +Roux +Rust SEMCYCLES SQUAREM SSH S_X2 Schmid +Singer SingleGroupClass +Stone Testlet UIRT +Wald +Willett X2 YAML Zh @@ -147,6 +173,7 @@ a0034177 aefa aefaInit aefaResults +atomistic bentlerQ bentlerT bifactorQ @@ -156,11 +183,14 @@ coefs covdata covdataClassifieder covdataFixedEffectComb +crossed customised cutoff demographical engineAEFA estModelTrials +estimand +estimands evaluateItemFit extractThetaOnly fitEMatUIRT @@ -186,14 +216,18 @@ knitr kwangwoon launchAEFA loadPercentage +longitudinal math maxExtraction maximising +membership minExtraction mirt mirtModel +mixedmirt modelSelectionCriteria ncat +nested newEngine nrow oblimax @@ -227,6 +261,7 @@ stylesheet summarise tandemI tandemII +tau00 testItemFit1 testMod1 testModel1 @@ -236,6 +271,7 @@ tracelines tryLCA turnOffMixedEst ubuntu +unidimensional v048 v1 z0 diff --git a/tests/FAST_TESTS.md b/tests/FAST_TESTS.md index 5629ada..5989696 100644 --- a/tests/FAST_TESTS.md +++ b/tests/FAST_TESTS.md @@ -13,6 +13,10 @@ testthat::test_file("tests/testthat/test-shiny-product-surface.R", reporter = reporter) testthat::test_file("tests/testthat/test-core-api-contract.R", reporter = reporter) +testthat::test_file("tests/testthat/test-aefa-parameter-recovery.R", + reporter = reporter) +testthat::test_file("tests/testthat/test-mixedmirt-parameter-recovery.R", + reporter = reporter) ``` `R-CMD-check` remains the package installation and multi-OS compatibility gate. diff --git a/tests/testthat/test-aefa-parameter-recovery.R b/tests/testthat/test-aefa-parameter-recovery.R new file mode 100644 index 0000000..1cd0cf3 --- /dev/null +++ b/tests/testthat/test-aefa-parameter-recovery.R @@ -0,0 +1,163 @@ +test_that("alignment and repeat helpers reject malformed inputs", { + .ensure_kaefa_namespace() + testthat::expect_error( + kaefa:::.alignIrtItemParameters(1:3, data.frame(a = 1, b = 1, row.names = "I1")), + "matrix or data.frame" + ) + testthat::expect_error( + kaefa:::.alignIrtItemParameters(data.frame(a = 1, b = 1, row.names = "I1"), 1:3), + "matrix or data.frame" + ) + testthat::expect_error(kaefa:::.summariseRecoveryRepeats(1:5), "data.frame") + testthat::expect_error( + kaefa:::.summariseRecoveryRepeats(data.frame(seed = 1L)), + "must contain columns" + ) + testthat::expect_error(kaefa:::.extractAefaIrtItems(list()), "aefa history or a mirt model") + testthat::expect_error( + kaefa:::.extractAefaIrtItems(list(estModelTrials = list())), + "no estimated model trials" + ) + testthat::expect_error( + kaefa:::.extractAefaIrtItems(structure(list(estModelTrials = list()), class = "aefa")), + "no estimated model trials" + ) + + testthat::skip_if_not_installed("mirt") + items <- data.frame(a = 1.2, b = -0.4, row.names = "Item1", stringsAsFactors = FALSE) + testthat::with_mocked_bindings( + coef = function(...) list(items = items), + .package = "mirt", + { + extracted <- kaefa:::.extractAefaIrtItems(methods::new("SingleGroupClass")) + testthat::expect_equal(extracted["Item1", "a"], 1.2) + testthat::expect_equal(extracted["Item1", "b"], -0.4) + from_history <- kaefa:::.extractAefaIrtItems( + structure( + list(estModelTrials = list(methods::new("SingleGroupClass"))), + class = "aefa" + ) + ) + testthat::expect_equal(from_history["Item1", "b"], -0.4) + } + ) +}) + +test_that("RMSE matches the Monte Carlo recovery definition", { + .ensure_kaefa_namespace() + estimated <- c(0.8, 1.1, 1.4) + truth <- c(1.0, 1.0, 1.0) + expected <- sqrt(mean((estimated - truth)^2)) + testthat::expect_equal(kaefa:::.parameterRecoveryRMSE(estimated, truth), expected) + testthat::expect_equal(expected, sqrt((0.04 + 0.01 + 0.16) / 3)) + testthat::expect_error(kaefa:::.parameterRecoveryRMSE(c(1, 2), 1), "same length") + testthat::expect_error(kaefa:::.parameterRecoveryRMSE(numeric(), numeric()), "empty") + testthat::expect_error(kaefa:::.parameterRecoveryRMSE(c(1, NA), c(1, 2)), "finite") +}) + +test_that("IRT parameter alignment is by item name and required columns", { + .ensure_kaefa_namespace() + estimated <- data.frame( + a = c(1.1, 0.9), + b = c(-0.2, 0.4), + row.names = c("Item2", "Item1") + ) + truth <- data.frame( + a = c(1.0, 1.2), + b = c(0.5, -0.1), + row.names = c("Item1", "Item2") + ) + aligned <- kaefa:::.alignIrtItemParameters(estimated, truth) + testthat::expect_identical(aligned$items, c("Item2", "Item1")) + testthat::expect_equal(aligned$estimated["Item1", "a"], 0.9) + testthat::expect_equal(aligned$truth["Item1", "b"], 0.5) + + unnamed <- estimated + rownames(unnamed) <- NULL + testthat::expect_error(kaefa:::.alignIrtItemParameters(unnamed, truth), "row names") + truth_unnamed <- truth + rownames(truth_unnamed) <- NULL + testthat::expect_error( + kaefa:::.alignIrtItemParameters(unnamed, truth_unnamed), + "row names" + ) + testthat::expect_error( + kaefa:::.alignIrtItemParameters(estimated[, "a", drop = FALSE], truth), + "Missing recovery columns" + ) + other <- truth + rownames(other) <- c("Q1", "Q2") + testthat::expect_error( + kaefa:::.alignIrtItemParameters(estimated, other), + "same item names" + ) + subset_truth <- truth[1, , drop = FALSE] + testthat::expect_error( + kaefa:::.alignIrtItemParameters(estimated, subset_truth), + "same item names" + ) +}) + +test_that("five-repeat recovery summary has a fixed output schema", { + .ensure_kaefa_namespace() + rmse_by_repeat <- data.frame( + seed = rep(c(11L, 22L, 33L, 44L, 55L), each = 2L), + parameter = rep(c("a", "b"), times = 5L), + rmse = c(0.10, 0.20, 0.12, 0.18, 0.11, 0.21, 0.09, 0.19, 0.13, 0.17), + stringsAsFactors = FALSE + ) + summarised <- kaefa:::.summariseRecoveryRepeats(rmse_by_repeat) + testthat::expect_identical(summarised$n_repeats, 5L) + testthat::expect_identical(names(summarised$summary), c("parameter", "n_repeats", "mean_rmse", "sd_rmse")) + testthat::expect_equal(summarised$summary$n_repeats, c(5L, 5L)) + a_rmse <- rmse_by_repeat$rmse[rmse_by_repeat$parameter == "a"] + testthat::expect_equal( + summarised$summary$mean_rmse[summarised$summary$parameter == "a"], + mean(a_rmse) + ) + testthat::expect_equal( + summarised$summary$sd_rmse[summarised$summary$parameter == "a"], + stats::sd(a_rmse) + ) + testthat::expect_error( + kaefa:::.summariseRecoveryRepeats(rmse_by_repeat[1:4, ]), + "exactly 5 repeats" + ) + incomplete_b <- rmse_by_repeat[rmse_by_repeat$parameter == "a" | + rmse_by_repeat$seed != 55L, ] + testthat::expect_error( + kaefa:::.summariseRecoveryRepeats(incomplete_b), + "exactly one RMSE value" + ) + duplicated_a <- rbind( + rmse_by_repeat, + data.frame(seed = 11L, parameter = "a", rmse = 0.99, stringsAsFactors = FALSE) + ) + testthat::expect_error( + kaefa:::.summariseRecoveryRepeats(duplicated_a), + "exactly one RMSE value" + ) + missing_rmse <- rmse_by_repeat + missing_rmse$rmse[1] <- NA_real_ + testthat::expect_error( + kaefa:::.summariseRecoveryRepeats(missing_rmse), + "missing values" + ) +}) + +test_that("recovery coverage log keeps MM and time-flow explicit exclusions", { + .ensure_kaefa_namespace() + coverage <- kaefa:::.recoveryCoverageExclusions() + testthat::expect_identical( + names(coverage), + c("surface", "status", "reason") + ) + testthat::expect_true(any(coverage$surface == "unidimensional 2PL via .mirt" & + coverage$status == "covered")) + testthat::expect_true(any(grepl("nested two-level random intercept", coverage$surface) & + coverage$status == "covered")) + testthat::expect_true(any(grepl("multiple-membership", coverage$surface) & + coverage$status == "excluded")) + testthat::expect_true(any(grepl("time-flow", coverage$surface) & + coverage$status == "excluded")) +}) diff --git a/tests/testthat/test-aefa-recovery-fits.R b/tests/testthat/test-aefa-recovery-fits.R new file mode 100644 index 0000000..e454b67 --- /dev/null +++ b/tests/testthat/test-aefa-recovery-fits.R @@ -0,0 +1,121 @@ +test_that(".mirt recovers known 2PL parameters with bounded RMSE", { + .ensure_kaefa_namespace() + testthat::skip_if_not_installed("mirt") + + set.seed(20260817) + true_a <- matrix(c(0.9, 1.1, 1.3, 1.5, 1.7), ncol = 1) + true_b <- c(-1.2, -0.6, 0, 0.6, 1.2) + true_d <- matrix(-true_a[, 1] * true_b, ncol = 1) + # Same design as test-fiifm-stability.R: N=1500 and SE=TRUE so the Hessian + # exists and leniency=FALSE can keep the fit. SE=FALSE leaves + # secondordertest NA/FALSE and .mirt() returns NULL. + response_data <- as.data.frame(mirt::simdata( + a = true_a, + d = true_d, + itemtype = "2PL", + N = 1500 + )) + names(response_data) <- paste0("Item", seq_len(ncol(response_data))) + truth <- data.frame( + a = true_a[, 1], + b = true_b, + row.names = names(response_data) + ) + + utils::capture.output( + fit <- suppressWarnings(kaefa::.mirt( + data = response_data, + model = 1, + method = "EM", + itemtype = "2PL", + SE = TRUE, + GenRandomPars = FALSE, + calcNull = FALSE, + leniency = FALSE, + NCYCLES = 400, + BURNIN = 100, + SEMCYCLES = 100 + )) + ) + if (!methods::is(fit, "SingleGroupClass")) { + testthat::fail( + "kaefa::.mirt did not return a single-group fit (converged second-order test required)" + ) + return(invisible(NULL)) + } + + estimated <- kaefa:::.extractAefaIrtItems(fit) + aligned <- kaefa:::.alignIrtItemParameters(estimated, truth) + rmse_a <- kaefa:::.parameterRecoveryRMSE(aligned$estimated$a, aligned$truth$a) + rmse_b <- kaefa:::.parameterRecoveryRMSE(aligned$estimated$b, aligned$truth$b) + testthat::expect_lt(rmse_a, 0.35) + testthat::expect_lt(rmse_b, 0.35) + testthat::expect_gt(stats::cor(aligned$estimated$a, aligned$truth$a), 0.90) +}) + +test_that("aefa five-seed 2PL recovery reports mean RMSE", { + .skip_expensive_ci_calls("aefa") + .ensure_kaefa_namespace() + testthat::skip_if_not_installed("mirt") + + true_a <- matrix(c(0.9, 1.1, 1.3, 1.5, 1.7), ncol = 1) + true_b <- c(-1.2, -0.6, 0, 0.6, 1.2) + true_d <- matrix(-true_a[, 1] * true_b, ncol = 1) + seeds <- c(20260817L, 20260818L, 20260819L, 20260820L, 20260821L) + rows <- list() + + for (seed in seeds) { + set.seed(seed) + response_data <- as.data.frame(mirt::simdata( + a = true_a, + d = true_d, + itemtype = "2PL", + N = 250 + )) + names(response_data) <- paste0("Item", seq_len(ncol(response_data))) + truth <- data.frame( + a = true_a[, 1], + b = true_b, + row.names = names(response_data) + ) + result <- try( + aefa( + response_data, + minExtraction = 1, + maxExtraction = 1, + turnOffMixedEst = TRUE, + skipggum = TRUE, + tryLCA = FALSE, + saveModelHistory = TRUE, + printItemFit = FALSE + ), + silent = TRUE + ) + if (inherits(result, "try-error") || is.null(result)) { + testthat::fail(paste("aefa recovery repeat failed for seed", seed)) + return(invisible(NULL)) + } + estimated <- kaefa:::.extractAefaIrtItems(result) + aligned <- kaefa:::.alignIrtItemParameters(estimated, truth) + rows[[length(rows) + 1L]] <- data.frame( + seed = seed, + parameter = c("a", "b"), + rmse = c( + kaefa:::.parameterRecoveryRMSE(aligned$estimated$a, aligned$truth$a), + kaefa:::.parameterRecoveryRMSE(aligned$estimated$b, aligned$truth$b) + ), + stringsAsFactors = FALSE + ) + } + + summarised <- kaefa:::.summariseRecoveryRepeats(do.call(rbind, rows)) + testthat::expect_identical(summarised$n_repeats, 5L) + testthat::expect_lt( + summarised$summary$mean_rmse[summarised$summary$parameter == "a"], + 0.45 + ) + testthat::expect_lt( + summarised$summary$mean_rmse[summarised$summary$parameter == "b"], + 0.45 + ) +}) diff --git a/tests/testthat/test-ci-action-pins.R b/tests/testthat/test-ci-action-pins.R new file mode 100644 index 0000000..4203052 --- /dev/null +++ b/tests/testthat/test-ci-action-pins.R @@ -0,0 +1,142 @@ +testthat::test_that("CI uses exactly the reviewed r-lib action references", { + reviewed_sha <- "d3c5be51b12e724e68f33216ca3c148b66d5f0b6" + workflow_paths <- c( + testthat::test_path("..", "..", ".github", "workflows", "R-CMD-check.yaml"), + testthat::test_path("..", "..", ".github", "workflows", "test-fast.yaml"), + testthat::test_path("..", "..", ".github", "workflows", "test-suite.yaml") + ) + expected_actions <- list( + "R-CMD-check.yaml" = c( + "setup-pandoc", + "setup-r", + "setup-r-dependencies", + "check-r-package" + ), + "test-fast.yaml" = c("setup-r", "setup-r-dependencies"), + "test-suite.yaml" = c("setup-r", "setup-r-dependencies") + ) + action_pattern <- paste0( + "r-lib/actions/", + "(setup-pandoc|setup-r-dependencies|setup-r|check-r-package)", + "@[^[:space:]#]+" + ) + + for (workflow_path in workflow_paths) { + workflow_text <- paste(readLines(workflow_path, warn = FALSE), collapse = "\n") + action_refs <- regmatches( + workflow_text, + gregexpr(action_pattern, workflow_text, perl = TRUE) + )[[1]] + expected_refs <- paste0( + "r-lib/actions/", + expected_actions[[basename(workflow_path)]], + "@", + reviewed_sha + ) + + testthat::expect_identical( + action_refs, + expected_refs, + info = paste("Unexpected r-lib action set in", basename(workflow_path)) + ) + } +}) + +testthat::test_that("R CMD check refreshes the reviewed dependency cache ABI", { + workflow_path <- testthat::test_path( + "..", "..", ".github", "workflows", "R-CMD-check.yaml" + ) + workflow_lines <- readLines(workflow_path, warn = FALSE) + dependency_step <- grep( + "r-lib/actions/setup-r-dependencies@", + workflow_lines, + fixed = TRUE + ) + testthat::expect_length(dependency_step, 1L) + dependency_block <- workflow_lines[ + dependency_step:min(dependency_step + 8L, length(workflow_lines)) + ] + + active_cache_version_pattern <- paste0( + "^[[:space:]]*cache-version:[[:space:]]*", + "['\\\"]2['\\\"][[:space:]]*(#.*)?$" + ) + testthat::expect_false( + grepl(active_cache_version_pattern, "# cache-version: '2'", perl = TRUE), + info = "A commented cache-version example must not satisfy the contract" + ) + testthat::expect_true( + any(grepl(active_cache_version_pattern, dependency_block, perl = TRUE)), + info = "The reviewed macOS TBB ABI cache refresh must remain explicit" + ) +}) + +testthat::test_that("CI rebuilds oneTBB dependents around dependency setup", { + workflow_paths <- c( + testthat::test_path("..", "..", ".github", "workflows", "R-CMD-check.yaml"), + testthat::test_path("..", "..", ".github", "workflows", "test-fast.yaml"), + testthat::test_path("..", "..", ".github", "workflows", "test-suite.yaml") + ) + + for (workflow_path in workflow_paths) { + workflow_lines <- readLines(workflow_path, warn = FALSE) + pre_step <- grep( + "name: Source-build stringfish before dependency resolution", + workflow_lines, + fixed = TRUE + ) + dependency_step <- grep( + "r-lib/actions/setup-r-dependencies@", + workflow_lines, + fixed = TRUE + ) + post_step <- grep( + "name: Rebuild RcppParallel-linked packages from source", + workflow_lines, + fixed = TRUE + ) + + testthat::expect_length(pre_step, 1L) + testthat::expect_length(dependency_step, 1L) + testthat::expect_length(post_step, 1L) + testthat::expect_true(pre_step < dependency_step) + testthat::expect_true(dependency_step < post_step) + + pre_block <- workflow_lines[ + pre_step:min(pre_step + 5L, length(workflow_lines)) + ] + post_block <- workflow_lines[ + post_step:min(post_step + 7L, length(workflow_lines)) + ] + testthat::expect_true(any(grepl( + 'install.packages(c("RcppParallel", "stringfish")', + pre_block, + fixed = TRUE + ))) + testthat::expect_true(any(grepl( + 'repos = "https://cloud.r-project.org"', + pre_block, + fixed = TRUE + ))) + testthat::expect_true(any(grepl( + 'type = "source"', + pre_block, + fixed = TRUE + ))) + testthat::expect_true(any(grepl( + 'dependsOnPkgs("RcppParallel", dependencies = "LinkingTo")', + post_block, + fixed = TRUE + ))) + testthat::expect_true(any(grepl( + 'install.packages("RcppParallel", type = "source"', + post_block, + fixed = TRUE + ))) + testthat::expect_true(any(grepl( + 'install.packages(linked, type = "source"', + post_block, + fixed = TRUE + ))) + } +}) diff --git a/tests/testthat/test-mixedmirt-parameter-recovery.R b/tests/testthat/test-mixedmirt-parameter-recovery.R new file mode 100644 index 0000000..4c4515b --- /dev/null +++ b/tests/testthat/test-mixedmirt-parameter-recovery.R @@ -0,0 +1,472 @@ +test_that("nested mixedmirt design rejects atomistic and undersized clusters", { + .ensure_kaefa_namespace() + testthat::skip_if_not_installed("mirt") + + testthat::expect_error( + kaefa:::.rejectAtomisticMixedmirtRecovery(NULL, ~ 1 | group), + "atomistic fallacy" + ) + testthat::expect_error( + kaefa:::.rejectAtomisticMixedmirtRecovery(data.frame(group = 1:40), NULL), + "nested random-effect formula" + ) + testthat::expect_error( + kaefa:::.rejectAtomisticMixedmirtRecovery( + data.frame(group = 1:40), + list(~ 1 | school, ~ 1 | rater) + ), + "Crossed random effects" + ) + testthat::expect_error( + kaefa:::.mixedmirtNestedRecoveryDesign(seed = 1L, n_groups = 10L), + "at least 30 groups" + ) + testthat::expect_error( + kaefa:::.mixedmirtNestedRecoveryDesign(seed = 1L, n_per_group = 2L), + "at least 5 persons" + ) + testthat::expect_error( + kaefa:::.mixedmirtNestedRecoveryDesign(seed = 1L, tau00 = 0), + "positive finite" + ) +}) + +test_that("nested mixedmirt design keeps cluster membership and known estimands", { + .ensure_kaefa_namespace() + testthat::skip_if_not_installed("mirt") + + design <- kaefa:::.mixedmirtNestedRecoveryDesign( + seed = 20260817L, + n_groups = 30L, + n_per_group = 5L + ) + testthat::expect_identical(nrow(design$response_data), 150L) + testthat::expect_identical(nrow(design$covdata), 150L) + testthat::expect_identical(nlevels(design$covdata$group), 30L) + testthat::expect_identical(as.character(design$estimands), c("b", "tau00")) + testthat::expect_identical(rownames(design$truth_items), names(design$response_data)) + testthat::expect_true(all(c("a", "b") %in% colnames(design$truth_items))) + testthat::expect_equal(design$truth_tau00, 0.5) + testthat::expect_false(is.null(design$random_formula)) + kaefa:::.rejectAtomisticMixedmirtRecovery(design$covdata, design$random_formula) +}) + +test_that("interval inclusion is exact and rejects incomplete bounds", { + .ensure_kaefa_namespace() + included <- kaefa:::.recoveryIntervalInclusion( + estimate = c(0.1, 0.9), + truth = c(0.0, 1.0), + lower = c(-0.2, 0.8), + upper = c(0.3, 0.85) + ) + testthat::expect_identical(included$covers, c(TRUE, FALSE)) + testthat::expect_error( + kaefa:::.recoveryIntervalInclusion(1, 1, 2, 0), + "upper bound" + ) + testthat::expect_error( + kaefa:::.recoveryIntervalInclusion(1, 1, NA_real_, 2), + "finite" + ) + testthat::expect_error( + kaefa:::.recoveryIntervalInclusion(1, c(1, 2), 0, 2), + "same length" + ) +}) + +test_that("mixedmirt extractors reject non-MixedClass and incomplete inputs", { + .ensure_kaefa_namespace() + testthat::expect_error( + kaefa:::.extractMixedmirtIrtItems(list()), + "MixedClass" + ) + testthat::expect_error( + kaefa:::.extractMixedmirtItemIntervals(list(), "Item1"), + "MixedClass" + ) + testthat::expect_error( + kaefa:::.extractMixedmirtGroupVariance(list(), "group"), + "MixedClass" + ) + testthat::expect_error( + kaefa:::.fitMixedmirtNestedRecovery(list()), + "mixedmirtNestedRecoveryDesign" + ) + testthat::expect_error( + kaefa:::.recoveryIntervalInclusion(numeric(), numeric(), numeric(), numeric()), + "empty" + ) + testthat::expect_error( + kaefa:::.rejectAtomisticMixedmirtRecovery(data.frame(group = 1:40), ~ 1), + "1 \\| G" + ) + testthat::expect_error( + kaefa:::.rejectAtomisticMixedmirtRecovery( + data.frame(group = 1:40), + ~ 1 | group:items + ), + "Crossed or interaction" + ) + testthat::expect_error( + kaefa:::.rejectAtomisticMixedmirtRecovery(data.frame(group = 1:40), "group"), + "nested random-effect formula" + ) + testthat::expect_true(isTRUE( + kaefa:::.rejectAtomisticMixedmirtRecovery( + data.frame(group = factor(rep(1:40, each = 2))), + list(~ 1 | group) + ) + )) + testthat::expect_error( + kaefa:::.mixedmirtNestedRecoveryDesign(seed = c(1, 2)), + "single finite" + ) + testthat::expect_error( + kaefa:::.mixedmirtNestedRecoveryDesign(seed = 1L, true_b = NA_real_), + "true_b" + ) + testthat::expect_error( + kaefa:::.mixedmirtNestedRecoveryDesign(seed = 1L, group_name = ""), + "group_name" + ) + testthat::expect_error( + kaefa:::.mixedmirtNestedRecoveryDesign(seed = 1L, residual_variance = 0), + "residual" + ) +}) + +.fake_mixedclass <- function(item_pars, + item_names = paste0("Item", seq_along(item_pars)), + item_se = NULL, + random_par = 0.5, + random_parnames = "COV_group", + random_se = 0.08, + random_name = "group", + include_random = TRUE, + group_se = 0.01) { + fit <- methods::new("MixedClass") + n_items <- length(item_pars) + dat <- as.data.frame( + matrix(0L, nrow = 2L, ncol = max(n_items, 0L)), + stringsAsFactors = FALSE + ) + if (n_items) { + names(dat) <- item_names + } + fit@Data <- list(K = rep(2L, n_items), data = dat) + + item_objects <- lapply(seq_len(n_items), function(i) { + item <- methods::new("dich") + item@par <- unname(item_pars[[i]]) + item@parnames <- names(item_pars[[i]]) + if (!is.null(item_se)) { + item@SEpar <- item_se[[i]] + } + item + }) + group_pars <- methods::new("GroupPars") + group_pars@par <- 0 + group_pars@parnames <- "MEAN_1" + if (!is.null(item_se) || !is.null(group_se)) { + group_pars@SEpar <- group_se + } + + random <- list() + if (isTRUE(include_random)) { + random_pars <- methods::new("RandomPars") + random_pars@par <- random_par + random_pars@parnames <- random_parnames + if (!is.null(random_se)) { + random_pars@SEpar <- random_se + } + random_pars@gdesign <- matrix( + 1, + nrow = 2L, + ncol = 1L, + dimnames = list(NULL, random_name) + ) + random <- list(random_pars) + } + + fit@ParObjects <- list( + pars = c(item_objects, list(group_pars)), + random = random, + lr.random = list() + ) + fit@Model <- list(lrPars = numeric(0)) + fit +} + +test_that("MixedClass IRTpars simplify is fail-closed on supported mirt", { + .ensure_kaefa_namespace() + testthat::skip_if_not_installed("mirt") + + fit <- .fake_mixedclass( + item_pars = list(c(a1 = 1, d = 0.6), c(a1 = 1, d = -0.3)), + item_se = list(c(0.05, 0.1), c(0.05, 0.1)) + ) + testthat::expect_true(methods::is(fit, "MixedClass")) + testthat::expect_null(mirt::coef(fit, IRTpars = TRUE, simplify = TRUE)$items) + testthat::expect_error( + kaefa:::.extractAefaIrtItems(fit), + "Could not extract IRT item parameters" + ) + testthat::expect_error( + kaefa:::.extractAefaIrtItems(list(estModelTrials = list(fit))), + "Could not extract IRT item parameters" + ) + testthat::expect_error( + kaefa:::.extractAefaIrtItems(structure(list(estModelTrials = list(fit)), class = "aefa")), + "Could not extract IRT item parameters" + ) +}) + +test_that("mixedmirt item extractor covers a/b, a1/d, and fail-closed shapes", { + .ensure_kaefa_namespace() + testthat::skip_if_not_installed("mirt") + + slope_intercept <- .fake_mixedclass( + item_pars = list(c(a1 = 1, d = 0.6), c(a1 = 2, d = -0.4)), + item_se = list(c(0.05, 0.1), c(0.05, 0.1)) + ) + from_ad <- kaefa:::.extractMixedmirtIrtItems(slope_intercept) + testthat::expect_equal(from_ad["Item1", "a"], 1) + testthat::expect_equal(from_ad["Item1", "b"], -0.6) + testthat::expect_equal(from_ad["Item2", "b"], 0.2) + + irt_ab <- .fake_mixedclass( + item_pars = list(c(a = 1.1, b = -0.5), c(a = 0.9, b = 0.4)), + item_se = list(c(0.05, 0.1), c(0.05, 0.1)) + ) + from_ab <- kaefa:::.extractMixedmirtIrtItems(irt_ab) + testthat::expect_equal(from_ab["Item1", "a"], 1.1) + testthat::expect_equal(from_ab["Item1", "b"], -0.5) + testthat::expect_equal(from_ab["Item2", "b"], 0.4) + + guessing_only <- .fake_mixedclass( + item_pars = list(c(g = 0.2), c(g = 0.1)), + item_se = list(0.01, 0.01) + ) + testthat::expect_error( + kaefa:::.extractMixedmirtIrtItems(guessing_only), + "lacks IRT a/b or a1/d" + ) + + empty_items <- .fake_mixedclass( + item_pars = list(), + include_random = FALSE, + item_se = NULL, + group_se = 0.01 + ) + testthat::expect_error( + kaefa:::.extractMixedmirtIrtItems(empty_items), + "Could not extract mixedmirt item parameters" + ) +}) + +test_that("mixedmirt interval extractor covers a/b, a1/d, and fail-closed inputs", { + .ensure_kaefa_namespace() + testthat::skip_if_not_installed("mirt") + + slope_intercept <- .fake_mixedclass( + item_pars = list(c(a1 = 1, d = 0.6)), + item_se = list(c(0.05, 0.1)) + ) + from_ad <- kaefa:::.extractMixedmirtItemIntervals(slope_intercept, "Item1") + testthat::expect_equal(from_ad$estimate, -0.6) + testthat::expect_lt(from_ad$lower, from_ad$estimate) + testthat::expect_gt(from_ad$upper, from_ad$estimate) + + irt_ab <- .fake_mixedclass( + item_pars = list(c(a = 1, b = -0.5)), + item_se = list(c(0.05, 0.1)) + ) + from_ab <- kaefa:::.extractMixedmirtItemIntervals(irt_ab, "Item1") + testthat::expect_equal(from_ab$estimate, -0.5) + testthat::expect_lt(from_ab$lower, from_ab$estimate) + testthat::expect_gt(from_ab$upper, from_ab$estimate) + + no_se <- .fake_mixedclass( + item_pars = list(c(a1 = 1, d = 0.6)), + item_se = NULL, + random_se = NULL, + group_se = NULL + ) + testthat::expect_error( + kaefa:::.extractMixedmirtItemIntervals(no_se, "Item1"), + "no Wald interval rows" + ) + testthat::expect_error( + kaefa:::.extractMixedmirtItemIntervals(slope_intercept, "MissingItem"), + "do not contain item" + ) + testthat::expect_error( + kaefa:::.extractMixedmirtItemIntervals(slope_intercept, character()), + "non-empty character" + ) + testthat::expect_error( + kaefa:::.extractMixedmirtItemIntervals(slope_intercept, 1), + "non-empty character" + ) + + guessing_only <- .fake_mixedclass( + item_pars = list(c(g = 0.2)), + item_se = list(0.01) + ) + testthat::expect_error( + kaefa:::.extractMixedmirtItemIntervals(guessing_only, "Item1"), + "lacks interval columns" + ) +}) + +test_that("mixedmirt group-variance extractor covers COV_ and fail-closed shapes", { + .ensure_kaefa_namespace() + testthat::skip_if_not_installed("mirt") + + with_ci <- .fake_mixedclass( + item_pars = list(c(a1 = 1, d = 0.6)), + item_se = list(c(0.05, 0.1)) + ) + tau <- kaefa:::.extractMixedmirtGroupVariance(with_ci, "group") + testthat::expect_equal(tau$estimate, 0.5) + testthat::expect_equal(tau$name, "COV_group") + testthat::expect_true(is.finite(tau$se)) + testthat::expect_true(is.finite(tau$lower)) + testthat::expect_true(is.finite(tau$upper)) + + no_ci <- .fake_mixedclass( + item_pars = list(c(a1 = 1, d = 0.6)), + item_se = NULL, + random_se = NULL, + group_se = NULL + ) + tau_na <- kaefa:::.extractMixedmirtGroupVariance(no_ci, "group") + testthat::expect_equal(tau_na$estimate, 0.5) + testthat::expect_true(is.na(tau_na$se)) + testthat::expect_true(is.na(tau_na$lower)) + testthat::expect_true(is.na(tau_na$upper)) + + no_random <- .fake_mixedclass( + item_pars = list(c(a1 = 1, d = 0.6)), + item_se = list(c(0.05, 0.1)), + include_random = FALSE + ) + testthat::expect_error( + kaefa:::.extractMixedmirtGroupVariance(no_random, "group"), + "do not contain the 'group' random effect" + ) + testthat::expect_error( + kaefa:::.extractMixedmirtGroupVariance(with_ci, "school"), + "do not contain the 'school' random effect" + ) + testthat::expect_error( + kaefa:::.extractMixedmirtGroupVariance(with_ci, ""), + "single non-empty string" + ) + testthat::expect_error( + kaefa:::.extractMixedmirtGroupVariance(with_ci, c("group", "school")), + "single non-empty string" + ) + + no_cov <- .fake_mixedclass( + item_pars = list(c(a1 = 1, d = 0.6)), + item_se = list(c(0.05, 0.1)), + random_par = 0.5, + random_parnames = "THETA_var", + random_se = 0.08 + ) + testthat::expect_error( + kaefa:::.extractMixedmirtGroupVariance(no_cov, "group"), + "no COV_ variance column" + ) + + crossed <- .fake_mixedclass( + item_pars = list(c(a1 = 1, d = 0.6)), + item_se = list(c(0.05, 0.1)), + random_par = c(0.5, 0.2), + random_parnames = c("COV_group", "COV_rater"), + random_se = c(0.08, 0.05) + ) + testthat::expect_error( + kaefa:::.extractMixedmirtGroupVariance(crossed, "group"), + "Crossed or multiple random-effect variances" + ) +}) + +test_that("mixedmirt par-row helper is fail-closed without a par row", { + .ensure_kaefa_namespace() + with_par <- matrix(c(1, -0.5), nrow = 1, dimnames = list("par", c("a", "b"))) + testthat::expect_identical(kaefa:::.mixedmirtParRow(with_par), "par") + testthat::expect_error( + kaefa:::.mixedmirtParRow(matrix(1, dimnames = list("estimate", "a"))), + "no 'par' row" + ) + testthat::expect_error( + kaefa:::.mixedmirtParRow(matrix(1)), + "no 'par' row" + ) +}) + +test_that("nested mixedmirt fit wrapper rejects incomplete designs and stays nested", { + .ensure_kaefa_namespace() + testthat::skip_if_not_installed("mirt") + + testthat::expect_error( + kaefa:::.fitMixedmirtNestedRecovery(list(response_data = data.frame(Item1 = 0:1))), + "mixedmirtNestedRecoveryDesign" + ) + testthat::expect_error( + kaefa:::.fitMixedmirtNestedRecovery(list(covdata = data.frame(group = 1:40))), + "mixedmirtNestedRecoveryDesign" + ) + testthat::expect_error( + kaefa:::.rejectAtomisticMixedmirtRecovery(data.frame(group = integer()), ~ 1 | group), + "atomistic fallacy" + ) + testthat::expect_error( + kaefa:::.rejectAtomisticMixedmirtRecovery( + data.frame(group = 1:40), + stats::as.formula("~ 1 | group | rater") + ), + "Crossed or interaction" + ) + + design <- kaefa:::.mixedmirtNestedRecoveryDesign( + seed = 20260817L, + n_groups = 30L, + n_per_group = 5L + ) + testthat::with_mocked_bindings( + .mixedmirt = function(...) methods::new("MixedClass"), + .package = "kaefa", + { + fit <- kaefa:::.fitMixedmirtNestedRecovery(design) + testthat::expect_true(methods::is(fit, "MixedClass")) + } + ) +}) + +test_that("mixedmirt coverage log covers nested RE and excludes MM/time-flow", { + .ensure_kaefa_namespace() + coverage <- kaefa:::.recoveryCoverageExclusions() + testthat::expect_identical(names(coverage), c("surface", "status", "reason")) + testthat::expect_true(any( + grepl("nested two-level random intercept", coverage$surface) & + coverage$status == "covered" + )) + testthat::expect_true(any( + grepl("lr.random", coverage$surface) & coverage$status == "excluded" + )) + testthat::expect_true(any( + grepl("multiple-membership", coverage$surface) & coverage$status == "excluded" + )) + testthat::expect_true(any( + grepl("crossed random effects", coverage$surface) & coverage$status == "excluded" + )) + testthat::expect_true(any( + grepl("time-flow", coverage$surface) & coverage$status == "excluded" + )) + testthat::expect_false(any( + grepl("mixedmirt nested", coverage$surface) & coverage$status == "excluded" + )) +}) diff --git a/tests/testthat/test-mixedmirt-recovery-fits.R b/tests/testthat/test-mixedmirt-recovery-fits.R new file mode 100644 index 0000000..41c4cf5 --- /dev/null +++ b/tests/testthat/test-mixedmirt-recovery-fits.R @@ -0,0 +1,116 @@ +.mixedmirt_recovery_once <- function(seed) { + design <- kaefa:::.mixedmirtNestedRecoveryDesign(seed = seed) + utils::capture.output( + fit <- suppressWarnings(kaefa:::.fitMixedmirtNestedRecovery(design)) + ) + if (!methods::is(fit, "MixedClass")) { + testthat::fail( + paste( + "kaefa::.mixedmirt did not return a MixedClass fit for seed", + seed, + "(converged second-order test required)" + ) + ) + return(invisible(NULL)) + } + + estimated <- kaefa:::.extractMixedmirtIrtItems(fit) + aligned <- kaefa:::.alignIrtItemParameters(estimated, design$truth_items) + tau <- kaefa:::.extractMixedmirtGroupVariance(fit, design$group_name) + item_intervals <- kaefa:::.extractMixedmirtItemIntervals(fit, aligned$items) + item_inclusion <- kaefa:::.recoveryIntervalInclusion( + estimate = item_intervals$estimate, + truth = aligned$truth$b, + lower = item_intervals$lower, + upper = item_intervals$upper + ) + tau_inclusion <- NULL + if (is.finite(tau$lower) && is.finite(tau$upper)) { + tau_inclusion <- kaefa:::.recoveryIntervalInclusion( + estimate = tau$estimate, + truth = design$truth_tau00, + lower = tau$lower, + upper = tau$upper + ) + } + + list( + design = design, + aligned = aligned, + rmse_b = kaefa:::.parameterRecoveryRMSE(aligned$estimated$b, aligned$truth$b), + rmse_tau00 = kaefa:::.parameterRecoveryRMSE(tau$estimate, design$truth_tau00), + item_inclusion = item_inclusion, + tau_inclusion = tau_inclusion, + tau = tau + ) +} + +test_that(".mixedmirt recovers nested Rasch b and tau00 with bounded RMSE", { + .ensure_kaefa_namespace() + testthat::skip_if_not_installed("mirt") + + recovered <- .mixedmirt_recovery_once(20260817L) + if (is.null(recovered)) { + return(invisible(NULL)) + } + + testthat::expect_lt(recovered$rmse_b, 0.35) + testthat::expect_lt(recovered$rmse_tau00, 0.35) + testthat::expect_gt(mean(recovered$item_inclusion$covers), 0.60) + testthat::expect_true(is.finite(recovered$tau$estimate)) + if (is.null(recovered$tau_inclusion)) { + testthat::fail("mixedmirt group-variance Wald interval was not estimable") + } else { + testthat::expect_true(is.finite(recovered$tau_inclusion$lower[1])) + testthat::expect_true(is.finite(recovered$tau_inclusion$upper[1])) + } +}) + +test_that("mixedmirt five-seed nested recovery reports mean RMSE", { + .skip_expensive_ci_calls("mixedmirt") + .ensure_kaefa_namespace() + testthat::skip_if_not_installed("mirt") + + seeds <- c(20260817L, 20260818L, 20260819L, 20260820L, 20260821L) + rows <- list() + item_cover_n <- 0L + item_interval_n <- 0L + tau_cover_n <- 0L + tau_interval_n <- 0L + + for (seed in seeds) { + recovered <- .mixedmirt_recovery_once(seed) + if (is.null(recovered)) { + return(invisible(NULL)) + } + rows[[length(rows) + 1L]] <- data.frame( + seed = seed, + parameter = c("b", "tau00"), + rmse = c(recovered$rmse_b, recovered$rmse_tau00), + stringsAsFactors = FALSE + ) + item_cover_n <- item_cover_n + sum(recovered$item_inclusion$covers) + item_interval_n <- item_interval_n + nrow(recovered$item_inclusion) + if (!is.null(recovered$tau_inclusion)) { + tau_cover_n <- tau_cover_n + sum(recovered$tau_inclusion$covers) + tau_interval_n <- tau_interval_n + nrow(recovered$tau_inclusion) + } + } + + summarised <- kaefa:::.summariseRecoveryRepeats(do.call(rbind, rows)) + testthat::expect_identical(summarised$n_repeats, 5L) + testthat::expect_lt( + summarised$summary$mean_rmse[summarised$summary$parameter == "b"], + 0.35 + ) + testthat::expect_lt( + summarised$summary$mean_rmse[summarised$summary$parameter == "tau00"], + 0.35 + ) + testthat::expect_gt(item_interval_n, 0L) + testthat::expect_gt(item_cover_n / item_interval_n, 0.60) + testthat::expect_identical(tau_interval_n, 5L) + # Five seeds cannot support a 95% coverage-rate claim (Harwell et al., 1996). + # Require only that the group-variance interval was estimable on every seed. + testthat::expect_true(tau_cover_n >= 0L) +})