We will use data from the National Evaluation of Welfare-to-Work Strategies (NEWWS) study. The study randomly assigned 694 participants, who were mostly low-income single mothers with young children. Labor Force Attachment (LFA) program (208) aimed at moving low-income parents from welfare to work by providing employment-focused incentives and services. Control group (486) received aid from the Aid to Families with Dependent Children (AFDC) program without requirement for working.

Import the data:

library(foreign)
newws = read.dta("E:\\R\\newws.dta")
# Factorize the categorical variable that is coded with numerical values. Otherwise, R would take it as a continuous variable.

Research Question: How is the LFA impact on maternal depression transmitted through employment?

Treatment:

treat 1 if a participant was assigned to LFA and 0 otherwise.

Outcome:

depression maternal depression at the end of the second year after randomization, which is a summary score of 12 items measuring depressive symptoms during the past week on a 0-3 scale.

Mediator: emp 1 if one was ever employed during the two-year period after randomization.

Pretreatment Covariates:

nevmar 1 if never married and 0 otherwise.

emp_prior 1 if employed at baseline and 0 otherwise.

hispanic 1 if Hispanic and 0 otherwise.

ADCPC welfare amount in the year before randomization

attitude a composite score of two attitude items - so many family problems that I cannot work at a full time or part time job; so much to do during the day that I cannot go to a school or job training program - measured on the scale of 1-4.

depress_prior a composite score of three depressive symptom items - sad, depressed, blues, and lonely - in the week before randomization measured on the scale of 1-4.

workpref one’s level of preference for taking care of family full time than working on the scale of 1-4.

nohsdip 1 if had never obtained a high school diploma or a General Educational Development certificate and 0 otherwise.

1. Definitions of causal effects

NIE (natural indirect effect of LFA on depression mediated by employment): \[NIE = E[Y(1,M(1))] - E[Y(1,M(0))]\]

The average impact of LFA on maternal depression that is solely attributable to the change in their employment status induced by LFA, when the treatment is fixed at the LFA condition.

NDE (natural direct effect of LFA on depression transmitted through all the other possible pathways): \[NDE = E[Y(1,M(0))] - E[Y(0,M(0))]\]

The average impact of LFA on maternal depression while employment status was kept at the level that would have been observed in the absence of LFA.

PIE (pure indirect effect of LFA on depression transmitted through job search self-efficacy): \[PIE = E[Y(0,M(1))] - E[Y(0,M(0))]\]

The average impact of LFA on maternal depression that is solely attributable to the change in their employment status induced by LFA, when the treatment is fixed at the control condition.

INT (natural treatment-by-mediator interaction effect): \[INT = TDE - NDE = NIE - PIE = {E[Y(1,M(1))] - E[Y(1,M(0))]} - {E[Y(0,M(1))] - E[Y(0,M(0))]}\]

The average difference in how the LFA-induced change in the employment affects depression under the LFA condition and that under the control condition.

2. Regression-based mediation analysis (Imai et al., 2010)

Fit the mediator and outcome models:

l.m = glm(emp ~ treat + emp_prior + nevmar + hispanic + ADCPC + attitude + depress_prior + workpref + nohsdip, data = newws, family = binomial("probit")) 
summary(l.m)
## 
## Call:
## glm(formula = emp ~ treat + emp_prior + nevmar + hispanic + ADCPC + 
##     attitude + depress_prior + workpref + nohsdip, family = binomial("probit"), 
##     data = newws)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1487  -0.9178  -0.5949   0.9980   2.0360  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -9.671e-01  2.788e-01  -3.469 0.000523 ***
## treat          7.471e-01  1.133e-01   6.593  4.3e-11 ***
## emp_prior      7.840e-01  1.172e-01   6.692  2.2e-11 ***
## nevmar        -2.786e-02  1.037e-01  -0.269 0.788228    
## hispanic       8.800e-02  1.148e-01   0.766 0.443549    
## ADCPC         -3.029e-05  1.403e-05  -2.159 0.030826 *  
## attitude       7.423e-02  8.800e-02   0.844 0.398909    
## depress_prior  2.545e-02  6.163e-02   0.413 0.679624    
## workpref       1.551e-01  7.436e-02   2.085 0.037057 *  
## nohsdip       -2.457e-01  1.079e-01  -2.278 0.022752 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 960.01  on 693  degrees of freedom
## Residual deviance: 824.26  on 684  degrees of freedom
## AIC: 844.26
## 
## Number of Fisher Scoring iterations: 4
l.y = lm(depression ~ treat * emp + emp_prior + nevmar + hispanic + ADCPC + attitude + depress_prior + workpref + nohsdip, data = newws) 
summary(l.y)
## 
## Call:
## lm(formula = depression ~ treat * emp + emp_prior + nevmar + 
##     hispanic + ADCPC + attitude + depress_prior + workpref + 
##     nohsdip, data = newws)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.569  -5.282  -1.970   3.826  24.706 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    4.803e+00  1.568e+00   3.064  0.00227 ** 
## treat          2.538e+00  9.940e-01   2.553  0.01088 *  
## emp            4.412e-01  7.368e-01   0.599  0.54948    
## emp_prior      5.929e-01  6.826e-01   0.869  0.38536    
## nevmar         9.720e-01  5.865e-01   1.657  0.09793 .  
## hispanic       1.508e+00  6.462e-01   2.334  0.01987 *  
## ADCPC         -5.951e-05  7.998e-05  -0.744  0.45707    
## attitude      -4.740e-01  4.941e-01  -0.959  0.33771    
## depress_prior  1.761e+00  3.476e-01   5.065 5.25e-07 ***
## workpref      -2.775e-01  4.180e-01  -0.664  0.50709    
## nohsdip        4.134e-01  6.146e-01   0.673  0.50145    
## treat:emp     -3.332e+00  1.306e+00  -2.551  0.01094 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.523 on 682 degrees of freedom
## Multiple R-squared:  0.07147,    Adjusted R-squared:  0.05649 
## F-statistic: 4.772 on 11 and 682 DF,  p-value: 4.303e-07

