Introduction

This document is an Appendix of the main presentation Measurement Invariance with Structural Equation Modeling. In particular, this Appendix may be useful for the interpretation of the section MI with Ordinal Variables

The items refer to the Italian version of the Short Negative Acts Questionnaire (SNAQ; Balducci et al., 2010; Notelaers et al., 2019)

It is important to specify that (due to the imbalance of the cells) two categories (4 and 5) were collapsed to one category (i.e., 4). Below, you can find the syntax with the relative item content (note that the value of 5 is collapsed in the value of 4)1.

# Proseguiamo l’intervista chiedendole di indicare con che frequenza, negli ultimi sei mesi, ha subìto ciascuno dei seguenti comportamenti nel suo luogo di lavoro
# 1 = mai, 2 = Una volta in tutto o di tanto in tanto, 3 = mensilmente, 4 = Settimanalmente, 5 = Quotidianamente 

dat <- mobbing_data2 %>%  mutate(  
  snaq1=replace(snaq1, snaq1==5, 4),
  # Le sono state nascoste informazioni che influenzano il suo lavoro
  snaq2=replace(snaq2, snaq2==5, 4),
  # Sono state diffuse chiacchiere e dicerie nei suoi riguardi
  snaq3=replace(snaq3, snaq3==5, 4),
  # Lei è stato/a ignorato/a, escluso/a o emarginato/a
  snaq4=replace(snaq4, snaq4==5, 4),
  # Sono state fatte osservazioni offensive sulla sua persona (ad es. sulle sue abitudini), sui suoi atteggiamenti o sulla sua vita privata
  snaq5=replace(snaq5, snaq5==5, 4),
  # Hanno alzato la voce con lei o è stato/a bersaglio di attacchi istintivi di rabbia
  snaq6=replace(snaq6, snaq6==5, 4),
  # Le hanno fatto notare i suoi errori
  snaq7=replace(snaq7, snaq7==5, 4),
  # Nel rivolgersi agli altri è stato/a ignorato/a o si è imbattuto in reazioni ostili
  snaq8=replace(snaq8, snaq8==5, 4),
  # Hanno criticato il suo lavoro ed il suo impegno
    snaq9=replace(snaq9, snaq9==5, 4)) %>%
  # Le sono stati fatti scherzi spiacevoli da persone con le quali non va d’accordo    
  select(gender, snaq1:snaq9, mob1_numeric) # the variable mob1_numeric is not necessary for the scope of this Appendix

