# NOT RUN {
mod.cat <- ' FU1 =~ u1 + u2 + u3 + u4
FU2 =~ u5 + u6 + u7 + u8 '
## the 2 factors are actually the same factor (FU) measured twice
longFacNames <- list(FU = c("FU1","FU2"))
## CONFIGURAL model: no constraints across groups or repeated measures
syntax.config <- measEq.syntax(configural.model = mod.cat,
# NOTE: data provides info about numbers of
# groups and thresholds
data = datCat,
ordered = paste0("u", 1:8),
parameterization = "theta",
ID.fac = "std.lv", ID.cat = "Wu.Estabrook.2016",
group = "g", longFacNames = longFacNames)
## print lavaan syntax to the Console
cat(as.character(syntax.config))
## print a summary of model features
summary(syntax.config)
## THRESHOLD invariance:
## only necessary to specify thresholds if you have no data
mod.th <- '
u1 | t1 + t2 + t3 + t4
u2 | t1 + t2 + t3 + t4
u3 | t1 + t2 + t3 + t4
u4 | t1 + t2 + t3 + t4
u5 | t1 + t2 + t3 + t4
u6 | t1 + t2 + t3 + t4
u7 | t1 + t2 + t3 + t4
u8 | t1 + t2 + t3 + t4
'
syntax.thresh <- measEq.syntax(configural.model = c(mod.cat, mod.th),
# NOTE: data not provided, so syntax must
# include thresholds, and number of
# groups == 2 is indicated by:
sample.nobs = c(1, 1),
parameterization = "theta",
ID.fac = "std.lv", ID.cat = "Wu.Estabrook.2016",
group = "g", group.equal = "thresholds",
longFacNames = longFacNames,
long.equal = "thresholds")
## notice that constraining 4 thresholds allows intercepts and residual
## variances to be freely estimated in all but the first group & occasion
cat(as.character(syntax.thresh))
## print a summary of model features
summary(syntax.thresh)
## Fit a model to the data either in a subsequent step (recommended):
mod.config <- as.character(syntax.config)
fit.config <- cfa(mod.config, data = datCat, group = "g",
ordered = paste0("u", 1:8), parameterization = "theta")
## or in a single step (not generally recommended):
fit.thresh <- measEq.syntax(configural.model = mod.cat, data = datCat,
ordered = paste0("u", 1:8),
parameterization = "theta",
ID.fac = "std.lv", ID.cat = "Wu.Estabrook.2016",
group = "g", group.equal = "thresholds",
longFacNames = longFacNames,
long.equal = "thresholds", return.fit = TRUE)
## compare their fit to test threshold invariance
anova(fit.config, fit.thresh)
## --------------------------------------------------------
## RECOMMENDED PRACTICE: fit one invariance model at a time
## --------------------------------------------------------
## - A downside of setting return.fit=TRUE is that if the model has trouble
## converging, you don't have the opportunity to investigate the syntax,
## or even to know whether an error resulted from the syntax-generator or
## from lavaan itself.
## - A downside of automatically fitting an entire set of invariance models
## (like the old measurementInvariance() function did) is that you might
## end up testing models that shouldn't even be fitted because less
## restrictive models already fail (e.g., don't test full scalar
## invariance if metric invariance fails! Establish partial metric
## invariance first, then test equivalent of intercepts ONLY among the
## indicators that have invariate loadings.)
## The recommended sequence is to (1) generate and save each syntax object,
## (2) print it to the screen to verify you are fitting the model you expect
## to (and potentially learn which identification constraints should be
## released when equality constraints are imposed), and (3) fit that model
## to the data, as you would if you had written the syntax yourself.
## Continuing from the examples above, after establishing invariance of
## thresholds, we proceed to test equivalence of loadings and intercepts
## (metric and scalar invariance, respectively)
## simultaneously across groups and repeated measures.
# }
# NOT RUN {
## metric invariance
syntax.metric <- measEq.syntax(configural.model = mod.cat, data = datCat,
ordered = paste0("u", 1:8),
parameterization = "theta",
ID.fac = "std.lv", ID.cat = "Wu.Estabrook.2016",
group = "g", longFacNames = longFacNames,
group.equal = c("thresholds","loadings"),
long.equal = c("thresholds","loadings"))
summary(syntax.metric) # summarize model features
mod.metric <- as.character(syntax.metric) # save as text
cat(mod.metric) # print/view lavaan syntax
## fit model to data
fit.metric <- cfa(mod.metric, data = datCat, group = "g",
ordered = paste0("u", 1:8), parameterization = "theta")
## test equivalence of loadings, given equivalence of thresholds
anova(fit.thresh, fit.metric)
## scalar invariance
syntax.scalar <- measEq.syntax(configural.model = mod.cat, data = datCat,
ordered = paste0("u", 1:8),
parameterization = "theta",
ID.fac = "std.lv", ID.cat = "Wu.Estabrook.2016",
group = "g", longFacNames = longFacNames,
group.equal = c("thresholds","loadings",
"intercepts"),
long.equal = c("thresholds","loadings",
"intercepts"))
summary(syntax.scalar) # summarize model features
mod.scalar <- as.character(syntax.scalar) # save as text
cat(mod.scalar) # print/view lavaan syntax
## fit model to data
fit.scalar <- cfa(mod.scalar, data = datCat, group = "g",
ordered = paste0("u", 1:8), parameterization = "theta")
## test equivalence of intercepts, given equal thresholds & loadings
anova(fit.metric, fit.scalar)
## For a single table with all results, you can pass the models to
## summarize to the compareFit() function
compareFit(fit.config, fit.thresh, fit.metric, fit.scalar)
## ------------------------------------------------------
## NOT RECOMMENDED: fit several invariance models at once
## ------------------------------------------------------
test.seq <- c("thresholds","loadings","intercepts","means","residuals")
meq.list <- list()
for (i in 0:length(test.seq)) {
if (i == 0L) {
meq.label <- "configural"
group.equal <- ""
long.equal <- ""
} else {
meq.label <- test.seq[i]
group.equal <- test.seq[1:i]
long.equal <- test.seq[1:i]
}
meq.list[[meq.label]] <- measEq.syntax(configural.model = mod.cat,
data = datCat,
ordered = paste0("u", 1:8),
parameterization = "theta",
ID.fac = "std.lv",
ID.cat = "Wu.Estabrook.2016",
group = "g",
group.equal = group.equal,
longFacNames = longFacNames,
long.equal = long.equal,
return.fit = TRUE)
}
compareFit(meq.list)
## -----------------
## Binary indicators
## -----------------
## borrow example data from Mplus user guide
myData <- read.table("http://www.statmodel.com/usersguide/chap5/ex5.16.dat")
names(myData) <- c("u1","u2","u3","u4","u5","u6","x1","x2","x3","g")
bin.mod <- '
FU1 =~ u1 + u2 + u3
FU2 =~ u4 + u5 + u6
'
## Must SIMULTANEOUSLY constrain thresholds, loadings, and intercepts
test.seq <- list(strong = c("thresholds","loadings","intercepts"),
means = "means",
strict = "residuals")
meq.list <- list()
for (i in 0:length(test.seq)) {
if (i == 0L) {
meq.label <- "configural"
group.equal <- ""
long.equal <- ""
} else {
meq.label <- names(test.seq)[i]
group.equal <- unlist(test.seq[1:i])
# long.equal <- unlist(test.seq[1:i])
}
meq.list[[meq.label]] <- measEq.syntax(configural.model = bin.mod,
data = myData,
ordered = paste0("u", 1:6),
parameterization = "theta",
ID.fac = "std.lv",
ID.cat = "Wu.Estabrook.2016",
group = "g",
group.equal = group.equal,
#longFacNames = longFacNames,
#long.equal = long.equal,
return.fit = TRUE)
}
compareFit(meq.list)
## ---------------------
## Multilevel Invariance
## ---------------------
## To test invariance across levels in a MLSEM, specify syntax as though
## you are fitting to 2 groups instead of 2 levels.
mlsem <- ' f1 =~ y1 + y2 + y3
f2 =~ y4 + y5 + y6 '
## metric invariance
syntax.metric <- measEq.syntax(configural.model = mlsem, meanstructure = TRUE,
ID.fac = "std.lv", sample.nobs = c(1, 1),
group = "cluster", group.equal = "loadings")
## by definition, Level-1 means must be zero, so fix them
syntax.metric <- update(syntax.metric,
change.syntax = paste0("y", 1:6, " ~ c(0, NA)*1"))
## save as a character string
mod.metric <- as.character(syntax.metric, groups.as.blocks = TRUE)
## convert from multigroup to multilevel
mod.metric <- gsub(pattern = "group:", replacement = "level:",
x = mod.metric, fixed = TRUE)
## fit model to data
fit.metric <- lavaan(mod.metric, data = Demo.twolevel, cluster = "cluster")
summary(fit.metric)
# }
Run the code above in your browser using DataLab