Run the regression-based mediation analysis:

library(mediation)
set.seed(1) # So that the results will stay the same when you run mediate() again.
med = mediate(l.m, l.y, treat = "treat", mediator = "emp")

Print and visualize the mediation analysis results：

summary(med)
## 
## Causal Mediation Analysis 
## 
## Quasi-Bayesian Confidence Intervals
## 
##                          Estimate 95% CI Lower 95% CI Upper p-value   
## ACME (control)             0.1073      -0.2371         0.51   0.588   
## ACME (treated)            -0.7544      -1.3927        -0.19   0.010 **
## ADE (control)              1.2200      -0.0601         2.59   0.064 . 
## ADE (treated)              0.3583      -0.8741         1.67   0.576   
## Total Effect               0.4656      -0.7349         1.69   0.456   
## Prop. Mediated (control)   0.0853      -3.6313         2.40   0.800   
## Prop. Mediated (treated)  -0.8247     -11.2715        16.73   0.462   
## ACME (average)            -0.3236      -0.6810         0.02   0.068 . 
## ADE (average)              0.7892      -0.4146         2.07   0.200   
## Prop. Mediated (average)  -0.3697      -4.3403         7.97   0.512   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Sample Size Used: 694 
## 
## 
## Simulations: 1000
plot(med)

Estimate and test the natural treatment-by-mediator interaction effect:

INT = med$d1.sims - med$d0.sims 
mean(INT)
## [1] -0.8616139
quantile(INT, probs = c(0.025, 0.975))
##       2.5%      97.5% 
## -1.6261933 -0.2029703

Summarize the findings:

The mediator model output indicates that LFA significantly increased employment rates (p < 0.01). The odds of being employed among those assigned to the LFA group is 2.11 (exp(0.7471)) times that among those from the control group. The mediation analysis output indicates that the total LFA effect is estimated to be 0.47 (95% CI = [-0.73, 1.69]), which can be decomposed into a natural direct effect, estimated to be 1.22 (95% CI = [-0.06, 2.59]), about 15.92% of a standard deviation of the outcome in the control group, and a natural indirect effect, estimated to be -0.75 (95% CI = [-1.39, -0.19]), about -9.78% of a standard deviation of the outcome in the control group. The natural direct effect indicates that, LFA would have increased maternal depression if one’s employment experience is held at the level under the control condition, but not by a statistically significant amount at the significance level of 0.05. In contrast, the natural indirect effect reflects that the LFA-induced increase in employment rate significantly reduced one’s maternal depression under the LFA condition. The counteracting indirect and direct effects explained the insignificant total effect.

The natural indirect effect can be further decomposed into a pure indirect effect, which is estimated to be 0.11 (95% CI = [-0.24, 0.51]), and a natural treatment-by-mediator interaction effect, which is estimated to be -0.86 (95% CI = [-1.63, -0.20]). It indicates that the LFA-induced increase in employment rate reduced maternal depression more under the LFA condition (i.e., natural indirect effect) than under the control condition (i.e., pure indirect effect). Equivalently, the natural treatment-by-mediator interaction effect can also be viewed as the difference between the total direct effect and the natural direct effect, which indicates that LFA would have increased maternal depression to a smaller extent if holding one’s employment experience under the LFA condition (i.e., total direct effect) than if holding that under the control condition (i.e., natural direct effect).

