Structural Equation Modeling | Exercises 4

1 What we are going to cover

  1. Ex.1 – Non-normal continuous data
  2. Ex.2 – Categorical data multi-group

2 Data

The data set used throughout is the European Social Survey ESS4-2008 Edition 4.5 was released on 1 December 2018. We will restrict the analysis to the Belgian case. Each line in the data set represents a Belgian respondent. The full dataset an documentation can be found on the ESS website

Codebook:

  • gvslvol Standard of living for the old, governments’ responsibility (0 Not governments’ responsibility at all - 10 Entirely governments’ responsibility)

  • gvslvue Standard of living for the unemployed, governments’ responsibility (0 Not governments’ responsibility at all - 10 Entirely governments’ responsibility)

  • gvhlthc Health care for the sick, governments’ responsibility (0 Not governments’ responsibility at all - 10 Entirely governments’ responsibility)

  • gvcldcr Child care services for working parents, governments’ responsibility (0 Not governments’ responsibility at all - 10 Entirely governments’ responsibility)

  • gvjbevn Job for everyone, governments’ responsibility (0 Not governments’ responsibility at all - 10 Entirely governments’ responsibility)

  • gvpdlwk Paid leave from work to care for sick family, governments’ responsibility (0 Not governments’ responsibility at all - 10 Entirely governments’ responsibility)

  • sbstrec Social benefits/services place too great strain on economy (1 Agree strongly - 5 Disagree strongly)

  • sbbsntx Social benefits/services cost businesses too much in taxes/charges (1 Agree strongly - 5 Disagree strongly)

  • sbprvpv Social benefits/services prevent widespread poverty (1 Agree strongly - 5 Disagree strongly)

  • sbeqsoc Social benefits/services lead to a more equal society (1 Agree strongly - 5 Disagree strongly)

  • sbcwkfm Social benefits/services make it easier to combine work and family (1 Agree strongly - 5 Disagree strongly)

  • sblazy Social benefits/services make people lazy (1 Agree strongly - 5 Disagree strongly)

  • sblwcoa Social benefits/services make people less willing care for one another (1 Agree strongly - 5 Disagree strongly)

  • sblwlka Social benefits/services make people less willing look after themselves/family (1 Agree strongly - 5 Disagree strongly)

  • agea Respondent’s age

  • eduyrs Years of full-time education completed

  • gndr Gender (1 Male, 2 Female)

  • hinctnta Household’s total net income, all sources (Deciles of the actual household income range in the given country.)

  • gincdif Government should reduce differences in income levels (1 Agree strongly - 5 Disagree strongly)

  • dfincac Large differences in income acceptable to reward talents and efforts (1 Agree strongly - 5 Disagree strongly)

  • smdfslv For fair society, differences in standard of living should be small (1 Agree strongly - 5 Disagree strongly)

3 Environment preparation

First, let’s load the necessary packages to load, manipulate, visualize and analyse the data.

# Uncomment this once if you need to install the packages on your system 

### DATA MANIPULATION ###
# install.packages("haven")                 # data import from spss
# install.packages("dplyr")                 # data manipulation
# install.packages("psych")                 # descriptives
# install.packages("stringr")               # string manipulation
# install.packages("purrr")                 # table manipulation 

### MODELING ###
# install.packages("lavaan")                # SEM modelling
# install.packages("MVN")                   # tests for multivariate normality
# install.packages("Amelia")                # performing multiple imputation

### VISUALIZATION ###
# install.packages("tidySEM")               # plotting SEM models
# install.packages("ggplot2")               # plotting 
# install.packages("patchwork")             # organization plots

# Load the packages 

### DATA MANIPULATION ###
library("haven")        
library("dplyr")      
library("psych")
library("stringr")
library("purrr")

### MODELING ###
library("lavaan")       
library("MVN")
library("Amelia")

### VISUALIZATION ###
library("tidySEM")
library("ggplot2")              
library("patchwork")    

4 Ex.1 – Non-normal continuous data

  1. Check the multivariate normality for the question gincdif, dfincac, smdfslv (egalitarianism) and gvslvol gvhlthc gvcldcr gvpdlwk (welfare support) .
  2. Estimate a simple CFA (2 latent variable), without any error covariances for egalitarianism (gincdif, dfincac, smdfslv) and welfare support (gvslvol gvhlthc gvcldcr gvpdlwk)
  3. Estimate the model using the ML estimation, then re-estimate using MLM.
  4. Compare fit statistics, loadings, and standard errors.
  5. Remove the covariance between egalitarianism and welfare support and re-estimate the model. Compare the fit (OPTIONAL).
ess_df <- haven::read_sav("https://github.com/albertostefanelli/SEM_labs/raw/master/data/ESS4_belgium.sav")

# select the egalitarianism and welfare support items
ess_df_eg <- ess_df[,c("gincdif", "dfincac", "smdfslv", "gvslvol", "gvhlthc", "gvcldcr","gvpdlwk")]
# remove NAs
ess_df_eg_na <- na.omit(ess_df_eg)

mvn_test <- mvn(data = ess_df_eg_na, # our data without NAs
                mvnTest = c("hz")    # type of normality test to perform
                )

mvn_test
$multivariateNormality
           Test       HZ p value MVN
1 Henze-Zirkler 6.155718       0  NO

$univariateNormality
              Test  Variable Statistic   p value Normality
1 Anderson-Darling  gincdif   105.3500  <0.001      NO    
2 Anderson-Darling  dfincac   102.9968  <0.001      NO    
3 Anderson-Darling  smdfslv   124.3559  <0.001      NO    
4 Anderson-Darling  gvslvol    47.8262  <0.001      NO    
5 Anderson-Darling  gvhlthc    51.7952  <0.001      NO    
6 Anderson-Darling  gvcldcr    37.7151  <0.001      NO    
7 Anderson-Darling  gvpdlwk    39.3915  <0.001      NO    

