Nothing
R/permuteMeasEq.R
In semTools: Useful Tools for Structural Equation Modeling
Defines functions
summ.mimic
summ.mgcfa
### Terrence D. Jorgensen
### Last updated: 12 March 2025
### permutation randomization test for measurement equivalence and DIF
## -----------------
## Class and Methods
## -----------------
##' Class for the Results of Permutation Randomization Tests of Measurement
##' Equivalence and DIF
##'
##' This class contains the results of tests of Measurement Equivalence and
##' Differential Item Functioning (DIF).
##'
##'
##' @name permuteMeasEq-class
##' @aliases permuteMeasEq-class show,permuteMeasEq-method
##' summary,permuteMeasEq-method hist,permuteMeasEq-method
##' @docType class
##'
##' @slot PT A `data.frame` returned by a call to
##' [lavaan::parTable()] on the constrained model
##' @slot modelType A character indicating the specified `modelType` in the
##' call to `permuteMeasEq`
##' @slot ANOVA A `numeric` vector indicating the results of the observed
##' (\eqn{\Delta})\eqn{\chi^2} test, based on the central \eqn{\chi^2}
##' distribution
##' @slot AFI.obs A vector of observed (changes in) user-selected fit measures
##' @slot AFI.dist The permutation distribution(s) of user-selected fit measures.
##' A `data.frame` with `n.Permutations` rows and one column for each
##' `AFI.obs`.
##' @slot AFI.pval A vector of *p* values (one for each element in slot
##' `AFI.obs`) calculated using slot `AFI.dist`, indicating the
##' probability of observing a change at least as extreme as `AFI.obs`
##' if the null hypothesis were true
##' @slot MI.obs A `data.frame` of observed Lagrange Multipliers
##' (modification indices) associated with the equality constraints or fixed
##' parameters specified in the `param` argument. This is a subset of the
##' output returned by a call to [lavaan::lavTestScore()] on the
##' constrained model.
##' @slot MI.dist The permutation distribution of the maximum modification index
##' (among those seen in slot `MI.obs$X2`) at each permutation of group
##' assignment or of `covariates`
##' @slot extra.obs If `permuteMeasEq` was called with an `extra`
##' function, the output when applied to the original data is concatenated
##' into this vector
##' @slot extra.dist A `data.frame`, each column of which contains the
##' permutation distribution of the corresponding statistic in slot
##' `extra.obs`
##' @slot n.Permutations An `integer` indicating the number of permutations
##' requested by the user
##' @slot n.Converged An `integer` indicating the number of permuation
##' iterations which yielded a converged solution
##' @slot n.nonConverged An `integer` vector of length
##' `n.Permutations` indicating how many times group assignment was
##' randomly permuted (at each iteration) before converging on a solution
##' @slot n.Sparse Only relevant with `ordered` indicators when
##' `modelType == "mgcfa"`. An `integer` vector of length
##' `n.Permutations` indicating how many times group assignment was
##' randomly permuted (at each iteration) before obtaining a sample with all
##' categories observed in all groups.
##' @slot oldSeed An `integer` vector storing the value of
##' `.Random.seed` before running `permuteMeasEq`. Only relevant
##' when using a parallel/multicore option and the original
##' `RNGkind() != "L'Ecuyer-CMRG"`. This enables users to restore their
##' previous `.Random.seed` state, if desired, by running:
##' `.Random.seed[-1] <- permutedResults@oldSeed[-1]`
##' @section Objects from the Class: Objects can be created via the
##' [semTools::permuteMeasEq()] function.
##'
##' @return
##' \itemize{
##' \item The `show` method prints a summary of the multiparameter
##' omnibus test results, using the user-specified AFIs. The parametric
##' (\eqn{\Delta})\eqn{\chi^2} test is also displayed.
##' \item The `summary` method prints the same information from the
##' `show` method, but when `extra = FALSE` (the default) it also
##' provides a table summarizing any requested follow-up tests of DIF using
##' modification indices in slot `MI.obs`. The user can also specify an
##' `alpha` level for flagging modification indices as significant, as
##' well as `nd` (the number of digits displayed). For each modification
##' index, the *p* value is displayed using a central \eqn{\chi^2}
##' distribution with the *df* shown in that column. Additionally, a
##' *p* value is displayed using the permutation distribution of the
##' maximum index, which controls the familywise Type I error rate in a manner
##' similar to Tukey's studentized range test. If any indices are flagged as
##' significant using the `tukey.p.value`, then a message is displayed for
##' each flagged index. The invisibly returned `data.frame` is the
##' displayed table of modification indices, unless
##' [semTools::permuteMeasEq()] was called with `param = NULL`,
##' in which case the invisibly returned object is `object`. If
##' `extra = TRUE`, the permutation-based *p* values for each
##' statistic returned by the `extra` function are displayed and returned
##' in a `data.frame` instead of the modification indices requested in the
##' `param` argument.
##' \item The `hist` method returns a list of `length == 2`,
##' containing the arguments for the call to `hist` and the arguments
##' to the call for `legend`, respectively. This list may facilitate
##' creating a customized histogram of `AFI.dist`, `MI.dist`, or
##' `extra.dist`
##' }
##'
##' @author Terrence D. Jorgensen (University of Amsterdam;
##' \email{TJorgensen314@@gmail.com})
##'
##' @seealso [semTools::permuteMeasEq()]
##'
##' @examples
##'
##' # See the example from the permuteMeasEq function
##'
setClass("permuteMeasEq", slots = c(PT = "data.frame",
modelType = "character",
ANOVA = "vector",
AFI.obs = "vector",
AFI.dist = "data.frame",
AFI.pval = "vector",
MI.obs = "data.frame",
MI.dist = "vector",
extra.obs = "vector",
extra.dist = "data.frame",
n.Permutations = "integer",
n.Converged = "integer",
n.nonConverged = "vector",
n.Sparse = "vector",
oldSeed = "integer"))
##' @rdname permuteMeasEq-class
##' @aliases show,permuteMeasEq-method
##' @export
setMethod("show", "permuteMeasEq", function(object) {
## print warning if there are nonConverged permutations
if (object@n.Permutations != object@n.Converged) {
warning(paste("Only", object@n.Converged, "out of",
object@n.Permutations, "models converged within",
max(object@n.nonConverged), "attempts per permutation.\n\n"))
}
## print ANOVA
cat("Omnibus p value based on parametric chi-squared difference test:\n\n")
print(round(object@ANOVA, digits = 3))
## print permutation results
cat("\n\nOmnibus p values based on nonparametric permutation method: \n\n")
AFI <- data.frame(AFI.Difference = object@AFI.obs, p.value = object@AFI.pval)
class(AFI) <- c("lavaan.data.frame","data.frame")
print(AFI, nd = 3)
invisible(object)
})
##' @rdname permuteMeasEq-class
##' @aliases summary,permuteMeasEq-method
##' @export
setMethod("summary", "permuteMeasEq", function(object, alpha = .05, nd = 3,
extra = FALSE) {
## print warning if there are nonConverged permutations
if (object@n.Permutations != object@n.Converged) {
warning(paste("Only", object@n.Converged, "out of",
object@n.Permutations, "models converged within",
max(object@n.nonConverged), "attempts per permutation.\n\n"))
}
## print ANOVA
cat("Omnibus p value based on parametric chi-squared difference test:\n\n")
print(round(object@ANOVA, digits = nd))
## print permutation results
cat("\n\nOmnibus p values based on nonparametric permutation method: \n\n")
AFI <- data.frame(AFI.Difference = object@AFI.obs, p.value = object@AFI.pval)
class(AFI) <- c("lavaan.data.frame","data.frame")
print(AFI, nd = nd)
## print extras or DIF test results, if any were requested
if (extra && length(object@extra.obs)) {
cat("\n\nUnadjusted p values of extra statistics,\n",
"based on permutation distribution of each statistic: \n\n")
MI <- data.frame(Statistic = object@extra.obs)
class(MI) <- c("lavaan.data.frame","data.frame")
MI$p.value <- sapply(names(object@extra.dist), function(nn) {
mean(abs(object@extra.dist[,nn]) >= abs(object@extra.obs[nn]), na.rm = TRUE)
})
MI$flag <- ifelse(MI$p.value < alpha, "* ", "")
print(MI, nd = nd)
} else if (length(object@MI.dist)) {
cat("\n\n Modification indices for equality constrained parameter estimates,\n",
"with unadjusted 'p.value' based on chi-squared distribution and\n",
"adjusted 'tukey.p.value' based on permutation distribution of the\n",
"maximum modification index per iteration: \n\n")
MI <- do.call(paste("summ", object@modelType, sep = "."),
args = list(object = object, alpha = alpha))
print(MI, nd = nd)
## print messages about potential DIF
if (all(MI$tukey.p.value > alpha)) {
cat("\n\n No equality constraints were flagged as significant.\n\n")
return(invisible(MI))
}
if (object@modelType == "mgcfa") {
cat("\n\nThe following equality constraints were flagged as significant:\n\n")
for (i in which(MI$tukey.p.value < alpha)) {
cat("Parameter '", MI$parameter[i], "' may differ between Groups '",
MI$group.lhs[i], "' and '", MI$group.rhs[i], "'.\n", sep = "")
}
cat("\nUse lavTestScore(..., epc = TRUE) on your constrained model to",
"display expected parameter changes for these equality constraints\n\n")
}
} else return(invisible(object))
invisible(MI)
})
summ.mgcfa <- function(object, alpha) {
MI <- object@MI.obs
class(MI) <- c("lavaan.data.frame","data.frame")
PT <- object@PT
eqPar <- rbind(PT[PT$plabel %in% MI$lhs, ], PT[PT$plabel %in% MI$rhs, ])
MI$flag <- ""
MI$parameter <- ""
MI$group.lhs <- ""
MI$group.rhs <- ""
for (i in 1:nrow(MI)) {
par1 <- eqPar$par[ eqPar$plabel == MI$lhs[i] ]
par2 <- eqPar$par[ eqPar$plabel == MI$rhs[i] ]
MI$parameter[i] <- par1
MI$group.lhs[i] <- eqPar$group.label[ eqPar$plabel == MI$lhs[i] ]
MI$group.rhs[i] <- eqPar$group.label[ eqPar$plabel == MI$rhs[i] ]
if (par1 != par2) {
myMessage <- paste0("Constraint '", MI$lhs[i], "==", MI$rhs[i],
"' refers to different parameters: \n'",
MI$lhs[i], "' is '", par1, "' in group '",
MI$group.lhs[i], "'\n'",
MI$rhs[i], "' is '", par2, "' in group '",
MI$group.rhs[i], "'\n")
warning(myMessage)
}
if (MI$tukey.p.value[i] < alpha) MI$flag[i] <- "* -->"
}
MI
}
summ.mimic <- function(object, alpha) {
MI <- object@MI.obs
class(MI) <- c("lavaan.data.frame","data.frame")
MI$flag <- ifelse(MI$tukey.p.value < alpha, "* ", "")
MI
}
##' @rdname permuteMeasEq-class
##' @aliases hist,permuteMeasEq-method
##' @importFrom stats qchisq dchisq quantile
##' @param object,x object of class `permuteMeasEq`
##' @param ... Additional arguments to pass to [graphics::hist()]
##' @param AFI `character` indicating the fit measure whose permutation
##' distribution should be plotted
##' @param alpha alpha level used to draw confidence limits in `hist` and
##' flag significant statistics in `summary` output
##' @param nd number of digits to display
##' @param extra `logical` indicating whether the `summary` output
##' should return permutation-based *p* values for each statistic returned
##' by the `extra` function. If `FALSE` (default), `summary`
##' will return permutation-based *p* values for each modification index.
##' @param printLegend `logical`. If `TRUE` (default), a legend will
##' be printed with the histogram
##' @param legendArgs `list` of arguments passed to the
##' [graphics::legend()] function. The default argument is a list
##' placing the legend at the top-left of the figure.
##' @export
setMethod("hist", "permuteMeasEq", function(x, ..., AFI, alpha = .05, nd = 3,
printLegend = TRUE,
legendArgs = list(x = "topleft")) {
histArgs <- list(...)
histArgs$x <- x@AFI.dist[[AFI]]
if (is.null(histArgs$col)) histArgs$col <- "grey69"
histArgs$freq <- !grepl("chi", AFI)
histArgs$ylab <- if (histArgs$freq) "Frequency" else "Probability Density"
if (printLegend) {
if (is.null(legendArgs$box.lty)) legendArgs$box.lty <- 0
if (nd < length(strsplit(as.character(1 / alpha), "")[[1]]) - 1) {
warning(paste0("The number of digits argument (nd = ", nd ,
") is too low to display your p value at the",
" same precision as your requested alpha level (alpha = ",
alpha, ")"))
}
if (x@AFI.pval[[AFI]] < (1 / 10^nd)) {
pVal <- paste(c("< .", rep(0, nd - 1),"1"), collapse = "")
} else {
pVal <- paste("=", round(x@AFI.pval[[AFI]], nd))
}
}
delta <- length(x@MI.dist) > 0L && x@modelType == "mgcfa"
if (grepl("chi", AFI)) { ####################################### Chi-squared
ChiSq <- x@AFI.obs[AFI]
DF <- x@ANOVA[2]
histArgs$xlim <- range(c(ChiSq, x@AFI.dist[[AFI]], qchisq(c(.01, .99), DF)))
xVals <- seq(histArgs$xlim[1], histArgs$xlim[2], by = .1)
theoDist <- dchisq(xVals, df = DF)
TheoCrit <- round(qchisq(p = alpha, df = DF, lower.tail = FALSE), 2)
Crit <- quantile(histArgs$x, probs = 1 - alpha)
if (ChiSq > histArgs$xlim[2]) histArgs$xlim[2] <- ChiSq
if (delta) {
histArgs$main <- expression(Permutation~Distribution~of~Delta*chi^2)
histArgs$xlab <- expression(Delta*chi^2)
if (printLegend) {
legendArgs$legend <- c(bquote(Theoretical~Delta*chi[Delta*.(paste("df =", DF))]^2 ~ Distribution),
bquote(Critical~chi[alpha~.(paste(" =", alpha))]^2 == .(round(TheoCrit, nd))),
bquote(.(paste("Permuted Critical Value =", round(Crit, nd)))),
bquote(Observed~Delta*chi^2 == .(round(ChiSq, nd))),
expression(paste("")),
bquote(Permuted~italic(p)~.(pVal)))
}
} else {
histArgs$main <- expression(Permutation~Distribution~of~chi^2)
histArgs$xlab <- expression(chi^2)
if (printLegend) {
legendArgs$legend <- c(bquote(Theoretical~chi[.(paste("df =", DF))]^2 ~ Distribution),
bquote(Critical~chi[alpha~.(paste(" =", alpha))]^2 == .(round(TheoCrit, nd))),
bquote(.(paste("Permuted Critical Value =", round(Crit, nd)))),
bquote(Observed~chi^2 == .(round(ChiSq, nd))),
expression(paste("")),
bquote(Permuted~italic(p)~.(pVal)))
}
}
H <- do.call(hist, c(histArgs["x"], plot = FALSE))
histArgs$ylim <- c(0, max(H$density, theoDist))
if (printLegend) {
legendArgs <- c(legendArgs, list(lty = c(2, 2, 1, 1, 0, 0),
lwd = c(2, 2, 2, 3, 0, 0),
col = c("black","black","black","red","","")))
}
} else { ################################################### other AFIs
badness <- grepl(pattern = "fmin|aic|bic|rmr|rmsea|cn|sic|hqc",
x = AFI, ignore.case = TRUE)
if (badness) {
Crit <- quantile(histArgs$x, probs = 1 - alpha)
} else {
Crit <- quantile(histArgs$x, probs = alpha)
}
histArgs$xlim <- range(histArgs$x, x@AFI.obs[AFI])
if (delta) {
histArgs$main <- bquote(~Permutation~Distribution~of~Delta*.(toupper(AFI)))
histArgs$xlab <- bquote(~Delta*.(toupper(AFI)))
if (printLegend) {
legendArgs$legend <- c(bquote(Critical~Delta*.(toupper(AFI))[alpha~.(paste(" =", alpha))] == .(round(Crit, nd))),
bquote(Observed~Delta*.(toupper(AFI)) == .(round(x@AFI.obs[AFI], nd))),
expression(paste("")),
bquote(Permuted~italic(p)~.(pVal)))
}
} else {
histArgs$main <- paste("Permutation Distribution of", toupper(AFI))
histArgs$xlab <- toupper(AFI)
if (printLegend) {
legendArgs$legend <- c(bquote(Critical~.(toupper(AFI))[alpha~.(paste(" =", alpha))] == .(round(Crit, nd))),
bquote(Observed~.(toupper(AFI)) == .(round(x@AFI.obs[AFI], nd))),
expression(paste("")),
bquote(Permuted~italic(p)~.(pVal)))
}
}
if (printLegend) {
legendArgs <- c(legendArgs, list(lty = c(1, 1, 0, 0),
lwd = c(2, 3, 0, 0),
col = c("black","red","","")))
}
}
## print histogram (and optionally, print legend)
suppressWarnings({
do.call(hist, histArgs)
if (grepl("chi", AFI)) {
lines(x = xVals, y = theoDist, lwd = 2, lty = 2)
abline(v = TheoCrit, col = "black", lwd = 2, lty = 2)
}
abline(v = Crit, col = "black", lwd = 2)
abline(v = x@AFI.obs[AFI], col = "red", lwd = 3)
if (printLegend) do.call(legend, legendArgs)
})
## return arguments to create histogram (and optionally, legend)
invisible(list(hist = histArgs, legend = legendArgs))
})
## --------------------
## Constructor Function
## --------------------
##' Permutation Randomization Tests of Measurement Equivalence and Differential
##' Item Functioning (DIF)
##'
##' The function `permuteMeasEq` provides tests of hypotheses involving
##' measurement equivalence, in one of two frameworks: multigroup CFA or MIMIC
##' models.
##'
##'
##' The function `permuteMeasEq` provides tests of hypotheses involving
##' measurement equivalence, in one of two frameworks:
##' \enumerate{
##' \item{1} For multiple-group CFA models, provide a pair of nested lavaan objects,
##' the less constrained of which (`uncon`) freely estimates a set of
##' measurement parameters (e.g., factor loadings, intercepts, or thresholds;
##' specified in `param`) in all groups, and the more constrained of which
##' (`con`) constrains those measurement parameters to equality across
##' groups. Group assignment is repeatedly permuted and the models are fit to
##' each permutation, in order to produce an empirical distribution under the
##' null hypothesis of no group differences, both for (a) changes in
##' user-specified fit measures (see `AFIs` and `moreAFIs`) and for
##' (b) the maximum modification index among the user-specified equality
##' constraints. Configural invariance can also be tested by providing that
##' fitted lavaan object to `con` and leaving `uncon = NULL`, in which
##' case `param` must be `NULL` as well.
##'
##' \item{2} In MIMIC models, one or a set of continuous and/or discrete
##' `covariates` can be permuted, and a constrained model is fit to each
##' permutation in order to provide a distribution of any fit measures (namely,
##' the maximum modification index among fixed parameters in `param`) under
##' the null hypothesis of measurement equivalence across levels of those
##' covariates.
##' }
##'
##' In either framework, modification indices for equality constraints or fixed
##' parameters specified in `param` are calculated from the constrained
##' model (`con`) using the function [lavaan::lavTestScore()].
##'
##' For multiple-group CFA models, the multiparameter omnibus null hypothesis of
##' measurement equivalence/invariance is that there are no group differences in
##' any measurement parameters (of a particular type). This can be tested using
##' the `anova` method on nested `lavaan` objects, as seen in the
##' output of [semTools::measurementInvariance()], or by inspecting
##' the change in alternative fit indices (AFIs) such as the CFI. The
##' permutation randomization method employed by `permuteMeasEq` generates
##' an empirical distribution of any `AFIs` under the null hypothesis, so
##' the user is not restricted to using fixed cutoffs proposed by Cheung &
##' Rensvold (2002), Chen (2007), or Meade, Johnson, & Braddy (2008).
##'
##' If the multiparameter omnibus null hypothesis is rejected, partial
##' invariance can still be established by freeing invalid equality constraints,
##' as long as equality constraints are valid for at least two indicators per
##' factor. Modification indices can be calculated from the constrained model
##' (`con`), but multiple testing leads to inflation of Type I error rates.
##' The permutation randomization method employed by `permuteMeasEq`
##' creates a distribution of the maximum modification index if the null
##' hypothesis is true, which allows the user to control the familywise Type I
##' error rate in a manner similar to Tukey's *q* (studentized range)
##' distribution for the Honestly Significant Difference (HSD) post hoc test.
##'
##' For MIMIC models, DIF can be tested by comparing modification indices of
##' regression paths to the permutation distribution of the maximum modification
##' index, which controls the familywise Type I error rate. The MIMIC approach
##' could also be applied with multiple-group models, but the grouping variable
##' would not be permuted; rather, the covariates would be permuted separately
##' within each group to preserve between-group differences. So whether
##' parameters are constrained or unconstrained across groups, the MIMIC
##' approach is only for testing null hypotheses about the effects of
##' `covariates` on indicators, controlling for common factors.
##'
##' In either framework, [lavaan::lavaan()]'s `group.label`
##' argument is used to preserve the order of groups seen in `con` when
##' permuting the data.
##'
##'
##' @importFrom lavaan lavInspect parTable
##'
##' @param nPermute An integer indicating the number of random permutations used
##' to form empirical distributions under the null hypothesis.
##' @param modelType A character string indicating type of model employed:
##' multiple-group CFA (`"mgcfa"`) or MIMIC (`"mimic"`).
##' @param con The constrained `lavaan` object, in which the parameters
##' specified in `param` are constrained to equality across all groups when
##' `modelType = "mgcfa"`, or which regression paths are fixed to zero when
##' `modelType = "mimic"`. In the case of testing *configural*
##' invariance when `modelType = "mgcfa"`, `con` is the configural
##' model (implicitly, the unconstrained model is the saturated model, so use
##' the defaults `uncon = NULL` and `param = NULL`). When
##' `modelType = "mimic"`, `con` is the MIMIC model in which the
##' covariate predicts the latent construct(s) but no indicators (unless they
##' have already been identified as DIF items).
##' @param uncon Optional. The unconstrained `lavaan` object, in which the
##' parameters specified in `param` are freely estimated in all groups.
##' When `modelType = "mgcfa"`, only in the case of testing
##' *configural* invariance should `uncon = NULL`. When
##' `modelType = "mimic"`, any non-`NULL uncon` is silently set to
##' `NULL`.
##' @param null Optional. A `lavaan` object, in which an alternative null
##' model is fit (besides the default independence model specified by
##' `lavaan`) for the calculation of incremental fit indices. See Widamin &
##' Thompson (2003) for details. If `NULL`, `lavaan`'s default
##' independence model is used.
##' @param param An optional character vector or list of character vectors
##' indicating which parameters the user would test for DIF following a
##' rejection of the omnibus null hypothesis tested using
##' (`more`)`AFIs`. Note that `param` does not guarantee certain
##' parameters *are* constrained in `con`; that is for the user to
##' specify when fitting the model. If users have any "anchor items" that they
##' would never intend to free across groups (or levels of a covariate), these
##' should be excluded from `param`; exceptions to a type of parameter can
##' be specified in `freeParam`. When `modelType = "mgcfa"`,
##' `param` indicates which parameters of interest are constrained across
##' groups in `con` and are unconstrained in `uncon`. Parameter names
##' must match those returned by `names(coef(con))`, but omitting any
##' group-specific suffixes (e.g., `"f1~1"` rather than `"f1~1.g2"`)
##' or user-specified labels (that is, the parameter names must follow the rules
##' of lavaan's [lavaan::model.syntax()]). Alternatively (or
##' additionally), to test all constraints of a certain type (or multiple types)
##' of parameter in `con`, `param` may take any combination of the
##' following values: `"loadings"`, `"intercepts"`,
##' `"thresholds"`, `"residuals"`, `"residual.covariances"`,
##' `"means"`, `"lv.variances"`, and/or `"lv.covariances"`. When
##' `modelType = "mimic"`, `param` must be a vector of individual
##' parameters or a list of character strings to be passed one-at-a-time to
##' `lavaan::lavTestScore(object = con, add = param[i])`,
##' indicating which (sets of) regression paths fixed to zero in `con` that
##' the user would consider freeing (i.e., exclude anchor items). If
##' `modelType = "mimic"` and `param` is a list of character strings,
##' the multivariate test statistic will be saved for each list element instead
##' of 1-*df* modification indices for each individual parameter, and
##' `names(param)` will name the rows of the `MI.obs` slot (see
##' [permuteMeasEq-class]). Set `param = NULL` (default) to avoid
##' collecting modification indices for any follow-up tests.
##' @param freeParam An optional character vector, silently ignored when
##' `modelType = "mimic"`. If `param` includes a type of parameter
##' (e.g., `"loadings"`), `freeParam` indicates exceptions (i.e.,
##' anchor items) that the user would *not* intend to free across groups
##' and should therefore be ignored when calculating *p* values adjusted
##' for the number of follow-up tests. Parameter types that are already
##' unconstrained across groups in the fitted `con` model (i.e., a
##' *partial* invariance model) will automatically be ignored, so they do
##' not need to be specified in `freeParam`. Parameter names must match
##' those returned by `names(coef(con))`, but omitting any group-specific
##' suffixes (e.g., `"f1~1"` rather than `"f1~1.g2"`) or
##' user-specified labels (that is, the parameter names must follow the rules of
##' lavaan [lavaan::model.syntax()]).
##' @param covariates An optional character vector, only applicable when
##' `modelType = "mimic"`. The observed data are partitioned into columns
##' indicated by `covariates`, and the rows are permuted simultaneously for
##' the entire set before being merged with the remaining data. Thus, the
##' covariance structure is preserved among the covariates, which is necessary
##' when (e.g.) multiple dummy codes are used to represent a discrete covariate
##' or when covariates interact. If `covariates = NULL` when
##' `modelType = "mimic"`, the value of `covariates` is inferred by
##' searching `param` for predictors (i.e., variables appearing after the
##' "`~`" operator).
##' @param AFIs A character vector indicating which alternative fit indices (or
##' chi-squared itself) are to be used to test the multiparameter omnibus null
##' hypothesis that the constraints specified in `con` hold in the
##' population. Any fit measures returned by [lavaan::fitMeasures()]
##' may be specified (including constants like `"df"`, which would be
##' nonsensical). If both `AFIs` and `moreAFIs` are `NULL`, only
##' `"chisq"` will be returned.
##' @param moreAFIs Optional. A character vector indicating which (if any)
##' alternative fit indices returned by [semTools::moreFitIndices()]
##' are to be used to test the multiparameter omnibus null hypothesis that the
##' constraints specified in `con` hold in the population.
##' @param maxSparse Only applicable when `modelType = "mgcfa"` and at
##' least one indicator is `ordered`. An integer indicating the maximum
##' number of consecutive times that randomly permuted group assignment can
##' yield a sample in which at least one category (of an `ordered`
##' indicator) is unobserved in at least one group, such that the same set of
##' parameters cannot be estimated in each group. If such a sample occurs, group
##' assignment is randomly permuted again, repeatedly until a sample is obtained
##' with all categories observed in all groups. If `maxSparse` is exceeded,
##' `NA` will be returned for that iteration of the permutation
##' distribution.
##' @param maxNonconv An integer indicating the maximum number of consecutive
##' times that a random permutation can yield a sample for which the model does
##' not converge on a solution. If such a sample occurs, permutation is
##' attempted repeatedly until a sample is obtained for which the model does
##' converge. If `maxNonconv` is exceeded, `NA` will be returned for
##' that iteration of the permutation distribution, and a warning will be
##' printed when using `show` or `summary`.
##' @param showProgress Logical. Indicating whether to display a progress bar
##' while permuting. Silently set to `FALSE` when using parallel options.
##' @param warn Sets the handling of warning messages when fitting model(s) to
##' permuted data sets. See [base::options()].
##' @param datafun An optional function that can be applied to the data
##' (extracted from `con`) after each permutation, but before fitting the
##' model(s) to each permutation. The `datafun` function must have an
##' argument named `data` that accepts a `data.frame`, and it must
##' return a `data.frame` containing the same column names. The column
##' order may differ, the values of those columns may differ (so be careful!),
##' and any additional columns will be ignored when fitting the model, but an
##' error will result if any column names required by the model syntax do not
##' appear in the transformed data set. Although available for any
##' `modelType`, `datafun` may be useful when using the MIMIC method
##' to test for nonuniform DIF (metric/weak invariance) by using product
##' indicators for a latent factor representing the interaction between a factor
##' and one of the `covariates`, in which case the product indicators would
##' need to be recalculated after each permutation of the `covariates`. To
##' access other R objects used within `permuteMeasEq`, the arguments to
##' `datafun` may also contain any subset of the following: `"con"`,
##' `"uncon"`, `"null"`, `"param"`, `"freeParam"`,
##' `"covariates"`, `"AFIs"`, `"moreAFIs"`, `"maxSparse"`,
##' `"maxNonconv"`, and/or `"iseed"`. The values for those arguments
##' will be the same as the values supplied to `permuteMeasEq`.
##' @param extra An optional function that can be applied to any (or all) of the
##' fitted lavaan objects (`con`, `uncon`, and/or `null`). This
##' function will also be applied after fitting the model(s) to each permuted
##' data set. To access the R objects used within `permuteMeasEq`, the
##' arguments to `extra` must be any subset of the following: `"con"`,
##' `"uncon"`, `"null"`, `"param"`, `"freeParam"`,
##' `"covariates"`, `"AFIs"`, `"moreAFIs"`, `"maxSparse"`,
##' `"maxNonconv"`, and/or `"iseed"`. The values for those arguments
##' will be the same as the values supplied to `permuteMeasEq`. The
##' `extra` function must return a named `numeric` vector or a named
##' `list` of scalars (i.e., a `list` of `numeric` vectors of
##' `length == 1`). Any unnamed elements (e.g., `""` or `NULL`)
##' of the returned object will result in an error.
##' @param parallelType The type of parallel operation to be used (if any). The
##' default is `"none"`. Forking is not possible on Windows, so if
##' `"multicore"` is requested on a Windows machine, the request will be
##' changed to `"snow"` with a message.
##' @param ncpus Integer: number of processes to be used in parallel operation.
##' If `NULL` (the default) and `parallelType %in%
##' c("multicore","snow")`, the default is one less than the maximum number of
##' processors detected by [parallel::detectCores()]. This default is
##' also silently set if the user specifies more than the number of processors
##' detected.
##' @param cl An optional \pkg{parallel} or \pkg{snow} cluster for use when
##' `parallelType = "snow"`. If `NULL`, a `"PSOCK"` cluster on
##' the local machine is created for the duration of the `permuteMeasEq`
##' call. If a valid [parallel::makeCluster()] object is supplied,
##' `parallelType` is silently set to `"snow"`, and `ncpus` is
##' silently set to `length(cl)`.
##' @param iseed Integer: Only used to set the states of the RNG when using
##' parallel options, in which case [base::RNGkind()] is set to
##' `"L'Ecuyer-CMRG"` with a message. See
##' [parallel::clusterSetRNGStream()] and Section 6 of
##' `vignette("parallel", "parallel")` for more details. If user supplies
##' an invalid value, `iseed` is silently set to the default (12345). To
##' set the state of the RNG when not using parallel options, call
##' [base::set.seed()] before calling `permuteMeasEq`.
##'
##' @return The [permuteMeasEq-class] object representing the results of
##' testing measurement equivalence (the multiparameter omnibus test) and DIF
##' (modification indices), as well as diagnostics and any `extra` output.
##'
##' @author Terrence D. Jorgensen (University of Amsterdam;
##' \email{TJorgensen314@@gmail.com})
##'
##' @seealso [stats::TukeyHSD()], [lavaan::lavTestScore()],
##' [semTools::measurementInvariance()],
##' [semTools::measurementInvarianceCat()]
##'
##' @references
##'
##' **Papers about permutation tests of measurement equivalence:**
##'
##' Jorgensen, T. D., Kite, B. A., Chen, P.-Y., & Short, S. D. (2018).
##' Permutation randomization methods for testing measurement equivalence and
##' detecting differential item functioning in multiple-group confirmatory
##' factor analysis. *Psychological Methods, 23*(4), 708--728.
##' \doi{10.1037/met0000152}
##'
##' Kite, B. A., Jorgensen, T. D., & Chen, P.-Y. (2018). Random permutation
##' testing applied to measurement invariance testing with ordered-categorical
##' indicators. *Structural Equation Modeling 25*(4), 573--587.
##' \doi{10.1080/10705511.2017.1421467}
##'
##' Jorgensen, T. D. (2017). Applying permutation tests and multivariate
##' modification indices to configurally invariant models that need
##' respecification. *Frontiers in Psychology, 8*(1455).
##' \doi{10.3389/fpsyg.2017.01455}
##'
##' **Additional reading:**
##'
##' Chen, F. F. (2007). Sensitivity of goodness of fit indexes to
##' lack of measurement invariance. *Structural Equation Modeling, 14*(3),
##' 464--504. \doi{10.1080/10705510701301834}
##'
##' Cheung, G. W., & Rensvold, R. B. (2002). Evaluating goodness-of-fit indexes
##' for testing measurement invariance. *Structural Equation Modeling,
##' 9*(2), 233--255. \doi{10.1207/S15328007SEM0902_5}
##'
##' Meade, A. W., Johnson, E. C., & Braddy, P. W. (2008). Power and sensitivity
##' of alternative fit indices in tests of measurement invariance. *Journal
##' of Applied Psychology, 93*(3), 568--592. \doi{10.1037/0021-9010.93.3.568}
##'
##' Widamin, K. F., & Thompson, J. S. (2003). On specifying the null model for
##' incremental fit indices in structural equation modeling. *Psychological
##' Methods, 8*(1), 16--37. \doi{10.1037/1082-989X.8.1.16}
##'
##' @examples
##'
##' \donttest{
##'
##' ########################
##' ## Multiple-Group CFA ##
##' ########################
##'
##' ## create 3-group data in lavaan example(cfa) data
##' HS <- lavaan::HolzingerSwineford1939
##' HS$ageGroup <- ifelse(HS$ageyr < 13, "preteen",
##' ifelse(HS$ageyr > 13, "teen", "thirteen"))
##'
##' ## specify and fit an appropriate null model for incremental fit indices
##' mod.null <- c(paste0("x", 1:9, " ~ c(T", 1:9, ", T", 1:9, ", T", 1:9, ")*1"),
##' paste0("x", 1:9, " ~~ c(L", 1:9, ", L", 1:9, ", L", 1:9, ")*x", 1:9))
##' fit.null <- cfa(mod.null, data = HS, group = "ageGroup")
##'
##' ## fit target model with varying levels of measurement equivalence
##' mod.config <- '
##' visual =~ x1 + x2 + x3
##' textual =~ x4 + x5 + x6
##' speed =~ x7 + x8 + x9
##' '
##' fit.config <- cfa(mod.config, data = HS, std.lv = TRUE, group = "ageGroup")
##' fit.metric <- cfa(mod.config, data = HS, std.lv = TRUE, group = "ageGroup",
##' group.equal = "loadings")
##' fit.scalar <- cfa(mod.config, data = HS, std.lv = TRUE, group = "ageGroup",
##' group.equal = c("loadings","intercepts"))
##'
##'
##' ####################### Permutation Method
##'
##' ## fit indices of interest for multiparameter omnibus test
##' myAFIs <- c("chisq","cfi","rmsea","mfi","aic")
##' moreAFIs <- c("gammaHat","adjGammaHat")
##'
##' ## Use only 20 permutations for a demo. In practice,
##' ## use > 1000 to reduce sampling variability of estimated p values
##'
##' ## test configural invariance
##' set.seed(12345)
##' out.config <- permuteMeasEq(nPermute = 20, con = fit.config)
##' out.config
##'
##' ## test metric equivalence
##' set.seed(12345) # same permutations
##' out.metric <- permuteMeasEq(nPermute = 20, uncon = fit.config, con = fit.metric,
##' param = "loadings", AFIs = myAFIs,
##' moreAFIs = moreAFIs, null = fit.null)
##' summary(out.metric, nd = 4)
##'
##' ## test scalar equivalence
##' set.seed(12345) # same permutations
##' out.scalar <- permuteMeasEq(nPermute = 20, uncon = fit.metric, con = fit.scalar,
##' param = "intercepts", AFIs = myAFIs,
##' moreAFIs = moreAFIs, null = fit.null)
##' summary(out.scalar)
##'
##' ## Not much to see without significant DIF.
##' ## Try using an absurdly high alpha level for illustration.
##' outsum <- summary(out.scalar, alpha = .50)
##'
##' ## notice that the returned object is the table of DIF tests
##' outsum
##'
##' ## visualize permutation distribution
##' hist(out.config, AFI = "chisq")
##' hist(out.metric, AFI = "chisq", nd = 2, alpha = .01,
##' legendArgs = list(x = "topright"))
##' hist(out.scalar, AFI = "cfi", printLegend = FALSE)
##'
##'
##' ####################### Extra Output
##'
##' ## function to calculate expected change of Group-2 and -3 latent means if
##' ## each intercept constraint were released
##' extra <- function(con) {
##' output <- list()
##' output["x1.vis2"] <- lavTestScore(con, release = 19:20, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[70]
##' output["x1.vis3"] <- lavTestScore(con, release = 19:20, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[106]
##' output["x2.vis2"] <- lavTestScore(con, release = 21:22, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[70]
##' output["x2.vis3"] <- lavTestScore(con, release = 21:22, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[106]
##' output["x3.vis2"] <- lavTestScore(con, release = 23:24, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[70]
##' output["x3.vis3"] <- lavTestScore(con, release = 23:24, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[106]
##' output["x4.txt2"] <- lavTestScore(con, release = 25:26, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[71]
##' output["x4.txt3"] <- lavTestScore(con, release = 25:26, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[107]
##' output["x5.txt2"] <- lavTestScore(con, release = 27:28, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[71]
##' output["x5.txt3"] <- lavTestScore(con, release = 27:28, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[107]
##' output["x6.txt2"] <- lavTestScore(con, release = 29:30, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[71]
##' output["x6.txt3"] <- lavTestScore(con, release = 29:30, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[107]
##' output["x7.spd2"] <- lavTestScore(con, release = 31:32, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[72]
##' output["x7.spd3"] <- lavTestScore(con, release = 31:32, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[108]
##' output["x8.spd2"] <- lavTestScore(con, release = 33:34, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[72]
##' output["x8.spd3"] <- lavTestScore(con, release = 33:34, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[108]
##' output["x9.spd2"] <- lavTestScore(con, release = 35:36, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[72]
##' output["x9.spd3"] <- lavTestScore(con, release = 35:36, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[108]
##' output
##' }
##'
##' ## observed EPC
##' extra(fit.scalar)
##'
##' ## permutation results, including extra output
##' set.seed(12345) # same permutations
##' out.scalar <- permuteMeasEq(nPermute = 20, uncon = fit.metric, con = fit.scalar,
##' param = "intercepts", AFIs = myAFIs,
##' moreAFIs = moreAFIs, null = fit.null, extra = extra)
##' ## summarize extra output
##' summary(out.scalar, extra = TRUE)
##'
##'
##' ###########
##' ## MIMIC ##
##' ###########
##'
##' ## Specify Restricted Factor Analysis (RFA) model, equivalent to MIMIC, but
##' ## the factor covaries with the covariate instead of being regressed on it.
##' ## The covariate defines a single-indicator construct, and the
##' ## double-mean-centered products of the indicators define a latent
##' ## interaction between the factor and the covariate.
##' mod.mimic <- '
##' visual =~ x1 + x2 + x3
##' age =~ ageyr
##' age.by.vis =~ x1.ageyr + x2.ageyr + x3.ageyr
##'
##' x1 ~~ x1.ageyr
##' x2 ~~ x2.ageyr
##' x3 ~~ x3.ageyr
##' '
##'
##' HS.orth <- indProd(var1 = paste0("x", 1:3), var2 = "ageyr", match = FALSE,
##' data = HS[ , c("ageyr", paste0("x", 1:3))] )
##' fit.mimic <- cfa(mod.mimic, data = HS.orth, meanstructure = TRUE)
##' summary(fit.mimic, stand = TRUE)
##'
##' ## Whereas MIMIC models specify direct effects of the covariate on an indicator,
##' ## DIF can be tested in RFA models by specifying free loadings of an indicator
##' ## on the covariate's construct (uniform DIF, scalar invariance) and the
##' ## interaction construct (nonuniform DIF, metric invariance).
##' param <- as.list(paste0("age + age.by.vis =~ x", 1:3))
##' names(param) <- paste0("x", 1:3)
##' # param <- as.list(paste0("x", 1:3, " ~ age + age.by.vis")) # equivalent
##'
##' ## test both parameters simultaneously for each indicator
##' do.call(rbind, lapply(param, function(x) lavTestScore(fit.mimic, add = x)$test))
##' ## or test each parameter individually
##' lavTestScore(fit.mimic, add = as.character(param))
##'
##'
##' ####################### Permutation Method
##'
##' ## function to recalculate interaction terms after permuting the covariate
##' datafun <- function(data) {
##' d <- data[, c(paste0("x", 1:3), "ageyr")]
##' indProd(var1 = paste0("x", 1:3), var2 = "ageyr", match = FALSE, data = d)
##' }
##'
##' set.seed(12345)
##' perm.mimic <- permuteMeasEq(nPermute = 20, modelType = "mimic",
##' con = fit.mimic, param = param,
##' covariates = "ageyr", datafun = datafun)
##' summary(perm.mimic)
##'
##' }
##'
##' @export
permuteMeasEq <- function(nPermute, modelType = c("mgcfa","mimic"),
con, uncon = NULL, null = NULL,
param = NULL, freeParam = NULL, covariates = NULL,
AFIs = NULL, moreAFIs = NULL,
maxSparse = 10, maxNonconv = 10, showProgress = TRUE,
warn = -1, datafun, extra,
parallelType = c("none","multicore","snow"),
ncpus = NULL, cl = NULL, iseed = 12345) {
## save arguments from call
availableArgs <- as.list(formals(permuteMeasEq))
argNames <- names(availableArgs)
if (missing(datafun)) argNames <- setdiff(argNames, "datafun")
if (missing(extra)) argNames <- setdiff(argNames, "extra")
for (aa in argNames) {
if (!is.null(eval(as.name(aa))))
suppressWarnings(availableArgs[[aa]] <- eval(as.name(aa)))
}
## check and return them
fullCall <- do.call(checkPermArgs, availableArgs)
## assign them to workspace (also adds old_RNG & oldSeed to workspace)
for (aa in names(fullCall)) assign(aa, fullCall[[aa]])
###################### SAVE OBSERVED RESULTS ##########################
AFI.obs <- do.call(getAFIs, fullCall)
## save modification indices if !is.null(param)
if (is.null(param)) {
MI.obs <- data.frame(NULL)
} else MI.obs <- do.call(getMIs, fullCall)
## anything extra?
if (!missing(extra)) {
extraArgs <- formals(extra)
neededArgs <- intersect(names(extraArgs), names(fullCall))
extraArgs <- do.call(c, lapply(neededArgs, function(nn) fullCall[nn]))
extraOut <- do.call(extra, extraArgs)
## check that extra() returns a named list of scalars
if (!is.list(extraOut)) extraOut <- as.list(extraOut)
wrongFormat <- paste('Function "extra" must return a numeric vector or a',
'list of scalars, with each element named.')
if (!all(sapply(extraOut, is.numeric))) stop(wrongFormat)
if (!all(sapply(extraOut, length) == 1L)) stop(wrongFormat)
if (is.null(names(extraOut)) | any(names(extraOut) == "")) stop(wrongFormat)
extra.obs <- do.call(c, extraOut)
} else extra.obs <- numeric(length = 0L)
######################### PREP DATA ##############################
argList <- fullCall[c("con","uncon","null","param","freeParam","covariates",
"AFIs","moreAFIs","maxSparse","maxNonconv","warn","iseed")]
argList$G <- lavInspect(con, "group")
## check for categorical variables
# catVars <- lavaan::lavNames(con, type = "ov.ord")
# numVars <- lavaan::lavNames(con, type = "ov.num")
# latentVars <- lavaan::lavNames(con, type = "lv.regular")
## assemble data to which the models were fit
if (length(argList$G)) {
dataList <- mapply(FUN = function(x, g, n) {
y <- data.frame(as.data.frame(x), g, stringsAsFactors = FALSE)
names(y) <- c(n, argList$G)
y
}, SIMPLIFY = FALSE,
x = lavInspect(con, "data"), g = lavInspect(con, "group.label"),
n = lavaan::lavNames(con, type = "ov",
group = seq_along(lavInspect(con, "group.label"))))
argList$d <- do.call(rbind, dataList)
} else {
argList$d <- as.data.frame(lavInspect(con, "data"))
names(argList$d) <- lavaan::lavNames(con, type = "ov")
}
## check that covariates are actual variables
if (modelType == "mimic") {
if (length(covariates) && !all(covariates %in% names(argList$d)))
stop('These specified covariates are not columns in the data.frame:\n',
paste(setdiff(covariates, names(argList$d)), collapse = ", "))
}
## anything extra?
if (!missing(extra)) argList$extra <- extra
if (!missing(datafun)) argList$datafun <- datafun
###################### PERMUTED RESULTS ###########################
## permute and return distributions of (delta)AFIs, largest MI, and extras
if (showProgress) {
mypb <- utils::txtProgressBar(min = 1, max = nPermute, initial = 1,
char = "=", width = 50, style = 3, file = "")
permuDist <- list()
for (j in 1:nPermute) {
permuDist[[j]] <- do.call(paste("permuteOnce", modelType, sep = "."),
args = c(argList, i = j))
utils::setTxtProgressBar(mypb, j)
}
close(mypb)
} else if (parallelType == "multicore") {
if (length(iseed)) set.seed(iseed)
argList$FUN <- as.name(paste("permuteOnce", modelType, sep = "."))
argList$X <- 1:nPermute
argList$mc.cores <- ncpus
argList$mc.set.seed <- TRUE
pmcl <- function(...) { parallel::mclapply(...) }
permuDist <- do.call(pmcl, args = argList)
## restore old RNG type
if (fullCall$old_RNG[1] != "L'Ecuyer-CMRG") RNGkind(fullCall$old_RNG[1])
} else if (parallelType == "snow") {
stopTheCluster <- FALSE
if (is.null(cl)) {
stopTheCluster <- TRUE
cl <- parallel::makePSOCKcluster(rep("localhost", ncpus))
}
parallel::clusterSetRNGStream(cl, iseed = iseed)
argList$cl <- cl
argList$X <- 1:nPermute
argList$fun <- paste("permuteOnce", modelType, sep = ".")
parallel::clusterExport(cl, varlist = c(argList$fun, "getAFIs","getMIs")) #FIXME: need update?
tempppl <- function(...) { parallel::parLapply(...) }
permuDist <- do.call(tempppl, args = argList)
if (stopTheCluster) parallel::stopCluster(cl)
## restore old RNG type
if (fullCall$old_RNG[1] != "L'Ecuyer-CMRG") RNGkind(fullCall$old_RNG[1])
} else {
argList$X <- 1:nPermute
argList$FUN <- paste("permuteOnce", modelType, sep = ".")
permuDist <- do.call(lapply, args = argList)
}
## extract AFI distribution
if (length(AFI.obs) > 1) {
AFI.dist <- as.data.frame(t(sapply(permuDist, function(x) x$AFI)))
}
if (length(AFI.obs) == 1L) {
AFI.dist <- data.frame(sapply(permuDist, function(x) x$AFI))
colnames(AFI.dist) <- names(AFI.obs)
}
## identify badness-of-fit measures
badness <- grepl(pattern = "fmin|chi|aic|bic|rmr|rmsea|cn|sic|hqc",
x = names(AFI.obs), ignore.case = TRUE)
## calculate all one-directional p-values
AFI.pval <- mapply(FUN = function(x, y, b) {
if (b) return(mean(x >= y, na.rm = TRUE))
mean(x <= y, na.rm = TRUE)
}, x = unclass(AFI.dist), y = AFI.obs, b = badness)
## extract distribution of maximum modification indices
MI.dist <- as.numeric(unlist(lapply(permuDist, function(x) x$MI)))
## calculate Tukey-adjusted p values for modification indices
if (!is.null(param)) {
MI.obs$tukey.p.value <- sapply(MI.obs$X2,
function(i) mean(i <= MI.dist, na.rm = TRUE))
MI.obs <- as.data.frame(unclass(MI.obs))
rownames(MI.obs) <- names(param)
}
## anything extra?
if (!missing(extra)) {
extra.dist <- do.call(rbind, lapply(permuDist, function(x) x$extra))
} else extra.dist <- data.frame(NULL)
## save parameter table for show/summary methods
PT <- as.data.frame(parTable(con))
PT$par <- paste0(PT$lhs, PT$op, PT$rhs)
if (length(lavInspect(con, "group")))
PT$group.label[PT$group > 0] <- lavInspect(con, "group.label")[PT$group[PT$group > 0] ]
## return observed results, permutation p values, and ANOVA results
if (is.null(uncon)) {
delta <- lavaan::anova(con)
} else {
delta <- lavaan::anova(uncon, con)
}
ANOVA <- sapply(delta[,c("Chisq diff","Df diff","Pr(>Chisq)")], function(x) x[2])
out <- new("permuteMeasEq", PT = PT, modelType = modelType, ANOVA = ANOVA,
AFI.obs = AFI.obs, AFI.dist = AFI.dist, AFI.pval = AFI.pval,
MI.obs = MI.obs, MI.dist = MI.dist,
extra.obs = extra.obs, extra.dist = extra.dist,
n.Permutations = nPermute, n.Converged = sum(!is.na(AFI.dist[,1])),
n.nonConverged = sapply(permuDist, function(x) x$n.nonConverged),
n.Sparse = sapply(permuDist, function(x) x$n.Sparse),
oldSeed = fullCall$oldSeed)
out
}
## ----------------
## Hidden Functions
## ----------------
## function to check validity of arguments to permuteMeasEq()
#' @importFrom lavaan lavInspect parTable
checkPermArgs <- function(nPermute, modelType, con, uncon, null,
param, freeParam, covariates, AFIs, moreAFIs,
maxSparse, maxNonconv, showProgress, warn,
datafun, extra, parallelType, ncpus, cl, iseed) {
fixedCall <- as.list(match.call())[-1]
fixedCall$nPermute <- as.integer(nPermute[1])
fixedCall$modelType <- modelType[1]
if (!fixedCall$modelType %in% c("mgcfa","mimic","long"))
stop('modelType must be one of c("mgcfa","mimic","long")')
if (fixedCall$modelType == "long") stop('modelType "long" is not yet available.')
if (fixedCall$modelType == "mgcfa" && lavInspect(con, "ngroups") == 1L)
stop('modelType = "mgcfa" applies only to multigroup models.')
if (fixedCall$modelType == "mimic") {
uncon <- NULL
fixedCall$uncon <- NULL
fixedCall <- c(fixedCall, list(uncon = NULL))
}
## strip white space
if (is.list(param)) {
fixedCall$param <- lapply(param, function(cc) gsub("[[:space:]]+", "", cc))
} else if (!is.null(param)) fixedCall$param <- gsub("[[:space:]]+", "", param)
if (!is.null(freeParam)) fixedCall$freeParam <- gsub("[[:space:]]+", "", freeParam)
if (fixedCall$modelType == "mimic") {
# PT <- lavaan::lavaanify(fixedCall$param)
# checkCovs <- unique(PT$rhs[PT$op == "~"])
# if (is.null(covariates)) covariates <- checkCovs
# if (length(setdiff(covariates, checkCovs)))
# warning('Argument "covariates" includes predictors not in argument "param"')
##### ordVars <- lavaan::lavNames(con, type = "ov.ord")
fixedCall$covariates <- as.character(covariates)
}
fixedCall$maxSparse <- as.integer(maxSparse[1])
fixedCall$maxNonconv <- as.integer(maxNonconv[1])
fixedCall$showProgress <- as.logical(showProgress[1])
fixedCall$warn <- as.integer(warn[1])
fixedCall$oldSeed <- as.integer(NULL)
parallelType <- as.character(parallelType[1])
if (!parallelType %in% c("none","multicore","snow")) parallelType <- "none"
if (!is.null(cl)) {
if (!is(cl, "cluster")) stop("Invalid cluster object. Check class(cl)")
parallelType <- "snow"
ncpus <- length(cl)
}
if (parallelType == "multicore" && .Platform$OS.type == "windows") {
parallelType <- "snow"
message("'multicore' option unavailable on Windows. Using 'snow' instead.")
}
## parallel settings, adapted from boot::boot()
if (parallelType != "none") {
if (is.null(ncpus) || ncpus > parallel::detectCores()) {
ncpus <- parallel::detectCores() - 1
}
if (ncpus <= 1L) {
parallelType <- "none"
} else {
fixedCall$showProgress <- FALSE
fixedCall$old_RNG <- RNGkind()
fixedCall$oldSeed <- .Random.seed
if (fixedCall$old_RNG[1] != "L'Ecuyer-CMRG") {
RNGkind("L'Ecuyer-CMRG")
message("Your RNGkind() was changed from ", fixedCall$old_RNG[1],
" to L'Ecuyer-CMRG, which is required for reproducibility ",
" in parallel jobs. Your RNGkind() has been returned to ",
fixedCall$old_RNG[1], " but the seed has not been set. ",
" The state of your previous RNG is saved in the slot ",
" named 'oldSeed', if you want to restore it using ",
" the syntax:\n",
".Random.seed[-1] <- permuteMeasEqObject@oldSeed[-1]")
}
fixedCall$iseed <- as.integer(iseed[1])
if (is.na(fixedCall$iseed)) fixedCall$iseed <- 12345
}
}
fixedCall$parallelType <- parallelType
if (is.null(ncpus)) {
fixedCall$ncpus <- NULL
fixedCall <- c(fixedCall, list(ncpus = NULL))
} else fixedCall$ncpus <- ncpus
## Check that "param" is NULL if uncon is NULL, and check for lavaan class.
## Also check that models are fitted to raw data, not summary stats.
notLavaan <- "Non-NULL 'con', 'uncon', or 'null' must be fitted lavaan object."
notRawData <- "lavaan models ('con', 'uncon', or 'null') must be fitted to raw data=, not summary statistics (e.g., sample.cov=)"
if (is.null(uncon)) {
if (!is.null(fixedCall$param) && fixedCall$modelType == "mgcfa") {
message(c(" When 'uncon = NULL', only configural invariance is tested.",
"\n So the 'param' argument was changed to NULL."))
fixedCall$param <- NULL
fixedCall <- c(fixedCall, list(param = NULL))
}
if (!inherits(con, "lavaan")) stop(notLavaan)
stopifnot(con@Data@data.type == "full")
} else {
if (!inherits(con, "lavaan")) stop(notLavaan)
if (!inherits(uncon, "lavaan")) stop(notLavaan)
stopifnot( con@Data@data.type == "full")
stopifnot(uncon@Data@data.type == "full")
}
if (!is.null(null)) {
if (!inherits(null, "lavaan")) stop(notLavaan)
stopifnot(null@Data@data.type == "full")
}
############ FIXME: check that lavInspect(con, "options")$conditional.x = FALSE (find defaults for continuous/ordered indicators)
if (!is.null(fixedCall$param)) {
## Temporarily warn about testing thresholds without necessary constraints. FIXME: check for binary indicators
if ("thresholds" %in% fixedCall$param | any(grepl("\\|", fixedCall$param))) {
warning(c("This function is not yet optimized for testing thresholds.\n",
"Necessary identification contraints might not be specified."))
}
## collect parameter types for "mgcfa"
if (fixedCall$modelType != "mimic") {
## save all estimates from constrained model
PT <- parTable(con)[ , c("lhs","op","rhs","group","plabel")]
## extract parameters of interest
paramTypes <- c("loadings","intercepts","thresholds","residuals","means",
"residual.covariances","lv.variances","lv.covariances")
params <- PT[paste0(PT$lhs, PT$op, PT$rhs) %in% setdiff(fixedCall$param,
paramTypes), ]
## add parameters by type, if any are specified
types <- intersect(fixedCall$param, paramTypes)
ov.names <- lavaan::lavNames(con, "ov")
isOV <- PT$lhs %in% ov.names
lv.names <- con@pta$vnames$lv[[1]]
isLV <- PT$lhs %in% lv.names & PT$rhs %in% lv.names
if ("loadings" %in% types) params <- rbind(params, PT[PT$op == "=~", ])
if ("intercepts" %in% types) {
params <- rbind(params, PT[isOV & PT$op == "~1", ])
}
if ("thresholds" %in% types) params <- rbind(params, PT[PT$op == "|", ])
if ("residuals" %in% types) {
params <- rbind(params, PT[isOV & PT$lhs == PT$rhs & PT$op == "~~", ])
}
if ("residual.covariances" %in% types) {
params <- rbind(params, PT[isOV & PT$lhs != PT$rhs & PT$op == "~~", ])
}
if ("means" %in% types) {
params <- rbind(params, PT[PT$lhs %in% lv.names & PT$op == "~1", ])
}
if ("lv.variances" %in% types) {
params <- rbind(params, PT[isLV & PT$lhs == PT$rhs & PT$op == "~~", ])
}
if ("lv.covariances" %in% types) {
params <- rbind(params, PT[isLV & PT$lhs != PT$rhs & PT$op == "~~", ])
}
## remove parameters specified by "freeParam" argument
params <- params[!paste0(params$lhs, params$op, params$rhs) %in% fixedCall$freeParam, ]
fixedCall$param <- paste0(params$lhs, params$op, params$rhs)
}
}
if (is.null(AFIs) & is.null(moreAFIs)) {
message("No AFIs were selected, so only chi-squared will be permuted.\n")
fixedCall$AFIs <- "chisq"
AFIs <- "chisq"
}
if ("ecvi" %in% AFIs & lavInspect(con, "ngroups") > 1L)
stop("ECVI is not available for multigroup models.")
## check estimators
leastSq <- grepl("LS", lavInspect(con, "options")$estimator)
if (!is.null(uncon)) {
if (uncon@Options$estimator != lavInspect(con, "options")$estimator)
stop("Models must be fit using same estimator.")
}
if (!is.null(null)) {
if (lavInspect(null, "options")$estimator != lavInspect(con, "options")$estimator)
stop("Models must be fit using same estimator.")
}
## check extra functions, if any
restrictedArgs <- c("con","uncon","null","param","freeParam","covariates",
"AFIs","moreAFIs","maxSparse","maxNonconv","iseed")
if (!missing(datafun)) {
if (!is.function(datafun)) stop('Argument "datafun" must be a function.')
extraArgs <- formals(datafun)
if (!all(names(extraArgs) %in% c(restrictedArgs, "data")))
stop('The user-supplied function "datafun" can only have any among the ',
'following arguments:\n', paste(restrictedArgs, collapse = ", "))
}
if (!missing(extra)) {
if (!is.function(extra)) stop('Argument "extra" must be a function.')
extraArgs <- formals(extra)
if (!all(names(extraArgs) %in% restrictedArgs))
stop('The user-supplied function "extra" can only have any among the ',
'following arguments:\n', paste(restrictedArgs, collapse = ", "))
}
## return evaluated list of other arguments
lapply(fixedCall, eval)
}
## function to extract fit measures
#' @importFrom lavaan lavInspect
getAFIs <- function(...) {
dots <- list(...)
AFI1 <- list()
AFI0 <- list()
leastSq <- grepl("LS", lavInspect(dots$con, "options")$estimator)
## check validity of user-specified AFIs, save output
if (!is.null(dots$AFIs)) {
IC <- grep("ic|logl", dots$AFIs, value = TRUE)
if (leastSq & length(IC)) {
stop(paste("Argument 'AFIs' includes invalid options:",
paste(IC, collapse = ", "),
"Information criteria unavailable for least-squares estimators.",
sep = "\n"))
}
if (!is.null(dots$uncon))
AFI1[[1]] <- lavaan::fitMeasures(dots$uncon, fit.measures = dots$AFIs,
baseline.model = dots$null)
AFI0[[1]] <- lavaan::fitMeasures(dots$con, fit.measures = dots$AFIs,
baseline.model = dots$null)
}
## check validity of user-specified moreAFIs
if (!is.null(dots$moreAFIs)) {
IC <- grep("ic|hqc", dots$moreAFIs, value = TRUE)
if (leastSq & length(IC)) {
stop(paste("Argument 'moreAFIs' includes invalid options:",
paste(IC, collapse = ", "),
"Information criteria unavailable for least-squares estimators.",
sep = "\n"))
}
if (!is.null(dots$uncon))
AFI1[[2]] <- moreFitIndices(dots$uncon, fit.measures = dots$moreAFIs)
AFI0[[2]] <- moreFitIndices(dots$con, fit.measures = dots$moreAFIs)
}
## save observed AFIs or delta-AFIs
if (is.null(dots$uncon)) {
AFI.obs <- unlist(AFI0)
} else {
AFI.obs <- unlist(AFI0) - unlist(AFI1)
}
AFI.obs
}
## Function to extract modification indices for equality constraints
#' @importFrom lavaan parTable
getMIs <- function(...) {
dots <- list(...)
if (dots$modelType == "mgcfa") {
## save all estimates from constrained model
PT <- parTable(dots$con)[ , c("lhs","op","rhs","group","plabel")]
## extract parameters of interest
params <- PT[paste0(PT$lhs, PT$op, PT$rhs) %in% dots$param, ]
## return modification indices for specified constraints (param)
MIs <- lavaan::lavTestScore(dots$con)$uni
MI.obs <- MIs[MIs$lhs %in% params$plabel, ]
} else if (dots$modelType == "mimic") {
if (is.list(dots$param)) {
MI <- lapply(dots$param, function(x) lavaan::lavTestScore(dots$con, add = x)$test)
MI.obs <- do.call(rbind, MI)
} else MI.obs <- lavaan::lavTestScore(dots$con, add = dots$param)$uni
} else if (dots$modelType == "long") {
## coming soon
}
MI.obs
}
## Functions to find delta-AFIs & maximum modification index in one permutation
permuteOnce.mgcfa <- function(i, d, G, con, uncon, null, param, freeParam,
covariates, AFIs, moreAFIs, maxSparse, maxNonconv,
iseed, warn, extra = NULL, datafun = NULL) {
old_warn <- options()$warn
options(warn = warn)
## save arguments from call
argNames <- names(formals(permuteOnce.mgcfa))
availableArgs <- lapply(argNames, function(x) eval(as.name(x)))
names(availableArgs) <- argNames
group.label <- lavaan::lavInspect(con, "group.label")
nSparse <- 0L
nTries <- 1L
while ( (nSparse <= maxSparse) & (nTries <= maxNonconv) ) {
## permute grouping variable
d[ , G] <- sample(d[ , G])
## transform data?
if (!is.null(datafun)) {
extraArgs <- formals(datafun)
neededArgs <- intersect(names(extraArgs), names(availableArgs))
extraArgs <- do.call(c, lapply(neededArgs, function(nn) availableArgs[nn]))
extraArgs$data <- d
originalNames <- colnames(d)
d <- do.call(datafun, extraArgs)
## coerce extraOut to data.frame
if (!is.data.frame(d)) stop('Argument "datafun" did not return a data.frame')
if (!all(originalNames %in% colnames(d)))
stop('The data.frame returned by argument "datafun" did not contain ',
'column names required by the model:\n',
paste(setdiff(originalNames, colnames(d)), collapse = ", "))
}
## for ordered indicators, check that groups have same observed categories
ordVars <- lavaan::lavNames(con, type = "ov.ord")
if (length(ordVars) > 0) {
try(onewayTables <- lavaan::lavTables(d, dimension = 1L,
categorical = ordVars, group = G),
silent = TRUE)
if (exists("onewayTables")) {
if (any(onewayTables$obs.prop == 1)) {
nSparse <- nSparse + 1L
next
}
} else {
## no "onewayTables" probably indicates empty categories in 1+ groups
nSparse <- nSparse + 1L
next
}
}
## fit null model, if it exists
if (!is.null(null)) {
out.null <- lavaan::update(null, data = d, group.label = group.label)
}
## fit constrained model, check for convergence
try(out0 <- lavaan::update(con, data = d, group.label = group.label))
if (!exists("out0")) {
nTries <- nTries + 1L
next
}
if (!lavaan::lavInspect(out0, "converged")) {
nTries <- nTries + 1L
next
}
## fit unconstrained model (unless NULL), check for convergence
if (!is.null(uncon)) {
try(out1 <- lavaan::update(uncon, data = d, group.label = group.label))
if (!exists("out1")) {
nTries <- nTries + 1L
next
}
if (!lavaan::lavInspect(out1, "converged")) {
nTries <- nTries + 1L
next
}
}
## If you get this far, everything converged, so break WHILE loop
break
}
## if WHILE loop ended before getting results, return NA
if ( (nSparse == maxSparse) | (nTries == maxNonconv) ) {
allAFIs <- c(AFIs, moreAFIs)
AFI <- rep(NA, sum(!is.na(allAFIs)))
names(AFI) <- allAFIs[!is.na(allAFIs)]
MI <- if (is.null(param)) NULL else NA
extra.obs <- NA
nTries <- nTries + 1L
} else {
availableArgs$con <- out0
if (exists("out1")) availableArgs$uncon <- out1
if (exists("out.null")) availableArgs$null <- out.null
AFI <- do.call(getAFIs, availableArgs)
## save max(MI) if !is.null(param)
if (is.null(param)) {
MI <- NULL
} else {
MI <- max(do.call(getMIs, c(availableArgs, modelType = "mgcfa"))$X2)
}
## anything extra?
if (!is.null(extra)) {
extraArgs <- formals(extra)
neededArgs <- intersect(names(extraArgs), names(availableArgs))
extraArgs <- do.call(c, lapply(neededArgs, function(nn) availableArgs[nn]))
extraOut <- do.call(extra, extraArgs)
## coerce extraOut to data.frame
if (!is.list(extraOut)) extraOut <- as.list(extraOut)
extra.obs <- data.frame(extraOut)
} else extra.obs <- data.frame(NULL)
}
options(warn = old_warn)
list(AFI = AFI, MI = MI, extra = extra.obs,
n.nonConverged = nTries - 1L, n.Sparse = nSparse)
}
permuteOnce.mimic <- function(i, d, G, con, uncon, null, param, freeParam,
covariates, AFIs, moreAFIs, maxSparse, maxNonconv,
iseed, warn, extra = NULL, datafun = NULL) {
old_warn <- options()$warn
options(warn = warn)
## save arguments from call
argNames <- names(formals(permuteOnce.mimic))
availableArgs <- lapply(argNames, function(x) eval(as.name(x)))
names(availableArgs) <- argNames
group.label <- lavaan::lavInspect(con, "group.label")
nTries <- 1L
while (nTries <= maxNonconv) {
## permute covariate(s) within each group
if (length(G)) {
for (gg in group.label) {
dG <- d[ d[[G]] == gg, ]
N <- nrow(dG)
newd <- dG[sample(1:N, N), covariates, drop = FALSE]
for (COV in covariates) d[d[[G]] == gg, COV] <- newd[ , COV]
}
} else {
N <- nrow(d)
newd <- d[sample(1:N, N), covariates, drop = FALSE]
for (COV in covariates) d[ , COV] <- newd[ , COV]
}
## transform data?
if (!is.null(datafun)) {
extraArgs <- formals(datafun)
neededArgs <- intersect(names(extraArgs), names(availableArgs))
extraArgs <- do.call(c, lapply(neededArgs, function(nn) availableArgs[nn]))
extraArgs$data <- d
originalNames <- colnames(d)
d <- do.call(datafun, extraArgs)
## coerce extraOut to data.frame
if (!is.data.frame(d)) stop('Argument "datafun" did not return a data.frame')
if (!all(originalNames %in% colnames(d)))
stop('The data.frame returned by argument "datafun" did not contain ',
'column names required by the model:\n',
paste(setdiff(originalNames, colnames(d)), collapse = ", "))
}
## fit null model, if it exists
if (!is.null(null)) {
out.null <- lavaan::update(null, data = d, group.label = group.label)
}
## fit constrained model
try(out0 <- lavaan::update(con, data = d, group.label = group.label))
## check for convergence
if (!exists("out0")) {
nTries <- nTries + 1L
next
}
if (!lavaan::lavInspect(out0, "converged")) {
nTries <- nTries + 1L
next
}
## If you get this far, everything converged, so break WHILE loop
break
}
## if WHILE loop ended before getting results, return NA
if (nTries == maxNonconv) {
allAFIs <- c(AFIs, moreAFIs)
AFI <- rep(NA, sum(!is.na(allAFIs)))
names(AFI) <- allAFIs[!is.na(allAFIs)]
MI <- if (is.null(param)) NULL else NA
extra.obs <- NA
nTries <- nTries + 1L
} else {
availableArgs$con <- out0
if (exists("out.null")) availableArgs$null <- out.null
AFI <- do.call(getAFIs, availableArgs)
if (is.null(param)) {
MI <- NULL
} else {
MI <- max(do.call(getMIs, c(availableArgs, modelType = "mimic"))$X2)
}
## anything extra?
if (!is.null(extra)) {
extraArgs <- formals(extra)
neededArgs <- intersect(names(extraArgs), names(availableArgs))
extraArgs <- do.call(c, lapply(neededArgs, function(nn) availableArgs[nn]))
extraOut <- do.call(extra, extraArgs)
## coerce extraOut to data.frame
if (!is.list(extraOut)) extraOut <- as.list(extraOut)
extra.obs <- data.frame(extraOut)
} else extra.obs <- data.frame(NULL)
}
options(warn = old_warn)
list(AFI = AFI, MI = MI, extra = extra.obs,
n.nonConverged = nTries - 1L, n.Sparse = integer(length = 0))
}
Try the semTools package in your browser
Any scripts or data that you put into this service are public.
semTools documentation built on April 3, 2025, 9:23 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions summ.mimic summ.mgcfa
### Terrence D. Jorgensen
### Last updated: 12 March 2025
### permutation randomization test for measurement equivalence and DIF
## -----------------
## Class and Methods
## -----------------
##' Class for the Results of Permutation Randomization Tests of Measurement
##' Equivalence and DIF
##'
##' This class contains the results of tests of Measurement Equivalence and
##' Differential Item Functioning (DIF).
##'
##'
##' @name permuteMeasEq-class
##' @aliases permuteMeasEq-class show,permuteMeasEq-method
##' summary,permuteMeasEq-method hist,permuteMeasEq-method
##' @docType class
##'
##' @slot PT A `data.frame` returned by a call to
##' [lavaan::parTable()] on the constrained model
##' @slot modelType A character indicating the specified `modelType` in the
##' call to `permuteMeasEq`
##' @slot ANOVA A `numeric` vector indicating the results of the observed
##' (\eqn{\Delta})\eqn{\chi^2} test, based on the central \eqn{\chi^2}
##' distribution
##' @slot AFI.obs A vector of observed (changes in) user-selected fit measures
##' @slot AFI.dist The permutation distribution(s) of user-selected fit measures.
##' A `data.frame` with `n.Permutations` rows and one column for each
##' `AFI.obs`.
##' @slot AFI.pval A vector of *p* values (one for each element in slot
##' `AFI.obs`) calculated using slot `AFI.dist`, indicating the
##' probability of observing a change at least as extreme as `AFI.obs`
##' if the null hypothesis were true
##' @slot MI.obs A `data.frame` of observed Lagrange Multipliers
##' (modification indices) associated with the equality constraints or fixed
##' parameters specified in the `param` argument. This is a subset of the
##' output returned by a call to [lavaan::lavTestScore()] on the
##' constrained model.
##' @slot MI.dist The permutation distribution of the maximum modification index
##' (among those seen in slot `MI.obs$X2`) at each permutation of group
##' assignment or of `covariates`
##' @slot extra.obs If `permuteMeasEq` was called with an `extra`
##' function, the output when applied to the original data is concatenated
##' into this vector
##' @slot extra.dist A `data.frame`, each column of which contains the
##' permutation distribution of the corresponding statistic in slot
##' `extra.obs`
##' @slot n.Permutations An `integer` indicating the number of permutations
##' requested by the user
##' @slot n.Converged An `integer` indicating the number of permuation
##' iterations which yielded a converged solution
##' @slot n.nonConverged An `integer` vector of length
##' `n.Permutations` indicating how many times group assignment was
##' randomly permuted (at each iteration) before converging on a solution
##' @slot n.Sparse Only relevant with `ordered` indicators when
##' `modelType == "mgcfa"`. An `integer` vector of length
##' `n.Permutations` indicating how many times group assignment was
##' randomly permuted (at each iteration) before obtaining a sample with all
##' categories observed in all groups.
##' @slot oldSeed An `integer` vector storing the value of
##' `.Random.seed` before running `permuteMeasEq`. Only relevant
##' when using a parallel/multicore option and the original
##' `RNGkind() != "L'Ecuyer-CMRG"`. This enables users to restore their
##' previous `.Random.seed` state, if desired, by running:
##' `.Random.seed[-1] <- permutedResults@oldSeed[-1]`
##' @section Objects from the Class: Objects can be created via the
##' [semTools::permuteMeasEq()] function.
##'
##' @return
##' \itemize{
##' \item The `show` method prints a summary of the multiparameter
##' omnibus test results, using the user-specified AFIs. The parametric
##' (\eqn{\Delta})\eqn{\chi^2} test is also displayed.
##' \item The `summary` method prints the same information from the
##' `show` method, but when `extra = FALSE` (the default) it also
##' provides a table summarizing any requested follow-up tests of DIF using
##' modification indices in slot `MI.obs`. The user can also specify an
##' `alpha` level for flagging modification indices as significant, as
##' well as `nd` (the number of digits displayed). For each modification
##' index, the *p* value is displayed using a central \eqn{\chi^2}
##' distribution with the *df* shown in that column. Additionally, a
##' *p* value is displayed using the permutation distribution of the
##' maximum index, which controls the familywise Type I error rate in a manner
##' similar to Tukey's studentized range test. If any indices are flagged as
##' significant using the `tukey.p.value`, then a message is displayed for
##' each flagged index. The invisibly returned `data.frame` is the
##' displayed table of modification indices, unless
##' [semTools::permuteMeasEq()] was called with `param = NULL`,
##' in which case the invisibly returned object is `object`. If
##' `extra = TRUE`, the permutation-based *p* values for each
##' statistic returned by the `extra` function are displayed and returned
##' in a `data.frame` instead of the modification indices requested in the
##' `param` argument.
##' \item The `hist` method returns a list of `length == 2`,
##' containing the arguments for the call to `hist` and the arguments
##' to the call for `legend`, respectively. This list may facilitate
##' creating a customized histogram of `AFI.dist`, `MI.dist`, or
##' `extra.dist`
##' }
##'
##' @author Terrence D. Jorgensen (University of Amsterdam;
##' \email{TJorgensen314@@gmail.com})
##'
##' @seealso [semTools::permuteMeasEq()]
##'
##' @examples
##'
##' # See the example from the permuteMeasEq function
##'
setClass("permuteMeasEq", slots = c(PT = "data.frame",
modelType = "character",
ANOVA = "vector",
AFI.obs = "vector",
AFI.dist = "data.frame",
AFI.pval = "vector",
MI.obs = "data.frame",
MI.dist = "vector",
extra.obs = "vector",
extra.dist = "data.frame",
n.Permutations = "integer",
n.Converged = "integer",
n.nonConverged = "vector",
n.Sparse = "vector",
oldSeed = "integer"))
##' @rdname permuteMeasEq-class
##' @aliases show,permuteMeasEq-method
##' @export
setMethod("show", "permuteMeasEq", function(object) {
## print warning if there are nonConverged permutations
if (object@n.Permutations != object@n.Converged) {
warning(paste("Only", object@n.Converged, "out of",
object@n.Permutations, "models converged within",
max(object@n.nonConverged), "attempts per permutation.\n\n"))
}
## print ANOVA
cat("Omnibus p value based on parametric chi-squared difference test:\n\n")
print(round(object@ANOVA, digits = 3))
## print permutation results
cat("\n\nOmnibus p values based on nonparametric permutation method: \n\n")
AFI <- data.frame(AFI.Difference = object@AFI.obs, p.value = object@AFI.pval)
class(AFI) <- c("lavaan.data.frame","data.frame")
print(AFI, nd = 3)
invisible(object)
})
##' @rdname permuteMeasEq-class
##' @aliases summary,permuteMeasEq-method
##' @export
setMethod("summary", "permuteMeasEq", function(object, alpha = .05, nd = 3,
extra = FALSE) {
## print warning if there are nonConverged permutations
if (object@n.Permutations != object@n.Converged) {
warning(paste("Only", object@n.Converged, "out of",
object@n.Permutations, "models converged within",
max(object@n.nonConverged), "attempts per permutation.\n\n"))
}
## print ANOVA
cat("Omnibus p value based on parametric chi-squared difference test:\n\n")
print(round(object@ANOVA, digits = nd))
## print permutation results
cat("\n\nOmnibus p values based on nonparametric permutation method: \n\n")
AFI <- data.frame(AFI.Difference = object@AFI.obs, p.value = object@AFI.pval)
class(AFI) <- c("lavaan.data.frame","data.frame")
print(AFI, nd = nd)
## print extras or DIF test results, if any were requested
if (extra && length(object@extra.obs)) {
cat("\n\nUnadjusted p values of extra statistics,\n",
"based on permutation distribution of each statistic: \n\n")
MI <- data.frame(Statistic = object@extra.obs)
class(MI) <- c("lavaan.data.frame","data.frame")
MI$p.value <- sapply(names(object@extra.dist), function(nn) {
mean(abs(object@extra.dist[,nn]) >= abs(object@extra.obs[nn]), na.rm = TRUE)
})
MI$flag <- ifelse(MI$p.value < alpha, "* ", "")
print(MI, nd = nd)
} else if (length(object@MI.dist)) {
cat("\n\n Modification indices for equality constrained parameter estimates,\n",
"with unadjusted 'p.value' based on chi-squared distribution and\n",
"adjusted 'tukey.p.value' based on permutation distribution of the\n",
"maximum modification index per iteration: \n\n")
MI <- do.call(paste("summ", object@modelType, sep = "."),
args = list(object = object, alpha = alpha))
print(MI, nd = nd)
## print messages about potential DIF
if (all(MI$tukey.p.value > alpha)) {
cat("\n\n No equality constraints were flagged as significant.\n\n")
return(invisible(MI))
}
if (object@modelType == "mgcfa") {
cat("\n\nThe following equality constraints were flagged as significant:\n\n")
for (i in which(MI$tukey.p.value < alpha)) {
cat("Parameter '", MI$parameter[i], "' may differ between Groups '",
MI$group.lhs[i], "' and '", MI$group.rhs[i], "'.\n", sep = "")
}
cat("\nUse lavTestScore(..., epc = TRUE) on your constrained model to",
"display expected parameter changes for these equality constraints\n\n")
}
} else return(invisible(object))
invisible(MI)
})
summ.mgcfa <- function(object, alpha) {
MI <- object@MI.obs
class(MI) <- c("lavaan.data.frame","data.frame")
PT <- object@PT
eqPar <- rbind(PT[PT$plabel %in% MI$lhs, ], PT[PT$plabel %in% MI$rhs, ])
MI$flag <- ""
MI$parameter <- ""
MI$group.lhs <- ""
MI$group.rhs <- ""
for (i in 1:nrow(MI)) {
par1 <- eqPar$par[ eqPar$plabel == MI$lhs[i] ]
par2 <- eqPar$par[ eqPar$plabel == MI$rhs[i] ]
MI$parameter[i] <- par1
MI$group.lhs[i] <- eqPar$group.label[ eqPar$plabel == MI$lhs[i] ]
MI$group.rhs[i] <- eqPar$group.label[ eqPar$plabel == MI$rhs[i] ]
if (par1 != par2) {
myMessage <- paste0("Constraint '", MI$lhs[i], "==", MI$rhs[i],
"' refers to different parameters: \n'",
MI$lhs[i], "' is '", par1, "' in group '",
MI$group.lhs[i], "'\n'",
MI$rhs[i], "' is '", par2, "' in group '",
MI$group.rhs[i], "'\n")
warning(myMessage)
}
if (MI$tukey.p.value[i] < alpha) MI$flag[i] <- "* -->"
}
MI
}
summ.mimic <- function(object, alpha) {
MI <- object@MI.obs
class(MI) <- c("lavaan.data.frame","data.frame")
MI$flag <- ifelse(MI$tukey.p.value < alpha, "* ", "")
MI
}
##' @rdname permuteMeasEq-class
##' @aliases hist,permuteMeasEq-method
##' @importFrom stats qchisq dchisq quantile
##' @param object,x object of class `permuteMeasEq`
##' @param ... Additional arguments to pass to [graphics::hist()]
##' @param AFI `character` indicating the fit measure whose permutation
##' distribution should be plotted
##' @param alpha alpha level used to draw confidence limits in `hist` and
##' flag significant statistics in `summary` output
##' @param nd number of digits to display
##' @param extra `logical` indicating whether the `summary` output
##' should return permutation-based *p* values for each statistic returned
##' by the `extra` function. If `FALSE` (default), `summary`
##' will return permutation-based *p* values for each modification index.
##' @param printLegend `logical`. If `TRUE` (default), a legend will
##' be printed with the histogram
##' @param legendArgs `list` of arguments passed to the
##' [graphics::legend()] function. The default argument is a list
##' placing the legend at the top-left of the figure.
##' @export
setMethod("hist", "permuteMeasEq", function(x, ..., AFI, alpha = .05, nd = 3,
printLegend = TRUE,
legendArgs = list(x = "topleft")) {
histArgs <- list(...)
histArgs$x <- x@AFI.dist[[AFI]]
if (is.null(histArgs$col)) histArgs$col <- "grey69"
histArgs$freq <- !grepl("chi", AFI)
histArgs$ylab <- if (histArgs$freq) "Frequency" else "Probability Density"
if (printLegend) {
if (is.null(legendArgs$box.lty)) legendArgs$box.lty <- 0
if (nd < length(strsplit(as.character(1 / alpha), "")[[1]]) - 1) {
warning(paste0("The number of digits argument (nd = ", nd ,
") is too low to display your p value at the",
" same precision as your requested alpha level (alpha = ",
alpha, ")"))
}
if (x@AFI.pval[[AFI]] < (1 / 10^nd)) {
pVal <- paste(c("< .", rep(0, nd - 1),"1"), collapse = "")
} else {
pVal <- paste("=", round(x@AFI.pval[[AFI]], nd))
}
}
delta <- length(x@MI.dist) > 0L && x@modelType == "mgcfa"
if (grepl("chi", AFI)) { ####################################### Chi-squared
ChiSq <- x@AFI.obs[AFI]
DF <- x@ANOVA[2]
histArgs$xlim <- range(c(ChiSq, x@AFI.dist[[AFI]], qchisq(c(.01, .99), DF)))
xVals <- seq(histArgs$xlim[1], histArgs$xlim[2], by = .1)
theoDist <- dchisq(xVals, df = DF)
TheoCrit <- round(qchisq(p = alpha, df = DF, lower.tail = FALSE), 2)
Crit <- quantile(histArgs$x, probs = 1 - alpha)
if (ChiSq > histArgs$xlim[2]) histArgs$xlim[2] <- ChiSq
if (delta) {
histArgs$main <- expression(Permutation~Distribution~of~Delta*chi^2)
histArgs$xlab <- expression(Delta*chi^2)
if (printLegend) {
legendArgs$legend <- c(bquote(Theoretical~Delta*chi[Delta*.(paste("df =", DF))]^2 ~ Distribution),
bquote(Critical~chi[alpha~.(paste(" =", alpha))]^2 == .(round(TheoCrit, nd))),
bquote(.(paste("Permuted Critical Value =", round(Crit, nd)))),
bquote(Observed~Delta*chi^2 == .(round(ChiSq, nd))),
expression(paste("")),
bquote(Permuted~italic(p)~.(pVal)))
}
} else {
histArgs$main <- expression(Permutation~Distribution~of~chi^2)
histArgs$xlab <- expression(chi^2)
if (printLegend) {
legendArgs$legend <- c(bquote(Theoretical~chi[.(paste("df =", DF))]^2 ~ Distribution),
bquote(Critical~chi[alpha~.(paste(" =", alpha))]^2 == .(round(TheoCrit, nd))),
bquote(.(paste("Permuted Critical Value =", round(Crit, nd)))),
bquote(Observed~chi^2 == .(round(ChiSq, nd))),
expression(paste("")),
bquote(Permuted~italic(p)~.(pVal)))
}
}
H <- do.call(hist, c(histArgs["x"], plot = FALSE))
histArgs$ylim <- c(0, max(H$density, theoDist))
if (printLegend) {
legendArgs <- c(legendArgs, list(lty = c(2, 2, 1, 1, 0, 0),
lwd = c(2, 2, 2, 3, 0, 0),
col = c("black","black","black","red","","")))
}
} else { ################################################### other AFIs
badness <- grepl(pattern = "fmin|aic|bic|rmr|rmsea|cn|sic|hqc",
x = AFI, ignore.case = TRUE)
if (badness) {
Crit <- quantile(histArgs$x, probs = 1 - alpha)
} else {
Crit <- quantile(histArgs$x, probs = alpha)
}
histArgs$xlim <- range(histArgs$x, x@AFI.obs[AFI])
if (delta) {
histArgs$main <- bquote(~Permutation~Distribution~of~Delta*.(toupper(AFI)))
histArgs$xlab <- bquote(~Delta*.(toupper(AFI)))
if (printLegend) {
legendArgs$legend <- c(bquote(Critical~Delta*.(toupper(AFI))[alpha~.(paste(" =", alpha))] == .(round(Crit, nd))),
bquote(Observed~Delta*.(toupper(AFI)) == .(round(x@AFI.obs[AFI], nd))),
expression(paste("")),
bquote(Permuted~italic(p)~.(pVal)))
}
} else {
histArgs$main <- paste("Permutation Distribution of", toupper(AFI))
histArgs$xlab <- toupper(AFI)
if (printLegend) {
legendArgs$legend <- c(bquote(Critical~.(toupper(AFI))[alpha~.(paste(" =", alpha))] == .(round(Crit, nd))),
bquote(Observed~.(toupper(AFI)) == .(round(x@AFI.obs[AFI], nd))),
expression(paste("")),
bquote(Permuted~italic(p)~.(pVal)))
}
}
if (printLegend) {
legendArgs <- c(legendArgs, list(lty = c(1, 1, 0, 0),
lwd = c(2, 3, 0, 0),
col = c("black","red","","")))
}
}
## print histogram (and optionally, print legend)
suppressWarnings({
do.call(hist, histArgs)
if (grepl("chi", AFI)) {
lines(x = xVals, y = theoDist, lwd = 2, lty = 2)
abline(v = TheoCrit, col = "black", lwd = 2, lty = 2)
}
abline(v = Crit, col = "black", lwd = 2)
abline(v = x@AFI.obs[AFI], col = "red", lwd = 3)
if (printLegend) do.call(legend, legendArgs)
})
## return arguments to create histogram (and optionally, legend)
invisible(list(hist = histArgs, legend = legendArgs))
})
## --------------------
## Constructor Function
## --------------------
##' Permutation Randomization Tests of Measurement Equivalence and Differential
##' Item Functioning (DIF)
##'
##' The function `permuteMeasEq` provides tests of hypotheses involving
##' measurement equivalence, in one of two frameworks: multigroup CFA or MIMIC
##' models.
##'
##'
##' The function `permuteMeasEq` provides tests of hypotheses involving
##' measurement equivalence, in one of two frameworks:
##' \enumerate{
##' \item{1} For multiple-group CFA models, provide a pair of nested lavaan objects,
##' the less constrained of which (`uncon`) freely estimates a set of
##' measurement parameters (e.g., factor loadings, intercepts, or thresholds;
##' specified in `param`) in all groups, and the more constrained of which
##' (`con`) constrains those measurement parameters to equality across
##' groups. Group assignment is repeatedly permuted and the models are fit to
##' each permutation, in order to produce an empirical distribution under the
##' null hypothesis of no group differences, both for (a) changes in
##' user-specified fit measures (see `AFIs` and `moreAFIs`) and for
##' (b) the maximum modification index among the user-specified equality
##' constraints. Configural invariance can also be tested by providing that
##' fitted lavaan object to `con` and leaving `uncon = NULL`, in which
##' case `param` must be `NULL` as well.
##'
##' \item{2} In MIMIC models, one or a set of continuous and/or discrete
##' `covariates` can be permuted, and a constrained model is fit to each
##' permutation in order to provide a distribution of any fit measures (namely,
##' the maximum modification index among fixed parameters in `param`) under
##' the null hypothesis of measurement equivalence across levels of those
##' covariates.
##' }
##'
##' In either framework, modification indices for equality constraints or fixed
##' parameters specified in `param` are calculated from the constrained
##' model (`con`) using the function [lavaan::lavTestScore()].
##'
##' For multiple-group CFA models, the multiparameter omnibus null hypothesis of
##' measurement equivalence/invariance is that there are no group differences in
##' any measurement parameters (of a particular type). This can be tested using
##' the `anova` method on nested `lavaan` objects, as seen in the
##' output of [semTools::measurementInvariance()], or by inspecting
##' the change in alternative fit indices (AFIs) such as the CFI. The
##' permutation randomization method employed by `permuteMeasEq` generates
##' an empirical distribution of any `AFIs` under the null hypothesis, so
##' the user is not restricted to using fixed cutoffs proposed by Cheung &
##' Rensvold (2002), Chen (2007), or Meade, Johnson, & Braddy (2008).
##'
##' If the multiparameter omnibus null hypothesis is rejected, partial
##' invariance can still be established by freeing invalid equality constraints,
##' as long as equality constraints are valid for at least two indicators per
##' factor. Modification indices can be calculated from the constrained model
##' (`con`), but multiple testing leads to inflation of Type I error rates.
##' The permutation randomization method employed by `permuteMeasEq`
##' creates a distribution of the maximum modification index if the null
##' hypothesis is true, which allows the user to control the familywise Type I
##' error rate in a manner similar to Tukey's *q* (studentized range)
##' distribution for the Honestly Significant Difference (HSD) post hoc test.
##'
##' For MIMIC models, DIF can be tested by comparing modification indices of
##' regression paths to the permutation distribution of the maximum modification
##' index, which controls the familywise Type I error rate. The MIMIC approach
##' could also be applied with multiple-group models, but the grouping variable
##' would not be permuted; rather, the covariates would be permuted separately
##' within each group to preserve between-group differences. So whether
##' parameters are constrained or unconstrained across groups, the MIMIC
##' approach is only for testing null hypotheses about the effects of
##' `covariates` on indicators, controlling for common factors.
##'
##' In either framework, [lavaan::lavaan()]'s `group.label`
##' argument is used to preserve the order of groups seen in `con` when
##' permuting the data.
##'
##'
##' @importFrom lavaan lavInspect parTable
##'
##' @param nPermute An integer indicating the number of random permutations used
##' to form empirical distributions under the null hypothesis.
##' @param modelType A character string indicating type of model employed:
##' multiple-group CFA (`"mgcfa"`) or MIMIC (`"mimic"`).
##' @param con The constrained `lavaan` object, in which the parameters
##' specified in `param` are constrained to equality across all groups when
##' `modelType = "mgcfa"`, or which regression paths are fixed to zero when
##' `modelType = "mimic"`. In the case of testing *configural*
##' invariance when `modelType = "mgcfa"`, `con` is the configural
##' model (implicitly, the unconstrained model is the saturated model, so use
##' the defaults `uncon = NULL` and `param = NULL`). When
##' `modelType = "mimic"`, `con` is the MIMIC model in which the
##' covariate predicts the latent construct(s) but no indicators (unless they
##' have already been identified as DIF items).
##' @param uncon Optional. The unconstrained `lavaan` object, in which the
##' parameters specified in `param` are freely estimated in all groups.
##' When `modelType = "mgcfa"`, only in the case of testing
##' *configural* invariance should `uncon = NULL`. When
##' `modelType = "mimic"`, any non-`NULL uncon` is silently set to
##' `NULL`.
##' @param null Optional. A `lavaan` object, in which an alternative null
##' model is fit (besides the default independence model specified by
##' `lavaan`) for the calculation of incremental fit indices. See Widamin &
##' Thompson (2003) for details. If `NULL`, `lavaan`'s default
##' independence model is used.
##' @param param An optional character vector or list of character vectors
##' indicating which parameters the user would test for DIF following a
##' rejection of the omnibus null hypothesis tested using
##' (`more`)`AFIs`. Note that `param` does not guarantee certain
##' parameters *are* constrained in `con`; that is for the user to
##' specify when fitting the model. If users have any "anchor items" that they
##' would never intend to free across groups (or levels of a covariate), these
##' should be excluded from `param`; exceptions to a type of parameter can
##' be specified in `freeParam`. When `modelType = "mgcfa"`,
##' `param` indicates which parameters of interest are constrained across
##' groups in `con` and are unconstrained in `uncon`. Parameter names
##' must match those returned by `names(coef(con))`, but omitting any
##' group-specific suffixes (e.g., `"f1~1"` rather than `"f1~1.g2"`)
##' or user-specified labels (that is, the parameter names must follow the rules
##' of lavaan's [lavaan::model.syntax()]). Alternatively (or
##' additionally), to test all constraints of a certain type (or multiple types)
##' of parameter in `con`, `param` may take any combination of the
##' following values: `"loadings"`, `"intercepts"`,
##' `"thresholds"`, `"residuals"`, `"residual.covariances"`,
##' `"means"`, `"lv.variances"`, and/or `"lv.covariances"`. When
##' `modelType = "mimic"`, `param` must be a vector of individual
##' parameters or a list of character strings to be passed one-at-a-time to
##' `lavaan::lavTestScore(object = con, add = param[i])`,
##' indicating which (sets of) regression paths fixed to zero in `con` that
##' the user would consider freeing (i.e., exclude anchor items). If
##' `modelType = "mimic"` and `param` is a list of character strings,
##' the multivariate test statistic will be saved for each list element instead
##' of 1-*df* modification indices for each individual parameter, and
##' `names(param)` will name the rows of the `MI.obs` slot (see
##' [permuteMeasEq-class]). Set `param = NULL` (default) to avoid
##' collecting modification indices for any follow-up tests.
##' @param freeParam An optional character vector, silently ignored when
##' `modelType = "mimic"`. If `param` includes a type of parameter
##' (e.g., `"loadings"`), `freeParam` indicates exceptions (i.e.,
##' anchor items) that the user would *not* intend to free across groups
##' and should therefore be ignored when calculating *p* values adjusted
##' for the number of follow-up tests. Parameter types that are already
##' unconstrained across groups in the fitted `con` model (i.e., a
##' *partial* invariance model) will automatically be ignored, so they do
##' not need to be specified in `freeParam`. Parameter names must match
##' those returned by `names(coef(con))`, but omitting any group-specific
##' suffixes (e.g., `"f1~1"` rather than `"f1~1.g2"`) or
##' user-specified labels (that is, the parameter names must follow the rules of
##' lavaan [lavaan::model.syntax()]).
##' @param covariates An optional character vector, only applicable when
##' `modelType = "mimic"`. The observed data are partitioned into columns
##' indicated by `covariates`, and the rows are permuted simultaneously for
##' the entire set before being merged with the remaining data. Thus, the
##' covariance structure is preserved among the covariates, which is necessary
##' when (e.g.) multiple dummy codes are used to represent a discrete covariate
##' or when covariates interact. If `covariates = NULL` when
##' `modelType = "mimic"`, the value of `covariates` is inferred by
##' searching `param` for predictors (i.e., variables appearing after the
##' "`~`" operator).
##' @param AFIs A character vector indicating which alternative fit indices (or
##' chi-squared itself) are to be used to test the multiparameter omnibus null
##' hypothesis that the constraints specified in `con` hold in the
##' population. Any fit measures returned by [lavaan::fitMeasures()]
##' may be specified (including constants like `"df"`, which would be
##' nonsensical). If both `AFIs` and `moreAFIs` are `NULL`, only
##' `"chisq"` will be returned.
##' @param moreAFIs Optional. A character vector indicating which (if any)
##' alternative fit indices returned by [semTools::moreFitIndices()]
##' are to be used to test the multiparameter omnibus null hypothesis that the
##' constraints specified in `con` hold in the population.
##' @param maxSparse Only applicable when `modelType = "mgcfa"` and at
##' least one indicator is `ordered`. An integer indicating the maximum
##' number of consecutive times that randomly permuted group assignment can
##' yield a sample in which at least one category (of an `ordered`
##' indicator) is unobserved in at least one group, such that the same set of
##' parameters cannot be estimated in each group. If such a sample occurs, group
##' assignment is randomly permuted again, repeatedly until a sample is obtained
##' with all categories observed in all groups. If `maxSparse` is exceeded,
##' `NA` will be returned for that iteration of the permutation
##' distribution.
##' @param maxNonconv An integer indicating the maximum number of consecutive
##' times that a random permutation can yield a sample for which the model does
##' not converge on a solution. If such a sample occurs, permutation is
##' attempted repeatedly until a sample is obtained for which the model does
##' converge. If `maxNonconv` is exceeded, `NA` will be returned for
##' that iteration of the permutation distribution, and a warning will be
##' printed when using `show` or `summary`.
##' @param showProgress Logical. Indicating whether to display a progress bar
##' while permuting. Silently set to `FALSE` when using parallel options.
##' @param warn Sets the handling of warning messages when fitting model(s) to
##' permuted data sets. See [base::options()].
##' @param datafun An optional function that can be applied to the data
##' (extracted from `con`) after each permutation, but before fitting the
##' model(s) to each permutation. The `datafun` function must have an
##' argument named `data` that accepts a `data.frame`, and it must
##' return a `data.frame` containing the same column names. The column
##' order may differ, the values of those columns may differ (so be careful!),
##' and any additional columns will be ignored when fitting the model, but an
##' error will result if any column names required by the model syntax do not
##' appear in the transformed data set. Although available for any
##' `modelType`, `datafun` may be useful when using the MIMIC method
##' to test for nonuniform DIF (metric/weak invariance) by using product
##' indicators for a latent factor representing the interaction between a factor
##' and one of the `covariates`, in which case the product indicators would
##' need to be recalculated after each permutation of the `covariates`. To
##' access other R objects used within `permuteMeasEq`, the arguments to
##' `datafun` may also contain any subset of the following: `"con"`,
##' `"uncon"`, `"null"`, `"param"`, `"freeParam"`,
##' `"covariates"`, `"AFIs"`, `"moreAFIs"`, `"maxSparse"`,
##' `"maxNonconv"`, and/or `"iseed"`. The values for those arguments
##' will be the same as the values supplied to `permuteMeasEq`.
##' @param extra An optional function that can be applied to any (or all) of the
##' fitted lavaan objects (`con`, `uncon`, and/or `null`). This
##' function will also be applied after fitting the model(s) to each permuted
##' data set. To access the R objects used within `permuteMeasEq`, the
##' arguments to `extra` must be any subset of the following: `"con"`,
##' `"uncon"`, `"null"`, `"param"`, `"freeParam"`,
##' `"covariates"`, `"AFIs"`, `"moreAFIs"`, `"maxSparse"`,
##' `"maxNonconv"`, and/or `"iseed"`. The values for those arguments
##' will be the same as the values supplied to `permuteMeasEq`. The
##' `extra` function must return a named `numeric` vector or a named
##' `list` of scalars (i.e., a `list` of `numeric` vectors of
##' `length == 1`). Any unnamed elements (e.g., `""` or `NULL`)
##' of the returned object will result in an error.
##' @param parallelType The type of parallel operation to be used (if any). The
##' default is `"none"`. Forking is not possible on Windows, so if
##' `"multicore"` is requested on a Windows machine, the request will be
##' changed to `"snow"` with a message.
##' @param ncpus Integer: number of processes to be used in parallel operation.
##' If `NULL` (the default) and `parallelType %in%
##' c("multicore","snow")`, the default is one less than the maximum number of
##' processors detected by [parallel::detectCores()]. This default is
##' also silently set if the user specifies more than the number of processors
##' detected.
##' @param cl An optional \pkg{parallel} or \pkg{snow} cluster for use when
##' `parallelType = "snow"`. If `NULL`, a `"PSOCK"` cluster on
##' the local machine is created for the duration of the `permuteMeasEq`
##' call. If a valid [parallel::makeCluster()] object is supplied,
##' `parallelType` is silently set to `"snow"`, and `ncpus` is
##' silently set to `length(cl)`.
##' @param iseed Integer: Only used to set the states of the RNG when using
##' parallel options, in which case [base::RNGkind()] is set to
##' `"L'Ecuyer-CMRG"` with a message. See
##' [parallel::clusterSetRNGStream()] and Section 6 of
##' `vignette("parallel", "parallel")` for more details. If user supplies
##' an invalid value, `iseed` is silently set to the default (12345). To
##' set the state of the RNG when not using parallel options, call
##' [base::set.seed()] before calling `permuteMeasEq`.
##'
##' @return The [permuteMeasEq-class] object representing the results of
##' testing measurement equivalence (the multiparameter omnibus test) and DIF
##' (modification indices), as well as diagnostics and any `extra` output.
##'
##' @author Terrence D. Jorgensen (University of Amsterdam;
##' \email{TJorgensen314@@gmail.com})
##'
##' @seealso [stats::TukeyHSD()], [lavaan::lavTestScore()],
##' [semTools::measurementInvariance()],
##' [semTools::measurementInvarianceCat()]
##'
##' @references
##'
##' **Papers about permutation tests of measurement equivalence:**
##'
##' Jorgensen, T. D., Kite, B. A., Chen, P.-Y., & Short, S. D. (2018).
##' Permutation randomization methods for testing measurement equivalence and
##' detecting differential item functioning in multiple-group confirmatory
##' factor analysis. *Psychological Methods, 23*(4), 708--728.
##' \doi{10.1037/met0000152}
##'
##' Kite, B. A., Jorgensen, T. D., & Chen, P.-Y. (2018). Random permutation
##' testing applied to measurement invariance testing with ordered-categorical
##' indicators. *Structural Equation Modeling 25*(4), 573--587.
##' \doi{10.1080/10705511.2017.1421467}
##'
##' Jorgensen, T. D. (2017). Applying permutation tests and multivariate
##' modification indices to configurally invariant models that need
##' respecification. *Frontiers in Psychology, 8*(1455).
##' \doi{10.3389/fpsyg.2017.01455}
##'
##' **Additional reading:**
##'
##' Chen, F. F. (2007). Sensitivity of goodness of fit indexes to
##' lack of measurement invariance. *Structural Equation Modeling, 14*(3),
##' 464--504. \doi{10.1080/10705510701301834}
##'
##' Cheung, G. W., & Rensvold, R. B. (2002). Evaluating goodness-of-fit indexes
##' for testing measurement invariance. *Structural Equation Modeling,
##' 9*(2), 233--255. \doi{10.1207/S15328007SEM0902_5}
##'
##' Meade, A. W., Johnson, E. C., & Braddy, P. W. (2008). Power and sensitivity
##' of alternative fit indices in tests of measurement invariance. *Journal
##' of Applied Psychology, 93*(3), 568--592. \doi{10.1037/0021-9010.93.3.568}
##'
##' Widamin, K. F., & Thompson, J. S. (2003). On specifying the null model for
##' incremental fit indices in structural equation modeling. *Psychological
##' Methods, 8*(1), 16--37. \doi{10.1037/1082-989X.8.1.16}
##'
##' @examples
##'
##' \donttest{
##'
##' ########################
##' ## Multiple-Group CFA ##
##' ########################
##'
##' ## create 3-group data in lavaan example(cfa) data
##' HS <- lavaan::HolzingerSwineford1939
##' HS$ageGroup <- ifelse(HS$ageyr < 13, "preteen",
##' ifelse(HS$ageyr > 13, "teen", "thirteen"))
##'
##' ## specify and fit an appropriate null model for incremental fit indices
##' mod.null <- c(paste0("x", 1:9, " ~ c(T", 1:9, ", T", 1:9, ", T", 1:9, ")*1"),
##' paste0("x", 1:9, " ~~ c(L", 1:9, ", L", 1:9, ", L", 1:9, ")*x", 1:9))
##' fit.null <- cfa(mod.null, data = HS, group = "ageGroup")
##'
##' ## fit target model with varying levels of measurement equivalence
##' mod.config <- '
##' visual =~ x1 + x2 + x3
##' textual =~ x4 + x5 + x6
##' speed =~ x7 + x8 + x9
##' '
##' fit.config <- cfa(mod.config, data = HS, std.lv = TRUE, group = "ageGroup")
##' fit.metric <- cfa(mod.config, data = HS, std.lv = TRUE, group = "ageGroup",
##' group.equal = "loadings")
##' fit.scalar <- cfa(mod.config, data = HS, std.lv = TRUE, group = "ageGroup",
##' group.equal = c("loadings","intercepts"))
##'
##'
##' ####################### Permutation Method
##'
##' ## fit indices of interest for multiparameter omnibus test
##' myAFIs <- c("chisq","cfi","rmsea","mfi","aic")
##' moreAFIs <- c("gammaHat","adjGammaHat")
##'
##' ## Use only 20 permutations for a demo. In practice,
##' ## use > 1000 to reduce sampling variability of estimated p values
##'
##' ## test configural invariance
##' set.seed(12345)
##' out.config <- permuteMeasEq(nPermute = 20, con = fit.config)
##' out.config
##'
##' ## test metric equivalence
##' set.seed(12345) # same permutations
##' out.metric <- permuteMeasEq(nPermute = 20, uncon = fit.config, con = fit.metric,
##' param = "loadings", AFIs = myAFIs,
##' moreAFIs = moreAFIs, null = fit.null)
##' summary(out.metric, nd = 4)
##'
##' ## test scalar equivalence
##' set.seed(12345) # same permutations
##' out.scalar <- permuteMeasEq(nPermute = 20, uncon = fit.metric, con = fit.scalar,
##' param = "intercepts", AFIs = myAFIs,
##' moreAFIs = moreAFIs, null = fit.null)
##' summary(out.scalar)
##'
##' ## Not much to see without significant DIF.
##' ## Try using an absurdly high alpha level for illustration.
##' outsum <- summary(out.scalar, alpha = .50)
##'
##' ## notice that the returned object is the table of DIF tests
##' outsum
##'
##' ## visualize permutation distribution
##' hist(out.config, AFI = "chisq")
##' hist(out.metric, AFI = "chisq", nd = 2, alpha = .01,
##' legendArgs = list(x = "topright"))
##' hist(out.scalar, AFI = "cfi", printLegend = FALSE)
##'
##'
##' ####################### Extra Output
##'
##' ## function to calculate expected change of Group-2 and -3 latent means if
##' ## each intercept constraint were released
##' extra <- function(con) {
##' output <- list()
##' output["x1.vis2"] <- lavTestScore(con, release = 19:20, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[70]
##' output["x1.vis3"] <- lavTestScore(con, release = 19:20, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[106]
##' output["x2.vis2"] <- lavTestScore(con, release = 21:22, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[70]
##' output["x2.vis3"] <- lavTestScore(con, release = 21:22, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[106]
##' output["x3.vis2"] <- lavTestScore(con, release = 23:24, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[70]
##' output["x3.vis3"] <- lavTestScore(con, release = 23:24, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[106]
##' output["x4.txt2"] <- lavTestScore(con, release = 25:26, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[71]
##' output["x4.txt3"] <- lavTestScore(con, release = 25:26, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[107]
##' output["x5.txt2"] <- lavTestScore(con, release = 27:28, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[71]
##' output["x5.txt3"] <- lavTestScore(con, release = 27:28, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[107]
##' output["x6.txt2"] <- lavTestScore(con, release = 29:30, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[71]
##' output["x6.txt3"] <- lavTestScore(con, release = 29:30, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[107]
##' output["x7.spd2"] <- lavTestScore(con, release = 31:32, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[72]
##' output["x7.spd3"] <- lavTestScore(con, release = 31:32, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[108]
##' output["x8.spd2"] <- lavTestScore(con, release = 33:34, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[72]
##' output["x8.spd3"] <- lavTestScore(con, release = 33:34, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[108]
##' output["x9.spd2"] <- lavTestScore(con, release = 35:36, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[72]
##' output["x9.spd3"] <- lavTestScore(con, release = 35:36, univariate = FALSE,
##' epc = TRUE, warn = FALSE)$epc$epc[108]
##' output
##' }
##'
##' ## observed EPC
##' extra(fit.scalar)
##'
##' ## permutation results, including extra output
##' set.seed(12345) # same permutations
##' out.scalar <- permuteMeasEq(nPermute = 20, uncon = fit.metric, con = fit.scalar,
##' param = "intercepts", AFIs = myAFIs,
##' moreAFIs = moreAFIs, null = fit.null, extra = extra)
##' ## summarize extra output
##' summary(out.scalar, extra = TRUE)
##'
##'
##' ###########
##' ## MIMIC ##
##' ###########
##'
##' ## Specify Restricted Factor Analysis (RFA) model, equivalent to MIMIC, but
##' ## the factor covaries with the covariate instead of being regressed on it.
##' ## The covariate defines a single-indicator construct, and the
##' ## double-mean-centered products of the indicators define a latent
##' ## interaction between the factor and the covariate.
##' mod.mimic <- '
##' visual =~ x1 + x2 + x3
##' age =~ ageyr
##' age.by.vis =~ x1.ageyr + x2.ageyr + x3.ageyr
##'
##' x1 ~~ x1.ageyr
##' x2 ~~ x2.ageyr
##' x3 ~~ x3.ageyr
##' '
##'
##' HS.orth <- indProd(var1 = paste0("x", 1:3), var2 = "ageyr", match = FALSE,
##' data = HS[ , c("ageyr", paste0("x", 1:3))] )
##' fit.mimic <- cfa(mod.mimic, data = HS.orth, meanstructure = TRUE)
##' summary(fit.mimic, stand = TRUE)
##'
##' ## Whereas MIMIC models specify direct effects of the covariate on an indicator,
##' ## DIF can be tested in RFA models by specifying free loadings of an indicator
##' ## on the covariate's construct (uniform DIF, scalar invariance) and the
##' ## interaction construct (nonuniform DIF, metric invariance).
##' param <- as.list(paste0("age + age.by.vis =~ x", 1:3))
##' names(param) <- paste0("x", 1:3)
##' # param <- as.list(paste0("x", 1:3, " ~ age + age.by.vis")) # equivalent
##'
##' ## test both parameters simultaneously for each indicator
##' do.call(rbind, lapply(param, function(x) lavTestScore(fit.mimic, add = x)$test))
##' ## or test each parameter individually
##' lavTestScore(fit.mimic, add = as.character(param))
##'
##'
##' ####################### Permutation Method
##'
##' ## function to recalculate interaction terms after permuting the covariate
##' datafun <- function(data) {
##' d <- data[, c(paste0("x", 1:3), "ageyr")]
##' indProd(var1 = paste0("x", 1:3), var2 = "ageyr", match = FALSE, data = d)
##' }
##'
##' set.seed(12345)
##' perm.mimic <- permuteMeasEq(nPermute = 20, modelType = "mimic",
##' con = fit.mimic, param = param,
##' covariates = "ageyr", datafun = datafun)
##' summary(perm.mimic)
##'
##' }
##'
##' @export
permuteMeasEq <- function(nPermute, modelType = c("mgcfa","mimic"),
con, uncon = NULL, null = NULL,
param = NULL, freeParam = NULL, covariates = NULL,
AFIs = NULL, moreAFIs = NULL,
maxSparse = 10, maxNonconv = 10, showProgress = TRUE,
warn = -1, datafun, extra,
parallelType = c("none","multicore","snow"),
ncpus = NULL, cl = NULL, iseed = 12345) {
## save arguments from call
availableArgs <- as.list(formals(permuteMeasEq))
argNames <- names(availableArgs)
if (missing(datafun)) argNames <- setdiff(argNames, "datafun")
if (missing(extra)) argNames <- setdiff(argNames, "extra")
for (aa in argNames) {
if (!is.null(eval(as.name(aa))))
suppressWarnings(availableArgs[[aa]] <- eval(as.name(aa)))
}
## check and return them
fullCall <- do.call(checkPermArgs, availableArgs)
## assign them to workspace (also adds old_RNG & oldSeed to workspace)
for (aa in names(fullCall)) assign(aa, fullCall[[aa]])
###################### SAVE OBSERVED RESULTS ##########################
AFI.obs <- do.call(getAFIs, fullCall)
## save modification indices if !is.null(param)
if (is.null(param)) {
MI.obs <- data.frame(NULL)
} else MI.obs <- do.call(getMIs, fullCall)
## anything extra?
if (!missing(extra)) {
extraArgs <- formals(extra)
neededArgs <- intersect(names(extraArgs), names(fullCall))
extraArgs <- do.call(c, lapply(neededArgs, function(nn) fullCall[nn]))
extraOut <- do.call(extra, extraArgs)
## check that extra() returns a named list of scalars
if (!is.list(extraOut)) extraOut <- as.list(extraOut)
wrongFormat <- paste('Function "extra" must return a numeric vector or a',
'list of scalars, with each element named.')
if (!all(sapply(extraOut, is.numeric))) stop(wrongFormat)
if (!all(sapply(extraOut, length) == 1L)) stop(wrongFormat)
if (is.null(names(extraOut)) | any(names(extraOut) == "")) stop(wrongFormat)
extra.obs <- do.call(c, extraOut)
} else extra.obs <- numeric(length = 0L)
######################### PREP DATA ##############################
argList <- fullCall[c("con","uncon","null","param","freeParam","covariates",
"AFIs","moreAFIs","maxSparse","maxNonconv","warn","iseed")]
argList$G <- lavInspect(con, "group")
## check for categorical variables
# catVars <- lavaan::lavNames(con, type = "ov.ord")
# numVars <- lavaan::lavNames(con, type = "ov.num")
# latentVars <- lavaan::lavNames(con, type = "lv.regular")
## assemble data to which the models were fit
if (length(argList$G)) {
dataList <- mapply(FUN = function(x, g, n) {
y <- data.frame(as.data.frame(x), g, stringsAsFactors = FALSE)
names(y) <- c(n, argList$G)
y
}, SIMPLIFY = FALSE,
x = lavInspect(con, "data"), g = lavInspect(con, "group.label"),
n = lavaan::lavNames(con, type = "ov",
group = seq_along(lavInspect(con, "group.label"))))
argList$d <- do.call(rbind, dataList)
} else {
argList$d <- as.data.frame(lavInspect(con, "data"))
names(argList$d) <- lavaan::lavNames(con, type = "ov")
}
## check that covariates are actual variables
if (modelType == "mimic") {
if (length(covariates) && !all(covariates %in% names(argList$d)))
stop('These specified covariates are not columns in the data.frame:\n',
paste(setdiff(covariates, names(argList$d)), collapse = ", "))
}
## anything extra?
if (!missing(extra)) argList$extra <- extra
if (!missing(datafun)) argList$datafun <- datafun
###################### PERMUTED RESULTS ###########################
## permute and return distributions of (delta)AFIs, largest MI, and extras
if (showProgress) {
mypb <- utils::txtProgressBar(min = 1, max = nPermute, initial = 1,
char = "=", width = 50, style = 3, file = "")
permuDist <- list()
for (j in 1:nPermute) {
permuDist[[j]] <- do.call(paste("permuteOnce", modelType, sep = "."),
args = c(argList, i = j))
utils::setTxtProgressBar(mypb, j)
}
close(mypb)
} else if (parallelType == "multicore") {
if (length(iseed)) set.seed(iseed)
argList$FUN <- as.name(paste("permuteOnce", modelType, sep = "."))
argList$X <- 1:nPermute
argList$mc.cores <- ncpus
argList$mc.set.seed <- TRUE
pmcl <- function(...) { parallel::mclapply(...) }
permuDist <- do.call(pmcl, args = argList)
## restore old RNG type
if (fullCall$old_RNG[1] != "L'Ecuyer-CMRG") RNGkind(fullCall$old_RNG[1])
} else if (parallelType == "snow") {
stopTheCluster <- FALSE
if (is.null(cl)) {
stopTheCluster <- TRUE
cl <- parallel::makePSOCKcluster(rep("localhost", ncpus))
}
parallel::clusterSetRNGStream(cl, iseed = iseed)
argList$cl <- cl
argList$X <- 1:nPermute
argList$fun <- paste("permuteOnce", modelType, sep = ".")
parallel::clusterExport(cl, varlist = c(argList$fun, "getAFIs","getMIs")) #FIXME: need update?
tempppl <- function(...) { parallel::parLapply(...) }
permuDist <- do.call(tempppl, args = argList)
if (stopTheCluster) parallel::stopCluster(cl)
## restore old RNG type
if (fullCall$old_RNG[1] != "L'Ecuyer-CMRG") RNGkind(fullCall$old_RNG[1])
} else {
argList$X <- 1:nPermute
argList$FUN <- paste("permuteOnce", modelType, sep = ".")
permuDist <- do.call(lapply, args = argList)
}
## extract AFI distribution
if (length(AFI.obs) > 1) {
AFI.dist <- as.data.frame(t(sapply(permuDist, function(x) x$AFI)))
}
if (length(AFI.obs) == 1L) {
AFI.dist <- data.frame(sapply(permuDist, function(x) x$AFI))
colnames(AFI.dist) <- names(AFI.obs)
}
## identify badness-of-fit measures
badness <- grepl(pattern = "fmin|chi|aic|bic|rmr|rmsea|cn|sic|hqc",
x = names(AFI.obs), ignore.case = TRUE)
## calculate all one-directional p-values
AFI.pval <- mapply(FUN = function(x, y, b) {
if (b) return(mean(x >= y, na.rm = TRUE))
mean(x <= y, na.rm = TRUE)
}, x = unclass(AFI.dist), y = AFI.obs, b = badness)
## extract distribution of maximum modification indices
MI.dist <- as.numeric(unlist(lapply(permuDist, function(x) x$MI)))
## calculate Tukey-adjusted p values for modification indices
if (!is.null(param)) {
MI.obs$tukey.p.value <- sapply(MI.obs$X2,
function(i) mean(i <= MI.dist, na.rm = TRUE))
MI.obs <- as.data.frame(unclass(MI.obs))
rownames(MI.obs) <- names(param)
}
## anything extra?
if (!missing(extra)) {
extra.dist <- do.call(rbind, lapply(permuDist, function(x) x$extra))
} else extra.dist <- data.frame(NULL)
## save parameter table for show/summary methods
PT <- as.data.frame(parTable(con))
PT$par <- paste0(PT$lhs, PT$op, PT$rhs)
if (length(lavInspect(con, "group")))
PT$group.label[PT$group > 0] <- lavInspect(con, "group.label")[PT$group[PT$group > 0] ]
## return observed results, permutation p values, and ANOVA results
if (is.null(uncon)) {
delta <- lavaan::anova(con)
} else {
delta <- lavaan::anova(uncon, con)
}
ANOVA <- sapply(delta[,c("Chisq diff","Df diff","Pr(>Chisq)")], function(x) x[2])
out <- new("permuteMeasEq", PT = PT, modelType = modelType, ANOVA = ANOVA,
AFI.obs = AFI.obs, AFI.dist = AFI.dist, AFI.pval = AFI.pval,
MI.obs = MI.obs, MI.dist = MI.dist,
extra.obs = extra.obs, extra.dist = extra.dist,
n.Permutations = nPermute, n.Converged = sum(!is.na(AFI.dist[,1])),
n.nonConverged = sapply(permuDist, function(x) x$n.nonConverged),
n.Sparse = sapply(permuDist, function(x) x$n.Sparse),
oldSeed = fullCall$oldSeed)
out
}
## ----------------
## Hidden Functions
## ----------------
## function to check validity of arguments to permuteMeasEq()
#' @importFrom lavaan lavInspect parTable
checkPermArgs <- function(nPermute, modelType, con, uncon, null,
param, freeParam, covariates, AFIs, moreAFIs,
maxSparse, maxNonconv, showProgress, warn,
datafun, extra, parallelType, ncpus, cl, iseed) {
fixedCall <- as.list(match.call())[-1]
fixedCall$nPermute <- as.integer(nPermute[1])
fixedCall$modelType <- modelType[1]
if (!fixedCall$modelType %in% c("mgcfa","mimic","long"))
stop('modelType must be one of c("mgcfa","mimic","long")')
if (fixedCall$modelType == "long") stop('modelType "long" is not yet available.')
if (fixedCall$modelType == "mgcfa" && lavInspect(con, "ngroups") == 1L)
stop('modelType = "mgcfa" applies only to multigroup models.')
if (fixedCall$modelType == "mimic") {
uncon <- NULL
fixedCall$uncon <- NULL
fixedCall <- c(fixedCall, list(uncon = NULL))
}
## strip white space
if (is.list(param)) {
fixedCall$param <- lapply(param, function(cc) gsub("[[:space:]]+", "", cc))
} else if (!is.null(param)) fixedCall$param <- gsub("[[:space:]]+", "", param)
if (!is.null(freeParam)) fixedCall$freeParam <- gsub("[[:space:]]+", "", freeParam)
if (fixedCall$modelType == "mimic") {
# PT <- lavaan::lavaanify(fixedCall$param)
# checkCovs <- unique(PT$rhs[PT$op == "~"])
# if (is.null(covariates)) covariates <- checkCovs
# if (length(setdiff(covariates, checkCovs)))
# warning('Argument "covariates" includes predictors not in argument "param"')
##### ordVars <- lavaan::lavNames(con, type = "ov.ord")
fixedCall$covariates <- as.character(covariates)
}
fixedCall$maxSparse <- as.integer(maxSparse[1])
fixedCall$maxNonconv <- as.integer(maxNonconv[1])
fixedCall$showProgress <- as.logical(showProgress[1])
fixedCall$warn <- as.integer(warn[1])
fixedCall$oldSeed <- as.integer(NULL)
parallelType <- as.character(parallelType[1])
if (!parallelType %in% c("none","multicore","snow")) parallelType <- "none"
if (!is.null(cl)) {
if (!is(cl, "cluster")) stop("Invalid cluster object. Check class(cl)")
parallelType <- "snow"
ncpus <- length(cl)
}
if (parallelType == "multicore" && .Platform$OS.type == "windows") {
parallelType <- "snow"
message("'multicore' option unavailable on Windows. Using 'snow' instead.")
}
## parallel settings, adapted from boot::boot()
if (parallelType != "none") {
if (is.null(ncpus) || ncpus > parallel::detectCores()) {
ncpus <- parallel::detectCores() - 1
}
if (ncpus <= 1L) {
parallelType <- "none"
} else {
fixedCall$showProgress <- FALSE
fixedCall$old_RNG <- RNGkind()
fixedCall$oldSeed <- .Random.seed
if (fixedCall$old_RNG[1] != "L'Ecuyer-CMRG") {
RNGkind("L'Ecuyer-CMRG")
message("Your RNGkind() was changed from ", fixedCall$old_RNG[1],
" to L'Ecuyer-CMRG, which is required for reproducibility ",
" in parallel jobs. Your RNGkind() has been returned to ",
fixedCall$old_RNG[1], " but the seed has not been set. ",
" The state of your previous RNG is saved in the slot ",
" named 'oldSeed', if you want to restore it using ",
" the syntax:\n",
".Random.seed[-1] <- permuteMeasEqObject@oldSeed[-1]")
}
fixedCall$iseed <- as.integer(iseed[1])
if (is.na(fixedCall$iseed)) fixedCall$iseed <- 12345
}
}
fixedCall$parallelType <- parallelType
if (is.null(ncpus)) {
fixedCall$ncpus <- NULL
fixedCall <- c(fixedCall, list(ncpus = NULL))
} else fixedCall$ncpus <- ncpus
## Check that "param" is NULL if uncon is NULL, and check for lavaan class.
## Also check that models are fitted to raw data, not summary stats.
notLavaan <- "Non-NULL 'con', 'uncon', or 'null' must be fitted lavaan object."
notRawData <- "lavaan models ('con', 'uncon', or 'null') must be fitted to raw data=, not summary statistics (e.g., sample.cov=)"
if (is.null(uncon)) {
if (!is.null(fixedCall$param) && fixedCall$modelType == "mgcfa") {
message(c(" When 'uncon = NULL', only configural invariance is tested.",
"\n So the 'param' argument was changed to NULL."))
fixedCall$param <- NULL
fixedCall <- c(fixedCall, list(param = NULL))
}
if (!inherits(con, "lavaan")) stop(notLavaan)
stopifnot(con@Data@data.type == "full")
} else {
if (!inherits(con, "lavaan")) stop(notLavaan)
if (!inherits(uncon, "lavaan")) stop(notLavaan)
stopifnot( con@Data@data.type == "full")
stopifnot(uncon@Data@data.type == "full")
}
if (!is.null(null)) {
if (!inherits(null, "lavaan")) stop(notLavaan)
stopifnot(null@Data@data.type == "full")
}
############ FIXME: check that lavInspect(con, "options")$conditional.x = FALSE (find defaults for continuous/ordered indicators)
if (!is.null(fixedCall$param)) {
## Temporarily warn about testing thresholds without necessary constraints. FIXME: check for binary indicators
if ("thresholds" %in% fixedCall$param | any(grepl("\\|", fixedCall$param))) {
warning(c("This function is not yet optimized for testing thresholds.\n",
"Necessary identification contraints might not be specified."))
}
## collect parameter types for "mgcfa"
if (fixedCall$modelType != "mimic") {
## save all estimates from constrained model
PT <- parTable(con)[ , c("lhs","op","rhs","group","plabel")]
## extract parameters of interest
paramTypes <- c("loadings","intercepts","thresholds","residuals","means",
"residual.covariances","lv.variances","lv.covariances")
params <- PT[paste0(PT$lhs, PT$op, PT$rhs) %in% setdiff(fixedCall$param,
paramTypes), ]
## add parameters by type, if any are specified
types <- intersect(fixedCall$param, paramTypes)
ov.names <- lavaan::lavNames(con, "ov")
isOV <- PT$lhs %in% ov.names
lv.names <- con@pta$vnames$lv[[1]]
isLV <- PT$lhs %in% lv.names & PT$rhs %in% lv.names
if ("loadings" %in% types) params <- rbind(params, PT[PT$op == "=~", ])
if ("intercepts" %in% types) {
params <- rbind(params, PT[isOV & PT$op == "~1", ])
}
if ("thresholds" %in% types) params <- rbind(params, PT[PT$op == "|", ])
if ("residuals" %in% types) {
params <- rbind(params, PT[isOV & PT$lhs == PT$rhs & PT$op == "~~", ])
}
if ("residual.covariances" %in% types) {
params <- rbind(params, PT[isOV & PT$lhs != PT$rhs & PT$op == "~~", ])
}
if ("means" %in% types) {
params <- rbind(params, PT[PT$lhs %in% lv.names & PT$op == "~1", ])
}
if ("lv.variances" %in% types) {
params <- rbind(params, PT[isLV & PT$lhs == PT$rhs & PT$op == "~~", ])
}
if ("lv.covariances" %in% types) {
params <- rbind(params, PT[isLV & PT$lhs != PT$rhs & PT$op == "~~", ])
}
## remove parameters specified by "freeParam" argument
params <- params[!paste0(params$lhs, params$op, params$rhs) %in% fixedCall$freeParam, ]
fixedCall$param <- paste0(params$lhs, params$op, params$rhs)
}
}
if (is.null(AFIs) & is.null(moreAFIs)) {
message("No AFIs were selected, so only chi-squared will be permuted.\n")
fixedCall$AFIs <- "chisq"
AFIs <- "chisq"
}
if ("ecvi" %in% AFIs & lavInspect(con, "ngroups") > 1L)
stop("ECVI is not available for multigroup models.")
## check estimators
leastSq <- grepl("LS", lavInspect(con, "options")$estimator)
if (!is.null(uncon)) {
if (uncon@Options$estimator != lavInspect(con, "options")$estimator)
stop("Models must be fit using same estimator.")
}
if (!is.null(null)) {
if (lavInspect(null, "options")$estimator != lavInspect(con, "options")$estimator)
stop("Models must be fit using same estimator.")
}
## check extra functions, if any
restrictedArgs <- c("con","uncon","null","param","freeParam","covariates",
"AFIs","moreAFIs","maxSparse","maxNonconv","iseed")
if (!missing(datafun)) {
if (!is.function(datafun)) stop('Argument "datafun" must be a function.')
extraArgs <- formals(datafun)
if (!all(names(extraArgs) %in% c(restrictedArgs, "data")))
stop('The user-supplied function "datafun" can only have any among the ',
'following arguments:\n', paste(restrictedArgs, collapse = ", "))
}
if (!missing(extra)) {
if (!is.function(extra)) stop('Argument "extra" must be a function.')
extraArgs <- formals(extra)
if (!all(names(extraArgs) %in% restrictedArgs))
stop('The user-supplied function "extra" can only have any among the ',
'following arguments:\n', paste(restrictedArgs, collapse = ", "))
}
## return evaluated list of other arguments
lapply(fixedCall, eval)
}
## function to extract fit measures
#' @importFrom lavaan lavInspect
getAFIs <- function(...) {
dots <- list(...)
AFI1 <- list()
AFI0 <- list()
leastSq <- grepl("LS", lavInspect(dots$con, "options")$estimator)
## check validity of user-specified AFIs, save output
if (!is.null(dots$AFIs)) {
IC <- grep("ic|logl", dots$AFIs, value = TRUE)
if (leastSq & length(IC)) {
stop(paste("Argument 'AFIs' includes invalid options:",
paste(IC, collapse = ", "),
"Information criteria unavailable for least-squares estimators.",
sep = "\n"))
}
if (!is.null(dots$uncon))
AFI1[[1]] <- lavaan::fitMeasures(dots$uncon, fit.measures = dots$AFIs,
baseline.model = dots$null)
AFI0[[1]] <- lavaan::fitMeasures(dots$con, fit.measures = dots$AFIs,
baseline.model = dots$null)
}
## check validity of user-specified moreAFIs
if (!is.null(dots$moreAFIs)) {
IC <- grep("ic|hqc", dots$moreAFIs, value = TRUE)
if (leastSq & length(IC)) {
stop(paste("Argument 'moreAFIs' includes invalid options:",
paste(IC, collapse = ", "),
"Information criteria unavailable for least-squares estimators.",
sep = "\n"))
}
if (!is.null(dots$uncon))
AFI1[[2]] <- moreFitIndices(dots$uncon, fit.measures = dots$moreAFIs)
AFI0[[2]] <- moreFitIndices(dots$con, fit.measures = dots$moreAFIs)
}
## save observed AFIs or delta-AFIs
if (is.null(dots$uncon)) {
AFI.obs <- unlist(AFI0)
} else {
AFI.obs <- unlist(AFI0) - unlist(AFI1)
}
AFI.obs
}
## Function to extract modification indices for equality constraints
#' @importFrom lavaan parTable
getMIs <- function(...) {
dots <- list(...)
if (dots$modelType == "mgcfa") {
## save all estimates from constrained model
PT <- parTable(dots$con)[ , c("lhs","op","rhs","group","plabel")]
## extract parameters of interest
params <- PT[paste0(PT$lhs, PT$op, PT$rhs) %in% dots$param, ]
## return modification indices for specified constraints (param)
MIs <- lavaan::lavTestScore(dots$con)$uni
MI.obs <- MIs[MIs$lhs %in% params$plabel, ]
} else if (dots$modelType == "mimic") {
if (is.list(dots$param)) {
MI <- lapply(dots$param, function(x) lavaan::lavTestScore(dots$con, add = x)$test)
MI.obs <- do.call(rbind, MI)
} else MI.obs <- lavaan::lavTestScore(dots$con, add = dots$param)$uni
} else if (dots$modelType == "long") {
## coming soon
}
MI.obs
}
## Functions to find delta-AFIs & maximum modification index in one permutation
permuteOnce.mgcfa <- function(i, d, G, con, uncon, null, param, freeParam,
covariates, AFIs, moreAFIs, maxSparse, maxNonconv,
iseed, warn, extra = NULL, datafun = NULL) {
old_warn <- options()$warn
options(warn = warn)
## save arguments from call
argNames <- names(formals(permuteOnce.mgcfa))
availableArgs <- lapply(argNames, function(x) eval(as.name(x)))
names(availableArgs) <- argNames
group.label <- lavaan::lavInspect(con, "group.label")
nSparse <- 0L
nTries <- 1L
while ( (nSparse <= maxSparse) & (nTries <= maxNonconv) ) {
## permute grouping variable
d[ , G] <- sample(d[ , G])
## transform data?
if (!is.null(datafun)) {
extraArgs <- formals(datafun)
neededArgs <- intersect(names(extraArgs), names(availableArgs))
extraArgs <- do.call(c, lapply(neededArgs, function(nn) availableArgs[nn]))
extraArgs$data <- d
originalNames <- colnames(d)
d <- do.call(datafun, extraArgs)
## coerce extraOut to data.frame
if (!is.data.frame(d)) stop('Argument "datafun" did not return a data.frame')
if (!all(originalNames %in% colnames(d)))
stop('The data.frame returned by argument "datafun" did not contain ',
'column names required by the model:\n',
paste(setdiff(originalNames, colnames(d)), collapse = ", "))
}
## for ordered indicators, check that groups have same observed categories
ordVars <- lavaan::lavNames(con, type = "ov.ord")
if (length(ordVars) > 0) {
try(onewayTables <- lavaan::lavTables(d, dimension = 1L,
categorical = ordVars, group = G),
silent = TRUE)
if (exists("onewayTables")) {
if (any(onewayTables$obs.prop == 1)) {
nSparse <- nSparse + 1L
next
}
} else {
## no "onewayTables" probably indicates empty categories in 1+ groups
nSparse <- nSparse + 1L
next
}
}
## fit null model, if it exists
if (!is.null(null)) {
out.null <- lavaan::update(null, data = d, group.label = group.label)
}
## fit constrained model, check for convergence
try(out0 <- lavaan::update(con, data = d, group.label = group.label))
if (!exists("out0")) {
nTries <- nTries + 1L
next
}
if (!lavaan::lavInspect(out0, "converged")) {
nTries <- nTries + 1L
next
}
## fit unconstrained model (unless NULL), check for convergence
if (!is.null(uncon)) {
try(out1 <- lavaan::update(uncon, data = d, group.label = group.label))
if (!exists("out1")) {
nTries <- nTries + 1L
next
}
if (!lavaan::lavInspect(out1, "converged")) {
nTries <- nTries + 1L
next
}
}
## If you get this far, everything converged, so break WHILE loop
break
}
## if WHILE loop ended before getting results, return NA
if ( (nSparse == maxSparse) | (nTries == maxNonconv) ) {
allAFIs <- c(AFIs, moreAFIs)
AFI <- rep(NA, sum(!is.na(allAFIs)))
names(AFI) <- allAFIs[!is.na(allAFIs)]
MI <- if (is.null(param)) NULL else NA
extra.obs <- NA
nTries <- nTries + 1L
} else {
availableArgs$con <- out0
if (exists("out1")) availableArgs$uncon <- out1
if (exists("out.null")) availableArgs$null <- out.null
AFI <- do.call(getAFIs, availableArgs)
## save max(MI) if !is.null(param)
if (is.null(param)) {
MI <- NULL
} else {
MI <- max(do.call(getMIs, c(availableArgs, modelType = "mgcfa"))$X2)
}
## anything extra?
if (!is.null(extra)) {
extraArgs <- formals(extra)
neededArgs <- intersect(names(extraArgs), names(availableArgs))
extraArgs <- do.call(c, lapply(neededArgs, function(nn) availableArgs[nn]))
extraOut <- do.call(extra, extraArgs)
## coerce extraOut to data.frame
if (!is.list(extraOut)) extraOut <- as.list(extraOut)
extra.obs <- data.frame(extraOut)
} else extra.obs <- data.frame(NULL)
}
options(warn = old_warn)
list(AFI = AFI, MI = MI, extra = extra.obs,
n.nonConverged = nTries - 1L, n.Sparse = nSparse)
}
permuteOnce.mimic <- function(i, d, G, con, uncon, null, param, freeParam,
covariates, AFIs, moreAFIs, maxSparse, maxNonconv,
iseed, warn, extra = NULL, datafun = NULL) {
old_warn <- options()$warn
options(warn = warn)
## save arguments from call
argNames <- names(formals(permuteOnce.mimic))
availableArgs <- lapply(argNames, function(x) eval(as.name(x)))
names(availableArgs) <- argNames
group.label <- lavaan::lavInspect(con, "group.label")
nTries <- 1L
while (nTries <= maxNonconv) {
## permute covariate(s) within each group
if (length(G)) {
for (gg in group.label) {
dG <- d[ d[[G]] == gg, ]
N <- nrow(dG)
newd <- dG[sample(1:N, N), covariates, drop = FALSE]
for (COV in covariates) d[d[[G]] == gg, COV] <- newd[ , COV]
}
} else {
N <- nrow(d)
newd <- d[sample(1:N, N), covariates, drop = FALSE]
for (COV in covariates) d[ , COV] <- newd[ , COV]
}
## transform data?
if (!is.null(datafun)) {
extraArgs <- formals(datafun)
neededArgs <- intersect(names(extraArgs), names(availableArgs))
extraArgs <- do.call(c, lapply(neededArgs, function(nn) availableArgs[nn]))
extraArgs$data <- d
originalNames <- colnames(d)
d <- do.call(datafun, extraArgs)
## coerce extraOut to data.frame
if (!is.data.frame(d)) stop('Argument "datafun" did not return a data.frame')
if (!all(originalNames %in% colnames(d)))
stop('The data.frame returned by argument "datafun" did not contain ',
'column names required by the model:\n',
paste(setdiff(originalNames, colnames(d)), collapse = ", "))
}
## fit null model, if it exists
if (!is.null(null)) {
out.null <- lavaan::update(null, data = d, group.label = group.label)
}
## fit constrained model
try(out0 <- lavaan::update(con, data = d, group.label = group.label))
## check for convergence
if (!exists("out0")) {
nTries <- nTries + 1L
next
}
if (!lavaan::lavInspect(out0, "converged")) {
nTries <- nTries + 1L
next
}
## If you get this far, everything converged, so break WHILE loop
break
}
## if WHILE loop ended before getting results, return NA
if (nTries == maxNonconv) {
allAFIs <- c(AFIs, moreAFIs)
AFI <- rep(NA, sum(!is.na(allAFIs)))
names(AFI) <- allAFIs[!is.na(allAFIs)]
MI <- if (is.null(param)) NULL else NA
extra.obs <- NA
nTries <- nTries + 1L
} else {
availableArgs$con <- out0
if (exists("out.null")) availableArgs$null <- out.null
AFI <- do.call(getAFIs, availableArgs)
if (is.null(param)) {
MI <- NULL
} else {
MI <- max(do.call(getMIs, c(availableArgs, modelType = "mimic"))$X2)
}
## anything extra?
if (!is.null(extra)) {
extraArgs <- formals(extra)
neededArgs <- intersect(names(extraArgs), names(availableArgs))
extraArgs <- do.call(c, lapply(neededArgs, function(nn) availableArgs[nn]))
extraOut <- do.call(extra, extraArgs)
## coerce extraOut to data.frame
if (!is.list(extraOut)) extraOut <- as.list(extraOut)
extra.obs <- data.frame(extraOut)
} else extra.obs <- data.frame(NULL)
}
options(warn = old_warn)
list(AFI = AFI, MI = MI, extra = extra.obs,
n.nonConverged = nTries - 1L, n.Sparse = integer(length = 0))
}
Try the semTools package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.