Mediation Analysis with Lavaan | MENAS

1 What we are going to cover

  1. MIMIC models
  2. Mediation Analysis

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)

In addition, we will use some other variables

  • 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 Belgium)

  • 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

# ### MODELING ###
# install.packages("lavaan")                # SEM modelling

# ### VISUALIZATION ###
# install.packages("tidySEM")               # plotting SEM models

# Load the packages 

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

### MODELING ###
library("lavaan")       

### VISUALIZATION ###
library("tidySEM")

4 Data exploration

It is always a good practice to check that everything is in order and make sense of the data that we are going to analyse.

ess_df <- haven::read_sav("https://github.com/albertostefanelli/SEM_labs/raw/master/data/ESS4_belgium.sav")


ess_df_selected <- ess_df %>% select(
                  ## Welfare support items ##
                  gvslvol, # the old
                  gvslvue, # the unemployed
                  gvhlthc, # the sick
                  gvcldcr, # working parents
                  gvjbevn, # job for everyone
                  gvpdlwk, # paid sick leave  
                  ##    Economic criticism items ##
                  sbstrec, # strain on economy
                  sbbsntx, # too much taxes
                  ##    Moral criticism items ##
                  sblazy,  # people lazy 
                  sblwcoa, # care for others
                  sblwlka,  # look after others
                  ## Egalitarianism ##
                  gincdif,
                  dfincac,
                  smdfslv,
                  ## Demographics ##
                  agea, 
                  eduyrs,
                  gndr,
                  hinctnta


)

descriptive_ess <- as.data.frame(psych::describe(ess_df_selected))

descriptive_ess <- dplyr::select(descriptive_ess, 
  n,
  mean,
  sd,
  median,
  min,
  max,
  skew,
  kurtosis)

descriptive_ess
            n      mean         sd median min max        skew    kurtosis
gvslvol  1759  7.871518  1.4811986      8   0  10 -0.78122430  1.37621274
gvslvue  1753  6.059897  1.9176408      6   0  10 -0.40989634  0.40781515
gvhlthc  1757  8.030734  1.4715101      8   0  10 -0.94788987  1.81361522
gvcldcr  1748  7.276316  1.7622804      7   0  10 -0.80981813  1.34667552
gvjbevn  1756  6.234624  2.2667339      7   0  10 -0.48699173 -0.08302998
gvpdlwk  1753  7.310325  1.7677965      8   0  10 -0.82684920  1.31394820
sbstrec  1735  2.963112  1.0193172      3   1   5  0.10310828 -0.92138250
sbbsntx  1714  2.531505  0.9970179      2   1   5  0.42532526 -0.51782206
sblazy   1751  2.894917  1.0557772      3   1   5  0.10529513 -0.92067158
sblwcoa  1751  2.919475  1.0352766      3   1   5  0.06846687 -1.05077060
sblwlka  1746  2.981672  1.0339269      3   1   5 -0.09401127 -1.07498593
gincdif  1751  2.233010  1.0590918      2   1   5  0.73154748 -0.24601836
dfincac  1756  2.625854  1.0544131      2   1   5  0.50378434 -0.57304293
smdfslv  1752  2.472603  0.9744553      2   1   5  0.65984232 -0.19176740
agea     1760 46.456818 18.7300429     46  15 105  0.20358225 -0.81002487
eduyrs   1759 12.666856  3.6579256     12   0  30  0.01431515  0.74369643
gndr     1760  1.509091  0.5000594      2   1   2 -0.03633866 -1.99981479
hinctnta 1567  7.456924  2.3668693      8   1  10 -0.70485395 -0.57439987

Q: Is everything ok ?

4.1 Sum score indicators

ess_df$welf_supp <- ess_df %>% 
  dplyr::select(starts_with("gv")) %>%  
  rowwise() %>% 
  dplyr::mutate(sum = sum(cur_data(),na.rm=T)/ncol(.)) %>%  
  select(sum) %>% 
  unlist() %>% as.vector()

ess_df$welf_crit <- ess_df %>% 
  dplyr::select(starts_with("sb")) %>%  
  rowwise() %>% 
  dplyr::mutate(sum = sum(cur_data(),na.rm=T)/ncol(.)) %>%  
  select(sum) %>% 
  unlist() %>% as.vector()

ess_df$egal <- ess_df %>% 
  dplyr::select(gincdif, dfincac,smdfslv ) %>%  
  rowwise() %>% 
  dplyr::mutate(sum = sum(cur_data(),na.rm=T)/ncol(.)) %>%  
  select(sum) %>% 
  unlist() %>% as.vector()

