Path Analysis and Mediation – Analysis Walkthrough
Every model, statistic, and console output used in the lecture slides
Frederick Anyan – PSY3100
13 October 2026
- 1 Load data
- 2 A first look at
lavaansyntax - 3 Getting bootstrap confidence intervals in R
- 4 Model 1: does provocation increase aggression through hostile cognition?
- 5 Model 2: does provocation also work through negative affect?
- 6 Model 3: does hostile cognition itself intensify negative affect?
- 7 Model 4: does trait irritability matter too?
- 8 Model fit: cutoff conventions, applied to our own models
- 9 Model Identification: Just, Over and Under Identified models
- 10 Fixed and constrained parameters
- 11 Multi-group analysis: does the pathway differ by gender?
- 12 The
manymomepackage - 13 Post-estimation: drawing the fitted model
- 14 Session info
This notebook works through the Path Analysis and Mediation lecture in R, in the same order as the slides. Every model, test and number shown in the lecture is computed here from the data, so you can run each step yourself and see where the numbers come from.
We use the gam_mediation.csv (Models 1-4, \(N=300\)) and
gam_multigroup.csv (the multi-group section, \(N=300\)). The theoretical framework
throughout is the General Aggression Model (Anderson & Bushman,
2002): provocation and trait irritability as person/situational inputs,
hostile cognition and negative affect as internal-state routes,
aggression as the outcome.
1 Load data
d <- read.csv("../source_material/data/gam_mediation.csv")
str(d)## 'data.frame': 300 obs. of 6 variables:
## $ id : int 1 2 3 4 5 6 7 8 9 10 ...
## $ provocation : num 19.1 20.1 17.8 25.4 24.4 19.9 23.5 10.8 11.5 12.6 ...
## $ trait_irritability: num 34.4 20.9 38.6 44.1 32.2 40.2 28.4 23.3 23.4 37.7 ...
## $ hostile_cognition : num 2.2 14 19.5 24.2 21.6 14.8 16.5 14.5 8.7 10 ...
## $ neg_affect : num 10 28.8 12.3 26.2 33.1 27.1 28.1 23.7 25.2 19.3 ...
## $ aggression : num 21.5 33.5 47.4 67.8 53.9 45 29.1 38.2 21.9 25.7 ...d %>%
select(-id) %>%
pivot_longer(everything(), names_to = "Variable", values_to = "value") %>%
group_by(Variable) %>%
summarise(N = n(), Mean = mean(value), SD = sd(value),
Min = min(value), Max = max(value), .groups = "drop") %>%
kable(digits = 2)| Variable | N | Mean | SD | Min | Max |
|---|---|---|---|---|---|
| aggression | 300 | 40.03 | 14.92 | 0.0 | 73.8 |
| hostile_cognition | 300 | 15.00 | 4.99 | 0.0 | 26.9 |
| neg_affect | 300 | 22.09 | 6.75 | 10.0 | 41.0 |
| provocation | 300 | 20.00 | 6.00 | 3.2 | 36.7 |
| trait_irritability | 300 | 30.00 | 7.00 | 11.0 | 50.0 |
cor(d[,-1]) %>% kable(digits = 2)| provocation | trait_irritability | hostile_cognition | neg_affect | aggression | |
|---|---|---|---|---|---|
| provocation | 1.00 | 0.28 | 0.50 | 0.43 | 0.49 |
| trait_irritability | 0.28 | 1.00 | 0.34 | 0.13 | 0.31 |
| hostile_cognition | 0.50 | 0.34 | 1.00 | 0.46 | 0.50 |
| neg_affect | 0.43 | 0.13 | 0.46 | 1.00 | 0.51 |
| aggression | 0.49 | 0.31 | 0.50 | 0.51 | 1.00 |
2 A first look at
lavaan syntax
model_demo <- '
Y ~ c*X # direct effect
M ~ a*X # X predicts the mediator
Y ~ b*M # M predicts the outcome
ab := a*b # indirect effect, defined
total := c + (a*b) # total effect, defined
'
cat(model_demo)##
## Y ~ c*X # direct effect
## M ~ a*X # X predicts the mediator
## Y ~ b*M # M predicts the outcome
## ab := a*b # indirect effect, defined
## total := c + (a*b) # total effect, defined3 Getting bootstrap confidence intervals in R
Every indirect and total-indirect effect in this lecture uses a
nonparametric bootstrap CI (percentile method), not the delta method,
since a product of two estimates is generally not symmetrically (i.e.,
normally) distributed. se = "bootstrap" plus
parameterEstimates(..., boot.ci.type = "perc") is the
pattern used throughout below.
4 Model 1: does provocation increase aggression through hostile cognition?
m1 <- '
hostile_cognition ~ a*provocation
aggression ~ cp*provocation + b*hostile_cognition
indirect := a*b
total := cp + (a*b)
'
fit1 <- sem(m1, data = d, se = "bootstrap", bootstrap = 5000, iseed = 4171)
summary(fit1, standardized = TRUE, fit.measures = TRUE, rsquare = TRUE)## lavaan 0.6-21 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 5
##
## Number of observations 300
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Model Test Baseline Model:
##
## Test statistic 205.759
## Degrees of freedom 3
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.000
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -2040.384
## Loglikelihood unrestricted model (H1) -2040.384
##
## Akaike (AIC) 4090.768
## Bayesian (BIC) 4109.287
## Sample-size adjusted Bayesian (SABIC) 4093.430
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.000
## P-value H_0: RMSEA <= 0.050 NA
## P-value H_0: RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.000
##
## Parameter Estimates:
##
## Standard errors Bootstrap
## Number of requested bootstrap draws 5000
## Number of successful bootstrap draws 4978
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## hostile_cognition ~
## provocatn (a) 0.416 0.043 9.577 0.000 0.416 0.501
## aggression ~
## provocatn (cp) 0.793 0.124 6.400 0.000 0.793 0.319
## hstl_cgnt (b) 1.022 0.151 6.763 0.000 1.022 0.342
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .hostile_cogntn 18.584 1.510 12.310 0.000 18.584 0.749
## .aggression 149.085 11.074 13.463 0.000 149.085 0.672
##
## R-Square:
## Estimate
## hostile_cogntn 0.251
## aggression 0.328
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## indirect 0.425 0.082 5.196 0.000 0.425 0.171
## total 1.218 0.115 10.561 0.000 1.218 0.490parameterEstimates(fit1, boot.ci.type = "perc", level = .95) %>%
filter(op == ":=") %>% kable(digits = 3)| lhs | op | rhs | label | est | se | z | pvalue | ci.lower | ci.upper |
|---|---|---|---|---|---|---|---|---|---|
| indirect | := | a*b | indirect | 0.425 | 0.082 | 5.196 | 0 | 0.274 | 0.600 |
| total | := | cp+(a*b) | total | 1.218 | 0.115 | 10.561 | 0 | 0.996 | 1.448 |
5 Model 2: does provocation also work through negative affect?
m2 <- '
hostile_cognition ~ a1*provocation
neg_affect ~ a2*provocation
aggression ~ cp*provocation + b1*hostile_cognition + b2*neg_affect
hostile_cognition ~~ neg_affect
ind1 := a1*b1
ind2 := a2*b2
total_ind := ind1 + ind2
'
fit2 <- sem(m2, data = d, se = "bootstrap", bootstrap = 5000, iseed = 4171)
summary(fit2, standardized = TRUE, fit.measures = TRUE, rsquare = TRUE)## lavaan 0.6-21 ended normally after 11 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 9
##
## Number of observations 300
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Model Test Baseline Model:
##
## Test statistic 327.093
## Degrees of freedom 6
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.000
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -2977.686
## Loglikelihood unrestricted model (H1) -2977.686
##
## Akaike (AIC) 5973.372
## Bayesian (BIC) 6006.706
## Sample-size adjusted Bayesian (SABIC) 5978.163
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.000
## P-value H_0: RMSEA <= 0.050 NA
## P-value H_0: RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.000
##
## Parameter Estimates:
##
## Standard errors Bootstrap
## Number of requested bootstrap draws 5000
## Number of successful bootstrap draws 4945
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## hostile_cognition ~
## provocatn (a1) 0.416 0.043 9.570 0.000 0.416 0.501
## neg_affect ~
## provocatn (a2) 0.486 0.058 8.338 0.000 0.486 0.433
## aggression ~
## provocatn (cp) 0.599 0.129 4.649 0.000 0.599 0.241
## hstl_cgnt (b1) 0.743 0.161 4.620 0.000 0.743 0.248
## neg_affct (b2) 0.638 0.122 5.212 0.000 0.638 0.289
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .hostile_cognition ~~
## .neg_affect 8.140 1.635 4.980 0.000 8.140 0.311
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .hostile_cogntn 18.584 1.509 12.314 0.000 18.584 0.749
## .neg_affect 36.886 2.623 14.061 0.000 36.886 0.813
## .aggression 135.521 9.992 13.563 0.000 135.521 0.611
##
## R-Square:
## Estimate
## hostile_cogntn 0.251
## neg_affect 0.187
## aggression 0.389
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## ind1 0.309 0.078 3.949 0.000 0.309 0.124
## ind2 0.310 0.073 4.263 0.000 0.310 0.125
## total_ind 0.619 0.095 6.500 0.000 0.619 0.249parameterEstimates(fit2, boot.ci.type = "perc", level = .95) %>%
filter(op %in% c(":=", "~~")) %>% kable(digits = 3)| lhs | op | rhs | label | est | se | z | pvalue | ci.lower | ci.upper |
|---|---|---|---|---|---|---|---|---|---|
| hostile_cognition | neg_affect | 8.140 | 1.635 | 4.980 | 0 | 4.931 | 11.325 | ||
| hostile_cognition | hostile_cognition | 18.584 | 1.509 | 12.314 | 0 | 15.566 | 21.491 | ||
| neg_affect | neg_affect | 36.886 | 2.623 | 14.061 | 0 | 31.569 | 41.763 | ||
| aggression | aggression | 135.521 | 9.992 | 13.563 | 0 | 114.631 | 153.647 | ||
| provocation | provocation | 35.935 | 0.000 | NA | NA | 35.935 | 35.935 | ||
| ind1 | := | a1*b1 | ind1 | 0.309 | 0.078 | 3.949 | 0 | 0.163 | 0.471 |
| ind2 | := | a2*b2 | ind2 | 0.310 | 0.073 | 4.263 | 0 | 0.176 | 0.458 |
| total_ind | := | ind1+ind2 | total_ind | 0.619 | 0.095 | 6.500 | 0 | 0.442 | 0.811 |
total_ind is the sum of the two specific indirect
effects, with its own percentile bootstrap CI in the table above. The
~~ row for hostile_cognition ~~ neg_affect is
the residual covariance between the two mediators; its standardized
value in the summary is the residual correlation.
6 Model 3: does hostile cognition itself intensify negative affect?
Model 3 replaces Model 2’s residual covariance between the mediators
with a directional path, neg_affect ~ hostile_cognition
(\(d_1\)). Model 2’s indirect effect
through negative affect splits exactly into Model 3’s ind2
(direct route to negative affect) and ind3 (the serial
route through hostile cognition first).
m3 <- '
hostile_cognition ~ a1*provocation
neg_affect ~ a2*provocation + d1*hostile_cognition
aggression ~ cp*provocation + b1*hostile_cognition + b2*neg_affect
ind1 := a1*b1
ind2 := a2*b2
ind3 := a1*d1*b2
total_ind := ind1 + ind2 + ind3
'
fit3 <- sem(m3, data = d, se = "bootstrap", bootstrap = 5000, iseed = 4171)
summary(fit3, standardized = TRUE, fit.measures = TRUE, rsquare = TRUE)## lavaan 0.6-21 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 9
##
## Number of observations 300
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Model Test Baseline Model:
##
## Test statistic 327.093
## Degrees of freedom 6
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.000
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -2977.686
## Loglikelihood unrestricted model (H1) -2977.686
##
## Akaike (AIC) 5973.372
## Bayesian (BIC) 6006.706
## Sample-size adjusted Bayesian (SABIC) 5978.163
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.000
## P-value H_0: RMSEA <= 0.050 NA
## P-value H_0: RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.000
##
## Parameter Estimates:
##
## Standard errors Bootstrap
## Number of requested bootstrap draws 5000
## Number of successful bootstrap draws 4962
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## hostile_cognition ~
## provocatn (a1) 0.416 0.043 9.570 0.000 0.416 0.501
## neg_affect ~
## provocatn (a2) 0.304 0.061 5.012 0.000 0.304 0.271
## hstl_cgnt (d1) 0.438 0.077 5.684 0.000 0.438 0.324
## aggression ~
## provocatn (cp) 0.599 0.129 4.643 0.000 0.599 0.241
## hstl_cgnt (b1) 0.743 0.161 4.608 0.000 0.743 0.248
## neg_affct (b2) 0.638 0.122 5.209 0.000 0.638 0.289
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .hostile_cogntn 18.584 1.511 12.303 0.000 18.584 0.749
## .neg_affect 33.321 2.427 13.730 0.000 33.321 0.734
## .aggression 135.521 10.005 13.545 0.000 135.521 0.611
##
## R-Square:
## Estimate
## hostile_cogntn 0.251
## neg_affect 0.266
## aggression 0.389
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## ind1 0.309 0.078 3.938 0.000 0.309 0.124
## ind2 0.194 0.055 3.504 0.000 0.194 0.078
## ind3 0.116 0.034 3.431 0.001 0.116 0.047
## total_ind 0.619 0.095 6.487 0.000 0.619 0.249parameterEstimates(fit3, boot.ci.type = "perc", level = .95) %>%
filter(op == ":=") %>% kable(digits = 3)| lhs | op | rhs | label | est | se | z | pvalue | ci.lower | ci.upper |
|---|---|---|---|---|---|---|---|---|---|
| ind1 | := | a1*b1 | ind1 | 0.309 | 0.078 | 3.938 | 0.000 | 0.163 | 0.471 |
| ind2 | := | a2*b2 | ind2 | 0.194 | 0.055 | 3.504 | 0.000 | 0.095 | 0.311 |
| ind3 | := | a1d1b2 | ind3 | 0.116 | 0.034 | 3.431 | 0.001 | 0.057 | 0.190 |
| total_ind | := | ind1+ind2+ind3 | total_ind | 0.619 | 0.095 | 6.487 | 0.000 | 0.441 | 0.812 |
# Model 2's ind2 = (Model 3's ind2) + (Model 3's ind3), exactly
coef(fit2)[["a2"]] * coef(fit2)[["b2"]]## [1] 0.3103coef(fit3)[["a2"]] * coef(fit3)[["b2"]] + coef(fit3)[["a1"]] * coef(fit3)[["d1"]] * coef(fit3)[["b2"]]## [1] 0.31037 Model 4: does trait irritability matter too?
m4 <- '
hostile_cognition ~ a1*provocation + a2*trait_irritability
aggression ~ c1*provocation + c2*trait_irritability + b1*hostile_cognition
provocation ~~ trait_irritability
ind_x1 := a1*b1
ind_x2 := a2*b1
'
fit4 <- sem(m4, data = d, se = "bootstrap", bootstrap = 5000, iseed = 4171)
summary(fit4, standardized = TRUE, fit.measures = TRUE, rsquare = TRUE)## lavaan 0.6-21 ended normally after 23 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 10
##
## Number of observations 300
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Model Test Baseline Model:
##
## Test statistic 253.402
## Degrees of freedom 6
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.000
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -3988.310
## Loglikelihood unrestricted model (H1) -3988.310
##
## Akaike (AIC) 7996.621
## Bayesian (BIC) 8033.659
## Sample-size adjusted Bayesian (SABIC) 8001.945
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.000
## P-value H_0: RMSEA <= 0.050 NA
## P-value H_0: RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.000
##
## Parameter Estimates:
##
## Standard errors Bootstrap
## Number of requested bootstrap draws 5000
## Number of successful bootstrap draws 4999
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## hostile_cognition ~
## provocatn (a1) 0.367 0.043 8.478 0.000 0.367 0.441
## trt_rrtbl (a2) 0.153 0.032 4.815 0.000 0.153 0.215
## aggression ~
## provocatn (c1) 0.748 0.125 5.991 0.000 0.748 0.301
## trt_rrtbl (c2) 0.270 0.105 2.562 0.010 0.270 0.126
## hstl_cgnt (b1) 0.922 0.156 5.906 0.000 0.922 0.308
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## provocation ~~
## trait_irrtblty 11.540 2.313 4.989 0.000 11.540 0.276
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .hostile_cogntn 17.522 1.475 11.879 0.000 17.522 0.707
## .aggression 145.993 10.885 13.412 0.000 145.993 0.658
## provocation 35.935 2.957 12.152 0.000 35.935 1.000
## trait_irrtblty 48.789 3.873 12.597 0.000 48.789 1.000
##
## R-Square:
## Estimate
## hostile_cogntn 0.293
## aggression 0.342
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## ind_x1 0.338 0.073 4.626 0.000 0.338 0.136
## ind_x2 0.141 0.037 3.797 0.000 0.141 0.066parameterEstimates(fit4, boot.ci.type = "perc", level = .95) %>%
filter(op %in% c(":=", "~~")) %>% kable(digits = 3)| lhs | op | rhs | label | est | se | z | pvalue | ci.lower | ci.upper |
|---|---|---|---|---|---|---|---|---|---|
| provocation | trait_irritability | 11.540 | 2.313 | 4.989 | 0 | 7.080 | 16.102 | ||
| hostile_cognition | hostile_cognition | 17.522 | 1.475 | 11.879 | 0 | 14.503 | 20.334 | ||
| aggression | aggression | 145.993 | 10.885 | 13.412 | 0 | 123.596 | 165.990 | ||
| provocation | provocation | 35.935 | 2.957 | 12.152 | 0 | 30.238 | 41.715 | ||
| trait_irritability | trait_irritability | 48.789 | 3.873 | 12.597 | 0 | 41.286 | 56.313 | ||
| ind_x1 | := | a1*b1 | ind_x1 | 0.338 | 0.073 | 4.626 | 0 | 0.204 | 0.490 |
| ind_x2 | := | a2*b1 | ind_x2 | 0.141 | 0.037 | 3.797 | 0 | 0.074 | 0.219 |
8 Model fit: cutoff conventions, applied to our own models
fitMeasures(fit4, c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr"))## chisq df pvalue cfi tli rmsea srmr
## 0 0 NA 1 1 0 0Model 4 as drawn is saturated (df = 0); the fixed-parameter and equality-constraint demonstrations below are what give these indices real work to do.
9 Model Identification: Just, Over and Under Identified models
9.1 Just-identified: Model 4 as drawn (df = 0, shown above)
9.2 Over-identified: fixing one direct effect, in the simplest model
m1_saturated <- '
hostile_cognition ~ a*provocation
aggression ~ cp*provocation + b*hostile_cognition
'
m1_over <- '
hostile_cognition ~ a*provocation
aggression ~ b*hostile_cognition # cp omitted: fixed to 0
'
fit1_sat <- sem(m1_saturated, data = d)
fit1_over <- sem(m1_over, data = d)
fitMeasures(fit1_sat, c("chisq","df","cfi","tli","rmsea"))## chisq df cfi tli rmsea
## 0 0 1 1 0fitMeasures(fit1_over, c("chisq","df","cfi","tli","rmsea"))## chisq df cfi tli rmsea
## 32.281 1.000 0.846 0.537 0.323anova(fit1_sat, fit1_over)##
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)
## fit1_sat 0 4091 4109 0.0
## fit1_over 1 4121 4136 32.3 32.3 0.323 1 1.3e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 19.3 Under-identified: reciprocal causation cannot be estimated
m_reciprocal <- '
aggression ~ b*hostile_cognition
hostile_cognition ~ bprime*aggression
'
fit_reciprocal <- sem(m_reciprocal, data = d)## 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.lavInspect(fit_reciprocal, "converged")## [1] TRUEsummary(fit_reciprocal)## lavaan 0.6-21 ended normally after 11 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 4
##
## Number of observations 300
##
## Model Test User Model:
##
## Test statistic NA
## Degrees of freedom -1
## P-value (Unknown) NA
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## aggression ~
## hstl_cg (b) 1.279 NA
## hostile_cognition ~
## aggrssn (bprm) 0.031 NA
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .aggression 167.230 NA
## .hostile_cogntn 22.679 NAThis chunk overrides the global warning=FALSE setting
(error=TRUE, warning=TRUE) so that lavaan’s
own warning is printed. With two observed variables there are \(p = 3\) and \(q =
4\) free parameters, so \(df_M =
-1\): lavaanreports no test statistic, and returns
NA for every standard error. The point estimates it prints
are one of infinitely many equally good solutions and carry no
meaning.
10 Fixed and constrained parameters
10.1 The complete-mediation model: fixing both direct effects to zero
m4_complete <- '
hostile_cognition ~ a1*provocation + a2*trait_irritability
aggression ~ b1*hostile_cognition # c1, c2 omitted = fixed to 0
provocation ~~ trait_irritability
'
fit4_complete <- sem(m4_complete, data = d)
anova(fit4, fit4_complete)##
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)
## fit4 0 7997 8034 0.0
## fit4_complete 2 8031 8061 38.6 38.6 0.247 2 4.2e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fitMeasures(fit4_complete, c("chisq","df","cfi","tli","rmsea","srmr"))## chisq df cfi tli rmsea srmr
## 38.568 2.000 0.852 0.557 0.247 0.08810.2 The equality-constrained model: constraining the two ‘a-paths’ to be equal
m4_equal <- '
hostile_cognition ~ a*provocation + a*trait_irritability
aggression ~ c1*provocation + c2*trait_irritability + b1*hostile_cognition
provocation ~~ trait_irritability
'
fit4_equal <- sem(m4_equal, data = d)
anova(fit4, fit4_equal)##
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)
## fit4 0 7997 8034 0.0
## fit4_equal 1 8006 8039 11.5 11.5 0.187 1 0.00071 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fitMeasures(fit4_equal, c("chisq","df","cfi","tli","rmsea","srmr"))## chisq df cfi tli rmsea srmr
## 11.458 1.000 0.958 0.746 0.187 0.047parameterEstimates(fit4_equal) %>% filter(op == "~") %>% kable(digits = 3)| lhs | op | rhs | label | est | se | z | pvalue | ci.lower | ci.upper |
|---|---|---|---|---|---|---|---|---|---|
| hostile_cognition | ~ | provocation | a | 0.247 | 0.024 | 10.424 | 0.000 | 0.201 | 0.294 |
| hostile_cognition | ~ | trait_irritability | a | 0.247 | 0.024 | 10.424 | 0.000 | 0.201 | 0.294 |
| aggression | ~ | provocation | c1 | 0.748 | 0.128 | 5.862 | 0.000 | 0.498 | 0.998 |
| aggression | ~ | trait_irritability | c2 | 0.270 | 0.111 | 2.419 | 0.016 | 0.051 | 0.488 |
| aggression | ~ | hostile_cognition | b1 | 0.922 | 0.164 | 5.638 | 0.000 | 0.601 | 1.242 |
11 Multi-group analysis: does the pathway differ by gender?
dmg <- read.csv("../source_material/data/gam_multigroup.csv")
dmg$gender <- factor(dmg$gender, levels = c(0,1), labels = c("Men","Women"))
table(dmg$gender)##
## Men Women
## 146 154lavaan numbers the groups in the order they first appear
in the data, so women are Group 1 and men Group 2 in the output
below.
mg_model_free leaves every path unlabeled, which keeps
every path free in each group under sem(..., group=). A
single label in a shared multi-group model string – as in
mg_model_constrained below – forces that path to be equal
across groups; lavaan prints a warning saying so.
11.1 The free model: every path differs by group
mg_model_free <- '
hostile_cognition ~ provocation
aggression ~ hostile_cognition + provocation
'
fit_free <- sem(mg_model_free, data = dmg, group = "gender")
summary(fit_free, standardized = TRUE, fit.measures = TRUE, rsquare = TRUE)## lavaan 0.6-21 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 14
##
## Number of observations per group:
## Women 154
## Men 146
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
## Test statistic for each group:
## Women 0.000
## Men 0.000
##
## Model Test Baseline Model:
##
## Test statistic 289.416
## Degrees of freedom 6
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.000
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1997.140
## Loglikelihood unrestricted model (H1) -1997.140
##
## Akaike (AIC) 4022.281
## Bayesian (BIC) 4074.134
## Sample-size adjusted Bayesian (SABIC) 4029.734
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.000
## P-value H_0: RMSEA <= 0.050 NA
## P-value H_0: RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.000
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
##
## Group 1 [Women]:
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## hostile_cognition ~
## provocation 0.493 0.058 8.536 0.000 0.493 0.567
## aggression ~
## hostile_cogntn 1.595 0.189 8.416 0.000 1.595 0.540
## provocation 0.781 0.165 4.735 0.000 0.781 0.304
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .hostile_cogntn 5.187 1.191 4.354 0.000 5.187 1.028
## .aggression -0.391 2.969 -0.132 0.895 -0.391 -0.026
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .hostile_cogntn 17.278 1.969 8.775 0.000 17.278 0.679
## .aggression 95.545 10.888 8.775 0.000 95.545 0.430
##
## R-Square:
## Estimate
## hostile_cogntn 0.321
## aggression 0.570
##
##
## Group 2 [Men]:
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## hostile_cognition ~
## provocation 0.247 0.063 3.914 0.000 0.247 0.308
## aggression ~
## hostile_cogntn 1.532 0.194 7.911 0.000 1.532 0.514
## provocation 0.700 0.155 4.508 0.000 0.700 0.293
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .hostile_cogntn 10.054 1.336 7.526 0.000 10.054 2.040
## .aggression 4.012 3.682 1.089 0.276 4.012 0.273
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .hostile_cogntn 21.981 2.573 8.544 0.000 21.981 0.905
## .aggression 120.349 14.086 8.544 0.000 120.349 0.557
##
## R-Square:
## Estimate
## hostile_cogntn 0.095
## aggression 0.44311.3 Testing equality of the ‘a-path’ across groups
anova(fit_free, fit_constrained)##
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)
## fit_free 0 4022 4074 0.00
## fit_constrained 1 4028 4077 8.13 8.13 0.218 1 0.0043 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1LL_free <- unname(fitMeasures(fit_free, "logl"))
LL_constrained <- unname(fitMeasures(fit_constrained, "logl"))
LR <- 2 * (LL_free - LL_constrained)
tibble(LL_free, LL_constrained, LR, p = pchisq(LR, df = 1, lower.tail = FALSE)) %>%
kable(digits = 4)| LL_free | LL_constrained | LR | p |
|---|---|---|---|
| -1997 | -2001 | 8.133 | 0.004347 |
# Constraining every regression path at once, as an alternative:
fit_constrained_all <- sem(mg_model_free, data = dmg, group = "gender",
group.equal = "regressions")
anova(fit_free, fit_constrained_all)##
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)
## fit_free 0 4022 4074 0.00
## fit_constrained_all 3 4025 4065 8.44 8.44 0.11 3 0.038 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 111.4 Per-group indirect effects, with bootstrap CIs
One label per group (c(a_w, a_m)) names each group’s
path without constraining it, so each group’s indirect effect, and their
difference, can be defined with :=. Point estimates are
identical to the free model above.
mg_model_labeled <- '
hostile_cognition ~ c(a_w, a_m)*provocation
aggression ~ c(b_w, b_m)*hostile_cognition + c(cp_w, cp_m)*provocation
ind_women := a_w*b_w
ind_men := a_m*b_m
ind_diff := ind_women - ind_men
'
fit_labeled <- sem(mg_model_labeled, data = dmg, group = "gender",
se = "bootstrap", bootstrap = 5000, iseed = 4171)
parameterEstimates(fit_labeled, boot.ci.type = "perc", level = .95) %>%
filter(op == ":=") %>% kable(digits = 3)| lhs | op | rhs | block | group | label | est | se | z | pvalue | ci.lower | ci.upper |
|---|---|---|---|---|---|---|---|---|---|---|---|
| ind_women | := | a_w*b_w | 0 | 0 | ind_women | 0.786 | 0.125 | 6.270 | 0.000 | 0.554 | 1.046 |
| ind_men | := | a_m*b_m | 0 | 0 | ind_men | 0.379 | 0.093 | 4.094 | 0.000 | 0.206 | 0.571 |
| ind_diff | := | ind_women-ind_men | 0 | 0 | ind_diff | 0.408 | 0.156 | 2.616 | 0.009 | 0.103 | 0.722 |
11.5 Cross-check: the equivalent interaction-term test
dmg$gender_dummy <- as.integer(dmg$gender == "Women")
mod_interact <- lm(hostile_cognition ~ provocation * gender_dummy, data = dmg)
summary(mod_interact)##
## Call:
## lm(formula = hostile_cognition ~ provocation * gender_dummy,
## data = dmg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.875 -2.993 0.187 2.989 10.105
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.0539 1.2688 7.92 4.7e-14 ***
## provocation 0.2473 0.0600 4.12 4.9e-05 ***
## gender_dummy -4.8669 1.7997 -2.70 0.0072 **
## provocation:gender_dummy 0.2458 0.0862 2.85 0.0047 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.45 on 296 degrees of freedom
## Multiple R-squared: 0.214, Adjusted R-squared: 0.206
## F-statistic: 26.8 on 3 and 296 DF, p-value: 2.21e-15The provocation:gender_dummy coefficient is the
difference between the women’s and the men’s \(a\) paths (\(0.493 - 0.247\)), tested in a single pooled
regression; it leads to the same conclusion as the multi-group \(\chi^2\) difference test.
12 The
manymome package
# Refit Model 3 (serial mediation) as an ordinary, unlabeled model --
# manymome does not require path labels in advance.
fit_m3_plain <- sem('
hostile_cognition ~ provocation
neg_affect ~ provocation + hostile_cognition
aggression ~ provocation + hostile_cognition + neg_affect
', data = d)12.1 Worked example 1: indirect effects along the serial mediation model
out1 <- indirect_effect(x = "provocation", y = "aggression",
m = "hostile_cognition", fit = fit_m3_plain,
boot_ci = TRUE, R = 5000, seed = 4171)
out1##
## == Indirect Effect ==
##
## Path: provocation -> hostile_cognition -> aggression
## Indirect Effect: 0.309
## 95.0% Bootstrap CI: [0.163 to 0.471]
##
## Computation Formula:
## (b.hostile_cognition~provocation)*(b.aggression~hostile_cognition)
##
## Computation:
## (0.41584)*(0.74286)
##
##
## Percentile confidence interval formed by nonparametric bootstrapping with 5000 bootstrap samples.
##
## Coefficients of Component Paths:
## Path Coefficient
## hostile_cognition~provocation 0.416
## aggression~hostile_cognition 0.743out2 <- indirect_effect(x = "provocation", y = "aggression",
m = c("hostile_cognition", "neg_affect"), fit = fit_m3_plain,
boot_ci = TRUE, R = 5000, seed = 4171)
out2##
## == Indirect Effect ==
##
## Path: provocation -> hostile_cognition -> neg_affect -> aggression
## Indirect Effect: 0.116
## 95.0% Bootstrap CI: [0.057 to 0.190]
##
## Computation Formula:
## (b.hostile_cognition~provocation)*(b.neg_affect~hostile_cognition)*(b.aggression~neg_affect)
##
## Computation:
## (0.41584)*(0.43801)*(0.63801)
##
##
## Percentile confidence interval formed by nonparametric bootstrapping with 5000 bootstrap samples.
##
## Coefficients of Component Paths:
## Path Coefficient
## hostile_cognition~provocation 0.416
## neg_affect~hostile_cognition 0.438
## aggression~neg_affect 0.63812.2 Worked example 2: conditional indirect effects, a continuous moderator
d_c <- d %>% mutate(
c_provocation = provocation - mean(provocation),
c_trait_irr = trait_irritability - mean(trait_irritability),
prov_x_irr = c_provocation * c_trait_irr
)
fit_mome <- sem('
hostile_cognition ~ c_provocation + c_trait_irr + prov_x_irr
aggression ~ hostile_cognition + c_provocation
', data = d_c)
parameterEstimates(fit_mome) %>% filter(op == "~") %>% kable(digits = 4)| lhs | op | rhs | est | se | z | pvalue | ci.lower | ci.upper |
|---|---|---|---|---|---|---|---|---|
| hostile_cognition | ~ | c_provocation | 0.3664 | 0.0419 | 8.7486 | 0.0000 | 0.2843 | 0.4485 |
| hostile_cognition | ~ | c_trait_irr | 0.1535 | 0.0359 | 4.2714 | 0.0000 | 0.0831 | 0.2240 |
| hostile_cognition | ~ | prov_x_irr | -0.0056 | 0.0061 | -0.9102 | 0.3627 | -0.0176 | 0.0064 |
| aggression | ~ | hostile_cognition | 1.0223 | 0.1635 | 6.2517 | 0.0000 | 0.7018 | 1.3428 |
| aggression | ~ | c_provocation | 0.7930 | 0.1358 | 5.8379 | 0.0000 | 0.5268 | 1.0593 |
w_levels <- mod_levels("c_trait_irr", fit = fit_mome)
out3 <- cond_indirect_effects(wlevels = w_levels, x = "c_provocation", y = "aggression",
m = "hostile_cognition", fit = fit_mome,
boot_ci = TRUE, R = 5000, seed = 4171)
out3##
## == Conditional indirect effects ==
##
## Path: c_provocation -> hostile_cognition -> aggression
## Conditional on moderator(s): c_trait_irr
## Moderator(s) represented by: c_trait_irr
##
## [c_trait_irr] (c_trait_irr) ind CI.lo CI.hi Sig hostile_cognition~c_provocation aggression~hostile_cognition
## 1 M+1.0SD 6.997 0.335 0.177 0.511 Sig 0.327 1.022
## 2 Mean -0.000 0.375 0.239 0.529 Sig 0.366 1.022
## 3 M-1.0SD -6.997 0.414 0.259 0.594 Sig 0.405 1.022
##
## - [CI.lo to CI.hi] are 95.0% percentile confidence intervals by nonparametric bootstrapping with 5000
## samples.
## - The 'ind' column shows the conditional indirect effects.
## - 'hostile_cognition~c_provocation','aggression~hostile_cognition' is/are the path coefficient(s) along the
## path conditional on the moderator(s).index_of_mome(x = "c_provocation", y = "aggression", m = "hostile_cognition",
w = "c_trait_irr", fit = fit_mome, boot_ci = TRUE, R = 5000, seed = 4171)##
## == Conditional indirect effects ==
##
## Path: c_provocation -> hostile_cognition -> aggression
## Conditional on moderator(s): c_trait_irr
## Moderator(s) represented by: c_trait_irr
##
## [c_trait_irr] (c_trait_irr) ind CI.lo CI.hi Sig hostile_cognition~c_provocation aggression~hostile_cognition
## 1 1 1 0.369 0.232 0.522 Sig 0.361 1.022
## 2 0 0 0.375 0.239 0.529 Sig 0.366 1.022
##
## == Index of Moderated Mediation ==
##
## Levels compared: Row 1 - Row 2
##
## x y Index CI.lo CI.hi
## Index c_provocation aggression -0.006 -0.018 0.006
##
## - [CI.lo, CI.hi]: 95% percentile confidence interval.13 Post-estimation: drawing the fitted model
semPaths(fit4, what = "std", layout = "tree", edge.label.cex = 0.9,
sizeMan = 10, edge.color = ntnublue)
lavaanPlot(model = fit4, coefs = TRUE, stand = TRUE, sig = .05)14 Session info
sessionInfo()## R version 4.4.2 (2024-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS 27.0.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Europe/Oslo
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] kableExtra_1.4.0 knitr_1.51 manymome_0.3.3 lavaanPlot_0.8.1 semPlot_1.1.8 lavaan_0.6-21
## [7] lubridate_1.9.5 forcats_1.0.1 stringr_1.6.0 dplyr_1.2.0 purrr_1.2.1 readr_2.2.0
## [13] tidyr_1.3.2 tibble_3.3.1 ggplot2_4.0.2 tidyverse_2.0.0 pacman_0.5.1
##
## loaded via a namespace (and not attached):
## [1] Rdpack_2.6.6 mnormt_2.1.2 pbapply_1.7-4 gridExtra_2.3 fdrtool_1.2.18
## [6] rlang_1.1.7 magrittr_2.0.4 otel_0.2.0 rockchalk_1.8.157 compiler_4.4.2
## [11] systemfonts_1.3.2 png_0.1-9 vctrs_0.7.2 reshape2_1.4.5 OpenMx_2.21.13
## [16] quadprog_1.5-8 pkgconfig_2.0.3 fastmap_1.2.0 arm_1.15-2 backports_1.5.0
## [21] pbivnorm_0.6.0 rmarkdown_2.30 tzdb_0.5.0 nloptr_2.2.1 xfun_0.57
## [26] cachem_1.1.0 kutils_1.73 jsonlite_2.0.0 jpeg_0.1-11 psych_2.6.1
## [31] parallel_4.4.2 cluster_2.1.8.2 R6_2.6.1 bslib_0.10.0 stringi_1.8.7
## [36] RColorBrewer_1.1-3 boot_1.3-32 rpart_4.1.24 jquerylib_0.1.4 Rcpp_1.1.1
## [41] base64enc_0.1-6 Matrix_1.7-5 splines_4.4.2 nnet_7.3-20 igraph_2.2.2
## [46] timechange_0.4.0 tidyselect_1.2.1 rstudioapi_0.18.0 abind_1.4-8 yaml_2.3.12
## [51] qgraph_1.9.8 lattice_0.22-9 plyr_1.8.9 withr_3.0.2 S7_0.2.1
## [56] coda_0.19-4.1 evaluate_1.0.5 foreign_0.8-91 RcppParallel_5.1.11-2 xml2_1.5.2
## [61] zip_2.3.3 pillar_1.11.1 DiagrammeR_1.0.11 carData_3.0-6 checkmate_2.3.4
## [66] stats4_4.4.2 reformulas_0.4.4 generics_0.1.4 hms_1.1.4 scales_1.4.0
## [71] minqa_1.2.8 gtools_3.9.5 xtable_1.8-8 glue_1.8.0 mi_1.2
## [76] Hmisc_5.2-5 tools_4.4.2 data.table_1.18.2.1 lme4_2.0-1 openxlsx_4.2.8.1
## [81] visNetwork_2.1.4 XML_3.99-0.23 grid_4.4.2 sem_3.1-16 rbibutils_2.4.1
## [86] colorspace_2.1-2 nlme_3.1-168 htmlTable_2.4.3 Formula_1.2-5 cli_3.6.5
## [91] textshaping_1.0.5 viridisLite_0.4.3 svglite_2.2.2 corpcor_1.6.10 glasso_1.11
## [96] gtable_0.3.6 sass_0.4.10 digest_0.6.39 htmlwidgets_1.6.4 farver_2.1.2
## [101] htmltools_0.5.9 lifecycle_1.0.5 lisrelToR_0.3 MASS_7.3-65