##Data Import
library(readxl)
stress_withdraw <- read_excel("Downloads/stress_withdraw.xlsx")
Local_groups <- read_excel("Downloads/Local_groups.xlsx")##Package Installation
#install.packages("mediation")
library(mediation)## Loading required package: MASS## Loading required package: Matrix## Loading required package: mvtnorm## Loading required package: sandwich## mediation: Causal Mediation Analysis
## Version: 4.5.0###Assignment
##1). Conduct causal mediation analysis for empirical example #2 through #4.
#Empirical Example #2 (Continuous M and Y with other Controls):
Ex2_m <- lm(depress~stress, data = stress_withdraw);summary(Ex2_m) ##
## Call:
## lm(formula = depress ~ stress, data = stress_withdraw)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.0095 -0.4195 -0.1609 0.2498 4.0278
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.79936 0.14331 5.578 6.11e-08 ***
## stress 0.17288 0.02965 5.831 1.63e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6819 on 260 degrees of freedom
## Multiple R-squared: 0.1156, Adjusted R-squared: 0.1122
## F-statistic: 34 on 1 and 260 DF, p-value: 1.63e-08Ex2_o <- lm(withdraw~stress+depress, data = stress_withdraw);summary(Ex2_o) ##
## Call:
## lm(formula = withdraw ~ stress + depress, data = stress_withdraw)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.1716 -0.9472 -0.2249 0.8490 2.9049
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.44706 0.25201 5.742 2.61e-08 ***
## stress -0.07685 0.05239 -1.467 0.144
## depress 0.76913 0.10306 7.463 1.29e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.133 on 259 degrees of freedom
## Multiple R-squared: 0.1804, Adjusted R-squared: 0.174
## F-statistic: 28.49 on 2 and 259 DF, p-value: 6.528e-12Ex2_mo <- mediate(Ex2_m,Ex2_o,treat="stress",mediator="depress", boot=T, sims=10000,boot.ci.type="bca")## Running nonparametric bootstrapsummary(Ex2_mo)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the BCa Method
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME 0.1330 0.0732 0.20 0.0002 ***
## ADE -0.0768 -0.1800 0.04 0.1770
## Total Effect 0.0561 -0.0739 0.17 0.3774
## Prop. Mediated 2.3694 -272.7353 2.16 0.3772
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 262
##
##
## Simulations: 10000#Empirical Example #3 (Binary M and Y, with continuous X)
Ex3_m <- glm(poverty ~age+groups, data = Local_groups, family = "binomial");summary(Ex3_m) ##
## Call:
## glm(formula = poverty ~ age + groups, family = "binomial", data = Local_groups)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.5745 -0.4738 -0.4079 -0.3204 3.5941
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.522441 0.119461 -12.744 < 2e-16 ***
## age -0.010859 0.002253 -4.820 1.44e-06 ***
## groups -0.368701 0.037558 -9.817 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 5743.2 on 9909 degrees of freedom
## Residual deviance: 5592.4 on 9907 degrees of freedom
## AIC: 5598.4
##
## Number of Fisher Scoring iterations: 6Ex3_o <- lm(SRH~poverty+groups+age, data = Local_groups);summary(Ex3_o) ##
## Call:
## lm(formula = SRH ~ poverty + groups + age, data = Local_groups)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.10875 0.06167 0.15693 0.22426 0.67949
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.9837700 0.0137680 71.45 <2e-16 ***
## poverty -0.2982256 0.0142843 -20.88 <2e-16 ***
## groups 0.0320033 0.0028152 11.37 <2e-16 ***
## age -0.0039251 0.0002448 -16.04 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3932 on 9906 degrees of freedom
## Multiple R-squared: 0.07695, Adjusted R-squared: 0.07667
## F-statistic: 275.3 on 3 and 9906 DF, p-value: < 2.2e-16Ex3_mo<- mediate(Ex3_m,Ex3_o,treat="age",mediator="poverty", boot=F, sims=1000)
summary(Ex3_mo)##
## Causal Mediation Analysis
##
## Quasi-Bayesian Confidence Intervals
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME 0.000395 -0.002413 0.00 0.83
## ADE -0.003928 -0.004419 0.00 <2e-16 ***
## Total Effect -0.003533 -0.006448 0.00 0.02 *
## Prop. Mediated -0.092565 -3.618450 0.41 0.85
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 9910
##
##
## Simulations: 1000#Empirical Example #4 (Poisson M, Binary Y, and Binary X):
Local_groups$pov <- as.factor(Local_groups$poverty)
Ex4_m <- glm(groups ~ poverty, data = Local_groups, family = "poisson");summary(Ex4_m) ##
## Call:
## glm(formula = groups ~ poverty, family = "poisson", data = Local_groups)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4932 -1.4932 -0.1107 0.7533 6.9838
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.108728 0.009943 10.94 <2e-16 ***
## poverty -0.605335 0.045383 -13.34 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 17226 on 9909 degrees of freedom
## Residual deviance: 17012 on 9908 degrees of freedom
## AIC: 30079
##
## Number of Fisher Scoring iterations: 6Ex4_o <- glm(SRH~poverty+groups, data = Local_groups, family = "binomial");summary(Ex4_o) ##
## Call:
## glm(formula = SRH ~ poverty + groups, family = "binomial", data = Local_groups)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8658 0.4610 0.6419 0.7138 1.2188
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.23741 0.03312 37.36 <2e-16 ***
## poverty -1.33432 0.07498 -17.80 <2e-16 ***
## groups 0.23770 0.02242 10.60 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 10254.8 on 9909 degrees of freedom
## Residual deviance: 9775.4 on 9907 degrees of freedom
## AIC: 9781.4
##
## Number of Fisher Scoring iterations: 4Ex4_mo <- mediate(Ex4_m,Ex4_o,treat="poverty",mediator="groups", boot=F, sims=1000)
summary(Ex4_mo)##
## Causal Mediation Analysis
##
## Quasi-Bayesian Confidence Intervals
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME (control) -0.0173 -0.0215 -0.01 <2e-16 ***
## ACME (treated) -0.0293 -0.0363 -0.02 <2e-16 ***
## ADE (control) -0.2750 -0.3089 -0.24 <2e-16 ***
## ADE (treated) -0.2870 -0.3205 -0.25 <2e-16 ***
## Total Effect -0.3044 -0.3381 -0.27 <2e-16 ***
## Prop. Mediated (control) 0.0568 0.0439 0.07 <2e-16 ***
## Prop. Mediated (treated) 0.0965 0.0750 0.12 <2e-16 ***
## ACME (average) -0.0233 -0.0289 -0.02 <2e-16 ***
## ADE (average) -0.2810 -0.3145 -0.25 <2e-16 ***
## Prop. Mediated (average) 0.0766 0.0596 0.10 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 9910
##
##
## Simulations: 1000Ex2_m <- lm(depress~stress, data = stress_withdraw);summary(Ex2_m) ##
## Call:
## lm(formula = depress ~ stress, data = stress_withdraw)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.0095 -0.4195 -0.1609 0.2498 4.0278
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.79936 0.14331 5.578 6.11e-08 ***
## stress 0.17288 0.02965 5.831 1.63e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6819 on 260 degrees of freedom
## Multiple R-squared: 0.1156, Adjusted R-squared: 0.1122
## F-statistic: 34 on 1 and 260 DF, p-value: 1.63e-08Ex2_o <- lm(withdraw~stress+depress, data = stress_withdraw);summary(Ex2_o) ##
## Call:
## lm(formula = withdraw ~ stress + depress, data = stress_withdraw)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.1716 -0.9472 -0.2249 0.8490 2.9049
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.44706 0.25201 5.742 2.61e-08 ***
## stress -0.07685 0.05239 -1.467 0.144
## depress 0.76913 0.10306 7.463 1.29e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.133 on 259 degrees of freedom
## Multiple R-squared: 0.1804, Adjusted R-squared: 0.174
## F-statistic: 28.49 on 2 and 259 DF, p-value: 6.528e-12Ex2_mo <- mediate(Ex2_m,Ex2_o,treat="stress",mediator="depress", boot=T, sims=10000,boot.ci.type="bca")## Running nonparametric bootstrapsummary(Ex2_mo)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the BCa Method
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME 0.1330 0.0740 0.21 0.0008 ***
## ADE -0.0768 -0.1780 0.04 0.1848
## Total Effect 0.0561 -0.0739 0.17 0.3682
## Prop. Mediated 2.3694 -248.1722 2.43 0.3674
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 262
##
##
## Simulations: 10000The effect of stress on our mediator (depression) is strongly significant. This indicates that stress is positively associated with depression and that a mediation can plausibly be suspected. We later see in our second model that our independent variable loses significance and that the only significant effect on the dependent variable is attributed to our mediator variable. This implies complete mediation. Our causal mediation analysis indicates that the effect on the dv is indirect and that our independent variable must work through depression in order to examine this effect. This is evident with the ACME significance.
Ex3_m2 <- glm(poverty ~age+groups, data = Local_groups, family = binomial(link = "probit"));summary(Ex3_m2) ##
## Call:
## glm(formula = poverty ~ age + groups, family = binomial(link = "probit"),
## data = Local_groups)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.5623 -0.4722 -0.4122 -0.3292 3.8088
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.956301 0.061319 -15.595 < 2e-16 ***
## age -0.005354 0.001134 -4.722 2.33e-06 ***
## groups -0.165291 0.017180 -9.621 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 5743.2 on 9909 degrees of freedom
## Residual deviance: 5599.8 on 9907 degrees of freedom
## AIC: 5605.8
##
## Number of Fisher Scoring iterations: 5Ex3_o2 <- lm(SRH~poverty+groups+age, data = Local_groups);summary(Ex3_o2) ##
## Call:
## lm(formula = SRH ~ poverty + groups + age, data = Local_groups)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.10875 0.06167 0.15693 0.22426 0.67949
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.9837700 0.0137680 71.45 <2e-16 ***
## poverty -0.2982256 0.0142843 -20.88 <2e-16 ***
## groups 0.0320033 0.0028152 11.37 <2e-16 ***
## age -0.0039251 0.0002448 -16.04 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3932 on 9906 degrees of freedom
## Multiple R-squared: 0.07695, Adjusted R-squared: 0.07667
## F-statistic: 275.3 on 3 and 9906 DF, p-value: < 2.2e-16Ex3_mo2 <- mediate(Ex3_m2,Ex3_o2,treat="age",mediator="poverty", boot=F, sims=100)
summary(Ex3_mo2)##
## Causal Mediation Analysis
##
## Quasi-Bayesian Confidence Intervals
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME 0.000353 -0.002709 0.00 0.90
## ADE -0.003884 -0.004335 0.00 <2e-16 ***
## Total Effect -0.003531 -0.006852 0.00 0.02 *
## Prop. Mediated -0.070498 -3.604698 0.43 0.92
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 9910
##
##
## Simulations: 100sens.cont <- medsens(Ex3_mo2, rho.by=.1, eps=.01, effect.type="both")## Warning in rho^2 * (1 - r.sq.m): Recycling array of length 1 in vector-array arithmetic is deprecated.
## Use c() or as.vector() instead.summary(sens.cont)##
## Mediation Sensitivity Analysis for Average Causal Mediation Effect
##
## Sensitivity Region
##
## Rho ACME 95% CI Lower 95% CI Upper R^2_M*R^2_Y* R^2_M~R^2_Y~
## [1,] -0.4 0 0 0 0.16 0.1388
##
## Rho at which ACME = 0: -0.4
## R^2_M*R^2_Y* at which ACME = 0: 0.16
## R^2_M~R^2_Y~ at which ACME = 0: 0.1388
##
##
## Mediation Sensitivity Analysis for Average Direct Effect
##
## Rho at which ADE = 0: -0.9
## R^2_M*R^2_Y* at which ADE = 0: 0.81
## R^2_M~R^2_Y~ at which ADE = 0: 0.7029