5 MIMIC model

Now that we constructed our indicators, we can apply our theoretical knowledge and test some simple hypotheses. We hypothesise that respondent’ structural characteristics influence their support for the welfare state. These type of models are called MIMIC models and stands for “Multiple Indicators, Multiple Causes”.

model_ws_mimic_1 <-'
welf_supp ~ welf_crit + gndr + eduyrs + hinctnta + agea
'

fit_ws_mimic_1 <- sem(model_ws_mimic_1, # model formula
                   data=ess_df 
  )

summary(fit_ws_mimic_1, 
        standardized=TRUE)
lavaan 0.6-10 ended normally after 1 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                         6
                                                      
                                                  Used       Total
  Number of observations                          1566        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

Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  welf_supp ~                                                           
    welf_crit         0.286    0.061    4.664    0.000    0.286    0.119
    gndr              0.112    0.062    1.794    0.073    0.112    0.045
    eduyrs           -0.002    0.009   -0.163    0.871   -0.002   -0.005
    hinctnta         -0.021    0.014   -1.461    0.144   -0.021   -0.040
    agea              0.003    0.002    1.917    0.055    0.003    0.050

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .welf_supp         1.507    0.054   27.982    0.000    1.507    0.981
r2 <- round(inspect(fit_ws_mimic_1, "r2"),3)
r2*100
welf_supp 
      1.9 
# residual variance Welfare support
sem_res_dv <- lavTech(fit_ws_mimic_1, "est", add.labels = TRUE)$psi["welf_supp","welf_supp"]
sem_res_dv
[1] 1.506806
# sample-implied variance Welfare Support
sample_var <- lavTech(fit_ws_mimic_1, "sampstat", add.labels = TRUE)[[1]]$cov["welf_supp","welf_supp"]
sample_var
[1] 1.535885

5.1 Comparison between the ML and least square estimator

Let’s compare the two different regression approaches

  1. Are the coefficients the same (magnitude and significance)?
  2. Is the number of degree of freedom the same across the two estimators?
  3. Are we estimating the same number of parameters?
  4. Are the residual variances equal?
lm_fit_ws_mimic_1 <- lm(model_ws_mimic_1, # model formula
                       data=ess_df, 
  )

summary(lm_fit_ws_mimic_1)

Call:
lm(formula = model_ws_mimic_1, data = ess_df)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.2669 -0.7337  0.0366  0.7607  3.0914 

Coefficients:
             Estimate Std. Error t value             Pr(>|t|)    
(Intercept)  6.206366   0.257738  24.080 < 0.0000000000000002 ***
welf_crit    0.286288   0.061495   4.655           0.00000351 ***
gndr         0.111576   0.062312   1.791               0.0736 .  
eduyrs      -0.001530   0.009410  -0.163               0.8709    
hinctnta    -0.021158   0.014507  -1.458               0.1449    
agea         0.003434   0.001795   1.913               0.0559 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.23 on 1560 degrees of freedom
  (194 observations deleted due to missingness)
Multiple R-squared:  0.01893,   Adjusted R-squared:  0.01579 
F-statistic: 6.021 on 5 and 1560 DF,  p-value: 0.0000158

Let’s compare the estimated coefficients.

sem_tidy_results <- table_results(fit_ws_mimic_1,             
  columns = c("label", 
              "est_sig"),
  digits = 3,
) %>% filter(str_detect(label, "welf_supp.ON") 
                        )
sem_tidy_results$est_sig_lm <- round(coef(lm_fit_ws_mimic_1)[-1],3)

sem_tidy_results
                   label  est_sig est_sig_lm
1 welf_supp.ON.welf_crit 0.286***      0.286
2      welf_supp.ON.gndr    0.112      0.112
3    welf_supp.ON.eduyrs   -0.002     -0.002
4  welf_supp.ON.hinctnta   -0.021     -0.021
5      welf_supp.ON.agea    0.003      0.003

The regression coefficients are identical (good!).

One thing to note is that we don’t have an intercept in the lavaan output. This highlights an important difference between ls and ML. SEM, and more in general ML, focuses on the covariance structure of the data. However, we can include an intercept (often called mean in SEM).

model_ws_mimic_1 <-'
welf_supp ~ 1 + welf_crit + gndr + eduyrs + hinctnta + agea
'

fit_ws_mimic_1 <- sem(model_ws_mimic_1, # model formula
                   data=ess_df 
  )

summary(fit_ws_mimic_1, 
        standardized=TRUE)
