SEM
Lucy Bryson
2025-06-09
Structural Equation Modelling
This document will define the model of resource access and its impact on eating behaviours, exploring both SES and Social Resources as mediators for eating behaviours in the Covid-19 pandemic.
library(lavaan)
library(MASS)
library(semPlot)
library(dagitty)
library(tidyverse)
library(haven)wave1 <- read_sav("covid-19_wave1_survey_cls.sav") |> mutate(wave = 1)
wave2 <- read_sav("covid-19_wave2_survey_cls.sav") |> mutate(wave = 2)
wave3 <- read_sav("covid-19_wave3_survey_cls.sav") |> mutate(wave = 3)
long_data1 <- bind_rows(wave1, wave2, wave3)
#Reverse coding SOCPROV_3
long_data1$CW1_SOCPROV_3 <- 4 - long_data1$CW1_SOCPROV_3
long_data1$CW2_SOCPROV_3 <- 4 - long_data1$CW2_SOCPROV_3
long_data1$CW3_SOCPROV_3 <- 4 - long_data1$CW3_SOCPROV_3
#converting labelled variables to numeric
long_data1 <- long_data1 |> mutate(across(where(is.labelled), as.numeric))
#converting CW1_EXCISESP to a consistent scale
long_data1 <- long_data1 |>
mutate(CW1_EXCISESP_scaled = scales::rescale(CW1_EXCISESP, to = c(0, 7)))
#removing unit errors or misunderstandings from CW1_FRTVEGSP
long_data1 <- long_data1 |>
mutate(
CW1_FRTVEGSP_trim = ifelse(CW1_FRTVEGSP > 20, NA, CW1_FRTVEGSP)
)
#reverse coding CW1_FOODAFFORD
long_data1 <- long_data1 |>
mutate(
ses_proxy = case_when(
CW1_FOODAFFORD == 1 ~ 4,
CW1_FOODAFFORD == 2 ~ 3,
CW1_FOODAFFORD == 3 ~ 2,
CW1_FOODAFFORD == 4 ~ 1,
TRUE ~ NA_real_
)
)model_sem <- '
# Measurement models
socialprov_t1 =~ CW1_SOCPROV_1 + CW1_SOCPROV_2 + CW1_SOCPROV_3
socialprov_t2 =~ CW2_SOCPROV_1 + CW2_SOCPROV_2 + CW2_SOCPROV_3
socialprov_t3 =~ CW3_SOCPROV_1 + CW3_SOCPROV_2 + CW3_SOCPROV_3
# Longitudinal autoregressive effects
socialprov_t2 ~ socialprov_t1 + ses_proxy
socialprov_t3 ~ socialprov_t2 + ses_proxy
CW2_FRTVEGSP ~ CW1_FRTVEGSP_trim + socialprov_t1 + ses_proxy
CW3_FRTVEGSP ~ CW2_FRTVEGSP + socialprov_t2 + ses_proxy
CW2_EXCISESP ~ CW1_EXCISESP_scaled + socialprov_t1 + ses_proxy
CW3_EXCISESP ~ CW2_EXCISESP + socialprov_t2 + ses_proxy
'fit_sem <- lavaan::sem(
model_sem,
data = long_data1,
missing = "fiml",
fixed.x = FALSE
)
summary(fit_sem, fit.measures = TRUE, standardized = TRUE)## lavaan 0.6-19 ended normally after 162 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 63
##
## Used Total
## Number of observations 58862 67518
## Number of missing patterns 59
##
## Model Test User Model:
##
## Test statistic 292.404
## Degrees of freedom 89
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## Test statistic 17288.816
## Degrees of freedom 117
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.988
## Tucker-Lewis Index (TLI) 0.984
##
## Robust Comparative Fit Index (CFI) NA
## Robust Tucker-Lewis Index (TLI) NA
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -294136.523
## Loglikelihood unrestricted model (H1) -293990.321
##
## Akaike (AIC) 588399.046
## Bayesian (BIC) 588964.972
## Sample-size adjusted Bayesian (SABIC) 588764.757
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.006
## 90 Percent confidence interval - lower 0.005
## 90 Percent confidence interval - upper 0.007
## P-value H_0: RMSEA <= 0.050 1.000
## P-value H_0: RMSEA >= 0.080 0.000
##
## Robust RMSEA NA
## 90 Percent confidence interval - lower NA
## 90 Percent confidence interval - upper NA
## P-value H_0: Robust RMSEA <= 0.050 NA
## P-value H_0: Robust RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.075
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Observed
## Observed information based on Hessian
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## socialprov_t1 =~
## CW1_SOCPROV_1 1.000 0.259 0.618
## CW1_SOCPROV_2 1.332 0.044 29.997 0.000 0.345 0.764
## CW1_SOCPROV_3 0.911 0.028 32.135 0.000 0.236 0.564
## socialprov_t2 =~
## CW2_SOCPROV_1 1.000 0.308 0.649
## CW2_SOCPROV_2 1.270 0.029 43.851 0.000 0.392 0.802
## CW2_SOCPROV_3 0.889 0.020 45.525 0.000 0.274 0.555
## socialprov_t3 =~
## CW3_SOCPROV_1 1.000 0.319 0.682
## CW3_SOCPROV_2 1.204 0.027 44.362 0.000 0.384 0.790
## CW3_SOCPROV_3 0.920 0.021 44.527 0.000 0.293 0.598
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## socialprov_t2 ~
## socialprov_t1 0.000 94.099 0.000 1.000 0.000 0.000
## ses_proxy 0.187 NA 0.606 0.269
## socialprov_t3 ~
## socialprov_t2 -0.001 31.234 -0.000 1.000 -0.001 -0.001
## ses_proxy 0.000 96.302 0.000 1.000 0.000 0.000
## CW2_FRTVEGSP ~
## CW1_FRTVEGSP_t 0.624 27.310 0.023 0.982 0.624 0.522
## socialprov_t1 0.000 405.887 0.000 1.000 0.000 0.000
## ses_proxy -2.259 18.290 -0.124 0.902 -2.259 -0.446
## CW3_FRTVEGSP ~
## CW2_FRTVEGSP 0.000 NA 0.000 0.000
## socialprov_t2 -0.000 NA -0.000 -0.000
## ses_proxy -0.000 NA -0.000 -0.000
## CW2_EXCISESP ~
## CW1_EXCISESP_s 1.281 73.780 0.017 0.986 1.281 0.506
## socialprov_t1 0.000 492.839 0.000 1.000 0.000 0.000
## ses_proxy -1.645 5.917 -0.278 0.781 -1.645 -0.335
## CW3_EXCISESP ~
## CW2_EXCISESP 0.000 108.914 0.000 1.000 0.000 0.000
## socialprov_t2 0.000 807.808 0.000 1.000 0.000 0.000
## ses_proxy -0.000 596.483 -0.000 1.000 -0.000 -0.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .socialprov_t3 ~~
## .CW3_FRTVEGSP -0.059 NA -0.185 -0.077
## .CW3_EXCISESP -0.075 0.077 -0.977 0.329 -0.235 -0.105
## .CW3_FRTVEGSP ~~
## .CW3_EXCISESP 0.745 0.061 12.248 0.000 0.745 0.140
## ses_proxy ~~
## CW1_FRTVEGSP_t 0.106 0.007 15.348 0.000 0.106 0.127
## CW1_EXCISESP_s 0.037 0.003 11.958 0.000 0.037 0.097
## CW1_FRTVEGSP_trim ~~
## CW1_EXCISESP_s 0.325 0.014 24.054 0.000 0.325 0.201
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .CW1_SOCPROV_1 1.184 0.005 231.634 0.000 1.184 2.825
## .CW1_SOCPROV_2 1.191 0.006 215.860 0.000 1.191 2.639
## .CW1_SOCPROV_3 1.145 0.005 223.296 0.000 1.145 2.741
## .CW2_SOCPROV_1 0.542 NA 0.542 1.142
## .CW2_SOCPROV_2 0.326 NA 0.326 0.667
## .CW2_SOCPROV_3 0.581 NA 0.581 1.178
## .CW3_SOCPROV_1 1.242 362.392 0.003 0.997 1.242 2.655
## .CW3_SOCPROV_2 1.221 436.446 0.003 0.998 1.221 2.511
## .CW3_SOCPROV_3 1.212 333.283 0.004 0.997 1.212 2.470
## .CW2_FRTVEGSP 9.980 NA 9.980 4.443
## .CW3_FRTVEGSP 3.938 NA 3.938 1.647
## .CW2_EXCISESP 7.657 NA 7.657 3.510
## .CW3_EXCISESP 3.056 2283.456 0.001 0.999 3.056 1.370
## ses_proxy 3.770 0.003 1077.686 0.000 3.770 8.497
## CW1_FRTVEGSP_t 3.839 0.015 251.231 0.000 3.839 2.041
## CW1_EXCISESP_s 1.199 0.007 173.212 0.000 1.199 1.391
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .CW1_SOCPROV_1 0.108 0.003 39.094 0.000 0.108 0.618
## .CW1_SOCPROV_2 0.085 0.004 21.625 0.000 0.085 0.416
## .CW1_SOCPROV_3 0.119 0.003 44.438 0.000 0.119 0.681
## .CW2_SOCPROV_1 0.130 0.003 50.356 0.000 0.130 0.578
## .CW2_SOCPROV_2 0.085 0.003 25.535 0.000 0.085 0.357
## .CW2_SOCPROV_3 0.169 0.003 62.552 0.000 0.169 0.692
## .CW3_SOCPROV_1 0.117 0.003 44.150 0.000 0.117 0.535
## .CW3_SOCPROV_2 0.089 0.003 27.679 0.000 0.089 0.376
## .CW3_SOCPROV_3 0.155 0.003 54.327 0.000 0.155 0.642
## .CW2_FRTVEGSP 2.963 129.843 0.023 0.982 2.963 0.587
## .CW3_FRTVEGSP 5.714 NA 5.714 1.000
## .CW2_EXCISESP 3.162 139.749 0.023 0.982 3.162 0.665
## .CW3_EXCISESP 4.978 0.164 30.436 0.000 4.978 1.000
## socialprov_t1 0.067 0.003 21.321 0.000 1.000 1.000
## .socialprov_t2 0.088 NA 0.928 0.928
## .socialprov_t3 0.102 0.008 13.242 0.000 1.000 1.000
## ses_proxy 0.197 0.002 89.680 0.000 0.197 1.000
## CW1_FRTVEGSP_t 3.537 0.041 86.899 0.000 3.537 1.000
## CW1_EXCISESP_s 0.743 0.008 88.022 0.000 0.743 1.000semPaths(
object = fit_sem, # your lavaan fit object
what = "std", # standardized estimates
whatLabels = "std", # display standardized values
style = "lisrel", # cleaner layout
layout = "tree", # top-down structure
edge.label.cex = 0.9, # size of the path coefficients
sizeMan = 5, # size of manifest variables
sizeLat = 7, # size of latent variables
edge.color = "blue", # path color
residuals = FALSE, # hide residuals for clarity
intercepts = FALSE, # hide intercepts
nCharNodes = 0, # show full variable names
optimizeLatRes = TRUE, # avoid overlap
rotation = 2, # flip layout if needed
title = FALSE, # remove default title
asize = 1.5 # arrow size
)
library(dagitty)
dag <- dagitty("
dag {
ses_proxy -> eating_t1
ses_proxy -> eating_t2
ses_proxy -> eating_t3
ses_proxy -> exercise_t1
ses_proxy -> exercise_t2
ses_proxy -> exercise_t3
socialprov_t1 -> eating_t1
socialprov_t1 -> exercise_t1
socialprov_t2 -> eating_t2
socialprov_t2 -> exercise_t2
socialprov_t3 -> eating_t3
socialprov_t3 -> exercise_t3
socialprov_t1 -> socialprov_t2
socialprov_t2 -> socialprov_t3
eating_t1 -> eating_t2
eating_t2 -> eating_t3
exercise_t1 -> exercise_t2
exercise_t2 -> exercise_t3
}
")
# plot
plot(dag)## Plot coordinates for graph not supplied! Generating coordinates, see ?coordinates for how to set your own.
library(DiagrammeR)## Warning: package 'DiagrammeR' was built under R version 4.4.3grViz("
digraph SEM {
# node styles
node [shape = box, style = filled, color = lightgray, fontname = Helvetica]
# SES and socialprov nodes
ses [label = 'SES (ses_proxy)', shape = ellipse, color = lightblue]
sp1 [label = 'Social Provision T1', shape = ellipse, color = lightgreen]
sp2 [label = 'Social Provision T2', shape = ellipse, color = lightgreen]
sp3 [label = 'Social Provision T3', shape = ellipse, color = lightgreen]
# Eating and exercise observed variables
eat1 [label = 'Eating T1']
ex1 [label = 'Exercise T1']
eat2 [label = 'Eating T2']
ex2 [label = 'Exercise T2']
eat3 [label = 'Eating T3']
ex3 [label = 'Exercise T3']
# Paths T1
sp1 -> eat1
ses -> eat1
sp1 -> ex1
ses -> ex1
# Paths T2
sp2 -> eat2
ses -> eat2
sp2 -> ex2
ses -> ex2
# Paths T3
sp3 -> eat3
ses -> eat3
sp3 -> ex3
ses -> ex3
# autoregressive paths
eat1 -> eat2
eat2 -> eat3
ex1 -> ex2
ex2 -> ex3
# socialprov continuity
sp1 -> sp2
sp2 -> sp3
}
")