%\VignetteIndexEntry{mediation-old}
\documentclass[11pt,letterpaper]{article}
%% === margins ===
\addtolength{\hoffset}{-0.75in} \addtolength{\voffset}{-0.75in}
\addtolength{\textwidth}{1.5in} \addtolength{\textheight}{1.5in}
\usepackage{Sweave}
%% === basic packages ===
\usepackage{latexsym}
\usepackage{amssymb,amsmath, bm}
\usepackage{graphicx}
%% === bibliography packages ===
\usepackage{natbib}
\bibliographystyle{natbib}
%% === hyperref options ===
\usepackage{color}
\usepackage[pdftex, bookmarksopen=true, bookmarksnumbered=true,
pdfstartview=FitH, breaklinks=true, urlbordercolor={0 1 0}, citebordercolor={0 0 1}]{hyperref}
% === dcolumn package ===
\usepackage{dcolumn}
\newcolumntype{.}{D{.}{.}{-1}}
\newcolumntype{d}[1]{D{.}{.}{#1}}
% === theorem package ===
\usepackage{theorem}
\theoremstyle{plain}
\theoremheaderfont{\scshape}
\newtheorem{theorem}{Theorem}
\newtheorem{algorithm}{Algorithm}
\newtheorem{assumption}{Assumption}
\newtheorem{lemma}{Lemma}
\newtheorem{remark}{Remark}
\newcommand{\qed}{\hfill \ensuremath{\Box}}
\newcommand\indep{\protect\mathpalette{\protect\independenT}{\perp}}
\DeclareMathOperator{\sgn}{sgn}
\DeclareMathOperator{\argmin}{arg\min}
\def\independenT#1#2{\mathrel{\rlap{$#1#2$}\mkern2mu{#1#2}}}
\providecommand{\norm}[1]{\lVert#1\rVert}
% ==== rotating package ===
\usepackage{rotating}
% == spacing between sections and subsections
\usepackage[compact]{titlesec}
% == multiple figure control
\usepackage{subfigure}
\def\subfigcapskip{-.2in}
\def\subfigtopskip{-.55in}
%\def\subfigbottomskip{-1in}
%\def\subfigcapmargin{1in}
%\usepackage{times}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% === new commands ===
\newcommand\ud{\mathrm{d}}
\newcommand\dist{\buildrel\rm d\over\sim}
\newcommand\ind{\stackrel{\rm indep.}{\sim}}
\newcommand\iid{\stackrel{\rm i.i.d.}{\sim}}
\newcommand\logit{{\rm logit}}
\renewcommand\r{\right}
\renewcommand\l{\left}
\newcommand\cA{\mathcal{A}}
\newcommand\E{\mathbb{E}}
\newcommand\cJ{\mathcal{J}}
\newcommand\bR{{\bf R}}
\newcommand\bmediation{{\bf mediation}}
\newcommand\bt{\mathbf{t}}
\newcommand\bone{\mathbf{1}}
\newcommand\bzero{\mathbf{0}}
\newcommand\var{{\rm Var}}
\newcommand\cov{{\rm Cov}}
\newcommand\tomega{\tilde\omega}
\begin{document}
\newcommand\spacingset[1]{\renewcommand{\baselinestretch}%
{#1}\small\normalsize}
\spacingset{1.2}
\newcommand{\blind}{0}
\newcommand{\polisci}{1}
\newcommand{\tit}{Causal Mediation Analysis Using R}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\if0\blind
{\title{\bf \tit\thanks{This is an {\bf \color{red} old} vignette for the {\bf mediation} package,
which itself is an updated version of the tutorial originally published in
an edited \citep{imai:etal:10}. The description is based on {\bf version 3.1}
of {\bf mediation} and {\bf \color{red} does not accurately describe the most current version}.
For the most up-to-date tutorial type {\tt vignette("mediation")} to read the
new vignette. Financial support from the National Science Foundation
(SES-0752050) is acknowledged.}}
\author{Kosuke Imai\thanks{Assistant Professor, Department of
Politics, Princeton University, Princeton NJ 08544. Phone:
609--258--6610, Email:
\href{mailto:kimai@princeton.edu}{kimai@princeton.edu}, URL:
\href{http://imai.princeton.edu}{http://imai.princeton.edu}}
\quad \quad Luke Keele\thanks{Associate Professor, Department of
Political Science, 2140 Derby Hall, Ohio State University,
Columbus, OH 43210 Phone: 614-247-4256, Email:
keele.4@polisci.osu.edu}
\quad \quad Dustin Tingley\thanks{Assistant Professor, Government Department, Harvard University. Email:
\href{dtingley@gov.harvard.edu}{dtingley@gov.harvard.edu}, URL:
\href{http://scholar.harvard.edu/dtingley}{http://scholar.harvard.edu/dtingley}}
\quad \quad Teppei Yamamoto\thanks{Assistant Professor, Department of
Political Science, Massachusetts Institute of Technology, 77 Massachusetts
Avenue, Cambridge, MA 02139. Phone: 617-253-6959, Email:
\href{mailto:teppei@mit.edu}{teppei@mit.edu}, URL:
\href{http://web.mit.edu/teppei/www}{http://web.mit.edu/teppei/www}}
}
\date{ \today
% First Draft: June 13, 2009 \\
% This Draft: July 13, 2009
}
\maketitle
}\fi
\if1\blind \title{\bf \tit} %\maketitle
\fi
\pdfbookmark[1]{Title Page}{Title Page}
\thispagestyle{empty}
\setcounter{page}{0}
\begin{abstract}
Causal mediation analysis is widely used across many disciplines to
investigate possible causal mechanisms. Such an analysis allows
researchers to explore various causal pathways, going beyond the
estimation of simple causal effects. Recently,
\citet*{imai:keel:yama:10} and \citet*{imai:keel:ting:10} developed
general algorithms to estimate causal mediation effects with the
variety of data types that are often encountered in practice. The
new algorithms can estimate causal mediation effects for linear and
nonlinear relationships, with parametric and nonparametric models,
with continuous and discrete mediators, and various types of outcome
variables. In this paper, we show how to implement these algorithms
in the statistical computing language \bR. Our easy-to-use
software, \bmediation, takes advantage of the object-oriented
programming nature of the \bR\ language and allows researchers to
estimate causal mediation effects in a straightforward manner.
Finally, \bmediation\ also implements sensitivity analyses which can
be used to formally assess the robustness of findings to the
potential violations of the key identifying assumption. After
describing the basic structure of the software, we illustrate its
use with several empirical examples.
\end{abstract}
\clearpage
\section{Introduction}
Causal mediation analysis is important for quantitative social science
research because it allows researchers to identify possible causal
mechanisms, thereby going beyond the simple estimation of causal
effects. As social scientists, we are often interested in empirically
testing a theoretical explanation of a particular causal phenomenon.
This is the primary goal of causal mediation analysis. Thus, causal
mediation analysis has a potential to overcome the common criticism of
quantitative social science research that it only provides a black-box
view of causality.
Recently, \citet*{imai:keel:yama:10} and \citet*{imai:keel:ting:10}
developed general algorithms for the estimation of causal mediation
effects with a wide variety of data that are often encountered in
practice. The new algorithms can estimate causal mediation effects
for linear and nonlinear relationships, with parametric and
nonparametric models, with continuous and discrete mediators, and
various types of outcome variables.
\citet{imai:keel:yama:10,imai:keel:ting:10} also develop sensitivity
analyses which can be used to formally assess the robustness of
findings to the potential violations of the key identifying
assumption. In this paper, we describe the easy-to-use software,
\bmediation, which allows researchers to conduct causal mediation
analysis within the statistical computing language \bR\ \citep{R:09}.
We illustrate the use of the software with some of the empirical
examples presented in \citet{imai:keel:ting:10}.
\paragraph{Installation and Updating.} Before we begin, we explain how
to install and update the software. First, researchers need to
install \bR\ which is available freely at the Comprehensive R Archive
Network (\href{http://cran.r-project.org}{http://cran.r-project.org}).
Next, open \bR\ and then type the following at the prompt,
\begin{verbatim}
R> install.packages("mediation")
\end{verbatim}
Once \bmediation\ is installed, the following command will load the package,
\begin{verbatim}
R> library("mediation")
\end{verbatim}
Finally, to update \bmediation\ to its latest version, try the
following command,
\begin{verbatim}
R> update.packages("mediation")
\end{verbatim}
\section{New Features in Versions 3.x}
This section describes the new features implemented since the release
of \bmediation\ version 3.0. Although the rest of this document has
also been updated to incorporate those additional features, we discuss
them separately here for those readers who are already familiar with
previous versions of the package and only want specific information
about the recent updates. Those who are new to this package can safely
skip this section and obtain the same information elsewhere in this
vignette or in the help documents.
\subsection{Features Added in Version 3.0}
We implemented many new features in version 3.0,
along with numerous small changes to the internal code for efficiency
improvement. The user-level changes from the previous version include:
\begin{itemize}
\item More model types have been added to the main mediation analysis function {\tt mediate}.
In addition to the models already implemented in the previous version, the {\tt mediate}
function can now be used for ordered response models fitted via the {\tt polr}
function, quantile regressions for mediators fitted via the {\tt qr} function in package
{\tt quantreg}, and the tobit model for censored outcome variables fitted via the
{\tt vglm} function in package {\tt VGAM}. These models have not been made available
for sensitivity analysis via {\tt medsens}, however.
\item Continuous treatment variables can now be used with the {\tt mediate} function
via the {\tt treat.value} and {\tt control.value} arguments. The estimated average
causal mediation effects (ACME) and average direct effects (ADE) will be equal to the
effects of changing treatment from {\tt control.value} to {\tt treat.value} while
holding the appropriate treatment status constant at either of these two values.
This functionality, however, has not been implemented for the {\tt medsens} function yet.
\item The sensitivity analysis via {\tt medsens} can now be conducted with respect
to the average direct effects. Users can choose from {\tt "indirect"}, {\tt "direct"}
or {\tt "both"} using the {\tt effect.type} argument. The {\tt summary} and {\tt plot}
methods are also added for these cases.
\item A {\tt plot} method function for {\tt mediate} objects is now included,
so that the results of mediation analysis using {\tt mediate} can be easily graphically
summarized. When the input includes treatment-mediator interaction and the effect estimates
therefore differ between the treatment and control conditions, the user can choose which
set of estimates to plot by setting the {\tt treatment} argument to {\tt "treated"},
{\tt "control"}, or {\tt "both"}.
\item The {\tt summary} output for {\tt mediate} objects has been updated to include
additional information and not to include unnecessary information (e.g. separate estimates
for $\bar\delta(0)$ and $\bar\delta(1)$ when they are equal).
\item The output of the {\tt mediate} function now by default includes the entire sets of
simulation draws ({\tt d0.sims}, {\tt z0.sims}, etc.) in addition to their means
and quantiles. This can be disabled by setting the {\tt long} argument to {\tt FALSE}.
\item The {\tt mediate} function can now run even when different data frames are used
for the mediator model ({\tt model.m}) and outcome model ({\tt model.y}). This is
done by internally refitting both models using a combined data frame created via {\tt merge}.
As a result, observations which contain missing values for either model will be
automatically listwise discarded. This feature is by default disabled and can be
enabled by setting the {\tt dropobs} option to {\tt TRUE}.
\item Confidence intervals can now be calculated using heteroskedasticity-consistent
standard errors for some parametric models (e.g. {\tt lm}, {\tt glm}). This
feature is controled by the new {\tt robustSE} argument of the {\tt mediate} function.
The standard errors are obtained via the {\tt vcovHC} function in package {\tt sandwich}
and users may also pass options for {\tt vcovHC} to use different types of
robust standard errors. However, these standard errors are currently not carried
through for sensitivity analysis; the output of {\tt medsens} will not be based
on robust standard errors regardless of the {\tt robustSE} argument.
\item Users no longer need to manually indicate whether their models include a
treatment-mediator interaction. The old {\tt INT} argument is therefore deprecated
(but still accepted with a warning).
\item Sampling weights can now be used for the {\tt mediate} and {\tt medsens}
functions via the {\tt weights} arguments in the input mediator and outcome models.
The estimated quantities will thus be weighted averages when the mediator and
outcome models are fitted with a non-{\tt NULL} {\tt weights} argument. This
will be useful if, for example, data come from a complex survey.
\item We have added a new function {\tt mediations} to easily conduct multiple
causal mediation analyses for different treatment-mediator-outcome combinations.
Both {\tt summary} and {\tt plot} methods have also been added for this function.
The function is intended for exploratory analyses only and users are cautioned against
its potential abuse. See Section~\ref{sec:mediations} for more discussion.
\end{itemize}
\subsection{Features Added in Version 3.1}
We made minor updates to the \bmediation\ package upon the release of
version 3.1.
\begin{itemize}
\item Confidence intervals can now be calculated using one way clustered
standard errors\footnote{See
\href{http://people.su.se/~ma/clustering.pdf}{http://people.su.se/$\sim$ma/clustering.pdf}
for a discussion of the implementation.}
for some parametric models (e.g.,
{\tt lm}, {\tt glm}). This feature is controlled by the new {\tt cluster}
argument of the {\tt mediate} function. These standard errors are currently
not carried through for sensitivity analysis. The {\tt cluster} option is
enabled by specifying the variable on which to cluster. The default is to not
cluster and hence is set to {\tt NULL}.
\item Ordered mediator and outcome models can now be fitted via the
{\tt bayespolr} function in the {\tt arm} package, along with the
{\tt polr} function.
\end{itemize}
\section{The Software}
In this section, we give an overview of the software by describing its
design and architecture. To avoid the duplication, we do not provide
the details of the methods that are implemented by \bmediation\ and
the assumptions that underline them. Readers are encouraged to read
\citet{imai:keel:yama:10,imai:keel:ting:10} for more information about
the methodology implemented in \bmediation.
\subsection{The Overview}
The methods implemented via \bmediation\ rely on the following
identification result obtained under the sequential ignorability
assumption of \citet{imai:keel:yama:10}, {\small \begin{eqnarray}
\bar\delta(t) & = & \int \int \E(Y_i \mid M_i = m, T_i=t, X_i = x)
\l\{ dF_{M_i \mid T_i = 1, X_i = x}(m) - dF_{M_i \mid T_i = 0, X_i = x}(m)\r\}\, dF_{X_i}(x), \label{eq:acme} \\
\bar\zeta(t) & = & \int \int \l\{ \E(Y_i \mid M_i = m, T_i = 1,
X_i = x) - \E(Y_i \mid M_i = m, T_i = 0, X_i = x) \r\}\, dF_{M_i
\mid T_i = t, X_i = x}(m)\, dF_{X_i}(x), \nonumber
\\ \label{eq:direct}
\end{eqnarray}}where $\bar\delta(t)$ and $\bar\zeta(t)$ are the
average causal mediation and average (natural) direct effects, respectively, and
$(Y_i,M_i,T_i,X_i)$ represents the observed outcome, mediator,
treatment, and pre-treatment covariates. The sequential ignorability
assumption states that the observed mediator status is as if randomly
assigned conditional on the randomized treatment variable and the pre-treatment covariates. Causal mediation analysis under this assumption requires two statistical
models; one for the mediator $f(M_i \mid T_i, X_i)$ and the other for
the outcome variable $f(Y_i \mid T_i, M_i, X_i)$. (Note that we use
the empirical distribution of $X_i$ to approximate $F_{X_i}$.) Once these
models are chosen and fitted by researchers, then \bmediation\ will
compute the estimated causal mediation and other relevant estimates
using the algorithms proposed in \citet{imai:keel:ting:10}. The
algorithms also produce uncertainty estimates such as standard errors
and confidence intervals, based on either a nonparametric bootstrap procedure (for parametric or nonparametric models) or a
quasi-Bayesian Monte Carlo approximation (for parametric models).
\begin{figure}[t]
\spacingset{1}
\begin{center}
\includegraphics[scale=.35]{mediation-chart.jpg}
\end{center}
\caption{The Diagram Illustrating the Use of the Software
\bmediation. Users first fit the mediator and outcome models.
Then, the function {\tt mediate} conducts causal mediation
analysis while {\tt medsens} implements sensitivity analysis. The
functions {\tt summary} and {\tt plot} help users interpret the
results of these analyses.} \label{fg:chart}
\end{figure}
Figure~\ref{fg:chart} graphically illustrates the three steps required
for a mediation analysis. The first step is to fit the mediator and
outcome models using, for example, regression models with the usual
\texttt{lm} or \texttt{glm} functions. In the second step, the
analysts takes the output objects from these models, which in
Figure~\ref{fg:chart} we call {\tt model.m} and {\tt model.y}, and use
them as inputs for the main function, {\tt mediate}. This function
then estimates the causal mediation effects, direct effects, and total
effect along with their uncertainty estimates. Finally, sensitivity
analysis can be conducted via the function {\tt medsens} which takes
the output of {\tt mediate} as an input. For these outputs,
there are both {\tt summary} and {\tt plot} methods to
display numerical and graphical summaries of the analyses,
respectively.
\subsection{Estimation of the Causal Mediation Effects}
Estimation of the causal mediation effects is based on
Algorithms~1~and~2 of \citet{imai:keel:ting:10}. These are general
algorithms in that they can be applied to any parametric
(Algorithm~1~or~2) or semi/nonparametric models (Algorithm~2) for the
mediator and outcome variables. Here, we briefly describe how these
algorithms have been implemented in \bmediation\ by taking advantage
of the object-oriented nature of \bR\ programming language.
\paragraph{Algorithm~1 for Parametric Models.}
We begin by explaining how to implement Algorithm~1 of
\citet{imai:keel:ting:10} for standard parametric models. First,
analysts fit parametric models for the mediator and outcome
variables. That is, we model the observed mediator $M_i$ given the
treatment $T_i$ and pre-treatment covariates $X_i$. Similarly, we
model the observed outcome $Y_i$ given the treatment, mediator, and
pre-treatment covariates. For example, to implement the Baron-Kenny
procedure in \bmediation, linear models are fitted for both the
mediator and outcome models using the {\tt lm} command.
The model objects from these two parametric models form the inputs for
the {\tt mediate} function. The user must also supply the names for
the mediator and outcome variables along with how many simulations
should be used for inference, and whether the mediator variable
interacts with the treatment variable in the outcome model. Given
these model objects, the estimation proceeds by simulating the model
parameters based on their approximate asymptotic distribution (i.e.,
the multivariate normal distribution with the mean equal to the
parameter estimates and the variance equal to the asymptotic variance
estimate), and then computing causal mediation effects of interest for
each parameter draw (e.g., using
equations~\eqref{eq:acme}~and~\eqref{eq:direct} for average causal
mediation and (natural) direct effects, respectively). This method of
inference can be viewed as an approximation to the Bayesian posterior
distribution due to the Bernstein-Von Mises Theorem
\citep{king:tomz:witt:00}. The advantage of this procedure is that it
is relatively computationally efficient (when compared to
Algorithm~2).
We take advantage of the object-oriented nature of \bR\ programming
language at several steps in the function {\tt mediate}. For
example, the functions like {\tt coef} and {\tt vcov} are useful
for extracting the point and uncertainty estimates from the model
objects to form the multivariate normal distribution from which the
parameter draws are sampled. In addition, the computation of the
estimated causal mediation effects of interest requires the prediction
of the mediator values under different treatment regimes as well as
the prediction of the outcome values under different treatment and
mediator values. This can be done by using {\tt model.frame} to set
the treatment and/or mediator values to specific levels while keeping
the values of the other variables unchanged. We then use the {\tt
model.matrix} and matrix multiplication with the distribution of
simulated parameters to compute the mediation and direct effects. The
main advantage of this approach is that it is applicable to a wide
range of parametric models and allows us to avoid coding a completely
separate function for different models.
\paragraph{Algorithm~2 for Non/Semiparametric Inference.}
The disadvantage of Algorithm~1 is that it cannot be easily applied to
non/semiparametric models. For such models, Algorithm~2, which is
based on nonparametric bootstrap, can be used although it is more
computationally intensive. Algorithm~2 may also be used for the usual
parametric models which may be useful since mediation effects are
known to have skewed distributions
\citep[e.g.,][]{MacKinnon:2004,Preacher:Hayes:08}. Specifically, in
Algorithm~2, we resample the observed data with replacement. Then,
for each of bootstrapped samples, we fit both the outcome and mediator
models and compute the quantities of interest. As before, the
computation requires the prediction of the mediator values under
different treatment regimes as well as the prediction of the outcome
values under different treatment and mediator values. To take
advantage of the object-oriented nature of the \bR\ language,
Algorithm~2 relies on the \texttt{predict} function to compute these
predictions, while we again manipulate the treatment and mediator
status using the \texttt{model.frame} function. This process is
repeated a large number of times and returns a bootstrap distribution
of the mediation, direct, and total effects. We use the percentiles
of the bootstrap distribution for confidence intervals. Thus,
Algorithm~2 allows analysts to estimate mediation effects with more
flexible model specifications or to estimate mediation effects for
quantiles of the distribution.
\subsection{Sensitivity Analysis}
Causal mediation analysis relies on the sequential ignorability
assumption that cannot be directly verified with the observed data.
The assumption implies that the treatment is ignorable given the
observed pre-treatment confounders and that the mediator is ignorable
given the observed treatment and the observed pre-treatment
covariates. In order to probe the plausibility of such a key
identification assumption, analysts must perform a sensitivity
analysis \citep{rose:02c}. Unfortunately, it is difficult to
construct a sensitivity analysis that is generally applicable to any
parametric or nonparametric model. Thus,
\citet{imai:keel:yama:10,imai:keel:ting:10} develop sensitivity
analyses for commonly used parametric models, which we implement in
\bmediation.
\paragraph{The Baron-Kenny Procedure.}
\citet{imai:keel:yama:10} develop a sensitivity analysis for the
Baron-Kenny procedure and \citet{imai:keel:ting:10} generalize it to
the linear structural equation model (LSEM) with an interaction term.
This general model is given by,
\begin{eqnarray}
M_i & = & \alpha_2 + \beta_2 T_i + \xi_2^\top X_i + \epsilon_{i2}, \label{eq:MgivenTX} \\
Y_i & = & \alpha_3 + \beta_3 T_i + \gamma M_i + \kappa T_i M_i + \xi_3^\top X_i + \epsilon_{i3}, \label{eq:YgivenTMX}
\end{eqnarray}
where the sensitivity parameter is the correlation between
$\epsilon_{i2}$ and $\epsilon_{i3}$, which we denote by $\rho$. Under
sequential ignorability, $\rho$ is equal to zero and thus the
magnitude of this correlation coefficient represents the departure
from the ignorability assumption (about the mediator). Note that the
treatment is assumed to be ignorable as it would be the case in
randomized experiments where the treatment is randomized but the
mediator is not. Theorem~2 of \citet{imai:keel:ting:10} shows how the
average causal mediation effects change as a function of $\rho$.
To obtain the confidence intervals for the sensitivity analysis, we
apply the following iterative algorithm to
equations~\eqref{eq:MgivenTX}~and~\eqref{eq:YgivenTMX} for a fixed
value of $\rho$. At the $t$th iteration, given the current values of
the coefficients, i.e., $\theta^{(t)}=(\alpha_2^{(t)}, \beta_2^{(t)},
\xi_2^{(t)}, \alpha_3^{(t)},
\beta_3^{(t)},\xi_3^{(t)},\gamma^{(t)},\kappa^{(t)})$, and a given
error correlation $\rho$, we compute the variance-covariance matrix of
$(\epsilon_{i2},\epsilon_{i3})$, which is denoted by $\Sigma^{(t)}$.
The matrix is computed by setting
${\sigma_j^{(t)}}^2=\norm{\hat{\epsilon}_j^{(t)}}^2/(n-L_j)$ and
$\sigma_{23}^{(t)}=\rho \sigma_2^{(t)}\sigma_3^{(t)}$, where
$\hat\epsilon_j^{(t)}$ is the residual vector and $L_j$ is the number
of coefficients for the mediator model ($j=2$) and the outcome model
($j=3$) at the $t$th iteration. We then update the parameters via
generalized least squares, i.e.,
\begin{eqnarray*}
\theta^{(t+1)} & = & \{V^\top ({\Sigma^{(t)}}^{-1} \otimes I_n)V\}^{-1}
V^\top ({\Sigma^{(t)}}^{-1} \otimes I_n)W
\end{eqnarray*}
where
$V = \l[\begin{array}{cccccccc} \mathbf{1} & T & X & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\
\mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{1} & T & M & TM & X
\end{array} \r]$,
$W \ = \ \l[\begin{array}{c} M \\
Y
\end{array} \r]$, $T=(T_1,\dots,T_n)^\top$, $M=(M_1,\dots,M_n)^\top$, and
$Y=(Y_1,\dots,Y_n)^\top$ are column vectors of length $n$, and $X = (X_1, \dots, X_n)^\top$ are the $(n \times K)$ matrix of observed pre-treatment covariates, and $\otimes$ represents the Kronecker product. We typically use
equation-by-equation least squares estimates as the starting values of $\theta$
and iterate these two steps until convergence. This is essentially an
application of the iterative feasible generalized least square (FGLS) algorithm of the
seemingly unrelated regression \citep{zell:62}, and thus the asymptotic variance
of $\hat\theta$ is given by $\var(\hat\theta) = \{V^\top (\Sigma^{-1} \otimes
I_n) V\}^{-1}$. Then, for a given value of $\rho$, the asymptotic variance of the estimated average causal mediation effects
is found, for example, by the Delta method and the confidence intervals can be constructed.
\paragraph{The Binary Outcome Case.}
The sensitivity analysis for binary outcomes parallels the case when
both the mediator and outcome are continuous. Here, we assume that
the model for the outcome is a probit regression. Using a probit
regression for the outcome allows us to assume the error terms are
jointly normal with a possibly non-zero correlation $\rho$.
\citet{imai:keel:ting:10} derive the average causal mediation effects
as a function of $\rho$ and a set of parameters that are identifiable
due to randomization of the treatment. This lets us use $\rho$ as a
sensitivity parameter in the same way as in the Baron-Kenny procedure.
For the calculation of confidence intervals, we rely on the
quasi-Bayesian approach of Algorithm~1 by approximating the posterior
distribution with the sampling distribution of the maximum likelihood
estimates.
% While it is possible to use a nonparametric bootstrap for the
% confidence intervals, it is much more computationally intensive.
% For a graphical summary of sensitivity analysis, we use a lowess
% smoother to smooth out sampling variability from the simulation when
% drawing the lines representing the point estimates and confidence
% intervals at a variety of $\rho$ values.
\paragraph{The Binary Mediator Case.}
Finally, a similar sensitivity analysis can also be conducted in a
situation where the mediator variable is dichotomous and the outcome
is continuous. In this case, we assume that the mediator can be
modeled as a probit regression where the error term is independently
and identically distributed as standard normal distribution. A linear
normal regression with error variance equal to $\sigma_3^2$ is used to
model the continuous outcome variable. We further assume that the two
error terms jointly follow a bivariate normal distribution with mean
zero and covariance $\rho\sigma_3$. Then, as in the other two cases,
we use the correlation between the two error terms $\rho$ as the
sensitivity parameter. \citet{imai:keel:ting:10} show that under this
setup, the causal mediation effects can be expressed as a function of
the model parameters that can be consistently estimated given a fixed
value of $\rho$. Uncertainty estimates are computed based on the
quasi-Bayesian approach, as in the binary outcome case. The results
can be graphically summarized via the \texttt{plot} function in a
manner similar to the other two cases.
\paragraph{Alternative Interpretations Based on $R^2$}
The main advantage of using $\rho$ as a sensitivity parameter is its
simplicity. However, applied researchers may find it difficult to
interpret the magnitude of this correlation coefficient. To overcome
this limitation, \citet{imai:keel:yama:10} proposed alternative
interpretations of $\rho$ based on the coefficients of determination
or $R^2$ and \citet{imai:keel:ting:10} extended them to the binary
mediator and binary outcome cases. In that formulation, it is assumed
that there exists a common unobserved pre-treament confounder in both
mediator and outcome models. Applied researchers are then required to
specify whether the coefficients of this unobserved confounder in the
two models have the same sign or not; i.e.,
$\sgn(\lambda_2\lambda_3)=1$ or $-1$ where $\lambda_2$ and $\lambda_3$
are the coefficients in the mediator and outcome models, respectively.
Once this information is provided, the average causal mediation effect
can be expressed as the function of ``the proportions of original
variances explained by the unobserved confounder'' where the original
variances refer to the variances of the mediator and the outcome (or
the variance of latent variable in the case of binary dependent
variable). This allows researchers to quantify how large the
unobserved confounder must be (relative to the observed pre-treatment
covariates in the model) in order for the original conclusions to be
reversed.
\paragraph{Sensitivity Analyses with Respect to the Average Direct Effects.}
In addition to the above procedures for the average causal mediation effects,
an analogous sensitivity analysis can also be conducted for the average direct
effects using the {\tt medsens} function. This is possible because it can be
shown that the average direct effects are also expressed as a function of $\rho$
in any of the above three cases. First, when both the mediator and outcome variables
are continuous and modelled using linear structural equations, the iterative FGLS algorithm
given above can be used to estimate the model coefficients and both the point and interval
estimates of the average direct effects can be obtained using the formulas given by
\citet[Section 4.1]{imai:keel:yama:10} and \citet[p.314 and Appendix A]{imai:keel:ting:10}.
Second, when the mediator is binary and modelled using the probit regression,
the average direct effect is given by $\bar\zeta(t) = \beta_3 + \kappa\E\{
\Phi(\alpha_2 + \beta_2 t + \xi_2^\top X_i)\}$ using the same notation as
\citet[footnote 12]{imai:keel:ting:10}, where all of the parameters, $(\alpha_2,
\beta_2, \xi_2, \beta_3, \kappa)$, are identified and consistently estimated
using the same procedure. Finally, for the binary outcome case,
the average direct effect is identified given a value of $\rho$ and expressed
as $\bar\zeta(t) = \E\{\Phi(\alpha_1 + \beta_1 t + \xi_1^\top X_i +
\beta_3(1-t)/\sqrt{\gamma^2\sigma_2^2 + 2\gamma\rho\sigma_2 + 1})
- \Phi(\alpha_1 + \beta_1 t + \xi_1^\top X_i - \beta_3 t/\sqrt{\gamma^2\sigma_2^2 + 2\gamma\rho\sigma_2 + 1}) \}$, which can be estimated via the same procedure as
described in Appendix H of \citet{imai:keel:ting:10}. The {\tt medsens} function
uses these procedures for sensitivity analyses for the average direct effects.
\subsection{Current Limitations}
\begin{table}[t]
\spacingset{1}
\begin{center}
\begin{tabular}{lcccccc}
\hline
&\multicolumn{6}{c}{\it Outcome Model Types} \\
\cline{2-7}
{\it Mediator Model Types} & Linear & GLM & Ordered & Censored & Quantile & Semiparametric \\
\hline
Linear ({\tt lm}) & \checkmark & \checkmark & \checkmark$^\ast$ & \checkmark &
\checkmark & \checkmark$^\ast$ \\
GLM (probit etc. via {\tt glm}) & \checkmark & \checkmark & \checkmark$^\ast$ & \checkmark &
\checkmark & \checkmark$^\ast$ \\
Ordered ({\tt polr}/{\tt bayespolr}) & \checkmark & \checkmark & \checkmark$^\ast$ & \checkmark &
\checkmark & \checkmark$^\ast$ \\
Censored (tobit via {\tt vglm}) & - & - & - & - & - & - \\
Quantile ({\tt rq}) & \checkmark$^\ast$ & \checkmark$^\ast$ & \checkmark$^\ast$ & \checkmark$^\ast$ &
\checkmark$^\ast$ & \checkmark$^\ast$ \\
Semiparametric ({\tt gam}) & \checkmark$^\ast$ & \checkmark$^\ast$ & \checkmark$^\ast$ & \checkmark$^\ast$ &
\checkmark$^\ast$ & \checkmark$^\ast$ \\
\hline
\end{tabular}
\caption{The Types of Models That Can be Handled by {\tt mediate} for the
Estimation of Causal Mediation Effects. Stars ($^\ast$) indicate the model combinations
that can only be estimated using the nonparametric bootstrap (i.e. with {\tt boot = TRUE}).}
\label{tab:MediateOptions}
\end{center}
\end{table}
Our software, \bmediation, is quite flexible and can handle many of
the model types that researchers are likely to use in practice.
Table~\ref{tab:MediateOptions} categorizes the types of the mediator
and outcome models and lists whether \bmediation\ can produce the
point and uncertainty estimates of causal mediation effects. For
example, while \bmediation\ can estimate average causal mediation
effects for censored outcomes with the tobit model, it
has not yet been extended to cases involving censored mediators.
In addition, {\tt mediate} can be used with common generalized linear models
(GLMs) such as binary logit and probit, Poisson and gamma regression models,
but it cannot be used with binomial mediator models with multiple trials.
In each situation handled by \bmediation, it is possible to have an
interaction term between treatment status and the mediator variable,
in which case the estimated quantities of interest will be reported
separately for the treatment and control groups.
\begin{table}[t]
\spacingset{1}
\begin{center}
\begin{tabular}{lcc}
\hline
&\multicolumn{2}{c}{\it Outcome Model Types} \\
\cline{2-3}
{\it Mediator Model Types} & Linear & Binary Probit \\
\hline
Linear & \checkmark & \checkmark \\
Binary Probit & \checkmark & - \\
\hline
\end{tabular}
\caption{The Types of Models That Can Be Handled by {\tt medsens} for
Sensitivity Analysis. } \label{tab:SensitivityOptions}
\end{center}
\end{table}
Our software provides a convenient way to probe the sensitivity of
results to potential violations of the ignorability assumption for
certain model types. The sensitivity analysis requires the specific
derivations for each combination of models, making it difficult to
develop a general sensitivity analysis method. As summarized in
Table~\ref{tab:SensitivityOptions}, \bmediation\ can handle several
cases that are likely to be encountered by applied researchers in
practice. When the mediator is continuous, then sensitivity analysis
can be conducted with continuous and binary outcome variables. In
addition, when the mediator is binary, sensitivity analysis is
available for continuous outcome variables. For sensitivity analyses
that combine binary or continuous mediators and outcomes, analysts
must use a probit regression model with a linear regression model.
It should be noted that unlike the
estimation of causal mediation effects, sensitivity analysis with
treatment/mediator interactions can only be done for the continuous
outcome/continuous mediator and continuous outcome/binary mediator cases.
In the future, we hope to expand the range of models that
are available for sensitivity analysis.
\section{Examples}
Next, we provide several examples to illustrate the use of
\bmediation\ for the estimation of causal mediation effects and
sensitivity analysis. The data used is available as part of the
package so that readers can replicate the results reported below. We
demonstrate the variety of models that can be used for the outcome and
mediating variables.
Before presenting our examples, we load the \bmediation\ library and
the example data set included with the library.
\begin{verbatim}
R> library("mediation")
mediation: R Package for Causal Mediation Analysis
Version: 3.0
R> data("jobs")
\end{verbatim}
This dataset is from the Job Search Intervention Study (JOBS II)
\citep{Vinokur:1997}. In the JOBS II field experiment, 1,801
unemployed workers received a pre-screening questionnaire and were
then randomly assigned to treatment and control groups. Those in the
treatment group participated in job-skills workshops. Those in the
control condition received a booklet describing job-search tips. In
follow-up interviews, two key outcome variables were measured; a
continuous measure of depressive symptoms based on the Hopkins Symptom
Checklist (\texttt{depress2}), and a binary variable representing
whether the respondent had become employed (\texttt{work1}). In the
JOBS II data, a continuous measure of job-search self-efficacy
represents a key mediating variable (\texttt{job\_seek}). In addition
to the outcome and mediators, the JOBS II data also include the
following list of baseline covariates that were measured prior to the
administration of the treatment: pre-treatment level of depression
(\texttt{depress1}), education (\texttt{educ}), income, race
(\texttt{nonwhite}), marital status (\texttt{marital}), age, sex,
previous occupation (\texttt{occp}), and the level of economic
hardship (\texttt{econ\_hard}).
\subsection{Estimation of Causal Mediation Effects}
\paragraph{The Baron-Kenny Procedure.}
We start with an example when both the mediator and the outcome are
continuous. In this instance, the results from either algorithm will
return point estimates essentially identical to the usual Baron and
Kenny procedure though the quasi-Bayesian or nonparametric bootstrap
approximation is used. Using the JOBS II data, we first estimate two
linear regressions for both the mediator and the outcome using the {\tt
lm} function.
\begin{verbatim}
R> model.m <- lm(job_seek ~ treat + depress1 + econ_hard + sex + age
+ occp + marital + nonwhite + educ + income, data = jobs)
R> model.y <- lm(depress2 ~ treat + job_seek + depress1 + econ_hard + sex + age
+ occp + marital + nonwhite + educ + income, data = jobs)
\end{verbatim}
These two model objects, {\tt model.m} and {\tt model.y}, become the
arguments for the \texttt{mediate} function. When different data frames
are used for these two models, the {\tt mediate} function internally refits
them using a merged dataset if the {\tt dropobs} argument is set to
{\tt TRUE}; otherwise an error message is returned.
In the first call to \texttt{mediate} below, we specify \texttt{boot
= TRUE} to call the nonparametric bootstrap with 1000 resamples
({\tt sims = 1000}). When this option is set to {\tt FALSE} in the second
call, inference proceeds via the quasi-Bayesian Monte Carlo
approximation using Algorithm~1 rather than Algorithm~2. In this case,
heteroskedasticity-robust standard errors can be used for the calculation of
confidence intervals for some types of models by setting the {\tt robustSE}
argument to {\tt TRUE}. Cluster-robust standard errors can also be
used by the {\tt cluster} argument. Finally, we must also
specify the variable names for the treatment indicator and the
mediator variable using {\tt treat} and {\tt mediator}, respectively.
\begin{verbatim}
R> out.1 <- mediate(model.m, model.y, sims = 1000, boot = TRUE, treat = "treat",
mediator = "job_seek")
R> out.2 <- mediate(model.m, model.y, sims = 1000, treat = "treat",
mediator = "job_seek")
\end{verbatim}
The objects from a call to \texttt{mediate}, i.e., {\tt out.1} and
{\tt out.2} above, are lists of class {\tt "mediate"} which contain several different
quantities from the analysis. For example, \texttt{out.1\$d0} returns
the point estimate for the average causal mediation effect based on
Algorithm~1. The help file contains a full list of values that are
contained in \texttt{mediate} objects.
The \texttt{summary}
function prints out the results of the analysis in tabular form:
\begin{verbatim}
R> summary(out.2)
Causal Mediation Analysis
Quasi-Bayesian Confidence Intervals
Mediation Effect: -0.01375 95 % CI -0.031812 0.003299
Direct Effect: -0.03873 95 % CI -0.11938 0.04668
Total Effect: -0.05247 95 % CI -0.13431 0.03150
Proportion of Total Effect via Mediation: 0.2028 95 % CI -1.383 1.934
Sample Size Used: 899
\end{verbatim}
The output from the {\tt summary} function displays the estimates
for the average causal mediation effect, direct effect, total effect,
and proportion of total effect mediated. The first column displays the
quantity of interest, the second column displays the point estimate,
and the other columns present the confidence intervals at the level
specified by the {\tt conf.level} argument which defaults to {\tt 0.95}.
The number of observations used for the estimation is then displayed at
the bottom. In this case, we
find that job search self-efficacy mediated the effect of the
treatment on depression in the negative direction. This effect,
however, was small with a point estimate of $-.014$ and the $95\%$
confidence interval $(-.032,-.003)$ barely covers $0$.
These results can be graphically presented by using the {\tt plot} function:
\begin{verbatim}
R> plot(out.2)
\end{verbatim}
The output is shown in Figure~\ref{fig:plot.mediate}. The standard arguments
for {\tt plot}, such as titles ({\tt main}, {\tt xlab}, etc.) and line widths
({\tt lwd}), can be used and behave as expected.
\begin{figure}[t]
\spacingset{1}
\vspace{-.5in}
\begin{center}
\includegraphics[scale=.8]{plot-mediate.pdf}
\end{center}
\vspace{-.5in}
\caption{Graphical Summary of Causal Mediation Analysis.
\label{fig:plot.mediate}}
\end{figure}
\paragraph{The Baron-Kenny Procedure with the Interaction Term.}
Analysts can also allow the causal mediation effect to vary with
treatment status. Here, the model for the outcome must be altered by
including an interaction term between the treatment indicator,
\texttt{treat}, and the mediator variable, \texttt{job\_seek}:
\begin{verbatim}
R> model.y <- lm(depress2 ~ treat + job_seek + treat:job_seek + depress1 + econ_hard
+ sex + age + occp + marital + nonwhite + educ + income, data = jobs)
R> out.3 <- mediate(model.m, model.y, sims = 1000, boot = TRUE,
treat = "treat", mediator = "job_seek")
R> summary(out.3)
Causal Mediation Analysis
Quasi-Bayesian Confidence Intervals
Mediation Effect_0: -0.01899 95 % CI -0.045271 0.003649
Mediation Effect_1: -0.01203 95 % CI -0.029640 0.002195
Direct Effect_0: -0.03596 95 % CI -0.1187 0.0426
Direct Effect_1: -0.02899 95 % CI -0.10879 0.05082
Total Effect: -0.04799 95 % CI -0.13495 0.03154
Proportion of Total Effect via Mediation_0: 0.2912 95 % CI -3.386 2.707
Proportion of Total Effect via Mediation_1: 0.1823 95 % CI -1.994 2.018
Mediation Effect (Average): -0.01551 95 % CI -0.035757 0.003138
Direct Effect (Average): -0.03248 95 % CI -0.1129 0.0464
Proportion of Total Effect via Mediation (Average): 0.2381 95 % CI -2.867 2.443
Sample Size Used: 899
\end{verbatim}
Again using the \texttt{summary} function provides a table of the
results. Now estimates for the mediation effects, direct effects and
proportion of totel effect mediated correspond
to the levels of the treatment and are printed as such in the tabular
summary. In this case, the mediation effect under the treatment
condition, listed as \texttt{Mediation Effect\_1} is estimated to be
$-.012$ while the mediation effect under the control condition,
\texttt{Mediation Effect\_0}, is $-.019$. The bottom of the table shows
the overall summary of these quantities averaged over the two treatment levels,
as well as the number of data points.
When using the {\tt plot} function on a {\tt mediate} object with interactions,
one can select which treatment condition to plot the estimated effects for by
specifying the {\tt treatment} argument. The default is to plot both estimates,
which is equivalent to setting {\tt treatment} to {\tt "both"}:
\begin{verbatim}
R> plot(out.3, treatment = "both")
\end{verbatim}
The resulting plot is in Figure~\ref{fig:plot.mediate.int}.
\begin{figure}[t]
\spacingset{1}
\vspace{-.5in}
\begin{center}
\includegraphics[scale=.8]{plot-mediate-int.pdf}
\end{center}
\vspace{-.5in}
\caption{Causal Mediation Analysis with Interaction
between Treatment and Mediator.
\label{fig:plot.mediate.int}}
\end{figure}
\paragraph{Use of Non/Semi-parametric Regression.}
The flexibility of \bmediation\ becomes readily apparent when we move
beyond standard linear regression models. For example, we might
suspect that the mediator has a nonlinear effect on the outcome.
Generalized Additive Models (GAMs) allow analysts to use splines for
flexible nonlinear fits. This presents no difficulties for the
\texttt{mediate} function. We model the mediator as before, but we
alter the outcome model using the \texttt{gam} function from the
\texttt{mgcv} library.
\begin{verbatim}
R> library(mgcv)
R> model.m <- lm(job_seek ~ treat + depress1 + econ_hard + sex + age + occp + marital
+ nonwhite + educ + income, data = jobs)
R> model.y <- gam(depress2 ~ treat + s(job_seek, bs = "cr") + depress1 + econ_hard
+ sex + age + occp + marital + nonwhite + educ + income, data = jobs)
\end{verbatim}
In this case we fit a Generalized Additive Model for the outcome
variable, but allow the effect of the \texttt{job\_seek} variable to
be non-linear and determined by the data. This is done by using the
\texttt{s} notation which allows the fit between the mediator and
the outcome to be modeled with a spline. Using the spline for the fit
allows the estimate for the mediator on the outcome to be a series of
piecewise polynomial regression fits. This semiparametric regression
model is a more general version of nonparametric regression models
such as lowess. The model above allows the estimate to vary across the
range of the predictor variable. Here, we specify the model with a
cubic basis function (\texttt{bs = "cr"}) for the smoothing spline and
leave the smoothing selection to be done at the program defaults which
is generalized cross-validation. Fully understanding how to fit such
models is beyond the scope here. Interested readers should consult
\cite{Wood:2006} for full technical details and \cite{Keele:2008}
provides coverage of these models from a social science perspective.
The call to \texttt{mediate} with a \texttt{gam} fit remains
unchanged except that when the outcome model is a semiparametric
regression only the nonparametric bootstrap is valid for calculating
uncertainty estimates, i.e., {\tt boot = TRUE}.
\begin{verbatim}
R> out.4 <- mediate(model.m, model.y, sims = 1000, boot = TRUE, treat = "treat",
mediator = "job_seek")
\end{verbatim}
The model for the mediator can be modeled with the \texttt{gam}
function as well. The \texttt{gam} function also allows analysts to
include interactions; thus analysts can still allow the mediation
effects to vary with treatment status. This simply requires altering
the model specification by using the \texttt{by} option in the
\texttt{gam} function and using two separate indicator variables for
treatment status. To fit this model we need one variable that
indicates whether the observation was in the treatment group and a
second variable that indicates whether the observation was in the
control group. To allow the mediation effect to vary with treatment
status, the call to \texttt{gam} takes the following form:
\begin{verbatim}
R> model.y <- gam(depress2 ~ treat + s(job_seek, by = treat) + s(job_seek, by = control)
+ depress1 + econ_hard + sex + age + occp + marital + nonwhite + educ + income,
data = jobs)
\end{verbatim}
In this case, we must also alter the options in the \texttt{mediate}
function by providing the variable
name for the control group indicator using the \texttt{control} option.
\begin{verbatim}
R> out.5 <- mediate(model.m, model.y, sims = 1000, boot = TRUE,
treat = "treat", mediator = "job_seek", control = "control")
R> summary(out.5)
Causal Mediation Analysis
Confidence Intervals Based on Nonparametric Bootstrap
Mediation Effect_0: -0.02254 95 % CI -0.051956 0.004447
Mediation Effect_1: -0.0171 95 % CI -0.038895 0.002281
Direct Effect_0: -0.009742 95 % CI -0.09718 0.07924
Direct Effect_1: -0.004302 95 % CI -0.08926 0.08285
Total Effect: -0.02684 95 % CI -0.11932 0.06449
Proportion of Total Effect via Mediation_0: 0.3892 95 % CI -4.411 9.161
Proportion of Total Effect via Mediation_1: 0.2768 95 % CI -3.614 6.960
Mediation Effect (Average): -0.01982 95 % CI -0.045831 0.003105
Direct Effect (Average): -0.007022 95 % CI -0.09207 0.08178
Proportion of Total Effect via Mediation (Average): 0.3398 95 % CI -3.870 8.105
Sample Size Used: 899
\end{verbatim}
\paragraph{Quantile Causal Mediation Effects.}
Researchers might also be interested in modeling mediation effects for
quantiles of the outcome. Quantile regression allows for a convenient
way to model the quantiles of the outcome distribution while adjusting
for a variety of covariates \citep{Koenker:2005}. For example, a
researcher might be interested in the 0.5 quantile (i.e., median)
of the distribution. This also presents no difficulties for the
\texttt{mediate} function. Again for these models, uncertainty
estimates are calculated using the nonparametric bootstrap. To use
quantile regression, we use the {\tt rq} function in package
\texttt{quantreg} and model
the median of the outcome, though other quantiles are also
permissible. Analysts can also relax the no-interaction assumption
for the quantile regression models as well. Below we estimate the
mediator with a standard linear regression, while for the outcome we
use \texttt{rq} to model the median.
\begin{verbatim}
R> library(quantreg)
R> model.m <- lm(job_seek ~ treat + depress1 + econ_hard + sex + age + occp + marital
+ nonwhite + educ + income, data = jobs)
R> model.y <- rq(depress2 ~ treat + job_seek + depress1 + econ_hard + sex
+ age + occp + marital + nonwhite + educ + income,
tau = 0.5, data = jobs)
R> out.6 <- mediate(model.m, model.y, sims = 1000, boot = TRUE, treat = "treat",
mediator = "job_seek")
R> summary(out.6)
Causal Mediation Analysis
Confidence Intervals Based on Nonparametric Bootstrap
Mediation Effect: -0.01424 95 % CI -0.031317 0.002109
Direct Effect: -0.04165 95 % CI -0.12413 0.03538
Total Effect: -0.05589 95 % CI -0.14305 0.02377
Proportion of Total Effect via Mediation: 0.2042 95 % CI -2.003 2.707
Sample Size Used: 899
\end{verbatim}
where the {\tt summary} command gives the estimated median causal
mediation effect along with the estimates for the other quantities of
interest.
It is also possible to estimate mediation effects for quantiles of the
outcome other than the median. This is done simply by specifying a
different outcome quantile in the quantile regression model. For
example, if the $10th$ percentile of the outcome were of interest,
then the user can change {\tt tau} option,
\begin{verbatim}
R> model.y <- rq(depress2 ~ treat + job_seek + depress1 + econ_hard + sex
+ age + occp + marital + nonwhite + educ + income,
tau = 0.1, data = jobs)
\end{verbatim}
Furthermore, it is straightforward to loop over a set of quantiles and
graph the mediation effects for a range of quantiles, as done in
\citet{imai:keel:ting:10}.
\paragraph{Discrete Mediator and Outcome Data.}
Often analysts use measures for the mediator and outcome that are
discrete. For standard methods, this has presented a number of
complications requiring individually tailored techniques. The
\bmediation\ software, however, can handle a number of different
discrete data types using the general algorithms developed in
\citet{imai:keel:ting:10}. For example, one outcome of interest in
the JOBS II study is a binary indicator (\texttt{work1}) for whether
the subject became employed after the training sessions. To estimate
the mediation effect, we simply use a probit regression instead of a
linear regression for the outcome and then call \texttt{mediate} as
before,
\begin{verbatim}
R> model.m <- lm(job_seek ~ treat + depress1 + econ_hard + sex + age
+ occp + marital + nonwhite + educ + income, data = jobs)
R> model.y <- glm(work1 ~ treat + job_seek + depress1 + econ_hard + sex + age
+ occp + marital + nonwhite + educ + income,
family = binomial(link = "probit"), data = jobs)
R> out.7 <- mediate(model.m, model.y, sims = 1000, treat = "treat",
mediator = "job_seek")
R> summary(out.7)
Causal Mediation Analysis
Quasi-Bayesian Confidence Intervals
Mediation Effect_0: 0.003290 95 % CI -0.0009537 0.0098747
Mediation Effect_1: 0.003538 95 % CI -0.001052 0.010728
Direct Effect_0: 0.05514 95 % CI -0.006427 0.113568
Direct Effect_1: 0.05538 95 % CI -0.006487 0.114191
Total Effect: 0.05867 95 % CI -0.003816 0.117091
Proportion of Total Effect via Mediation_0: 0.0455 95 % CI -0.0768 0.4070
Proportion of Total Effect via Mediation_1: 0.05042 95 % CI -0.0777 0.4110
Mediation Effect (Average): 0.003414 95 % CI -0.001003 0.010269
Direct Effect (Average): 0.05526 95 % CI -0.006457 0.113879
Proportion of Total Effect via Mediation (Average): 0.04802 95 % CI -0.07725 0.41005
Sample Size Used: 899
\end{verbatim}
In the table printed by the \texttt{summary} function, the estimated
average causal mediation effect along with the quasi-Bayesian
confidence interval are printed on the first line followed by the
direct and total effects, and the proportion of the total effect due
to mediation. Even though the outcome model does not include an interaction
term between the treatment and mediator, the estimated effects slightly differ
between the treatment and control conditions. This difference, however,
is solely due to the nonlinearity in the outcome model and should be
extremely small. It is also possible to use a logit model for the
outcome instead of a probit model. However, we recommend the use of a
probit model because our implementation of the sensitivity analysis
below requires a probit model for analytical tractability.
The mediator can also be binary or ordered measures as well. This
simply requires modeling the mediator with either a probit or ordered
probit model. For demonstration purposes, the \texttt{jobs} data
contains two variables, \texttt{job\_dich} and \texttt{job\_disc},
which are recoded versions of \texttt{job\_seek}. The first measure
is simply the continuous scale divided at the median into a binary
variable. The second measure, \texttt{job\_disc}, recodes the
continuous scale into a discrete four point scale. We emphasize that
this is for demonstration purposes only, and analysts in general
should not recode continuous measures into discrete measures.
Estimating mediation effects with a binary mediator is quite similar
to the case above with a binary outcome. We simply now use a probit
model for the mediator and a linear regression for the outcome,
\begin{verbatim}
R> model.m <- glm(job_dich ~ treat + depress1 + econ_hard + sex + age
+ occp + marital + nonwhite + educ + income, data = job,
family = binomial(link = "probit"))
R> model.y <- lm(depress2 ~ treat + job_dich + treat:job_dich + depress1
+ econ_hard + sex + age + occp + marital + nonwhite
+ educ + income, data = jobs)
\end{verbatim}
In this example we include the interaction term between the treatment and
mediator to allow the effect of the mediator to vary with
treatment status. The user now calls \texttt{mediate} and can use
either the quasi-Bayesian approximation or nonparametric bootstrap.
\begin{verbatim}
R> out.8 <- mediate(model.m, model.y, sims = 1000, treat = "treat",
mediator = "job_dich")
R> summary(out.8)
Causal Mediation Analysis
Quasi-Bayesian Confidence Intervals
Mediation Effect_0: -0.01823 95 % CI -0.038944 -0.002616
Mediation Effect_1: -0.02040 95 % CI -0.038965 -0.003155
Direct Effect_0: -0.02823 95 % CI -0.10786 0.05154
Direct Effect_1: -0.03039 95 % CI -0.11164 0.04764
Total Effect: -0.04863 95 % CI -0.13146 0.02839
Proportion of Total Effect via Mediation_0: 0.2782 95 % CI -3.055 3.325
Proportion of Total Effect via Mediation_1: 0.3162 95 % CI -3.443 3.628
Mediation Effect (Average): -0.01932 95 % CI -0.037878 -0.003201
Direct Effect (Average): -0.02931 95 % CI -0.11052 0.04807
Proportion of Total Effect via Mediation (Average): 0.301 95 % CI -3.016 3.477
Sample Size Used: 899
\end{verbatim}
\noindent In the table, we see that \texttt{Mediation Effect\_0} is
the mediation effect under the control condition, while
\texttt{Mediation Effect\_1} is the mediation effect under the
treatment condition. The same notation applies to the direct effects.
As the reader can see, the output also indicates which algorithm is
used for the 95\% confidence intervals.
When the mediator is an ordered variable, we switch to an ordered
probit model for the mediator. In \bR, the \texttt{polr} function
in the \texttt{MASS} library provides this functionality. (Alternatively,
the \texttt{bayespolr} function in the \texttt{arm} library can also be used.)
Thus, we fit the outcome and mediator models below,
\begin{verbatim}
R> model.m <- polr(job_disc ~ treat + depress1 + econ_hard + sex + age
+ occp + marital + nonwhite + educ + income,
data = jobs, method = "probit", Hess = TRUE)
R> model.y <- lm(depress2 ~ treat + job_disc + depress1 + econ_hard + sex + age
+ occp + marital + nonwhite + educ + income, data = jobs)
\end{verbatim}
The reader should note that in the call to \texttt{polr} the
\texttt{Hess = TRUE} should be specified to compute the asymptotic variance
covariance matrix for the estimated model coefficients in the quasi-Bayesian
approximation. Otherwise, the Hessian of the log-likelihood will be recalculated
within the {\tt mediate} call.
Once we have estimated these two models, analysis proceeds as before,
\begin{verbatim}
R> out.9 <- mediate(model.m, model.y, sims = 1000, treat = "treat",
mediator = "job_disc")
R> summary(out.9)
Causal Mediation Analysis
Quasi-Bayesian Confidence Intervals
Mediation Effect: -0.01048 95 % CI -0.03188 0.01008
Direct Effect: -0.03311 95 % CI -0.10759 0.04848
Total Effect: -0.04360 95 % CI -0.12376 0.03895
Proportion of Total Effect via Mediation: 0.1726 95 % CI -2.752 2.077
Sample Size Used: 899
\end{verbatim}
Again, for any of these data types, analysts can relax the
no-interaction assumption as before by including the interaction
between treatment and the mediator variable in the outcome model.
\subsection{Sensitivity Analysis}
Once analysts have estimated mediation effects, they should explore how
robust their finding is to the ignorability assumption, whenever possible.
The \texttt{medsens} function allows analysts to conduct sensitivity
analyses for mediation effects as well as direct effects.
In this section we provide a demonstration of these functionalities.
Currently, \bmediation\ can conduct sensitivity analyses when both mediator
and outcome are continuous and modelled with the {\tt lm} function, as well
as when either mediator or outcome (but not both) is binary and modelled using
a binary probit regression via {\tt glm}.
\paragraph{The Baron-Kenny Procedure.}
As before, one must first fit models for the mediator and outcome and then pass these model objects through the mediate function:
\begin{verbatim}
R> model.m <- lm(job_seek ~ treat + depress1 + econ_hard + sex + age + occp
+ marital + nonwhite + educ + income, data = jobs)
R> model.y <- lm(depress2 ~ treat + job_seek + depress1 + econ_hard + sex + age + occp
+ marital + nonwhite + educ + income, data = jobs)
R> med.cont <- mediate(model.m, model.y, sims=1000, treat = "treat",
mediator = "job_seek")
\end{verbatim}
Once the analyst estimates the mediation effects, the output from the
{\tt mediate} function becomes the argument for \texttt{medsens},
which is the workhorse function. The {\tt medsens} function
recognizes main arguments specified in the {\tt mediate} function and
thus there is no need to specify options such as \texttt{treat} or
\texttt{mediator}.
\begin{verbatim}
R> sens.cont <- medsens(med.cont, rho.by = 0.05)
\end{verbatim}
The \texttt{rho.by} option specifies how finely incremented the parameter $\rho$
is for the sensitivity analysis. When $\rho$ is on a coarser grid, this speeds
up estimation considerably but comes at the cost of estimating the
robustness of the original conclusion only imprecisely.
After running the sensitivity analysis via \texttt{medsens}, the
\texttt{summary} function can be used to produce a table with the
values of $\rho$ for which the confidence interval contains zero.
This allows the analyst to immediately see the approximate range of
$\rho$ where the sign of the causal mediation effect is indeterminate.
The second section of the table contains the value of $\rho$ for which
the mediation effect is exactly zero, which in this application is
$-0.19$. The table also presents coefficients of determination that
correspond to the critical value of $\rho$ where the mediation effect
is zero. First, ${R^\ast}^2_M{R^\ast}^2_Y$ is the product of
coefficients of determination which represents the proportion of the
\emph{previously unexplained} variance in the mediator and outcome
variables that is explained by an unobservable pretreatment
unconfounder. An alternative formulation is in terms of the
proportion of the \emph{original} variance explained by an unobserved
confounder, which we denote as $\widetilde{R}^2_M\widetilde{R}^2_Y$ .
\begin{verbatim}
R> summary(sens.cont)
Mediation Sensitivity Analysis for Average Causal Mediation Effect
Sensitivity Region
Rho Med. Eff. 95% CI Lower 95% CI Upper R^2_M*R^2_Y* R^2_M~R^2_Y~
[14,] -0.30 0.0064 -0.0025 0.0153 0.0900 0.0588
[15,] -0.25 0.0028 -0.0026 0.0082 0.0625 0.0408
[16,] -0.20 -0.0007 -0.0049 0.0036 0.0400 0.0261
[17,] -0.15 -0.0040 -0.0105 0.0025 0.0225 0.0147
Rho at which ACME = 0: -0.2
R^2_M*R^2_Y* at which ACME = 0: 0.04
R^2_M~R^2_Y~ at which ACME = 0: 0.0261
\end{verbatim}
The table above presents the estimated mediation effect along with its
confidence interval for each value of $\rho$. The reader can verify
that when $\rho$ is equal to zero, the reported mediation effect
matches the estimate produced by the \texttt{mediate} function. For
other values of $\rho$, the mediation effect is calculated under
different levels of unobserved confounding.
The information from the sensitivity analysis can also be summarized
graphically using the \texttt{plot} function. First, passing the
\texttt{medsens} object to \texttt{plot} and specifying the
\texttt{sens.par} option to \texttt{"rho"}, i.e.,
\begin{verbatim}
R> plot(sens.cont, sens.par = "rho")
\end{verbatim}
produces the left hand side of Figure~\ref{jobplot1}. In the plot, the dashed
horizontal line represents the estimated mediation effect under the sequential
ignorability assumption, and the solid line represents the mediation
effect under various values of $\rho$. The grey region represents the
95\% confidence bands.
Similarly, we can also plot the sensitivity analysis in terms of the
coefficients of determination as discussed above. Here we specify \texttt{sens.par}
option to \texttt{"R2"}. We also need to specify two additional pieces of information.
First, \texttt{r.type} option tells the plot function whether to plot
${R^\ast}^2_M{R^\ast}^2_Y$ or $\widetilde{R}^2_M\widetilde{R}^2_Y$. To plot the
former \texttt{r.type} is set to \texttt{"residual"} and to plot the latter
\texttt{r.type} is set
to \texttt{"total"}. Finally, the \texttt{sign.prod} option specifies the sign of the
product of the coefficients of the unobserved confounder in the mediator and outcome models.
This product indicates whether the unobserved confounder affects both mediator and outcome
variables in the same direction (\texttt{"positive"}) or different directions
(\texttt{"negative"}), thereby
reflecting the analyst's expectation about the nature of confounding.
For example, the following command produces the plot representing the sensitivity
of estimates with respect to the proportion of the original variances explained by
the unobserved confounder when the confounder is hypothesized to affect the mediator
and outcome variables in opposite directions.
\begin{verbatim}
R> plot(sens.cont, sens.par = "R2", r.type = "total", sign.prod = "negative")
\end{verbatim}
The resulting plot is shown on the right hand side of Figure~\ref{jobplot1}. Each contour
line represents the mediation effect for the corresponding values of
$\widetilde{R}^2_M$ and $\widetilde{R}^2_Y$. For example, the 0 contour line
corresponds to values of the product $\widetilde{R}^2_M\widetilde{R}^2_Y$ such
that the ACME is 0. As reported in the table, even a small proportion of
original variance unexplained by the confounder, $.02\%$, produces mediation
effects of 0. Accordingly, the right hand side of Figure~\ref{jobplot1} shows
how increases in $\widetilde{R}^2_M\widetilde{R}^2_Y$ (moving from the lower
left to upper right) produce \emph{positive} mediation effects.
For both types of sensitivity plots, the user can specify additional options
available in the plot function such as alternative
title (\texttt{main}) and axis labels (\texttt{xlab}, \texttt{ylab}) or
manipulate common graphical options (e.g., \texttt{xlim}).
\begin{figure}[t]
\spacingset{1}
\vspace{-.5in}
\begin{center}
\includegraphics[scale=.8]{Rho-R2-ContCont-JOBS.pdf}
\end{center}
\vspace{-.5in}
\caption{Sensitivity Analysis with Continuous Outcome and Mediator.
\label{jobplot1}}
\end{figure}
\paragraph{Binary Outcome.}
The \texttt{medsens} function also extends to analyses where the
mediator is binary and the outcome is continuous, as well as when the
mediator is continuous and the outcome is binary. If either variable
is binary, \texttt{medsens} takes an additional argument. For
example, recall the binary outcome model estimated earlier:
\begin{verbatim}
R> model.y <- glm(work1 ~ treat + job_seek + depress1 + econ_hard + sex + age
+ occp + marital + nonwhite + educ + income,
family = binomial(link = "probit"), data = jobs)
R> med.bout <- mediate(model.m, model.y, sims = 1000, treat = "treat",
mediator = "job_seek")
\end{verbatim}
The call to \texttt{medsens} works as before, with the output from the \texttt{mediate}
function passed through \texttt{medsens}.
\begin{verbatim}
R> sens.bout <- medsens(med.bout, rho.by = 0.05, sims = 1000)
\end{verbatim}
The \texttt{sims} option provides control over the number of Monte Carlo draws
used to compute confidence bands.
When either the mediator or outcome is binary, the exact values of
sensitivity parameters where the mediation effects are zero cannot be analytically
obtained as in the fully continuous case \citep[see][Section 4]{imai:keel:yama:10}.
Thus, this information is reported based on the signs of the estimated mediation
effects under various values of $\rho$ and corresponding coefficients of determination.
The usage of the \texttt{summary} function, however, remains
identical to the fully continuous case
in that the output table contains the estimated mediation
effects and the corresponding values of $\rho$ for which the
confidence region contains zero.
As in the case with continuous mediator and outcome variables, we can plot the
results of the sensitivity analysis. The following code produces Figure~\ref{SensBinOut}.
\begin{verbatim}
R> plot(sens.bout, sens.par = "rho")
R> plot(sens.bout, sens.par = "R2", r.type = "total", sign.prod = "positive")
\end{verbatim}
On the left hand side we plot the ACME in terms of $\rho$, while we plot the ACME
in terms of $\widetilde{R}^2_M$ and $\widetilde{R}^2_Y$ on the right hand side.
In the $\rho$ plot, the dashed line represents the estimated mediation effect
under sequential
ignorability, and the solid line represents the mediation effect under
various values of $\rho$. The grey region represents the 95\%
confidence bands. In the $\widetilde{R}^2$ plot the ACME is plotted against various values of $\widetilde{R}^2_M$ and $\widetilde{R}^2_Y$ and is interpreted in the same way as above.
\begin{figure}[t]
\spacingset{1}
\vspace{-.5in}
\begin{center}
\includegraphics[scale=.8]{Rho-R2-JobsBinOutcome.pdf}
\end{center}
\vspace{-.5in}
\caption{Sensitivity Analysis with Continuous Outcome and Binary Mediator.
\label{SensBinOut}}
\end{figure}
When the outcome is binary, the proportion of the the total effect due to
mediation can also be calculated as function of the sensitivity parameter $\rho$.
The \texttt{pr.plot} option in the plot command (in conjunction with the
\texttt{sens.par = "rho"} option) allows users to plot a summary
of the sensitivity analysis for the proportion mediated. For example, the
following call would provide a plot of this quantity:
\begin{verbatim}
R> plot(sens.bout, sens.par = "rho", pr.plot = TRUE)
\end{verbatim}
\paragraph{Binary Mediator.}
The final form of sensitivity analysis deals with the case where the
outcome variable is continuous but the mediator is binary. For the
purpose of illustration, we simply dichotomize the \texttt{job\_seek}
variable to produce a binary measure \texttt{job\_dich}. We fit a
probit model for the mediator and linear regression for the outcome
variable.
\begin{verbatim}
R> model.m <- glm(job_dich ~ treat + depress1 + econ_hard + sex + age
+ occp + marital + nonwhite + educ + income,
data = jobs, family = binomial(link = "probit"))
R> model.y <- lm(depress2 ~ treat + job_dich+ depress1 + econ_hard + sex + age
+ occp + marital + nonwhite + educ + income, data = jobs)
R> med.bmed <- mediate(model.m, model.y, sims = 1000, treat = "treat",
mediator = "job_dich")
R> sens.bmed <- medsens(med.bmed, rho.by = 0.05, sims = 1000)
\end{verbatim}
Again we can pass the output of the \texttt{medsens} function through the \texttt{plot} function,
\begin{verbatim}
R> plot(sens.bmed, sens.par = "rho")
\end{verbatim}
producing Figure \ref{SensBinMed}. The plot is interpreted in the same
way as the above cases. The user also has the option to plot sensitivity results
in terms of the coefficients of determination just as in the case with
continuous outcome and mediator variables.
\begin{figure}[t]
\spacingset{1}
\vspace{-.5in}
\begin{center}
\includegraphics[scale=.8]{BinaryMediator-Sensitivity-Rpaper.pdf}
\end{center}
\vspace{-.5in}
\caption{Sensitivity Analysis with Continuous Outcome and Binary Mediator.
\label{SensBinMed}}
\end{figure}
When the mediator variable is binary, the plotted values of the mediation effect
and their confidence bands may not be perfectly smooth curves due to simulation errors.
This is especially likely when the number of simulations (\texttt{sims}) is set
to a small value. In such situations, the user can choose to set the
\texttt{smooth.effect} and \texttt{smooth.ci} options to \texttt{TRUE} in the
\texttt{plot} function so that the corresponding values become smoothed
out via a lowess smoother before being plotted. Although this option
often makes the produced graph look nicer, the user should be cautious as
the adjustment could affect one's substantive conclusions in a
significant way. A recommended alternative is to increase the number of
simulations.
\section{Multiple Causal Mediation Analyses via {\tt mediations}}
\label{sec:mediations}
Researchers may face a situation where they want to conduct multiple causal
mediation analyses. For example, the researcher may have run several different
experiments that had different treatment conditions but had the same mediator
and outcome variables. Or, a researcher may measure several different mediating
variables in their experiment and wish to estimate causal mediation effects for
each of these cases. In these situations the researcher would need to write
separate lines of code for the different models. A more efficient way is to
loop over the different combinations of interest. The {\tt mediations} function
facilitates this process.
The {\tt mediations} function can be used to conduct causal mediation analyses
over the range of treatment, mediator, and outcome combinations that could be
present within a study. For example, if there were two treatment conditions, two
measured mediators, and two outcomes, mediations would run eight causal
mediation analyses. The results can then be analyzed with the {\tt
summary} and {\tt plot} functions. If there is only one
treatment variable, the user needs only a single data set that contains
all mediators and outcomes of interest. If there are multiple treatment
variables, then the researcher must first create separate data frames, each
with all mediator/outcome variables to be considered, and named using a
convention where the name of the data frame object begins with the name of the
treatment variable contained within it.
Use of the {\tt mediations} function is straightforward. First, the user creates
vectors of the treatment, mediator, outcome, and covariates names that are to be used.
The treatment, mediator, and outcome variables are looped over to produce
all possible combinations of mediator and outcome models. The set of covariates
is held constant across the models. The different model combinations are then
entered into the {\tt mediate} function.
The user can specify several different
options for {\tt mediations}. First, the type of model to be used for the
mediator and outcome equations must be specified. Currently, {\tt mediations} can
handle a variety of different model combinations. Either the mediator or outcome
equation can be modeled with linear regression, probit, ordered probit, and
quantile regression. The {\tt mediations} function also handles the case where tobit is used for
the outcome. Next, the user must also specify the standard features that are necessary for the
{\tt mediate} function, such as the number of simulations ({\tt sims}). In addition, if
quantile regression or tobit regression is used, the user needs to specify
the appropriate parameters, such as the quantiles to be estimated ({\tt tau.m},
{\tt tau.y}) and cutpoints ({\tt LowerY}, {\tt UpperY}).
If the user wishes for estimates to include weights
(e.g., survey weights) then the name of the weight variable needs to be passed
via the {\tt weights} argument.
The output of {\tt mediations} contains a collection of {\tt mediate} objects. The naming
system uses the name of the particular outcome variable, mediator, and
treatment, each separated by a dot. For example, if the output of {\tt mediations}
is {\tt output} and the outcome, mediator and treatment variable names
{\tt O}, {\tt M1} or {\tt M2}, and {\tt T}, respectively,
then {\tt output\$O.M2.T} would be the {\tt mediate} object from using mediator
{\tt M2}.
There are several limitations to the code. First, it works only with a subset of
the model types that will be accommodated if {\tt mediate} is used individually.
Second, one cannot specify separate sets of covariates for different
treatment/mediator/outcome combinations. Users should use {\tt mediate} separately
for individual models if more flexibility is required in their specific
applications.
While {\tt mediations} is a helpful function, it is also easy to abuse.
First, researchers must be aware of the assumptions being made at every step of
causal mediation analysis. For example, if there are multiple mediators, a
strong assumption must be made for the estimates to be interpreted as causal
quantities. One such assumption is that the mediators are not causally related.
Second, care must be taken against the risks of multiple testing when a
large number of treatment/mediator/outcome combinations are involved. The
{\tt mediations} function simply runs independent {\tt mediate} processes and
does not use any correction procedures for multiple testing, so the confidence
intervals may be too narrow. We recommend that the use of this function be limited
for exploratory analyses.
An illustration of the {\tt mediations} function can be found in the {\tt Examples}
section of its help document; see {\tt ?mediations} or {\tt example(mediations)}.
\section{Concluding Remarks}
Causal mediation analysis is a key tool for social scientific
research. In this paper, we describe our easy-to-use software for
causal mediation analysis, \bmediation, that implements the new
methods and algorithms introduced by \citet*{imai:keel:yama:10} and
\citet*{imai:keel:ting:10}. The software provides a flexible, unified
approach to causal mediation analysis in various situations
encountered by applied researchers. The object-oriented nature of the
\bR\ programming made it possible for us to implement these algorithms
in a fairly general way. In addition to the estimation of causal
mediation effects, \bmediation\ implements formal sensitivity analyses
so that researchers can assess the robustness of their findings to the
potential violations of the key identifying assumption. This is an
important contribution for at least two reasons. First, even in
randomized experiments, causal mediation analysis requires an
additional assumption that is not directly testable from the observed
data. Thus, researchers must evaluate the consequences of potential
violations of the assumption via sensitivity analysis. Second, the
accumulation of such sensitivity analyses is essential for
interpreting the relative degree of robustness across different
studies. Thus, the development of easy-to-use software, such as
\bmediation, facilitates causal mediation analysis in applied social
science research in several critical directions.
\clearpage
\pdfbookmark[1]{References}{References}
\bibliography{mediation-old}
\end{document}