Furthermore, given that this document is an extract of a wider article (Perinelli, Balducci, & Fraccaroli, under review, see also here http://dx.doi.org/10.13140/RG.2.2.26701.10720), keep in mind that several packages listed below are not necessary for the invariance analysis reported subsequently.

Note that I have used the procedure and the scripts described in Svetina, Rutkowski, & Rutkowski (2020, Multiple-group invariance with categorical outcomes using updated guidelines: An illustration using Mplus and the lavaan/semTools packages. Structural Equation Modeling, 27, 111–130. https://doi.org/10.1080/10705511.2019.1602776)

Load Libraries and Environment

library(dplyr)
library(readxl)
library(psych)
library(apaTables)
library(formattable)
library(lavaan)
library(semTools)
library(MplusAutomation)
library(pROC)
library(randomForest)
library(WVPlots)
library(gridExtra)

# added for this appendix
library(summarytools)

load("../.RData")

Gender \(x\) Age Stats

The sample consists of 357 Italian workers

summarytools::st_css(bootstrap = FALSE) # set to FALSE, otherwise too wide tables
mobbing_data2 %$% summarytools::ctable(age_range,
                                       gender, style = 'rmarkdown',
                                       chisq = TRUE,
                                       headings = FALSE) %>%
  print(method="render")
gender
age_range male female Total
18-34 34 ( 51.5% ) 32 ( 48.5% ) 66 ( 100.0% )
35-54 90 ( 45.0% ) 110 ( 55.0% ) 200 ( 100.0% )
55-oltre 53 ( 58.2% ) 38 ( 41.8% ) 91 ( 100.0% )
Total 177 ( 49.6% ) 180 ( 50.4% ) 357 ( 100.0% )
 Χ2 = 4.5082   df = 2   p = .1050

Generated by summarytools 1.0.1 (R version 4.2.0)
2022-06-16

Kendall rank correlations (in APA style)

The package apaTables only formats tables for Pearson correlations.

Luckily, here you can find how to use that package also for Kendall rank correlations: https://github.com/dstanley4/apaTables/issues/18
Note that you have to install a apa_tables_mod.zip file in a proper folder, and then call it with the source function.

source("../analysis.R/apa_tables_mod/apaCorrelationTable_mod.R")
source("../analysis.R/apa_tables_mod//apaTables.R")
source("../analysis.R/apa_tables_mod/rtfMakeTable.R")
source("../analysis.R/apa_tables_mod/rtfMakeDocument.R")

library(tidyverse)

Whole sample

mobbing_data2 %>%
  select(snaq1, snaq2, snaq3, snaq4, snaq5, snaq6, snaq7, snaq8,snaq9) %>%
  apa.cor.table(., table.number=1, landscape = F, cor.method = "kendall") %>% 
  pluck(., 3) %>%
  knitr::kable(.)
Variable M SD 1 2 3 4 5 6 7 8
1. snaq1 1.92 1.13
2. snaq2 1.64 0.93 .43**
[.34, .51]
3. snaq3 1.69 1.04 .49** .44**
[.41, .57] [.35, .52]
4. snaq4 1.38 0.77 .38** .56** .43**
[.28, .46] [.49, .63] [.34, .51]
5. snaq5 1.58 0.82 .33** .29** .29** .46**
[.23, .42] [.20, .39] [.20, .39] [.37, .54]
6. snaq6 2.10 0.89 .25** .32** .19** .32** .29**
[.15, .35] [.23, .41] [.08, .28] [.22, .41] [.19, .38]
7. snaq7 1.72 0.91 .50** .45** .48** .44** .40** .31**
[.41, .57] [.36, .53] [.40, .56] [.35, .52] [.31, .48] [.21, .40]
8. snaq8 1.57 0.82 .39** .46** .31** .44** .41** .38** .45**
[.30, .47] [.37, .54] [.22, .40] [.36, .52] [.32, .49] [.28, .46] [.36, .53]
9. snaq9 1.14 0.52 .30** .35** .32** .41** .25** .22** .26** .32**
[.20, .39] [.26, .44] [.23, .41] [.32, .50] [.15, .35] [.12, .32] [.16, .36] [.23, .41]

Male sample

mobbing_data2 %>%
  filter (gender == "male") %>%
  select(snaq1, snaq2, snaq3, snaq4, snaq5, snaq6, snaq7, snaq8,snaq9) %>%
  apa.cor.table(., table.number=1, landscape = F, cor.method = "kendall") %>% 
  pluck(., 3) %>%
  knitr::kable(.)
Variable M SD 1 2 3 4 5 6 7 8
1. snaq1 1.94 1.17
2. snaq2 1.72 1.05 .48**
[.36, .59]
3. snaq3 1.67 1.12 .56** .48**
[.45, .66] [.36, .59]
4. snaq4 1.45 0.88 .43** .64** .43**
[.30, .54] [.55, .72] [.30, .54]
5. snaq5 1.58 0.86 .28** .37** .22** .44**
[.13, .41] [.24, .49] [.07, .35] [.31, .55]
6. snaq6 2.20 1.06 .27** .33** .24** .32** .31**
[.12, .40] [.19, .45] [.09, .37] [.18, .45] [.17, .44]
7. snaq7 1.73 0.96 .51** .49** .44** .46** .33** .27**
[.39, .61] [.37, .59] [.31, .55] [.33, .57] [.19, .45] [.13, .40]
8. snaq8 1.59 0.91 .38** .48** .28** .55** .42** .36** .41**
[.25, .50] [.36, .59] [.14, .41] [.44, .65] [.29, .54] [.22, .48] [.28, .53]
9. snaq9 1.16 0.59 .30** .36** .38** .41** .26** .24** .30** .40**
[.16, .42] [.22, .48] [.25, .50] [.28, .52] [.12, .39] [.10, .37] [.16, .43] [.27, .52]

Female sample

mobbing_data2 %>%
  filter (gender == "female") %>%
  select(snaq1, snaq2, snaq3, snaq4, snaq5, snaq6, snaq7, snaq8,snaq9) %>%
  apa.cor.table(., table.number=1, landscape = F, cor.method = "kendall") %>% 
  pluck(., 3) %>%
  knitr::kable(.)
Variable M SD 1 2 3 4 5 6 7 8
1. snaq1 1.90 1.09
2. snaq2 1.56 0.79 .38**
[.25, .50]
3. snaq3 1.71 0.97 .43** .40**
[.30, .54] [.26, .51]
4. snaq4 1.30 0.63 .32** .46** .45**
[.18, .44] [.34, .57] [.33, .56]
5. snaq5 1.58 0.79 .38** .22** .37** .49**
[.25, .50] [.07, .35] [.24, .49] [.37, .59]
6. snaq6 1.99 0.67 .23** .31** .13 .31** .28**
[.08, .36] [.17, .44] [-.01, .27] [.17, .43] [.14, .41]
7. snaq7 1.71 0.86 .49** .41** .53** .42** .48** .36**
[.37, .59] [.28, .53] [.42, .63] [.30, .54] [.35, .58] [.23, .48]
8. snaq8 1.54 0.73 .40** .44** .34** .33** .40** .41** .49**
[.27, .51] [.31, .55] [.21, .47] [.19, .45] [.27, .51] [.28, .52] [.37, .60]
9. snaq9 1.11 0.43 .30** .34** .26** .42** .24** .21** .23** .24**
[.16, .43] [.21, .47] [.12, .39] [.29, .54] [.10, .38] [.06, .34] [.08, .36] [.09, .37]

Reliability for ordinal data

Computation of the Nonlinear SEM Reliability Coefficient (Yang & Green, 2015) and AVE for S-NAQ, through reliability() function of semTools package (omega3 and avevar output, respectively)

fit_snaq_total <- cfa(model_snaq, data=dat, ordered=c("snaq1","snaq2","snaq3","snaq4", "snaq5",
                                                      "snaq6","snaq7", "snaq8", "snaq9"),
                      estimator="WLSMV")

fit_snaq_gender <- cfa(model_snaq, data=dat, ordered=c("snaq1","snaq2","snaq3","snaq4", "snaq5",
                                                       "snaq6","snaq7", "snaq8", "snaq9"),
                       estimator="WLSMV", group = "gender")

semTools::reliability(fit_snaq_total)
## For constructs with categorical indicators, Zumbo et al.`s (2007) "ordinal alpha" is calculated in addition to the standard alpha, which treats ordinal variables as numeric. See Chalmers (2018) for a critique of "alpha.ord" and the response by Zumbo & Kroc (2019). Likewise, average variance extracted is calculated from polychoric (polyserial) not Pearson correlations.
##             mobbing
## alpha     0.8786073
## alpha.ord 0.9231591
## omega     0.8827226
## omega2    0.8827226
## omega3    0.8960171
## avevar    0.5905185
semTools::reliability(fit_snaq_gender)
## For constructs with categorical indicators, Zumbo et al.`s (2007) "ordinal alpha" is calculated in addition to the standard alpha, which treats ordinal variables as numeric. See Chalmers (2018) for a critique of "alpha.ord" and the response by Zumbo & Kroc (2019). Likewise, average variance extracted is calculated from polychoric (polyserial) not Pearson correlations.
## $male
##             mobbing
## alpha     0.8848273
## alpha.ord 0.9268769
## omega     0.8909220
## omega2    0.8909220
## omega3    0.9117575
## avevar    0.6133607
## 
## $female
##             mobbing
## alpha     0.8706102
## alpha.ord 0.9200622
## omega     0.8795105
## omega2    0.8795105
## omega3    0.8977970
## avevar    0.5824758

Summary of reliability results

knitr::kable(data.frame
             (Sample = c("Whole",
                         "Male",
                         "Female"),
              rho_nl = c(round(semTools::reliability(fit_snaq_total, what = "omega3")[1,1], 2),
                         round(semTools::reliability(fit_snaq_gender, what = "omega3")$male[1,1],2),
                         round(semTools::reliability(fit_snaq_gender, what = "omega3")$female[1,1],2)),
              AVE = c(round(semTools::reliability(fit_snaq_total, what = "ave")[1,1], 2), 
                      round(semTools::reliability(fit_snaq_gender, what = "ave")$male[1,1],2),
                      round(semTools::reliability(fit_snaq_gender, what = "ave")$female[1,1],2))
               ),
             col.names = c("Sample", "$\\rho_{NL}$", "$AVE$"),
             format = "simple"
             )
Sample \(\rho_{NL}\) \(AVE\)
Whole 0.90 0.59
Male 0.91 0.61
Female 0.90 0.58

Scripts and Results

In the following tabsets you can find the three models for the invariance analysis with ordinal data (see Svetina, Rutkowski, & Rutkowski, 2020, for more info and Mplus equivalent scripts)

Baseline Model

#### INVARIANCE

# empty matrix in which we will put results
all.results<-matrix(NA, nrow = 3, ncol = 6)

# Specifying the baseline model with four items
mod.cat <- 'F1 =~ snaq1 + snaq2 + snaq3 + snaq4 + snaq5 + snaq6 + snaq7 + snaq8 + snaq9'

# Baseline model: no constraints across groups or repeated measures
baseline <- measEq.syntax(configural.model = mod.cat,
                          data = dat,
                          ordered = c("snaq1", "snaq2", "snaq3", "snaq4", "snaq5", "snaq6",
                                      "snaq7", "snaq8", "snaq9"),
                          parameterization = "delta",
                          ID.fac = "std.lv",
                          ID.cat = "Wu.Estabrook.2016",
                          group = "gender",
                          group.equal = "configural" )

# For a little bit of orientation/instructions in what model looks like
summary(baseline)
## This lavaan model syntax specifies a CFA with 9 manifest indicators (9 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 9 ordinal indicators were identified using the "delta" 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:
## 
##       F1 
## snaq1 ord
## snaq2 ord
## snaq3 ord
## snaq4 ord
## snaq5 ord
## snaq6 ord
## snaq7 ord
## snaq8 ord
## snaq9 ord
## 
## The following types of parameter were constrained to equality across groups:
## 
##  configural
# To see all of the constraints in the model
cat(as.character(baseline))
## ## LOADINGS:
## 
## F1 =~ c(NA, NA)*snaq1 + c(lambda.1_1.g1, lambda.1_1.g2)*snaq1
## F1 =~ c(NA, NA)*snaq2 + c(lambda.2_1.g1, lambda.2_1.g2)*snaq2
## F1 =~ c(NA, NA)*snaq3 + c(lambda.3_1.g1, lambda.3_1.g2)*snaq3
## F1 =~ c(NA, NA)*snaq4 + c(lambda.4_1.g1, lambda.4_1.g2)*snaq4
## F1 =~ c(NA, NA)*snaq5 + c(lambda.5_1.g1, lambda.5_1.g2)*snaq5
## F1 =~ c(NA, NA)*snaq6 + c(lambda.6_1.g1, lambda.6_1.g2)*snaq6
## F1 =~ c(NA, NA)*snaq7 + c(lambda.7_1.g1, lambda.7_1.g2)*snaq7
## F1 =~ c(NA, NA)*snaq8 + c(lambda.8_1.g1, lambda.8_1.g2)*snaq8
## F1 =~ c(NA, NA)*snaq9 + c(lambda.9_1.g1, lambda.9_1.g2)*snaq9
## 
## ## THRESHOLDS:
## 
## snaq1 | c(NA, NA)*t1 + c(snaq1.thr1.g1, snaq1.thr1.g2)*t1
## snaq1 | c(NA, NA)*t2 + c(snaq1.thr2.g1, snaq1.thr2.g2)*t2
## snaq1 | c(NA, NA)*t3 + c(snaq1.thr3.g1, snaq1.thr3.g2)*t3
## snaq2 | c(NA, NA)*t1 + c(snaq2.thr1.g1, snaq2.thr1.g2)*t1
## snaq2 | c(NA, NA)*t2 + c(snaq2.thr2.g1, snaq2.thr2.g2)*t2
## snaq2 | c(NA, NA)*t3 + c(snaq2.thr3.g1, snaq2.thr3.g2)*t3
## snaq3 | c(NA, NA)*t1 + c(snaq3.thr1.g1, snaq3.thr1.g2)*t1
## snaq3 | c(NA, NA)*t2 + c(snaq3.thr2.g1, snaq3.thr2.g2)*t2
## snaq3 | c(NA, NA)*t3 + c(snaq3.thr3.g1, snaq3.thr3.g2)*t3
## snaq4 | c(NA, NA)*t1 + c(snaq4.thr1.g1, snaq4.thr1.g2)*t1
## snaq4 | c(NA, NA)*t2 + c(snaq4.thr2.g1, snaq4.thr2.g2)*t2
## snaq4 | c(NA, NA)*t3 + c(snaq4.thr3.g1, snaq4.thr3.g2)*t3
## snaq5 | c(NA, NA)*t1 + c(snaq5.thr1.g1, snaq5.thr1.g2)*t1
## snaq5 | c(NA, NA)*t2 + c(snaq5.thr2.g1, snaq5.thr2.g2)*t2
## snaq5 | c(NA, NA)*t3 + c(snaq5.thr3.g1, snaq5.thr3.g2)*t3
## snaq6 | c(NA, NA)*t1 + c(snaq6.thr1.g1, snaq6.thr1.g2)*t1
## snaq6 | c(NA, NA)*t2 + c(snaq6.thr2.g1, snaq6.thr2.g2)*t2
## snaq6 | c(NA, NA)*t3 + c(snaq6.thr3.g1, snaq6.thr3.g2)*t3
## snaq7 | c(NA, NA)*t1 + c(snaq7.thr1.g1, snaq7.thr1.g2)*t1
## snaq7 | c(NA, NA)*t2 + c(snaq7.thr2.g1, snaq7.thr2.g2)*t2
## snaq7 | c(NA, NA)*t3 + c(snaq7.thr3.g1, snaq7.thr3.g2)*t3
## snaq8 | c(NA, NA)*t1 + c(snaq8.thr1.g1, snaq8.thr1.g2)*t1
## snaq8 | c(NA, NA)*t2 + c(snaq8.thr2.g1, snaq8.thr2.g2)*t2
## snaq8 | c(NA, NA)*t3 + c(snaq8.thr3.g1, snaq8.thr3.g2)*t3
## snaq9 | c(NA, NA)*t1 + c(snaq9.thr1.g1, snaq9.thr1.g2)*t1
## snaq9 | c(NA, NA)*t2 + c(snaq9.thr2.g1, snaq9.thr2.g2)*t2
## snaq9 | c(NA, NA)*t3 + c(snaq9.thr3.g1, snaq9.thr3.g2)*t3
## 
## ## INTERCEPTS:
## 
## snaq1 ~ c(0, 0)*1 + c(nu.1.g1, nu.1.g2)*1
## snaq2 ~ c(0, 0)*1 + c(nu.2.g1, nu.2.g2)*1
## snaq3 ~ c(0, 0)*1 + c(nu.3.g1, nu.3.g2)*1
## snaq4 ~ c(0, 0)*1 + c(nu.4.g1, nu.4.g2)*1
## snaq5 ~ c(0, 0)*1 + c(nu.5.g1, nu.5.g2)*1
## snaq6 ~ c(0, 0)*1 + c(nu.6.g1, nu.6.g2)*1
## snaq7 ~ c(0, 0)*1 + c(nu.7.g1, nu.7.g2)*1
## snaq8 ~ c(0, 0)*1 + c(nu.8.g1, nu.8.g2)*1
## snaq9 ~ c(0, 0)*1 + c(nu.9.g1, nu.9.g2)*1
## 
## ## SCALING FACTORS:
## 
## snaq1 ~*~ c(1, 1)*snaq1
## snaq2 ~*~ c(1, 1)*snaq2
## snaq3 ~*~ c(1, 1)*snaq3
## snaq4 ~*~ c(1, 1)*snaq4
## snaq5 ~*~ c(1, 1)*snaq5
## snaq6 ~*~ c(1, 1)*snaq6
## snaq7 ~*~ c(1, 1)*snaq7
## snaq8 ~*~ c(1, 1)*snaq8
## snaq9 ~*~ c(1, 1)*snaq9
## 
## 
## ## LATENT MEANS/INTERCEPTS:
## 
## F1 ~ c(0, 0)*1 + c(alpha.1.g1, alpha.1.g2)*1
## 
## ## COMMON-FACTOR VARIANCES:
## 
## F1 ~~ c(1, 1)*F1 + c(psi.1_1.g1, psi.1_1.g2)*F1
# Have to specify as.character to submit to lavaan
model.baseline <- as.character(baseline)
# Fitting baseline model in lavaan via cfa function
fit.baseline <- cfa(model.baseline, data = dat, group = "gender",
                    ordered = c("snaq1", "snaq2", "snaq3", "snaq4", "snaq5", "snaq6",
                                "snaq7", "snaq8", "snaq9"))



# Obtaining results from baseline model
summary(fit.baseline, fit.measures=TRUE, rsquare=TRUE)
## lavaan 0.6-11 ended normally after 16 iterations
## 
##   Estimator                                       DWLS
##   Optimization method                           NLMINB
##   Number of model parameters                        72
##                                                       
##   Number of observations per group:                   
##     male                                           177
##     female                                         180
##                                                       
## Model Test User Model:
##                                               Standard      Robust
##   Test Statistic                                87.526     150.014
##   Degrees of freedom                                54          54
##   P-value (Chi-square)                           0.003       0.000
##   Scaling correction factor                                  0.634
##   Shift parameter for each group:                                 
##       male                                                   5.889
##       female                                                 5.989
##        simple second-order correction                             
##   Test statistic for each group:
##     male                                        35.941      62.612
##     female                                      51.585      87.402
## 
## Model Test Baseline Model:
## 
##   Test statistic                              6242.272    3266.319
##   Degrees of freedom                                72          72
##   P-value                                        0.000       0.000
##   Scaling correction factor                                  1.932
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.995       0.970
##   Tucker-Lewis Index (TLI)                       0.993       0.960
##                                                                   
##   Robust Comparative Fit Index (CFI)                            NA
##   Robust Tucker-Lewis Index (TLI)                               NA
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.059       0.100
##   90 Percent confidence interval - lower         0.035       0.081
##   90 Percent confidence interval - upper         0.081       0.119
##   P-value RMSEA <= 0.05                          0.241       0.000
##                                                                   
##   Robust RMSEA                                                  NA
##   90 Percent confidence interval - lower                        NA
##   90 Percent confidence interval - upper                        NA
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.071       0.071
## 
## Parameter Estimates:
## 
##   Standard errors                           Robust.sem
##   Information                                 Expected
##   Information saturated (h1) model        Unstructured
## 
## 
## Group 1 [male]:
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   F1 =~                                               
##     snaq1   (l.1_)    0.792    0.037   21.425    0.000
##     snaq2   (l.2_)    0.880    0.028   31.364    0.000
##     snaq3   (l.3_)    0.793    0.040   19.979    0.000
##     snaq4   (l.4_)    0.924    0.026   35.328    0.000
##     snaq5   (l.5_)    0.657    0.054   12.122    0.000
##     snaq6   (l.6_)    0.534    0.066    8.121    0.000
##     snaq7   (l.7_)    0.757    0.045   16.743    0.000
##     snaq8   (l.8_)    0.808    0.041   19.841    0.000
##     snaq9   (l.9_)    0.832    0.060   13.924    0.000
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .snaq1   (n.1.)    0.000                           
##    .snaq2   (n.2.)    0.000                           
##    .snaq3   (n.3.)    0.000                           
##    .snaq4   (n.4.)    0.000                           
##    .snaq5   (n.5.)    0.000                           
##    .snaq6   (n.6.)    0.000                           
##    .snaq7   (n.7.)    0.000                           
##    .snaq8   (n.8.)    0.000                           
##    .snaq9   (n.9.)    0.000                           
##     F1      (a.1.)    0.000                           
## 
## Thresholds:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     snq1|t1 (s1.1)   -0.092    0.095   -0.974    0.330
##     snq1|t2 (s1.2)    0.790    0.106    7.454    0.000
##     snq1|t3 (s1.3)    1.182    0.123    9.622    0.000
##     snq2|t1 (s2.1)    0.121    0.095    1.274    0.203
##     snq2|t2 (s2.2)    1.100    0.118    9.286    0.000
##     snq2|t3 (s2.3)    1.338    0.133   10.088    0.000
##     snq3|t1 (s3.1)    0.280    0.096    2.919    0.004
##     snq3|t2 (s3.2)    1.182    0.123    9.622    0.000
##     snq3|t3 (s3.3)    1.374    0.135   10.161    0.000
##     snq4|t1 (s4.1)    0.559    0.100    5.587    0.000
##     snq4|t2 (s4.2)    1.338    0.133   10.088    0.000
##     snq4|t3 (s4.3)    1.693    0.165   10.287    0.000
##     snq5|t1 (s5.1)    0.178    0.095    1.873    0.061
##     snq5|t2 (s5.2)    1.374    0.135   10.161    0.000
##     snq5|t3 (s5.3)    1.693    0.165   10.287    0.000
##     snq6|t1 (s6.1)   -0.733    0.104   -7.032    0.000
##     snq6|t2 (s6.2)    0.679    0.103    6.603    0.000
##     snq6|t3 (s6.3)    1.154    0.121    9.513    0.000
##     snq7|t1 (s7.1)    0.050    0.095    0.525    0.600
##     snq7|t2 (s7.2)    0.979    0.113    8.669    0.000
##     snq7|t3 (s7.3)    1.537    0.149   10.342    0.000
##     snq8|t1 (s8.1)    0.236    0.095    2.471    0.013
##     snq8|t2 (s8.2)    1.272    0.128    9.919    0.000
##     snq8|t3 (s8.3)    1.492    0.145   10.317    0.000
##     snq9|t1 (s9.1)    1.272    0.128    9.919    0.000
##     snq9|t2 (s9.2)    1.907    0.193    9.884    0.000
##     snq9|t3 (s9.3)    2.003    0.209    9.599    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     F1      (p.1_)    1.000                           
##    .snaq1             0.372                           
##    .snaq2             0.226                           
##    .snaq3             0.371                           
##    .snaq4             0.147                           
##    .snaq5             0.568                           
##    .snaq6             0.714                           
##    .snaq7             0.427                           
##    .snaq8             0.347                           
##    .snaq9             0.307                           
## 
## Scales y*:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     snaq1             1.000                           
##     snaq2             1.000                           
##     snaq3             1.000                           
##     snaq4             1.000                           
##     snaq5             1.000                           
##     snaq6             1.000                           
##     snaq7             1.000                           
##     snaq8             1.000                           
##     snaq9             1.000                           
## 
## R-Square:
##                    Estimate
##     snaq1             0.628
##     snaq2             0.774
##     snaq3             0.629
##     snaq4             0.853
##     snaq5             0.432
##     snaq6             0.286
##     snaq7             0.573
##     snaq8             0.653
##     snaq9             0.693
## 
## 
## Group 2 [female]:
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   F1 =~                                               
##     snaq1   (l.1_)    0.734    0.048   15.427    0.000
##     snaq2   (l.2_)    0.765    0.042   18.168    0.000
##     snaq3   (l.3_)    0.782    0.039   20.087    0.000
##     snaq4   (l.4_)    0.857    0.042   20.633    0.000
##     snaq5   (l.5_)    0.711    0.049   14.369    0.000
##     snaq6   (l.6_)    0.601    0.056   10.821    0.000
##     snaq7   (l.7_)    0.861    0.034   25.261    0.000
##     snaq8   (l.8_)    0.786    0.040   19.653    0.000
##     snaq9   (l.9_)    0.738    0.094    7.889    0.000
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .snaq1   (n.1.)    0.000                           
##    .snaq2   (n.2.)    0.000                           
##    .snaq3   (n.3.)    0.000                           
##    .snaq4   (n.4.)    0.000                           
##    .snaq5   (n.5.)    0.000                           
##    .snaq6   (n.6.)    0.000                           
##    .snaq7   (n.7.)    0.000                           
##    .snaq8   (n.8.)    0.000                           
##    .snaq9   (n.9.)    0.000                           
##     F1      (a.1.)    0.000                           
## 
## Thresholds:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     snq1|t1 (s1.1)   -0.154    0.094   -1.635    0.102
##     snq1|t2 (s1.2)    0.924    0.110    8.417    0.000
##     snq1|t3 (s1.3)    1.221    0.124    9.842    0.000
##     snq2|t1 (s2.1)    0.182    0.094    1.931    0.053
##     snq2|t2 (s2.2)    1.348    0.132   10.194    0.000
##     snq2|t3 (s2.3)    1.764    0.172   10.279    0.000
##     snq3|t1 (s3.1)    0.028    0.094    0.297    0.766
##     snq3|t2 (s3.2)    1.164    0.121    9.634    0.000
##     snq3|t3 (s3.3)    1.593    0.153   10.434    0.000
##     snq4|t1 (s4.1)    0.746    0.104    7.189    0.000
##     snq4|t2 (s4.2)    1.593    0.153   10.434    0.000
##     snq4|t3 (s4.3)    2.287    0.268    8.526    0.000
##     snq5|t1 (s5.1)    0.140    0.094    1.486    0.137
##     snq5|t2 (s5.2)    1.282    0.128   10.030    0.000
##     snq5|t3 (s5.3)    1.915    0.192    9.948    0.000
##     snq6|t1 (s6.1)   -0.924    0.110   -8.417    0.000
##     snq6|t2 (s6.2)    1.164    0.121    9.634    0.000
##     snq6|t3 (s6.3)    1.764    0.172   10.279    0.000
##     snq7|t1 (s7.1)   -0.056    0.094   -0.595    0.552
##     snq7|t2 (s7.2)    1.192    0.122    9.740    0.000
##     snq7|t3 (s7.3)    1.593    0.153   10.434    0.000
##     snq8|t1 (s8.1)    0.182    0.094    1.931    0.053
##     snq8|t2 (s8.2)    1.314    0.130   10.115    0.000
##     snq8|t3 (s8.3)    2.010    0.208    9.656    0.000
##     snq9|t1 (s9.1)    1.383    0.135   10.264    0.000
##     snq9|t2 (s9.2)    2.128    0.231    9.219    0.000
##     snq9|t3 (s9.3)    2.539    0.350    7.258    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     F1      (p.1_)    1.000                           
##    .snaq1             0.461                           
##    .snaq2             0.414                           
##    .snaq3             0.389                           
##    .snaq4             0.266                           
##    .snaq5             0.494                           
##    .snaq6             0.638                           
##    .snaq7             0.258                           
##    .snaq8             0.382                           
##    .snaq9             0.455                           
## 
## Scales y*:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     snaq1             1.000                           
##     snaq2             1.000                           
##     snaq3             1.000                           
##     snaq4             1.000                           
##     snaq5             1.000                           
##     snaq6             1.000                           
##     snaq7             1.000                           
##     snaq8             1.000                           
##     snaq9             1.000                           
## 
## R-Square:
##                    Estimate
##     snaq1             0.539
##     snaq2             0.586
##     snaq3             0.611
##     snaq4             0.734
##     snaq5             0.506
##     snaq6             0.362
##     snaq7             0.742
##     snaq8             0.618
##     snaq9             0.545
# Extracting fit indices into the first row of all results matrix
all.results[1,]<-round(data.matrix(fitmeasures(fit.baseline,
                                               fit.measures = c("chisq.scaled","df.scaled","pvalue.scaled",
                                                                "cfi.scaled", "tli.scaled", "rmsea.scaled"))),
                                                                digits=3)

Equal Thresholds Model

# To remain consistent with Wu and Estabrook's (2016) notation, we call this step as prop4 to indicate the
# alignment with Proposition 4 in Wu and Estabrook's article.
prop4 <- measEq.syntax(configural.model = mod.cat,
                       data = dat,
                       ordered = c("snaq1", "snaq2", "snaq3", "snaq4", "snaq5", "snaq6",
                                   "snaq7", "snaq8", "snaq9"),
                                   parameterization = "delta",
                                   ID.fac = "std.lv",
                                   ID.cat = "Wu.Estabrook.2016",
                                   group = "gender",
                                   group.equal = c("thresholds"))

model.prop4 <- as.character(prop4)

# Fitting thresholds invariance model in lavaan via cfa function
fit.prop4 <- cfa(model.prop4, data = dat, group ="gender",
                 ordered = c("snaq1", "snaq2", "snaq3", "snaq4", "snaq5", "snaq6",
                             "snaq7", "snaq8", "snaq9"))

# Obtaining results from thresholds invariance model
summary(fit.prop4, fit.measures=TRUE, rsquare=TRUE)
## lavaan 0.6-11 ended normally after 81 iterations
## 
##   Estimator                                       DWLS
##   Optimization method                           NLMINB
##   Number of model parameters                        90
##   Number of equality constraints                    27
##                                                       
##   Number of observations per group:                   
##     male                                           177
##     female                                         180
##                                                       
## Model Test User Model:
##                                               Standard      Robust
##   Test Statistic                                92.166     165.795
##   Degrees of freedom                                63          63
##   P-value (Chi-square)                           0.010       0.000
##   Scaling correction factor                                  0.609
##   Shift parameter for each group:                                 
##       male                                                   7.136
##       female                                                 7.257
##        simple second-order correction                             
##   Test statistic for each group:
##     male                                        38.429      70.264
##     female                                      53.737      95.531
## 
## Model Test Baseline Model:
## 
##   Test statistic                              6242.272    3266.319
##   Degrees of freedom                                72          72
##   P-value                                        0.000       0.000
##   Scaling correction factor                                  1.932
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.995       0.968
##   Tucker-Lewis Index (TLI)                       0.995       0.963
##                                                                   
##   Robust Comparative Fit Index (CFI)                            NA
##   Robust Tucker-Lewis Index (TLI)                               NA
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.051       0.096
##   90 Percent confidence interval - lower         0.026       0.078
##   90 Percent confidence interval - upper         0.073       0.114
##   P-value RMSEA <= 0.05                          0.448       0.000
##                                                                   
##   Robust RMSEA                                                  NA
##   90 Percent confidence interval - lower                        NA
##   90 Percent confidence interval - upper                        NA
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.071       0.071
## 
## Parameter Estimates:
## 
##   Standard errors                           Robust.sem
##   Information                                 Expected
##   Information saturated (h1) model        Unstructured
## 
## 
## Group 1 [male]:
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   F1 =~                                               
##     snaq1   (l.1_)    0.792    0.037   21.425    0.000
##     snaq2   (l.2_)    0.880    0.028   31.364    0.000
##     snaq3   (l.3_)    0.793    0.040   19.979    0.000
##     snaq4   (l.4_)    0.924    0.026   35.328    0.000
##     snaq5   (l.5_)    0.657    0.054   12.122    0.000
##     snaq6   (l.6_)    0.534    0.066    8.121    0.000
##     snaq7   (l.7_)    0.757    0.045   16.743    0.000
##     snaq8   (l.8_)    0.808    0.041   19.841    0.000
##     snaq9   (l.9_)    0.832    0.060   13.924    0.000
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .snaq1   (n.1.)    0.000                           
##    .snaq2   (n.2.)    0.000                           
##    .snaq3   (n.3.)    0.000                           
##    .snaq4   (n.4.)    0.000                           
##    .snaq5   (n.5.)    0.000                           
##    .snaq6   (n.6.)    0.000                           
##    .snaq7   (n.7.)    0.000                           
##    .snaq8   (n.8.)    0.000                           
##    .snaq9   (n.9.)    0.000                           
##     F1      (a.1.)    0.000                           
## 
## Thresholds:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     snq1|t1 (s1.1)   -0.099    0.095   -1.051    0.293
##     snq1|t2 (s1.2)    0.825    0.101    8.182    0.000
##     snq1|t3 (s1.3)    1.147    0.119    9.662    0.000
##     snq2|t1 (s2.1)    0.124    0.095    1.314    0.189
##     snq2|t2 (s2.2)    1.076    0.115    9.360    0.000
##     snq2|t3 (s2.3)    1.362    0.132   10.303    0.000
##     snq3|t1 (s3.1)    0.286    0.096    2.985    0.003
##     snq3|t2 (s3.2)    1.142    0.118    9.670    0.000
##     snq3|t3 (s3.3)    1.411    0.136   10.379    0.000
##     snq4|t1 (s4.1)    0.569    0.100    5.682    0.000
##     snq4|t2 (s4.2)    1.293    0.128   10.091    0.000
##     snq4|t3 (s4.3)    1.737    0.164   10.591    0.000
##     snq5|t1 (s5.1)    0.187    0.095    1.976    0.048
##     snq5|t2 (s5.2)    1.308    0.130   10.035    0.000
##     snq5|t3 (s5.3)    1.763    0.163   10.819    0.000
##     snq6|t1 (s6.1)   -0.737    0.104   -7.101    0.000
##     snq6|t2 (s6.2)    0.695    0.098    7.058    0.000
##     snq6|t3 (s6.3)    1.137    0.117    9.692    0.000
##     snq7|t1 (s7.1)    0.038    0.095    0.399    0.690
##     snq7|t2 (s7.2)    1.035    0.107    9.663    0.000
##     snq7|t3 (s7.3)    1.470    0.139   10.595    0.000
##     snq8|t1 (s8.1)    0.250    0.095    2.626    0.009
##     snq8|t2 (s8.2)    1.185    0.124    9.548    0.000
##     snq8|t3 (s8.3)    1.571    0.147   10.670    0.000
##     snq9|t1 (s9.1)    1.278    0.128    9.965    0.000
##     snq9|t2 (s9.2)    1.856    0.190    9.788    0.000
##     snq9|t3 (s9.3)    2.048    0.214    9.559    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     F1      (p.1_)    1.000                           
##    .snaq1             0.372                           
##    .snaq2             0.226                           
##    .snaq3             0.371                           
##    .snaq4             0.147                           
##    .snaq5             0.568                           
##    .snaq6             0.714                           
##    .snaq7             0.427                           
##    .snaq8             0.347                           
##    .snaq9             0.307                           
## 
## Scales y*:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     snaq1             1.000                           
##     snaq2             1.000                           
##     snaq3             1.000                           
##     snaq4             1.000                           
##     snaq5             1.000                           
##     snaq6             1.000                           
##     snaq7             1.000                           
##     snaq8             1.000                           
##     snaq9             1.000                           
## 
## R-Square:
##                    Estimate
##     snaq1             0.628
##     snaq2             0.774
##     snaq3             0.629
##     snaq4             0.853
##     snaq5             0.432
##     snaq6             0.286
##     snaq7             0.573
##     snaq8             0.653
##     snaq9             0.693
## 
## 
## Group 2 [female]:
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   F1 =~                                               
##     snaq1   (l.1_)    0.654    0.107    6.131    0.000
##     snaq2   (l.2_)    0.610    0.104    5.878    0.000
##     snaq3   (l.3_)    0.573    0.099    5.786    0.000
##     snaq4   (l.4_)    0.687    0.153    4.478    0.000
##     snaq5   (l.5_)    0.661    0.118    5.602    0.000
##     snaq6   (l.6_)    0.416    0.062    6.697    0.000
##     snaq7   (l.7_)    0.726    0.109    6.662    0.000
##     snaq8   (l.8_)    0.606    0.103    5.867    0.000
##     snaq9   (l.9_)    0.532    0.223    2.384    0.017
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .snaq1   (n.1.)    0.032    0.119    0.269    0.788
##    .snaq2   (n.2.)   -0.018    0.130   -0.142    0.887
##    .snaq3   (n.3.)    0.268    0.119    2.261    0.024
##    .snaq4   (n.4.)   -0.023    0.201   -0.114    0.909
##    .snaq5   (n.5.)    0.065    0.136    0.481    0.631
##    .snaq6   (n.6.)   -0.100    0.103   -0.969    0.333
##    .snaq7   (n.7.)    0.076    0.121    0.632    0.527
##    .snaq8   (n.8.)    0.118    0.128    0.918    0.358
##    .snaq9   (n.9.)    0.285    0.424    0.671    0.502
##     F1      (a.1.)    0.000                           
## 
## Thresholds:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     snq1|t1 (s1.1)   -0.099    0.095   -1.051    0.293
##     snq1|t2 (s1.2)    0.825    0.101    8.182    0.000
##     snq1|t3 (s1.3)    1.147    0.119    9.662    0.000
##     snq2|t1 (s2.1)    0.124    0.095    1.314    0.189
##     snq2|t2 (s2.2)    1.076    0.115    9.360    0.000
##     snq2|t3 (s2.3)    1.362    0.132   10.303    0.000
##     snq3|t1 (s3.1)    0.286    0.096    2.985    0.003
##     snq3|t2 (s3.2)    1.142    0.118    9.670    0.000
##     snq3|t3 (s3.3)    1.411    0.136   10.379    0.000
##     snq4|t1 (s4.1)    0.569    0.100    5.682    0.000
##     snq4|t2 (s4.2)    1.293    0.128   10.091    0.000
##     snq4|t3 (s4.3)    1.737    0.164   10.591    0.000
##     snq5|t1 (s5.1)    0.187    0.095    1.976    0.048
##     snq5|t2 (s5.2)    1.308    0.130   10.035    0.000
##     snq5|t3 (s5.3)    1.763    0.163   10.819    0.000
##     snq6|t1 (s6.1)   -0.737    0.104   -7.101    0.000
##     snq6|t2 (s6.2)    0.695    0.098    7.058    0.000
##     snq6|t3 (s6.3)    1.137    0.117    9.692    0.000
##     snq7|t1 (s7.1)    0.038    0.095    0.399    0.690
##     snq7|t2 (s7.2)    1.035    0.107    9.663    0.000
##     snq7|t3 (s7.3)    1.470    0.139   10.595    0.000
##     snq8|t1 (s8.1)    0.250    0.095    2.626    0.009
##     snq8|t2 (s8.2)    1.185    0.124    9.548    0.000
##     snq8|t3 (s8.3)    1.571    0.147   10.670    0.000
##     snq9|t1 (s9.1)    1.278    0.128    9.965    0.000
##     snq9|t2 (s9.2)    1.856    0.190    9.788    0.000
##     snq9|t3 (s9.3)    2.048    0.214    9.559    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     F1      (p.1_)    1.000                           
##    .snaq1             0.366                           
##    .snaq2             0.263                           
##    .snaq3             0.209                           
##    .snaq4             0.171                           
##    .snaq5             0.427                           
##    .snaq6             0.306                           
##    .snaq7             0.183                           
##    .snaq8             0.227                           
##    .snaq9             0.236                           
## 
## Scales y*:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     snaq1             1.123    0.151    7.431    0.000
##     snaq2             1.254    0.186    6.740    0.000
##     snaq3             1.364    0.208    6.568    0.000
##     snaq4             1.247    0.252    4.956    0.000
##     snaq5             1.076    0.155    6.949    0.000
##     snaq6             1.445    0.144   10.054    0.000
##     snaq7             1.186    0.158    7.505    0.000
##     snaq8             1.297    0.198    6.545    0.000
##     snaq9             1.388    0.509    2.728    0.006
## 
## R-Square:
##                    Estimate
##     snaq1             0.539
##     snaq2             0.586
##     snaq3             0.611
##     snaq4             0.734
##     snaq5             0.506
##     snaq6             0.362
##     snaq7             0.742
##     snaq8             0.618
##     snaq9             0.545
# Extracting fit indices into the second row of all results matrix
all.results[2,]<-round(data.matrix(fitmeasures(fit.prop4,
                                          fit.measures = c("chisq.scaled","df.scaled","pvalue.scaled",
                                                            "cfi.scaled", "tli.scaled", "rmsea.scaled"))),
                       digits=3)

lavTestLRT(fit.baseline,fit.prop4)
## Scaled Chi-Squared Difference Test (method = "satorra.2000")
## 
## lavaan NOTE:
##     The "Chisq" column contains standard test statistics, not the
##     robust test that should be reported per model. A robust difference
##     test is a function of two standard (not robust) statistics.
##  
##              Df AIC BIC  Chisq Chisq diff Df diff Pr(>Chisq)
## fit.baseline 54         87.526                              
## fit.prop4    63         92.166     13.514       9     0.1407

Equal Thresholds and Loadings Model

# Proposition 7 per Wu and Estabrook (2016)
prop7 <- measEq.syntax(configural.model = mod.cat,
                       data = dat,
                       ordered = c("snaq1", "snaq2", "snaq3", "snaq4", "snaq5", "snaq6",
                                   "snaq7", "snaq8", "snaq9"),
                                   parameterization = "delta",
                                   ID.fac = "std.lv",
                                   ID.cat = "Wu.Estabrook.2016",
                                   group = "gender",
                                   group.equal = c("thresholds","loadings"))
model.prop7 <- as.character(prop7)

fit.prop7 <- cfa(model.prop7, data = dat, group =
                   "gender", ordered = c("snaq1", "snaq2", "snaq3", "snaq4", "snaq5", "snaq6",
                                         "snaq7", "snaq8", "snaq9"))
summary(fit.prop7, fit.measures=TRUE, rsquare=TRUE)
## lavaan 0.6-11 ended normally after 42 iterations
## 
##   Estimator                                       DWLS
##   Optimization method                           NLMINB
##   Number of model parameters                        91
##   Number of equality constraints                    36
##                                                       
##   Number of observations per group:                   
##     male                                           177
##     female                                         180
##                                                       
## Model Test User Model:
##                                               Standard      Robust
##   Test Statistic                                97.147     158.085
##   Degrees of freedom                                71          71
##   P-value (Chi-square)                           0.021       0.000
##   Scaling correction factor                                  0.678
##   Shift parameter for each group:                                 
##       male                                                   7.378
##       female                                                 7.503
##        simple second-order correction                             
##   Test statistic for each group:
##     male                                        41.211      68.128
##     female                                      55.935      89.957
## 
## Model Test Baseline Model:
## 
##   Test statistic                              6242.272    3266.319
##   Degrees of freedom                                72          72
##   P-value                                        0.000       0.000
##   Scaling correction factor                                  1.932
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.996       0.973
##   Tucker-Lewis Index (TLI)                       0.996       0.972
##                                                                   
##   Robust Comparative Fit Index (CFI)                            NA
##   Robust Tucker-Lewis Index (TLI)                               NA
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.046       0.083
##   90 Percent confidence interval - lower         0.018       0.066
##   90 Percent confidence interval - upper         0.067       0.101
##   P-value RMSEA <= 0.05                          0.611       0.001
##                                                                   
##   Robust RMSEA                                                  NA
##   90 Percent confidence interval - lower                        NA
##   90 Percent confidence interval - upper                        NA
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.071       0.071
## 
## Parameter Estimates:
## 
##   Standard errors                           Robust.sem
##   Information                                 Expected
##   Information saturated (h1) model        Unstructured
## 
## 
## Group 1 [male]:
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   F1 =~                                               
##     snaq1   (l.1_)    0.794    0.035   22.636    0.000
##     snaq2   (l.2_)    0.875    0.028   31.712    0.000
##     snaq3   (l.3_)    0.789    0.038   20.677    0.000
##     snaq4   (l.4_)    0.924    0.025   36.244    0.000
##     snaq5   (l.5_)    0.671    0.051   13.266    0.000
##     snaq6   (l.6_)    0.530    0.057    9.296    0.000
##     snaq7   (l.7_)    0.765    0.043   17.918    0.000
##     snaq8   (l.8_)    0.805    0.039   20.495    0.000
##     snaq9   (l.9_)    0.830    0.059   14.094    0.000
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .snaq1   (n.1.)    0.000                           
##    .snaq2   (n.2.)    0.000                           
##    .snaq3   (n.3.)    0.000                           
##    .snaq4   (n.4.)    0.000                           
##    .snaq5   (n.5.)    0.000                           
##    .snaq6   (n.6.)    0.000                           
##    .snaq7   (n.7.)    0.000                           
##    .snaq8   (n.8.)    0.000                           
##    .snaq9   (n.9.)    0.000                           
##     F1      (a.1.)    0.000                           
## 
## Thresholds:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     snq1|t1 (s1.1)   -0.092    0.093   -0.984    0.325
##     snq1|t2 (s1.2)    0.822    0.098    8.363    0.000
##     snq1|t3 (s1.3)    1.139    0.113   10.080    0.000
##     snq2|t1 (s2.1)    0.083    0.094    0.886    0.376
##     snq2|t2 (s2.2)    1.107    0.112    9.865    0.000
##     snq2|t3 (s2.3)    1.403    0.128   10.935    0.000
##     snq3|t1 (s3.1)    0.254    0.091    2.783    0.005
##     snq3|t2 (s3.2)    1.166    0.114   10.229    0.000
##     snq3|t3 (s3.3)    1.444    0.131   11.047    0.000
##     snq4|t1 (s4.1)    0.552    0.098    5.663    0.000
##     snq4|t2 (s4.2)    1.309    0.123   10.629    0.000
##     snq4|t3 (s4.3)    1.756    0.159   11.063    0.000
##     snq5|t1 (s5.1)    0.249    0.094    2.646    0.008
##     snq5|t2 (s5.2)    1.242    0.120   10.381    0.000
##     snq5|t3 (s5.3)    1.675    0.155   10.780    0.000
##     snq6|t1 (s6.1)   -0.749    0.104   -7.235    0.000
##     snq6|t2 (s6.2)    0.700    0.098    7.149    0.000
##     snq6|t3 (s6.3)    1.145    0.116    9.857    0.000
##     snq7|t1 (s7.1)    0.093    0.092    1.012    0.311
##     snq7|t2 (s7.2)    1.002    0.102    9.854    0.000
##     snq7|t3 (s7.3)    1.388    0.127   10.890    0.000
##     snq8|t1 (s8.1)    0.231    0.092    2.495    0.013
##     snq8|t2 (s8.2)    1.204    0.120   10.070    0.000
##     snq8|t3 (s8.3)    1.591    0.141   11.301    0.000
##     snq9|t1 (s9.1)    1.245    0.131    9.482    0.000
##     snq9|t2 (s9.2)    1.900    0.176   10.771    0.000
##     snq9|t3 (s9.3)    2.082    0.200   10.390    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     F1      (p.1_)    1.000                           
##    .snaq1             0.370                           
##    .snaq2             0.235                           
##    .snaq3             0.378                           
##    .snaq4             0.146                           
##    .snaq5             0.550                           
##    .snaq6             0.719                           
##    .snaq7             0.414                           
##    .snaq8             0.351                           
##    .snaq9             0.311                           
## 
## Scales y*:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     snaq1             1.000                           
##     snaq2             1.000                           
##     snaq3             1.000                           
##     snaq4             1.000                           
##     snaq5             1.000                           
##     snaq6             1.000                           
##     snaq7             1.000                           
##     snaq8             1.000                           
##     snaq9             1.000                           
## 
## R-Square:
##                    Estimate
##     snaq1             0.630
##     snaq2             0.765
##     snaq3             0.622
##     snaq4             0.854
##     snaq5             0.450
##     snaq6             0.281
##     snaq7             0.586
##     snaq8             0.649
##     snaq9             0.689
## 
## 
## Group 2 [female]:
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   F1 =~                                               
##     snaq1   (l.1_)    0.794    0.035   22.636    0.000
##     snaq2   (l.2_)    0.875    0.028   31.712    0.000
##     snaq3   (l.3_)    0.789    0.038   20.677    0.000
##     snaq4   (l.4_)    0.924    0.025   36.244    0.000
##     snaq5   (l.5_)    0.671    0.051   13.266    0.000
##     snaq6   (l.6_)    0.530    0.057    9.296    0.000
##     snaq7   (l.7_)    0.765    0.043   17.918    0.000
##     snaq8   (l.8_)    0.805    0.039   20.495    0.000
##     snaq9   (l.9_)    0.830    0.059   14.094    0.000
## 
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .snaq1   (n.1.)    0.042    0.117    0.363    0.716
##    .snaq2   (n.2.)   -0.114    0.130   -0.879    0.379
##    .snaq3   (n.3.)    0.216    0.112    1.933    0.053
##    .snaq4   (n.4.)   -0.099    0.148   -0.673    0.501
##    .snaq5   (n.5.)    0.183    0.118    1.546    0.122
##    .snaq6   (n.6.)   -0.106    0.104   -1.020    0.308
##    .snaq7   (n.7.)    0.156    0.110    1.425    0.154
##    .snaq8   (n.8.)    0.078    0.119    0.653    0.514
##    .snaq9   (n.9.)   -0.026    0.297   -0.089    0.929
##     F1      (a.1.)    0.000                           
## 
## Thresholds:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     snq1|t1 (s1.1)   -0.092    0.093   -0.984    0.325
##     snq1|t2 (s1.2)    0.822    0.098    8.363    0.000
##     snq1|t3 (s1.3)    1.139    0.113   10.080    0.000
##     snq2|t1 (s2.1)    0.083    0.094    0.886    0.376
##     snq2|t2 (s2.2)    1.107    0.112    9.865    0.000
##     snq2|t3 (s2.3)    1.403    0.128   10.935    0.000
##     snq3|t1 (s3.1)    0.254    0.091    2.783    0.005
##     snq3|t2 (s3.2)    1.166    0.114   10.229    0.000
##     snq3|t3 (s3.3)    1.444    0.131   11.047    0.000
##     snq4|t1 (s4.1)    0.552    0.098    5.663    0.000
##     snq4|t2 (s4.2)    1.309    0.123   10.629    0.000
##     snq4|t3 (s4.3)    1.756    0.159   11.063    0.000
##     snq5|t1 (s5.1)    0.249    0.094    2.646    0.008
##     snq5|t2 (s5.2)    1.242    0.120   10.381    0.000
##     snq5|t3 (s5.3)    1.675    0.155   10.780    0.000
##     snq6|t1 (s6.1)   -0.749    0.104   -7.235    0.000
##     snq6|t2 (s6.2)    0.700    0.098    7.149    0.000
##     snq6|t3 (s6.3)    1.145    0.116    9.857    0.000
##     snq7|t1 (s7.1)    0.093    0.092    1.012    0.311
##     snq7|t2 (s7.2)    1.002    0.102    9.854    0.000
##     snq7|t3 (s7.3)    1.388    0.127   10.890    0.000
##     snq8|t1 (s8.1)    0.231    0.092    2.495    0.013
##     snq8|t2 (s8.2)    1.204    0.120   10.070    0.000
##     snq8|t3 (s8.3)    1.591    0.141   11.301    0.000
##     snq9|t1 (s9.1)    1.245    0.131    9.482    0.000
##     snq9|t2 (s9.2)    1.900    0.176   10.771    0.000
##     snq9|t3 (s9.3)    2.082    0.200   10.390    0.000
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     F1      (p.1_)    0.648    0.150    4.333    0.000
##    .snaq1             0.351                           
##    .snaq2             0.335                           
##    .snaq3             0.248                           
##    .snaq4             0.197                           
##    .snaq5             0.304                           
##    .snaq6             0.315                           
##    .snaq7             0.142                           
##    .snaq8             0.255                           
##    .snaq9             0.367                           
## 
## Scales y*:
##                    Estimate  Std.Err  z-value  P(>|z|)
##     snaq1             1.148    0.118    9.691    0.000
##     snaq2             1.097    0.112    9.761    0.000
##     snaq3             1.239    0.123   10.034    0.000
##     snaq4             1.155    0.129    8.965    0.000
##     snaq5             1.296    0.137    9.464    0.000
##     snaq6             1.419    0.131   10.848    0.000
##     snaq7             1.385    0.138   10.045    0.000
##     snaq8             1.217    0.130    9.334    0.000
##     snaq9             1.109    0.175    6.334    0.000
## 
## R-Square:
##                    Estimate
##     snaq1             0.538
##     snaq2             0.597
##     snaq3             0.619
##     snaq4             0.738
##     snaq5             0.489
##     snaq6             0.366
##     snaq7             0.728
##     snaq8             0.623
##     snaq9             0.549
# Extracting fit indices into the third row of all results matrix
all.results[3,]<-round(data.matrix(fitmeasures(fit.prop7,
                                               fit.measures = c("chisq.scaled","df.scaled","pvalue.scaled",
                                                                "cfi.scaled", "tli.scaled", "rmsea.scaled"))), digits=3)
lavTestLRT(fit.prop4, fit.prop7)
## Scaled Chi-Squared Difference Test (method = "satorra.2000")
## 
## lavaan NOTE:
##     The "Chisq" column contains standard test statistics, not the
##     robust test that should be reported per model. A robust difference
##     test is a function of two standard (not robust) statistics.
##  
##           Df AIC BIC  Chisq Chisq diff Df diff Pr(>Chisq)
## fit.prop4 63         92.166                              
## fit.prop7 71         97.147     5.1401       8     0.7425

Summary of Results

Indices of Fit

summary_sem <- cbind(c("Baseline", "Equal Tresholds", "Equal Tresholds and Loadings"), all.results) %>%
  as.data.frame() %>% 
  rename(Model = V1,
         WLSMV_based_chi2 = V2,
         df = V3,
         p_value = V4,
         CFI = V5,
         TLI = V6,
         RMSEA = V7
         ) 

summary_sem %>% 
  knitr::kable()
Model WLSMV_based_chi2 df p_value CFI TLI RMSEA
Baseline 150.014 54 0 0.97 0.96 0.1
Equal Tresholds 165.795 63 0 0.968 0.963 0.096
Equal Tresholds and Loadings 158.085 71 0 0.973 0.972 0.083

Model comparison

summary_2 <- rbind(
(lavTestLRT(fit.baseline,fit.prop4)),
(lavTestLRT(fit.prop4, fit.prop7))
)[c(1,2,4), c(5,6,7)] %>%
  as_tibble() %>%
  round(., 3)

summary_2 %>% 
  knitr::kable()
Chisq diff Df diff Pr(>Chisq)
NA NA NA
13.514 9 0.141
5.140 8 0.743

Merge results

options(knitr.kable.NA = '-') # This command allows to report "-" instead of "NA" when printing a Kable

cbind(summary_sem, summary_2) %>% 
    knitr::kable(., "simple",
               col.names = c("Model", "WLSMV-based$\\chi^2$", "$df$", "$p$", "CFI", "TLI", "RMSEA",
                             "Scaled-$\\Delta\\chi^2$", "$\\Delta df$", "$p$"),
               align = c("l", "c", "c", "c", "c", "c", "c", "c", "c", "c"))
Model WLSMV-based\(\chi^2\) \(df\) \(p\) CFI TLI RMSEA Scaled-\(\Delta\chi^2\) \(\Delta df\) \(p\)
Baseline 150.014 54 0 0.97 0.96 0.1 - - -
Equal Tresholds 165.795 63 0 0.968 0.963 0.096 13.514 9 0.141
Equal Tresholds and Loadings 158.085 71 0 0.973 0.972 0.083 5.140 8 0.743

  1. See DiStefano, Shi, & Morgan (2021) and Rutkowski, Svetina, & Liaw (2019) for a debate on this issue↩︎