New semTools function for testing measurement equivalence / invariance

[Terrence Jorgensen]

For the past few years, there have been 3 separate functions for automatically testing the omnibus null hypotheses of measurement invariance, each under very constrained conditions (e.g., simple structure, first-order measurement). 

• measurementInvariance(): only for numeric indicators, only tests across multiple groups
• measurementInvarianceCat(): only for ordered indicators, only tests across multiple groups
• longInvariance(): only for numeric indicators of 1 (repeatedly measured) factor, only tests across repeated measures

In the next version of semTools, these functions will still available, but they have been deprecated.  

The new function measEq.syntax() addresses all of the above limitations by writing lavaan model syntax and (optionally) fitting the model to data.  By specifying a configural model with minimal syntax (i.e., the pattern of factor loadings), users can set arguments to specify models with invariance constraints imposed on any measurement and/or structural parameters across any combination of multiple groups and/or repeated measures.  The configural model may have complex structure (i.e., cross-loadings and residual covariances are allowed), higher-order factors, and any combination of numeric and ordered indicators.  Different repeatedly measured factors do not need to all be measured on the same number of occasions, nor do the same indicators even need to be measured on all occasions (e.g., if some test items are replaced when measuring abilities at later ages).  Beyond the choice between std.lv = TRUE or FALSE in the previous functions, users can also choose effects-coding to identify common-factor distributions (except in models with cross-loadings, such as bifactor models).  Users can also choose between 4 different methods for identifying the latent item-response distributions underlying any ordered indicators, no longer being confined to use the Millsap & Tein (2004) recommendations implemented in measurementInvarianceCat().  The measEq.syntax() function (and now the compareFit() function) can also be applied to lavaan.mi objects, if users fit lavaan models to multiple imputations.  If a user still has a more complex situation than can be automated using the measEq.syntax() function, it is still possible to specify a model with similar needs, then print the lavaan model syntax and edit it manually to fit the desired model.

I encourage anyone who has had difficulty using the previous invariance-testing functions to try using this, and please let me know if there are any difficulties.  Users can install the development version of semTools and find examples on the help page to get started:

devtools::install_github("simsem/semTools/semTools")
library(semTools)
?measEq.syntax

The help page examples show how the user can mimic the Console output of the previous functions by looping over a sequence of constraints, saving the fitted models in a list, and comparing their fit using the compareFit() function.  Eventually, I will write an interactive wrapper function to do this, integrating follow-up tests rather than only providing omnibus tests that assume all previous levels of invariance are tenable.  But that will only be possible after updating the existing partialInvariance() function to work in this new framework, and I probably will not have time to work on that soon. 

Thanks for any feedback (please post here so others can benefit),

Terrence D. Jorgensen

Postdoctoral Researcher, Methods and Statistics

Research Institute for Child Development and Education, the University of Amsterdam

UvA web page: http://www.uva.nl/profile/t.d.jorgensen

[Giacomo Mason]

Thank you for the hard work Terrence, I'll try this out soon.

Giacomo

[Giacomo Mason]

I can't seem to be able to install semTools from GitHub using devtools. Is that possible?

Best,

G

On Sunday, August 26, 2018 at 11:02:43 PM UTC+1, Terrence Jorgensen wrote:

[Terrence Jorgensen]

I can't seem to be able to install semTools from GitHub using devtools

