2.1 Software used: R

- This is not an R course!
- We will learn some R as we go along
- I will use RStudio
- Many packages or libraries exist to do specific analyses

7.1 1-factor model: 3 indicators

Now that we have an idea of what is in our data we apply our theoretical knowledge to the our empirical data. Our measurement model will be formed by 3 items that loads on the latent factor “Welfare Support.”

7.2 Degrees of freedom

Sorts of identification:

1. non-identified: \(df < 0\), impossible to estimate parameters
2. over-identified: \(df > 0\), we should strive for this
3. just-identified: \(df = 0\), it is ok, but the fit cannot be assessed
4. empirical under-identification: it can happen when two or more indicators are highly correlated (multicollinearity)

We can quickly check if our model has positive degrees of freedom using the following formula.

Degrees of freedom = Number of parameters to estimate - Pieces of information

Pieces of information = \(\dfrac{p(p+1)}{2}=\dfrac{3(3+1)}{2}=6\).

where \(p\) is the number of “manifest” indicators \(i\) (also called “measured variables” or “observed indicators”), in this case gvslvol gvslvue, gvhlthc.

Number of parameters to estimate = \(\psi{f} + \lambda_{ji} + \theta_i= 1 + 2 + 3=6\)

where \(\psi{f}\) is the latent factor \(f\) variance, \(\lambda_{fi}\) is the loading on the latent factor \(f\) for each observed indicator \(i\), and \(\theta_{i}\) is the residual variance for each observed indicator \(i\).