lavaan 0.6-10 ended normally after 21 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                         7
                                                      
                                                  Used       Total
  Number of observations                          1566        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

Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  welf_supp ~                                                           
    welf_crit         0.286    0.061    4.664    0.000    0.286    0.119
    gndr              0.112    0.062    1.794    0.073    0.112    0.045
    eduyrs           -0.002    0.009   -0.163    0.871   -0.002   -0.005
    hinctnta         -0.021    0.014   -1.461    0.144   -0.021   -0.040
    agea              0.003    0.002    1.917    0.055    0.003    0.050

Intercepts:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .welf_supp         6.206    0.257   24.126    0.000    6.206    5.008

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .welf_supp         1.507    0.054   27.982    0.000    1.507    0.981

Typically we include intercepts only when it’s relevant to our scientific questions. For example, when we fit multi-group models. It is worth noting that the number of free parameters increased by 1 because we were estimating an additional intercept but this does not change the degrees of freedom because we are simply adding another term to the total parameters. However, we have 7 estimated parameters and 21 piece of informations in the variance-covariance matrix so we should have a total of 8 degree of freedom instead of 0 as reported in the lavaan output.

Let’s try to understand what it is happening by looking at the model- and sample-implied covariance matrices.

# Returns the observed covariance matrix.
lavInspect(fit_ws_mimic_1, "sampstat")$cov
          wlf_sp  wlf_cr  gndr    eduyrs  hnctnt  agea   
welf_supp   1.536                                        
welf_crit   0.068   0.267                                
gndr        0.030   0.002   0.250                        
eduyrs     -0.074   0.274  -0.083  13.325                
hinctnta   -0.136   0.106  -0.077   3.391   5.589        
agea        0.916  -1.625   0.102 -14.753 -10.002 330.596
# Returns the model estimated covariance matrix.
fitted(fit_ws_mimic_1, "theta")$cov
          wlf_sp  wlf_cr  gndr    eduyrs  hnctnt  agea   
welf_supp   1.536                                        
welf_crit   0.068   0.267                                
gndr        0.030   0.002   0.250                        
eduyrs     -0.074   0.274  -0.083  13.325                
hinctnta   -0.136   0.106  -0.077   3.391   5.589        
agea        0.916  -1.625   0.102 -14.753 -10.002 330.596

There is no difference between the sample implied and the model-implied covariance matrices. This means that the fit function is equal to 0 and we are able to perfectly reproduce the structure of data. Given the fact that we have 0 degree of freedom, it is likely that we are fitting a fully saturated model where all the exogenous variables are correlated to each other. Let’s check this by looking at what parameters are being estimated.

lavTech(fit_ws_mimic_1, "est", add.labels = TRUE)
$lambda
          welf_supp welf_crit gndr eduyrs hinctnta agea
welf_supp         1         0    0      0        0    0
welf_crit         0         1    0      0        0    0
gndr              0         0    1      0        0    0
eduyrs            0         0    0      1        0    0
hinctnta          0         0    0      0        1    0
agea              0         0    0      0        0    1

$theta
          welf_supp welf_crit gndr eduyrs hinctnta agea
welf_supp         0         0    0      0        0    0
welf_crit         0         0    0      0        0    0
gndr              0         0    0      0        0    0
eduyrs            0         0    0      0        0    0
hinctnta          0         0    0      0        0    0
agea              0         0    0      0        0    0

$psi
          welf_supp    welf_crit         gndr      eduyrs     hinctnta        agea
welf_supp  1.506806  0.000000000  0.000000000   0.0000000   0.00000000   0.0000000
welf_crit  0.000000  0.266691092  0.001544433   0.2741784   0.10602824  -1.6253434
gndr       0.000000  0.001544433  0.249996330  -0.0825394  -0.07702471   0.1021638
eduyrs     0.000000  0.274178382 -0.082539403  13.3252664   3.39055504 -14.7525628
hinctnta   0.000000  0.106028244 -0.077024706   3.3905550   5.58937772 -10.0020698
agea       0.000000 -1.625343445  0.102163797 -14.7525628 -10.00206985 330.5964093

$beta
          welf_supp welf_crit     gndr      eduyrs    hinctnta        agea
welf_supp         0 0.2862882 0.111576 -0.00152983 -0.02115809 0.003434201
welf_crit         0 0.0000000 0.000000  0.00000000  0.00000000 0.000000000
gndr              0 0.0000000 0.000000  0.00000000  0.00000000 0.000000000
eduyrs            0 0.0000000 0.000000  0.00000000  0.00000000 0.000000000
hinctnta          0 0.0000000 0.000000  0.00000000  0.00000000 0.000000000
agea              0 0.0000000 0.000000  0.00000000  0.00000000 0.000000000