$Descriptives
           n     Mean   Std.Dev Median Min Max 25th 75th       Skew   Kurtosis
gincdif 1728 2.230903 1.0612045      2   1   5    1    3  0.7362838 -0.2431670
dfincac 1728 2.626736 1.0550139      2   1   5    2    3  0.5063052 -0.5763850
smdfslv 1728 2.471065 0.9776167      2   1   5    2    3  0.6607562 -0.1964199
gvslvol 1728 7.884838 1.4786813      8   0  10    7    9 -0.7864489  1.3981895
gvhlthc 1728 8.042824 1.4717630      8   0  10    7    9 -0.9593240  1.8591224
gvcldcr 1728 7.287616 1.7607401      7   0  10    6    8 -0.8186374  1.3878985
gvpdlwk 1728 7.328704 1.7606153      8   0  10    6    8 -0.8356154  1.3712611
model_eg_ws <- '
## Egalitarianism ##
egual =~  gincdif + dfincac + smdfslv
## Welfare support ##
ws =~   gvslvol + gvhlthc + gvcldcr + gvpdlwk   
'

fit_eg_ws_ml <- cfa(model_eg_ws,    # model formula
                  data = ess_df,    # data frame
                  estimator = "ML"  # select the estimator 
)

summary(fit_eg_ws_ml,
        standardized=TRUE, 
        fit.measures = TRUE
)
lavaan 0.6.15 ended normally after 33 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                        15

                                                  Used       Total
  Number of observations                          1728        1760

Model Test User Model:
                                                      
  Test statistic                               291.138
  Degrees of freedom                                13
  P-value (Chi-square)                           0.000

Model Test Baseline Model:

  Test statistic                              2434.089
  Degrees of freedom                                21
  P-value                                        0.000

User Model versus Baseline Model:

  Comparative Fit Index (CFI)                    0.885
  Tucker-Lewis Index (TLI)                       0.814

Loglikelihood and Information Criteria:

  Loglikelihood user model (H0)             -19543.347
  Loglikelihood unrestricted model (H1)             NA
                                                      
  Akaike (AIC)                               39116.693
  Bayesian (BIC)                             39198.514
  Sample-size adjusted Bayesian (SABIC)      39150.861

Root Mean Square Error of Approximation:

  RMSEA                                          0.111
  90 Percent confidence interval - lower         0.100
  90 Percent confidence interval - upper         0.123
  P-value H_0: RMSEA <= 0.050                    0.000
  P-value H_0: RMSEA >= 0.080                    1.000

Standardized Root Mean Square Residual:

  SRMR                                           0.051

Parameter Estimates:

  Standard errors                             Standard
  Information                                 Expected
  Information saturated (h1) model          Structured

Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  egual =~                                                              
    gincdif           1.000                               0.688    0.648
    dfincac          -0.615    0.059  -10.362    0.000   -0.423   -0.401
    smdfslv           0.870    0.082   10.628    0.000    0.598    0.612
  ws =~                                                                 
    gvslvol           1.000                               1.161    0.785
    gvhlthc           0.912    0.038   23.711    0.000    1.059    0.720
    gvcldcr           0.847    0.043   19.896    0.000    0.983    0.559
    gvpdlwk           0.882    0.043   20.618    0.000    1.024    0.582

Covariances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  egual ~~                                                              
    ws               -0.209    0.030   -6.942    0.000   -0.262   -0.262

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .gincdif           0.653    0.049   13.424    0.000    0.653    0.580
   .dfincac           0.934    0.036   25.622    0.000    0.934    0.839
   .smdfslv           0.597    0.039   15.464    0.000    0.597    0.625
   .gvslvol           0.838    0.054   15.568    0.000    0.838    0.383
   .gvhlthc           1.043    0.053   19.752    0.000    1.043    0.482
   .gvcldcr           2.131    0.083   25.622    0.000    2.131    0.688
   .gvpdlwk           2.049    0.082   25.110    0.000    2.049    0.661
    egual             0.473    0.053    8.867    0.000    1.000    1.000
    ws                1.348    0.082   16.346    0.000    1.000    1.000
fit_eg_ws_mlr <- cfa(model_eg_ws,       # model formula
                  data = ess_df,        # data frame
                  estimator = "MLR"     # select the estimator 
)

summary(fit_eg_ws_mlr,
        standardized=TRUE, 
        fit.measures = TRUE
)
lavaan 0.6.15 ended normally after 33 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                        15

                                                  Used       Total
  Number of observations                          1728        1760

Model Test User Model:
                                              Standard      Scaled
  Test Statistic                               291.138     230.235
  Degrees of freedom                                13          13
  P-value (Chi-square)                           0.000       0.000
  Scaling correction factor                                  1.265
    Yuan-Bentler correction (Mplus variant)                       

Model Test Baseline Model:

  Test statistic                              2434.089    1588.615
  Degrees of freedom                                21          21
  P-value                                        0.000       0.000
  Scaling correction factor                                  1.532

User Model versus Baseline Model:

  Comparative Fit Index (CFI)                    0.885       0.861
  Tucker-Lewis Index (TLI)                       0.814       0.776
                                                                  
  Robust Comparative Fit Index (CFI)                         0.886
  Robust Tucker-Lewis Index (TLI)                            0.815

Loglikelihood and Information Criteria:

  Loglikelihood user model (H0)             -19543.347  -19543.347
  Scaling correction factor                                  1.684
      for the MLR correction                                      
  Loglikelihood unrestricted model (H1)             NA          NA
  Scaling correction factor                                  1.489
      for the MLR correction                                      
                                                                  
  Akaike (AIC)                               39116.693   39116.693
  Bayesian (BIC)                             39198.514   39198.514
  Sample-size adjusted Bayesian (SABIC)      39150.861   39150.861