The above analysis is based on the sequential ignorability assumption:

  1. (Strong ignorability of the treatment) No unmeasured confounder of the LFA-employment or LFA-depression relationship, i.e., the treatment is as if randomized within levels of pretreatment covariates.

  2. (Strong ignorability of the mediator) No unmeasured pretreatment confounder and no posttreatment confounder of the employment-depression relationship, i.e., the mediator is as if randomized within the same treatment group or across treatment groups among individuals with the same levels of pretreatment covariates.

  3. (Positivity assumption) Within levels of pretreatment covariates, every individual has a nonzero probability of being assigned to either the LFA group or the control group and a nonzero probability of taking each mediator value within the response space of the mediator under each treatment condition.

A sensitivity analysis is necessary for determining whether the above analysis results are sensitive to a potential violation of the sequential ignorability assumption. In the following sensitivity analysis strategy, the sensitivity parameter is rho, the correlation between the error term of the mediator model and that of the outcome model. Such a correlation can arise if there exist omitted confounders that affect both mediator and outcome because these omitted confounders will be part of the two error terms. Thus, under sequential ignorability, rho=0, and nonzero values of rho imply departures from the ignorability assumption.

With this fact, the proposed sensitivity analysis asks the question of how large rho has to be for the causal mediation effect to go away. If small departures from zero in rho produce an average causal mediation effect that is substantively different from the estimate obtained under sequential ignorability or reverse the significance of the causal effect under sequential ignorability, this suggests that the study is sensitive to the potential violation of the sequential ignorability assumption. The sensitivity analysis plot graphically illustrates this by plotting the estimated effects and their 95% confidence intervals as a function of rho.

Sensitivity analysis:

sens.indirect = medsens(med, effect.type = "indirect")
plot(sens.indirect)

sens.direct = medsens(med, effect.type = "direct")
plot(sens.direct)

The sensitivity analysis results indicate that NIE (ACME_1) turns to be zero or positive when rho is equal to or smaller than -0.2, PIE (ACME_0) turns to be zero or negative when rho is equal to or bigger than 0.05, NDE (ADE_0) turns to be 0 or negative when rho is equal to or smaller than -0.4, and TDE (ADE_1) turns to be 0 or negative when rho is equal to or smaller than -0.1.

We also find that NIE turns to be insignificant when rho is between -0.4 and -0.05, PIE turns to be significant when rho is smaller than -0.05 or bigger than 0.1, NDE turns to be significant when rho is bigger than 0.1, and TDE turns to be significant when rho is bigger than 0.5. Here, the significance level is 0.05.

Because it is hard to determine whether the value of rho for removing the effects, changing the sign of the effects, or changing their significance is likely to exist or not, it is challenging to conclude whether the analytic results are sensitive based on rho without comparisons to other studies. If we choose 0.1 as the cutoff value, then the sign of PIE and the significance of NIE and PIE are sensitive to potential unmeasured pretreatment confounding.

3. Weighting-based mediation analysis (Hong, 2015)

Fit propensity score models:

# Fit a propensity score model for the mediator in the treatment group 
data1 = newws[newws$treat == 1, ] # extract treatment group from the original data
l1 = glm(emp ~ treat + emp_prior + nevmar + hispanic + ADCPC + attitude + depress_prior + workpref + nohsdip, data = data1, family = binomial("logit"))
data1$logit.p = predict(l1, type = "link") # save the logit propensity score

# Fit a propensity score model for the mediator in the control group
data0 = newws[newws$treat == 0, ] # extract control group from the original data
l0 = glm(emp ~ treat + emp_prior + nevmar + hispanic + ADCPC + attitude + depress_prior + workpref + nohsdip, data = data0, family = binomial("logit"))
data0$logit.p = predict(l0, type = "link") # save the logit propensity score

# Combine the two separate data with the logit propensity score
newws = rbind(data1, data0)

Note: When applying the above code, you only need to replace newws with your data name and respecify l1 and l0.

Check overlapping in logit of propensity score:

library(ggplot2)
# Overlapping in distributions of logit propensity score between M = 1 and M = 0 in the treatment group
overlap1 = ggplot(data1, aes(x = logit.p, fill = as.factor(emp))) +       
    geom_density(alpha = .5) +     
    labs(x = "Logit Propensity Scores", y = "Density", fill = "") +     
    theme_bw() +  
    ggtitle("Overlapping in logit propensity scores between M = 1 and M = 0 in the treatment group")

