Path Diagrams

library(lavaan)

## This is lavaan 0.6-19
## lavaan is FREE software! Please report any bugs.

library(tidygraph)

## 
## Attaching package: 'tidygraph'

## The following object is masked from 'package:stats':
## 
##     filter

library(ggraph)

## Loading required package: ggplot2

library(igraph)

## 
## Attaching package: 'igraph'

## The following object is masked from 'package:tidygraph':
## 
##     groups

## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum

## The following object is masked from 'package:base':
## 
##     union

library(grid) 
library(semPlot)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:igraph':
## 
##     as_data_frame, groups, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

### Simulating data to fit lavaan models

set.seed(2025)
dummy_data <- as.data.frame(matrix(rnorm(900 * 20), ncol = 20))
colnames(dummy_data) <- c(
  "MRelQual","WRelQual",
  "MDistress","MEsteem","MLifeSat",
  "WDistress","WEsteem","WLifeSat",
  "MHS", "MBS", "WHS", "WBS", "MRelDep", "WRelDep", "MSES", "WSES", "RelSES" 
)

names(dummy_data) <- make.names(names(dummy_data), unique = TRUE)
# mean centering data to reduce multicollinearity
center_colmeans <- function(x) {
    xcenter = colMeans(x)
    x - rep(xcenter, rep.int(nrow(x), ncol(x)))
}

dummy_data<- center_colmeans(dummy_data)

# creating factor variable for employment where 0 = unemployed and 1 = employed
set.seed(123)   
n <- 900             
p_emp <- 0.95         # probability of being employed, nz unemployment rate is currently 5%
Memp_num <- rbinom(n, size = 1, prob = p_emp)
dummy_data$M_Emp <- factor(Memp_num, ordered=TRUE,
                     levels = c(0, 1)) # 0 = unemployed, 1 = employed
Wemp_num <- rbinom(n, size = 1, prob = p_emp)
dummy_data$W_Emp<- factor(Wemp_num, ordered=TRUE,
                     levels = c(0, 1)) # 0 = unemployed, 1 = employed

prop.table(table(dummy_data$M_Emp)) #checking that this works, it did, 95% of peoople in our simulated dataset are employed

## 
##          0          1 
## 0.05333333 0.94666667

# creating new variable for within-couple employment
dummy_data <- dummy_data %>%
  mutate(
    Couple_Emp = case_when(
      M_Emp == 1  & W_Emp == 1  ~ 1,
      M_Emp == 0 & W_Emp == 0 ~ 2,
      M_Emp == 1  & W_Emp == 0 ~ 3,
      M_Emp == 0 & W_Emp == 1  ~ 4 # man unemployed, woman employed
    ),
    Couple_Emp = factor(Couple_Emp, ordered = TRUE, levels = c(1:4
    ))
  )

dummy_data %>% count(Couple_Emp) %>% mutate(freq = n / sum(n))

##   Couple_Emp   n        freq
## 1          1 809 0.898888889
## 2          2   2 0.002222222
## 3          3  43 0.047777778
## 4          4  46 0.051111111

# if assuming that 95% of the population is employed, very small numbers of unemployed men and employed women. 

Model 1. Baseline Relative Deprivation

model_syntax <- '
  MWellbeing =~ MDistress + MEsteem + MLifeSat 
  WWellbeing =~ WDistress + WEsteem + WLifeSat

  MWellbeing ~ MRelDep + WRelDep
  WWellbeing ~ MRelDep + WRelDep

  MRelQual   ~ MRelDep + WRelDep
  WRelQual   ~ MRelDep + WRelDep

  MRelDep ~~ WRelDep
  MWellbeing ~~ WWellbeing
  MRelQual   ~~ WRelQual
  MWellbeing ~~ MRelQual
  WWellbeing ~~ WRelQual
'

# =~ is factor loadings
# ~ is structual equations
# ~~ is covariance 

fit1 <- lavaan::lavaan(model_syntax, data = dummy_data, 
                       fixed.x=FALSE,    
                       auto.var=TRUE,    
                       auto.fix.first=FALSE)

## Warning: lavaan->lav_model_vcov():  
##    Could not compute standard errors! The information matrix could not be 
##    inverted. This may be a symptom that the model is not identified.

## Warning: lavaan->lav_object_post_check():  
##    some estimated lv variances are negative

#Plot
semPaths(fit1,                   
         layout = "tree2",        
         rotation=2,             
         edge.label.cex = 0.8,    
         sizeMan = 5,              
         sizeLat = 7,              
         nCharNodes = 0,           
         residuals = FALSE,        
         intercepts = FALSE)  