Root Mean Square Error of Approximation:

  RMSEA                                          0.111       0.098
  90 Percent confidence interval - lower         0.100       0.089
  90 Percent confidence interval - upper         0.123       0.108
  P-value H_0: RMSEA <= 0.050                    0.000       0.000
  P-value H_0: RMSEA >= 0.080                    1.000       0.999
                                                                  
  Robust RMSEA                                               0.111
  90 Percent confidence interval - lower                     0.098
  90 Percent confidence interval - upper                     0.123
  P-value H_0: Robust RMSEA <= 0.050                         0.000
  P-value H_0: Robust RMSEA >= 0.080                         1.000

Standardized Root Mean Square Residual:

  SRMR                                           0.051       0.051

Parameter Estimates:

  Standard errors                             Sandwich
  Information bread                           Observed
  Observed information based on                Hessian

Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  egual =~                                                              
    gincdif           1.000                               0.688    0.648
    dfincac          -0.615    0.059  -10.340    0.000   -0.423   -0.401
    smdfslv           0.870    0.088    9.881    0.000    0.598    0.612
  ws =~                                                                 
    gvslvol           1.000                               1.161    0.785
    gvhlthc           0.912    0.034   26.668    0.000    1.059    0.720
    gvcldcr           0.847    0.081   10.453    0.000    0.983    0.559
    gvpdlwk           0.882    0.081   10.889    0.000    1.024    0.582

Covariances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  egual ~~                                                              
    ws               -0.209    0.030   -6.884    0.000   -0.262   -0.262

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .gincdif           0.653    0.055   11.828    0.000    0.653    0.580
   .dfincac           0.934    0.035   26.623    0.000    0.934    0.839
   .smdfslv           0.597    0.047   12.668    0.000    0.597    0.625
   .gvslvol           0.838    0.095    8.794    0.000    0.838    0.383
   .gvhlthc           1.043    0.091   11.499    0.000    1.043    0.482
   .gvcldcr           2.131    0.164   13.008    0.000    2.131    0.688
   .gvpdlwk           2.049    0.171   11.992    0.000    2.049    0.661
    egual             0.473    0.056    8.480    0.000    1.000    1.000
    ws                1.348    0.119   11.316    0.000    1.000    1.000
tidy_results_ml <- table_results(fit_eg_ws_ml,        
  columns = c("label", 
              "est_sig", 
              "se"),
  digits = 2
)

tidy_results_mlr <- table_results(fit_eg_ws_mlr,        
  columns = c("label", 
              "est_sig", 
              "se"),
  digits = 2
)



data.frame(Parameters = tidy_results_mlr$label, 
           "ML Model" = paste(tidy_results_ml$est_sig,"(",tidy_results_ml$se, ")"),
           "MLR Model" = paste(tidy_results_mlr$est_sig,"(",tidy_results_mlr$se, ")")
)
          Parameters          ML.Model         MLR.Model
1   egual.BY.gincdif     1.00 ( 0.00 )     1.00 ( 0.00 )
2   egual.BY.dfincac -0.61*** ( 0.06 ) -0.61*** ( 0.06 )
3   egual.BY.smdfslv  0.87*** ( 0.08 )  0.87*** ( 0.09 )
4      ws.BY.gvslvol     1.00 ( 0.00 )     1.00 ( 0.00 )
5      ws.BY.gvhlthc  0.91*** ( 0.04 )  0.91*** ( 0.03 )
6      ws.BY.gvcldcr  0.85*** ( 0.04 )  0.85*** ( 0.08 )
7      ws.BY.gvpdlwk  0.88*** ( 0.04 )  0.88*** ( 0.08 )
8  Variances.gincdif  0.65*** ( 0.05 )  0.65*** ( 0.06 )
9  Variances.dfincac  0.93*** ( 0.04 )  0.93*** ( 0.04 )
10 Variances.smdfslv  0.60*** ( 0.04 )  0.60*** ( 0.05 )
11 Variances.gvslvol  0.84*** ( 0.05 )  0.84*** ( 0.10 )
12 Variances.gvhlthc  1.04*** ( 0.05 )  1.04*** ( 0.09 )
13 Variances.gvcldcr  2.13*** ( 0.08 )  2.13*** ( 0.16 )
14 Variances.gvpdlwk  2.05*** ( 0.08 )  2.05*** ( 0.17 )
15   Variances.egual  0.47*** ( 0.05 )  0.47*** ( 0.06 )
16      Variances.ws  1.35*** ( 0.08 )  1.35*** ( 0.12 )
17     egual.WITH.ws -0.21*** ( 0.03 ) -0.21*** ( 0.03 )
model_eg_ws_cov <- '
## Egalitarianism ##
egual =~  gincdif + dfincac + smdfslv
## Welfare support ##
ws =~   gvslvol + gvhlthc + gvcldcr + gvpdlwk   
ws ~~ 0*egual
'

fit_eg_ws_mlr_cov <- cfa(model_eg_ws_cov,    # model formula
                  data = ess_df,             # data frame
                  estimator = "MLR"          # select the estimator 
)

summary(fit_eg_ws_mlr_cov,
        standardized=TRUE, 
        fit.measures = TRUE
)
lavaan 0.6.15 ended normally after 33 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                        14

                                                  Used       Total
  Number of observations                          1728        1760

Model Test User Model:
                                              Standard      Scaled
  Test Statistic                               347.870     279.437
  Degrees of freedom                                14          14
  P-value (Chi-square)                           0.000       0.000
  Scaling correction factor                                  1.245
    Yuan-Bentler correction (Mplus variant)                       

Model Test Baseline Model:

  Test statistic                              2434.089    1588.615
  Degrees of freedom                                21          21
  P-value                                        0.000       0.000
  Scaling correction factor                                  1.532

User Model versus Baseline Model:

  Comparative Fit Index (CFI)                    0.862       0.831
  Tucker-Lewis Index (TLI)                       0.792       0.746
                                                                  
  Robust Comparative Fit Index (CFI)                         0.862
  Robust Tucker-Lewis Index (TLI)                            0.794

Loglikelihood and Information Criteria:

  Loglikelihood user model (H0)             -19571.712  -19571.712
  Scaling correction factor                                  1.733
      for the MLR correction                                      
  Loglikelihood unrestricted model (H1)             NA          NA
  Scaling correction factor                                  1.489
      for the MLR correction                                      
                                                                  
  Akaike (AIC)                               39171.425   39171.425
  Bayesian (BIC)                             39247.791   39247.791
  Sample-size adjusted Bayesian (SABIC)      39203.314   39203.314

Root Mean Square Error of Approximation:

  RMSEA                                          0.117       0.105
  90 Percent confidence interval - lower         0.107       0.095
  90 Percent confidence interval - upper         0.128       0.114
  P-value H_0: RMSEA <= 0.050                    0.000       0.000
  P-value H_0: RMSEA >= 0.080                    1.000       1.000
                                                                  
  Robust RMSEA                                               0.117
  90 Percent confidence interval - lower                     0.105
  90 Percent confidence interval - upper                     0.129
  P-value H_0: Robust RMSEA <= 0.050                         0.000
  P-value H_0: Robust RMSEA >= 0.080                         1.000

Standardized Root Mean Square Residual:

  SRMR                                           0.082       0.082

Parameter Estimates:

  Standard errors                             Sandwich
  Information bread                           Observed
  Observed information based on                Hessian

Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  egual =~                                                              
    gincdif           1.000                               0.702    0.662
    dfincac          -0.612    0.064   -9.493    0.000   -0.429   -0.407
    smdfslv           0.827    0.089    9.275    0.000    0.581    0.595
  ws =~                                                                 
    gvslvol           1.000                               1.158    0.783
    gvhlthc           0.921    0.035   26.480    0.000    1.066    0.725
    gvcldcr           0.847    0.082   10.373    0.000    0.981    0.557
    gvpdlwk           0.881    0.081   10.851    0.000    1.020    0.580

Covariances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  egual ~~                                                              
    ws                0.000                               0.000    0.000

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .gincdif           0.632    0.059   10.726    0.000    0.632    0.562
   .dfincac           0.928    0.036   25.920    0.000    0.928    0.834
   .smdfslv           0.618    0.047   13.005    0.000    0.618    0.647
   .gvslvol           0.844    0.096    8.813    0.000    0.844    0.386
   .gvhlthc           1.028    0.092   11.212    0.000    1.028    0.475
   .gvcldcr           2.137    0.165   12.920    0.000    2.137    0.690
   .gvpdlwk           2.057    0.171   12.004    0.000    2.057    0.664
    egual             0.493    0.060    8.161    0.000    1.000    1.000
    ws                1.341    0.119   11.316    0.000    1.000    1.000
