# NOT RUN {
# }
# 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
'
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)
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab