New semTools function for testing measurement equivalence / invariance

3425 views
Longitudinal
announcement
categorical
invariance
measurement
multigroup
reproducibility
semTools
syntax
Skip to first unread message

Terrence Jorgensen

unread,
Aug 26, 2018, 6:02:43 PM8/26/18
to lavaan
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

Giacomo Mason

unread,
Aug 28, 2018, 5:51:56 AM8/28/18
to lavaan
Thank you for the hard work Terrence, I'll try this out soon.

Giacomo

Giacomo Mason

unread,
Aug 28, 2018, 11:24:25 AM8/28/18
to lavaan
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

unread,
Aug 29, 2018, 8:35:36 AM8/29/18
to lavaan
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?


And is devtools installed?

library(devtools)

Carlos Lemos

unread,
Aug 30, 2018, 1:17:16 PM8/30/18
to lav...@googlegroups.com
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

unread,
Sep 10, 2018, 7:06:32 AM9/10/18
to lav...@googlegroups.com
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

unread,
Sep 12, 2018, 3:58:58 PM9/12/18
to lavaan

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

unread,
Sep 12, 2018, 4:29:30 PM9/12/18
to lavaan

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.

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

Terrence Jorgensen

unread,
Sep 12, 2018, 4:37:53 PM9/12/18
to lavaan

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

Carlos Lemos

unread,
Sep 13, 2018, 11:08:24 AM9/13/18
to lav...@googlegroups.com
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

unread,
Oct 12, 2018, 11:56:23 AM10/12/18
to lavaan


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

unread,
Oct 18, 2018, 7:04:39 AM10/18/18
to lavaan
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

unread,
Oct 19, 2018, 5:37:58 AM10/19/18
to lavaan

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

unread,
Oct 20, 2018, 5:17:41 PM10/20/18
to lavaan
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

unread,
Oct 25, 2018, 11:51:40 AM10/25/18
to lav...@googlegroups.com
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:


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

unread,
Oct 26, 2018, 6:20:15 AM10/26/18
to lavaan
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:


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)


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.


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:
 

Terrence D. Jorgensen
Assistant Professor, Methods and Statistics

Carlos Lemos

unread,
Oct 27, 2018, 5:25:19 PM10/27/18
to lav...@googlegroups.com
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

unread,
Oct 31, 2018, 10:42:59 AM10/31/18
to lavaan
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

unread,
Nov 2, 2018, 9:48:03 AM11/2/18
to lavaan
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

Johanna Schüller

unread,
Nov 4, 2018, 10:59:26 AM11/4/18
to lav...@googlegroups.com
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"+ 
        "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)

This is the warning I receive back: 

lavaan 0.6-3 did not run (perhaps do.fit = FALSE)?
** WARNING ** Estimates below are simply the starting values

  Optimization method                           NLMINB
  Number of free parameters                          0

  Number of observations                             0

  Estimator                                         ML


Thanks again for your help,
Johanna

Jessica Fritz

unread,
Nov 7, 2018, 9:22:43 AM11/7/18
to lavaan
Hi Terry,

Many thanks for this absolutely brilliant tool!

I am trying to apply the Liu, Y., et al. (2017; doi:10.1037/met0000075) measurement invariance procedure for longitudinal, categorical measures to my data.

I have used the measEq.syntax() function as depicted below.

I compared the resulting syntax to the syntax in the Liu et al.'s (2017) manuscript (i.e. Appendix 3A). As far as I can see the syntax performs exactly as expected, with one exception:
Liu et al. specify the following for the mean structure:

# Latent common factor means

# Fix T1 common factor mean to zero

C1FAMo ~ 0*1

# Freely estimate T2-T4 common factor means

C2FAMo ~ 1

C3FAMo ~ 1

C4FAMo ~ 1


However, when I try to replicate the syntax with measEq.syntax() (as shown below) the first mean is not fixed to zero but is freely estimated:

## LATENT MEANS/INTERCEPTS:

FS_early ~ alpha.1*1 + NA*1
FS_late ~ alpha.2*1 + NA*1

I of course know that I manually can adapt the syntax and can run my model manually in congruence with the Liu et al (2017) recommendations.

However I was wondering whether this discrepancy is (a) a mistake I made by misspecifying the arguments in the measEq.syntax() function, or whether there is (b) indeed a discrepancy between the original Liu specifications and the measEq.syntax(ID.cat = "millsap.tein.2004") specifications?

If (b) is the case, is the discrepancy on purpose? And if yes, could you perhas help me understand why?

I'll attach my complete syntax so the interested reader can see that all other aspects of my syntax, created with the measEq.syntax(ID.cat = "millsap.tein.2004") function, seem to be congruent with the original syntax Liu et al (2017) suggested.

Very many thanks for any help,

Jes

F_Model <- 'FS_early =~ RF1FQ1_occ1 + RF1FQ2_occ1 + RF1FQ3_occ1 + RF1FQ4_occ1 + RF1FQ8_occ1
                    FS_late  =~  RF1FQ1_occ2 + RF1FQ2_occ2 + RF1FQ3_occ2 + RF1FQ4_occ2 + RF1FQ8_occ2'
longFacNames <- list(FS = c("FS_early","FS_late")); longFacNames
longIndNames <- list(RF1FQ1 = c("RF1FQ1_occ1","RF1FQ1_occ2"),
                     RF1FQ2 = c("RF1FQ2_occ1","RF1FQ2_occ2"),
                     RF1FQ3 = c("RF1FQ3_occ1","RF1FQ3_occ2"),
                     RF1FQ4 = c("RF1FQ4_occ1","RF1FQ4_occ2"),
                     RF1FQ8 = c("RF1FQ8_occ1","RF1FQ8_occ2")); longIndNames
#Run configural model
syntax.fit.F_Model_T <- measEq.syntax(configural.model = F_Model,
                                           data = F_Data,
                                           estimator="WLSMV",
                                           ordered=c("RF1FQ1_occ1", "RF1FQ2_occ1", "RF1FQ3_occ1", "RF1FQ4_occ1", "RF1FQ8_occ1",
                                                     "RF1FQ1_occ2", "RF1FQ2_occ2", "RF1FQ3_occ2", "RF1FQ4_occ2", "RF1FQ8_occ2"),
                                           parameterization = "theta",
                                           ID.fac = "UL",
                                           ID.cat = "millsap.tein.2004",
                                           longFacNames = longFacNames,
                                           longIndNames = longIndNames,
                                           auto = "all")
cat(as.character(syntax.fit.F_Model_T))
summary(syntax.fit.F_Model_T)



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

Terrence Jorgensen

unread,
Nov 8, 2018, 5:06:12 AM11/8/18
to lavaan
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"+ 
        "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"'

Is there still a problem if you remove the unnecessary quotation marks around each variable name in the syntax?

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+
            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
'

FYI, lavaan accepts a character vector of any length, so the syntax does not need to be a single string.  (My function actually exploits that fact.)  So with this many variables, it might be simpler to write your syntax automatically like this:

retest.model.FPP_lavaan <- c(paste("FPP_T1 =~", var1), paste("FPP_T2 =~", var2))

And you can still add other parameters to that character vector if necessary, for example, a residual covariance:

retest.model.FPP_lavaan <- c(paste("FPP_T1 =~", var1), paste("FPP_T2 =~", var2),
                             "PPQ24_A_Macht ~~ PPQ36_A_Macht")

Hope this helps,

Terrence Jorgensen

unread,
Nov 8, 2018, 3:37:29 PM11/8/18
to lavaan
when I try to replicate the syntax with measEq.syntax() (as shown below) the first mean is not fixed to zero but is freely estimated:

## LATENT MEANS/INTERCEPTS:

FS_early ~ alpha.1*1 + NA*1
FS_late ~ alpha.2*1 + NA*1

Woops!  You're right, thanks for finding this mistake.  I spent so much time making sure all the nuances would work to implement Wu & Estabrook's recommendations (because I am not a fan of Millsap & Tein's) that I hadn't checked the ID.cat="millsap" situation carefully enough.  It is now fixed in the development version, which you can install using:

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

Please let me know if you find any other discrepancies.  Thanks,

Jessica Fritz

unread,
Nov 9, 2018, 5:33:54 AM11/9/18
to lavaan
Hi Terry,

Thank you so much for looking at this so quickly!

I've installed the new version and as far as I can see everything in the measEq.syntax(ID.cat = "millsap.tein.2004") seems to be corret!

Thousand thanks for your help.

All the best,

Jes

Pavneet Kaur

unread,
Dec 1, 2018, 5:47:34 PM12/1/18
to lavaan

Hello Dr. Terrence.

First of all, Thanks for this measEq.syntax() function as it is really helpful.

I have a doubt regarding the group.equal ( ), as stated underneath:

 

type5 <- 't5=~ g4q1.5+ g4q3.5+ g4q5.5 + g4q7.5 + g5q1.5+ g5q3.5+ g5q5.5 + g5q7.5

 

g5q1.5 ~~ g5q5.5'

 

After testing invariance for threshold and loadings, using group.partial ( ) for one indicator, the fit indices for latent-varaible invariance  is given as:

 

test_th.lo.var<-measEq.syntax(type5, data=data,  ID.cat="Wu.Estabrook.2016",  group="group", group.equal=c("thresholds","loadings","lv.variances"), group.partial = c("type5=~g5q3.5","g5q3.5|t1"), warn=T, return.fit=T,  ordered=c("g4q1.5", "g4q3.5", "g4q5.5", "g4q7.5", "g5q1.5", "g5q3.5", "g5q5.5", "g5q7.5"),  parameterization="theta", estimator="wlsmv",    meanstructure=TRUE)

chisq.scaled     df.scaled pvalue.scaled    cfi.scaled  rmsea.scaled    tli.scaled

#63.048        46.000         0.048         0.955         0.020         0.945

Now, I am interested  in testing the invariance of latent-variable means , but it is yielding the same fit-indices as above:

test_th.lo.var.means<-measEq.syntax(type5, data=data,  ID.cat="Wu.Estabrook.2016",  group="group", group.equal=c("thresholds","loadings","lv.variances","means"), warn=T, return.fit=T, ordered = c("g4q1.5", "g4q3.5", "g4q5.5", "g4q7.5", "g5q1.5", "g5q3.5", "g5q5.5", "g5q7.5"), parameterization = "theta", estimator="wlsmv", meanstructure=TRUE)

 

#chisq.scaled     df.scaled pvalue.scaled    cfi.scaled  rmsea.scaled    tli.scaled

#63.048        46.000         0.048         0.955         0.020         0.945

I am unsure why these two fit-indices are same, though the “test_th.lo.var.means” has means set to be equal on both groups? Thanks in advance,

Pavneet Kaur

 

Terrence Jorgensen

unread,
Dec 3, 2018, 5:38:13 AM12/3/18
to lavaan

group.equal=c("thresholds","loadings","lv.variances","means")

You cannot validly test equality of means unless intercepts are equivalent.

I am unsure why these two fit-indices are same, though the “test_th.lo.var.means” has means set to be equal on both groups? 

Means were probably already held equal because intercepts were not held equal.

FYI, if you only have binary indicators, you need to simultaneously constrain group.equal=c("thresholds","loadings","intercepts") to equality after configural invariance, because there is not enough information to differentiate between those 3 potential sources of non-invariance.

Terrence D. Jorgensen
Assistant Professor, Methods and Statistics
Research Institute for Child Development and Education, the University of Amsterdam

Pavneet Bharaj

unread,
Dec 3, 2018, 2:51:26 PM12/3/18
to lav...@googlegroups.com

Thanks Professor Terrence for your reply back. That was useful.

I have fixed the “intercepts” to be equal in group.equal argument, but even them I am having the same issue.

I am wondering as group.equal for “means” provided a different fit index for model “t2”, but not “t5”. Even if all the indices remain same, I would expect a change in degrees of freedom. Please let me know if I am doing something wrong in the given code:

 

type5 <- 't5=~ g4q1.5+ g4q3.5+ g4q5.5 + g4q7.5 + g5q1.5+ g5q3.5+ g5q5.5 + g5q7.5 

          g5q1.5 ~~ g5q5.5'

#=====  CHECKING EQUAL VARIANCES ========#

test_th.lo.var<-measEq.syntax(type5, data=data,  ID.cat="Wu.Estabrook.2016",

                              group="group", group.equal=c("loadings","thresholds", "intercepts", "lv.variances"),

                               warn=T, return.fit=T,

                              ordered=c("g4q1.5", "g4q3.5", "g4q5.5", "g4q7.5",

                                        "g5q1.5", "g5q3.5", "g5q5.5", "g5q7.5"),

                               parameterization="theta", estimator="wlsmv",

                              meanstructure=TRUE)

summary(test_th.lo.var, fit.measures=T, rsquare=T)

#chisq.scaled     df.scaled pvalue.scaled    cfi.scaled  rmsea.scaled    tli.scaled

#62.004        46.000         0.058         0.957         0.019         0.948

 

#=====  CHECKING EQUAL MEANS ========#

test_th.lo.var.means<-measEq.syntax(type5, data=data,  ID.cat="Wu.Estabrook.2016", group="group", group.equal=c("loadings","thresholds", "intercepts", "lv.variances","means"),   warn=T, return.fit=T,

                                               ordered=c("g4q1.5", "g4q3.5", "g4q5.5", "g4q7.5",

                                              "g5q1.5", "g5q3.5", "g5q5.5", "g5q7.5"),

                                    parameterization="theta", estimator="wlsmv",

                                    meanstructure=TRUE)

summary(test_th.lo.var.means, fit.measures=T, rsquare=T)

#chisq.scaled     df.scaled pvalue.scaled    cfi.scaled  rmsea.scaled    tli.scaled

# 62.004        46.000         0.058         0.957         0.019         0.948

#==============================#

type2 <- 'type2=~ g4q1.2+ g4q3.2+ g4q5.2+ g4q7.2+ g5q1.2+ g5q3.2+ g5q5.2+  g5q7.2 

g4q3.2 ~~ g5q5.2

g5q1.2 ~~ g5q5.2'

#=====  CHECKING EQUAL VARIANCES ========#

test_partial_th.lo2.var<-measEq.syntax(type2, data=data,  ID.cat="Wu.Estabrook.2016",

                                       group="group", group.equal=c("thresholds","loadings", "intercepts","lv.variances"),

                                       group.partial=c("type2=~g5q3.2","g5q3.2|t1"),warn=T, return.fit=T,

                                       ordered=c("g4q1.2", "g4q3.2", "g4q5.2", "g4q7.2",

                                                 "g5q1.2", "g5q3.2", "g5q5.2", "g5q7.2"),

                                        parameterization="theta", estimator="wlsmv",

                                       meanstructure=TRUE)

summary(test_partial_th.lo2.var, fit.measures=T, rsquare=T)

#chisq.scaled     df.scaled pvalue.scaled    cfi.scaled  rmsea.scaled    tli.scaled

#  48.741        43.000         0.253         0.954         0.012         0.940

 

#=====  CHECKING EQUAL MEANS ========#

test_partial_th.lo2_var.means<-measEq.syntax(type2, data=data,  ID.cat="Wu.Estabrook.2016",

                                             group="group", group.equal=c("thresholds","loadings","intercepts", "lv.variances",  "means"),  warn=T, return.fit=T,

                                                      ordered=c("g4q1.2", "g4q3.2", "g4q5.2", "g4q7.2",

                                                       "g5q1.2", "g5q3.2", "g5q5.2", "g5q7.2"),

                                          parameterization="theta", estimator="wlsmv",  meanstructure=TRUE)

summary(test_partial_th.lo2_var.means, fit.measures=T, rsquare=T)

#chisq.scaled     df.scaled pvalue.scaled    cfi.scaled  rmsea.scaled    tli.scaled

#  55.472        44.000         0.115         0.909         0.017         0.884


Terrence Jorgensen

unread,
Dec 4, 2018, 8:08:23 AM12/4/18
to lavaan

I have fixed the “intercepts” to be equal in group.equal argument, but even them I am having the same issue.

Do you have the issue after installing the latest development version of semTools?

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

If so, I suggest you leave return.fit=FALSE so you can inspect the model syntax to see what is happening.  For example, does adding "means" to group.equal= not lose a df because the means are already constrained in the model without "means" in group.equal=?  Or because adding "means" to group.equal= does not give them the same labels to add the equality constraint, or because the means are not free parameters?  This will help me figure out if I need to update the function.

I am wondering as group.equal for “means” provided a different fit index for model “t2”, but not “t5”. 

The change occurred for type2 because you only included group.partial=c("type2=~g5q3.2","g5q3.2|t1") in the model without means constrained, not when means were constrained. 

Pavneet Bharaj

unread,
Dec 4, 2018, 11:10:11 AM12/4/18
to lav...@googlegroups.com

I have fixed the “intercepts” to be equal in group.equal argument, but even them I am having the same issue.

Do you have the issue after installing the latest development version of semTools?

 

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

 

 

Yeah, the issue is still there, even though I have installed semTools again using above code.

 

If so, I suggest you leave return.fit=FALSE so you can inspect the model syntax to see what is happening.  For example, does adding "means" to group.equal= not lose a df because the means are already constrained in the model without "means" in group.equal=?  Or because adding "means" to group.equal= does not give them the same labels to add the equality constraint, or because the means are not free parameters?  This will help me figure out if I need to update the function.

 

 

So, I tried to use fit.measure=F and to clarify I am providing the syntax and results underneath:

 

Return.fit=FALSE

test_th.lo.var <- measEq.syntax(type5, data=data,  ID.cat="Wu.Estabrook.2016",

                   group="group", group.equal=c("thresholds","loadings", "lv.variances"),

  warn=T, return.fit=F,

                                      ordered=c("g4q1.5", "g4q3.5", "g4q5.5", "g4q7.5",

                                                "g5q1.5", "g5q3.5", "g5q5.5", "g5q7.5"),

                                      parameterization="theta", estimator="wlsmv",

                                      meanstructure=TRUE)

 

summary(test_th.lo.var)

 

 

output:

This lavaan model syntax specifies a CFA with 8 manifest indicators (8 of which are ordinal) of 1 common factor(s).

 

To identify the location and scale of each common factor, the factor means and variances were fixed to 0 and 1, respectively, unless equality constraints on measurement parameters allow them to be freed.

 

The location and scale of each latent item-response underlying 8 ordinal indicators were identified using the "theta" parameterization, and the identification constraints recommended by Wu & Estabrook (2016). For details, read:

 

        https://doi.org/10.1007/s11336-016-9506-0

 

Pattern matrix indicating num(eric), ord(ered), and lat(ent) indicators per factor:

 

       t5

g4q1.5 ord

g4q3.5 ord

g4q5.5 ord

g4q7.5 ord

g5q1.5 ord

g5q3.5 ord

g5q5.5 ord

g5q7.5 ord

 

The following types of parameter were constrained to equality across groups:

 

        thresholds

        loadings

        lv.variances

 

 

test_th.lo.var.means <- measEq.syntax(type5, data=data,  ID.cat="Wu.Estabrook.2016",

                                      group="group", group.equal=c("thresholds","loadings","lv.variances","means"),

                                      warn=T, return.fit=F,

                                      ordered=c("g4q1.5", "g4q3.5", "g4q5.5", "g4q7.5",

                                                "g5q1.5", "g5q3.5", "g5q5.5", "g5q7.5"),

                                     parameterization="theta", estimator="wlsmv",

                                      meanstructure=TRUE)

summary(test_th.lo.var.means)

 

 

OUTPUT:

This lavaan model syntax specifies a CFA with 8 manifest indicators (8 of which are ordinal) of 1 common factor(s).
 
To identify the location and scale of each common factor, the factor means and variances were fixed to 0 and 1, respectively, unless equality constraints on measurement parameters allow them to be freed.
 
The location and scale of each latent item-response underlying 8 ordinal indicators were identified using the "theta" parameterization, and the identification constraints recommended by Wu & Estabrook (2016). For details, read:
 
        https://doi.org/10.1007/s11336-016-9506-0 
 
Pattern matrix indicating num(eric), ord(ered), and lat(ent) indicators per factor:
 
       t5 
g4q1.5 ord
g4q3.5 ord
g4q5.5 ord
g4q7.5 ord
g5q1.5 ord
g5q3.5 ord
g5q5.5 ord
g5q7.5 ord
 
The following types of parameter were constrained to equality across groups:
 
        thresholds
        loadings
        lv.variances
        means

 

 

Syntax and output with fit.measure=T

 

test_th.lo.var<-measEq.syntax(type5, data=data,  ID.cat="Wu.Estabrook.2016",

                                      group="group", group.equal=c("thresholds","loadings", "lv.variances"),

                                     

                                      warn=T, return.fit=T,

                                      ordered=c("g4q1.5", "g4q3.5", "g4q5.5", "g4q7.5",

                                                "g5q1.5", "g5q3.5", "g5q5.5", "g5q7.5"),

                                     

                                      parameterization="theta", estimator="wlsmv",

                                      meanstructure=TRUE)

 

summary(test_th.lo.var, fit.measures=T, rsquare=T)

 

lavaan 0.6-3 ended normally after 57 iterations
 
  Optimization method                           NLMINB
  Number of free parameters                         42
  Number of equality constraints                    16
 
  Number of observations per group         
  0                                                919
  1                                                926
 
  Estimator                                       DWLS      Robust
  Model Fit Test Statistic                      59.510      63.048
  Degrees of freedom                                46          46
  P-value (Chi-square)                           0.087       0.048
  Scaling correction factor                                  0.943
  Shift parameter for each group:
    0                                                       -0.018
    1                                                       -0.018
    for simple second-order correction (Mplus variant)
 
Chi-square for each group:
 
  0                                             34.682      36.748
  1                                             24.827      26.301
 
Model test baseline model:
 
  Minimum Function Test Statistic              466.605     432.226
  Degrees of freedom                                56          56
  P-value                                        0.000       0.000
 
User model versus baseline model:
 
  Comparative Fit Index (CFI)                    0.967       0.955
  Tucker-Lewis Index (TLI)                       0.960       0.945
 
  Robust Comparative Fit Index (CFI)                            NA
  Robust Tucker-Lewis Index (TLI)                               NA
 
Root Mean Square Error of Approximation:
 
  RMSEA                                          0.018       0.020
  90 Percent Confidence Interval          0.000  0.030       0.002  0.031
  P-value RMSEA <= 0.05                          1.000       1.000
 
  Robust RMSEA                                                  NA
  90 Percent Confidence Interval                                NA     NA
 
Standardized Root Mean Square Residual:
 
  SRMR                                           0.073       0.073
 
Parameter Estimates:
 
  Information                                 Expected
  Information saturated (h1) model        Unstructured
  Standard Errors                           Robust.sem
 
 
Group 1 [0]:
 
Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)
  t5 =~                                               
    g4q1.5  (l.1_)    0.555    0.065    8.508    0.000
    g4q3.5  (l.2_)    1.013    0.143    7.076    0.000
    g4q5.5  (l.3_)    0.601    0.071    8.412    0.000
    g4q7.5  (l.4_)    0.273    0.110    2.471    0.013
    g5q1.5  (l.5_)    0.264    0.048    5.540    0.000
    g5q3.5  (l.6_)   -0.019    0.043   -0.452    0.652
    g5q5.5  (l.7_)    0.377    0.065    5.803    0.000
    g5q7.5  (l.8_)   -0.149    0.044   -3.389    0.001
 
Covariances:
                   Estimate  Std.Err  z-value  P(>|z|)
 .g5q1.5 ~~                                           
   .g5q5.5  (t.7_)    0.173    0.077    2.255    0.024
 
Intercepts:
                   Estimate  Std.Err  z-value  P(>|z|)
   .g4q1.5  (n.1.)    0.000                           
   .g4q3.5  (n.2.)    0.000                           
   .g4q5.5  (n.3.)    0.000                           
   .g4q7.5  (n.4.)    0.000                           
   .g5q1.5  (n.5.)    0.000                           
   .g5q3.5  (n.6.)    0.000                           
   .g5q5.5  (n.7.)    0.000                           
   .g5q7.5  (n.8.)    0.000                           
    t5      (a.1.)    0.000                           
 
Thresholds:
                   Estimate  Std.Err  z-value  P(>|z|)
    g41.5|1 (g41.)    0.638    0.053   12.151    0.000
    g43.5|1 (g43.)    0.581    0.073    7.932    0.000
    g45.5|1 (g45.)   -0.364    0.050   -7.239    0.000
    g47.5|1 (g47.)    2.188    0.118   18.611    0.000
    g51.5|1 (g51.)   -0.132    0.043   -3.064    0.002
    g53.5|1 (g53.)   -0.094    0.041   -2.275    0.023
    g55.5|1 (g55.)    1.349    0.066   20.578    0.000
    g57.5|1 (g57.)   -0.087    0.042   -2.077    0.038
 
Variances:
                   Estimate  Std.Err  z-value  P(>|z|)
   .g4q1.5  (t.1_)    1.000                           
   .g4q3.5  (t.2_)    1.000                           
   .g4q5.5  (t.3_)    1.000                           
   .g4q7.5  (t.4_)    1.000                           
   .g5q1.5  (t.5_)    1.000                           
   .g5q3.5  (t.6_)    1.000                           
   .g5q5.5  (t.7_)    1.000                           
   .g5q7.5  (t.8_)    1.000                           
    t5      (p.1_)    1.000                           
 
Scales y*:
                   Estimate  Std.Err  z-value  P(>|z|)
    g4q1.5            0.875                           
    g4q3.5            0.703                           
    g4q5.5            0.857                           
    g4q7.5            0.965                           
    g5q1.5            0.967                           
    g5q3.5            1.000                           
    g5q5.5            0.936                           
    g5q7.5            0.989                           
 
R-Square:
                   Estimate
    g4q1.5            0.235
    g4q3.5            0.506
    g4q5.5            0.265
    g4q7.5            0.069
    g5q1.5            0.065
    g5q3.5            0.000
    g5q5.5            0.125
    g5q7.5            0.022
 
 
Group 2 [1]:
 
Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)
  t5 =~                                               
    g4q1.5  (l.1_)    0.555    0.065    8.508    0.000
    g4q3.5  (l.2_)    1.013    0.143    7.076    0.000
    g4q5.5  (l.3_)    0.601    0.071    8.412    0.000
    g4q7.5  (l.4_)    0.273    0.110    2.471    0.013
    g5q1.5  (l.5_)    0.264    0.048    5.540    0.000
    g5q3.5  (l.6_)   -0.019    0.043   -0.452    0.652
    g5q5.5  (l.7_)    0.377    0.065    5.803    0.000
    g5q7.5  (l.8_)   -0.149    0.044   -3.389    0.001
 
