## Not run:
#
# ########################
# ## 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
# '
# miout <- measurementInvariance(mod.config, data = HS, std.lv = TRUE,
# group = "ageGroup")
#
# (fit.config <- miout[["fit.configural"]])
# (fit.metric <- miout[["fit.loadings"]])
# (fit.scalar <- miout[["fit.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 the interaction terms after permuting the covariate
# datafun <- function(data) {
# d <- data[, !names(data) %in% 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)
#
# ## End(Not run)
Run the code above in your browser using DataLab