title("Model 1. Baseline Relative Deprivation Model", 
      cex.main = 1.2,   
      line = 2)

Model 2. Baseline SES Model

\[ \text{Relative SES} \;=\; \frac{\text{Man's SES}}{\text{Man's SES} \;+\; \text{Woman's SES}} \]

model_syntax.2 <- '
  MWellbeing =~ MDistress + MEsteem + MLifeSat
  WWellbeing =~ WDistress + WEsteem + WLifeSat

  MWellbeing  ~ MSES + WSES + RelSES
  WWellbeing  ~ MSES + WSES + RelSES
  MRelQual    ~ MSES + WSES + RelSES
  WRelQual    ~ MSES + WSES + RelSES

  MSES   ~~ WSES
  MSES   ~~ RelSES
  WSES   ~~ RelSES

  MWellbeing ~~ WWellbeing
  MWellbeing ~~ MRelQual
  WWellbeing ~~ WRelQual
  MRelQual   ~~ WRelQual
'

fit2 <- lavaan::lavaan(model_syntax.2, data = dummy_data, 
                       fixed.x=FALSE,    
                       auto.var=TRUE,    
                       auto.fix.first=FALSE)

## Warning: lavaan->lav_model_vcov():  
##    Could not compute standard errors! The information matrix could not be 
##    inverted. This may be a symptom that the model is not identified.

## Warning: lavaan->lav_object_post_check():  
##    some estimated ov variances are negative

#Plot
semPaths(fit2,                   
         layout = "tree2",        
         rotation=2,             
         edge.label.cex = 0.8,    
         sizeMan = 5,              
         sizeLat = 7,              
         nCharNodes = 0,           
         residuals = FALSE,        
         intercepts = FALSE)  
title("Model 2. Baseline SES Model", 
      cex.main = 1.2,   
      line = 2)

Model 3. Baseline Employment Model.

model_syntax.3 <- '
  MWellbeing =~ MDistress + MEsteem + MLifeSat
  WWellbeing =~ WDistress + WEsteem + WLifeSat

  MWellbeing ~ M_Emp + W_Emp + Couple_Emp
  WWellbeing ~ M_Emp + W_Emp + Couple_Emp

  MRelQual   ~ M_Emp + W_Emp + Couple_Emp
  WRelQual   ~ M_Emp + W_Emp + Couple_Emp

  M_Emp   ~~ W_Emp
  M_Emp   ~~ Couple_Emp
  W_Emp   ~~ Couple_Emp

  MWellbeing ~~ WWellbeing
  MRelQual   ~~ WRelQual
  MWellbeing ~~ MRelQual
  WWellbeing ~~ WRelQual
'

fit3 <- lavaan::lavaan(model_syntax.3, data = dummy_data, 
                       fixed.x=FALSE,    
                       auto.var=TRUE,    
                       auto.fix.first=FALSE)

## Warning: lavaan->lav_partable_check():  
##    automatically added intercepts are set to zero: ("MRelQual", "WRelQual")

## Warning: lavaan->lav_samplestats_step2():  
##    correlation between variables Couple_Emp and M_Emp is (nearly) 1.0

## Warning: lavaan->lav_model_vcov():  
##    Could not compute standard errors! The information matrix could not be 
##    inverted. This may be a symptom that the model is not identified.

## Warning: lavaan->lav_test_satorra_bentler():  
##    could not invert information matrix needed for robust test statistic

#Plot
semPaths(fit3,                   
         layout = "tree2",        
         rotation=2,             
         edge.label.cex = 0.9,
         label.cex=1.2,
         sizeMan = 5,              
         sizeLat = 7,              
         nCharNodes = 0,           
         residuals = FALSE,        
         intercepts = FALSE)    
title("Model 3. Baseline Employment Model", 
      cex.main = 1.2,   
      line = 2)

Model 4. Relative Deprivation with Sexism as a Moderator

# Creating interaction terms
dummy_data$MInt_BS <- dummy_data$MRelDep * dummy_data$MBS
dummy_data$MInt_HS <- dummy_data$MRelDep * dummy_data$MHS
dummy_data$WInt_BS<- dummy_data$WRelDep * dummy_data$WBS
dummy_data$WInt_HS<- dummy_data$WRelDep * dummy_data$WHS

model_syntax.4 <- '
  MWellbeing =~ MDistress + MEsteem + MLifeSat
  WWellbeing =~ WDistress + WEsteem + WLifeSat

  # Relative‐deprivation paths (actor + partner)
  MWellbeing ~ MRelDep + WRelDep
  WWellbeing ~ MRelDep + WRelDep

  MRelQual   ~ MRelDep + WRelDep
  WRelQual   ~ MRelDep + WRelDep

  # Sexism actor paths
  MWellbeing ~ MBS + MHS
  MRelQual   ~ MBS + MHS

  WWellbeing ~ WBS + WHS
  WRelQual   ~ WBS + WHS

  # Interactions actor
  MWellbeing ~ MInt_BS + MInt_HS
  MRelQual   ~ MInt_BS + MInt_HS

  WWellbeing ~ WInt_BS + WInt_HS
  WRelQual   ~ WInt_BS + WInt_HS

  # Covariances
  MRelDep    ~~ WRelDep
  MWellbeing ~~ WWellbeing
  MRelQual   ~~ WRelQual
  MWellbeing ~~ MRelQual
  WWellbeing ~~ WRelQual
  MBS ~~ MHS
  MBS ~~ WBS
  MBS ~~ WHS
  MHS ~~ WBS
  MHS ~~ WHS
  WBS ~~ WHS
  MInt_BS ~~ MRelDep
  MInt_BS~~ MBS
  MInt_HS ~~ MRelDep
  MInt_HS~~ MHS
  WInt_BS ~~ WRelDep
  WInt_BS~~ WBS
  WInt_HS ~~ WRelDep
  WInt_HS~~ WHS
'


fit4 <- lavaan::lavaan(model_syntax.4, data = dummy_data, 
                       fixed.x=FALSE,    
                       auto.var=TRUE,    
                       auto.fix.first=FALSE)

## Warning: lavaan->lav_model_vcov():  
##    Could not compute standard errors! The information matrix could not be 
##    inverted. This may be a symptom that the model is not identified.

## Warning: lavaan->lav_object_post_check():  
##    some estimated lv variances are negative

#Plot
semPaths(fit4,                   
         layout = "tree2",        
         rotation=2,             
         edge.label.cex = 0.8,    
         sizeMan = 5,              
         sizeLat = 7,              
         nCharNodes = 0,           
         residuals = FALSE,        
         intercepts = FALSE)    
title("Model 4. Relative Deprivation with Sexism as Moderator", 
      cex.main = 1.2,   
      line = 2)

lavInspect(fit4, "npar") # no. of parameters

## [1] 69

Model 5. SES with Sexism as Moderator

dummy_data$SESxMBS <- dummy_data$RelSES * dummy_data$MBS
dummy_data$SESxMHS <- dummy_data$RelSES * dummy_data$MHS
dummy_data$SESxWBS<- dummy_data$RelSES * dummy_data$WBS
dummy_data$SESxWHS<- dummy_data$RelSES * dummy_data$WHS

model_syntax.5 <- '
  MWellbeing =~ MDistress + MEsteem + MLifeSat
  WWellbeing =~ WDistress + WEsteem + WLifeSat

  # SES regressions (actor + partner)
  MWellbeing ~ MSES + WSES + RelSES
  WWellbeing ~ MSES + WSES + RelSES

  MRelQual   ~ MSES + WSES + RelSES
  WRelQual   ~ MSES + WSES + RelSES

  # Sexism actor‐only
  MWellbeing ~ MBS + MHS
  MRelQual   ~ MBS + MHS

  WWellbeing ~ WBS + WHS
  WRelQual   ~ WBS + WHS

  # Interactions (RelSES × Sexism), actor‐only
  MWellbeing ~ SESxMBS + SESxMHS
  MRelQual   ~ SESxMBS + SESxMHS

  WWellbeing ~ SESxWBS + SESxWHS
  WRelQual   ~ SESxWBS + SESxWHS

  # Covariances among SES IVs
  MSES  ~~ WSES
  MSES  ~~ RelSES
  WSES  ~~ RelSES
  MBS ~~ MHS
  MBS ~~ WBS
  MBS ~~ WHS
  MHS ~~ WBS
  MHS ~~ WHS
  WBS ~~ WHS
  SESxMBS ~~ RelSES
  SESxMBS ~~ MBS
  SESxMHS ~~ RelSES
  SESxMHS ~~ MHS
  SESxWBS ~~ RelSES
  SESxWBS ~~ WBS
  SESxWHS ~~ RelSES
  SESxWHS ~~ WHS
  MWellbeing ~~ WWellbeing
  MRelQual   ~~ WRelQual
  MWellbeing ~~ MRelQual
  WWellbeing ~~ WRelQual
'

fit5 <- lavaan::lavaan(model_syntax.5, data = dummy_data, 
                       fixed.x=FALSE,    
                       auto.var=TRUE,    
                       auto.fix.first=FALSE)

## Warning: lavaan->lav_model_vcov():  
##    Could not compute standard errors! The information matrix could not be 
##    inverted. This may be a symptom that the model is not identified.

## Warning: lavaan->lav_object_post_check():  
##    some estimated ov variances are negative

## Warning: lavaan->lav_object_post_check():  
##    some estimated lv variances are negative

#Plot
semPaths(fit5,                   
         layout = "tree2",        
         rotation=2,             
         edge.label.cex = 0.8,    
         sizeMan = 5,              
         sizeLat = 7,              
         nCharNodes = 0,           
         residuals = FALSE,        
         intercepts = FALSE)  
title("Model 5. SES with Sexism as Moderator", 
      cex.main = 1.2,   
      line = 2)

lavInspect(fit5, "npar") # no. of parameters

## [1] 76

(ninfo.fit5<- (17*(17+1))/2)

## [1] 153

Model 6. Employment with Sexism as Moderator

# creating interaction terms (need to decide how to code employment status, so for now I'm just going to use men's employment as a 0-1 dummy variable)

dummy_data <- dummy_data %>%
  mutate(
    M_Emp.n      = as.numeric(as.character(M_Emp)),
    W_Emp.n      = as.numeric(as.character(W_Emp)),
    Couple_Emp.n = as.numeric(as.character(Couple_Emp))
  )
dummy_data$EmpxMBS <- dummy_data$M_Emp.n * dummy_data$MBS
dummy_data$EmpxMHS <- dummy_data$M_Emp.n * dummy_data$MHS
dummy_data$EmpxWBS<- dummy_data$M_Emp.n * dummy_data$WBS
dummy_data$EmpxWHS<- dummy_data$M_Emp.n * dummy_data$WHS


model_syntax.6 <- '
  # Measurement models
  MWellbeing =~ MDistress + MEsteem + MLifeSat
  WWellbeing =~ WDistress + WEsteem + WLifeSat

  # Regressions on employment
  MWellbeing ~ M_Emp.n + W_Emp.n + Couple_Emp.n
  WWellbeing ~ M_Emp.n + W_Emp.n + Couple_Emp.n

  MRelQual   ~ M_Emp.n + W_Emp.n + Couple_Emp.n
  WRelQual   ~ M_Emp.n + W_Emp.n + Couple_Emp.n

  # Sexism actor‐only effects
  MWellbeing ~ MBS + MHS
  MRelQual   ~ MBS + MHS

  WWellbeing ~ WBS + WHS
  WRelQual   ~ WBS + WHS

  # Employment×Sexism interactions (actor only)
  MWellbeing ~ EmpxMBS + EmpxMHS
  MRelQual   ~ EmpxMBS + EmpxMHS

  WWellbeing ~ EmpxWBS + EmpxWHS
  WRelQual   ~ EmpxWBS + EmpxWHS

  M_Emp.n ~~ W_Emp.n
  M_Emp.n ~~ Couple_Emp.n
  W_Emp.n ~~ Couple_Emp.n
  EmpxMBS ~~ M_Emp.n
  EmpxMBS ~~ MBS
  EmpxMHS ~~ M_Emp.n
  EmpxMHS ~~ MHS
  EmpxWBS ~~ M_Emp.n
  EmpxWBS ~~ WBS
  EmpxWHS ~~ M_Emp.n
  EmpxWHS ~~ WHS
  MBS ~~ MHS
  MBS ~~ WBS
  MBS ~~ WHS
  MHS ~~ WBS
  MHS ~~ WHS
  WBS ~~ WHS

  # Latent DV covariances
  MWellbeing ~~ WWellbeing
  MRelQual   ~~ WRelQual
  MWellbeing ~~ MRelQual
  WWellbeing ~~ WRelQual
'


fit6 <- lavaan::lavaan(
  model          = model_syntax.6,
  data           = dummy_data,
  fixed.x        = FALSE,
  auto.var       = TRUE,
  auto.fix.first = FALSE
)

## Warning: lavaan->lav_lavaan_step11_estoptim():  
##    Model estimation FAILED! Returning starting values.

semPaths(fit6,
         what       = "path",
         layout     = "tree2",
         rotation   = 2,
         reorder    = FALSE,
         sizeMan    = 5,
         sizeLat    = 7,
         nCharNodes = 0,
         residuals  = FALSE,
         intercepts = FALSE,
         label.cex  = 1.4)

title("Model 6. Employment with Sexism as Moderator", 
      cex.main = 1.2,   
      line = 2)

lavInspect(fit6, "npar") # no. of parameters

## [1] 76

```