Covariances:
                   Estimate  Std.Err  z-value  P(>|z|)
 .g5q1.5 ~~                                           
   .g5q5.5  (t.7_)    0.154    0.067    2.289    0.022
 
Intercepts:
                   Estimate  Std.Err  z-value  P(>|z|)
   .g4q1.5  (n.1.)   -0.044    0.071   -0.619    0.536
   .g4q3.5  (n.2.)    0.085    0.086    0.997    0.319
   .g4q5.5  (n.3.)    0.046    0.070    0.655    0.513
   .g4q7.5  (n.4.)    0.022    0.145    0.154    0.878
   .g5q1.5  (n.5.)   -0.025    0.061   -0.412    0.680
   .g5q3.5  (n.6.)    0.025    0.059    0.430    0.667
   .g5q5.5  (n.7.)    0.328    0.080    4.114    0.000
   .g5q7.5  (n.8.)   -0.040    0.059   -0.684    0.494
    t5      (a.1.)    0.000                           
 
Thresholds:
                   Estimate  Std.Err  z-value  P(>|z|)
    g41.5|1 (g41.)    0.638    0.053   12.151    0.000
    g43.5|1 (g43.)    0.581    0.073    7.932    0.000
    g45.5|1 (g45.)   -0.364    0.050   -7.239    0.000
    g47.5|1 (g47.)    2.188    0.118   18.611    0.000
    g51.5|1 (g51.)   -0.132    0.043   -3.064    0.002
    g53.5|1 (g53.)   -0.094    0.041   -2.275    0.023
    g55.5|1 (g55.)    1.349    0.066   20.578    0.000
    g57.5|1 (g57.)   -0.087    0.042   -2.077    0.038
 