Q: Why \(\lambda_{fi}=2\) ? (Pro Tip: (do not count parameters we fixed for scaling).

Degrees of freedom = Number of parameters to estimate - Pieces of informations =\(6-6=0\)

Q: What are the implication of having 0 degrees of freedom?

Remember:

• Factor variance (\(\psi{j}\)): how much individuals differ on the factor
• Factor loadings (\(\lambda_{ji}\)): the effect of the factor or latent variable on the measure
• Residual (error) variance (\(\theta_{i}\)) : variance in the measure not attributable to the factor

Let’s now fit our CFA model using gvslvol, gvslvue, gvhlthc.

model_ws_3 <-'welf_supp =~ gvslvol + gvslvue + gvhlthc'

fit_ws_3 <- cfa(model_ws_3,             # model formula
                data = ess_df_selected  # data frame
  )

summary(fit_ws_3, 
        standardized = TRUE)

lavaan 0.6-10 ended normally after 26 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                         6
                                                      
                                                  Used       Total
  Number of observations                          1752        1760
                                                                  
Model Test User Model:
                                                      
  Test statistic                                 0.000
  Degrees of freedom                                 0

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
  welf_supp =~                                                          
    gvslvol           1.000                               1.233    0.832
    gvslvue           0.593    0.048   12.395    0.000    0.730    0.381
    gvhlthc           0.890    0.062   14.292    0.000    1.097    0.746

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .gvslvol           0.677    0.103    6.555    0.000    0.677    0.308
   .gvslvue           3.143    0.112   28.083    0.000    3.143    0.855
   .gvhlthc           0.962    0.086   11.171    0.000    0.962    0.444
    welf_supp         1.519    0.123   12.359    0.000    1.000    1.000

7.3 Lavaan syntax

Let’s take a closer look at the lavaan syntax

• latent variable definition =~ (is measured by)
• regression ~ (is regressed on)
• (residual) (co)variance ~~ (is correlated with)
• intercept ~1 (intercept)
• new model parameter := (is equal to)

So, if we want, we can freely estimate the first factor loading and fix the factor variance to 1:

7.4 Alternative model specification

Remember: Every factor in CFA should be given a metric to be identified.

1. Fix one factor loading to 1 (“marker variable” method)
2. Fix factor variance to 1
3. Fix the average loading to 1 (“effect coding” method).

model_ws_3_alt <-'
welf_supp =~ NA*gvslvol + gvslvue + gvhlthc
welf_supp~~1*welf_supp
'

fit_ws_3_alt <- cfa(model_ws_3_alt,            # model formula
                    data = ess_df_selected     # data frame
  )

summary(fit_ws_3_alt, 
        standardized = TRUE)

lavaan 0.6-10 ended normally after 21 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                         6
                                                      
                                                  Used       Total
  Number of observations                          1752        1760
                                                                  
Model Test User Model:
                                                      
  Test statistic                                 0.000
  Degrees of freedom                                 0

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
  welf_supp =~                                                          
    gvslvol           1.233    0.050   24.718    0.000    1.233    0.832
    gvslvue           0.730    0.050   14.516    0.000    0.730    0.381
    gvhlthc           1.097    0.047   23.321    0.000    1.097    0.746

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
    welf_supp         1.000                               1.000    1.000
   .gvslvol           0.677    0.103    6.555    0.000    0.677    0.308
   .gvslvue           3.143    0.112   28.083    0.000    3.143    0.855
   .gvhlthc           0.962    0.086   11.171    0.000    0.962    0.444

Q: What is the difference between these two models?

7.5 1-factor model: 6 indicators

Since we have 3 more indicators that measure our latent construct welfare support, let’s fit a model with 6 items instead of 3. This will allows us to estimate fit indices, one of the most valuable metric in a SEM model.

model_ws_6 <-'welf_supp =~ gvslvol + gvslvue + gvhlthc + gvcldcr + gvjbevn + gvpdlwk'


fit_ws_6 <- cfa(model_ws_6,             # model formula
                data = ess_df_selected  # data frame
  )

summary(fit_ws_6, 
        standardized = TRUE)

lavaan 0.6-10 ended normally after 28 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                        12
                                                      
                                                  Used       Total
  Number of observations                          1740        1760
                                                                  
Model Test User Model:
                                                      
  Test statistic                               326.956
  Degrees of freedom                                 9
  P-value (Chi-square)                           0.000

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
  welf_supp =~                                                          
    gvslvol           1.000                               1.129    0.761
    gvslvue           0.727    0.046   15.908    0.000    0.821    0.428
    gvhlthc           0.945    0.037   25.657    0.000    1.067    0.724
    gvcldcr           0.895    0.042   21.086    0.000    1.010    0.573
    gvjbevn           1.042    0.054   19.184    0.000    1.176    0.519
    gvpdlwk           0.913    0.042   21.526    0.000    1.031    0.586

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .gvslvol           0.923    0.048   19.088    0.000    0.923    0.420
   .gvslvue           3.001    0.108   27.859    0.000    3.001    0.817
   .gvhlthc           1.031    0.049   21.109    0.000    1.031    0.476
   .gvcldcr           2.087    0.081   25.883    0.000    2.087    0.672
   .gvjbevn           3.755    0.140   26.790    0.000    3.755    0.731
   .gvpdlwk           2.032    0.079   25.627    0.000    2.032    0.657
    welf_supp         1.274    0.077   16.546    0.000    1.000    1.000

Q: Why do we have 9 degrees of freedom in this model?

7.6 Model diagnostics

Important things to check after running the model:

1. Estimates
   • Do the indicators load well on the factor(s)?
2. Model fit [see next]
   • Are the fit indices good?
3. Heywood cases/ reasonable solution
   • Are variances positive?
   • Are \(r^2\) is below 1?
   • Are there any extreme standard errors?

# Returns the observed covariance matrix.
lavInspect(fit_ws_6, "sampstat")

$cov
        gvslvl gvslvu gvhlth gvcldc gvjbvn gvpdlw
gvslvol 2.197                                    
gvslvue 0.894  3.675                             
gvhlthc 1.353  0.800  2.169                      
gvcldcr 1.034  0.785  0.886  3.107               
gvjbevn 1.197  1.476  1.285  1.226  5.139        
gvpdlwk 1.060  0.806  0.946  1.689  1.156  3.095 

# Returns the model estimated covariance matrix.
fitted(fit_ws_6)

$cov
        gvslvl gvslvu gvhlth gvcldc gvjbvn gvpdlw
gvslvol 2.197                                    
gvslvue 0.926  3.675                             
gvhlthc 1.204  0.875  2.169                      
gvcldcr 1.140  0.829  1.077  3.107               
gvjbevn 1.328  0.966  1.255  1.188  5.139        
gvpdlwk 1.163  0.846  1.099  1.041  1.212  3.095 

# Returns the difference between the observed and estimated covariance matrix
residuals(fit_ws_6)

$type
[1] "raw"

$cov
        gvslvl gvslvu gvhlth gvcldc gvjbvn gvpdlw
gvslvol  0.000                                   
gvslvue -0.033  0.000                            
gvhlthc  0.149 -0.075  0.000                     
gvcldcr -0.106 -0.044 -0.192  0.000              
gvjbevn -0.131  0.510  0.030  0.037  0.000       
gvpdlwk -0.103 -0.040 -0.154  0.648 -0.056  0.000

# Returns the standard errors for the free parameters in the model.
lavTech(fit_ws_6, "se")

$lambda
           [,1]
[1,] 0.00000000
[2,] 0.04571461
[3,] 0.03682999
[4,] 0.04244399
[5,] 0.05432472
[6,] 0.04241931

$theta
          [,1]      [,2]       [,3]       [,4]      [,5]       [,6]
[1,] 0.0483618 0.0000000 0.00000000 0.00000000 0.0000000 0.00000000
[2,] 0.0000000 0.1077257 0.00000000 0.00000000 0.0000000 0.00000000
[3,] 0.0000000 0.0000000 0.04886228 0.00000000 0.0000000 0.00000000
[4,] 0.0000000 0.0000000 0.00000000 0.08062486 0.0000000 0.00000000
[5,] 0.0000000 0.0000000 0.00000000 0.00000000 0.1401613 0.00000000
[6,] 0.0000000 0.0000000 0.00000000 0.00000000 0.0000000 0.07930721

$psi
           [,1]
[1,] 0.07699622

# Returns the parameter estimates with confidence intervals.
parameterEstimates(fit_ws_6, standardized=TRUE)

         lhs op       rhs   est    se      z pvalue ci.lower ci.upper std.lv std.all std.nox
1  welf_supp =~   gvslvol 1.000 0.000     NA     NA    1.000    1.000  1.129   0.761   0.761
2  welf_supp =~   gvslvue 0.727 0.046 15.908      0    0.638    0.817  0.821   0.428   0.428
3  welf_supp =~   gvhlthc 0.945 0.037 25.657      0    0.873    1.017  1.067   0.724   0.724
4  welf_supp =~   gvcldcr 0.895 0.042 21.086      0    0.812    0.978  1.010   0.573   0.573
5  welf_supp =~   gvjbevn 1.042 0.054 19.184      0    0.936    1.149  1.176   0.519   0.519
6  welf_supp =~   gvpdlwk 0.913 0.042 21.526      0    0.830    0.996  1.031   0.586   0.586
7    gvslvol ~~   gvslvol 0.923 0.048 19.088      0    0.828    1.018  0.923   0.420   0.420
8    gvslvue ~~   gvslvue 3.001 0.108 27.859      0    2.790    3.212  3.001   0.817   0.817
9    gvhlthc ~~   gvhlthc 1.031 0.049 21.109      0    0.936    1.127  1.031   0.476   0.476
10   gvcldcr ~~   gvcldcr 2.087 0.081 25.883      0    1.929    2.245  2.087   0.672   0.672
11   gvjbevn ~~   gvjbevn 3.755 0.140 26.790      0    3.480    4.030  3.755   0.731   0.731
12   gvpdlwk ~~   gvpdlwk 2.032 0.079 25.627      0    1.877    2.188  2.032   0.657   0.657
13 welf_supp ~~ welf_supp 1.274 0.077 16.546      0    1.123    1.425  1.000   1.000   1.000

Lavaan offers a quick way to check model fit and statistics

summary(fit_ws_6,            # fitted model 
        fit.measures = TRUE, # returns commonly used fit measures 
        standardized = TRUE  # indicates that we want standardized results
        )

lavaan 0.6-10 ended normally after 28 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                        12
                                                      
                                                  Used       Total
  Number of observations                          1740        1760
                                                                  
Model Test User Model:
                                                      
  Test statistic                               326.956
  Degrees of freedom                                 9
  P-value (Chi-square)                           0.000

Model Test Baseline Model:

  Test statistic                              2621.015
  Degrees of freedom                                15
  P-value                                        0.000

User Model versus Baseline Model:

  Comparative Fit Index (CFI)                    0.878
  Tucker-Lewis Index (TLI)                       0.797

Loglikelihood and Information Criteria:

  Loglikelihood user model (H0)             -19550.568
  Loglikelihood unrestricted model (H1)             NA
                                                      
  Akaike (AIC)                               39125.136
  Bayesian (BIC)                             39190.676
  Sample-size adjusted Bayesian (BIC)        39152.553

Root Mean Square Error of Approximation:

  RMSEA                                          0.142
  90 Percent confidence interval - lower         0.129
  90 Percent confidence interval - upper         0.156
  P-value RMSEA <= 0.05                          0.000

Standardized Root Mean Square Residual:

  SRMR                                           0.061

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
  welf_supp =~                                                          
    gvslvol           1.000                               1.129    0.761
    gvslvue           0.727    0.046   15.908    0.000    0.821    0.428
    gvhlthc           0.945    0.037   25.657    0.000    1.067    0.724
    gvcldcr           0.895    0.042   21.086    0.000    1.010    0.573
    gvjbevn           1.042    0.054   19.184    0.000    1.176    0.519
    gvpdlwk           0.913    0.042   21.526    0.000    1.031    0.586

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .gvslvol           0.923    0.048   19.088    0.000    0.923    0.420
   .gvslvue           3.001    0.108   27.859    0.000    3.001    0.817
   .gvhlthc           1.031    0.049   21.109    0.000    1.031    0.476
   .gvcldcr           2.087    0.081   25.883    0.000    2.087    0.672
   .gvjbevn           3.755    0.140   26.790    0.000    3.755    0.731
   .gvpdlwk           2.032    0.079   25.627    0.000    2.032    0.657
    welf_supp         1.274    0.077   16.546    0.000    1.000    1.000

# we can also just extract the factor loadings.
# using lavaan terminology, the factors loadings are the lambdas

inspect(fit_ws_6,           # fitted model 
        what="std")$lambda  # standardized factors loadings

        wlf_sp
gvslvol  0.761
gvslvue  0.428
gvhlthc  0.724
gvcldcr  0.573
gvjbevn  0.519
gvpdlwk  0.586

# more info on the factors loading can be obtained using the tidy SEM package
# first we extract all the estimated parameters 

tidy_results <- table_results(fit_ws_6,             
  columns = c("label", 
              "est_sig", 
              "se", 
              "confint"),
  digits = 2
)

tidy_results %>% filter(str_detect(label, "welf_supp."))

                 label est_sig   se      confint
1 welf_supp.BY.gvslvol    1.00 0.00 [1.00, 1.00]
2 welf_supp.BY.gvslvue 0.73*** 0.05 [0.64, 0.82]
3 welf_supp.BY.gvhlthc 0.94*** 0.04 [0.87, 1.02]
4 welf_supp.BY.gvcldcr 0.89*** 0.04 [0.81, 0.98]
5 welf_supp.BY.gvjbevn 1.04*** 0.05 [0.94, 1.15]
6 welf_supp.BY.gvpdlwk 0.91*** 0.04 [0.83, 1.00]

# we can also take a look at residual variances 
# using lavaan terminology, the residual variances are the thetas

theta <- round(inspect(fit_ws_6, "est")$theta,3)
theta.std <- round(inspect(fit_ws_6, "std")$theta,3) 
r2 <- round(inspect(fit_ws_6, "r2"),3)

data.frame(row.names = c(),                       # empty the columns names 
           Variables = colnames(theta),           # variable names 
           "Residuals" = diag(theta),             # diagonal theta
           "Std. Residuals" = diag(theta.std),    # diagonal std. theta
           "R Squared" = r2                       # R-squared
        )

  Variables Residuals Std..Residuals R.Squared
1   gvslvol     0.923          0.420     0.580
2   gvslvue     3.001          0.817     0.183
3   gvhlthc     1.031          0.476     0.524
4   gvcldcr     2.087          0.672     0.328
5   gvjbevn     3.755          0.731     0.269
6   gvpdlwk     2.032          0.657     0.343

7.7 Model Fit statistics

We have two different type of fit statistics.

1. Global fit measures:
   • They take into account how the entire entire model fit the data
     
   • (some) rules of thumb: \(CFI/TLI > 0.95, RMSEA < 0.05, SRMR < 0.06\)
     
   • current practice is: chi-square value + df + pvalue, RMSEA, CFI and SRMR
   • DO NOT cherry pick your fit indices
   • Check Hu and Bentler (1999)
2. “Local” fit measures:
   • looking at just one part of the model
   • usually involves the use of Modification Indexes
   • Check Thoemmes, Rosseel, and Textor (2018)

## Global fit measures ##
# we can select which global fit measures to extract
fitMeasures(fit_ws_6, c("logl","AIC", "BIC", "chisq", "df", "pvalue", "cfi", "tli","rmsea"), output = "matrix")

                 
logl   -19550.568
aic     39125.136
bic     39190.676
chisq     326.956
df          9.000
pvalue      0.000
cfi         0.878
tli         0.797
rmsea       0.142

Q: How our model perform in terms of global fit measures?

We can have a better grasp of what is happening taking a look at the local fit of our model.

## Local fit measures: modification indices ##
mi <- inspect(fit_ws_6,"mi")
mi.sorted <- mi[order(-mi$mi),] # sort from high to low mi.sorted[1:5,] # only display some large MI values
mi.sorted[1:5,] # only display some large MI values

       lhs op     rhs      mi    epc sepc.lv sepc.all sepc.nox
27 gvcldcr ~~ gvpdlwk 243.943  0.918   0.918    0.446    0.446
15 gvslvol ~~ gvhlthc 168.789  0.621   0.621    0.637    0.637
23 gvhlthc ~~ gvcldcr  58.187 -0.375  -0.375   -0.256   -0.256
21 gvslvue ~~ gvjbevn  47.839  0.607   0.607    0.181    0.181
25 gvhlthc ~~ gvpdlwk  39.471 -0.310  -0.310   -0.214   -0.214

# let's plot the modification indices 
plot(mi.sorted$mi) # plot the MI values
abline(h=3.84) # add a horizontal reference line (chisq value for 1 df where p=0.05)

Q: What the plot is suggesting ?

1. In certain contexts, we can modify the model based on a review of:
   • MI’s in combination with EPC’s (Expected Value Change). Both need to be “substantial”
   • Theory or the source of the data (e.g. review the content of the test items)
2. However, modifying a CFA moves it away from a strictly confirmatory model
   • The more modifications, the more exploratory the model becomes
   • Maybe this model was not ready for a confirmatory modelling strategy?