Mediator effects

library(mediation)

Fit a model to estimate the total effect. treat is the treatment and depress2 the outcome. The other variables are covariates.

a <- lm(depress2 ~ treat + econ_hard + sex + age, data=jobs)

Now the models for the mediation test (job_seek is the mediator):

b <- lm(job_seek ~ treat + econ_hard + sex + age, data=jobs)
c <- lm(depress2 ~ treat + job_seek + econ_hard + sex + age, data=jobs)

Estimate the CIs for the various effects by nonparametric bootstrap:

contcont.boot <- mediate(b, c, boot=TRUE, sims=1000, treat="treat", mediator="job_seek")

Running nonparametric bootstrap

summary(contcont.boot)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

               Estimate 95% CI Lower 95% CI Upper p-value
ACME            -0.0157      -0.0400         0.01    0.21
ADE             -0.0403      -0.1237         0.04    0.40
Total Effect    -0.0560      -0.1452         0.03    0.24
Prop. Mediated   0.2811      -2.0516         2.53    0.37

Sample Size Used: 899 


Simulations: 1000 

Average direct effect (ADE):

coef(c)["treat"]

     treat 
-0.0402647 

Average causal mediation effect (ACME):

coef(b)["treat"] * coef(c)["job_seek"]

      treat 
-0.01574465 

Total effect:

coef(a)["treat"]

      treat 
-0.05600935 

Or the ADE + ACME:

coef(c)["treat"] + coef(b)["treat"] * coef(c)["job_seek"]

      treat 
-0.05600935 

Proportion mediated is ACME / total effect:

coef(b)["treat"] * coef(c)["job_seek"] / coef(a)["treat"]

    treat 
0.2811075 

Summary of all model coefficients:

summary(a)

Call:
lm(formula = depress2 ~ treat + econ_hard + sex + age, data = jobs)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.0584 -0.4935 -0.1346  0.3738  2.8561 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.3268664  0.1098292  12.081  < 2e-16 ***
treat       -0.0560093  0.0451637  -1.240   0.2152    
econ_hard    0.1357868  0.0215838   6.291 4.93e-10 ***
sex          0.1086375  0.0427532   2.541   0.0112 *  
age         -0.0004517  0.0020347  -0.222   0.8244    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6369 on 894 degrees of freedom
Multiple R-squared:  0.04924,   Adjusted R-squared:  0.04498 
F-statistic: 11.57 on 4 and 894 DF,  p-value: 3.635e-09

summary(b)

Call:
lm(formula = job_seek ~ treat + econ_hard + sex + age, data = jobs)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.1658 -0.4092  0.0577  0.5819  1.1121 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.670585   0.125169  29.325   <2e-16 ***
treat        0.065615   0.051472   1.275   0.2027    
econ_hard    0.053162   0.024598   2.161   0.0309 *  
sex         -0.007637   0.048725  -0.157   0.8755    
age          0.004586   0.002319   1.978   0.0482 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7259 on 894 degrees of freedom
Multiple R-squared:  0.01147,   Adjusted R-squared:  0.007048 
F-statistic: 2.593 on 4 and 894 DF,  p-value: 0.03529

summary(c)

Call:
lm(formula = depress2 ~ treat + job_seek + econ_hard + sex + 
    age, data = jobs)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.5192 -0.4461 -0.1503  0.3597  2.7568 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.2076414  0.1480548  14.911  < 2e-16 ***
treat       -0.0402647  0.0435059  -0.926   0.3550    
job_seek    -0.2399550  0.0282433  -8.496  < 2e-16 ***
econ_hard    0.1485434  0.0208269   7.132 2.04e-12 ***
sex          0.1068049  0.0411471   2.596   0.0096 ** 
age          0.0006489  0.0019625   0.331   0.7410    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.613 on 893 degrees of freedom
Multiple R-squared:  0.1203,    Adjusted R-squared:  0.1154 
F-statistic: 24.43 on 5 and 893 DF,  p-value: < 2.2e-16

library(jtools)

And some pics:

plot_summs(a, c, model.names = c("a", "b"))

plot_summs(b)