Variances:
                   Estimate  Std.Err  z-value  P(>|z|)
   .g4q1.5  (t.1_)    1.000                           
   .g4q3.5  (t.2_)    1.000                           
   .g4q5.5  (t.3_)    1.000                           
   .g4q7.5  (t.4_)    1.000                           
   .g5q1.5  (t.5_)    1.000                           
   .g5q3.5  (t.6_)    1.000                           
   .g5q5.5  (t.7_)    1.000                           
   .g5q7.5  (t.8_)    1.000                           
    t5      (p.1_)    1.000                           
 
Scales y*:
                   Estimate  Std.Err  z-value  P(>|z|)
    g4q1.5            0.875                           
    g4q3.5            0.703                           
    g4q5.5            0.857                           
    g4q7.5            0.965                           
    g5q1.5            0.967                           
    g5q3.5            1.000                           
    g5q5.5            0.936                           
    g5q7.5            0.989                           
 
R-Square:
                   Estimate
    g4q1.5            0.235
    g4q3.5            0.506
    g4q5.5            0.265
    g4q7.5            0.069
    g5q1.5            0.065
    g5q3.5            0.000
    g5q5.5            0.125
    g5q7.5            0.022

 

 

 

test_th.lo.var.means<-measEq.syntax(type5, data=data,  ID.cat="Wu.Estabrook.2016",

                                      group="group", group.equal=c("thresholds","loadings","lv.variances","means"),

                                      warn=T, return.fit=T,

                                   ordered=c("g4q1.5", "g4q3.5", "g4q5.5", "g4q7.5",

                                                "g5q1.5", "g5q3.5", "g5q5.5", "g5q7.5"),

                                  parameterization="theta", estimator="wlsmv",

                                      meanstructure=TRUE)

summary(test_th.lo.var.means, fit.measures=T, rsquare=T)

 

lavaan 0.6-3 ended normally after 57 iterations
 
  Optimization method                           NLMINB
  Number of free parameters                         42
  Number of equality constraints                    16
 
  Number of observations per group         
  0                                                919
  1                                                926
 
  Estimator                                       DWLS      Robust
  Model Fit Test Statistic                      59.510      63.048
  Degrees of freedom                                46          46
  P-value (Chi-square)                           0.087       0.048
  Scaling correction factor                                  0.943
  Shift parameter for each group:
    0                                                       -0.018
    1                                                       -0.018
    for simple second-order correction (Mplus variant)
 
Chi-square for each group:
 
  0                                             34.682      36.748
  1                                             24.827      26.301
 
Model test baseline model:
 
  Minimum Function Test Statistic