My colleague was able to install it.  Do you have the latest R (3.5.1) and development version of lavaan (0.6-3.1295 installed?

http://lavaan.ugent.be/development.html

And is devtools installed?

library(devtools)

[Carlos Lemos]

I can't thank you enough for this, Terrence!

I had a hard time implementing the Wu & Estabrook procedure following your advice, and I can imagine the amount of expertise and work required to write this function.

I will certainly try it and give feedback.

Best regards,

Carlos

--
You received this message because you are subscribed to the Google Groups "lavaan" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lavaan+un...@googlegroups.com.
To post to this group, send email to lav...@googlegroups.com.
Visit this group at https://groups.google.com/group/lavaan.
For more options, visit https://groups.google.com/d/optout.

[Carlos Lemos]

Hi again Terrence,

I’m trying to use the new measEq.Syntax() function to study measurement invariance on the ISSP religion data. Here is the structure of the data frame:

> str(df)

'data.frame':   35198 obs. of  20 variables:

 $ V28         : Ord.factor w/ 3 levels "I don't believe in God"<..: 3 3 2 2 2 2 2 2 2 2 ...

 $ V30         : Ord.factor w/ 4 levels "No, definitely not"<..: 2 4 2 2 4 NA 1 NA 2 1 ...

 $ V31         : Ord.factor w/ 4 levels "No, definitely not"<..: 2 4 2 1 NA 1 1 NA 1 3 ...

 $ V32         : Ord.factor w/ 4 levels "No, definitely not"<..: 2 2 2 1 NA 1 1 1 2 3 ...

 $ V33         : Ord.factor w/ 4 levels "No, definitely not"<..: 3 4 2 1 1 3 1 2 2 3 ...

 $ V35         : Ord.factor w/ 5 levels "Strongly disagree"<..: 2 5 3 1 2 NA 1 3 3 NA ...

 $ V37         : Ord.factor w/ 5 levels "Strongly disagree"<..: 3 2 3 1 1 1 1 3 3 NA ...

 $ V46         : Ord.factor w/ 4 levels "Never"<"Yearly"<..: 4 3 3 1 2 2 3 2 1 3 ...

 $ V47         : Ord.factor w/ 4 levels "Never"<"Yearly"<..: 3 NA 2 1 2 NA 1 1 1 3 ...

 $ V48         : Ord.factor w/ 4 levels "Never"<"Yearly"<..: 4 4 2 1 2 4 3 1 2 3 ...

 $ V49         : Ord.factor w/ 5 levels "Never"<"Yearly"<..: 4 1 2 1 1 1 1 1 1 3 ...

 $ V50         : Ord.factor w/ 5 levels "Never"<"Yearly"<..: 2 1 1 2 1 1 1 1 1 2 ...

 $ V51         : Ord.factor w/ 5 levels "Highly non-religious"<..: 4 1 2 1 1 4 2 2 2 2 ...

 $ ATTEND      : Ord.factor w/ 4 levels "Never"<"Yearly"<..: 4 1 2 1 1 2 1 1 1 3 ...

 $ AGE         : Ord.factor w/ 5 levels "0-24"<"25-44"<..: 4 2 2 2 2 3 3 3 3 3 ...

 $ SEX         : Factor w/ 2 levels "Male","Female": 2 2 1 1 1 1 2 1 1 2 ...

 $ DEGREE      : Ord.factor w/ 6 levels "No formal qualification"<..: 2 4 4 6 4 6 3 3 4 3 ...

 $ RELIGGRP    : Factor w/ 5 levels "No religion",..: 2 2 4 1 1 1 2 1 1 2 ...

 $ COUNTRY.NAME: Factor w/ 26 levels "AU-Australia",..: 2 3 4 5 7 8 9 12 13 17 ...

 $ YEAR        : Ord.factor w/ 1 level "2008": 1 1 1 1 1 1 1 1 1 1 ...

>

> 

 

I defined the following items as ordinal and converted to numeric, and for simplicity I used “SEX” as the group factor (just two levels):

> indicators <- c("V28", "V30", "V31", "V32", "V33",        

+                 "V35", "V37",         

+                 "V46", "V47", "V48", "V49", "V50",        

+                 "V51", "ATTEND")

> for (ind in indicators) {

+   df[,ind] <- as.numeric(df[,ind])

+ }

> group <- "SEX"

> 

I now tried to compute a configural model, first with “cfa” and then with “measEq.syntax”, for a single factor:

> model <- "afterlife =~ V30 + V31 + V32 + V33"

> # CONFIGURAL

> config <- cfa(model, data = df, ordered = indicators, parameterization = "theta",

+               estimator = "WLSMV", group = group, group.equal = "")

> config <- measEq.syntax(model, data = df, ordered = indicators, parameterization = "theta",

+                         ID.fac = "std.lv", ID.cat = "Wu.Estabrook.2016", estimator = "WLSMV",

+                         group = group, group.equal = "", return.fit = TRUE)

Error in specify[[1]][RR, CC] : subscript out of bounds

> 

The first runs, the second gives an error. However, if I define this 2-factor model:

model <- "afterlife =~ V30 + V31 + V32 + V33

          belief.God =~ V28 + V35 + V37"

 

it runs OK:

> config <- cfa(model, data = df, ordered = indicators, parameterization = "theta",

+               estimator = "WLSMV", group = group, group.equal = "")

> config <- measEq.syntax(model, data = df, ordered = indicators, parameterization = "theta",

+                         ID.fac = "std.lv", ID.cat = "Wu.Estabrook.2016", estimator = "WLSMV",

+                         group = group, group.equal = "", return.fit = TRUE)

> 

What am I doing wrong?

I proceeded to compute a model with constrained thresholds and loadings, which tests for metric invariance (scales of both latent response variables and factors fixed across groups):

> metric <- measEq.syntax(model, data = df, ordered = indicators, parameterization = "theta",

+                         ID.fac = "std.lv", ID.cat = "Wu.Estabrook.2016", estimator = "WLSMV",

+                         group = group, group.equal = c("thresholds","loadings"),

                          return.fit = TRUE)

> 

and it runs OK and has good fit (cfi = 0.99936, rmsea = 0.03707745, srmr = 0.0183). I inspected the parameter table and checked that the identification conditions in eq. (19) of Wu and Estabrook are nicely in place (thank you very much for this!!!). However, when I check the estimator:

> unlist(lavInspect(metric, what = "Options"))["estimator"]

estimator

   "DWLS"

> 

I don’t get “WLSMV”, which is recommended for MGCFA with ordinal variables.

Finally, I have another doubt. I’m trying to run multiple imputed files with measEq.syntax(), but since this takes a long time with “mice”, I thought of saving the imputed data frames in a “mids” object and then running the analysis. However, I could not figure out how to do this with the help page of measEq.syntax().

Thank you very much once again for writing this wonderful function and for your help and attention.

Sincerely,

 

Carlos M. Lemos

Em seg, 27 de ago de 2018 às 00:02, Terrence Jorgensen <tjorge...@gmail.com> escreveu:

[dvillarre...@gmail.com]

Hi Terrence,

Thank you very much for this function.

I have a problem running my model, it is a longitudinal invariance analysis with categorical variables.

They are a total of seven measurements (DEP0 to DEP6), with nine items in each measurement. The items score from 0 to 3, however, in some follow-ups the items only reach scores from 0 to 2.

 

 

model <- "

DEP0 =~ f0V1 + f0V2 + f0V3 + f0V4 + f0V5 + f0V6 + f0V7 + f0V8 + f0V9
DEP1 =~ f1V1 + f1V2 + f1V3 + f1V4 + f1V5 + f1V6 + f1V7 + f1V8 + f1V9
DEP2 =~ f2V1 + f2V2 + f2V3 + f2V4 + f2V5 + f2V6 + f2V7 + f2V8 + f2V9
DEP3 =~ f3V1 + f3V2 + f3V3 + f3V4 + f3V5 + f3V6 + f3V7 + f3V8 + f3V9
DEP4 =~ f4V1 + f4V2 + f4V3 + f4V4 + f4V5 + f4V6 + f4V7 + f4V8 + f4V9
DEP5 =~ f5V1 + f5V2 + f5V3 + f5V4 + f5V5 + f5V6 + f5V7 + f5V8 + f5V9
DEP6 =~ f6V1 + f6V2 + f6V3 + f6V4 + f6V5 + f6V6 + f6V7 + f6V8 + f6V9
"


var0 <- c("f0V1", "f0V2", "f0V3", "f0V4", "f0V5", "f0V6", "f0V7", "f0V8", "f0V9")
var1 <- c("f1V1", "f1V2", "f1V3", "f1V4", "f1V5", "f1V6", "f1V7", "f1V8", "f1V9")
var2 <- c("f2V1", "f2V2", "f2V3", "f2V4", "f2V5", "f2V6", "f2V7", "f2V8", "f2V9")    
var3 <- c("f3V1", "f3V2", "f3V3", "f3V4", "f3V5", "f3V6", "f3V7", "f3V8", "f3V9")
var4 <- c("f4V1", "f4V2", "f4V3", "f4V4", "f4V5", "f4V6", "f4V7", "f4V8", "f4V9")
var5 <- c("f5V1", "f5V2", "f5V3", "f5V4", "f5V5", "f5V6", "f5V7", "f5V8", "f5V9")
var6 <- c("f6V1", "f6V2", "f6V3", "f6V4", "f6V5", "f6V6", "f6V7", "f6V8", "f6V9")
 

longFacNames <- list(FU = c("DEP0","DEP1","DEP2","DEP3","DEP4","DEP5","DEP6"))



syntax.config <- measEq.syntax(configural.model = model, data = DEP_SAL, parameterization = "theta",
                        ID.fac = "std.lv", ID.cat = "millsap", estimator = "WLSMV",
                        return.fit = TRUE, longFacNames=longFacNames,
                        ordered = c("f0V1", "f0V2", "f0V3", "f0V4", "f0V5", "f0V6", "f0V7", "f0V8", "f0V9",
                                    "f1V1", "f1V2", "f1V3", "f1V4", "f1V5", "f1V6", "f1V7", "f1V8", "f1V9",
                                    "f2V1", "f2V2", "f2V3", "f2V4", "f2V5", "f2V6", "f2V7", "f2V8", "f2V9",
                                    "f3V1", "f3V2", "f3V3", "f3V4", "f3V5", "f3V6", "f3V7", "f3V8", "f3V9",
                                    "f4V1", "f4V2", "f4V3", "f4V4", "f4V5", "f4V6", "f4V7", "f4V8", "f4V9",
                                    "f5V1", "f5V2", "f5V3", "f5V4", "f5V5", "f5V6", "f5V7", "f5V8", "f5V9",
                                    "f6V1", "f6V2", "f6V3", "f6V4","f6V5", "f6V6", "f6V7", "f6V8", "f6V9"))

 

cat(as.character(syntax.config))


summary(syntax.config, fit.measures=TRUE, standardized=TRUE)  

When I mail the analysis I get this error

Error in lav_samplestats_step1(Y = Data, ov.names = ov.names, ov.types = ov.types,  :
  lavaan ERROR: some categories of variable `f0V8' are empty in group 1; frequencies are [1670 45 43 0]

[Terrence Jorgensen]

I defined the following items as ordinal and converted to numeric

That is not necessary. If the ordered variables are already declared as ordered in the data.frame, then you don't even need to use the ordered= argument.

http://lavaan.ugent.be/tutorial/cat.html

Error in specify[[1]][RR, CC] : subscript out of bounds

I'm glad you tried it with both a 1-factor and 2-factor model. That helped me locate the source of the error and add a check so that syntax is only written for factor covariances if there are any.  If you install the development version of semTools again, it should work

when I check the estimator:

> unlist(lavInspect(metric, what = "Options"))["estimator"]

estimator

   "DWLS"

> 

I don’t get “WLSMV”, which is recommended for MGCFA with ordinal variables.

WLSMV is not an estimator.  That is just the shortcut for setting estimator = "DWLS", test = "scaled.shifted", and se = "robust.sem".  Nothing to do with this function; check the result for config when fitted with cfa().

Finally, I have another doubt. I’m trying to run multiple imputed files with measEq.syntax(), but since this takes a long time with “mice”, I thought of saving the imputed data frames in a “mids” object and then running the analysis. However, I could not figure out how to do this with the help page of measEq.syntax().

Yes, you should impute your data first so you are not wasting time repeating the imputation process before each analysis.  If you fit the configural model to the imputed data with runMI(), you can submit that object to measEq.syntax().  Here is an example.

set.seed(12345)
HSMiss <- HolzingerSwineford1939[ , paste("x", 1:9, sep = "")]
HSMiss$x5 <- ifelse(HSMiss$x1 <= quantile(HSMiss$x1, .3), NA, HSMiss$x5)
HSMiss$x9 <- ifelse(is.na(HSMiss$x5), NA, HSMiss$x9)
HSMiss$school <- HolzingerSwineford1939$school
HS.amelia <- amelia(HSMiss, m = 20, noms = "school")
imps <- HS.amelia$imputations

HS.model <- '
visual  =~ x1 + a*x2 + b*x3
textual =~ x4 + x5 + x6
speed   =~ x7 + x8 + x9
ab := a*b
'
mgfit1 <- cfa.mi(HS.model, data = imps, std.lv = TRUE, group = "school")

syntax.mi <- measEq.syntax(configural.model = mgfit1, group = "school",
                           group.equal = c("loadings","intercepts"))
cat(as.character(syntax.mi))
summary(syntax.mi)

## return fit
fit.mi <- measEq.syntax(configural.model = mgfit1, group = "school",
                        group.equal = c("loadings","intercepts"),
                        return.fit = TRUE)
summary(fit.mi)

Terrence D. Jorgensen

Postdoctoral Researcher, Methods and Statistics

Research Institute for Child Development and Education, the University of Amsterdam

http://www.uva.nl/profile/t.d.jorgensen

[Terrence Jorgensen]

When I mail the analysis I get this error

Error in lav_samplestats_step1(Y = Data, ov.names = ov.names, ov.types = ov.types,  :
  lavaan ERROR: some categories of variable `f0V8' are empty in group 1; frequencies are [1670 45 43 0]

The source of the error is clearly noted in the error message, and it has nothing to do with the syntax-writing function.  If you fit the configural model yourself with cfa(), you will notice the same error occurs.  

Terrence D. Jorgensen

Postdoctoral Researcher, Methods and Statistics

Research Institute for Child Development and Education, the University of Amsterdam

http://www.uva.nl/profile/t.d.jorgensen

[Carlos Lemos]

Hi again Terrence,

I'm glad that my feedback was somehow useful to you. I will reinstall semTools and retry.

Thank you very much for also clarifying my confusion about WLSMV and for sending the example on combining cfa.mi() with measEq.syntax(). I will make sure to set std.lv = TRUE in the call to cfa.mi (I missed that in my example) and then ID.fac = "std.lv" and ID.cat = "Wu.Estabrook.2016" in the call to measEq.syntax().

Best regards,

Carlos

[Una]

Thank you for all the work, and making our work simpler. I am now trying the function to test longitudinal variance, and I'm having problems with the long.partial argument.

This is the syntax:

mi.i<-measEq.syntax(ep.m, data= data,estimator="MLM", meanstructure=TRUE, int.lv.free=TRUE, auto.fix.first=T , orthogonal=T,
                    ID.fac="auto.fix.first",  return.fit = TRUE,
                    longFacNames =  list(X = c("X1","X2")), auto=T,
                    long.equal=c("loadings","intercepts"), long.partial = c("x2.1 ~ 1","x2.2 ~ 1"))

And this is the warning message:

Warning message: In if (long.partial == "") { : the condition has length > 1 and only the first element will be used

If I use only long.partial="x2 ~ 1", it doesn't give a warning message, but it still equals the intercepts.

Hope this feedback helps...

Kind regards,

Una

[Terrence Jorgensen]

long.partial = c("x2.1 ~ 1","x2.2 ~ 1"))

And this is the warning message:

Warning message: In if (long.partial == "") { : the condition has length > 1 and only the first element will be used

If I use only long.partial="x2 ~ 1", it doesn't give a warning message, but it still equals the intercepts.

 

I fixed the warning by checking the length too, but it is not something that would have actually led to a problem for anyone.

The problem (I am guessing) is that you are using actual observed-variable names in the vector: Are x2.1 and x2.2 the same variable measured at times 1 and 2, respectively?  

When you use long.partial="x2 ~ 1", that seems to refer to a name that you give to the observed variables x2.1 and x2.2 that are repeatedly measured (i.e., x2 itself is not a variable in your data set, right?).  In order to do that, you need to tell the function that "x2" refers to c("x2.1","x2.2") using the longIndNames= argument, which works the same way longFacNames= does.  As described in the "Repeated Measures" part of the Details section in the help page, if you don't make your own longIndNames, the default will be to assign automated labels based on the order, using the factor name.  So, since your longitudinal factor is named "X", the automated indicator name would probably be "._X_ind.2" (assuming it is the second indicator specified after =~ both times).  But of course I would recommend assigning your own longIndNames so you can use sensible names in long.partial=.

[Una]

Thank you for the feedback, I did judge wrongly how the automated labels would look like. If I understood correctly, in long.partial I can not use actual observed-variables names (which were in my case  x2.1  and  x2.2, the same variable measured two times, as you supposed) and I have to use a name refered to in longIndNames. But if I only define one indicator name in longIndNames, it seems to me it does not test the invariance of other indicators (those not defined in longIndNames). Some of my models are quite small, so this is not a problem to name them one by one. But for the bigger CFA models I was quite thankful for the automated indicator label formation, but it seems it is not compatible with partial invariance. I tried a few options of forming the label to use in long.partial based on ""._factor_ind.1", but did not find the right one (the one that would free the parameter). Is there maybe a way to access the automated indicator names, and then use those names when testing partial invariance?

Thank you for all the effort put in SEM tools,

Una

[Terrence Jorgensen]

If I understood correctly, in long.partial I can not use actual observed-variables names

Correct

I have to use a name refered to in longIndNames.

Correct

Is there maybe a way to access the automated indicator names, and then use those names when testing partial invariance?

Yes, and I should probably add this to the help page.  I store it in the @call slot, even if it was not part of your actual call.  From the help-page example:

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, data = datCat,
                               ordered = paste0("u", 1:8),
                               parameterization = "theta",

                               ID.fac = "std.lv", ID.cat = "Wu.Estabrook.2016",

                               group = "g", longFacNames = longFacNames)

## extract default longIndNames
syntax.config@call$longIndNames

Even when you fit the model, I store the measEq.syntax object in lavaan's @external slot.  So from your posted example, here is how you can find the names:

mi.i@external$measEq.syntax@call$longIndNames

Thank you for all the effort put in SEM tools,

Happy to help :-)

[Carlos Lemos]

Hi Terrence,

I got a warning with group.partial that looks similar to that reported by Una. I'm testing a 4-factor model with 10 ordinal indicators, and I'm trying to establish partial scalar invariance across countries (group) by freeing one indicator in each factor (Wu & Estabrook parameterization, no marker indicators) and I get the following:

> # SCALAR.PARTIAL
> scalar.partial <- measEq.syntax(model, data = df, ordered = indicators, parameterization = "theta",

+                         ID.fac = "std.lv", ID.cat = "Wu.Estabrook.2016", estimator = "WLSMV",

+                         group = group, group.equal = c("thresholds","loadings","intercepts"),
+                         group.partial = c("V33~1","V37~1","V48~1","ATTEND~1"), return.fit = TRUE)
Warning message:
In if (group.partial == "") { :

  the condition has length > 1 and only the first element will be used

It seems that the function is working OK because using lavInspect(scalar.partial, what = "free") I get (only two countries shown):

...
$`NO-Norway`$nu
       intrcp
V30         0
V31         0
V32         0
V33      1579
V28         0
V35         0
V37      1580
V46         0
V47         0
V48      1581
V49         0
V51         0
ATTEND   1582

...

$`IT-Italy`$nu
       intrcp
V30         0
V31         0
V32         0
V33      1751
V28         0
V35         0
V37      1752
V46         0
V47         0
V48      1753
V49         0
V51         0
ATTEND   1754

i.e. the intercepts I want are free and different for each country. I don't know whether I'm doing some error and/or the warning is serious.

I also have some questions not directly related to measEq.syntax(). If you can answer them without spending much time I would be very grateful:

- Is there any function in lavaan to test the equality of the group variance-covariance matrices (general and/or multiple pairwise)?

- How can I get Chi-square contributions for each group other than using summary()?

- Some authors state that 0.05 <= RMSEA < 0.08 means acceptable fit, but others (particularly Hu & Bentler (1999)) recommend 0.06 as cutoff for acceptable fit. I'm using the more restrictive limit, particularly because a recent paper warned on the use of criteria based on ML estimation for assessing fit of models with ordinal items obtained using WLSMV:

https://www.tandfonline.com/doi/abs/10.1080/10705511.2014.882658

Do you know about any recent study that recommends criteria for assessing fit in measurement invariance tests of complex models with large samples and ordinal items?

Thank you once more for your attention, and for creating measEq.syntax()and making it available!

Carlos M. Lemos

[Terrence Jorgensen]

It seems that the function is working OK

Correct.  As I told Una above, that warning is not something that will cause a problem for anyone, so you can ignore it.  But in the development version, the warning is gone.

devtools::install_github("simsem/semTools/semTools")

Is there any function in lavaan to test the equality of the group variance-covariance matrices (general and/or multiple pairwise)?

No, I've always just written my own saturated-model syntax using outer() with paste().

varnames <- paste0("x", 1:6)

## specify model parameters
covstruc <- outer(varnames, varnames, function(x, y) paste(x, "~~", y))
satMod <- c(paste(varnames, "~ 1"), # mean structure
            covstruc[lower.tri(covstruc, diag = TRUE)]) # covariance structure
cat(satMod, sep = "\n")

## fit model, unconstrained across groups
fit1 <- lavaan(satMod, data = HolzingerSwineford1939, group = "school")

## fit model, all parameters constrained across groups
fit0 <- lavaan(satMod, data = HolzingerSwineford1939, group = "school",
               group.equal = c("intercepts","residuals","residual.covariances"))

Note that "intercepts" constrains all observed-variable intercepts, even when the intercepts happen to be means because there are no predictors.  Likewise, "residuals" and "residual.covariances" constrain all observed-variable total (co)variances when there happen to be no predictors. 

Because the saturated model has df = 0, the test of equality is merely the chi-squared from fit0.

How can I get Chi-square contributions for each group other than using summary()?

Careful, it isn't really valid to compare those across multigroup models once any parameters are constrained across groups.  But they are stored in the object's @test slot:

fit@test[[1]]$stat.group # naïve
fit@test[[2]]$stat.group # robust

- Some authors state that 0.05 <= RMSEA < 0.08 means acceptable fit, but others (particularly Hu & Bentler (1999)) recommend 0.06 as cutoff for acceptable fit. I'm using the more restrictive limit, particularly because a recent paper warned on the use of criteria based on ML estimation for assessing fit of models with ordinal items obtained using WLSMV:

https://www.tandfonline.com/doi/abs/10.1080/10705511.2014.882658

So you are using a fixed cutoff even though the article said it's not a good idea?

Do you know about any recent study that recommends criteria for assessing fit in measurement invariance tests of complex models with large samples and ordinal items?

Yes (warning: shameless self-promotion)

https://www.tandfonline.com/doi/abs/10.1080/10705511.2017.1421467

You can use the semTools function permuteMeasEq() to do this.  This paper is specifically about getting better control of Type I errors with categorical indicators, but see our first paper showing why fixed cutoffs give inconsistent results even under ideal circumstances.

http://dx.doi.org/10.1037/met0000152

It also shows that permutation provides an actual test of configural invariance, even if your model does not fit perfectly.  If your configural model does not fit well, consider this paper too:

 

https://www.frontiersin.org/articles/10.3389/fpsyg.2017.01455/full 

Terrence D. Jorgensen

Assistant Professor, Methods and Statistics

[Carlos Lemos]

Hi Terrence,

I'm once again extremely grateful for your help and guidance!

Congratulations also for your papers - they are great, and likely to be useful not only to me but to many people in this forum. I hope I get our paper published, to be able to cite yours'!

The problem of the cutoff criteria is very important, because in the model for the countries choosing the cutoff 0.06 instead of 0.08 for RMSEA makes an enormous difference - the RMSEA for the scalar-invariant model is 0.077 (0.075,0.078) [90%].

In fact, two of my questions are related. In a test comparing a metric-invariant model for religious groups (equal thresholds and loadings across groups) with a model with invariant thresholds, loadings, and factor variance-covariance matrices (= lv.variances and = lv.covariances), I get RMSEA = 0.072, so I'm rejecting the hypothesis of invariant \Phi across groups. I thought that one indirect way of confirming the rejection criteria would be to test the \Phi_{g} matrices for equality in the scalar invariant model (metric+intercepts fixed, but \Phi again free across groups), which holds OK (good fit indices, small enough delta(CFI), etc.).

So, I will have to study your papers and your code, which will take me some time...

Best regards, and many, many thanks,

Carlos M. Lemos

[Johanna Schüller]

Dear Terrence,

I tried to use the measEq.syntax() function but I receive the error message:  

"lavaan 0.6-3 did not run (perhaps do.fit = FALSE)?

** WARNING ** Estimates below are simply the starting values"

Would you help me out and tell me where I went wrong? 

 Here is my code:

retest.model.FPP <- 'FPP_T1 =~ "PPQ7_Emp"+ "PPQ19_Emp"+ "PPQ25_Emp"+ "PPQ31_Emp"+ "PPQ37_Emp"+ 

        "PPQ2_Furchtl"+ "PPQ8_Furchtl"+ "PPQ14_Furchtl"+  "PPQ26_Furchtl"+ "PPQ38_Furchtl"+ 

        "PPQ3_Ego"+ "PPQ9_Ego"+ "PPQ21_Ego"+ "PPQ27_Ego"+ "PPQ39_Ego"+ 

        "PPQ4_Impul"+ "PPQ10_Impul"+ "PPQ16_Impul"+ "PPQ22_Impul"+ "PPQ40_Impul"+ 

        "PPQ11_Soz"+ "PPQ17_Soz"+ "PPQ23_Soz"+ "PPQ29_Soz"+ "PPQ35_Soz" + 

        "PPQ6_Macht" + "PPQ12_Macht"+ "PPQ18_Macht"+ "PPQ24_Macht"+ "PPQ36_Macht"

        FPP_T2 =~ "PPQ7_A_Emp" + "PPQ19_A_Emp" + "PPQ25_A_Emp"+ "PPQ31_A_Emp"+ "PPQ37_A_Emp" + 

        "PPQ2_A_Furchtl"+ "PPQ8_A_Furchtl"+ "PPQ14_A_Furchtl"+  "PPQ26_A_Furchtl"+ "PPQ38_A_Furchtl"+ 

        "PPQ3_A_Ego"+ "PPQ9_A_Ego"+ "PPQ21_A_Ego"+ "PPQ27_A_Ego"+ "PPQ39_A_Ego"+ 

        "PPQ4_A_Impul"+ "PPQ10_A_Impul"+ "PPQ16_A_Impul"+ "PPQ22_A_Impul"+ "PPQ40_A_Impul"+ 

        "PPQ11_A_Soz"+ "PPQ17_A_Soz" + "PPQ23_A_Soz"+ "PPQ29_A_Soz"+ "PPQ35_A_Soz"+

        "PPQ6_A_Macht"+ "PPQ12_A_Macht" + "PPQ18_A_Macht"+ "PPQ24_A_Macht"+ "PPQ36_A_Macht"'

var1 <- c("PPQ7_Emp", "PPQ19_Emp", "PPQ25_Emp", "PPQ31_Emp", "PPQ37_Emp", 

            "PPQ2_Furchtl", "PPQ8_Furchtl", "PPQ14_Furchtl",  "PPQ26_Furchtl", "PPQ38_Furchtl", 

            "PPQ3_Ego", "PPQ9_Ego", "PPQ21_Ego", "PPQ27_Ego", "PPQ39_Ego", 

            "PPQ4_Impul", "PPQ10_Impul", "PPQ16_Impul", "PPQ22_Impul", "PPQ40_Impul", 

            "PPQ11_Soz", "PPQ17_Soz", "PPQ23_Soz", "PPQ29_Soz", "PPQ35_Soz" , 

            "PPQ6_Macht" , "PPQ12_Macht", "PPQ18_Macht", "PPQ24_Macht", "PPQ36_Macht")

var2 <- c("PPQ7_A_Emp" , "PPQ19_A_Emp" , "PPQ25_A_Emp", "PPQ31_A_Emp", "PPQ37_A_Emp" , 

            "PPQ2_A_Furchtl", "PPQ8_A_Furchtl", "PPQ14_A_Furchtl",  "PPQ26_A_Furchtl", "PPQ38_A_Furchtl", 

            "PPQ3_A_Ego", "PPQ9_A_Ego", "PPQ21_A_Ego", "PPQ27_A_Ego", "PPQ39_A_Ego", 

            "PPQ4_A_Impul", "PPQ10_A_Impul", "PPQ16_A_Impul", "PPQ22_A_Impul", "PPQ40_A_Impul", 

            "PPQ11_A_Soz", "PPQ17_A_Soz" , "PPQ23_A_Soz", "PPQ29_A_Soz", "PPQ35_A_Soz",

            "PPQ6_A_Macht", "PPQ12_A_Macht" , "PPQ18_A_Macht", "PPQ24_A_Macht", "PPQ36_A_Macht")

 

longFacNames <- list(FU = c("FPP_T1", "FPP_T2"))

measEq.syntax(retest.model.FPP_lavaan, #spezifiziertes Modell über Lavaan

              ID.fac = "std.lv", #standardisierter latenter Faktor

              ID.cat = "Wu", #default Guideline

              longFacNames = longFacNames,

              return.fit = TRUE)

              

Thank you very much in advance!

Johanna 

[Terrence Jorgensen]

I receive the error message:  

"lavaan 0.6-3 did not run (perhaps do.fit = FALSE)?

** WARNING ** Estimates below are simply the starting values"

That is a warning, not an error.  

Would you help me out and tell me where I went wrong? 

 Here is my code:

retest.model.FPP <- 'FPP_T1 =~ "PPQ7_Emp"+ "PPQ19_Emp"+ "PPQ25_Emp"+ "PPQ31_Emp"+ "PPQ37_Emp"+ 

        "PPQ2_Furchtl"+ "PPQ8_Furchtl"+ "PPQ14_Furchtl"+  "PPQ26_Furchtl"+ "PPQ38_Furchtl"+ 

        "PPQ3_Ego"+ "PPQ9_Ego"+ "PPQ21_Ego"+ "PPQ27_Ego"+ "PPQ39_Ego"+ 

        "PPQ4_Impul"+ "PPQ10_Impul"+ "PPQ16_Impul"+ "PPQ22_Impul"+ "PPQ40_Impul"+ 

        "PPQ11_Soz"+ "PPQ17_Soz"+ "PPQ23_Soz"+ "PPQ29_Soz"+ "PPQ35_Soz" + 

        "PPQ6_Macht" + "PPQ12_Macht"+ "PPQ18_Macht"+ "PPQ24_Macht"+ "PPQ36_Macht"

        FPP_T2 =~ "PPQ7_A_Emp" + "PPQ19_A_Emp" + "PPQ25_A_Emp"+ "PPQ31_A_Emp"+ "PPQ37_A_Emp" + 

        "PPQ2_A_Furchtl"+ "PPQ8_A_Furchtl"+ "PPQ14_A_Furchtl"+  "PPQ26_A_Furchtl"+ "PPQ38_A_Furchtl"+ 

        "PPQ3_A_Ego"+ "PPQ9_A_Ego"+ "PPQ21_A_Ego"+ "PPQ27_A_Ego"+ "PPQ39_A_Ego"+ 

        "PPQ4_A_Impul"+ "PPQ10_A_Impul"+ "PPQ16_A_Impul"+ "PPQ22_A_Impul"+ "PPQ40_A_Impul"+ 

        "PPQ11_A_Soz"+ "PPQ17_A_Soz" + "PPQ23_A_Soz"+ "PPQ29_A_Soz"+ "PPQ35_A_Soz"+

        "PPQ6_A_Macht"+ "PPQ12_A_Macht" + "PPQ18_A_Macht"+ "PPQ24_A_Macht"+ "PPQ36_A_Macht"'

measEq.syntax(retest.model.FPP_lavaan, #spezifiziertes Modell über Lavaan

              ID.fac = "std.lv", #standardisierter latenter Faktor

              ID.cat = "Wu", #default Guideline

              longFacNames = longFacNames,

              return.fit = TRUE) 

It can't be all of your code, because you are not submitting the syntax object to measEq.syntax().  Did you fit the model first?  Is that when the warning message was printed?

Terrence D. Jorgensen

Assistant Professor, Methods and Statistics

Research Institute for Child Development and Education, the University of Amsterdam

http://www.uva.nl/profile/t.d.jorgensen

[Johanna Schüller]

Thanks for the answer. This is my correct code:

library(semTools)

retest.model.FPP_lavaan <- 'FPP_T1 =~ "PPQ7_Emp"+ "PPQ19_Emp"+ "PPQ25_Emp"+ "PPQ31_Emp"+ "PPQ37_Emp"+ 

        "PPQ2_Furchtl"+ "PPQ8_Furchtl"+ "PPQ14_Furchtl"+  "PPQ26_Furchtl"+ "PPQ38_Furchtl"+ 

        "PPQ3_Ego"+ "PPQ9_Ego"+ "PPQ21_Ego"+ "PPQ27_Ego"+ "PPQ39_Ego"+ 

        "PPQ4_Impul"+ "PPQ10_Impul"+ "PPQ16_Impul"+ "PPQ22_Impul"+ "PPQ40_Impul"+ 

        "PPQ11_Soz"+ "PPQ17_Soz"+ "PPQ23_Soz"+ "PPQ29_Soz"+ "PPQ35_Soz" + 

        "PPQ6_Macht" + "PPQ12_Macht"+ "PPQ18_Macht"+ "PPQ24_Macht"+ "PPQ36_Macht"

        FPP_T2 =~ "PPQ7_A_Emp" + "PPQ19_A_Emp" + "PPQ25_A_Emp"+ "PPQ31_A_Emp"+ "PPQ37_A_Emp" + 

        "PPQ2_A_Furchtl"+ "PPQ8_A_Furchtl"+ "PPQ14_A_Furchtl"+  "PPQ26_A_Furchtl"+ "PPQ38_A_Furchtl"+ 

        "PPQ3_A_Ego"+ "PPQ9_A_Ego"+ "PPQ21_A_Ego"+ "PPQ27_A_Ego"+ "PPQ39_A_Ego"+