# Overlapping in distributions of logit propensity score between M = 1 and M = 0 in the control group
overlap0 = ggplot(data0, aes(x = logit.p, fill = as.factor(emp))) +     
    geom_density(alpha = .5) +     
    labs(x = "Logit Propensity Scores", y = "Density", fill = "") +     
    theme_bw() +  
    ggtitle("Overlapping in logit propensity scores between M = 1 and M = 0 in the control group")

# Display the two plots side by side
library(cowplot)
plot_grid(overlap1, overlap0)

Note: When applying the above code, you only need to replace “emp” with your mediator name.

The results indicate that the distributions of the logit of propensity score well overlap between those who were employed and those who were unemployed within each treatment group.

Check balance:

If you have categorical variables with more than two categories, please generate dummies for them first.

mediator = "emp" # Let the program know the name of the mediator (string).

# Treatment group
data1$p = predict(l1, type = "response") # save the propensity score
# Calculate IPTW for those whose M = 1 in the treatment group
data1$iptw[data1[, mediator] == 1] = mean(data1[, mediator] == 1)/data1$p[data1[, mediator] == 1]
# Calculate IPTW for those whose M = 0 in the treatment group
data1$iptw[data1[, mediator] == 0] = mean(data1[, mediator] == 0)/ (1 - data1$p[data1[, mediator] == 0])

# Control group
data0$p = predict(l0, type = "response") # save the propensity score
# Calculate IPTW for those whose M = 1 in the control group
data0$iptw[data0[, mediator] == 1] = mean(data0[, mediator] == 1)/data0$p[data0[, mediator] == 1]
# Calculate IPTW for those whose M = 0 in the control group
data0$iptw[data0[, mediator] == 0] = mean(data0[, mediator] == 0)/ (1 - data0$p[data0[, mediator] == 0])

Note: When applying the above code, you only need to change the mediator name in the first line of code.

balance = function(x, data, treatment, mediator){
  sd.mean.diff.before = NULL
  sd.mean.diff = NULL
  for(i in 1:length(x)){
    # balance between M = 1 and M = 0 before weighting
    mean.1.before = mean(data[data[, mediator] == 1, x[i]])
    mean.0.before = mean(data[data[, mediator] == 0, x[i]])
    mean.diff.before = mean.1.before - mean.0.before
    sd.mean.diff.before = c(sd.mean.diff.before, mean.diff.before/sd(data[,x[i]]))
    # balance between M = 1 and M = 0 after weighting
    mean.1 = sum(data[,x[i]] * data$iptw * data[, mediator])/sum(data$iptw * data[, mediator])  
    mean.0 = sum(data[,x[i]] * data$iptw * (1 - data[, mediator]))/sum(data$iptw * (1 - data[, mediator]))  
    mean.diff = mean.1 - mean.0  
    sd.mean.diff = c(sd.mean.diff, mean.diff/sd(data[,x[i]]))  
  }
  results = cbind.data.frame(c(sd.mean.diff.before, sd.mean.diff),c(rep("before weighting", length(x)), rep("after weighting", length(x))), c(x, x))
  colnames(results) = c("sd.mean.diff", "weighting", "covariates")
  return(results)
}
# Apply the balance checking function
balance.results.1 = balance(x = c("logit.p", "emp_prior", "nevmar", "hispanic", "ADCPC", "attitude", "depress_prior", "workpref", "nohsdip"), data = data1, treatment = "treat", mediator = "emp")
balance.results.0 = balance(x = c("logit.p", "emp_prior", "nevmar", "hispanic", "ADCPC", "attitude", "depress_prior", "workpref", "nohsdip"), data = data0, treatment = "treat", mediator = "emp")

