R implementation for causal moderated mediation analysis
This document illustrates the implementation of the R package
moderate.mediation. The package is for researchers to
estimate and test the conditional and moderated mediation effects,
assess their sensitivity to unmeasured pretreatment confounding, and
visualize the results. We built the package based on the Monte Carlo
method (default) and the bootstrapping method. When the sample size is
small (e.g., n = 50 or 100 under the studied conditions in the
simulations), and the mediator or the outcome is binary, we recommend
the bootstrapping method, because it generates point estimates with
smaller bias than the Monte Carlo method. Otherwise, we recommend the
Monte Carlo method, because its running speed is much faster and its
statistical power is slightly higher at a small sample size when both
the mediator and outcome are continuous. Because a relatively large
sample size is usually required to detect a significant moderated
mediation effect, we set the Monte Carlo method to be the default
method. The package is applicable to a binary or continuous treatment, a
binary or continuous mediator, a binary or continuous outcome, and one
or more moderators of any scale.
Loading the package:
library(moderate.mediation)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:
data(newws)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. Estimation and inference
The estimation and inference are conducted via the
modmed function. A challenge of moderated mediation
analysis lies in the specification of the mediator and outcome models,
which can be very complex as more interactions are involved. To ease
model building, the function allows users to specify the main model
predictors and the moderators of each main model coefficient following
the hierarchical form. In the NEWWS example, with a binary treatment, a
binary mediator, and a continuous outcome, we assume that paths T➔M,
T➔Y, M➔Y, and TM➔Y are moderated by the number of children and welfare
amount at baseline, while X➔M and X➔Y are not moderated. We are
interested in the difference in the mediation mechanism between those
who had three or more children and those who had two children at
baseline, while baseline welfare amount is set at its median, 5050.
Hence, we specified the function as follows:
data(newws)
results = modmed(data = newws,
treatment = "treat",
mediator = "emp",
outcome = "depression",
covariates.disc = c("emp_prior", "nevmar", "hispanic", "nohsdip"),
covariates.cont = c("workpref", "attitude", "depress_prior"),
moderators.disc = "CHCNT",
moderators.cont = "ADCPC",
m.model = list(intercept = c("ADCPC", "CHCNT"), treatment =
c("ADCPC", "CHCNT"), emp_prior = NULL,
nevmar = NULL, hispanic = NULL, nohsdip = NULL, workpref = NULL,
attitude = NULL, depress_prior = NULL),
y.model = list(intercept = c("ADCPC", "CHCNT"), treatment =
c("ADCPC", "CHCNT"), mediator = c("ADCPC", "CHCNT"), tm =
c("ADCPC", "CHCNT"), emp_prior = NULL, nevmar = NULL, hispanic =
NULL, nohsdip = NULL, workpref = NULL, attitude = NULL,
depress_prior = NULL),
comp.treatment.value = 1,
ref.treatment.value = 0,
comp.mod.disc.values = 3,
ref.mod.disc.values = 2,
comp.mod.cont.values = 5050,
ref.mod.cont.values = 5050,
m.scale = "binary",
y.scale = "continuous",
method = "mc",
nmc = 1000,
conf.level = 0.95,
seed = 1) Variable names Users are required to specify the
names of the data, treatment, mediator, and outcome via the arguments,
data, treatment, mediator, and
outcome. If there are any missing values in the data, users
need to impute them before running the function. Whenever applicable,
names of the discrete and continuous moderators \(\bf{W}\) and other pretreatment covariates
\(\bf{X}\) are specified separately via
moderators.disc, moderators.cont,
covariates.disc, and covariates.cont, whose
default values are NULL. For example, if there are no
discrete moderators, users do not need to specify
moderators.disc. Users do not need to reformat discrete
variables. No matter if the discrete variables are coded as string or
numerical values, the program can automatically factorize them. This
eases data cleaning and avoids mistakes.
\(\bf{W}\) contains pretreatment covariates that moderate the mediation mechanism. As clarified in the definition and identification sections, \(\bf{W}\) is chosen based on theoretical interest. \(\bf{X}\) contains other pretreatment covariates that are confounders within levels of \(\bf{W}\). If treatment is randomized, \(\bf{X}\) should contain confounders of the mediator–outcome relationship within levels of \(\bf{W}\). If treatment is not randomized, \(\bf{X}\) should also contain confounders of the treatment–mediator and treatment– outcome relationships within levels of \(\bf{W}\). If \(\bf{W}\) is not specified, only the population average effects are estimated and tested.
Model specification Users can specify the mediator
and outcome models via m.model and y.model,
each as a list. The names of the elements in each list include the
intercept and the predictors in the corresponding main model.
Specifically, in m.model, the names must include
intercept, treatment, and the covariates in
covariates.disc and covariates.cont that
predict the mediator (i.e., a covariate that only confounds the
treatmentoutcome relationship does not need to be included). In
y.model, the names must include intercept,
treatment, mediator, and the covariates in
covariates.disc and covariates.cont that
predict the outcome (i.e., a covariate that only confounds the
treatment-mediator relationship does not need to be included). If the
treatment is assumed to interact with the mediator when affecting the
outcome, an additional element should be added to y.model,
named as tm. We suggested users always include
tm unless it barely changes the final estimates. As
VanderWeele (2015) wrote on page 46, “An investigator might be tempted
to only include such exposure-mediator interactions in the model if the
interaction is statistically significant. This approach is problematic.
It is problematic because power to detect interaction tends to be very
low unless the sample size is very large. Such exposure-mediator
interaction may be important in capturing the dynamics of mediation. A
better approach is perhaps to include them by default and only exclude
them if they do not seem to change the estimates of the direct and
indirect effects very much”.
Each element of each list is equal to a vector of moderators in the
function of the corresponding main model coefficient. Each moderator
specified in moderators.disc and
moderators.cont must moderate at least one slope in either
the main mediator model or the main outcome model. The vector
corresponding to the intercept must contain all the moderators in the
corresponding model because their coefficients represent the main
effects of the moderators. If a main model coefficient is not moderated,
the corresponding vector of moderators needs to be specified as
NULL. The set of moderators in m.model and
that in y.model are not necessarily the same. Specifically,
for a moderator that moderates coefficient(s) in one model, it is
possible that it does not occur in the other model; it is also possible
that it is included in the other model as a predictor rather than a
moderator, in which case users will only include it in the vector
corresponding to the intercept when specifying the list for the other
model. The moderators that are included in both the mediator model and
the outcome model essentially also confound the mediator–outcome
relationship. The union of the two sets of moderators in
m.model and y.model should be the same as the
union of moderators.disc and
moderators.cont.
In the NEWWS example, it is assumed that the paths T➔M, T➔Y, M➔Y, and
TM➔Y are moderated by both "CHCNT" and
"ADCPC".
Variable values and scales Users can specify t and
t’, respectively, via comp.treatment.value and
ref.treatment.value, where "comp" stands for
comparison, while "ref" stands for reference. t = 1 and t’
= 0 by default. In other words, if treatment is binary, one does not
need to specify t. Similarly, users can specify vectors \(\bf{w}_1\) and \(\bf{w}_2\), respectively, through
comp.mod.disc.values (and/or
comp.mod.cont.values) and ref.mod.disc.values
(and/or ref.mod.cont.values), whose default values are
NULL. To be specific, if the focal interest is in the
conditional effects rather than the moderated effects, users do not need
to specify comp.mod.disc.values or
comp.mod.cont.values; if there are no discrete moderators,
users do not need to specify comp.mod.disc.values or
ref.mod.disc.values; if there are no continuous moderators,
users do not need to specify comp.mod.cont.values or
ref.mod.cont.values. If not NULL, the length
and order of each value vector should be consistent with those of the
corresponding name vector.
If one does not want to condition some moderators on specific values,
one may specify their values to be NA. In the NEWWS
example, if we want to check how each effect varies by CHCNT over all
the individuals rather than those whose ADCPC is at 5050, we could
replace 5050 with NA in the above syntax.
The mediator and the outcome are continuous by default. If the
mediator (outcome) is binary, m.scale (y.scale) needs to be
specified as "binary", and a probit regression is
fitted.
In the NEWWS example, we evaluated the difference in the mediation
mechanism between CHCNT = 3 vs. CHCNT= 2 while
ADCPC is set at its median (5050).
Other arguments There are additional arguments that
users do not need to specify unless they want to change the default.
Users can specify the estimation method via method, which
is set to "mc" if the Monte Carlo method is chosen and to
"boot" if the bootstrapping method is chosen (default =
"mc"). If method = "mc", one can specify the
number of simulations that the Monte Carlo algorithm takes via
nmc (default = 1000); If method = "boot", one
can specify the number of bootstrapped samples via nboot
(default = 1000). conf.level, which indicates the
confidence level, is set to 0.95 by default. Users can also specify
seed (default = NULL).
2. Numerical summary of analysis results
The output of the modmed function can be passed via
object to the summary_modmed function for
numerical summary of the analysis results.
summary_modmed(object = results)## Treatment: treat
##
## Mediator: emp
##
## Outcome: depression
##
## Pre-treatment confounders: emp_prior, nevmar, hispanic, nohsdip, workpref, attitude, depress_prior
##
## Moderators: CHCNT, ADCPC
##
## Compare values of the treatment: 1
##
## Reference values of the treatment: 0
##
## Compare values of the moderators: 3, 5050
##
## Reference values of the moderators: 2, 5050
##
## Estimation method: Monte Carlo Method
##
## Causal Effects:
## Estimate Std. Error 95% CI Lower 2.5% 95% CI Upper 2.5%
## TE 0.47925508 0.6716924 -0.85981724 1.77596645
## TIE -0.76498887 0.3527762 -1.47218136 -0.09806794
## PIE -0.01629142 0.2270166 -0.44934204 0.45481519
## PDE 1.24424395 0.7224715 -0.19591531 2.61923999
## TDE 0.49554650 0.7115971 -0.91393794 1.90552911
## INT -0.74869745 0.4205389 -1.61948024 0.02738178
## TE.ref 2.19848837 1.1361880 0.05964914 4.44545134
## TIE.ref -0.95614534 0.6924313 -2.49873699 0.30463453
## PIE.ref 0.16624820 0.4421947 -0.72998825 1.05244606
## PDE.ref 3.15463371 1.2675048 0.70577443 5.58663299
## TDE.ref 2.03224017 1.2354353 -0.40493093 4.44795983
## INT.ref -1.12239354 0.8372979 -2.89434696 0.50936546
## TE.dif -3.26938152 1.5823285 -6.37740944 -0.06492951
## TIE.dif 0.43697787 0.7531679 -0.95055948 2.10809748
## PIE.dif -0.16362524 0.4688132 -1.05735349 0.74375372
## PDE.dif -3.70635939 1.6847718 -6.92068720 -0.33509981
## TDE.dif -3.10575628 1.6622477 -6.44243372 0.12034522
## INT.dif 0.60060311 0.9031720 -1.09897932 2.53426161
## ---
## CI is confidence interval constructed based on simulation of mediator and outcome model parameters (number of simulations is 1000)
##
## Mediator Model:
## glm(formula = as.formula(paste(mediator, "~", formula.m)), family = binomial(link = "probit"),
## data = data)
## main moderation Estimate Std. Error z value
## (Intercept) Intercept Intercept -8.994616e-01 2.997011e-01 -3.0011954
## ADCPC ADCPC -5.120716e-05 1.770930e-05 -2.8915407
## CHCNT2 CHCNT2 -1.034358e-01 1.519718e-01 -0.6806254
## CHCNT3 CHCNT3 1.115883e-01 1.635839e-01 0.6821473
## treat treat Intercept 5.695776e-01 2.275527e-01 2.5030574
## treat:ADCPC ADCPC 6.269787e-05 2.981701e-05 2.1027554
## treat:CHCNT2 CHCNT2 7.564875e-02 2.861609e-01 0.2643574
## treat:CHCNT3 CHCNT3 -3.904822e-01 2.890170e-01 -1.3510698
## emp_prior1 emp_prior1 Intercept 7.927838e-01 1.171389e-01 6.7678976
## nevmar1 nevmar1 Intercept -4.622593e-02 1.116981e-01 -0.4138471
## hispanic1 hispanic1 Intercept 7.791587e-02 1.156756e-01 0.6735721
## nohsdip1 nohsdip1 Intercept -2.443537e-01 1.084837e-01 -2.2524461
## workpref workpref Intercept 1.583547e-01 7.509759e-02 2.1086518
## attitude attitude Intercept 7.699394e-02 8.905917e-02 0.8645257
## depress_prior depress_prior Intercept 3.298730e-02 6.207356e-02 0.5314228
## Pr(>|z|)
## (Intercept) 2.689220e-03 **
## ADCPC 3.833579e-03 **
## CHCNT2 4.961086e-01
## CHCNT3 4.951458e-01
## treat 1.231256e-02 *
## treat:ADCPC 3.548716e-02 *
## treat:CHCNT2 7.915046e-01
## treat:CHCNT3 1.766731e-01
## emp_prior1 1.306672e-11 ***
## nevmar1 6.789861e-01
## hispanic1 5.005834e-01
## nohsdip1 2.429410e-02 *
## workpref 3.497465e-02 *
## attitude 3.872992e-01
## depress_prior 5.951258e-01
## ---
## Signif.codes:
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Null deviance: 960.01 on 693 degrees of freedom
## Residual deviance: 817.6 on 679 degrees of freedom
## AIC: 847.5997
##
## Outcome Model:
## lm(formula = as.formula(paste(outcome, "~", formula.y)), data = data)
## main moderation Estimate Std. Error
## (Intercept) Intercept Intercept 4.058645e+00 1.7907773457
## ADCPC ADCPC -2.571579e-05 0.0001247530
## CHCNT2 CHCNT2 1.246529e-01 1.0925244717
## CHCNT3 CHCNT3 6.166837e-01 1.1734409467
## treat treat Intercept 2.337643e+00 2.1719542091
## treat:ADCPC ADCPC 2.370110e-05 0.0002562197
## treat:CHCNT2 CHCNT2 2.010385e+00 2.6829663548
## treat:CHCNT3 CHCNT3 -1.274618e+00 2.6346287059
## emp emp Intercept 7.131112e-01 1.3615226098
## emp:ADCPC ADCPC -3.045652e-04 0.0002063319
## emp:CHCNT2 CHCNT2 1.472583e+00 1.7034188039
## emp:CHCNT3 CHCNT3 1.440190e+00 1.8078479481
## treat:emp treat:emp Intercept -4.054790e+00 2.6922253282
## treat:emp:ADCPC ADCPC 3.711774e-04 0.0003537832
## treat:emp:CHCNT2 CHCNT2 -1.265219e+00 3.4022130215
## treat:emp:CHCNT3 CHCNT3 -1.654510e+00 3.4129707843
## emp_prior1 emp_prior1 Intercept 6.827357e-01 0.6888712767
## nevmar1 nevmar1 Intercept 1.232385e+00 0.6295330118
## hispanic1 hispanic1 Intercept 1.332704e+00 0.6503718578
## nohsdip1 nohsdip1 Intercept 4.537165e-01 0.6174878435
## workpref workpref Intercept -2.050631e-01 0.4218136727
## attitude attitude Intercept -5.015208e-01 0.5017808235
## depress_prior depress_prior Intercept 1.817711e+00 0.3497288888
## t value Pr(>|t|)
## (Intercept) 2.26641544 2.374355e-02 *
## ADCPC -0.20613356 8.367491e-01
## CHCNT2 0.11409619 9.091957e-01
## CHCNT3 0.52553451 5.993853e-01
## treat 1.07628557 2.821864e-01
## treat:ADCPC 0.09250305 9.263260e-01
## treat:CHCNT2 0.74931432 4.539304e-01
## treat:CHCNT3 -0.48379404 6.286899e-01
## emp 0.52376009 6.006184e-01
## emp:ADCPC -1.47609339 1.403879e-01
## emp:CHCNT2 0.86448699 3.876294e-01
## emp:CHCNT3 0.79663237 4.259463e-01
## treat:emp -1.50611089 1.325093e-01
## treat:emp:ADCPC 1.04916636 2.944793e-01
## treat:emp:CHCNT2 -0.37188122 7.100987e-01
## treat:emp:CHCNT3 -0.48477115 6.279969e-01
## emp_prior1 0.99109323 3.219973e-01
## nevmar1 1.95761844 5.068891e-02 .
## hispanic1 2.04914105 4.083669e-02 *
## nohsdip1 0.73477809 4.627314e-01
## workpref -0.48614618 6.270223e-01
## attitude -0.99948178 3.179217e-01
## depress_prior 5.19748615 2.684134e-07 ***
## ---
## Signif.codes:
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.52 on 671 degrees of freedom
## Multiple R-squared: 0.0871 Adjusted R-squared: 0.05715
## F-statistic: 2.909 on 22 and 671 DF, p-value: 0The output includes the estimation and inference results of the population average, conditional, and/or moderated effects. In addition, it also includes the fitting results of the mediator and outcome models, which allow one to assess the moderation of not only each mediation effect but also each specific path. The model fitting results are shown in the hierarchical form and thus better demonstrate how each specific path is moderated.
The results for the NEWWS data show that the total LFA effect is estimated to be 0.48 (SE = 0.67, 95% CI = [-0.86, 1.78]), which can be decomposed into a pure (natural) direct effect, estimated to be 1.24 (SE = 0.72, 95% CI = [–0.20, 2.62]), about 16.07% of a standard deviation of the outcome in the control group; and a total (natural) indirect effect, estimated to be -0.76 (SE = 0.35, 95% CI = [-1.47, -0.10]), about -9.88% of a standard deviation of the outcome in the control group. The pure 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. In contrast, the total indirect effect reflects that the LFA-induced increase in employment rate significantly reduced one’s maternal depression, when the treatment is held at the LFA condition. The counteracting indirect and direct effects explained the insignificant total effect.
The total indirect effect can be further decomposed into a pure indirect effect, which is estimated to be -0.02 (SE = 0.23, 95% CI = [-0.45, 0.45]); and a natural treatment-bymediator interaction effect, which is estimated to be -0.75 (SE = 0.42, 95% CI = [-1.62, 0.03]). It indicates that the LFA-induced increase in employment rate reduced maternal depression more under the LFA condition (i.e., total 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 pure 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 than if holding that under the control condition.
Through the moderated mediation analysis, we further detected a significantly positive pure direct effect among those who had two children at baseline and received median level ($5050) of welfare in the year prior to randomization, which is estimated to be 3.15 (SE = 1.27, 95% CI = [0.71, 5.59]). This conditional pure direct effect is significantly higher than that of those with three children and $5,050 welfare in the year prior to randomization. The magnitude of the difference is estimated to be 3.71 (SE = 1.68, 95% CI = [0.34, 6.92]).
The above numerical summary provides important information about the conditional effects at a given set of values of the moderators and how they differ between two given sets of values of the moderators. To further assess how the effects vary across the whole value range of a moderator, the R package offers a graphical tool.
3. Graphical summary of analysis results
The output of the modmed function can be passed via
object to the modmed.plot function for
generating a whole picture of how the mediation mechanism differs by a
moderator.
modmed.plot(object = results,
effect = "PDE",
moderator = "CHCNT",
other.mod.cont.values = 5050)
Plot of a conditional causal effect versus a
moderator Each time users choose one focal effect (specified
via effect, which can be "TE",
"TIE", "PDE", "PIE",
"TDE", or "INT") to be plotted against a focal
moderator (specified via moderator), while conditioning on
given values of the other moderators (specified via
other.mod.disc.values for discrete moderators and/or
other.mod.cont.values for continuous moderators). The
length and order of other.mod.disc.values must match those
of moderators.disc (specified in the modmed
function) with the focal moderator removed if it is discrete. If one
does not want to condition some moderators on specific values, one may
specify their values to be NA. The same is true for
other.mod.cont.values.
Applying the above syntax to the NEWWS data, we obtained the above
figure, which represents how the pure direct effect varied with the
number of children that participants had at baseline, among those who
received $5050 welfare in the year before randomization. For each
category of the moderator, a boxplot is generated to represent the Q
estimates of the conditional PDE obtained from the Monte Carlo or
bootstrapping algorithm. Each blue dot in the middle represents the
final effect estimate, while each blue interval denotes the confidence
interval at the confidence level specified in the modmed
function.
Among those who received $5050 welfare in the year before
randomization, the pure direct effect of LFA on maternal depression is
estimated to be positive among those with one or two children at
baseline, while significant only among those with two children. In
contrast, among those with three or more children, LFA would have
reduced maternal depression if not changing one’s employment experience,
but not by a statistically significant amount. This matches and
supplements the numerical summary obtained from the
summary_modmed function.
Conditional distribution of a continuous moderator
If the focal moderator is continuous, one can choose to add its sample
distribution given values of the other moderators via
is.dist.moderator = TRUE. It will not be plotted if the
subsample within the given levels of the other moderators has fewer than
30 observations. Five percentile lines are added to the distribution of
the moderator, which facilitates the evaluation of conditional effects
at given percentiles of each continuous moderator. The percentiles can
be specified via probs. The default is
c(0.1, 0.25, 0.5, 0.75, 0.9).
By running the following syntax, we further assessed how the pure direct effect varied with welfare amount in the year before randomization and how the welfare amount is distributed, among those who had two children at baseline.
modmed.plot(object = results,
effect = "PDE",
moderator = "ADCPC",
other.mod.disc.values = 2,
is.dist.moderator = TRUE,
probs = c(0.1, 0.25, 0.5, 0.75, 0.9),
ncore = 6)
As shown in the above figure, the estimates of the conditional PDE are represented by the blue curve. To obtain such a curve, the function estimates the conditional PDE at each of 15 evenly spaced values of the moderator and fits a loess line to the 15 estimates. Similarly, we obtained the grey band, which represents the confidence band. The dashed vertical lines divide the significant and insignificant regions. The limit of the x-axis is determined by the range of the focal moderator in the whole sample.
From the figure, among those who had two children at baseline, the pure direct effect of LFA on maternal depression increased as one received more welfare amount in the year prior to randomization. The effect is significant when the welfare amount is larger than $2526. The sample distribution of the welfare amount indicates that about 70% of the participants with two children are in the significant region.
By changing “PDE” to the other effects, one at a time, we can obtain
a figure for each effect. These figures provide researchers with a
general picture of how each effect changes with a moderator given values
of the other moderators, as well as the moderator values at which the
sign or significance of the effect is flipped. If particularly
interested in a visualized difference in an effect between two values of
the focal moderator, one may use the modmed function to
estimate the difference and test its significance. If the figure shows
that an effect is linearly associated with a continuous moderator, one
may estimate and test the difference in the effect between two
arbitrarily chosen moderator values that are one unit apart, which
implies the average change in the effect with one unit increase in the
moderator.
4. Sensitivity analysis
In a causal mediation analysis, it is important to theoretically
discuss which unmeasured variables may be confounders, assess how
plausible it is, and evaluate how additional adjustment for these
potential unmeasured confounders would have changed the sign or
significance of the effects. The sensitivity of the analysis results to
unmeasured pretreatment confounding can be assessed by passing the
output of the modmed function via object to
the modmed.sens function.
sens.results = modmed.sens(object = results,
sens.effect = c("TIE", "PIE",
"TDE", "PDE","INT", "TIE.ref",
"PIE.ref", "PDE.ref", "TDE.ref",
"INT.ref", "TIE.dif", "PIE.dif",
"PDE.dif", "TDE.dif","INT.dif"),
range.b.m = NULL,
range.b.y = NULL,
grid.b.m = 10,
grid.b.y = 10,
U.scale = "binary",
p.u = 0.5,
t.rand = TRUE,
t.model = NULL,
t.scale = "binary",
b.t = NULL,
iter = 10,
nsim = 5,
ncore = 6)The effects whose sensitivity will be evaluated can be specified via
sens.effect. It is a vector,
c("TIE", "PIE", "TDE", "PDE", "INT", "TIE.ref", "PIE.ref", "PDE.ref", "TDE.ref", "INT.ref", "TIE.dif", "PIE.dif", "PDE.dif","TDE.dif", "INT.dif"),
by default (i.e., if sens.effect is not specified).
sens.effect can also be specified as a subvector of the
default. It does not matter how the effects are ordered. For example, if
a researcher is mainly interested in the sensitivity of moderated
indirect effects, sens.effect can be specified as
c("PIE.dif", "TIE.dif"). Following the algorithm for
sensitivity analysis, users can specify value ranges for the sensitivity
parameters via range.b.m and range.b.y (e.g.,
c(-2, 2)); the values are automatically generated if NULL
is specified. They are divided into a
grid.b.m-by-grid.b.y grid, which is 10-by-10
by default. The finer the grid, the smoother the sensitivity analysis
plot as shown later. The sensitivity parameters are the slopes of U in
the standardized mediator and outcome models, in which both the
independent and dependent variables are standardized. If the dependent
variable is binary, its latent index is standardized instead (Grace et
al., 2018). The variance of the latent index is the sum of the predictor
variance and the assumed error variance (McKelvey & Zavoina,
1975).
As clarified in the sensitivity analysis section, the conditional
distribution of U also relies on the marginal distribution of U, which
is assumed to be binary (U.scale = "binary") with Pr(U = 1)
= 0.5 (p.u = 0.5) by default. One can change
U.scale to be "continuous" for a continuous U,
which is assumed to follow a standard normal distribution
(sigma.u = 1) by default. p.u and
sigma.u can be adjusted.
By default, the treatment is randomized (t.rand = TRUE).
In other words, there is no confounder of the treatment– mediator or
treatment–outcome relationship. Hence, by default,
t.model = NULL. If treatment is not randomized
(t.rand = FALSE), one needs to specify t.model
in the same way as specifying m.model and
y.model, because a treatment model is required for the
derivation of the conditional distribution of the unmeasured
pretreatment confounder U. The names of the elements in
t.model must include intercept and the
covariates in covariates.disc and
covariates.cont that predict the treatment (i.e., a
covariate that only confounds the mediator-outcome relationship does not
need to be included). Each element of the list is equal to a vector of
moderators in the function of the main treatment model coefficient. The
moderators of the intercept must contain all the moderators in the
treatment model because their slopes represent the main effects of the
moderators. If a main model coefficient is not moderated, the
corresponding vector of moderators needs to be specified as
NULL. A moderator in the mediator and outcome models does
not necessarily moderate any coefficient in the treatment model. If it
is only a predictor rather than a moderator in the treatment model,
users will only include it in the vector corresponding to the intercept
when specifying t.model. For example, suppose the treatment
is not randomized, we may specify that t.model in the
modmed.sens function as
t.model = list(intercept = c("ADCPC", "CHCNT"), emp_prior = "ADCPC", nevmar = "CHCNT", hispanic = NULL, nohsdip = NULL, workpref = NULL, attitude = NULL, depress_prior = NULL).
Besides, a value for sensitivity parameter \(\beta^t_u\) needs to be specified via
b.t. One may use the standardized coefficient estimates of
the observed covariates in the treatment model as referent values for
the specification of b.t. The treatment is binary by
default. If the treatment is continuous, t.scale needs to
be specified as "continuous". The sensitivity analysis is
conducted given one single value of \(\beta^t_u\). If there are multiple values
of \(\beta^t_u\) to be assessed, one
may conduct the sensitivity analysis multiple times to assess the
influence of each value of \(\beta^t_u\), one at a time.
In addition, users can specify the number of iterations that the
stochastic EM algorithm takes for simulating U (i.e., S) and the number
of simulations of U (i.e., K) respectively via iter and
nsim, which are 10 and 5 by default. Each run of the
modmed function only takes seconds, while the
modmed.sens function involves nsim runs of the
modmed function, each of which involves iter
iterations, for each cell of the grid. To speed up calculation, the
package allows parallel computing by setting the number of CPU cores via
ncore. Its default value is 2. A user may use more cores to
speed up the program. detectCores() can be used to detect
the number of all cores available on a computer. It is not recommended
to use up all the cores. At least one core should be saved for users to
run other programs on the computer while running the R program. A
progress bar is displayed as the program runs. An application of the
above syntax to the NEWWS data took about 1 hour. It takes less time if
sens.effect is specified as a subvector of all the effects.
For example, it takes 15 minutes if
sens.effect = "TIE.ref". We reduced the size of the grid to
5-by-5 and obtained similar results for all the effects in 15 minutes,
while the visualized results are less smooth. We would suggest users use
a coarse grid (e.g., 5-by- 5) to get a preliminary sense of sensitivity
and then use a fine grid (e.g., 10-by-10) for the final report of the
sensitivity analysis results. In the cases when Step 2.5 does not
involve a spline regression, it takes less than 10 minutes for a 10-by-
10 grid with only one core.
The output includes a grid.b.m-by-grid.b.y
table separately for the point estimates, CI lower bounds, and CI upper
bounds of each effect after adjusting for U at the strengths reflected
by the column and row names. Each column (row) name corresponds to one
of the grid.b.y (grid.b.m) sensitivity
parameter values that evenly divide range.b.y
(range.b.m). The table can be extracted via
sens.results$results.new. The output also includes slopes
of each observed pretreatment confounder in the standardized mediator
and outcome models, which can be used as referent values to calibrate
the strength of unmeasured pretreatment confounding. The slopes can be
extracted via sens.results$X.coef.plot.
In the NEWWS example, as illustrated in the identification section,
one’s salary expectation at baseline is a potential unmeasured
confounder of the relationship between employment and depression within
levels of baseline welfare amount and number of children. To assess how
an additional adjustment for such an unmeasured confounder would affect
the initial results, one may specify the sensitivity parameter values
based on prior knowledge about its association with employment and
depression or a comparable observed pretreatment covariate. For example,
it may be reasoned that the unique confounding role of baseline salary
expectation is similar to that of baseline employment. Hence, one may
use the standardized coefficient estimates of baseline employment in the
mediator and outcome models as referent values when specifying the
sensitivity parameters, by setting range.b.m and
range.b.m to the values and setting grid.b.m
and grid.b.m to 1:
modmed.sens(object = results,
range.b.m = 0.31,
range.b.y = 0.04,
grid.b.m = 1,
grid.b.y = 1)## $b.t
## NULL
##
## $X.coef.plot
## X.coef.plot.m X.coef.plot.y
## emp_prior1 0.31030821 0.04068706
## nevmar1 -0.01940212 0.07875448
## hispanic1 0.03084250 0.08031976
## nohsdip1 -0.10357677 0.02928150
## workpref 0.11993452 -0.02364645
## attitude 0.04843231 -0.04803222
## depress_prior 0.02339182 0.19624892
##
## $range.b.m
## [1] 0.31
##
## $range.b.y
## [1] 0.04
##
## $b.y.all
## [1] 0.04
##
## $b.m.all
## [1] 0.31
##
## $results.new
## $results.new$TIE
## 0.04
## 0.31 -0.8429615
##
## $results.new$CIL.TIE
## 0.04
## 0.31 -1.5418
##
## $results.new$CIU.TIE
## 0.04
## 0.31 -0.2113415
##
## $results.new$PIE
## 0.04
## 0.31 -0.05096788
##
## $results.new$CIL.PIE
## 0.04
## 0.31 -0.4858237
##
## $results.new$CIU.PIE
## 0.04
## 0.31 0.3811576
##
## $results.new$PDE
## 0.04
## 0.31 1.326268
##
## $results.new$CIL.PDE
## 0.04
## 0.31 -0.08569929
##
## $results.new$CIU.PDE
## 0.04
## 0.31 2.722105
##
## $results.new$TDE
## 0.04
## 0.31 0.5342744
##
## $results.new$CIL.TDE
## 0.04
## 0.31 -0.797396
##
## $results.new$CIU.TDE
## 0.04
## 0.31 1.881271
##
## $results.new$INT
## 0.04
## 0.31 -0.7919937
##
## $results.new$CIL.INT
## 0.04
## 0.31 -1.626356
##
## $results.new$CIU.INT
## 0.04
## 0.31 -0.03102925
##
## $results.new$TIE.ref
## 0.04
## 0.31 -1.012629
##
## $results.new$CIL.TIE.ref
## 0.04
## 0.31 -2.441112
##
## $results.new$CIU.TIE.ref
## 0.04
## 0.31 0.2252989
##
## $results.new$PIE.ref
## 0.04
## 0.31 0.1141786
##
## $results.new$CIL.PIE.ref
## 0.04
## 0.31 -0.7406788
##
## $results.new$CIU.PIE.ref
## 0.04
## 0.31 0.9980871
##
## $results.new$PDE.ref
## 0.04
## 0.31 3.203819
##
## $results.new$CIL.PDE.ref
## 0.04
## 0.31 0.7700738
##
## $results.new$CIU.PDE.ref
## 0.04
## 0.31 5.683864
##
## $results.new$TDE.ref
## 0.04
## 0.31 2.077011
##
## $results.new$CIL.TDE.ref
## 0.04
## 0.31 -0.315912
##
## $results.new$CIU.TDE.ref
## 0.04
## 0.31 4.489824
##
## $results.new$INT.ref
## 0.04
## 0.31 -1.126808
##
## $results.new$CIL.INT.ref
## 0.04
## 0.31 -2.855432
##
## $results.new$CIU.INT.ref
## 0.04
## 0.31 0.3395653
##
## $results.new$TIE.dif
## 0.04
## 0.31 0.3899642
##
## $results.new$CIL.TIE.dif
## 0.04
## 0.31 -1.087823
##
## $results.new$CIU.TIE.dif
## 0.04
## 0.31 1.922661
##
## $results.new$PIE.dif
## 0.04
## 0.31 -0.1220958
##
## $results.new$CIL.PIE.dif
## 0.04
## 0.31 -1.111169
##
## $results.new$CIU.PIE.dif
## 0.04
## 0.31 0.8530015
##
## $results.new$PDE.dif
## 0.04
## 0.31 -3.613411
##
## $results.new$CIL.PDE.dif
## 0.04
## 0.31 -6.899852
##
## $results.new$CIU.PDE.dif
## 0.04
## 0.31 -0.2404738
##
## $results.new$TDE.dif
## 0.04
## 0.31 -3.101351
##
## $results.new$CIL.TDE.dif
## 0.04
## 0.31 -6.317415
##
## $results.new$CIU.TDE.dif
## 0.04
## 0.31 0.1468521
##
## $results.new$INT.dif
## 0.04
## 0.31 0.5120601
##
## $results.new$CIL.INT.dif
## 0.04
## 0.31 -1.300382
##
## $results.new$CIU.INT.dif
## 0.04
## 0.31 2.388272Correspondingly, we obtained the adjusted estimates of the causal effects, as shown in the right panel of the above table. A comparison of the adjusted results with the original results in the left panel of the above table reveals little influence of the additional adjustment on the estimation and inference results. Therefore, the original conclusions may be robust to the omission of baseline salary expectation from the analysis.
5. Visualization of sensitivity analysis results
The output of the modmed function and the output of the
modmed.sens function can be passed respectively via
object and sens.results to the
sens.plot function for visualizing the sensitivity analysis
results for each effect specified via effect, which must be
included in sens.effect of the modmed.sens
function.
sens.plot(object = results,
sens.results = sens.results,
effect = "PIE")
Applying the above syntax to the NEWWS data, we obtained the above
figure. The X and Y axes are respectively the slopes of U in the
standardized mediator and outcome models. The blue number on the origin
indicates the original estimate of the pure indirect effect, without
adjustment for U. Each black contour represents the combinations of
sensitivity parameters that result in the adjusted pure indirect effect
estimate as indicated by the number at the end of the contour. The
sensitivity parameters on the red dashed contours reduce the PIE
estimate to zero. Each blue dotted contour indicates the boundary at
which the significance of the PIE is changed at the confidence level
specified in the modmed function. The adjusted PIE is
insignificant in the region that contains the zero lines. The results
are more sensitive if the magnitudes of the sensitivity parameters that
change the sign or significance of PIE are smaller. The strengths of the
sensitivity parameters can be calibrated by the black dots. Each dot is
plotted based on the slopes of each observed pretreatment confounder in
the standardized mediator and outcome models and thus represents an
unmeasured confounder comparable to the observed confounder. A message
of which observed confounder each dot is comparable to is displayed
after each run of the function.
The above sensitivity analysis plot shows that an additional
adjustment of a binary unmeasured pretreatment confounder U with Pr(U =
1) = 0.5 would reverse the sign of PIE as long as the confounding role
of U is as strong as one’s baseline educational attainment
(nohsdip) or work preference (workpref).
Nevertheless, for the significance of PIE to be altered, U must be much
stronger than the strongest observed pretreatment confounder. Because we
have controlled for the most important pretreatment confounders in
theory, it is almost impossible for the remaining confounders to reverse
the significance of PIE even collectively. Hence, the sign of PIE may be
sensitive but its significance is more robust. By changing “PIE” to the
other effects, one at a time, we can obtain a sensitivity analysis plot
for each effect. The plots show that the other effects are mostly robust
to unmeasured pretreatment confounding.
In addition, we assessed sensitivity to a continuous U that follows a standard normal distribution and obtained very similar results. It is more computationally intensive because the derivation of the conditional distribution of a continuous U relies on adaptive numerical integration. Hence, if U is expected to follow a standard normal distribution, one may assess its influence by specifying U to be binary with Pr(U = 1) = 0.5 to reduce computing time.
It is worth noting that, if the treatment is not randomized and if \(\beta_u^t\) u is specified to be nonzero, the contours would cross the X and Y axes, because even though U is not associated with M (Y), it can still cause bias as long as U is associated with Y (M), given that U is associated with T. In addition, as shown in Fig. 3, the sensitivity analysis plot is symmetric with respect to the origin, indicating that the bias caused by unmeasured confounder \(U_1\) and that caused by \(U_2\) would be the same when \(\beta_{u1}^m = -\beta_{u2}^m\) and \(\beta_{u1}^y = -\beta_{u2}^y\). Similarly, when the treatment is not randomized, the bias would be the same only if \(\beta_{u1}^t = -\beta_{u2}^t\) in addition to \(\beta_{u1}^m = -\beta_{u2}^m\) and \(\beta_{u1}^y = -\beta_{u2}^y\). Therefore, at a given nonzero value of \(\beta_u^t\), the sensitivity analysis plot is asymmetric about the origin.