The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

Getting Started with RobustMediate

Overview

RobustMediate estimates natural direct and indirect effects (NDE/NIE) for continuous treatments by g-computation from the fitted mediator and outcome models. It adds balance diagnostics, two sensitivity analyses, and a reporting helper (diagnose()).

Simulated data

set.seed(42)
n  <- 600
Z1 <- rnorm(n)
Z2 <- rbinom(n, 1, 0.5)
X  <- 0.5 * Z1 + 0.3 * Z2 + rnorm(n)  # continuous treatment
M  <- 0.4 * X  + 0.2 * Z1 + rnorm(n)  # mediator
Y  <- 0.3 * X  + 0.5 * M  + 0.1 * Z1 + rnorm(n)  # outcome
dat <- data.frame(Y, X, M, Z1, Z2)

Fitting the model

library(RobustMediate)

fit <- robustmediate(
  treatment_formula = X ~ Z1 + Z2,
  mediator_formula  = M ~ X + Z1 + Z2,
  outcome_formula   = Y ~ X + M + Z1 + Z2,
  data    = dat,
  R       = 50,    # use >= 500 in practice
  verbose = FALSE
)
print(fit)
#> -- RobustMediate fit ------------------------------------------
#>   Treatment: X  |  Mediator: M  |  Outcome: Y
#>   N = 600  |  Ref dose = 0.113  |  Bootstrap reps = 50
#> 
#>   Effects at focal dose (1.778 vs ref 0.113):
#>     NDE   0.5469  [0.3761, 0.6731]
#>     NIE   0.3070  [0.2259, 0.3816]
#>     TE    0.8539  [0.6898, 0.9728]
#>     % mediated: 36.0%
#> 
#>   Balance (max |SMD| after weighting):
#>     Treatment pathway: 0.122  (1 covariate(s) above 0.10)
#>     Mediator pathway:  0.058  (0 covariate(s) above 0.10)
#> 
#>   Sensitivity: indirect effect is zero at rho = 0.421 (see sensitivity_rho())
#> ---------------------------------------------------------------

1. Love Plot — plot_balance()

The love plot shows standardised mean differences (SMDs) before and after weighting by stabilised treatment-model weights. The mediator panel is a rough screening check. The dashed lines mark the common |SMD| < 0.10 rule of thumb.

plot_balance(fit)

2. Dose-Response Curve — plot_mediation()

NDE and NIE as smooth curves over the full treatment range, with pointwise bootstrap confidence bands. The dotted vertical line marks the reference dose.

plot_mediation(fit, estimands = c("NDE", "NIE", "TE"))

3. Sensitivity — sensitivity_rho() and plot_sensitivity()

The indirect effect is identified only if no unmeasured variable confounds the mediator-outcome relationship. Following Imai, Keele, and Yamamoto (2010), rho is the correlation between the errors of the mediator and outcome models; rho = 0 under that assumption. The plot shows the indirect effect as rho moves away from zero, and the dashed line marks where it reaches zero.

sensitivity_rho(fit)
#> -- Sensitivity to mediator-outcome confounding (Imai et al., 2010) --
#>   Indirect effect per unit of treatment (rho = 0): 0.1844
#>   Indirect effect for contrast of 1.665 units:       0.3070
#>   Indirect effect is zero at rho = 0.421
#>   Equivalent R2*_M x R2*_Y = 0.1771
#>   Point estimates only; assumes linear models and ignorable treatment.
plot_sensitivity(fit)

The analysis uses a closed-form result for linear models. It is skipped when the mediator or outcome model is not a Gaussian identity-link model, or when the treatment or mediator enters through an interaction or transformation.

4. Diagnose — diagnose()

Prints a formatted report with a ready-to-paste Results paragraph.

diagnose(fit)
#> ==============================================================
#>            RobustMediate: Diagnostics Report
#> ==============================================================
#> 
#> -- Balance diagnostics ---------------------------------------
#>   Treatment pathway: max |SMD| = 0.122  (1 covariate(s) above 0.10)
#>   Mediator pathway:  max |SMD| = 0.058  (0 covariate(s) above 0.10)
#>   WARNING: Balance exceeds 0.10 on at least one covariate.
#>   Inspect plot_balance() carefully.
#> 
#> -- Mediation effects (focal dose = 1.778 vs ref = 0.113 ) ----------------
#>   NDE   +0.5469  [+0.3761, +0.6731]  (95% CI)
#>   NIE   +0.3070  [+0.2259, +0.3816]  (95% CI)
#>   TE    +0.8539  [+0.6898, +0.9728]  (95% CI)
#>   % mediated: 36.0%
#> 
#> -- Sensitivity to mediator-outcome confounding --------------
#>   Indirect effect is zero at rho = 0.421
#>   Equivalent R2*_M x R2*_Y = 0.1771
#> 
#>   Suggested text for Results section:
#>   "Following Imai, Keele, and Yamamoto (2010), the estimated indirect
#>    effect would be zero if the correlation between the errors of the
#>    mediator and outcome models were 0.42. Equivalently, an unmeasured
#>    confounder would need to explain proportions of the residual variance
#>    of the mediator and outcome whose product is 0.177."
#>   Use plot(fit, type = "sensitivity") for visualization.
#> 
#> --------------------------------------------------------------

Path-wise robustness — sensitivity_paths()

sensitivity_paths() gives robustness values in the partial R-squared metric of Cinelli and Hazlett (2020) for the a-path, the b-path, and the indirect effect. For the indirect effect it reports two values: one for the joint significance test (equal to the smaller path value) and one for the Sobel test, which is always smaller because one confounder weakens both paths.

sp <- sensitivity_paths(fit)
sp
#> -- Path-wise and indirect-effect robustness values --------------
#>               quantity estimate    se   stat  df rv_q1 rv_alpha
#>                 a-path    0.396 0.041  9.652 596 0.325    0.269
#>                 b-path    0.466 0.041 11.316 595 0.369    0.317
#>  indirect (joint test)    0.184    NA     NA  NA 0.325    0.269
#>  indirect (Sobel test)    0.184 0.025  7.344  NA 0.325    0.260
#> 
#>   RV(joint) / RV(Sobel) = 1.03
#>   rv_q1: confounder strength that reduces the estimate to zero.
#>   rv_alpha: strength that makes the test non-significant (alpha = 0.05).
#>   If the indirect effect is tested as a product (Sobel, Monte Carlo,
#>   bootstrap), report the Sobel value, not the smaller path value.
plot_paths(sp)

Clustered data

Supply cluster_var and the bootstrap resamples whole clusters:

data(sim_mediation)
fit_clus <- robustmediate(
  treatment_formula = X ~ Z1 + Z2 + Z3,
  mediator_formula  = M ~ X + Z1 + Z2 + Z3,
  outcome_formula   = Y ~ X + M + Z1 + Z2 + Z3,
  data        = sim_mediation,
  cluster_var = "school",
  R           = 500
)

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.