# let's compare the fit of the different models
model_fit <-  function(lavobject) {
  vars <- c("chisq", "df", "cfi", "tli", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper", "rmsea.pvalue", "srmr")
  return(fitmeasures(lavobject)[vars] %>% data.frame() %>% round(2) %>% t())
}

table_fit <- 
  list(model_fit(fit_eg_ws_mlr), 
       model_fit(fit_eg_ws_mlr_cov)) %>% 
  reduce(rbind)
 

rownames(table_fit) <- c("eg_ws_mlr", "eg_ws_mlr_cov")

table_fit
               chisq df  cfi  tli rmsea rmsea.ci.lower rmsea.ci.upper rmsea.pvalue srmr
eg_ws_mlr     291.14 13 0.88 0.81  0.11           0.10           0.12            0 0.05
eg_ws_mlr_cov 347.87 14 0.86 0.79  0.12           0.11           0.13            0 0.08
anova(fit_eg_ws_mlr, fit_eg_ws_mlr_cov)

Scaled Chi-Squared Difference Test (method = "satorra.bentler.2001")

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_eg_ws_mlr     13 39117 39199 291.14                                           
fit_eg_ws_mlr_cov 14 39171 39248 347.87     57.321       1 0.00000000000003702 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5 Ex.2 – Measurement invariance with categorical data

  1. Re-estimate a CFA (2 latent variable) for egalitarianism (gincdif, dfincac, smdfslv) and social criticism (bprvpv, sbeqsoc, sbcwkfm) treating the items as ordinal
  2. Test if the two latent constructs reach scalar invariance between male and feamle (free/fixed factor loading and thresholds).
  3. Compare and interpret the fit.
  4. Even if not necessary, relax certain equality constraints, re-fit the model, compare the fit.
  5. Compare the metric model fitted with continuous data with the one fitted with ordinal data (OPTIONAL)

Lavaan syntax for relaxing equality constraints:

y | c(label1,label1)*t1 # set threshold equal across 2 groups
y | c(label1,label2)*t1 # set threshold free across 2 groups
u3 ~*~ c(1,1)*u3 # fix scale of the endogenous variable "u3" to 1 in both groups
ess_df$gndr <- factor(ess_df$gndr,
                      levels = c("1", "2"),         # levels 
                      labels = c("Male", "Female")) # labels 

model_eg_sc <- '
## Egalitarianism ##
egual =~  gincdif + dfincac + smdfslv
## Social Criticism ##
wc_socia =~ sbprvpv + sbeqsoc + sbcwkfm
'

fit_eg_sc_conf<- cfa(model_eg_sc,           # model formula
                     data = ess_df,         # data frame
                     estimator = "WLSMV",   # select the estimator 
                     group = "gndr",
                     ordered = c(
                     "gincdif",
                     "dfincac",
                     "smdfslv",
                     "sbprvpv",
                     "sbeqsoc",
                     "sbcwkfm" )
          )


fit_eg_sc_metric <- cfa(model_eg_sc,           # model formula
                     data = ess_df,            # data frame
                     estimator = "WLSMV",      # select the estimator 
                     group = "gndr",
                     ordered = c(
                     "gincdif",
                     "dfincac",
                     "smdfslv",
                     "sbprvpv",
                     "sbeqsoc",
                     "sbcwkfm" ),
                     group.equal=c("loadings")
)

fit_eg_sc_scalar <- cfa(model_eg_sc,       # model formula
                     data = ess_df,        # data frame
                     estimator = "WLSMV",  # select the estimator 
                     group = "gndr",
                     ordered = c(
                     "gincdif",
                     "dfincac",
                     "smdfslv",
                     "sbprvpv",
                     "sbeqsoc",
                     "sbcwkfm" ),
                     group.equal=c("loadings",
                                "thresholds")
)

anova(fit_eg_sc_conf,fit_eg_sc_metric,fit_eg_sc_scalar)

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_eg_sc_conf   16          87.126                                 
fit_eg_sc_metric 20          91.971      5.770       4    0.21699   
fit_eg_sc_scalar 36         121.815     39.163      16    0.00103 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# According to the Lagrange multiplier test it seems that  an improvement in fit over the metric model is difficult 
# However, for didactic purpose we are going to free the the equality constraints with the highest chi-square contribution 
lavTestScore(fit_eg_sc_scalar)
$test

total score test:

   test     X2 df p.value
1 score 34.435 28   0.187

$uni

univariate score tests:

     lhs op   rhs     X2 df p.value
1   .p2. == .p55.  0.399  1   0.528
2   .p3. == .p56.  0.015  1   0.903
3   .p5. == .p58.  0.081  1   0.777
4   .p6. == .p59.  2.986  1   0.084
5   .p7. == .p60.  1.276  1   0.259
6   .p8. == .p61.  3.413  1   0.065
7   .p9. == .p62.  0.075  1   0.785
8  .p10. == .p63.  0.001  1   0.981
9  .p11. == .p64. 12.723  1   0.000
10 .p12. == .p65.  7.049  1   0.008
11 .p13. == .p66.  0.930  1   0.335
12 .p14. == .p67.  2.246  1   0.134
13 .p15. == .p68.  0.760  1   0.383
14 .p16. == .p69.  0.093  1   0.760
15 .p17. == .p70.  0.225  1   0.635
16 .p18. == .p71.  0.055  1   0.815
17 .p19. == .p72.  0.035  1   0.851
18 .p20. == .p73.  0.005  1   0.942
19 .p21. == .p74.  0.135  1   0.714
20 .p22. == .p75.  0.313  1   0.576
21 .p23. == .p76.  2.192  1   0.139
22 .p24. == .p77.  1.190  1   0.275
23 .p25. == .p78.  0.015  1   0.903
24 .p26. == .p79.  0.399  1   0.528
25 .p27. == .p80.  2.955  1   0.086
26 .p28. == .p81.  5.497  1   0.019
27 .p29. == .p82.  1.005  1   0.316
28 .p30. == .p83.  0.111  1   0.739
parTable(fit_eg_sc_scalar)
     id      lhs  op      rhs user block group free ustart exo label plabel  start    est    se
1     1    egual  =~  gincdif    1     1     1    0      1   0         .p1.  1.000  1.000 0.000
2     2    egual  =~  dfincac    1     1     1    1     NA   0  .p2.   .p2. -0.648 -0.617 0.048
3     3    egual  =~  smdfslv    1     1     1    2     NA   0  .p3.   .p3.  0.820  0.986 0.076
4     4 wc_socia  =~  sbprvpv    1     1     1    0      1   0         .p4.  1.000  1.000 0.000
5     5 wc_socia  =~  sbeqsoc    1     1     1    3     NA   0  .p5.   .p5.  1.523  1.247 0.092
6     6 wc_socia  =~  sbcwkfm    1     1     1    4     NA   0  .p6.   .p6.  0.753  0.692 0.045
7     7  gincdif   |       t1    0     1     1    5     NA   0  .p7.   .p7. -0.673 -0.701 0.042
8     8  gincdif   |       t2    0     1     1    6     NA   0  .p8.   .p8.  0.513  0.459 0.040
9     9  gincdif   |       t3    0     1     1    7     NA   0  .p9.   .p9.  0.922  0.915 0.046
10   10  gincdif   |       t4    0     1     1    8     NA   0 .p10.  .p10.  1.885  1.886 0.077
11   11  dfincac   |       t1    0     1     1    9     NA   0 .p11.  .p11. -1.118 -1.217 0.048
12   12  dfincac   |       t2    0     1     1   10     NA   0 .p12.  .p12.  0.259  0.178 0.034
13   13  dfincac   |       t3    0     1     1   11     NA   0 .p13.  .p13.  0.765  0.735 0.039
14   14  dfincac   |       t4    0     1     1   12     NA   0 .p14.  .p14.  1.760  1.691 0.065
15   15  smdfslv   |       t1    0     1     1   13     NA   0 .p15.  .p15. -1.212 -1.237 0.052
16   16  smdfslv   |       t2    0     1     1   14     NA   0 .p16.  .p16.  0.268  0.259 0.038
17   17  smdfslv   |       t3    0     1     1   15     NA   0 .p17.  .p17.  0.852  0.838 0.044
18   18  smdfslv   |       t4    0     1     1   16     NA   0 .p18.  .p18.  1.835  1.824 0.073
19   19  sbprvpv   |       t1    0     1     1   17     NA   0 .p19.  .p19. -1.175 -1.170 0.051
20   20  sbprvpv   |       t2    0     1     1   18     NA   0 .p20.  .p20.  0.537  0.539 0.038
21   21  sbprvpv   |       t3    0     1     1   19     NA   0 .p21.  .p21.  1.118  1.106 0.047
22   22  sbprvpv   |       t4    0     1     1   20     NA   0 .p22.  .p22.  2.224  2.183 0.092
23   23  sbeqsoc   |       t1    0     1     1   21     NA   0 .p23.  .p23. -1.346 -1.386 0.060
24   24  sbeqsoc   |       t2    0     1     1   22     NA   0 .p24.  .p24.  0.540  0.509 0.041
25   25  sbeqsoc   |       t3    0     1     1   23     NA   0 .p25.  .p25.  1.169  1.165 0.050
26   26  sbeqsoc   |       t4    0     1     1   24     NA   0 .p26.  .p26.  2.224  2.177 0.091
27   27  sbcwkfm   |       t1    0     1     1   25     NA   0 .p27.  .p27. -1.493 -1.428 0.057
28   28  sbcwkfm   |       t2    0     1     1   26     NA   0 .p28.  .p28.  0.472  0.542 0.036
29   29  sbcwkfm   |       t3    0     1     1   27     NA   0 .p29.  .p29.  1.230  1.264 0.048
30   30  sbcwkfm   |       t4    0     1     1   28     NA   0 .p30.  .p30.  2.301  2.327 0.100
31   31  gincdif  ~~  gincdif    0     1     1    0      1   0        .p31.  1.000  0.531 0.000
32   32  dfincac  ~~  dfincac    0     1     1    0      1   0        .p32.  1.000  0.822 0.000
33   33  smdfslv  ~~  smdfslv    0     1     1    0      1   0        .p33.  1.000  0.544 0.000
34   34  sbprvpv  ~~  sbprvpv    0     1     1    0      1   0        .p34.  1.000  0.587 0.000
35   35  sbeqsoc  ~~  sbeqsoc    0     1     1    0      1   0        .p35.  1.000  0.358 0.000
36   36  sbcwkfm  ~~  sbcwkfm    0     1     1    0      1   0        .p36.  1.000  0.802 0.000
37   37    egual  ~~    egual    0     1     1   29     NA   0        .p37.  0.050  0.469 0.044
38   38 wc_socia  ~~ wc_socia    0     1     1   30     NA   0        .p38.  0.050  0.413 0.039
39   39    egual  ~~ wc_socia    0     1     1   31     NA   0        .p39.  0.000  0.142 0.023
40   40  gincdif ~*~  gincdif    0     1     1    0      1   0        .p40.  1.000  1.000 0.000
41   41  dfincac ~*~  dfincac    0     1     1    0      1   0        .p41.  1.000  1.000 0.000
42   42  smdfslv ~*~  smdfslv    0     1     1    0      1   0        .p42.  1.000  1.000 0.000
43   43  sbprvpv ~*~  sbprvpv    0     1     1    0      1   0        .p43.  1.000  1.000 0.000
44   44  sbeqsoc ~*~  sbeqsoc    0     1     1    0      1   0        .p44.  1.000  1.000 0.000
45   45  sbcwkfm ~*~  sbcwkfm    0     1     1    0      1   0        .p45.  1.000  1.000 0.000
46   46  gincdif  ~1             0     1     1    0      0   0        .p46.  0.000  0.000 0.000
47   47  dfincac  ~1             0     1     1    0      0   0        .p47.  0.000  0.000 0.000
48   48  smdfslv  ~1             0     1     1    0      0   0        .p48.  0.000  0.000 0.000
49   49  sbprvpv  ~1             0     1     1    0      0   0        .p49.  0.000  0.000 0.000
50   50  sbeqsoc  ~1             0     1     1    0      0   0        .p50.  0.000  0.000 0.000
51   51  sbcwkfm  ~1             0     1     1    0      0   0        .p51.  0.000  0.000 0.000
52   52    egual  ~1             0     1     1    0      0   0        .p52.  0.000  0.000 0.000
53   53 wc_socia  ~1             0     1     1    0      0   0        .p53.  0.000  0.000 0.000
54   54    egual  =~  gincdif    1     2     2    0      1   0        .p54.  1.000  1.000 0.000
55   55    egual  =~  dfincac    1     2     2   32     NA   0  .p2.  .p55. -0.610 -0.617 0.048
56   56    egual  =~  smdfslv    1     2     2   33     NA   0  .p3.  .p56.  1.042  0.986 0.076
57   57 wc_socia  =~  sbprvpv    1     2     2    0      1   0        .p57.  1.000  1.000 0.000
58   58 wc_socia  =~  sbeqsoc    1     2     2   34     NA   0  .p5.  .p58.  1.109  1.247 0.092
59   59 wc_socia  =~  sbcwkfm    1     2     2   35     NA   0  .p6.  .p59.  0.610  0.692 0.045
60   60  gincdif   |       t1    0     2     2   36     NA   0  .p7.  .p60. -0.636 -0.701 0.042
61   61  gincdif   |       t2    0     2     2   37     NA   0  .p8.  .p61.  0.516  0.459 0.040
62   62  gincdif   |       t3    0     2     2   38     NA   0  .p9.  .p62.  1.020  0.915 0.046
63   63  gincdif   |       t4    0     2     2   39     NA   0 .p10.  .p63.  2.015  1.886 0.077
64   64  dfincac   |       t1    0     2     2   40     NA   0 .p11.  .p64. -1.422 -1.217 0.048
65   65  dfincac   |       t2    0     2     2   41     NA   0 .p12.  .p65.  0.042  0.178 0.034
66   66  dfincac   |       t3    0     2     2   42     NA   0 .p13.  .p66.  0.654  0.735 0.039
67   67  dfincac   |       t4    0     2     2   43     NA   0 .p14.  .p67.  1.594  1.691 0.065
68   68  smdfslv   |       t1    0     2     2   44     NA   0 .p15.  .p68. -1.209 -1.237 0.052
69   69  smdfslv   |       t2    0     2     2   45     NA   0 .p16.  .p69.  0.366  0.259 0.038
70   70  smdfslv   |       t3    0     2     2   46     NA   0 .p17.  .p70.  0.964  0.838 0.044
71   71  smdfslv   |       t4    0     2     2   47     NA   0 .p18.  .p71.  1.993  1.824 0.073
72   72  sbprvpv   |       t1    0     2     2   48     NA   0 .p19.  .p72. -1.259 -1.170 0.051
73   73  sbprvpv   |       t2    0     2     2   49     NA   0 .p20.  .p73.  0.523  0.539 0.038
74   74  sbprvpv   |       t3    0     2     2   50     NA   0 .p21.  .p74.  1.102  1.106 0.047
75   75  sbprvpv   |       t4    0     2     2   51     NA   0 .p22.  .p75.  2.201  2.183 0.092
76   76  sbeqsoc   |       t1    0     2     2   52     NA   0 .p23.  .p76. -1.498 -1.386 0.060
77   77  sbeqsoc   |       t2    0     2     2   53     NA   0 .p24.  .p77.  0.435  0.509 0.041
78   78  sbeqsoc   |       t3    0     2     2   54     NA   0 .p25.  .p78.  1.123  1.165 0.050
79   79  sbeqsoc   |       t4    0     2     2   55     NA   0 .p26.  .p79.  2.112  2.177 0.091
80   80  sbcwkfm   |       t1    0     2     2   56     NA   0 .p27.  .p80. -1.447 -1.428 0.057
81   81  sbcwkfm   |       t2    0     2     2   57     NA   0 .p28.  .p81.  0.601  0.542 0.036
82   82  sbcwkfm   |       t3    0     2     2   58     NA   0 .p29.  .p82.  1.312  1.264 0.048
83   83  sbcwkfm   |       t4    0     2     2   59     NA   0 .p30.  .p83.  2.405  2.327 0.100
84   84  gincdif  ~~  gincdif    0     2     2    0      1   0        .p84.  1.000  0.525 0.000
85   85  dfincac  ~~  dfincac    0     2     2    0      1   0        .p85.  1.000  0.807 0.000
86   86  smdfslv  ~~  smdfslv    0     2     2    0      1   0        .p86.  1.000  0.483 0.000
87   87  sbprvpv  ~~  sbprvpv    0     2     2    0      1   0        .p87.  1.000  0.506 0.000
88   88  sbeqsoc  ~~  sbeqsoc    0     2     2    0      1   0        .p88.  1.000  0.341 0.000
89   89  sbcwkfm  ~~  sbcwkfm    0     2     2    0      1   0        .p89.  1.000  0.740 0.000
90   90    egual  ~~    egual    0     2     2   60     NA   0        .p90.  0.050  0.449 0.056
91   91 wc_socia  ~~ wc_socia    0     2     2   61     NA   0        .p91.  0.050  0.410 0.047
92   92    egual  ~~ wc_socia    0     2     2   62     NA   0        .p92.  0.000  0.057 0.021
93   93  gincdif ~*~  gincdif    0     2     2   63     NA   0        .p93.  1.000  1.013 0.044
94   94  dfincac ~*~  dfincac    0     2     2   64     NA   0        .p94.  1.000  1.011 0.040
95   95  smdfslv ~*~  smdfslv    0     2     2   65     NA   0        .p95.  1.000  1.043 0.042
96   96  sbprvpv ~*~  sbprvpv    0     2     2   66     NA   0        .p96.  1.000  1.045 0.041
97   97  sbeqsoc ~*~  sbeqsoc    0     2     2   67     NA   0        .p97.  1.000  1.011 0.041
98   98  sbcwkfm ~*~  sbcwkfm    0     2     2   68     NA   0        .p98.  1.000  1.033 0.039
99   99  gincdif  ~1             0     2     2    0      0   0        .p99.  0.000  0.000 0.000
100 100  dfincac  ~1             0     2     2    0      0   0       .p100.  0.000  0.000 0.000
101 101  smdfslv  ~1             0     2     2    0      0   0       .p101.  0.000  0.000 0.000
102 102  sbprvpv  ~1             0     2     2    0      0   0       .p102.  0.000  0.000 0.000
103 103  sbeqsoc  ~1             0     2     2    0      0   0       .p103.  0.000  0.000 0.000
104 104  sbcwkfm  ~1             0     2     2    0      0   0       .p104.  0.000  0.000 0.000
105 105    egual  ~1             0     2     2   69     NA   0       .p105.  0.000 -0.101 0.043
106 106 wc_socia  ~1             0     2     2   70     NA   0       .p106.  0.000  0.040 0.040
107 107     .p2.  ==    .p55.    2     0     0    0     NA   0               0.000  0.000 0.000
108 108     .p3.  ==    .p56.    2     0     0    0     NA   0               0.000  0.000 0.000
109 109     .p5.  ==    .p58.    2     0     0    0     NA   0               0.000  0.000 0.000
110 110     .p6.  ==    .p59.    2     0     0    0     NA   0               0.000  0.000 0.000
111 111     .p7.  ==    .p60.    2     0     0    0     NA   0               0.000  0.000 0.000
112 112     .p8.  ==    .p61.    2     0     0    0     NA   0               0.000  0.000 0.000
113 113     .p9.  ==    .p62.    2     0     0    0     NA   0               0.000  0.000 0.000
114 114    .p10.  ==    .p63.    2     0     0    0     NA   0               0.000  0.000 0.000
115 115    .p11.  ==    .p64.    2     0     0    0     NA   0               0.000  0.000 0.000
116 116    .p12.  ==    .p65.    2     0     0    0     NA   0               0.000  0.000 0.000
117 117    .p13.  ==    .p66.    2     0     0    0     NA   0               0.000  0.000 0.000
118 118    .p14.  ==    .p67.    2     0     0    0     NA   0               0.000  0.000 0.000
119 119    .p15.  ==    .p68.    2     0     0    0     NA   0               0.000  0.000 0.000
120 120    .p16.  ==    .p69.    2     0     0    0     NA   0               0.000  0.000 0.000
121 121    .p17.  ==    .p70.    2     0     0    0     NA   0               0.000  0.000 0.000
122 122    .p18.  ==    .p71.    2     0     0    0     NA   0               0.000  0.000 0.000
123 123    .p19.  ==    .p72.    2     0     0    0     NA   0               0.000  0.000 0.000
124 124    .p20.  ==    .p73.    2     0     0    0     NA   0               0.000  0.000 0.000
125 125    .p21.  ==    .p74.    2     0     0    0     NA   0               0.000  0.000 0.000
126 126    .p22.  ==    .p75.    2     0     0    0     NA   0               0.000  0.000 0.000
127 127    .p23.  ==    .p76.    2     0     0    0     NA   0               0.000  0.000 0.000
128 128    .p24.  ==    .p77.    2     0     0    0     NA   0               0.000  0.000 0.000
129 129    .p25.  ==    .p78.    2     0     0    0     NA   0               0.000  0.000 0.000
130 130    .p26.  ==    .p79.    2     0     0    0     NA   0               0.000  0.000 0.000
131 131    .p27.  ==    .p80.    2     0     0    0     NA   0               0.000  0.000 0.000
132 132    .p28.  ==    .p81.    2     0     0    0     NA   0               0.000  0.000 0.000
133 133    .p29.  ==    .p82.    2     0     0    0     NA   0               0.000  0.000 0.000
134 134    .p30.  ==    .p83.    2     0     0    0     NA   0               0.000  0.000 0.000
# according to lavTestScore the most problematic parameter is .p11. ==  .p64. which correspond to dfincac   |   t1

model_eg_sc_p <- '
## Egalitarianism ##
egual =~  gincdif + dfincac + smdfslv
## Social Criticism ##
wc_socia =~ sbprvpv + sbeqsoc + sbcwkfm
# fix scale of variable "dfincac" to 1 in both groups
dfincac ~*~ c(1,1)*dfincac

'

fit_eg_sc_scalar_p <- cfa(model_eg_sc_p,       # model formula
                     data = ess_df,            # data frame
                     estimator = "WLSMV",      # select the estimator 
                     group = "gndr",
                     ordered = c(
                     "gincdif",
                     "dfincac",
                     "smdfslv",
                     "sbprvpv",
                     "sbeqsoc",
                     "sbcwkfm" ),
                     group.equal=c("loadings","thresholds"),
                     # relax the threshold t1 of the variable dfincac
                     group.partial=c("equal =~ dfincac","dfincac | t1")) 

# we do not reach partial invariance
anova(fit_eg_sc_conf,fit_eg_sc_metric,fit_eg_sc_scalar_p)

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_eg_sc_conf     16          87.126                                
fit_eg_sc_metric   20          91.971     5.7701       4    0.21699  
fit_eg_sc_scalar_p 36         113.037    27.8118      16    0.03329 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit_eg_sc_metric_cont <- cfa(model_eg_sc,           # model formula
                     data = ess_df,            # data frame
                     estimator = "ML",         # select the estimator 
                     group = "gndr",
                     group.equal=c("loadings")
)

table_fit <- 
  list(model_fit(fit_eg_sc_metric_cont), 
       model_fit(fit_eg_sc_metric)) %>% 
  reduce(rbind)
 

rownames(table_fit) <- c("Continous", "Ordinal")

table_fit
          chisq df  cfi  tli rmsea rmsea.ci.lower rmsea.ci.upper rmsea.pvalue srmr
Continous 64.01 20 0.96 0.94  0.05           0.04           0.06         0.44 0.04
Ordinal   91.97 20 0.97 0.96  0.06           0.05           0.08         0.03 0.05
tidy_results_continous <- table_results(fit_eg_sc_metric_cont,        
  columns = c("label", 
              "est_sig", 
              "se"),
  digits = 2
)

tidy_results_ordinal <- table_results(fit_eg_sc_metric,        
  columns = c("label", 
              "est_sig", 
              "se"),
  digits = 2
)



data.frame(Parameters = tidy_results_continous$label[1:6], 
           "Continous Model" = paste(tidy_results_continous$est_sig[1:6],"(",tidy_results_continous$se[1:6], ")"),
           "Ordinal Model" = paste(tidy_results_ordinal$est_sig[1:6],"(",tidy_results_ordinal$se[1:6], ")")
)
                Parameters   Continous.Model     Ordinal.Model
1    egual.BY.gincdif.Male     1.00 ( 0.00 )     1.00 ( 0.00 )
2    egual.BY.dfincac.Male -0.59*** ( 0.06 ) -0.60*** ( 0.05 )
3    egual.BY.smdfslv.Male  0.84*** ( 0.08 )  0.99*** ( 0.08 )
4 wc_socia.BY.sbprvpv.Male     1.00 ( 0.00 )     1.00 ( 0.00 )
5 wc_socia.BY.sbeqsoc.Male  1.24*** ( 0.12 )  1.22*** ( 0.09 )
6 wc_socia.BY.sbcwkfm.Male  0.62*** ( 0.05 )  0.69*** ( 0.04 )

6 !!Support Ukraine!!