[Mediation-information] possible bug in 'mediation-0.4'

[Teppei Yamamoto]

Hi Liviu,

The current version of the mediate function uses the bootstrap for the 
construction of percentile confidence intervals, not for bias 
correction. So for the linear models, the point estimate of the 
mediation effect will be exactly the same as the product of 
coefficients, but the confidence interval will slightly differ from the 
standard one based on the t table.

If you do worry about the bias, you can estimate it using the bootstrap 
draws contained in the output (d0.sims, etc) as

d0.bias <- mean(med_a$d0.sim) - med_a$d0

Then the simplest bias-corrected estimate will be

med_a$d0 - d0.bias

The previous version of the function computed the estimate differently, 
but it was NOT a bias-corrected estimate.

Teppei


Liviu Andronic wrote:
> Hello
> I am wondering whether I spotted a bug in latest 'mediation'. Consider
> the following:
>> require(mediation)
>> summary(mod.m<- lm(EFFICACY ~ LEADER + AGE + SEX + DIPLOMA + DEPART + SUPTENURE, indir))  ##estimate m model
>
> Call:
> lm(formula = EFFICACY ~ LEADER + AGE + SEX + DIPLOMA + DEPART +
>      SUPTENURE, data = indir)
>
> Residuals:
>       Min       1Q   Median       3Q      Max
> -2.92860 -0.43824 -0.06034  0.49659  1.93870
>
> Coefficients:
>                Estimate Std. Error t value Pr(>|t|)
> (Intercept)  4.6369421  0.8721587   5.317 7.92e-07 ***
> LEADER       0.2229302  0.0944945   2.359   0.0205 *
> AGE         -0.0224524  0.0107333  -2.092   0.0393 *
> SEX          0.2862008  0.1745420   1.640   0.1046
> DIPLOMA      0.1807359  0.0977275   1.849   0.0678 .
> DEPART      -0.0154174  0.0102971  -1.497   0.1379
> SUPTENURE   -0.0007209  0.0020428  -0.353   0.7250
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> Residual standard error: 0.7816 on 88 degrees of freedom
>    (17 observations deleted due to missingness)
> Multiple R-squared: 0.1959,	Adjusted R-squared: 0.1411
> F-statistic: 3.574 on 6 and 88 DF,  p-value: 0.003242
>
>> summary(mod.y<- lm(CREATIV ~ EFFICACY + LEADER + AGE + SEX + DIPLOMA + DEPART + SUPTENURE, indir))  ##estimate y model
>
> Call:
> lm(formula = CREATIV ~ EFFICACY + LEADER + AGE + SEX + DIPLOMA +
>      DEPART + SUPTENURE, data = indir)
>
> Residuals:
>       Min       1Q   Median       3Q      Max
> -1.20977 -0.44844 -0.09279  0.37973  1.29044
>
> Coefficients:
>               Estimate Std. Error t value Pr(>|t|)
> (Intercept)  2.184317   0.784891   2.783 0.006606 **
> EFFICACY     0.141007   0.083462   1.689 0.094707 .
> LEADER       0.203025   0.076287   2.661 0.009268 **
> AGE         -0.008247   0.008610  -0.958 0.340814
> SEX          0.246680   0.138728   1.778 0.078871 .
> DIPLOMA      0.037031   0.077987   0.475 0.636094
> DEPART      -0.030152   0.008164  -3.693 0.000386 ***
> SUPTENURE   -0.002973   0.001601  -1.857 0.066644 .
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> Residual standard error: 0.612 on 87 degrees of freedom
>    (17 observations deleted due to missingness)
> Multiple R-squared: 0.3114,	Adjusted R-squared: 0.256
> F-statistic:  5.62 on 7 and 87 DF,  p-value: 2.275e-05
>
>> med_a<- mediate(mod.m, mod.y, dropobs=T, treat='LEADER', mediator='EFFICACY', boot=T, sims=5000)  ##supply m and y models for mediation analysis
>
>> summary(med_a)
>
> Causal Mediation Analysis
>
> Confidence Intervals Based on Nonparametric Bootstrap
>
>                           Estimate 95% CI Lower 95% CI Upper p-value
> Mediation Effect           0.0314      -0.0228       0.0842    0.23
> Direct Effect              0.2030       0.0424       0.3672    0.01
> Total Effect               0.2345       0.0765       0.3834    0.00
> Proportion via Mediation   0.1341      -0.1035       0.4998    0.23
>
> Sample Size Used: 95
>
>
> Now the trouble is, I think, with the 'Mediation Effect' estimate: it
> doesn't seem to be the bootstrap estimate, but the estimate obtained
> directly from hte regressions. If you use the regression coefficients
> you get this:
> #a*b=0.0314
>> (med_axb<- as.numeric(coef(mod.m)['LEADER'] * coef(mod.y)['EFFICACY'])) #data
> [1] 0.03143479
>
> And the 'medaite' object contains this:
>> med_a$d0 #??bug
> [1] 0.03143479
>
> If, however, you compute the average over the bootstrap estimates as
> stored in the object you get a different estimate, and this is what I
> would expect printed in the summary() output above:
>> mean(med_a$d0.sims) #boot/acme
> [1] 0.02689732
>
> Computing the bias worked fine with the previosu version of
> 'mediation', but now it yields zero:
>> med_a$d0 - med_axb #bias
> [1] 4.857226e-17
>
>
> I may be wrong but I think there is a bug, and I am not sure if it
> affects only the 'Mediation Effect' or some other coefficients (or the
> CI).
>
> Regards
> Liviu
>
>

-- 
====================================
Teppei Yamamoto
Assistant Professor
Department of Political Science
Massachusetts Institute of Technology
http://web.mit.edu/teppei/www/
====================================