diff --git a/semTools/NAMESPACE b/semTools/NAMESPACE index fd39ce8..6fcae1d 100644 --- a/semTools/NAMESPACE +++ b/semTools/NAMESPACE @@ -1,5 +1,8 @@ # Generated by roxygen2: do not edit by hand +S3method(print,epcEquivCheckStd) +S3method(print,summaryEpcEquivFit) +S3method(summary,epcequivfit.data.frame) export(AVE) export(PAVranking) export(SSpower) @@ -17,6 +20,8 @@ export(compareFit) export(discriminantValidity) export(efa.ekc) export(efaUnrotate) +export(epcEquivCheck) +export(epcEquivFit) export(findRMSEApower) export(findRMSEApowernested) export(findRMSEAsamplesize) @@ -47,7 +52,6 @@ export(maximalRelia) export(measEq.syntax) export(measurementInvariance) export(measurementInvarianceCat) -export(miPowerFit) export(modificationIndices.mi) export(modificationindices.mi) export(modindices.mi) @@ -120,6 +124,7 @@ importFrom(lavaan,lavTestLRT) importFrom(lavaan,lavaan) importFrom(lavaan,lavaanList) importFrom(lavaan,lavaanify) +importFrom(lavaan,modificationIndices) importFrom(lavaan,parTable) importFrom(methods,as) importFrom(methods,getMethod) diff --git a/semTools/NEWS.md b/semTools/NEWS.md index 4d50562..361074f 100644 --- a/semTools/NEWS.md +++ b/semTools/NEWS.md @@ -12,7 +12,14 @@ - `return.df=` argument replaced by `simplify=` - Coefficients now returned (by default) in a list, with headers indicating how to appropriately interpret each coefficient. The verbose headers can be suppressed (mimicking old behavior) by setting `simplify = -1L` -- TODO: Sunthud's update +- `miPowerFit()` has been renamed to `epcEquivFit()` and its functionality has been substantially extended. + - The output of `epcEquivFit()` now has its own class, `epcequivfit.data.frame`. A dedicated `summary()` method is provided, which reports both local classifications of fixed parameters and the resulting global equivalence decision. + - By default, `epcEquivFit()` focuses on equivalence testing based on expected parameter changes (EPCs). The Saris, Satorra, and van der Veld (2009) decision rule based on modification indices and power is retained but is no longer treated as the primary result, as simulation results favor the EPC equivalence-testing approach. + - The equivalence-testing procedure has been extended to include an underpowered classification when the EPC confidence interval (CI) is wider than the smallest effect size of interest (SESOI) region. + - The arguments and the resulting column names with `target.epc` has been renamed to reflect the smallest effect size of interest (SESOI), consistent with terminology in the equivalence-testing literature. + - Standardized EPCs for indicator residual covariances are now computed by scaling residual covariances by residual variances, rather than by total indicator variances as in `lavaan` and other SEM packages, providing a more appropriate effect-size representation for residual relationships. + - A new function, `epcEquivCheck()`, is introduced to assess whether trivial misspecification is difficult to establish for a given model under the EPC equivalence-testing framework. + - Clipboard utilities have been updated to accommodate these changes, and all relevant documentation has been revised to reflect the transition from `miPowerFit()` to `epcEquivFit()`. ## Bug Fixes: diff --git a/semTools/R/clipboard.R b/semTools/R/clipboard.R index 6da46fa..2370458 100644 --- a/semTools/R/clipboard.R +++ b/semTools/R/clipboard.R @@ -18,8 +18,8 @@ ##' [FitDiff-class]. ##' @param what The attributes of the `lavaan` object to be copied in the ##' clipboard. `"summary"` is to copy the screen provided from the -##' `summary` function. `"mifit"` is to copy the result from the -##' [miPowerFit()] function. Other attributes listed in the +##' `summary` function. `"epceqfit"` is to copy the result from the +##' [epcEquivFit()] function. Other attributes listed in the ##' `inspect` method in the [lavaan::lavaan-class] could also be ##' used, such as `"coef"`, `"se"`, `"fit"`, `"samp"`, and ##' so on. Ignored for [FitDiff-class]-class objects. @@ -33,7 +33,7 @@ ##' @param writeArgs `list` of additional arguments to be passed to ##' [utils::write.table()] ##' @param \dots Additional arguments when passing a `lavaan` object to the -##' `summary` or [miPowerFit()] function. +##' `summary` or [epcEquivFit()] function. ##' ##' @return The resulting output will be saved into a clipboard or a file. If ##' using the `clipboard` function, users may paste it in the other @@ -60,8 +60,8 @@ ##' # pass additional arguments to summary() method for class?lavaan ##' clipboard(fit, rsquare = TRUE, standardized = TRUE, fit.measures = TRUE) ##' -##' # Copy modification indices and fit stats from the miPowerFit() function -##' clipboard(fit, "mifit") +##' # Copy the EPC equivalence testing results from the epcEquivFit() function +##' clipboard(fit, "epceqfit") ##' ##' # Copy the parameter estimates ##' clipboard(fit, "coef") @@ -78,8 +78,8 @@ ##' # Save the summary of the lavaan object ##' saveFile(fit, "out.txt") ##' -##' # Save modification indices and fit stats from the miPowerFit() function -##' saveFile(fit, "out.txt", "mifit") +##' # Save the EPC equivalence testing results from the epcEquivFit() function +##' saveFile(fit, "out.txt", "epceqfit") ##' ##' # Save the parameter estimates ##' saveFile(fit, "out.txt", "coef") @@ -162,13 +162,13 @@ saveFileLavaan <- function(object, file, what = "summary", tableFormat = FALSE, writeArgs$x <- paste(utils::capture.output(summary(object, ...)), collapse = "\n") } - } else if (what == "mifit") { + } else if (what %in% c("epceqfit", "mifit")) { # "mifit" retained for backward compatibility if (tableFormat) { - writeArgs$x <- miPowerFit(object, ...) + writeArgs$x <- epcEquivFit(object, ...) if (is.null(writeArgs$row.names)) writeArgs$row.names <- FALSE if (is.null(writeArgs$col.names)) writeArgs$col.names <- TRUE } else { - writeArgs$x <- paste(utils::capture.output(miPowerFit(object, ...)), + writeArgs$x <- paste(utils::capture.output(epcEquivFit(object, ...)), collapse = "\n") } } else { diff --git a/semTools/R/fitIndices.R b/semTools/R/fitIndices.R index 2df5f15..50740f9 100644 --- a/semTools/R/fitIndices.R +++ b/semTools/R/fitIndices.R @@ -152,8 +152,8 @@ ##' ##' @seealso ##' \itemize{ -##' \item [miPowerFit()] For the modification indices and their -##' power approach for model fit evaluation +##' \item [epcEquivFit()] For the equivalence testing based on expected +##' parameter changes for model fit evaluation ##' \item [nullRMSEA()] For RMSEA of the default independence model ##' } ##' diff --git a/semTools/R/miPowerFit.R b/semTools/R/miPowerFit.R index 8833f2d..b1515d5 100644 --- a/semTools/R/miPowerFit.R +++ b/semTools/R/miPowerFit.R @@ -1,390 +1,1538 @@ -### Sunthud Pornprasertmanit and Terrence D. Jorgensen (added ... argument) -### Last updated: 2 September 2021 - - -##' Modification indices and their power approach for model fit evaluation -##' -##' The model fit evaluation approach using modification indices and expected -##' parameter changes. -##' -##' To decide whether a parameter should be freed, one can inspect its -##' modification index (MI) and expected parameter change (EPC). -##' Those values can be used to evaluate model fit by 2 methods. -##' -##' Method 1: Saris, Satorra, and van der Veld (2009, pp. 570--573) used -##' power (probability of detecting a significant MI) and EPC to decide whether -##' to free a parametr. First, one should evaluate whether a parameter's MI -##' is significant. Second, one should evaluate whether the power to detect a -##' target EPC is high enough. The combination of criteria leads to the -##' so-called "JRule" first implemented with LISREL (van der Veld et al., 2008): -##' -##' \itemize{ -##' \item If the MI is not significant and the power is low, -##' the test is inconclusive. -##' \item If the MI is not significant and the power is high, -##' there is no misspecification. -##' \item If the MI is significant and the power is low, -##' the fixed parameter is misspecified. -##' \item If the MI is significant and the power is high, -##' the EPC is investigated. If the EPC is large (greater than the -##' the target EPC), the parameter is misspecified. If the EPC is low -##' (lower than the target EPC), the parameter is not misspecificied. -##' } -##' -##' Method 2: The confidence interval (CI) of an EPC is calculated. -##' These CIs are compared with the range of trivial -##' misspecification, which could be (-`delta`, `delta`) or (0, -##' `delta`) for nonnegative parameters. -##' -##' \itemize{ -##' \item If a CI overlaps with the range of trivial misspecification, -##' the test is inconclusive. -##' \item If a CI completely exceeds the range of trivial misspecification, -##' the fixed parameters are severely misspecified. -##' \item If a CI is completely within the range of trivial misspecification, -##' the fixed parameters are trivially misspecified. -##' } -##' -##' -##' @aliases miPowerFit miPowerFit -##' @importFrom lavaan lavInspect -##' @importFrom stats qnorm qchisq pchisq -##' -##' @param lavaanObj The lavaan model object used to evaluate model fit -##' @param stdLoad The amount of standardized factor loading that one would like -##' to be detected (rejected). The default value is 0.4, which is suggested by -##' Saris and colleagues (2009, p. 571). -##' @param cor The amount of factor or error correlations that one would like to -##' be detected (rejected). The default value is 0.1, which is suggested by -##' Saris and colleagues (2009, p. 571). -##' @param stdBeta The amount of standardized regression coefficients that one -##' would like to be detected (rejected). The default value is 0.1, which is -##' suggested by Saris and colleagues (2009, p. 571). -##' @param intcept The amount of standardized intercept (similar to Cohen's -##' *d* that one would like to be detected (rejected). The default value -##' is 0.2, which is equivalent to a low effect size proposed by Cohen (1988, -##' 1992). -##' @param stdDelta The vector of the standardized parameters that one would -##' like to be detected (rejected). If this argument is specified, the value -##' here will overwrite the other arguments above. The order of the vector -##' must be the same as the row order from modification indices from the -##' `lavaan` object. If a single value is specified, the value will be -##' applied to all parameters. -##' @param delta The vector of the unstandardized parameters that one would like -##' to be detected (rejected). If this argument is specified, the value here -##' will overwrite the other arguments above. The order of the vector must be -##' the same as the row order from modification indices from the `lavaan` -##' object. If a single value is specified, the value will be applied to all -##' parameters. -##' @param cilevel The confidence level of the confidence interval of expected -##' parameter changes. The confidence intervals are used in the equivalence -##' testing. -##' @param \dots arguments passed to [lavaan::modificationIndices()], -##' except for `delta`, which is already an argument (which can be -##' substituted for `stdDelta` or specific sets of parameters using -##' `stdLoad`, `cor`, `stdBeta`, and `intcept`). -##' -##' @return A data frame with these variables: -##' \enumerate{ -##' \item `lhs`: The left-hand side variable, with respect to the operator in -##' in the lavaan [lavaan::model.syntax()] -##' \item `op`: The lavaan syntax operator: "~~" represents covariance, -##' "=~" represents factor loading, "~" represents regression, and -##' "~1" represents intercept. -##' \item `rhs`: The right-hand side variable -##' \item `group`: The level of the group variable for the parameter in question -##' \item `mi`: The modification index of the fixed parameter -##' \item `epc`: The EPC if the parameter is freely estimated -##' \item `target.epc`: The target EPC that represents the minimum size -##' of misspecification that one would like to be detected -##' by the test with a high power -##' \item `std.epc`: The standardized EPC if the parameter is freely estimated -##' \item `std.target.epc`: The standardized target expected parameter change -##' \item `significant.mi`: Represents whether the modification index value is -##' significant -##' \item `high.power`: Represents whether the power is enough to detect the -##' target expected parameter change -##' \item `decision.pow`: The decision whether the parameter is misspecified -##' or not based on Saris et al's method: `"M"` represents the parameter -##' is misspecified, `"NM"` represents the parameter is not misspecified, -##' `"EPC:M"` represents the parameter is misspecified decided by -##' checking the expected parameter change value, `"EPC:NM"` represents -##' the parameter is not misspecified decided by checking the expected -##' parameter change value, and `"I"` represents the decision is -##' inconclusive. -##' \item `se.epc`: The standard errors of the expected parameter changes. -##' \item `lower.epc`: The lower bound of the confidence interval of expected -##' parameter changes. -##' \item `upper.epc`: The upper bound of the confidence interval of expected -##' parameter changes. -##' \item `lower.std.epc`: Lower confidence limit of standardized EPCs -##' \item `upper.std.epc`: Upper confidence limit of standardized EPCs -##' \item `decision.ci`: Decision whether the parameter is misspecified -##' based on the CI method: `"M"` represents the -##' parameter is misspecified, `"NM"` represents the parameter is not -##' misspecified, and `"I"` represents the decision is inconclusive. -##' } -##' -##' The row numbers matches with the results obtained from the -##' `inspect(object, "mi")` function. -##' -##' @author Sunthud Pornprasertmanit (\email{psunthud@@gmail.com}) -##' -##' @seealso [moreFitIndices()] For the additional fit indices -##' information -##' -##' @references -##' Cohen, J. (1988). *Statistical power analysis for the -##' behavioral sciences* (2nd ed.). Hillsdale, NJ: Erlbaum. -##' -##' Cohen, J. (1992). A power primer. *Psychological Bulletin, 112*(1), -##' 155--159. \doi{10.1037/0033-2909.112.1.155} -##' -##' Saris, W. E., Satorra, A., & van der Veld, W. M. (2009). Testing structural -##' equation models or detection of misspecifications? *Structural Equation -##' Modeling, 16*(4), 561--582. \doi{10.1080/10705510903203433} -##' -##' van der Veld, W. M., Saris, W. E., & Satorra, A. (2008). -##' *JRule 3.0 Users Guide*. \doi{10.13140/RG.2.2.13609.90729} -##' -##' @examples -##' -##' library(lavaan) -##' -##' HS.model <- ' visual =~ x1 + x2 + x3 ' -##' fit <- cfa(HS.model, data = HolzingerSwineford1939, -##' group = "sex", group.equal = c("loadings","intercepts")) -##' miPowerFit(fit, free.remove = FALSE, op = "=~") # loadings -##' miPowerFit(fit, free.remove = FALSE, op = "~1") # intercepts -##' -##' model <- ' -##' # latent variable definitions -##' ind60 =~ x1 + x2 + x3 -##' dem60 =~ y1 + a*y2 + b*y3 + c*y4 -##' dem65 =~ y5 + a*y6 + b*y7 + c*y8 -##' -##' # regressions -##' dem60 ~ ind60 -##' dem65 ~ ind60 + dem60 -##' -##' # residual correlations -##' y1 ~~ y5 -##' y2 ~~ y4 + y6 -##' y3 ~~ y7 -##' y4 ~~ y8 -##' y6 ~~ y8 -##' ' -##' fit2 <- sem(model, data = PoliticalDemocracy, meanstructure = TRUE) -##' miPowerFit(fit2, stdLoad = 0.3, cor = 0.2, stdBeta = 0.2, intcept = 0.5) -##' -##' @export -miPowerFit <- function(lavaanObj, stdLoad = 0.4, cor = 0.1, stdBeta = 0.1, - intcept = 0.2, stdDelta = NULL, delta = NULL, +### Sunthud Pornprasertmanit; with contributions by Terrence D. Jorgensen +### Last updated: 3 February 2026 + + +#' EPC Equivalence Fit Evaluation Using Modification Indices +#' +#' Evaluates model fit from an equivalence-testing perspective by +#' aggregating local EPC-based diagnostics into a global, fit-style +#' assessment. The procedure combines modification indices (MI), +#' expected parameter changes (EPC), statistical power, and confidence +#' intervals relative to a smallest effect size of interest (SESOI). +#' +#' Two complementary local decision rules are implemented: +#' +#' \strong{Method 1 (Power-based; Saris, Satorra, & van der Veld, 2009).} +#' Modification indices, statistical power, and EPC magnitude are jointly +#' evaluated (the J-rule) to classify fixed parameters as misspecified, +#' not misspecified, or inconclusive. +#' +#' \strong{Method 2 (CI-based equivalence testing).} +#' Confidence intervals of EPCs are compared against a trivial +#' misspecification region defined by the SESOI to determine whether +#' fixed parameters are severely misspecified, trivially misspecified, +#' underpowered, or inconclusive. +#' +#' The resulting local classifications are returned in a single data +#' frame and can be summarized to yield a global equivalence-style fit +#' evaluation. +#' +#' @param lavaanObj A fitted \code{lavaan} object used to evaluate model fit. +#' @param stdLoad Standardized factor loading defining the SESOI for +#' loading misspecifications. Default is 0.4. +#' @param cor Default standardized correlation defining the SESOI for +#' covariance misspecifications. Used for both latent and residual +#' covariances unless overridden. +#' @param corLatent Standardized latent factor correlation defining the +#' SESOI for latent covariance misspecifications. If \code{NULL}, +#' defaults to \code{cor}. +#' @param corResidual Standardized residual correlation defining the +#' SESOI for residual covariance misspecifications. If \code{NULL}, +#' defaults to \code{cor}. +#' @param stdBeta Standardized regression coefficient defining the SESOI +#' for structural misspecifications. Default is 0.1. +#' @param stdIntcept Standardized intercept (Cohen's \emph{d}) defining +#' the SESOI for intercept misspecifications. Default is 0.2. +#' @param stdSesoi Optional vector of standardized SESOI values. If +#' provided, overrides operator-specific SESOI definitions. +#' @param sesoi Optional vector of unstandardized SESOI values. If +#' provided, overrides \code{stdSesoi} and all operator-specific SESOI +#' arguments. +#' @param cilevel Confidence level for EPC confidence intervals used in +#' CI-based equivalence testing. +#' @param \dots Additional arguments passed to +#' \code{\link[lavaan]{modificationIndices}}. +#' +#' @details +#' This function provides a local-to-global equivalence-based alternative +#' to traditional exact-fit evaluation. It is designed to assess whether +#' fixed parameters are substantively misspecified relative to a SESOI, +#' rather than whether a model fits exactly. +#' +#' Models with categorical indicators or unsupported constraints may +#' not be fully supported. +#' +#' @return A data frame with one row per fixed parameter, containing: +#' \enumerate{ +#' \item Parameter identifiers: \code{lhs}, \code{op}, \code{rhs}, and \code{group}. +#' \item Modification index (\code{mi}) and expected parameter change estimates (\code{epc}). +#' \item Unstandardized and standardized smallest effect size of interest values +#' (\code{sesoi}, \code{std.sesoi}). +#' \item Power-based decision (\code{decision.pow}) and related diagnostics, including +#' whether the modification index is statistically significant +#' (\code{significant.mi}) and whether the misfit at the SESOI has power greater than 0.80 +#' (\code{high.power}). +#' \item EPC-related statistics, including the standard error of the EPC +#' (\code{se.epc}), confidence interval bounds for the EPC +#' (\code{lower.epc}, \code{upper.epc}), and confidence interval bounds for the +#' standardized EPC (\code{lower.std.epc}, \code{upper.std.epc}). +#' \item Confidence-interval–based equivalence decision (\code{decision.ci}). +#' } +#' +#' @references +#' Saris, W. E., Satorra, A., & van der Veld, W. M. (2009). +#' Testing structural equation models or detection of misspecifications? +#' \emph{Structural Equation Modeling, 16}(4), 561--582. +#' +#' @seealso \code{\link{epcEquivCheck}} +#' +#' @importFrom stats qchisq pchisq qnorm +#' @importFrom lavaan modificationIndices lavNames fitMeasures +#' +#' @aliases epcEquivFit miPowerFit +#' +#' @examples +#' +#' library(lavaan) +#' +#' one.model <- ' onefactor =~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 ' +#' fit <- cfa(one.model, data = HolzingerSwineford1939) +#' out <- epcEquivFit(fit) +#' out +#' summary(out) +#' +#' @export +epcEquivFit <- function(lavaanObj, + stdLoad = 0.4, + cor = 0.1, + corLatent = NULL, + corResidual = NULL, + stdBeta = 0.1, + stdIntcept = 0.2, + stdSesoi = NULL, + sesoi = NULL, cilevel = 0.90, ...) { - mi <- lavaan::modificationIndices(lavaanObj, ...) - mi <- mi[mi$op != "==",] - sigma <- mi[,"epc"] / sqrt(mi[,"mi"]) - if (is.null(delta)) { - if (is.null(stdDelta)) - stdDelta <- getTrivialEpc(mi, stdLoad = stdLoad, cor = cor, - stdBeta = stdBeta, intcept = intcept) - if (length(stdDelta) == 1) stdDelta <- rep(stdDelta, nrow(mi)) - delta <- unstandardizeEpc(mi, stdDelta, findTotalVar(lavaanObj)) - } - if (length(delta) == 1) delta <- rep(delta, nrow(mi)) - ncp <- (delta / sigma)^2 - alpha <- 0.05 - desiredPow <- 0.80 - cutoff <- qchisq(1 - alpha, df = 1) - pow <- 1 - pchisq(cutoff, df = 1, ncp = ncp) - sigMI <- mi[,"mi"] > cutoff - highPow <- pow > desiredPow - group <- rep(1, nrow(mi)) - if ("group" %in% colnames(mi)) group <- mi[ , "group"] - decision <- mapply(decisionMIPow, sigMI = sigMI, highPow = highPow, - epc = mi[ , "epc"], trivialEpc = delta) - if (is.null(stdDelta)) stdDelta <- standardizeEpc(mi, findTotalVar(lavaanObj), - delta = delta) - result <- cbind(mi[ , 1:3], group, as.numeric(mi[ , "mi"]), mi[ , "epc"], - delta, standardizeEpc(mi, findTotalVar(lavaanObj)), - stdDelta, sigMI, highPow, decision) - # New method - crit <- abs(qnorm((1 - cilevel)/2)) - seepc <- abs(result[,6]) / sqrt(abs(result[,5])) - lowerepc <- result[,6] - crit * seepc - upperepc <- result[,6] + crit * seepc - stdlowerepc <- standardizeEpc(mi, findTotalVar(lavaanObj), delta = lowerepc) - stdupperepc <- standardizeEpc(mi, findTotalVar(lavaanObj), delta = upperepc) - isVar <- mi[,"op"] == "~~" & mi[,"lhs"] == mi[,"rhs"] - decisionci <- mapply(decisionCIEpc, targetval = as.numeric(stdDelta), - lower = stdlowerepc, upper = stdupperepc, - positiveonly = isVar) - result <- cbind(result, seepc, lowerepc, upperepc, stdlowerepc, - stdupperepc, decisionci) - result <- result[!is.na(decision), ] - colnames(result) <- c("lhs","op","rhs","group","mi","epc","target.epc", - "std.epc","std.target.epc","significant.mi", - "high.power","decision.pow","se.epc","lower.epc", - "upper.epc","lower.std.epc","upper.std.epc","decision.ci") - class(result) <- c("lavaan.data.frame","data.frame") - return(result) + dots <- list(...) + + df_model <- lavaan::fitMeasures(lavaanObj, "df") + + if (is.na(df_model) || df_model <= 0) { + stop( + "epcEquivFit() requires a model with positive degrees of freedom (df > 0).\n", + "The supplied lavaan model has df = ", df_model, ".\n", + "EPC-based equivalence testing is not defined for just-identified or ", + "under-identified models.", + call. = FALSE + ) + } + + # deprecated aliases + if (!is.null(dots$intcept)) { + warning( + "'intcept' is deprecated; please use 'stdIntcept' instead.", + call. = FALSE + ) + stdIntcept <- dots$intcept + } + + if (!is.null(dots$stdDelta)) { + warning( + "'stdDelta' is deprecated; please use 'stdSesoi' instead.", + call. = FALSE + ) + stdSesoi <- dots$stdDelta + } + + if (!is.null(dots$delta)) { + warning( + "'delta' is deprecated; please use 'sesoi' instead.", + call. = FALSE + ) + sesoi <- dots$delta + } + + lv_names <- lavaan::lavNames(lavaanObj, type="lv") + mi <- lavaan::modificationIndices(lavaanObj) + mi <- mi[mi$op != "==",] + sigma <- mi[,"epc"] / sqrt(mi[,"mi"]) + sigma[!is.finite(sigma)] <- NA_real_ + if (is.null(corLatent)) corLatent <- cor + if (is.null(corResidual)) corResidual <- cor + if (is.null(sesoi)) { + if (is.null(stdSesoi)) + stdSesoi <- getTrivialEpc(mi, lv_names=lv_names, stdLoad = stdLoad, + corLatent = corLatent, corResidual = corResidual, + stdBeta = stdBeta, stdIntcept = stdIntcept) + if (length(stdSesoi) == 1) stdSesoi <- rep(stdSesoi, nrow(mi)) + sesoi <- unstandardizeEpc(mi, stdSesoi, lavInspectTotalVar(lavaanObj), lavInspectResidualVar(lavaanObj)) + } + if (length(sesoi) == 1) sesoi <- rep(sesoi, nrow(mi)) + ncp <- (sesoi / sigma)^2 + alpha <- 0.05 + desiredPow <- 0.80 + cutoff <- stats::qchisq(1 - alpha, df = 1) + pow <- 1 - stats::pchisq(cutoff, df = 1, ncp = ncp) + sigMI <- mi[,"mi"] > cutoff + highPow <- pow > desiredPow + group <- rep(1, nrow(mi)) + if ("group" %in% colnames(mi)) group <- mi[ , "group"] + decision <- mapply(decisionMIPow, sigMI = sigMI, highPow = highPow, + epc = mi[ , "epc"], trivialEpc = sesoi) + if (is.null(stdSesoi)) stdSesoi <- standardizeEpc(mi, lavInspectTotalVar(lavaanObj), + lavInspectResidualVar(lavaanObj), + sesoi = sesoi) + result <- cbind(mi[ , 1:3], group, as.numeric(mi[ , "mi"]), mi[ , "epc"], + sesoi, mi[ , "sepc.all"], + stdSesoi, sigMI, highPow, decision) + # New method + crit <- abs(stats::qnorm((1 - cilevel)/2)) + seepc <- abs(result[,6]) / sqrt(abs(result[,5])) + lowerepc <- result[,6] - crit * seepc + upperepc <- result[,6] + crit * seepc + stdlowerepc <- standardizeEpc(mi, lavInspectTotalVar(lavaanObj), + lavInspectResidualVar(lavaanObj), sesoi = lowerepc) + stdupperepc <- standardizeEpc(mi, lavInspectTotalVar(lavaanObj), + lavInspectResidualVar(lavaanObj), sesoi = upperepc) + isVar <- mi[,"op"] == "~~" & mi[,"lhs"] == mi[,"rhs"] + decisionci <- mapply(decisionCIEpc, targetval = as.numeric(stdSesoi), + lower = stdlowerepc, upper = stdupperepc, + positiveonly = isVar) + result <- cbind(result, pow, seepc, lowerepc, upperepc, stdlowerepc, + stdupperepc, decisionci) + result <- result[!is.na(decision), ] + colnames(result) <- c("lhs","op","rhs","group","mi","epc","sesoi", + "std.epc","std.sesoi","significant.mi", + "high.power","decision.pow","pow","se.epc","lower.epc", + "upper.epc","lower.std.epc","upper.std.epc","decision.ci") + class(result) <- c("epcequivfit.data.frame","lavaan.data.frame","data.frame") + + # backward compatibility alias + result$target.epc <- result$sesoi + result$std.target.epc <- result$std.sesoi + + return(result) } +#FIXME: Change to .Defunct after a few version updates +miPowerFit <- function(...) { + .Deprecated("epcEquivFit", + msg = "miPowerFit() has been replaced by epcEquivFit().") + epcEquivFit(...) +} + +#' EPC Equivalence Feasibility Check for Standardized Parameters +#' +#' Performs an EPC-based feasibility check to assess whether a set of +#' standardized population parameters defines a valid population +#' covariance matrix and whether trivially misspecified parameters +#' remain within a user-defined smallest effect size of interest (SESOI). +#' Feasibility is evaluated by constructing implied population models +#' under targeted parameter perturbations and examining EPC behavior +#' using \code{\link{epcEquivFit}}. +#' +#' This function focuses on standardized parameters and supports +#' recursive SEMs with continuous indicators only. +#' +#' @param lavaanObj A fitted \code{lavaan} object representing the target model. +#' @param minRelEffect A scalar in (0, 1) specifying the minimum relative +#' magnitude of the standardized perturbation to be evaluated. The +#' default value of 0.75 indicates that perturbations equal to 75\% of +#' the SESOI are treated as trivial. If EPCs exceed the SESOI under +#' such perturbations, EPC equivalence testing is not recommended. +#' @param stdLoad Standardized factor loading used to define the SESOI +#' for loading misspecifications. +#' @param cor Standardized correlation used as a default SESOI for +#' covariance misspecifications. This value is used for both latent +#' and residual covariances unless overridden by +#' \code{corLatent} or \code{corResidual}. +#' @param corLatent Standardized latent factor correlation used to +#' define the SESOI for latent covariance misspecifications. If +#' \code{NULL}, defaults to \code{cor}. +#' @param corResidual Standardized residual correlation used to define +#' the SESOI for indicator residual covariance misspecifications. If +#' \code{NULL}, defaults to \code{cor}. +#' @param stdBeta Standardized regression coefficient used to define +#' the SESOI for structural misspecifications. +#' +#' @details +#' The procedure first checks whether the standardized parameters imply +#' a positive definite population covariance matrix. It then evaluates +#' EPC behavior under both positive and negative trivial +#' misspecifications by repeatedly constructing implied population +#' covariance matrices with perturbed parameters +#' (\code{minRelEffect} \eqn{\times} SESOI), refitting the model, and +#' re-evaluating EPCs. +#' +#' Models with categorical indicators, formative indicators, or +#' multiple-group structures are not supported. +#' +#' @return An object of class \code{"epcEquivCheckStd"} containing: +#' \itemize{ +#' \item \code{feasible}: Logical indicator of whether a valid +#' standardized population model exists. +#' \item \code{any_M}: Logical indicator of whether any EPC exceeded +#' the SESOI under the evaluated misspecifications. +#' \item \code{recommendation}: Character string summarizing feasibility +#' (e.g., \code{"RECOMMENDED"}, \code{"NOT RECOMMENDED"}). +#' \item \code{M_table}: Data frame summarizing EPCs exceeding the SESOI, +#' if any. +#' \item \code{testeffect}: Data frame reporting the smallest tested +#' standardized perturbations in each direction. +#' } +#' +#' @importFrom lavaan lavaan +#' +#' @seealso \code{\link{epcEquivFit}} +#' +#' @examples +#' +#' library(lavaan) +#' +#' one.model <- ' onefactor =~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 ' +#' fit <- cfa(one.model, data = HolzingerSwineford1939) +#' \donttest{ +#' epcEquivCheck(fit) +#' } +#' +#' @export +epcEquivCheck <- function(lavaanObj, + minRelEffect = 0.75, + stdLoad = 0.4, + cor = 0.1, + corLatent = NULL, + corResidual = NULL, + stdBeta = 0.1) { + + tol <- 1e-10 + + .new_epcEquivCheckStd_infeasible <- function(reason) { + out <- list( + feasible = FALSE, + any_M = NA, + recommendation = "NOT APPLICABLE", + reason = reason, + M_table = NULL, + testeffect = NULL + ) + class(out) <- "epcEquivCheckStd" + out + } + + + if (minRelEffect <= 0 || minRelEffect >= 1) { + stop("minRelEffect must be between 0 and 1.") + } + + if (is.null(corLatent)) corLatent <- cor + if (is.null(corResidual)) corResidual <- cor + + # Candidate misspecifications and EPC-based feasibility outputs + miout <- epcEquivFit(lavaanObj, stdLoad = stdLoad, corLatent = corLatent, + corResidual = corResidual, stdBeta = stdBeta) + + # Scope checks: features not supported in this standardized-parameter feasibility check + if (any(miout$op == "|")) { + stop("Models with categorical indicators are not supported.") + } + if (any(miout$op == "~1")) { + stop("Models with a mean structure are not supported.") + } + if (any(miout$op == "<~")) { + stop("Models with formative indicators are not supported.") + } + if (lavaan::lavInspect(lavaanObj, "ngroups") > 1) { + stop("Multiple-group models are not supported yet.") + } + + # Extract standardized parameter matrices + std <- lavaan::lavInspect(lavaanObj, "std.all") + lambda <- std$lambda + corpsi <- cov2cor_safe(std$psi) + cortheta <- cov2cor_safe(std$theta) + stdbeta <- std$beta + if (is.null(stdbeta)) stdbeta <- matrix(0, nrow(corpsi), ncol(corpsi), + dimnames = dimnames(corpsi)) + + # ---- Existence check: do standardized parameters define ANY PD population Sigma? ---- + residVarPsi0 <- findFactorResidualVar( + beta = stdbeta, + corPsi = corpsi, + totalVarPsi = rep(1, nrow(stdbeta)) + ) + if (any(!is.finite(residVarPsi0)) || any(residVarPsi0 < -tol)) { + return(.new_epcEquivCheckStd_infeasible("std_params_not_generative")) + } + residVarPsi0[abs(residVarPsi0) < tol] <- 0 + + Phi0 <- findFactorTotalCov( + beta = stdbeta, + corPsi = corpsi, + errorVarPsi = residVarPsi0 + ) + + Sigma_y0 <- lambda %*% Phi0 %*% t(lambda) + + residVarTheta0 <- findIndResidualVar( + lambda = lambda, + totalFactorCov = Phi0, + totalVarTheta = rep(1, nrow(lambda)) + ) + if (any(!is.finite(residVarTheta0)) || any(residVarTheta0 < -tol)) { + return(.new_epcEquivCheckStd_infeasible("std_params_not_generative")) + } + residVarTheta0[abs(residVarTheta0) < tol] <- 0 + + Theta0 <- cor2cov_safe(cortheta, sqrt(residVarTheta0)) + Sigma0 <- Sigma_y0 + Theta0 + + if (!isPD(Sigma0)) { + return(.new_epcEquivCheckStd_infeasible("std_params_not_generative")) + } + + # ---- Search feasibility for positive perturbations ---- + result <- matrix(NA, nrow(miout), nrow(miout)) + sepc_mat <- matrix(NA, nrow(miout), nrow(miout)) + testeffect <- rep(NA_real_, nrow(miout)) + kseq <- seq(minRelEffect, 0.1, length.out = 10) # decreasing misspecification + + for (i in 1:nrow(miout)) { + row <- miout[i,] + found <- FALSE + k_found <- NA_real_ + + for (k in kseq) { + + # Start from base standardized parameter matrices each time + tlambda <- lambda + tcorpsi <- corpsi + tcortheta <- cortheta + tstdbeta <- stdbeta + tempRelEffect <- NA_real_ + + # Apply a single targeted perturbation (scaled by k) + if (row$op == "=~") { + tlambda[row$rhs, row$lhs] <- stdLoad * k + tempRelEffect <- stdLoad + } else if (row$op == "~~") { + + # latent (psi) vs observed residual (theta) determined by membership + if (row$lhs %in% rownames(corpsi)) { + tcorpsi[row$lhs, row$rhs] <- corLatent * k + tcorpsi[row$rhs, row$lhs] <- corLatent * k + tempRelEffect <- corLatent + } else if (row$lhs %in% rownames(lambda)) { + tcortheta[row$lhs, row$rhs] <- corResidual * k + tcortheta[row$rhs, row$lhs] <- corResidual * k + tempRelEffect <- corResidual + } + + } else if (row$op == "~") { + tstdbeta[row$lhs, row$rhs] <- stdBeta * k + tempRelEffect <- stdBeta + } + + # Derive implied factor residual variances needed for total factor covariance + residVarPsi <- findFactorResidualVar( + beta = tstdbeta, + corPsi = tcorpsi, + totalVarPsi = rep(1, nrow(tstdbeta)) + ) + if (any(!is.finite(residVarPsi)) || any(residVarPsi < -tol)) next + residVarPsi[abs(residVarPsi) < tol] <- 0 + + Phi <- tryCatch( + findFactorTotalCov( + beta = tstdbeta, + corPsi = tcorpsi, + errorVarPsi = residVarPsi + ), + error = function(e) NULL + ) + if (is.null(Phi) || any(!is.finite(Phi))) next + + # Implied indicator covariance (excluding residuals) + Sigma_y_noTheta <- tlambda %*% Phi %*% t(tlambda) + + # Derive indicator residual variances under standardized total variances (=1) + residVarTheta <- findIndResidualVar( + lambda = tlambda, + totalFactorCov = Phi, + totalVarTheta = rep(1, nrow(tlambda)) + ) + if (any(!is.finite(residVarTheta)) || any(residVarTheta < -tol)) next + residVarTheta[abs(residVarTheta) < tol] <- 0 + + Theta <- cor2cov_safe(tcortheta, sqrt(residVarTheta)) + Sigma_implied <- Sigma_y_noTheta + Theta + + # Only proceed if Sigma is PD and the model is estimable under that population + if (isPD(Sigma_implied)) { + tempout <- tryCatch( + suppressWarnings(lavaan::lavaan(lavaanObj, + sample.cov = Sigma_implied, + sample.nobs = 1000000L, + std.lv = TRUE)), + error = function(e) NULL + ) + if (is.null(tempout)) next + + if (lavCheckAdmissibleFit(tempout)) { + tempEqTest <- epcEquivFit(tempout, + stdLoad = stdLoad, + corLatent = corLatent, + corResidual = corResidual, + stdBeta = stdBeta) + result[i, ] <- tempEqTest[, "decision.ci"] + sepc_mat[i, ] <- tempEqTest[,"std.epc"] + found <- TRUE + k_found <- k + } + } + + if (found) break + } + testeffect[i] <- k_found * tempRelEffect + } + + # ---- Search feasibility for negative perturbations ---- + result2 <- matrix(NA, nrow(miout), nrow(miout)) + sepc_mat2 <- matrix(NA, nrow(miout), nrow(miout)) + testeffect2 <- rep(NA_real_, nrow(miout)) + kseq2 <- seq(minRelEffect, 0.1, length.out = 10) # decreasing misspecification + + for (i in 1:nrow(miout)) { + row <- miout[i,] + found2 <- FALSE + k_found2 <- NA_real_ + + for (k2 in kseq2) { + + tlambda2 <- lambda + tcorpsi2 <- corpsi + tcortheta2 <- cortheta + tstdbeta2 <- stdbeta + tempRelEffect2 <- NA_real_ + if (row$op == "=~") { + tlambda2[row$rhs, row$lhs] <- -stdLoad * k2 + tempRelEffect2 <- -stdLoad + } else if (row$op == "~~") { + + if (row$lhs %in% rownames(corpsi)) { + tcorpsi2[row$lhs, row$rhs] <- -corLatent * k2 + tcorpsi2[row$rhs, row$lhs] <- -corLatent * k2 + tempRelEffect2 <- -corLatent + } else if (row$lhs %in% rownames(lambda)) { + tcortheta2[row$lhs, row$rhs] <- -corResidual * k2 + tcortheta2[row$rhs, row$lhs] <- -corResidual * k2 + tempRelEffect2 <- -corResidual + } + + } else if (row$op == "~") { + tstdbeta2[row$lhs, row$rhs] <- -stdBeta * k2 + tempRelEffect2 <- -stdBeta + } + + residVarPsi2 <- findFactorResidualVar( + beta = tstdbeta2, + corPsi = tcorpsi2, + totalVarPsi = rep(1, nrow(tstdbeta2)) + ) + if (any(!is.finite(residVarPsi2)) || any(residVarPsi2 < -tol)) next + residVarPsi2[abs(residVarPsi2) < tol] <- 0 + + Phi2 <- tryCatch( + findFactorTotalCov( + beta = tstdbeta2, + corPsi = tcorpsi2, + errorVarPsi = residVarPsi2 + ), + error = function(e) NULL + ) + if (is.null(Phi2) || any(!is.finite(Phi2))) next + + Sigma_y_noTheta2 <- tlambda2 %*% Phi2 %*% t(tlambda2) + + residVarTheta2 <- findIndResidualVar( + lambda = tlambda2, + totalFactorCov = Phi2, + totalVarTheta = rep(1, nrow(tlambda2)) + ) + if (any(!is.finite(residVarTheta2)) || any(residVarTheta2 < -tol)) next + residVarTheta2[abs(residVarTheta2) < tol] <- 0 + + Theta2 <- cor2cov_safe(tcortheta2, sqrt(residVarTheta2)) + Sigma_implied2 <- Sigma_y_noTheta2 + Theta2 + + if (isPD(Sigma_implied2)) { + tempout2 <- tryCatch( + suppressWarnings(lavaan(lavaanObj, + sample.cov = Sigma_implied2, + sample.nobs = 1000000L, + std.lv = TRUE)), + error = function(e) NULL + ) + if (is.null(tempout2)) next + + if (lavCheckAdmissibleFit(tempout2)) { + tempEqTest2 <- epcEquivFit(tempout2, + stdLoad = stdLoad, + corLatent = corLatent, + corResidual = corResidual, + stdBeta = stdBeta) + result2[i, ] <- tempEqTest2[, "decision.ci"] + sepc_mat2[i, ] <- tempEqTest2[,"std.epc"] + found2 <- TRUE + k_found2 <- k2 + } + } + + if (found2) break + } + + testeffect2[i] <- tempRelEffect2 * k_found2 + } + + resultall <- cbind(result, result2) + M_pos <- extract_M_table( + result_mat = result, + miout = miout, + sepc_mat = sepc_mat, + direction = "positive" + ) + + M_neg <- extract_M_table( + result_mat = result2, + miout = miout, + sepc_mat = sepc_mat2, + direction = "negative" + ) + + M_all <- rbind(M_pos, M_neg) + + T_all <- data.frame( + lhs = miout$lhs, + op = miout$op, + rhs = miout$rhs, + effect_positive = testeffect, + effect_negative = testeffect2, + stringsAsFactors = FALSE + ) + feasible <- TRUE + any_M <- any(resultall == "M", na.rm = TRUE) + + recommendation <- if (!feasible) { + "NOT APPLICABLE" + } else if (any_M) { + "NOT RECOMMENDED" + } else { + "RECOMMENDED" + } + + out <- list( + feasible = feasible, + any_M = any_M, + recommendation = recommendation, + M_table = M_all, + testeffect = T_all + ) + + class(out) <- "epcEquivCheckStd" + return(out) +} ## ---------------- ## Hidden Functions ## ---------------- -## totalFacVar: Find total factor variances when regression coeffient matrix -## and factor residual covariance matrix are specified +# lavInspectTotalVar() +# ------------------------------------------------------------------ +# Internal lavaan helper used by epcEquivFit() and EPC-related utilities. +# Extracts total variances of observed and latent variables from a +# fitted lavaan object using the implied covariance matrix. The +# function supports single-group and multi-group models and returns +# group-specific variance vectors for use in EPC standardization and +# variance-based transformations. No EPC estimation or statistical +# inference is performed. -totalFacVar <- function(beta, psi) { - ID <- diag(nrow(psi)) - total <- solve(ID - beta) %*% psi %*% t(solve(ID - beta)) - return(diag(total)) +##' @importFrom lavaan lavInspect +lavInspectTotalVar <- function(lavaanObj) { + result <- list() + nGroups <- lavInspect(lavaanObj, "ngroups") + cov.all <- lavInspect(lavaanObj, "cov.all") + if (nGroups == 1) cov.all <- list(cov.all) + for (i in 1:nGroups) { + temp <- diag(cov.all[[i]]) + names(temp) <- rownames(cov.all[[i]]) + result[[i]] <- temp + } + return(result) } -## findTotalVar: find the total indicator and factor variances +# lavInspectResidualVar() +# ------------------------------------------------------------------ +# Internal lavaan helper used by epcEquivFit() and EPC-related utilities. +# Extracts estimated residual variances for indicators (theta) and +# latent variables (psi) from a fitted lavaan object, handling both +# single-group and multi-group models. ##' @importFrom lavaan lavInspect -findTotalVar <- function(lavaanObj) { - result <- list() - nGroups <- lavInspect(lavaanObj, "ngroups") - cov.all <- lavInspect(lavaanObj, "cov.all") - if (nGroups == 1) cov.all <- list(cov.all) - for (i in 1:nGroups) { - temp <- diag(cov.all[[i]]) - names(temp) <- rownames(cov.all[[i]]) - result[[i]] <- temp - } - return(result) +lavInspectResidualVar <- function(lavaanObj) { + result <- list() + nGroups <- lavInspect(lavaanObj, "ngroups") + est <- lavInspect(lavaanObj, "est") + if (nGroups == 1) est <- list(est) + for (i in 1:nGroups) { + temppsi <- diag(est[[i]]$psi) + names(temppsi) <- rownames(est[[i]]$psi) + templambda <- est[[i]]$lambda + if(ncol(templambda) < nrow(templambda)) { + temptheta <- diag(est[[i]]$theta) + names(temptheta) <- rownames(est[[i]]$theta) + temp <- c(temptheta, temppsi) + } else { + temp <- temppsi + } + result[[i]] <- temp + } + return(result) } -## getTrivialEpc: find the trivial misspecified expected parameter changes -## given the type of parameters in each row of modification indices +# getTrivialEpc() +# ------------------------------------------------------------------ +# Internal utility used by epcEquivFit() and EPC equivalence diagnostics. +# Assigns operator-specific smallest effect size of interest (SESOI) +# values to each fixed parameter based on its type (e.g., factor +# loadings, structural regressions, latent covariances, residual +# covariances, intercepts). The resulting vector defines the magnitude +# of trivial misspecification used for EPC evaluation. +getTrivialEpc <- function( + mi, + lv_names, + stdLoad = 0.4, + corLatent = 0.1, + corResidual = 0.1, + stdBeta = 0.1, + stdIntcept = 0.2 +) { -getTrivialEpc <- function(mi, stdLoad=0.4, cor=0.1, stdBeta=0.1, intcept=0.2) { - op <- mi[,"op"] - result <- gsub("=~", stdLoad, op) - result <- gsub("~~", cor, result) - result <- gsub("~1", intcept, result) - result <- gsub("~", stdBeta, result) - return(result) -} + result <- numeric(nrow(mi)) + + for (i in seq_len(nrow(mi))) { + op <- mi[i, "op"] + lhs <- mi[i, "lhs"] + rhs <- mi[i, "rhs"] + + if (op == "=~") { + result[i] <- stdLoad + + } else if (op == "~~") { + if (lhs %in% lv_names && rhs %in% lv_names) { + result[i] <- corLatent + } else { + result[i] <- corResidual + } + + } else if (op == "~1") { + result[i] <- stdIntcept + + } else if (op == "~") { + result[i] <- stdBeta + + } else { + result[i] <- NA_real_ + } + } -## unstandardizeEpc: Transform from standardized EPC to unstandardized EPC - -unstandardizeEpc <- function(mi, delta, totalVar) { - name <- names(totalVar[[1]]) - lhsPos <- match(mi[,"lhs"], name) - rhsPos <- match(mi[,"rhs"], name) - group <- rep(1, nrow(mi)) - if("group" %in% colnames(mi)) group <- mi[,"group"] - getVar <- function(pos, group) totalVar[[group]][pos] - lhsVar <- mapply(getVar, pos=lhsPos, group=group) - rhsVar <- mapply(getVar, pos=rhsPos, group=group) - FUN <- function(op, lhsVar, rhsVar, delta) { - if(op == "|") return(NA) - lhsSD <- sqrt(lhsVar) - rhsSD <- sqrt(rhsVar) - if(!is.numeric(delta)) delta <- as.numeric(delta) - if(op == "=~") { - return((rhsSD * delta) / lhsSD) - } else if (op == "~~") { - return(lhsSD * delta * rhsSD) - } else if (op == "~1") { - return(lhsSD * delta) - } else if (op == "~") { - return((lhsSD * delta) / rhsSD) - } else { - return(NA) - } - } - unstdDelta <- mapply(FUN, op=mi[,"op"], lhsVar=lhsVar, rhsVar=rhsVar, delta=delta) - return(unstdDelta) + result } -## unstandardizeEpc: Transform from unstandardized EPC to standardized EPC. -## If delta is null, the unstandardized epc from the modification indices -## data.frame are used - -standardizeEpc <- function(mi, totalVar, delta = NULL) { - if(is.null(delta)) delta <- mi[,"epc"] - name <- names(totalVar[[1]]) - lhsPos <- match(mi[,"lhs"], name) - rhsPos <- match(mi[,"rhs"], name) - group <- rep(1, nrow(mi)) - if("group" %in% colnames(mi)) group <- mi[,"group"] - getVar <- function(pos, group) totalVar[[group]][pos] - lhsVar <- mapply(getVar, pos=lhsPos, group=group) - rhsVar <- mapply(getVar, pos=rhsPos, group=group) - FUN <- function(op, lhsVar, rhsVar, delta) { - lhsSD <- sqrt(lhsVar) - rhsSD <- sqrt(rhsVar) - if(!is.numeric(delta)) delta <- as.numeric(delta) - if(op == "=~") { - #stdload = beta * sdlatent / sdindicator = beta * lhs / rhs - return((delta / rhsSD) * lhsSD) - } else if (op == "~~") { - #r = cov / (sd1 * sd2) - return(delta / (lhsSD * rhsSD)) - } else if (op == "~1") { - #d = meanDiff/sd - return(delta / lhsSD) - } else if (op == "~") { - #beta = b * sdX / sdY = b * rhs / lhs - return((delta / lhsSD) * rhsSD) - } else { - return(NA) - } - } - stdDelta <- mapply(FUN, op=mi[,"op"], lhsVar=lhsVar, rhsVar=rhsVar, delta=delta) - return(stdDelta) + +# unstandardizeEpc() +# ------------------------------------------------------------------ +# Internal utility used by epcEquivFit() and related EPC diagnostics. +# Converts standardized effect-size thresholds (SESOI) back to the +# unstandardized EPC scale using total and residual variances of the +# involved variables. The transformation is operator-specific +# (e.g., loadings, regressions, covariances, intercepts) and provides +# unstandardized quantities required for EPC evaluation. +unstandardizeEpc <- function(mi, sesoi, totalVar, residualVar) { + name <- names(totalVar[[1]]) + lhsPos <- match(mi[,"lhs"], name) + rhsPos <- match(mi[,"rhs"], name) + group <- rep(1, nrow(mi)) + if("group" %in% colnames(mi)) group <- mi[,"group"] + getVar <- function(pos, group) totalVar[[group]][pos] + getVarRes <- function(pos, group) residualVar[[group]][pos] + lhsVar <- mapply(getVar, pos=lhsPos, group=group) + rhsVar <- mapply(getVar, pos=rhsPos, group=group) + lhsVarRes <- mapply(getVarRes, pos=lhsPos, group=group) + rhsVarRes <- mapply(getVarRes, pos=rhsPos, group=group) + + FUN <- function(op, lhsVar, rhsVar, lhsVarRes, rhsVarRes, sesoi) { + if(op == "|") return(NA) + lhsSD <- sqrt(lhsVar) + rhsSD <- sqrt(rhsVar) + lhsSDRes <- sqrt(lhsVarRes) + rhsSDRes <- sqrt(rhsVarRes) + if(!is.numeric(sesoi)) sesoi <- as.numeric(sesoi) + if(op == "=~") { + return((rhsSD * sesoi) / lhsSD) + } else if (op == "~~") { + return(lhsSDRes * sesoi * rhsSDRes) + } else if (op == "~1") { + return(lhsSD * sesoi) + } else if (op == "~") { + return((lhsSD * sesoi) / rhsSD) + } else { + return(NA) + } + } + sesoi <- mapply(FUN, op=mi[,"op"], lhsVar=lhsVar, rhsVar=rhsVar, + lhsVarRes=lhsVarRes, rhsVarRes=rhsVarRes, sesoi=sesoi) + return(sesoi) } -## decisionMIPow: provide the decision given the significance of modification -## indices and power to detect trivial misspecification +# standardizeEpc() +# ------------------------------------------------------------------ +# Internal utility used by epcEquivFit() and related EPC diagnostics. +# Transforms unstandardized EPCs or SESOI values into standardized +# effect-size metrics based on total and residual variances of the +# involved variables. The standardization is operator-specific +# (e.g., loadings, regressions, covariances, intercepts) and produces +# standardized quantities suitable for comparison against standardized +# SESOI thresholds. +standardizeEpc <- function(mi, totalVar, residualVar, sesoi = NULL) { + if(is.null(sesoi)) sesoi <- mi[,"epc"] + name <- names(totalVar[[1]]) + lhsPos <- match(mi[,"lhs"], name) + rhsPos <- match(mi[,"rhs"], name) + group <- rep(1, nrow(mi)) + if("group" %in% colnames(mi)) group <- mi[,"group"] + getVar <- function(pos, group) totalVar[[group]][pos] + getVarRes <- function(pos, group) residualVar[[group]][pos] + lhsVar <- mapply(getVar, pos=lhsPos, group=group) + rhsVar <- mapply(getVar, pos=rhsPos, group=group) + lhsVarRes <- mapply(getVarRes, pos=lhsPos, group=group) + rhsVarRes <- mapply(getVarRes, pos=rhsPos, group=group) + FUN <- function(op, lhsVar, rhsVar, lhsVarRes, rhsVarRes, sesoi) { + lhsSD <- sqrt(lhsVar) + rhsSD <- sqrt(rhsVar) + lhsSDRes <- sqrt(lhsVarRes) + rhsSDRes <- sqrt(rhsVarRes) + if(!is.numeric(sesoi)) sesoi <- as.numeric(sesoi) + if(op == "=~") { + #stdload = beta * sdlatent / sdindicator = beta * lhs / rhs + return((sesoi / rhsSD) * lhsSD) + } else if (op == "~~") { + #r = cov / (sd1 * sd2) + return(sesoi / (lhsSDRes * rhsSDRes)) + } else if (op == "~1") { + #d = meanDiff/sd + return(sesoi / lhsSD) + } else if (op == "~") { + #beta = b * sdX / sdY = b * rhs / lhs + return((sesoi / lhsSD) * rhsSD) + } else { + return(NA) + } + } + stdSesoi <- mapply(FUN, op=mi[,"op"], lhsVar=lhsVar, rhsVar=rhsVar, + lhsVarRes=lhsVarRes, rhsVarRes=rhsVarRes, sesoi=sesoi) + return(stdSesoi) +} +# decisionMIPow() +# ------------------------------------------------------------------ +# Core decision rule for power-based (SSV-style) EPC classification. +# Given the significance of the modification index, an indicator of +# high statistical power, the EPC value, and a smallest effect size +# of interest (SESOI), this function classifies the parameter into +# categories reflecting severity and substantive relevance +# (e.g., M, NM, I, EPC:M, EPC:NM). decisionMIPow <- function(sigMI, highPow, epc, trivialEpc) { - if(is.na(sigMI) | is.na(highPow)) return(NA) - if(sigMI & highPow) { - if(abs(epc) > abs(trivialEpc)) { - return("EPC:M") - } else { - return("EPC:NM") - } - } else if (sigMI & !highPow) { - return("M") - } else if (!sigMI & highPow) { - return("NM") - } else if (!sigMI & !highPow) { - return("I") - } else { - return(NA) - } + if(is.na(sigMI) | is.na(highPow)) return(NA) + if(sigMI & highPow) { + if(abs(epc) > abs(trivialEpc)) { + return("EPC:M") + } else { + return("EPC:NM") + } + } else if (sigMI & !highPow) { + return("M") + } else if (!sigMI & highPow) { + return("NM") + } else if (!sigMI & !highPow) { + return("I") + } else { + return(NA) + } } +# decisionCIEpc() +# ------------------------------------------------------------------ +# Core decision rule for CI-based EPC classification. +# Given an EPC confidence interval and a smallest effect size of +# interest (SESOI), this function classifies the parameter as +# Severe (M), Not Misspecified (NM), Inconclusive (I), or Underpowered (U). decisionCIEpc <- function(targetval, lower, upper, positiveonly = FALSE) { - if(is.na(lower) | is.na(upper)) return(NA) - if(positiveonly) { - if(lower > targetval) { - return("M") - } else if (upper < targetval) { - return("NM") - } else { - return("I") - } - } else { - negtargetval <- -targetval - if(lower > targetval | upper < negtargetval) { - return("M") - } else if (upper < targetval & negtargetval < lower) { - return("NM") - } else { - return("I") - } - } + if(is.na(lower) | is.na(upper)) return(NA) + if(positiveonly) { + ciwidth <- upper - lower + trivialwidth <- targetval + if (ciwidth > trivialwidth) { + return("U") + } else if (lower > targetval) { + return("M") + } else if (upper < targetval) { + return("NM") + } else { + return("I") + } + } else { + negtargetval <- -targetval + ciwidth <- upper - lower + trivialwidth <- 2*targetval + if(ciwidth > trivialwidth) { + return("U") + } else if(lower > targetval | upper < negtargetval) { + return("M") + } else if (upper < targetval & negtargetval < lower) { + return("NM") + } else { + return("I") + } + } +} + +# isPD() +# ------------------------------------------------------------------ +# Internal numerical utility used by epcEquivCheck() and related helpers. +# Checks whether a symmetric matrix is positive definite by verifying +# that all eigenvalues exceed a small tolerance. This function is used +# to screen implied population covariance matrices for admissibility +# before refitting models under perturbed parameters. +isPD <- function(M, tol = 1e-8) { + if (is.list(M)) { + return(all(vapply(M, isPD, logical(1), tol = tol))) + } + ev <- eigen(M, symmetric = TRUE, only.values = TRUE)$values + all(ev > tol) +} + +# cov2cor_safe() +# ------------------------------------------------------------------ +# Internal numerical utility used by epcEquivCheck() and related helpers. +# Converts a covariance matrix to a correlation matrix while safely +# handling zero or non-positive variances. The function computes +# correlations only for valid sub-blocks, enforces unit diagonals, +# and suppresses numerical artifacts to maintain stability in +# downstream standardized-parameter computations. +cov2cor_safe <- function(S) { + S <- as.matrix(S) + p <- nrow(S) + if (p != ncol(S)) stop("Input must be a square matrix.") + + v <- diag(S) + + # initialize correlation matrix + R <- matrix(0, p, p) + dimnames(R) <- dimnames(S) + + # indices with positive variance + pos <- which(v > 0) + + # normal cov2cor where possible + if (length(pos) > 0) { + sd <- sqrt(v[pos]) + R[pos, pos] <- S[pos, pos] / (sd %o% sd) + } + + # enforce symmetry and unit diagonal + R[!is.finite(R)] <- 0 + diag(R) <- 1 + + R +} + +# cor2cov_safe() +# ------------------------------------------------------------------ +# Internal numerical utility used by epcEquivCheck() and related helpers. +# Converts a correlation matrix to a covariance matrix given a vector +# of standard deviations, while safely handling zero, negative, or +# non-finite standard deviations. The function operates only on valid +# sub-blocks, enforces symmetry, and removes numerical artifacts to +# ensure stability in subsequent covariance-based computations. +cor2cov_safe <- function(R, sd) { + R <- as.matrix(R) + p <- nrow(R) + if (p != ncol(R)) stop("R must be square.") + if (length(sd) != p) stop("Length of sd must match dimension of R.") + + # initialize covariance matrix + S <- matrix(0, p, p) + dimnames(S) <- dimnames(R) + + # indices with positive sd + pos <- which(sd > 0 & is.finite(sd)) + + if (length(pos) > 0) { + D <- diag(sd[pos], length(pos)) + S[pos, pos] <- D %*% R[pos, pos, drop = FALSE] %*% D + } + + # enforce symmetry and clean numerical noise + S[!is.finite(S)] <- 0 + S +} + +# extract_M_table() +# ------------------------------------------------------------------ +# Internal helper used by epcEquivCheck(). +# Extracts and formats combinations of perturbed and resulting fixed +# parameters that yield Severe (M) EPC decisions. The function maps +# row–column indices from a decision matrix back to the corresponding +# parameter labels in the \code{epcEquivFit} output and returns a tidy summary +# table for reporting purposes only. +extract_M_table <- function(result_mat, miout, sepc_mat, direction) { + + idx <- which(result_mat == "M", arr.ind = TRUE) + if (nrow(idx) == 0) return(NULL) + + data.frame( + perturbed_lhs = miout$lhs[idx[, "row"]], + perturbed_op = miout$op[idx[, "row"]], + perturbed_rhs = miout$rhs[idx[, "row"]], + + resulting_lhs = miout$lhs[idx[, "col"]], + resulting_op = miout$op[idx[, "col"]], + resulting_rhs = miout$rhs[idx[, "col"]], + + resulting_sepc = sepc_mat[cbind(idx[, "row"], idx[, "col"])], + + direction = direction, + stringsAsFactors = FALSE + ) +} + +# print.epcEquivCheckStd() +# ------------------------------------------------------------------ +# Internal print method for epcEquivCheckStd objects. +# Formats and displays feasibility and recommendation results from +# standardized-parameter EPC equivalence checks. + +#' @export +print.epcEquivCheckStd <- function(x, ...) { + cat("EPC Equivalence Feasibility (Standardized Parameters)\n") + cat("----------------------------------------------------\n") + + cat("Feasible standardized population:", x$feasible, "\n") + cat("Any EPC exceeding SESOI:", x$any_M, "\n") + cat("Recommendation:", x$recommendation, "\n\n") + + if (x$recommendation == "NOT RECOMMENDED") { + cat("Non-equivalent EPCs detected (summary):\n") + print(x$M_table) + } else if (x$recommendation == "RECOMMENDED") { + cat("No EPC exceeded the SESOI under tested misspecifications.\n") + } else { + cat("Standardized parameters do not define a valid population model.\n") + } + + invisible(x) } +# summary.epcequivfit.data.frame() +# ------------------------------------------------------------------ +# Internal summary constructor used by epcEquivFit(). +# Aggregates and ranks fixed-parameter EPC results from a data-frame +# representation into a structured summaryEpcEquivFit object for +# downstream printing and global decision reporting. This function +# performs data processing only and does not conduct EPC estimation +# or statistical inference. + +#' @method summary epcequivfit.data.frame +#' @rdname epcEquivFit +#' @param object An object returned by \code{\link{epcEquivFit}}. +#' @param top Number of top-ranked EPCs to display. +#' @param ssv Logical; whether to include power-based diagnostics. +#' @export +summary.epcequivfit.data.frame <- function(object, ..., top = 5, ssv = FALSE) { + miout <- object + stopifnot(is.data.frame(miout)) + num_cols <- c( + "std.epc", "std.sesoi", + "lower.std.epc", "upper.std.epc", + "se.epc", "pow" + ) + + for (nm in num_cols) { + if (!is.null(miout[[nm]])) { + miout[[nm]] <- as.numeric(as.character(miout[[nm]])) + } + } + + miout$std.sesoi[miout$std.sesoi == 0] <- NA + miout$se.epc[miout$se.epc == 0] <- NA + + + # ---- Derived quantities ---- + miout$severity <- abs(miout$std.epc / miout$std.sesoi) + + miout$ci_gap <- with(miout, pmin( + abs(lower.std.epc - std.sesoi), + abs(upper.std.epc - std.sesoi), + abs(lower.std.epc + std.sesoi), + abs(upper.std.epc + std.sesoi), + na.rm = TRUE + )) + + # ---- EPC equivalence testing (CI-based; primary) ---- + epc_equiv <- list( + n_M = sum(miout$decision.ci == "M", na.rm = TRUE), + n_I = sum(miout$decision.ci == "I", na.rm = TRUE), + n_U = sum(miout$decision.ci == "U", na.rm = TRUE), + n_NM = sum(miout$decision.ci == "NM", na.rm = TRUE), + any_M = any(miout$decision.ci == "M", na.rm = TRUE) + ) + + # Severe EPCs + top_M <- miout[miout$decision.ci == "M", ] + top_M <- top_M[is.finite(top_M$severity), ] + top_M <- top_M[order(top_M$severity, decreasing = TRUE), ] + if (nrow(top_M) > top) top_M <- top_M[1:top, ] + + # CI-inconclusive EPCs (exclude NM) + top_I_ci <- miout[miout$decision.ci == "I", ] + top_I_ci <- top_I_ci[is.finite(top_I_ci$ci_gap), ] + top_I_ci <- top_I_ci[order(top_I_ci$ci_gap, decreasing = TRUE), ] + if (nrow(top_I_ci) > top) top_I_ci <- top_I_ci[1:top, ] + + # ---- SSV / power-based diagnostics (secondary) ---- + ssv_out <- list( + n_M = sum(miout$decision.pow %in% c("M", "EPC:M"), na.rm = TRUE), + n_I = sum(miout$decision.pow == "I", na.rm = TRUE), + n_NM = sum(miout$decision.pow %in% c("NM", "EPC:NM"), na.rm = TRUE), + any_M = any(miout$decision.pow %in% c("M", "EPC:M"), na.rm = TRUE) + ) + + # Severe SSV + top_M_pow <- miout[miout$decision.pow == "M", ] + top_M_pow <- top_M_pow[is.finite(top_M_pow$severity), ] + top_M_pow <- top_M_pow[order(top_M_pow$severity, decreasing = TRUE), ] + if (nrow(top_M_pow) > top) top_M_pow <- top_M_pow[1:top, ] + + # Underpowered EPCs (exclude NM) + top_I_pow <- miout[miout$decision.pow == "I", ] + top_I_pow <- top_I_pow[is.finite(top_I_pow$pow), ] + top_I_pow <- top_I_pow[order(top_I_pow$pow), ] # lowest power first + if (nrow(top_I_pow) > top) top_I_pow <- top_I_pow[1:top, ] + + out <- list( + epc_equivalence = epc_equiv, + ssv = ssv_out, + top_non_equiv = top_M, + top_inconclusive_ci = top_I_ci, + top_non_pow = top_M_pow, + top_inconclusive_pow = top_I_pow, + show_ssv = ssv + ) + + class(out) <- "summaryEpcEquivFit" + out +} + +# global_localfit_decision() +# ------------------------------------------------------------------ +# Internal decision utility used by print.summaryEpcEquivFit(). +# Aggregates local fixed-parameter classifications (e.g., Severe, +# Inconclusive, Underpowered, Not Misspecified) into a single global +# decision using a conservative priority rule. +global_localfit_decision <- function(n_M, n_I, n_U, n_NM) { + + if (n_U > 0) { + return("UNDERPOWERED") + } + + if (n_M > 0) { + return("SEVERE MISSPECIFICATION") + } + + if (n_I > 0) { + return("INCONCLUSIVE") + } + + if (n_NM > 0) { + return("EQUIVALENT") + } + + return("NO PARAMETERS EVALUATED") +} + +# print.summaryEpcEquivFit() +# ------------------------------------------------------------------ +# Internal print method for summaryEpcEquivFit objects. +# Formats and displays results from global EPC equivalence testing +# and the Saris, Satorra, and van der Veld (2009) power-based method. +# This function summarizes aggregate fixed-parameter diagnostics +# for display purposes only and is not part of the analytical +# implementation of epcEquivFit(). + +#' @export +print.summaryEpcEquivFit <- function(x, ...) { + + cat("Global EPC Evaluation Summary\n") + cat("--------------------------------\n\n") + + ## ---- EPC equivalence testing (primary) ---- + cat("[1. EPC Equivalence Testing: CI-based]\n") + cat("Severe (M):", x$epc_equivalence$n_M, "\n") + cat("Inconclusive (I):", x$epc_equivalence$n_I, "\n") + cat("CI-Underpowered (U):", x$epc_equivalence$n_U, "\n") + cat("Trivial / Not Misspecified (NM):", x$epc_equivalence$n_NM, "\n\n") + + global_epc <- global_localfit_decision( + n_M = x$epc_equivalence$n_M, + n_I = x$epc_equivalence$n_I, + n_U = x$epc_equivalence$n_U, + n_NM = x$epc_equivalence$n_NM + ) + + cat("Global EPC Equivalence Decision:", global_epc, "\n\n") + + if (x$epc_equivalence$n_M > 0) { + cat("1.1 Top Severe EPCs (ranked by |std.epc / SESOI|):\n") + print( + x$top_non_equiv[ + , c("lhs","op","rhs","std.epc","std.sesoi","severity"), + drop = FALSE + ] + ) + cat("\n") + } + + if (x$epc_equivalence$n_I > 0) { + cat("1.2 Top CI-inconclusive EPCs\n") + cat("(ranked by distance to equivalence bounds; larger = needs narrower CI):\n") + print( + x$top_inconclusive_ci[ + , c("lhs","op","rhs","lower.std.epc","upper.std.epc","ci_gap"), + drop = FALSE + ] + ) + cat("\n") + } + + if(isTRUE(x$show_ssv)) { + ## ---- SSV / power-based diagnostics (secondary) ---- + cat("[2. Saris, Satorra, Van der Veld (2009) / Power-based Diagnostics]\n") + cat("Severe (M, EPC:M):", x$ssv$n_M, "\n") + cat("Inconclusive (I):", x$ssv$n_I, "\n") + cat("Trivial / Not Misspecified (NM, EPC:NM):", x$ssv$n_NM, "\n\n") + + global_ssv <- global_localfit_decision( + n_M = x$ssv$n_M, + n_I = x$ssv$n_I, + n_U = 0, + n_NM = x$ssv$n_NM + ) + + cat("Global SSV / Power-based Decision:", global_ssv, "\n\n") + + if (x$ssv$n_M > 0) { + cat("2.1 Top severe EPCs (ranked by |std.epc / SESOI|):\n") + print( + x$top_non_pow[ + , c("lhs","op","rhs","std.epc","std.sesoi","severity"), + drop = FALSE + ] + ) + cat("\n") + } + + if (x$ssv$n_I > 0) { + cat("2.2 Top inconclusive EPCs\n") + cat("(ranked by lowest approximate power):\n") + print( + x$top_inconclusive_pow[ + , c("lhs","op","rhs","pow","se.epc"), + drop = FALSE + ] + ) + } + } + + invisible(x) +} + +# augmentBetaWithGamma() +# ------------------------------------------------------------------ +# Copied from simsem. Internal utility used by epcEquivCheck(). +# Augments the factor regression coefficient matrix (beta) with +# regression paths from exogenous observed variables (gamma) so that +# endogenous and exogenous predictors can be handled within a single +# block-recursive system. +# NOTE: This function is prepared for future support of exogenous covariates +augmentBetaWithGamma <- function(beta, gamma) { + nf <- nrow(beta) + nz <- ncol(gamma) + result <- matrix(0, nf + nz, nf + nz) + result[(nz + 1):(nz + nf), (nz + 1):(nz + nf)] <- beta + result[(nz + 1):(nz + nf), 1:nz] <- gamma + result +} + +# augmentPsiWithExogenousCov() +# ------------------------------------------------------------------ +# Copied from simsem. Internal utility used by epcEquivCheck(). +# Augments the factor residual covariance matrix (psi) with the +# covariance matrix of exogenous observed variables (sigmaxx), +# producing a joint covariance structure compatible with the +# augmented regression system. +# NOTE: This function is prepared for future support of exogenous covariates +augmentPsiWithExogenousCov <- function(psi, sigmaxx) { + nf <- nrow(psi) + nz <- ncol(sigmaxx) + result <- matrix(0, nf + nz, nf + nz) + result[(nz + 1):(nz + nf), (nz + 1):(nz + nf)] <- psi + result[1:nz, 1:nz] <- sigmaxx + result +} + +# findRecursiveSet() +# ------------------------------------------------------------------ +# Copied from simsem. Internal helper for epcEquivCheck(). +# Identifies a recursive (block-triangular) ordering of variables +# implied by the regression coefficient matrix (beta). This ordering +# is required to compute analytical residual variances in recursive +# SEMs. +findRecursiveSet <- function(beta) { + if (any(diag(beta) != 0)) + stop("Diagonal elements of beta must be zero.") + result <- list() + ni <- nrow(beta) + fix.variable <- rep(FALSE, ni) + ni.sofar <- 0 + i <- 1 + while (ni.sofar < ni) { + temp <- findRowZero(beta, fix.variable) + if (is.null(temp)) + stop("The matrix is not recursive.") + fix.variable[temp] <- TRUE + result[[i]] <- temp + i <- i + 1 + ni.sofar <- ni.sofar + length(temp) + } + return(result) +} + +# findRowZero() +# ------------------------------------------------------------------ +# Internal helper for recursive SEM utilities (copied/adapted from simsem). +# Identifies rows of a square coefficient matrix whose unfixed entries +# are all zero. This is used when constructing a recursive (block-triangular) +# ordering of variables implied by a regression matrix (e.g., beta). +# +# Rows already marked as fixed are ignored, and the function returns the +# indices of rows that contain only zeros among the remaining unfixed columns. +# If no such rows exist, NULL is returned. +# +# This function is used by findRecursiveSet() to iteratively determine +# which variables can be treated as exogenous at each recursion step. +findRowZero <- function(square.matrix, is.row.fixed = FALSE) { + ni <- nrow(square.matrix) + if (length(is.row.fixed) == 1) { + if (is.row.fixed == FALSE) + is.row.fixed <- rep(FALSE, ni) + } + result <- NULL + desired.zero <- sum(!is.row.fixed) + for (i in 1:ni) { + if (is.row.fixed[i] == FALSE) { + temp <- sum(square.matrix[i, !is.row.fixed] == 0, na.rm = TRUE) + if (temp == desired.zero) + result <- c(result, i) + } + } + return(result) +} + +# findIndResidualVar() +# ------------------------------------------------------------------ +# Copied from simsem. Internal utility used by epcEquivCheck(). +# Computes indicator residual variances given factor loadings, total +# factor covariance, and total indicator variances. The resulting +# implied residual variances are used when constructing EPC values +# under trivial misspecification based on a fitted lavaan object. +findIndResidualVar <- function(lambda, totalFactorCov, totalVarTheta = NULL, + kappa = NULL, covcov = NULL) { + ni <- nrow(lambda) + if (is.null(totalVarTheta)) totalVarTheta <- rep(1, ni) + factor.part <- lambda %*% totalFactorCov %*% t(lambda) + error.var <- totalVarTheta - pmax(diag(factor.part), 0) + if (!is.null(kappa)) error.var <- error.var - diag(kappa %*% covcov %*% t(kappa)) + error.var[(error.var < 0) & (totalVarTheta == 0)] <- 0 + return(as.vector(error.var)) +} + +# findFactorTotalCov() +# ------------------------------------------------------------------ +# Copied from simsem. Internal utility used by epcEquivCheck(). +# Computes the total covariance matrix of latent factors implied by +# a recursive regression structure and factor residual covariances. +# The resulting implied factor covariance matrix is used as an +# intermediate quantity for evaluating EPCs under trivial +# misspecification based on a fitted lavaan object. +findFactorTotalCov <- function(beta, psi = NULL, corPsi = NULL, + totalVarPsi = NULL, errorVarPsi = NULL, + gamma = NULL, covcov = NULL) { + if (is.null(psi)) { + if (is.null(errorVarPsi)) + errorVarPsi <- findFactorResidualVar(beta, corPsi, totalVarPsi) + psi <- suppressWarnings(cor2cov_safe(as.matrix(corPsi), sqrt(errorVarPsi))) + } + iden <- diag(nrow(beta)) + inv <- solve(iden - beta) + facTotalCov <- inv %*% psi %*% t(inv) + if (!is.null(gamma)) { + facTotalCov <- facTotalCov + (inv %*% gamma %*% covcov %*% t(gamma) %*% t(inv)) + } + return(facTotalCov) +} + +# findFactorResidualVar() +# ------------------------------------------------------------------ +# Copied from simsem. Core analytical routine used by epcEquivCheck(). +# Computes factor residual variances implied by a recursive SEM given +# total variances, factor correlations, and regression coefficients. +# These analytically derived residual variances provide the basis +# for evaluating EPCs under trivial misspecification based on a +# fitted lavaan object. +findFactorResidualVar <- function(beta, corPsi, totalVarPsi = NULL, + gamma = NULL, covcov = NULL) { + stopifnot(nrow(beta) == ncol(beta)) + stopifnot(nrow(corPsi) == ncol(corPsi)) + stopifnot(nrow(beta) == nrow(corPsi)) + if (!is.null(gamma)) { + if (is.null(totalVarPsi)) totalVarPsi <- rep(1, nrow(beta)) + beta <- augmentBetaWithGamma(beta, gamma) + corPsi <- augmentPsiWithExogenousCov(corPsi, cov2cor(covcov)) + totalVarPsi <- c(diag(covcov), totalVarPsi) + } + if (sum(diag(corPsi)) == 0) diag(corPsi) <- 1 + ni <- nrow(beta) + set <- findRecursiveSet(beta) + errorVar <- rep(1, ni) + if (is.null(totalVarPsi)) totalVarPsi <- rep(1, ni) + errorVar[set[[1]]] <- totalVarPsi[set[[1]]] + iv <- NULL + ivCor <- corPsi[set[[1]], set[[1]]] + startVar <- totalVarPsi[set[[1]]] + ivCov <- suppressWarnings(cor2cov_safe(as.matrix(ivCor), sqrt(startVar))) + for (i in seq_len(length(set) - 1)) { + iv <- c(iv, set[[i]]) + dv <- set[[i + 1]] + tempBeta <- matrix(beta[dv, iv], nrow = length(dv), ncol = length(iv)) + var.reg <- (tempBeta %*% ivCov %*% t(tempBeta)) + tempPsi <- corPsi[dv, dv] + tempPsiSd <- rep(0, length(dv)) + for (j in 1:length(dv)) { + errorVar[dv[j]] <- totalVarPsi[dv[j]] - var.reg[j, j] + if (is.na(errorVar[dv[j]]) || errorVar[dv[j]] < 0) { + tempPsiSd[j] <- NA_real_ + } else { + tempPsiSd[j] <- sqrt(errorVar[dv[j]]) + } + } + if (i < (length(set) - 1)) { + tempPsi <- suppressWarnings(cor2cov_safe(as.matrix(tempPsi), as.matrix(tempPsiSd))) #### FLAG CHANGES + real.tempPsi <- matrix(0, length(iv) + length(dv), length(iv) + length(dv)) + real.tempPsi[1:length(iv), 1:length(iv)] <- ivCov + real.tempPsi[(length(iv) + 1):(length(iv) + length(dv)), (length(iv) + + 1):(length(iv) + length(dv))] <- tempPsi + agg <- c(iv, dv) + blank.path <- matrix(0, nrow = length(iv), ncol = length(agg)) + temp.path2 <- beta[dv, agg] + temp.path2 <- rbind(blank.path, temp.path2) + ID <- matrix(0, length(agg), length(agg)) + diag(ID) <- 1 + ivCov <- solve(ID - temp.path2) %*% real.tempPsi %*% t(solve(ID - temp.path2)) + } + } + if (!is.null(gamma)) { + errorVar <- errorVar[(nrow(covcov) + 1):length(errorVar)] + } + return(as.vector(errorVar)) +} + +# lavCheckAdmissibleFit() +# ------------------------------------------------------------------ +# Internal admissibility and numerical-sanity check for fitted lavaan +# objects. This function verifies that a fitted model is suitable for +# downstream analytical use (e.g., EPC evaluation, population-based +# refitting) by screening for common estimation pathologies. +# +# The checks include: +# - model convergence +# - positive definiteness of observed, implied, residual, and latent +# covariance matrices +# - strictly positive residual and latent variances +# - positive definiteness of the parameter vcov matrix (if available) +# - positive definiteness of the Hessian matrix (if available) +# +# The function is intentionally conservative: failure of any check +# returns FALSE. It is not intended as a model-fit diagnostic, but as a +# numerical admissibility gate for internal algorithms. +lavCheckAdmissibleFit <- function(fit, tol = 1e-8) { + cov.ov <- lavaan::lavInspect(fit, "cov.ov") # sample cov / polychoric + sigma <- lavaan::lavInspect(fit, "sigma.hat") # model-implied cov + theta <- lavaan::lavInspect(fit, "theta") # residual covariances + psi <- lavaan::lavInspect(fit, "cov.lv") # latent covariances + vc <- lavaan::lavInspect(fit, "vcov") # Variance-Covariance Matrix + # hessian may not exist (e.g., se = "none") + hess <- try(lavaan::lavInspect(fit, "hessian"), silent = TRUE) + hess_pd <- if (inherits(hess, "try-error") || is.null(hess)) NA else isPD(hess, tol) + + vccheck <- TRUE + if(!is.null(vc)) { + if(all(round(vc,4) == 0)) { + vccheck <- TRUE + } else { + vccheck <- isPD(vc, tol) + } + } + all( + converged = lavaan::lavInspect(fit, "converged"), + + cov.ov.pd = isPD(cov.ov, tol), + sigma.pd = isPD(sigma, tol), + + theta.pd = isPD(theta, tol), + psi.pd = isPD(psi, tol), + + vc = vccheck, + + theta.var.pos = diag_pos(theta, tol), # no negative residual variances + psi.var.pos = diag_pos(psi, tol), # no negative latent variances + + hessian.pd = hess_pd + ) +} + +# diag_pos() +# ------------------------------------------------------------------ +# Internal numerical utility for model-admissibility checks. +# Tests whether all diagonal elements of a matrix are strictly positive +# up to a numerical tolerance. This is used to screen for negative or +# zero variances in residual (theta) and latent (psi) covariance +# matrices returned by lavaan. +# +# If a list of matrices is supplied (e.g., multigroup models), the +# check is applied recursively to each element. +diag_pos <- function(M, tol = 1e-8) { + if (is.list(M)) { + return(all(vapply(M, diag_pos, logical(1), tol = tol))) + } + all(diag(M) > tol) +} diff --git a/semTools/R/semTools.R b/semTools/R/semTools.R index bd8575f..49e4881 100644 --- a/semTools/R/semTools.R +++ b/semTools/R/semTools.R @@ -15,7 +15,7 @@ ##' [moreFitIndices()], ##' [nullRMSEA()], ##' [singleParamTest()], -##' [miPowerFit()], and +##' [epcEquivFit()], and ##' [chisqSmallN()]} ##' \item{Measurement Invariance: ##' [measEq.syntax()], diff --git a/semTools/man/clipboard.Rd b/semTools/man/clipboard.Rd index 9054965..91a0d3c 100644 --- a/semTools/man/clipboard.Rd +++ b/semTools/man/clipboard.Rd @@ -17,14 +17,14 @@ saveFile(object, file, what = "summary", tableFormat = FALSE, \item{what}{The attributes of the \code{lavaan} object to be copied in the clipboard. \code{"summary"} is to copy the screen provided from the -\code{summary} function. \code{"mifit"} is to copy the result from the -\code{\link[=miPowerFit]{miPowerFit()}} function. Other attributes listed in the +\code{summary} function. \code{"epceqfit"} is to copy the result from the +\code{\link[=epcEquivFit]{epcEquivFit()}} function. Other attributes listed in the \code{inspect} method in the \link[lavaan:lavaan-class]{lavaan::lavaan} could also be used, such as \code{"coef"}, \code{"se"}, \code{"fit"}, \code{"samp"}, and so on. Ignored for \linkS4class{FitDiff}-class objects.} \item{\dots}{Additional arguments when passing a \code{lavaan} object to the -\code{summary} or \code{\link[=miPowerFit]{miPowerFit()}} function.} +\code{summary} or \code{\link[=epcEquivFit]{epcEquivFit()}} function.} \item{file}{A file name used for saving the result.} @@ -66,8 +66,8 @@ clipboard(fit) # pass additional arguments to summary() method for class?lavaan clipboard(fit, rsquare = TRUE, standardized = TRUE, fit.measures = TRUE) -# Copy modification indices and fit stats from the miPowerFit() function -clipboard(fit, "mifit") +# Copy the EPC equivalence testing results from the epcEquivFit() function +clipboard(fit, "epceqfit") # Copy the parameter estimates clipboard(fit, "coef") @@ -84,8 +84,8 @@ clipboard(fit, "fit") # Save the summary of the lavaan object saveFile(fit, "out.txt") -# Save modification indices and fit stats from the miPowerFit() function -saveFile(fit, "out.txt", "mifit") +# Save the EPC equivalence testing results from the epcEquivFit() function +saveFile(fit, "out.txt", "epceqfit") # Save the parameter estimates saveFile(fit, "out.txt", "coef") diff --git a/semTools/man/epcEquivCheck.Rd b/semTools/man/epcEquivCheck.Rd new file mode 100644 index 0000000..f5505aa --- /dev/null +++ b/semTools/man/epcEquivCheck.Rd @@ -0,0 +1,90 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/miPowerFit.R +\name{epcEquivCheck} +\alias{epcEquivCheck} +\title{EPC Equivalence Feasibility Check for Standardized Parameters} +\usage{ +epcEquivCheck(lavaanObj, minRelEffect = 0.75, stdLoad = 0.4, cor = 0.1, + corLatent = NULL, corResidual = NULL, stdBeta = 0.1) +} +\arguments{ +\item{lavaanObj}{A fitted \code{lavaan} object representing the target model.} + +\item{minRelEffect}{A scalar in (0, 1) specifying the minimum relative +magnitude of the standardized perturbation to be evaluated. The +default value of 0.75 indicates that perturbations equal to 75\\% of +the SESOI are treated as trivial. If EPCs exceed the SESOI under +such perturbations, EPC equivalence testing is not recommended.} + +\item{stdLoad}{Standardized factor loading used to define the SESOI +for loading misspecifications.} + +\item{cor}{Standardized correlation used as a default SESOI for +covariance misspecifications. This value is used for both latent +and residual covariances unless overridden by +\code{corLatent} or \code{corResidual}.} + +\item{corLatent}{Standardized latent factor correlation used to +define the SESOI for latent covariance misspecifications. If +\code{NULL}, defaults to \code{cor}.} + +\item{corResidual}{Standardized residual correlation used to define +the SESOI for indicator residual covariance misspecifications. If +\code{NULL}, defaults to \code{cor}.} + +\item{stdBeta}{Standardized regression coefficient used to define +the SESOI for structural misspecifications.} +} +\value{ +An object of class \code{"epcEquivCheckStd"} containing: +\itemize{ +\item \code{feasible}: Logical indicator of whether a valid +standardized population model exists. +\item \code{any_M}: Logical indicator of whether any EPC exceeded +the SESOI under the evaluated misspecifications. +\item \code{recommendation}: Character string summarizing feasibility +(e.g., \code{"RECOMMENDED"}, \code{"NOT RECOMMENDED"}). +\item \code{M_table}: Data frame summarizing EPCs exceeding the SESOI, +if any. +\item \code{testeffect}: Data frame reporting the smallest tested +standardized perturbations in each direction. +} +} +\description{ +Performs an EPC-based feasibility check to assess whether a set of +standardized population parameters defines a valid population +covariance matrix and whether trivially misspecified parameters +remain within a user-defined smallest effect size of interest (SESOI). +Feasibility is evaluated by constructing implied population models +under targeted parameter perturbations and examining EPC behavior +using \code{\link{epcEquivFit}}. +} +\details{ +This function focuses on standardized parameters and supports +recursive SEMs with continuous indicators only. + +The procedure first checks whether the standardized parameters imply +a positive definite population covariance matrix. It then evaluates +EPC behavior under both positive and negative trivial +misspecifications by repeatedly constructing implied population +covariance matrices with perturbed parameters +(\code{minRelEffect} \eqn{\times} SESOI), refitting the model, and +re-evaluating EPCs. + +Models with categorical indicators, formative indicators, or +multiple-group structures are not supported. +} +\examples{ + +library(lavaan) + +one.model <- ' onefactor =~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 ' +fit <- cfa(one.model, data = HolzingerSwineford1939) +\donttest{ +epcEquivCheck(fit) +} + +} +\seealso{ +\code{\link{epcEquivFit}} +} diff --git a/semTools/man/epcEquivFit.Rd b/semTools/man/epcEquivFit.Rd new file mode 100644 index 0000000..13af1f9 --- /dev/null +++ b/semTools/man/epcEquivFit.Rd @@ -0,0 +1,127 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/miPowerFit.R +\name{epcEquivFit} +\alias{epcEquivFit} +\alias{miPowerFit} +\alias{summary.epcequivfit.data.frame} +\title{EPC Equivalence Fit Evaluation Using Modification Indices} +\usage{ +epcEquivFit(lavaanObj, stdLoad = 0.4, cor = 0.1, corLatent = NULL, + corResidual = NULL, stdBeta = 0.1, stdIntcept = 0.2, stdSesoi = NULL, + sesoi = NULL, cilevel = 0.9, ...) + +\method{summary}{epcequivfit.data.frame}(object, ..., top = 5, ssv = FALSE) +} +\arguments{ +\item{lavaanObj}{A fitted \code{lavaan} object used to evaluate model fit.} + +\item{stdLoad}{Standardized factor loading defining the SESOI for +loading misspecifications. Default is 0.4.} + +\item{cor}{Default standardized correlation defining the SESOI for +covariance misspecifications. Used for both latent and residual +covariances unless overridden.} + +\item{corLatent}{Standardized latent factor correlation defining the +SESOI for latent covariance misspecifications. If \code{NULL}, +defaults to \code{cor}.} + +\item{corResidual}{Standardized residual correlation defining the +SESOI for residual covariance misspecifications. If \code{NULL}, +defaults to \code{cor}.} + +\item{stdBeta}{Standardized regression coefficient defining the SESOI +for structural misspecifications. Default is 0.1.} + +\item{stdIntcept}{Standardized intercept (Cohen's \emph{d}) defining +the SESOI for intercept misspecifications. Default is 0.2.} + +\item{stdSesoi}{Optional vector of standardized SESOI values. If +provided, overrides operator-specific SESOI definitions.} + +\item{sesoi}{Optional vector of unstandardized SESOI values. If +provided, overrides \code{stdSesoi} and all operator-specific SESOI +arguments.} + +\item{cilevel}{Confidence level for EPC confidence intervals used in +CI-based equivalence testing.} + +\item{\dots}{Additional arguments passed to +\code{\link[lavaan]{modificationIndices}}.} + +\item{object}{An object returned by \code{\link{epcEquivFit}}.} + +\item{top}{Number of top-ranked EPCs to display.} + +\item{ssv}{Logical; whether to include power-based diagnostics.} +} +\value{ +A data frame with one row per fixed parameter, containing: +\enumerate{ +\item Parameter identifiers: \code{lhs}, \code{op}, \code{rhs}, and \code{group}. +\item Modification index (\code{mi}) and expected parameter change estimates (\code{epc}). +\item Unstandardized and standardized smallest effect size of interest values +(\code{sesoi}, \code{std.sesoi}). +\item Power-based decision (\code{decision.pow}) and related diagnostics, including +whether the modification index is statistically significant +(\code{significant.mi}) and whether the misfit at the SESOI has power greater than 0.80 +(\code{high.power}). +\item EPC-related statistics, including the standard error of the EPC +(\code{se.epc}), confidence interval bounds for the EPC +(\code{lower.epc}, \code{upper.epc}), and confidence interval bounds for the +standardized EPC (\code{lower.std.epc}, \code{upper.std.epc}). +\item Confidence-interval–based equivalence decision (\code{decision.ci}). +} +} +\description{ +Evaluates model fit from an equivalence-testing perspective by +aggregating local EPC-based diagnostics into a global, fit-style +assessment. The procedure combines modification indices (MI), +expected parameter changes (EPC), statistical power, and confidence +intervals relative to a smallest effect size of interest (SESOI). +} +\details{ +Two complementary local decision rules are implemented: + +\strong{Method 1 (Power-based; Saris, Satorra, & van der Veld, 2009).} +Modification indices, statistical power, and EPC magnitude are jointly +evaluated (the J-rule) to classify fixed parameters as misspecified, +not misspecified, or inconclusive. + +\strong{Method 2 (CI-based equivalence testing).} +Confidence intervals of EPCs are compared against a trivial +misspecification region defined by the SESOI to determine whether +fixed parameters are severely misspecified, trivially misspecified, +underpowered, or inconclusive. + +The resulting local classifications are returned in a single data +frame and can be summarized to yield a global equivalence-style fit +evaluation. + +This function provides a local-to-global equivalence-based alternative +to traditional exact-fit evaluation. It is designed to assess whether +fixed parameters are substantively misspecified relative to a SESOI, +rather than whether a model fits exactly. + +Models with categorical indicators or unsupported constraints may +not be fully supported. +} +\examples{ + +library(lavaan) + +one.model <- ' onefactor =~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 ' +fit <- cfa(one.model, data = HolzingerSwineford1939) +out <- epcEquivFit(fit) +out +summary(out) + +} +\references{ +Saris, W. E., Satorra, A., & van der Veld, W. M. (2009). +Testing structural equation models or detection of misspecifications? +\emph{Structural Equation Modeling, 16}(4), 561--582. +} +\seealso{ +\code{\link{epcEquivCheck}} +} diff --git a/semTools/man/miPowerFit.Rd b/semTools/man/miPowerFit.Rd deleted file mode 100644 index d59f56b..0000000 --- a/semTools/man/miPowerFit.Rd +++ /dev/null @@ -1,192 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/miPowerFit.R -\name{miPowerFit} -\alias{miPowerFit} -\title{Modification indices and their power approach for model fit evaluation} -\usage{ -miPowerFit(lavaanObj, stdLoad = 0.4, cor = 0.1, stdBeta = 0.1, - intcept = 0.2, stdDelta = NULL, delta = NULL, cilevel = 0.9, ...) -} -\arguments{ -\item{lavaanObj}{The lavaan model object used to evaluate model fit} - -\item{stdLoad}{The amount of standardized factor loading that one would like -to be detected (rejected). The default value is 0.4, which is suggested by -Saris and colleagues (2009, p. 571).} - -\item{cor}{The amount of factor or error correlations that one would like to -be detected (rejected). The default value is 0.1, which is suggested by -Saris and colleagues (2009, p. 571).} - -\item{stdBeta}{The amount of standardized regression coefficients that one -would like to be detected (rejected). The default value is 0.1, which is -suggested by Saris and colleagues (2009, p. 571).} - -\item{intcept}{The amount of standardized intercept (similar to Cohen's -\emph{d} that one would like to be detected (rejected). The default value -is 0.2, which is equivalent to a low effect size proposed by Cohen (1988, -1992).} - -\item{stdDelta}{The vector of the standardized parameters that one would -like to be detected (rejected). If this argument is specified, the value -here will overwrite the other arguments above. The order of the vector -must be the same as the row order from modification indices from the -\code{lavaan} object. If a single value is specified, the value will be -applied to all parameters.} - -\item{delta}{The vector of the unstandardized parameters that one would like -to be detected (rejected). If this argument is specified, the value here -will overwrite the other arguments above. The order of the vector must be -the same as the row order from modification indices from the \code{lavaan} -object. If a single value is specified, the value will be applied to all -parameters.} - -\item{cilevel}{The confidence level of the confidence interval of expected -parameter changes. The confidence intervals are used in the equivalence -testing.} - -\item{\dots}{arguments passed to \code{\link[lavaan:modificationIndices]{lavaan::modificationIndices()}}, -except for \code{delta}, which is already an argument (which can be -substituted for \code{stdDelta} or specific sets of parameters using -\code{stdLoad}, \code{cor}, \code{stdBeta}, and \code{intcept}).} -} -\value{ -A data frame with these variables: -\enumerate{ -\item \code{lhs}: The left-hand side variable, with respect to the operator in -in the lavaan \code{\link[lavaan:model.syntax]{lavaan::model.syntax()}} -\item \code{op}: The lavaan syntax operator: "~~" represents covariance, -"=~" represents factor loading, "~" represents regression, and -"~1" represents intercept. -\item \code{rhs}: The right-hand side variable -\item \code{group}: The level of the group variable for the parameter in question -\item \code{mi}: The modification index of the fixed parameter -\item \code{epc}: The EPC if the parameter is freely estimated -\item \code{target.epc}: The target EPC that represents the minimum size -of misspecification that one would like to be detected -by the test with a high power -\item \code{std.epc}: The standardized EPC if the parameter is freely estimated -\item \code{std.target.epc}: The standardized target expected parameter change -\item \code{significant.mi}: Represents whether the modification index value is -significant -\item \code{high.power}: Represents whether the power is enough to detect the -target expected parameter change -\item \code{decision.pow}: The decision whether the parameter is misspecified -or not based on Saris et al's method: \code{"M"} represents the parameter -is misspecified, \code{"NM"} represents the parameter is not misspecified, -\code{"EPC:M"} represents the parameter is misspecified decided by -checking the expected parameter change value, \code{"EPC:NM"} represents -the parameter is not misspecified decided by checking the expected -parameter change value, and \code{"I"} represents the decision is -inconclusive. -\item \code{se.epc}: The standard errors of the expected parameter changes. -\item \code{lower.epc}: The lower bound of the confidence interval of expected -parameter changes. -\item \code{upper.epc}: The upper bound of the confidence interval of expected -parameter changes. -\item \code{lower.std.epc}: Lower confidence limit of standardized EPCs -\item \code{upper.std.epc}: Upper confidence limit of standardized EPCs -\item \code{decision.ci}: Decision whether the parameter is misspecified -based on the CI method: \code{"M"} represents the -parameter is misspecified, \code{"NM"} represents the parameter is not -misspecified, and \code{"I"} represents the decision is inconclusive. -} - -The row numbers matches with the results obtained from the -\code{inspect(object, "mi")} function. -} -\description{ -The model fit evaluation approach using modification indices and expected -parameter changes. -} -\details{ -To decide whether a parameter should be freed, one can inspect its -modification index (MI) and expected parameter change (EPC). -Those values can be used to evaluate model fit by 2 methods. - -Method 1: Saris, Satorra, and van der Veld (2009, pp. 570--573) used -power (probability of detecting a significant MI) and EPC to decide whether -to free a parametr. First, one should evaluate whether a parameter's MI -is significant. Second, one should evaluate whether the power to detect a -target EPC is high enough. The combination of criteria leads to the -so-called "JRule" first implemented with LISREL (van der Veld et al., 2008): - -\itemize{ -\item If the MI is not significant and the power is low, -the test is inconclusive. -\item If the MI is not significant and the power is high, -there is no misspecification. -\item If the MI is significant and the power is low, -the fixed parameter is misspecified. -\item If the MI is significant and the power is high, -the EPC is investigated. If the EPC is large (greater than the -the target EPC), the parameter is misspecified. If the EPC is low -(lower than the target EPC), the parameter is not misspecificied. -} - -Method 2: The confidence interval (CI) of an EPC is calculated. -These CIs are compared with the range of trivial -misspecification, which could be (-\code{delta}, \code{delta}) or (0, -\code{delta}) for nonnegative parameters. - -\itemize{ -\item If a CI overlaps with the range of trivial misspecification, -the test is inconclusive. -\item If a CI completely exceeds the range of trivial misspecification, -the fixed parameters are severely misspecified. -\item If a CI is completely within the range of trivial misspecification, -the fixed parameters are trivially misspecified. -} -} -\examples{ - -library(lavaan) - -HS.model <- ' visual =~ x1 + x2 + x3 ' -fit <- cfa(HS.model, data = HolzingerSwineford1939, - group = "sex", group.equal = c("loadings","intercepts")) -miPowerFit(fit, free.remove = FALSE, op = "=~") # loadings -miPowerFit(fit, free.remove = FALSE, op = "~1") # intercepts - -model <- ' - # latent variable definitions - ind60 =~ x1 + x2 + x3 - dem60 =~ y1 + a*y2 + b*y3 + c*y4 - dem65 =~ y5 + a*y6 + b*y7 + c*y8 - - # regressions - dem60 ~ ind60 - dem65 ~ ind60 + dem60 - - # residual correlations - y1 ~~ y5 - y2 ~~ y4 + y6 - y3 ~~ y7 - y4 ~~ y8 - y6 ~~ y8 -' -fit2 <- sem(model, data = PoliticalDemocracy, meanstructure = TRUE) -miPowerFit(fit2, stdLoad = 0.3, cor = 0.2, stdBeta = 0.2, intcept = 0.5) - -} -\references{ -Cohen, J. (1988). \emph{Statistical power analysis for the -behavioral sciences} (2nd ed.). Hillsdale, NJ: Erlbaum. - -Cohen, J. (1992). A power primer. \emph{Psychological Bulletin, 112}(1), -155--159. \doi{10.1037/0033-2909.112.1.155} - -Saris, W. E., Satorra, A., & van der Veld, W. M. (2009). Testing structural -equation models or detection of misspecifications? \emph{Structural Equation -Modeling, 16}(4), 561--582. \doi{10.1080/10705510903203433} - -van der Veld, W. M., Saris, W. E., & Satorra, A. (2008). -\emph{JRule 3.0 Users Guide}. \doi{10.13140/RG.2.2.13609.90729} -} -\seealso{ -\code{\link[=moreFitIndices]{moreFitIndices()}} For the additional fit indices -information -} -\author{ -Sunthud Pornprasertmanit (\email{psunthud@gmail.com}) -} diff --git a/semTools/man/moreFitIndices.Rd b/semTools/man/moreFitIndices.Rd index 7e88d01..e2882bd 100644 --- a/semTools/man/moreFitIndices.Rd +++ b/semTools/man/moreFitIndices.Rd @@ -188,8 +188,8 @@ structural equation modeling} (pp. 209--231). New York, NY: Guilford. } \seealso{ \itemize{ -\item \code{\link[=miPowerFit]{miPowerFit()}} For the modification indices and their -power approach for model fit evaluation +\item \code{\link[=epcEquivFit]{epcEquivFit()}} For the equivalence testing based on expected +parameter changes for model fit evaluation \item \code{\link[=nullRMSEA]{nullRMSEA()}} For RMSEA of the default independence model } } diff --git a/semTools/man/semTools.Rd b/semTools/man/semTools.Rd index b3ecf61..188e0fe 100644 --- a/semTools/man/semTools.Rd +++ b/semTools/man/semTools.Rd @@ -15,7 +15,7 @@ functions a \emph{suite}. Our suites include: \code{\link[=moreFitIndices]{moreFitIndices()}}, \code{\link[=nullRMSEA]{nullRMSEA()}}, \code{\link[=singleParamTest]{singleParamTest()}}, -\code{\link[=miPowerFit]{miPowerFit()}}, and +\code{\link[=epcEquivFit]{epcEquivFit()}}, and \code{\link[=chisqSmallN]{chisqSmallN()}}} \item{Measurement Invariance: \code{\link[=measEq.syntax]{measEq.syntax()}},