$nu
          intercept
welf_supp         0
welf_crit         0
gndr              0
eduyrs            0
hinctnta          0
agea              0

$alpha
          intercept
welf_supp  6.206366
welf_crit  2.639527
gndr       1.498084
eduyrs    12.747765
hinctnta   7.459770
agea      46.996169

Let’s go through the output

  1. \(\Lambda\) factor loading matrix. In this case, it is empty because we are not estimating any latent factor.
  2. \(\Theta\) variance-covariance matrix of the residuals. Since the sample implied and model implied covariance matrices are the same, this is equal to 0.
  3. \(\Phi\) model implied variance-covariance matrix. We can observe that both variances and covariances of our predictors are estimated! 4 \(\Beta\) matrix represents the regression path between endogenous and exogenous variables.
  4. \(\nu\) matrix indicates whether or not we estimate the intercept for the observed variables,
  5. \(\alpha\) matrix is the sample-implied intercept matrix.

Let’s try to include the covariances among the predictors in the model syntax and check the estimated number of free parameters.

model_ws_mimic_2 <-'
welf_supp ~ 1 + welf_crit + gndr + eduyrs + hinctnta + agea

# covariance 
welf_crit ~~ gndr
welf_crit ~~ eduyrs
welf_crit ~~ hinctnta
welf_crit ~~ agea
gndr ~~ eduyrs
gndr ~~ hinctnta
gndr ~~ agea
eduyrs ~~ hinctnta
eduyrs ~~ agea
hinctnta ~~ agea

## res. variances 
welf_supp ~~ welf_supp
welf_crit ~~ welf_crit
gndr  ~~ gndr
eduyrs ~~ eduyrs
hinctnta ~~ hinctnta
agea ~~ agea

'

fit_ws_mimic_2 <- sem(model_ws_mimic_2, # model formula
                   data=ess_df      # data frame
  )


summary(fit_ws_mimic_2, standardized=TRUE)
lavaan 0.6-10 ended normally after 84 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                        27
                                                      
                                                  Used       Total
  Number of observations                          1566        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

Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  welf_supp ~                                                           
    welf_crit         0.286    0.061    4.664    0.000    0.286    0.119
    gndr              0.112    0.062    1.794    0.073    0.112    0.045
    eduyrs           -0.002    0.009   -0.163    0.871   -0.002   -0.005
    hinctnta         -0.021    0.014   -1.461    0.144   -0.021   -0.040
    agea              0.003    0.002    1.917    0.055    0.003    0.050

Covariances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  welf_crit ~~                                                          
    gndr              0.002    0.007    0.237    0.813    0.002    0.006
    eduyrs            0.274    0.048    5.696    0.000    0.274    0.145
    hinctnta          0.106    0.031    3.424    0.001    0.106    0.087
    agea             -1.625    0.241   -6.750    0.000   -1.625   -0.173
  gndr ~~                                                               
    eduyrs           -0.083    0.046   -1.788    0.074   -0.083   -0.045
    hinctnta         -0.077    0.030   -2.573    0.010   -0.077   -0.065
    agea              0.102    0.230    0.445    0.657    0.102    0.011
  eduyrs ~~                                                             
    hinctnta          3.391    0.234   14.470    0.000    3.391    0.393
    agea            -14.753    1.718   -8.586    0.000  -14.753   -0.222
  hinctnta ~~                                                           
    agea            -10.002    1.115   -8.968    0.000  -10.002   -0.233

Intercepts:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .welf_supp         6.206    0.257   24.126    0.000    6.206    5.008
    welf_crit         2.640    0.013  202.264    0.000    2.640    5.111
    gndr              1.498    0.013  118.567    0.000    1.498    2.996
    eduyrs           12.748    0.092  138.195    0.000   12.748    3.492
    hinctnta          7.460    0.060  124.865    0.000    7.460    3.155
    agea             46.996    0.459  102.284    0.000   46.996    2.585

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .welf_supp         1.507    0.054   27.982    0.000    1.507    0.981
    welf_crit         0.267    0.010   27.982    0.000    0.267    1.000
    gndr              0.250    0.009   27.982    0.000    0.250    1.000
    eduyrs           13.325    0.476   27.982    0.000   13.325    1.000
    hinctnta          5.589    0.200   27.982    0.000    5.589    1.000
    agea            330.596   11.815   27.982    0.000  330.596    1.000