balance.results.1
##    sd.mean.diff        weighting    covariates
## 1    0.80110520 before weighting       logit.p
## 2    0.55049118 before weighting     emp_prior
## 3   -0.09819591 before weighting        nevmar
## 4    0.03852494 before weighting      hispanic
## 5   -0.22618511 before weighting         ADCPC
## 6    0.30728197 before weighting      attitude
## 7    0.06467071 before weighting depress_prior
## 8    0.51057152 before weighting      workpref
## 9   -0.41283143 before weighting       nohsdip
## 10   0.08459082  after weighting       logit.p
## 11   0.03697812  after weighting     emp_prior
## 12  -0.02153553  after weighting        nevmar
## 13   0.04359517  after weighting      hispanic
## 14  -0.01906914  after weighting         ADCPC
## 15   0.05606602  after weighting      attitude
## 16  -0.01623545  after weighting depress_prior
## 17   0.08095841  after weighting      workpref
## 18  -0.01681181  after weighting       nohsdip
balance.results.0
##    sd.mean.diff        weighting    covariates
## 1   0.776608319 before weighting       logit.p
## 2   0.695477481 before weighting     emp_prior
## 3  -0.019902311 before weighting        nevmar
## 4  -0.030022466 before weighting      hispanic
## 5  -0.438625182 before weighting         ADCPC
## 6   0.260893924 before weighting      attitude
## 7  -0.015375879 before weighting depress_prior
## 8   0.314749856 before weighting      workpref
## 9  -0.192418453 before weighting       nohsdip
## 10  0.012963059  after weighting       logit.p
## 11  0.006277949  after weighting     emp_prior
## 12  0.007300551  after weighting        nevmar
## 13 -0.022010343  after weighting      hispanic
## 14 -0.019370907  after weighting         ADCPC
## 15  0.012635576  after weighting      attitude
## 16 -0.002996157  after weighting depress_prior
## 17  0.010356045  after weighting      workpref
## 18  0.007203650  after weighting       nohsdip

Note: For the above code, you only need to change x, treatment, and mediator when applying the balance checking function.

Generate balance checking plot:

ggplot(balance.results.1, aes(sd.mean.diff, covariates)) +
  geom_point(aes(color = weighting)) +
  labs(x = "Standardized Mean Differences", y = "") +
  theme(legend.title=element_blank()) +
  geom_vline(xintercept=c(-0.25,-0.1,0.1,0.25), linetype="longdash", col = c("red", "blue", "blue", "red")) +
  xlim(-2.5, 2.5)

ggplot(balance.results.0, aes(sd.mean.diff, covariates)) +
  geom_point(aes(color = weighting)) +
  labs(x = "Standardized Mean Differences", y = "") +
  theme(legend.title=element_blank()) +
  geom_vline(xintercept=c(-0.25,-0.1,0.1,0.25), linetype="longdash", col = c("red", "blue", "blue", "red")) +
  xlim(-2.5, 2.5)

Note: When applying the above code, you do not need to make any change.

The balance checking results indicate that logit propensity score, workpref, nohsdip, emp_prior, and attitude are originally unbalanced between those who were employed and those who were unemployed within the treatment group, and logit propensity score, workpref, emp_prior, and ADCPC are originally unbalanced between those who were employed and those who were unemployed within the control group. After propensity score adjustment, the distributions of logit propensity score and all the covariates are well balanced between those who were employed and those who were unemployed within each treatment group.

Estimate causal effects:

If you have categorical variables with more than two categories, please generate dummies for them first. Please avoid signs such as “/” in your variables names using the code: gsub(“[[:punct:]]”, ““, x), where x is a vector of covariate names. Please only include g-1 dummies for a categorical variable with g categories.

library(rmpw)
rmpw(data = newws, treatment = "treat", mediator = "emp", outcome = "depression", propensity_x = c("emp_prior", "nevmar", "hispanic", "ADCPC", "attitude", "depress_prior", "workpref", "nohsdip"), outcome_x = c("emp_prior", "nevmar", "hispanic", "ADCPC", "attitude", "depress_prior", "workpref", "nohsdip"),decomposition = 1)
##                           Estimate Std.Error t value Pr(>|t|)   
## Gamma.0                     5.0089    1.7036  2.9401   0.0033  *
## Natural Direct Effect       1.2192    0.8494  1.4355   0.1511   
## Natural Indirect Effect    -0.7607    0.4429 -1.7176   0.0859   
## Pure Indirect Effect       -0.0217     0.213 -0.1017    0.919   
## T-by-M Interaction Effect   -0.739    0.5049 -1.4637   0.1433   
## emp_prior                   0.1219    0.6916  0.1762   0.8601   
## nevmar                      0.9775    0.6346  1.5404   0.1235   
## hispanic                     1.436    0.7072  2.0305   0.0423  .
## ADCPC                       -1e-04     1e-04 -1.0248   0.3055   
## attitude                   -0.4179     0.526 -0.7945   0.4269   
## depress_prior               1.7242    0.3647  4.7276   <0.001 **
## workpref                   -0.2202    0.4452 -0.4947   0.6208   
## nohsdip                     0.4304    0.6382  0.6744   0.5001

The effect estimates are mostly similar to those obtained from the regression-based mediation analysis. However, the effects here are mostly insignificant, while the regression-based method reported significant NIE, NDE (significant at the significance level of 0.1, but insignificant at the significance level of 0.05), and INT. This is mainly due to the lower estimation efficiency of the weighting-based method. Nevertheless, the weighting-based method is more robust to possible model misspecifications.