kosukeimai / mediation

R package mediation
61 stars 29 forks source link

Precision of upper confidence bound #36

Closed mspeekenbrink closed 3 years ago

mspeekenbrink commented 3 years ago

When calling summary on a mediate object, the underlying code uses the printCoefmat method. By default, this assumes the results of the second-last column (i.e. the upper-confidence interval here) are a test statistic, which is reported with less precision than other values such the lower confidence interval. That is inconsistent between the bounds, but can also lead to rounding to 0, which can lead to inconsistencies between the confidence intervals and the p-values. For instance

> dat <- read.csv("https://mspeekenbrink.github.io/sdam-jasp-companion/data/redist2015.csv")
> mod2 <- lm(SC_mean ~ income, data=dat)
> mod3 <- lm(redist ~ SC_mean + income, data=dat)
> set.seed(20201027)
> med <- mediate(model.m = mod2, model.y = mod3, sims=1000, boot = TRUE, boot.ci.type = "bca", treat = "income", mediator = "SC_mean")
Running nonparametric bootstrap

> summary(med)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the BCa Method

               Estimate 95% CI Lower 95% CI Upper p-value    
ACME           -0.00229     -0.00463         0.00   0.014 *  
ADE            -0.00292     -0.00650         0.00   0.078 .  
Total Effect   -0.00521     -0.00801         0.00  <2e-16 ***
Prop. Mediated  0.43953      0.12895         1.29   0.014 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Sample Size Used: 301 

Simulations: 1000 

It would be better if the upper confidence interval had the same precision as the lower CI. An easy fix (I believe) is to indicate that there is no test statistic, replacing all calls to printCoefmat to

printCoefmat(smat, tst.ind=NULL)

Removing the digits=3 argument that is used at the moment also allows the digits to be controlled by the general options, rather than hardcoded. E.g.

> x <- med
> smat <- c(x$d1, x$d1.ci, x$d1.p)
> smat <- rbind(smat, c(x$z0, x$z0.ci, x$z0.p))
> smat <- rbind(smat, c(x$tau.coef, x$tau.ci, x$tau.p))
> smat <- rbind(smat, c(x$n0, x$n0.ci, x$n0.p))
> rownames(smat) <- c("ACME", "ADE",
+                     "Total Effect", "Prop. Mediated")
> clp <- "95"
> colnames(smat) <- c("Estimate", paste(clp, "% CI Lower", sep=""),
+                     paste(clp, "% CI Upper", sep=""), "p-value")
> printCoefmat(smat, tst.ind=NULL)
                  Estimate 95% CI Lower 95% CI Upper p-value    
ACME           -0.00228812  -0.00462947  -0.00046455   0.014 *  
ADE            -0.00291772  -0.00649948   0.00024205   0.078 .  
Total Effect   -0.00520584  -0.00801499  -0.00287072  <2e-16 ***
Prop. Mediated  0.43953008   0.12894917   1.29021831   0.014 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
kosukeimai commented 3 years ago

@mspeekenbrink This is a reasonable suggestion. If you like to contribute, a PR is welcome.

mspeekenbrink commented 3 years ago

@kosukeimai Done

kosukeimai commented 3 years ago

Thanks!