So, we are fitting a fully saturated model where both the residual variances and the covariances are estimated.

semPlot::semPaths(fit_ws_mimic_1, title = FALSE, curvePivot = TRUE)

5.2 Obtain the same residual variance between the LM and SEM model [OPTIONAL]

There is a discrepancy in what is reported in the lm() and sem() output concerning the residuals between. First we need to convert the residual standard error in residual variances to be allowed to compare them.

# model rank 
k <- length(lm_fit_ws_mimic_1$coefficients)-1 
# residual sum of square
RSS <- sum(lm_fit_ws_mimic_1$residuals**2)
# sample size 
n <- length(lm_fit_ws_mimic_1$residuals)
# residual standard error
res_se <- sqrt(RSS/(n-(1+k))) 
# residual variance
lm_res_dv <- round(res_se^2,4) 
lm_res_dv
[1] 1.5126
data.frame(row.names = c(),                       # empty the columns names 
           models = c("MIMIC","lm"),
           res = c(sem_res_dv, lm_res_dv)
)
  models      res
1  MIMIC 1.506806
2     lm 1.512600

We can notice a discrepancy in the residuals between the two estimation procedures. Are the SEM and the LM model equivalent? Let’s try to fix the covariances among the predictors to zero and see if the residual variance of Welfare Support changes.

model_ws_mimic_3 <-'

welf_supp ~ 1 + welf_crit + gndr + eduyrs + hinctnta + agea

# intercepts 
welf_crit ~ 0
gndr ~ 0
eduyrs ~ 0
hinctnta ~ 0
agea ~ 0


# covariance 
gndr ~~ 0*eduyrs
gndr ~~ 0*hinctnta
gndr ~~ 0*agea
eduyrs ~~ 0*hinctnta
eduyrs ~~ 0*agea
hinctnta ~~ 0*agea


'

fit_ws_mimic_3 <- sem(model_ws_mimic_3, # model formula
                   data=ess_df, # data frame
  )


summary(fit_ws_mimic_3, 
        standardized=TRUE
        )
lavaan 0.6-10 ended normally after 66 iterations

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

Parameter Estimates:

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

Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  welf_supp ~                                                           
    welf_crit         0.286    0.012   24.823    0.000    0.286    0.520
    gndr              0.112    0.020    5.681    0.000    0.112    0.119
    eduyrs           -0.002    0.002   -0.654    0.513   -0.002   -0.014
    hinctnta         -0.021    0.004   -5.338    0.000   -0.021   -0.112
    agea              0.003    0.001    5.579    0.000    0.003    0.117

Covariances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  gndr ~~                                                               
    eduyrs            0.000                               0.000    0.000
    hinctnta          0.000                               0.000    0.000
    agea              0.000                               0.000    0.000
  eduyrs ~~                                                             
    hinctnta          0.000                               0.000    0.000
    agea              0.000                               0.000    0.000
  hinctnta ~~                                                           
    agea              0.000                               0.000    0.000

Intercepts:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .welf_supp         6.206    0.031  200.080    0.000    6.206    4.195
    welf_crit         0.000                               0.000    0.000
    gndr              0.000                               0.000    0.000
    eduyrs            0.000                               0.000    0.000
    hinctnta          0.000                               0.000    0.000
    agea              0.000                               0.000    0.000

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .welf_supp         1.507    0.054   27.982    0.000    1.507    0.689
    welf_crit         7.234    0.259   27.982    0.000    7.234    1.000
    gndr              2.494    0.089   27.982    0.000    2.494    1.000
    eduyrs          175.831    6.284   27.982    0.000  175.831    1.000
    hinctnta         61.238    2.188   27.982    0.000   61.238    1.000
    agea           2539.236   90.745   27.982    0.000 2539.236    1.000

The residuals are still different. This is not due to model specification but rather to how the residual variance is calculated with ML estimators.

For OLS: \(\hat{\sigma}^2_{OLS} = \frac{\sum_{i=1}^{N} \hat{\zeta_i}^2}{N-k}R\)

For ML: \(\hat{\sigma}^2_{ML} = \frac{\sum_{i=1}^{N} \hat{\zeta_i}^2}{N}\).

So to make the two variance comparable: \(\hat{\sigma}^2_{ML} = \left(\frac{N-k}{n}\right) \hat{\sigma}^2_{OLS}\)

sample_size <- nrow(lm_fit_ws_mimic_1$model)
# number of estimated parameters 
rank <- 6
((sample_size-rank)/sample_size)*lm_res_dv
[1] 1.506805

