From 30a7e0a4af7c3c6d9b62172bd4462137e24e166d Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Mon, 13 Jul 2026 16:59:25 +0900 Subject: [PATCH 1/2] fix(aefa): preserve valid model-selection formulas --- R/kaefa-package.r | 10 +- R/kaefa.R | 362 ++++++++++++++---- README.Rmd | 4 +- README.md | 12 +- docs/papers/README.md | 41 +- inst/shiny-app/README.md | 6 +- inst/shiny-app/app.R | 6 +- man/aefa.Rd | 9 +- man/kaefa.Rd | 10 +- tests/testthat/helper-test-data.R | 11 +- .../testthat/test-aefa-advanced-parameters.R | 19 +- tests/testthat/test-aefa-greedy-algorithm.R | 4 +- .../testthat/test-model-selection-criteria.R | 66 ++++ tests/testthat/test-shiny-app.R | 2 +- tests/testthat/test-test-bootstrap-contract.R | 14 + 15 files changed, 476 insertions(+), 100 deletions(-) create mode 100644 tests/testthat/test-model-selection-criteria.R diff --git a/R/kaefa-package.r b/R/kaefa-package.r index 98c8f6a..cfbb49c 100644 --- a/R/kaefa-package.r +++ b/R/kaefa-package.r @@ -22,10 +22,12 @@ #' #' The kaefa may inspect this issues from the MMMM or MM #' in statistical learning theory perspectives using model -#' selection criteria like the DIC (Kang, 2008; Kang, Cohen, & Sung, 2009; -#' Jiao, Kamata, Wang, & Jin, 2012; Jiao & Zhang, 2015) with maximising -#' generalisability of the number of factor decisions in every calibration -#' (Kang, 2008; Preacher, Zhang, Kim, & Mels, 2013). +#' likelihood-based model-selection criteria such as AIC, AICc, BIC, and +#' sample-size-adjusted BIC. Posterior DIC (Spiegelhalter et al., 2002) is +#' accepted only when the fitted model actually supplies it; DIC is not +#' reconstructed from AIC. These criteria support generalisable factor-count +#' decisions in each calibration (Kang, 2008; Preacher, Zhang, Kim, & Mels, +#' 2013). #' #' If researcher provide of demographical information in kaefa, #' kaefa will inspect the optimal number of factor and optimal IRT model, diff --git a/R/kaefa.R b/R/kaefa.R index 3e4b134..e0ab631 100644 --- a/R/kaefa.R +++ b/R/kaefa.R @@ -1,5 +1,24 @@ # kaefa.R +#' Choose a bounded local worker count for AEFA +#' +#' @param cpu_idle Optional scalar percentage of idle CPU capacity. +#' @param physical_cores Detected physical core count. +#' @return A positive integer worker count. +#' @keywords internal +#' @noRd +.aefaParallelProcessorCount <- function(cpu_idle, physical_cores) { + if (length(physical_cores) != 1L || !is.numeric(physical_cores) || + !is.finite(physical_cores) || physical_cores < 1) { + return(1L) + } + if (length(cpu_idle) == 1L && is.numeric(cpu_idle) && + is.finite(cpu_idle) && cpu_idle < 10) { + return(1L) + } + as.integer(max(1, round(physical_cores / 3))) +} + #' Initalize aefa engine, #' This function initalise the aefa cluster. #' If someone have Remote Cluster informaiton with SSH, put the information in the argument. @@ -214,52 +233,38 @@ aefaInit <- function(RemoteClusters = getOption("kaefaServers"), debug = F, sshK # connList <- connList[sample(x = 1:length(connList), size = length(connList))] return(connList) } - - - - if (is.null(suppressWarnings(NCmisc::top()$CPU$idle))) { - parallelProcessors <- round(parallel::detectCores(all.tests = FALSE, logical = FALSE)/3) - if (1 >= parallelProcessors) { - parallelProcessors <- 1 - } - } else if (suppressWarnings(NCmisc::top()$CPU$idle) > 50) { - parallelProcessors <- round(parallel::detectCores(all.tests = FALSE, logical = FALSE)/3) - if (1 >= parallelProcessors) { - parallelProcessors <- 1 - } - } else if (suppressWarnings(NCmisc::top()$CPU$idle) <= 50) { - parallelProcessors <- round(parallel::detectCores(all.tests = FALSE, logical = FALSE)/3) - if (1 >= parallelProcessors) { - parallelProcessors <- 1 + cpuIdle <- tryCatch( + suppressWarnings(NCmisc::top()$CPU$idle), + error = function(e) { + if (isTRUE(debug)) { + message("CPU idle sampling unavailable; using detected core count: ", + conditionMessage(e)) } - } else if (suppressWarnings(NCmisc::top()$CPU$idle) < 30) { - parallelProcessors <- round(parallel::detectCores(all.tests = FALSE, logical = FALSE)/3) - if (1 >= parallelProcessors) { - parallelProcessors <- 1 - } - } else if (suppressWarnings(NCmisc::top()$CPU$idle) < 10) { - parallelProcessors <- 1 - } + NULL + } + ) + physicalCores <- parallel::detectCores(all.tests = FALSE, logical = FALSE) + parallelProcessors <- .aefaParallelProcessorCount(cpuIdle, physicalCores) # setting up cluster if (!is.null(RemoteClusters)) { halfCores <- function() { max(1, round(0.3 * future::availableCores()))} try(future::plan(list( future::tweak(future::cluster, workers = assignClusterNodes(RemoteClusters), gc = T, homogeneous = FALSE, earlySignal = T), - future::tweak(future::multiprocess, workers = halfCores, gc = T, earlySignal = T) + future::tweak(future::multisession, workers = halfCores, gc = T, earlySignal = T) ))) # try(future::plan(future::tweak("future::cluster", # workers = tryCatch(assignClusterNodes(RemoteClusters), error = function(e){assignClusterNodes(RemoteClusters)}), gc = T # ))) } else if (NROW(future::plan("list")) == 1) { if (length(grep("openblas|microsoft", extSoftVersion()["BLAS"])) > 0) { - options(aefaConn = future::plan("future::multiprocess", workers = parallelProcessors), + options(aefaConn = future::plan(future::multisession, workers = parallelProcessors), gc = T) } else if (parallel::detectCores(logical = F) == 1) { options(aefaConn = future::plan(future::sequential), gc = T) } else { options(aefaConn = (tryCatch(future::plan(strategy = list(future::tweak(future::cluster(workers = parallelProcessors)), - future::multiprocess(workers = parallelProcessors)), gc = T), error = function(e) { + future::tweak(future::multisession, workers = parallelProcessors)), gc = T), error = function(e) { }))) } } @@ -432,6 +437,142 @@ evaluateItemFit <- function(mirtModel, RemoteClusters = NULL, rotate = "bifactor } } +#' Extract a scientifically valid information criterion from a mirt fit +#' +#' DIC is a posterior-deviance criterion (Spiegelhalter et al., 2002) and +#' therefore cannot be reconstructed from a maximum-likelihood AIC value. +#' AICc is reconstructed only from its original small-sample correction +#' (Hurvich and Tsai, 1989): AIC + 2 k (k + 1) / (n - k - 1). +#' +#' @param fit The model's named fit-statistic list. +#' @param criterion Requested information criterion. +#' @param sample_size Number of response patterns used to fit the model. +#' @return One finite information-criterion value. +#' @keywords internal +#' @noRd +.aefaFitCriterionValue <- function(fit, criterion, sample_size) { + criterion <- toupper(criterion) + fit_value <- function(name) { + value <- fit[[name]] + if (length(value) == 1L && is.numeric(value) && is.finite(value)) { + return(as.numeric(value)) + } + NULL + } + + if (criterion == "DIC") { + value <- fit_value("DIC") + if (is.null(value)) { + stop( + paste0( + "DIC is unavailable: DIC requires posterior deviance draws and an ", + "effective parameter count, but this mirt ML/MAP fit does not supply ", + "a DIC value. Use AIC, AICc, BIC, or saBIC; DIC is never replaced by AIC." + ), + call. = FALSE + ) + } + return(value) + } + + if (criterion == "CAIC") { + stop( + "CAIC is ambiguous and is not an alias for corrected AIC; use 'AICc' explicitly.", + call. = FALSE + ) + } + + if (criterion == "AICC") { + value <- fit_value("AICc") + if (!is.null(value)) { + return(value) + } + + aic <- fit_value("AIC") + log_likelihood <- fit_value("logLik") + if (is.null(aic) || is.null(log_likelihood)) { + stop( + "AICc is unavailable because the fit supplies neither AICc nor both AIC and logLik.", + call. = FALSE + ) + } + if (length(sample_size) != 1L || !is.numeric(sample_size) || + !is.finite(sample_size) || sample_size <= 0) { + stop("AICc requires one finite positive sample size.", call. = FALSE) + } + + parameter_count <- (aic + 2 * log_likelihood) / 2 + if (!is.finite(parameter_count) || parameter_count < 0) { + stop("AICc could not recover a valid parameter count from AIC and logLik.", call. = FALSE) + } + if (sample_size <= parameter_count + 1) { + stop( + paste0( + "AICc is undefined because n (", sample_size, + ") must be greater than k + 1 (", parameter_count + 1, ")." + ), + call. = FALSE + ) + } + + return(aic + (2 * parameter_count * (parameter_count + 1)) / + (sample_size - parameter_count - 1)) + } + + fit_name <- switch( + criterion, + AIC = "AIC", + BIC = "BIC", + SABIC = "SABIC", + NULL + ) + if (is.null(fit_name)) { + stop( + "Unsupported model selection criterion. Use AIC, AICc, BIC, saBIC, or model-supplied DIC.", + call. = FALSE + ) + } + + value <- fit_value(fit_name) + if (is.null(value)) { + stop( + paste0(criterion, " is unavailable because the fitted model does not supply a finite value."), + call. = FALSE + ) + } + value +} + +#' Select the original candidate index with the lowest finite score +#' +#' @param scores Numeric vector aligned one-to-one with model candidates. +#' @return The original candidate index, or code{NA_integer_} when none is valid. +#' @keywords internal +#' @noRd +.aefaBestScoreIndex <- function(scores) { + finite_indices <- which(is.finite(scores)) + if (length(finite_indices) == 0L) { + return(NA_integer_) + } + finite_indices[[which.min(scores[finite_indices])]] +} + +#' Return proposed item exclusions that can make calibration progress +#' +#' @param excluded Item names already excluded from calibration. +#' @param proposed Item names proposed by the current fit diagnostic. +#' @param available Item names present in the fitted model. +#' @return Unique, available item names not previously excluded. +#' @keywords internal +#' @noRd +.aefaValidNewItemExclusions <- function(excluded, proposed, available) { + excluded <- unique(as.character(excluded)) + proposed <- unique(as.character(proposed)) + available <- unique(as.character(available)) + proposed <- proposed[!is.na(proposed) & nzchar(proposed)] + setdiff(intersect(proposed, available), excluded) +} + #' doing automated exploratory factor analysis (aefa) for research capability to identify unexplained factor structure with complexly cross-classified multilevel structured data in R environment #' #' This function implements a greedy search algorithm to efficiently explore the model space and find improved model configurations. @@ -464,7 +605,12 @@ evaluateItemFit <- function(mirtModel, RemoteClusters = NULL, rotate = "bifactor #' @param resampling Do you want to do resampling with replace? default is TRUE and activate nrow is over samples argument. #' @param samples Specify the number samples with resampling. default is 5000. #' @param printDebugMsg Do you want to see the debugging messeages? default is FALSE -#' @param modelSelectionCriteria Which critera want to use model selection work? 'DIC' (default), 'AIC', 'AICc', 'BIC', 'saBIC' available. AIC and DIC will be identical when no prior parameter distributions are included. +#' @param modelSelectionCriteria Information criterion used for model selection. +#' code{"AIC"} is the default for the maximum-likelihood/MAP models fitted by +#' current versions of pkg{mirt}. code{"AICc"}, code{"BIC"}, and +#' code{"saBIC"} are also available. code{"DIC"} is accepted only when the +#' fitted model actually supplies a posterior-based DIC value; it is never +#' substituted with AIC. #' @param saveRawEstModels Do you want to save raw estimated models before model selection work? default is FALSE #' @param fitEMatUIRT Do you want to fit the model with EM at UIRT? default is FALSE #' @param ranefautocomb Do you want to find global-optimal random effect combination? default is TRUE @@ -514,7 +660,7 @@ aefa <- efa <- function(data, model = NULL, minExtraction = 1, maxExtraction = i filename = "aefa.RDS", printItemFit = T, rotate = c("bifactorQ","geominQ", "geominT", "bentlerQ", "bentlerT", "oblimin", "simplimax", "tandemII", "tandemI", "entropy", "quartimax"), resampling = T, samples = 5000, - printDebugMsg = F, modelSelectionCriteria = "DIC", saveRawEstModels = F, fitEMatUIRT = F, + printDebugMsg = F, modelSelectionCriteria = "AIC", saveRawEstModels = F, fitEMatUIRT = F, ranefautocomb = T, PV_Q1 = T, tryLCA = F, forcingQMC = F, turnOffMixedEst = F, fitIndicesCutOff = 0.005, anchor = colnames(data), skipggum = F, powertest = F, idling = 0, leniency = F) { @@ -623,13 +769,19 @@ aefa <- efa <- function(data, model = NULL, minExtraction = 1, maxExtraction = i tryCatch(rm(estModel), error = function(e) {NULL}) } modelDONE <- FALSE + modelAttempt <- 0L + lastModelError <- "engineAEFA returned no model" while (!modelDONE) { + modelAttempt <- modelAttempt + 1L if(workDirectory != getwd()){ setwd('/tmp') setwd(workDirectory) } tryCatch(aefaInit(RemoteClusters = RemoteClusters, debug = printDebugMsg, - sshKeyPath = sshKeyPath), error = function(e) {NULL}) + sshKeyPath = sshKeyPath), error = function(e) { + message("AEFA cluster initialization was unavailable; continuing locally: ", + conditionMessage(e)) + }) # general condition estModel <- tryCatch(engineAEFA(data = data.frame(data[, !colnames(data) %in% badItemNames]), model = model, GenRandomPars = GenRandomPars, NCYCLES = NCYCLES, @@ -639,9 +791,24 @@ aefa <- efa <- function(data, model = NULL, minExtraction = 1, maxExtraction = i fitEMatUIRT = fitEMatUIRT, ranefautocomb = ranefautocomb, tryLCA = tryLCA, forcingQMC = forcingQMC, turnOffMixedEst = turnOffMixedEst, anchor = anchor[!anchor %in% DIFitems], skipggumInternal = skipggum, powertest = powertest, idling = idling, leniency = leniency), error = function(e) { + lastModelError <<- conditionMessage(e) + NULL }) - if (exists("estModel")) { + if (!is.null(estModel)) { modelDONE <- TRUE + } else if (modelAttempt >= 3L) { + stop( + paste0( + "AEFA model estimation failed after ", modelAttempt, + " attempts. Last error: ", lastModelError + ), + call. = FALSE + ) + } else { + message( + "AEFA model estimation attempt ", modelAttempt, + " failed; retrying. Reason: ", lastModelError + ) } } @@ -714,47 +881,36 @@ aefa <- efa <- function(data, model = NULL, minExtraction = 1, maxExtraction = i if (class(estModel) == "list") { # model fit evaluation - modModelFit <- vector() + # Keep scores aligned with their original candidates. A skipped + # Heywood solution must not shift the selected model index. + modModelFit <- rep(Inf, NROW(estModel)) for (i in 1:NROW(estModel)) { if (sum(c("MixedClass", "SingleGroupClass", "DiscreteClass", "MultipleGroupClass") %in% class(estModel[[i]])) > 0) { # heywood case filter if(class(estModel[[i]]) %in% "MixedClass"){ - if(sum(invisible(.exportParmsEME(estModel[[i]], quiet = T))@Fit$h2 > 1) > 0){ + if(any(invisible(.exportParmsEME(estModel[[i]], quiet = T))@Fit$h2 > 1, + na.rm = TRUE)){ next() } } else if(class(estModel[[i]]) %in% "SingleGroupClass"){ - if(sum(estModel[[i]]@Fit$h2 > 1) > 0){ + if(any(estModel[[i]]@Fit$h2 > 1, na.rm = TRUE)){ next() } } - if (toupper(modelSelectionCriteria) %in% c("DIC")) { - modModelFit[[length(modModelFit) + 1]] <- estModel[[i]]@Fit$DIC - } else if (toupper(modelSelectionCriteria) %in% c("AIC")) { - modModelFit[[length(modModelFit) + 1]] <- estModel[[i]]@Fit$AIC - } else if (toupper(modelSelectionCriteria) %in% c("AICC", "CAIC")) { - modModelFit[[length(modModelFit) + 1]] <- estModel[[i]]@Fit$AICc - } else if (toupper(modelSelectionCriteria) %in% c("BIC")) { - modModelFit[[length(modModelFit) + 1]] <- estModel[[i]]@Fit$BIC - } else if (toupper(modelSelectionCriteria) %in% c("SABIC")) { - modModelFit[[length(modModelFit) + 1]] <- estModel[[i]]@Fit$SABIC - } else { - stop("please specify model fit type correctly: DIC (default), AIC, BIC, AICc (aka cAIC), saBIC") - } + modModelFit[[i]] <- .aefaFitCriterionValue( + fit = estModel[[i]]@Fit, + criterion = modelSelectionCriteria, + sample_size = nrow(estModel[[i]]@Data$data) + ) } } # select model - if (exists("modModelFit")) { - if (length(which(modModelFit == min(modModelFit[is.finite(modModelFit)], - na.rm = T))[1]) > 0) { - estModel <- estModel[[which(modModelFit == min(modModelFit[is.finite(modModelFit)], - na.rm = T))[1]]] - } else { - message("Can not find any optimal model") - STOP <- T - } + selectedModelIndex <- .aefaBestScoreIndex(modModelFit) + if (!is.na(selectedModelIndex)) { + estModel <- estModel[[selectedModelIndex]] } else { message("Can not find any optimal model") STOP <- T @@ -776,7 +932,7 @@ aefa <- efa <- function(data, model = NULL, minExtraction = 1, maxExtraction = i } } } - # save model history of DIC evaluated model + # save the selected model history if (saveModelHistory) { if(workDirectory != getwd()){ setwd('/tmp') @@ -792,7 +948,10 @@ aefa <- efa <- function(data, model = NULL, minExtraction = 1, maxExtraction = i }) } fitDONE <- FALSE + fitAttempt <- 0L + lastFitError <- "evaluateItemFit returned no fit statistics" while (!fitDONE) { + fitAttempt <- fitAttempt + 1L if(workDirectory != getwd()){ setwd('/tmp') setwd(workDirectory) @@ -893,7 +1052,10 @@ aefa <- efa <- function(data, model = NULL, minExtraction = 1, maxExtraction = i estItemFit <- estItemFitRotationSearchTmp[[paste0(rotateCandidates)]] } else { estItemFit <- tryCatch(evaluateItemFit(estModel, RemoteClusters = RemoteClusters, - rotate = rotateCandidates, PV_Q1 = PV_Q1), error = function(e) {NULL}) + rotate = rotateCandidates, PV_Q1 = PV_Q1), error = function(e) { + lastFitError <<- conditionMessage(e) + NULL + }) } @@ -901,7 +1063,10 @@ aefa <- efa <- function(data, model = NULL, minExtraction = 1, maxExtraction = i # estimate item fit measures rotateCandidates <- 'none' estItemFit <- tryCatch(evaluateItemFit(estModel, RemoteClusters = RemoteClusters, - rotate = rotateCandidates, PV_Q1 = PV_Q1), error = function(e) {NULL}) + rotate = rotateCandidates, PV_Q1 = PV_Q1), error = function(e) { + lastFitError <<- conditionMessage(e) + NULL + }) } if (exists("estItemFit")) { @@ -909,6 +1074,20 @@ aefa <- efa <- function(data, model = NULL, minExtraction = 1, maxExtraction = i fitDONE <- TRUE } } + if (!fitDONE && fitAttempt >= 3L) { + stop( + paste0( + "AEFA item-fit evaluation failed after ", fitAttempt, + " attempts. Last error: ", lastFitError + ), + call. = FALSE + ) + } else if (!fitDONE) { + message( + "AEFA item-fit evaluation attempt ", fitAttempt, + " failed; retrying. Reason: ", lastFitError + ) + } } if (printItemFit) { tryCatch(print(estItemFit), error = function(e) {NULL}) @@ -982,8 +1161,22 @@ aefa <- efa <- function(data, model = NULL, minExtraction = 1, maxExtraction = i S_X2Cond4 <- FALSE } - # flagging bad item - if (ZhCond) { + itemFitRemovalRequired <- any(c( + ZhCond, PVCond1, PVCond2, PVCond3, + S_X2Cond1, S_X2Cond2, S_X2Cond3, S_X2Cond4 + )) + badItemNamesBefore <- unique(as.character(badItemNames)) + + # A factor model cannot keep deleting indicators indefinitely. + # Retain the current fitted model at the established three-item + # floor and make the residual misfit visible in the log. + if (itemFitRemovalRequired && length(estItemFit$item) <= 3L) { + message( + "AEFA item exclusion stopped at the three-item floor; ", + "the retained model still has item-fit diagnostics above the configured cutoff." + ) + STOP <- TRUE + } else if (ZhCond) { badItemNames <- c(badItemNames, as.character(estItemFit$item[which(estItemFit$Zh == min(estItemFit$Zh[is.finite(estItemFit$Zh)], na.rm = T))])) } else if (PVCond1) { @@ -1009,8 +1202,15 @@ aefa <- efa <- function(data, model = NULL, minExtraction = 1, maxExtraction = i max(estItemFit$RMSEA.S_X2[is.finite(estItemFit$RMSEA.S_X2)], na.rm = T))])) } else if(class(estModel) %in% "MultipleGroupClass" && checkDIF){ # DIF part (fixed effect) - if (toupper(modelSelectionCriteria) %in% c("DIC")) { # find out - seqStat <- 'AIC' + if (toupper(modelSelectionCriteria) %in% c("DIC")) { + stop( + paste0( + "DIC cannot be substituted with AIC during sequential DIF selection; ", + "mirt::DIF does not provide a posterior-DIC sequence statistic. ", + "Use AIC, AICc, BIC, or saBIC." + ), + call. = FALSE + ) } else if (toupper(modelSelectionCriteria) %in% c("AIC")) { seqStat <- 'AIC' } else if (toupper(modelSelectionCriteria) %in% c("AICC", "CAIC")) { @@ -1081,6 +1281,40 @@ aefa <- efa <- function(data, model = NULL, minExtraction = 1, maxExtraction = i STOP <- TRUE } + if (itemFitRemovalRequired && !STOP) { + proposedBadItems <- setdiff( + unique(as.character(badItemNames)), + badItemNamesBefore + ) + availableItems <- colnames(estModel@Data$data) + validNewBadItems <- .aefaValidNewItemExclusions( + excluded = badItemNamesBefore, + proposed = proposedBadItems, + available = availableItems + ) + invalidProposals <- setdiff(proposedBadItems, availableItems) + badItemNames <- unique(c(badItemNamesBefore, validNewBadItems)) + + if (length(invalidProposals) > 0L) { + message( + "AEFA item diagnostic proposed names absent from the fitted model: ", + paste(invalidProposals, collapse = ", ") + ) + } + if (length(validNewBadItems) == 0L) { + message( + "AEFA item exclusion stopped because the diagnostic did not identify ", + "a new available item; repeated calibration would make no progress." + ) + STOP <- TRUE + } else { + message( + "AEFA item-fit exclusion: ", + paste(validNewBadItems, collapse = ", ") + ) + } + } + # adjust model if supplied model is confirmatory model if (!is.null(model) && (!is.numeric(model) | !is.integer(model)) && "Zh" %in% colnames(estItemFit)) { diff --git a/README.Rmd b/README.Rmd index 271af9a..f646ed6 100644 --- a/README.Rmd +++ b/README.Rmd @@ -22,7 +22,7 @@ The goal of kaefa is to improve researchers' ability to identify unexplained fac The automated exploratory factor analysis (aefa) framework implements a **greedy search algorithm** to efficiently explore the model space and find improved model configurations. The algorithm iteratively: 1. Evaluates multiple model candidates with different factor structures and item response models -2. Selects the best model based on information criteria (DIC, AIC, BIC, etc.) +2. Selects the best model using AIC by default, with AICc, BIC, and sample-size-adjusted BIC available. DIC is used only when a fitted model actually supplies posterior DIC; it is never approximated with AIC. 3. Assesses item fit and removes poorly fitting items one at a time 4. Re-estimates the model until convergence to a locally optimal solution @@ -32,6 +32,8 @@ This greedy approach enables efficient exploration of the model space while seek - Preacher, K. J., Zhang, G., Kim, C., & Mels, G. (2013). Choosing the optimal number of factors in exploratory factor analysis: A model selection perspective. Multivariate Behavioral Research, 48(1), 28-56. https://doi.org/10.1080/00273171.2012.710386 - Jennrich, R. I., & Bentler, P. M. (2011). Exploratory bi-factor analysis. Psychometrika, 76(4), 537-549. https://doi.org/10.1007/s11336-011-9218-4 +- Hurvich, C. M., & Tsai, C.-L. (1989). Regression and time series model selection in small samples. Biometrika, 76(2), 297-307. https://doi.org/10.1093/biomet/76.2.297 +- Spiegelhalter, D. J., Best, N. G., Carlin, B. P., & van der Linde, A. (2002). Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society: Series B, 64(4), 583-639. https://doi.org/10.1111/1467-9868.00353 ## Installation diff --git a/README.md b/README.md index 18c9164..76d0386 100644 --- a/README.md +++ b/README.md @@ -18,8 +18,9 @@ find improved model configurations. The algorithm iteratively: 1. Evaluates multiple model candidates with different factor structures and item response models -2. Selects the best model based on information criteria (DIC, AIC, BIC, - etc.) +2. Selects the best model using AIC by default, with AICc, BIC, and + sample-size-adjusted BIC available. DIC is used only when a fitted model + actually supplies posterior DIC; it is never approximated with AIC. 3. Assesses item fit and removes poorly fitting items one at a time 4. Re-estimates the model until convergence to a locally optimal solution @@ -38,6 +39,13 @@ research (Preacher, Zhang, Kim, & Mels, 2013; Jennrich & Bentler, 2011). - Jennrich, R. I., & Bentler, P. M. (2011). Exploratory bi-factor analysis. Psychometrika, 76(4), 537-549. +- Hurvich, C. M., & Tsai, C.-L. (1989). Regression and time series model + selection in small samples. Biometrika, 76(2), 297-307. + +- Spiegelhalter, D. J., Best, N. G., Carlin, B. P., & van der Linde, A. + (2002). Bayesian measures of model complexity and fit. Journal of the + Royal Statistical Society: Series B, 64(4), 583-639. + ## Installation diff --git a/docs/papers/README.md b/docs/papers/README.md index 86be622..cc2caa8 100644 --- a/docs/papers/README.md +++ b/docs/papers/README.md @@ -86,10 +86,43 @@ cited with its DOI. Open-access / preprint links are noted where available. Both have expectation 1 under good fit. +## 5. Corrected Akaike information criterion (`AICc`) + +- **Source:** Hurvich, C. M., & Tsai, C.-L. (1989). Regression and time series + model selection in small samples. *Biometrika, 76*(2), 297-307. + DOI: [10.1093/biomet/76.2.297](https://doi.org/10.1093/biomet/76.2.297) +- **Canonical equation.** For maximized log-likelihood `logLik`, parameter count + `k`, and sample size `n`: + + AIC = -2 * logLik + 2 * k + AICc = AIC + 2 * k * (k + 1) / (n - k - 1) + + Current `mirt` fits supply `AIC` and `logLik` but do not consistently expose + `AICc`. Kaefa therefore recovers `k = (AIC + 2 * logLik) / 2` and applies the + correction exactly. The statistic is undefined when `n <= k + 1`; kaefa + reports that reason rather than returning a fabricated finite score. + +## 6. Deviance information criterion (`DIC`) boundary + +- **Source:** Spiegelhalter, D. J., Best, N. G., Carlin, B. P., & van der Linde, + A. (2002). Bayesian measures of model complexity and fit. *Journal of the + Royal Statistical Society: Series B, 64*(4), 583-639. + DOI: [10.1111/1467-9868.00353](https://doi.org/10.1111/1467-9868.00353) +- **Canonical equation.** With posterior mean deviance `Dbar`, deviance at the + posterior mean parameters `D(theta_bar)`, and effective parameter count `pD`: + + pD = Dbar - D(theta_bar) + DIC = Dbar + pD + + DIC is a posterior-deviance criterion. The maximum-likelihood/MAP models + produced by current `mirt` versions do not expose the posterior quantities + needed to reconstruct it. Kaefa accepts DIC only when the fitted model + supplies a finite DIC value and never relabels AIC as DIC. + ## 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 are therefore correct by -construction (subject to `mirt`'s own validation). The only package-local numeric -formulas are the fit-based **decision rules** documented above; the `Zh` rule is -the one that had drifted out of internal consistency and has been restored. +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. diff --git a/inst/shiny-app/README.md b/inst/shiny-app/README.md index bac2581..e4c5f49 100644 --- a/inst/shiny-app/README.md +++ b/inst/shiny-app/README.md @@ -29,7 +29,7 @@ launchAEFA() 2. **Configure Model**: - Set minimum and maximum number of factors to explore - Choose rotation method (bifactorQ recommended for most cases) - - Select model selection criteria (DIC is default) + - Select model selection criteria (AIC is default for mirt ML/MAP fits) - Optionally enable model history saving to inspect candidate models 3. **Run Analysis**: Click the "Run Analysis" button @@ -56,7 +56,9 @@ following inputs and pass them to `aefa()`: - **Rotation method**: select input mapped to `rotate` (e.g., `bifactorQ`, `geominQ`, `oblimin`). - **Model selection criteria**: select input mapped to - `modelSelectionCriteria` (e.g., `DIC`, `AIC`, `AICc`, `BIC`, `saBIC`). + `modelSelectionCriteria` (e.g., `AIC`, `AICc`, `BIC`, `saBIC`). DIC is an + API-only option and requires a fitted model that actually supplies posterior + DIC; the application never substitutes AIC for it. - **Model history toggle**: checkbox mapped to `saveModelHistory` (recommended for inspecting candidate models). diff --git a/inst/shiny-app/app.R b/inst/shiny-app/app.R index 155725a..8a7f244 100644 --- a/inst/shiny-app/app.R +++ b/inst/shiny-app/app.R @@ -53,8 +53,8 @@ ui <- fluidPage( selected = "bifactorQ"), selectInput("modelSelection", "Model Selection Criteria:", - choices = c("DIC", "AIC", "AICc", "BIC", "SABIC"), - selected = "DIC"), + choices = c("AIC", "AICc", "BIC", "SABIC"), + selected = "AIC"), checkboxInput("saveHistory", "Save Model History", TRUE), @@ -139,7 +139,7 @@ ui <- fluidPage( tags$ul( tags$li(strong("Minimum/Maximum Factors:"), "Set the range of factors to explore. The algorithm will test models with different numbers of factors and select the best one."), tags$li(strong("Rotation Method:"), "Choose the rotation method for factor loadings. 'bifactorQ' is recommended for most cases."), - tags$li(strong("Model Selection Criteria:"), "Choose the criteria for selecting the best model. 'DIC' (Deviance Information Criterion) is the default.") + tags$li(strong("Model Selection Criteria:"), "Choose the criterion for selecting the best model. AIC is the default for the maximum-likelihood/MAP models fitted by mirt; AICc applies the finite-sample correction.") ), h4("Step 3: Run the Analysis"), diff --git a/man/aefa.Rd b/man/aefa.Rd index 74d0e08..7583314 100644 --- a/man/aefa.Rd +++ b/man/aefa.Rd @@ -44,7 +44,7 @@ aefa( resampling = T, samples = 5000, printDebugMsg = F, - modelSelectionCriteria = "DIC", + modelSelectionCriteria = "AIC", saveRawEstModels = F, fitEMatUIRT = F, ranefautocomb = T, @@ -107,7 +107,12 @@ aefa( \item{printDebugMsg}{Do you want to see the debugging messages? default is FALSE} -\item{modelSelectionCriteria}{Which criteria want to use model selection work? 'DIC' (default), 'AIC', 'AICc', 'BIC', 'saBIC' available. AIC and DIC will be identical when no prior parameter distributions are included.} +\item{modelSelectionCriteria}{Information criterion used for model selection. +\code{"AIC"} is the default for the maximum-likelihood/MAP models fitted by +current versions of \pkg{mirt}. \code{"AICc"}, \code{"BIC"}, and +\code{"saBIC"} are also available. \code{"DIC"} is accepted only when the +fitted model actually supplies a posterior-based DIC value; it is never +substituted with AIC.} \item{saveRawEstModels}{Do you want to save raw estimated models before model selection work? default is FALSE} diff --git a/man/kaefa.Rd b/man/kaefa.Rd index 270c091..6bcf5b1 100644 --- a/man/kaefa.Rd +++ b/man/kaefa.Rd @@ -29,10 +29,12 @@ that they are possible improper solutions. The kaefa may inspect this issues from the MMMM or MM in statistical learning theory perspectives using model -selection criteria like the DIC (Kang, 2008; Kang, Cohen, & Sung, 2009; -Jiao, Kamata, Wang, & Jin, 2012; Jiao & Zhang, 2015) with maximising -generalisability of the number of factor decisions in every calibration -(Kang, 2008; Preacher, Zhang, Kim, & Mels, 2013). +likelihood-based model-selection criteria such as AIC, AICc, BIC, and +sample-size-adjusted BIC. Posterior DIC (Spiegelhalter et al., 2002) is +accepted only when the fitted model actually supplies it; DIC is not +reconstructed from AIC. These criteria support generalisable factor-count +decisions in each calibration (Kang, 2008; Preacher, Zhang, Kim, & Mels, +2013). If researcher provide of demographical information in kaefa, kaefa will inspect the optimal number of factor and optimal IRT model, diff --git a/tests/testthat/helper-test-data.R b/tests/testthat/helper-test-data.R index bb57262..f99ab6a 100644 --- a/tests/testthat/helper-test-data.R +++ b/tests/testthat/helper-test-data.R @@ -80,7 +80,16 @@ create_binary_test_data <- function(n_items = 10, n_obs = 100) { } .merge_test_defaults <- function(defaults, dots) { - utils::modifyList(defaults, dots) + dot_names <- names(dots) + if (is.null(dot_names)) { + dot_names <- rep("", length(dots)) + } + + is_named <- !is.na(dot_names) & nzchar(dot_names) + positional <- dots[!is_named] + named <- dots[is_named] + + c(positional, utils::modifyList(defaults, named)) } .skip_expensive_ci_calls <- function(function_name) { diff --git a/tests/testthat/test-aefa-advanced-parameters.R b/tests/testthat/test-aefa-advanced-parameters.R index 3b40ec8..c2ea499 100644 --- a/tests/testthat/test-aefa-advanced-parameters.R +++ b/tests/testthat/test-aefa-advanced-parameters.R @@ -15,16 +15,15 @@ expect_valid_aefa <- function(result) { # Test Suite 1: Model Selection Criteria Variations # ============================================================ -test_that("aefa respects different model selection criteria - DIC", { - test_data <- create_test_data(n_items = 5, n_obs = 100) - - result <- try(aefa(test_data, - minExtraction = 1, - maxExtraction = 1, - modelSelectionCriteria = "DIC"), - silent = TRUE) - - expect_valid_aefa(result) +test_that("DIC is not fabricated for mirt ML/MAP fits", { + expect_error( + kaefa:::.aefaFitCriterionValue( + list(AIC = 120, logLik = -50), + criterion = "DIC", + sample_size = 100 + ), + "posterior deviance draws.*never replaced by AIC" + ) }) test_that("aefa respects different model selection criteria - AIC", { diff --git a/tests/testthat/test-aefa-greedy-algorithm.R b/tests/testthat/test-aefa-greedy-algorithm.R index 676e0fe..0814c07 100644 --- a/tests/testthat/test-aefa-greedy-algorithm.R +++ b/tests/testthat/test-aefa-greedy-algorithm.R @@ -90,7 +90,7 @@ test_that("aefa greedy search evaluates model candidates", { expect_true(is.list(result$itemFitTrials)) last_model <- result$estModelTrials[[length(result$estModelTrials)]] if (isS4(last_model) && "Fit" %in% slotNames(last_model)) { - fit_values <- c(last_model@Fit$AIC, last_model@Fit$BIC, last_model@Fit$DIC) + fit_values <- c(last_model@Fit$AIC, last_model@Fit$BIC, last_model@Fit$SABIC) expect_true(any(is.finite(fit_values)), info = "Selected model should include information criteria values") } @@ -107,7 +107,7 @@ test_that("aefa greedy search evaluates model candidates", { test_that("aefa uses information criteria for model selection", { test_data <- create_test_data(n_items = 6, n_obs = 100) - # The function should select best model based on DIC, AIC, BIC, etc. + # The function should select the best model using a criterion supplied by mirt. result <- try(aefa(test_data, minExtraction = 1, maxExtraction = 2), silent = TRUE) if (!inherits(result, "try-error") && !is.null(result)) { diff --git a/tests/testthat/test-model-selection-criteria.R b/tests/testthat/test-model-selection-criteria.R new file mode 100644 index 0000000..6b387b5 --- /dev/null +++ b/tests/testthat/test-model-selection-criteria.R @@ -0,0 +1,66 @@ +test_that("AICc preserves the Hurvich-Tsai small-sample correction", { + fit <- list(AIC = 120, logLik = -50) + parameter_count <- 10 + expected <- 120 + (2 * parameter_count * (parameter_count + 1)) / + (100 - parameter_count - 1) + + expect_equal( + kaefa:::.aefaFitCriterionValue(fit, "AICc", sample_size = 100), + expected, + tolerance = 1e-12 + ) +}) + +test_that("model-supplied criteria are preserved exactly", { + fit <- list(AIC = 120, AICc = 123.5, BIC = 140, SABIC = 130, DIC = 117.25) + + expect_identical(kaefa:::.aefaFitCriterionValue(fit, "AIC", 100), 120) + expect_identical(kaefa:::.aefaFitCriterionValue(fit, "AICc", 100), 123.5) + expect_identical(kaefa:::.aefaFitCriterionValue(fit, "BIC", 100), 140) + expect_identical(kaefa:::.aefaFitCriterionValue(fit, "saBIC", 100), 130) + expect_identical(kaefa:::.aefaFitCriterionValue(fit, "DIC", 100), 117.25) +}) + +test_that("invalid information-criterion substitutions fail visibly", { + fit <- list(AIC = 120, logLik = -50) + + expect_error( + kaefa:::.aefaFitCriterionValue(fit, "DIC", 100), + "DIC is unavailable" + ) + expect_error( + kaefa:::.aefaFitCriterionValue(fit, "CAIC", 100), + "not an alias for corrected AIC" + ) + expect_error( + kaefa:::.aefaFitCriterionValue(fit, "AICc", 11), + "AICc is undefined" + ) +}) + +test_that("candidate selection retains original indices after skipped fits", { + expect_identical(kaefa:::.aefaBestScoreIndex(c(Inf, 20, 10)), 3L) + expect_identical(kaefa:::.aefaBestScoreIndex(c(NA, Inf)), NA_integer_) +}) + +test_that("item exclusions require a new name present in the fitted model", { + expect_identical( + kaefa:::.aefaValidNewItemExclusions( + excluded = "Item1", + proposed = c("Item1", "Item2", "missing", NA_character_), + available = c("Item1", "Item2", "Item3") + ), + "Item2" + ) + expect_identical( + kaefa:::.aefaValidNewItemExclusions("Item1", "Item1", c("Item1", "Item2")), + character() + ) +}) + +test_that("local worker selection survives unavailable CPU telemetry", { + expect_identical(kaefa:::.aefaParallelProcessorCount(NULL, 12), 4L) + expect_identical(kaefa:::.aefaParallelProcessorCount(c(10, 20), 12), 4L) + expect_identical(kaefa:::.aefaParallelProcessorCount(5, 12), 1L) + expect_identical(kaefa:::.aefaParallelProcessorCount(80, NA_real_), 1L) +}) diff --git a/tests/testthat/test-shiny-app.R b/tests/testthat/test-shiny-app.R index ffc8855..1e5cb9c 100644 --- a/tests/testthat/test-shiny-app.R +++ b/tests/testthat/test-shiny-app.R @@ -191,7 +191,7 @@ test_that("Model selection criteria choices are valid", { ui_html <- as.character(ui) # Check for expected model selection criteria - expected_criteria <- c("DIC", "AIC", "AICc", "BIC", "SABIC") + expected_criteria <- c("AIC", "AICc", "BIC", "SABIC") for (criterion in expected_criteria) { expect_true( diff --git a/tests/testthat/test-test-bootstrap-contract.R b/tests/testthat/test-test-bootstrap-contract.R index 78626f6..14c0e61 100644 --- a/tests/testthat/test-test-bootstrap-contract.R +++ b/tests/testthat/test-test-bootstrap-contract.R @@ -31,3 +31,17 @@ test_that("test file lookup uses explicit base package helpers", { perl = TRUE )) }) + +test_that("test defaults preserve positional data and override named options", { + data_arg <- data.frame(item = 1:3) + + merged <- .merge_test_defaults( + list(NCYCLES = 120, BURNIN = 40), + list(data_arg, BURNIN = 10) + ) + + expect_identical(merged[[1]], data_arg) + expect_identical(names(merged)[[1]], "") + expect_identical(merged$NCYCLES, 120) + expect_identical(merged$BURNIN, 10) +}) From 10bece38849b4b80660fe104eb09e6e116650e04 Mon Sep 17 00:00:00 2001 From: Seongho Bae Date: Mon, 13 Jul 2026 17:08:35 +0900 Subject: [PATCH 2/2] fix(aefa): satisfy current-head review and coverage gates --- R/kaefa.R | 11 ++++++----- inst/WORDLIST | 2 ++ man/aefa.Rd | 3 ++- tests/testthat/test-aefa-advanced-parameters.R | 2 ++ tests/testthat/test-model-selection-criteria.R | 2 ++ 5 files changed, 14 insertions(+), 6 deletions(-) diff --git a/R/kaefa.R b/R/kaefa.R index e0ab631..d0680d9 100644 --- a/R/kaefa.R +++ b/R/kaefa.R @@ -546,7 +546,7 @@ evaluateItemFit <- function(mirtModel, RemoteClusters = NULL, rotate = "bifactor #' Select the original candidate index with the lowest finite score #' #' @param scores Numeric vector aligned one-to-one with model candidates. -#' @return The original candidate index, or code{NA_integer_} when none is valid. +#' @return The original candidate index, or \code{NA_integer_} when none is valid. #' @keywords internal #' @noRd .aefaBestScoreIndex <- function(scores) { @@ -604,11 +604,12 @@ evaluateItemFit <- function(mirtModel, RemoteClusters = NULL, rotate = "bifactor #' #' @param resampling Do you want to do resampling with replace? default is TRUE and activate nrow is over samples argument. #' @param samples Specify the number samples with resampling. default is 5000. -#' @param printDebugMsg Do you want to see the debugging messeages? default is FALSE +#' @param printDebugMsg Show additional debugging messages. Default is FALSE. +#' Nonfatal fallback reasons and terminal errors remain visible regardless. #' @param modelSelectionCriteria Information criterion used for model selection. -#' code{"AIC"} is the default for the maximum-likelihood/MAP models fitted by -#' current versions of pkg{mirt}. code{"AICc"}, code{"BIC"}, and -#' code{"saBIC"} are also available. code{"DIC"} is accepted only when the +#' \code{"AIC"} is the default for the maximum-likelihood/MAP models fitted by +#' current versions of \pkg{mirt}. \code{"AICc"}, \code{"BIC"}, and +#' \code{"saBIC"} are also available. \code{"DIC"} is accepted only when the #' fitted model actually supplies a posterior-based DIC value; it is never #' substituted with AIC. #' @param saveRawEstModels Do you want to save raw estimated models before model selection work? default is FALSE diff --git a/inst/WORDLIST b/inst/WORDLIST index 00d99a3..f03f48c 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -25,6 +25,7 @@ 18637 20 200 +2002 2008 2009 2011 @@ -102,6 +103,7 @@ Kang LAN LCA Leiman +MAP MH MHRM MIRT diff --git a/man/aefa.Rd b/man/aefa.Rd index 7583314..b89cf1e 100644 --- a/man/aefa.Rd +++ b/man/aefa.Rd @@ -105,7 +105,8 @@ aefa( \item{samples}{Specify the number samples with resampling. default is 5000.} -\item{printDebugMsg}{Do you want to see the debugging messages? default is FALSE} +\item{printDebugMsg}{Show additional debugging messages. Default is FALSE. +Nonfatal fallback reasons and terminal errors remain visible regardless.} \item{modelSelectionCriteria}{Information criterion used for model selection. \code{"AIC"} is the default for the maximum-likelihood/MAP models fitted by diff --git a/tests/testthat/test-aefa-advanced-parameters.R b/tests/testthat/test-aefa-advanced-parameters.R index c2ea499..6323489 100644 --- a/tests/testthat/test-aefa-advanced-parameters.R +++ b/tests/testthat/test-aefa-advanced-parameters.R @@ -4,6 +4,8 @@ context("AEFA Advanced Parameter Interactions and Combinations") +.ensure_kaefa_namespace() + expect_valid_aefa <- function(result) { expect_false(inherits(result, "try-error")) expect_s3_class(result, "aefa") diff --git a/tests/testthat/test-model-selection-criteria.R b/tests/testthat/test-model-selection-criteria.R index 6b387b5..723d3d8 100644 --- a/tests/testthat/test-model-selection-criteria.R +++ b/tests/testthat/test-model-selection-criteria.R @@ -1,3 +1,5 @@ +.ensure_kaefa_namespace() + test_that("AICc preserves the Hurvich-Tsai small-sample correction", { fit <- list(AIC = 120, logLik = -50) parameter_count <- 10