The lm estimate is now equal to the one we obtain from the ML.

Take home:

  1. The magnitude and the significance of the fixed effects of an ordinary least squares regression is equivalent to maximum likelihood.
  2. The number of degree of freedom is calculated using the number of observations - 1 in lm. In SEM, we use the unique pieces of info in the observed covariance matrix - the number of estimated parameters.
  3. In lavaan, regression models are always fully saturated. That is, we are also estimating the covariances between all the included predictors
  4. ML and ordinary least squares estimators calculate the residual variance in different ways.

6 Mediation analysis

In his most simple form, mediation analysis (or path analysis) tests whether the relationship between two variables is explained by a third intermediate variable. It can have a casual interpretation such as the extent to which a variable (mediator) participates in the transmittance of change from a cause to its effect. Consider a classical mediation setup with three variables:

  • Y is the dependent variable (Welfare support)
  • X is the predictor (Income)
  • M is a mediator (Egalitarianism)

This results in different paths

  1. a path: Test whether X and M are significantly associated
  2. b path: Test whether M and Y are significantly associated
  3. c path: Test whether X and Y are significantly associated (Direct Effect)
  4. c’ path: Test whether Y from X are significantly associated after controlling for M (Indirect Effect). This is usually called “the amount of mediation”.

Note that the Total Effect is equal to Direct Effect + Indirect Effect or \(c= ab +c'\). It can be interpreted as the pairwise correlation between the two variables.

Before getting into mediation analysis, let’s take a closer look at the lavaan syntax for path modelling:

  • regression ~ (is regressed on)
  • (residual) (co)variance ~~ (is correlated with)
  • intercept ~1 (intercept)
  • new model parameter := (is equal to)
model_mediation <- '

## Direct effect ##
welf_supp ~ c*hinctnta

## Mediator ##
egal ~ a*hinctnta
welf_supp ~ b*egal

## Indirect effect (a*b) ##
ab := a*b
## Total effect ##
total := c + (a*b)

'

fit_mediation <- sem(model_mediation, # model formula
           data=ess_df                # data frame
  )

summary(fit_mediation, standardized=TRUE)
lavaan 0.6-10 ended normally after 9 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                         5
                                                      
                                                  Used       Total
  Number of observations                          1567        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

Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  welf_supp ~                                                           
    hinctnta   (c)   -0.017    0.013   -1.305    0.192   -0.017   -0.033
  egal ~                                                                
    hinctnta   (a)    0.024    0.006    3.977    0.000    0.024    0.100
  welf_supp ~                                                           
    egal       (b)   -0.321    0.054   -5.904    0.000   -0.321   -0.148

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .welf_supp         1.499    0.054   27.991    0.000    1.499    0.976
   .egal              0.323    0.012   27.991    0.000    0.323    0.990

Defined Parameters:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
    ab               -0.008    0.002   -3.298    0.001   -0.008   -0.015
    total            -0.025    0.013   -1.884    0.060   -0.025   -0.048

The indirect effect is significant and negative. We can say that egalitarianism mediates the effect between income and welfare support. However, the total effect is still not significant.

Let’s plot our model to get better grasp of what it is happening.

# let's organize our plot on 4 rows 
# this help our readers by having a more comprehensible plot

lay <- get_layout(
"gincdif", "dfincac", "smdfslv", "", 
"", "egal", "", "", 
"hinctnta", "", "welf_supp", "", 
"",  "gvslvol", "gvslvue", "gvhlthc",
rows = 4)

plot_mediation <- graph_sem(model = fit_mediation,   # model fit
          layout = lay,        # layout
          angle = 170          # adjust the arrows 
          #label = "est_std"   # get standardized results (not rounded)
          )

plot_mediation

Q: The path between welfare support and income (direct effect) is not significant. Nor the total effect. What does that mean ?

A: We have a so called indirect-only mediation such that the indirect effect pathways fully account for the overall impact of IV on DV with the direct effect being insignificant

Zhao, Lynch and Chen (2010) classify mediation effects as following:

  • Complementary mediation: Mediated effect (a x b) and direct effect (c) both exist and point at the same direction.
  • Competitive mediation: Mediated effect (a x b) and direct effect (c) both exist and point in opposite directions.
  • Indirect-only mediation: Mediated effect (a x b) exists, but no direct effect (c).
  • Direct-only non-mediation: Direct effect (c) exists, but no indirect effect.
  • No-effect non-mediation: Nether direct effect (c), nor indirect effect exists.

6.1 More complex path modelling [serial mediation]

There are situation where we might want to test chain linking of the mediators, with a specified direction flow. For instance, we might hypothesise that the effect of income on Welfare Support is mediated by both Egalitarianism and Welfare Criticism.

model_mediation_serial <- '

## Direct effect ##
welf_supp ~ c*hinctnta
welf_supp ~ d*welf_crit

## Mediator ##
egal ~ a*hinctnta
#welf_supp ~ b*egal
welf_crit ~ e*egal

## Indirect effect ##
#ab := a*b
aed := a*e*d
## Total effect ##
total := c + (a*e*d) 



'

fit_mediation_serial <- sem(model_mediation_serial, # model formula
           data = ess_df                # data frame
  )

summary(fit_mediation_serial, standardized=TRUE)
lavaan 0.6-10 ended normally after 9 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                         7
                                                      
                                                  Used       Total
  Number of observations                          1567        1760
                                                                  
Model Test User Model:
                                                      
  Test statistic                                50.201
  Degrees of freedom                                 2
  P-value (Chi-square)                           0.000

Parameter Estimates:

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

Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  welf_supp ~                                                           
    hinctnta   (c)   -0.030    0.013   -2.270    0.023   -0.030   -0.057
    welf_crit  (d)    0.270    0.060    4.503    0.000    0.270    0.113
  egal ~                                                                
    hinctnta   (a)    0.024    0.006    3.977    0.000    0.024    0.100
  welf_crit ~                                                           
    egal       (e)    0.114    0.023    4.998    0.000    0.114    0.125

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .welf_supp         1.513    0.054   27.991    0.000    1.513    0.984
   .egal              0.323    0.012   27.991    0.000    0.323    0.990
   .welf_crit         0.264    0.009   27.991    0.000    0.264    0.984

Defined Parameters:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
    aed               0.001    0.000    2.560    0.010    0.001    0.001
    total            -0.029    0.013   -2.214    0.027   -0.029   -0.055

Income has a negative total effect on Welfare Support (\(\beta_{Total}=-0.029, -2.214, p<.05\)). As theorized, this effect was serially mediated by Egalitarianism and Welfare Criticism. The indirect pathway of the effect of Income on Welfare Support via Egalitarianism and Welfare Criticism is significant and positive (\(\beta_{aed}= .001, z=2.560, p<.05\)). On the contrary, the direct effect of Income on Welfare support is negative (\(\beta_{c}=-0.030, -2.270, p<.05\)) suggesting the presence of a competitive mediation.

semPlot::semPaths(fit_mediation_serial, 
                  title = FALSE, 
                  what= 'col', 
                  whatLabels = 'std',
                  layout= "tree2", 
                  curvePivot = TRUE, 
                  #color=c('black'),
                  edge.color=c('black'))

?semPaths

6.2 Comparing nested models

In some cases, we want to test whether adding a path between 2 predictors leads to a improvement in the fit of the model. Remember that fully saturated models have a perfect fit so it makes no sense to compare them.

model_mediation_serial_2 <- '

## Direct effect ##
welf_supp ~ c*hinctnta
welf_supp ~ d*welf_crit
welf_supp ~ f*egal

## Mediator ##
egal ~ a*hinctnta
#welf_supp ~ b*egal
welf_crit ~ e*egal

## Indirect effect ##
#ab := a*b
aed := a*e*d
fa := f*a
## Total effect ##
total := c + (a*e*d) 

'

fit_mediation_serial_2 <- sem(model_mediation_serial_2, # model formula
           data = ess_df                # data frame
  )

summary(fit_mediation_serial_2, 
        standardized=TRUE
        )
lavaan 0.6-10 ended normally after 9 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                         8
                                                      
                                                  Used       Total
  Number of observations                          1567        1760
                                                                  
Model Test User Model:
                                                      
  Test statistic                                 8.023
  Degrees of freedom                                 1
  P-value (Chi-square)                           0.005

Parameter Estimates:

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

Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  welf_supp ~                                                           
    hinctnta   (c)   -0.022    0.013   -1.696    0.090   -0.022   -0.042
    welf_crit  (d)    0.316    0.060    5.301    0.000    0.316    0.132
    egal       (f)   -0.355    0.054   -6.533    0.000   -0.355   -0.164
  egal ~                                                                
    hinctnta   (a)    0.024    0.006    3.977    0.000    0.024    0.100
  welf_crit ~                                                           
    egal       (e)    0.114    0.023    4.998    0.000    0.114    0.125

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .welf_supp         1.473    0.053   27.991    0.000    1.473    0.958
   .egal              0.323    0.012   27.991    0.000    0.323    0.990
   .welf_crit         0.264    0.009   27.991    0.000    0.264    0.984

Defined Parameters:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
    aed               0.001    0.000    2.684    0.007    0.001    0.002
    fa               -0.009    0.003   -3.397    0.001   -0.009   -0.016
    total            -0.021    0.013   -1.629    0.103   -0.021   -0.040

Q: Why we do have 1 degree of freedom ?

Let’s extract the fit indices of these two model and see whether there is an improvement in the fit.

fitm_model_serial <-  fitMeasures(fit_mediation_serial, c("logl",
                                                "AIC", 
                                                "BIC", 
                                                "chisq",
                                                "df",
                                                "pvalue",
                                                "cfi",
                                                "srmr"), 
                                                output = "matrix")

fitm_model_serial_2 <- fitMeasures(fit_mediation_serial_2, c("logl",
                                                   "AIC", 
                                                   "BIC", 
                                                   "chisq",
                                                   "df",
                                                   "pvalue",
                                                   "cfi",
                                                   "srmr"), 
                                                   output = "matrix")



data.frame(Fit=rownames(fitm_model_serial),
           "serial" = round(fitm_model_serial[,1],2),
           "plus_f" = round(fitm_model_serial_2[,1],2)

  
)
          Fit   serial   plus_f
logl     logl -5064.39 -5043.30
aic       aic 10142.78 10102.60
bic       bic 10180.28 10145.46
chisq   chisq    50.20     8.02
df         df     2.00     1.00
pvalue pvalue     0.00     0.00
cfi       cfi     0.55     0.94
srmr     srmr     0.06     0.02
lavaan::anova(fit_mediation_serial, fit_mediation_serial_2)
Chi-Squared Difference Test

                       Df   AIC   BIC   Chisq Chisq diff Df diff       Pr(>Chisq)    
fit_mediation_serial_2  1 10103 10146  8.0234                                        
fit_mediation_serial    2 10143 10180 50.2012     42.178       1 0.00000000008334 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Adding and removing path changes the model and, consequently, the interpretation of the results. Make sure to think through this.

6.3 Boostrapping

The delta method used for testing mediation is known to be problematic because the sampling distribution of the indirect path product terms is not normal. Bootstrapping is a common workaround for this problem that does not make strong assumptions about the distribution of the coefficient of interest (i.e., the sampling distributions of the two mediated paths). We can implement this using the argument se = "bootstrap".

fit_mediation_serial_bs <- sem(model_mediation_serial_2, # model formula
           data = ess_df,
           se = "bootstrap"
  )

summary(fit_mediation_serial_bs, standardized=TRUE)
lavaan 0.6-10 ended normally after 9 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of model parameters                         8
                                                      
                                                  Used       Total
  Number of observations                          1567        1760
                                                                  
Model Test User Model:
                                                      
  Test statistic                                 8.023
  Degrees of freedom                                 1
  P-value (Chi-square)                           0.005

Parameter Estimates:

  Standard errors                            Bootstrap
  Number of requested bootstrap draws             1000
  Number of successful bootstrap draws            1000

Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  welf_supp ~                                                           
    hinctnta   (c)   -0.022    0.014   -1.619    0.105   -0.022   -0.042
    welf_crit  (d)    0.316    0.081    3.883    0.000    0.316    0.132
    egal       (f)   -0.355    0.058   -6.154    0.000   -0.355   -0.164
  egal ~                                                                
    hinctnta   (a)    0.024    0.006    4.051    0.000    0.024    0.100
  welf_crit ~                                                           
    egal       (e)    0.114    0.024    4.665    0.000    0.114    0.125

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .welf_supp         1.473    0.063   23.300    0.000    1.473    0.958
   .egal              0.323    0.012   26.511    0.000    0.323    0.990
   .welf_crit         0.264    0.012   22.382    0.000    0.264    0.984

Defined Parameters:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
    aed               0.001    0.000    2.250    0.024    0.001    0.002
    fa               -0.009    0.003   -3.392    0.001   -0.009   -0.016
    total            -0.021    0.014   -1.554    0.120   -0.021   -0.040

7 !!Support Ukraine!!

References

Zhao, Xinshu, John G Lynch Jr, and Qimei Chen. 2010. “Reconsidering Baron and Kenny: Myths and Truths about Mediation Analysis.” Journal of Consumer Research 37 (2): 197–206.

  1. Alberto Stefanelli↩︎