Causal Mediation Analysis with Multiple Treatments and Latent Confounders

Let $\mathbf{Z}_{i}$ denote the observed treatments assigned to individual $i$ which is a vector with $J\geq 2$ treatment variables, i.e., $\mathbf{Z}_{i}=(Z_{i1},Z_{i2},\ldots,Z_{iJ})^{\top}$. Let $Y_{i}$ denote the observed outcome for individual $i$ and $M_{i}$ denote some observed intermediate variable on the causal path from the treatments to the outcome. To streamline notation of the random variables, we suppress subscript $i$ for individual below.

The components of $\mathbf{Z}$ can be correlated with each other through an unobserved common cause of them. For each $j=1,\ldots,J$, assume that $\Pr(Z_{j}=z_{j})>0$ for any $z_{j}$th treatment level with $z_{j}\in\{1,\ldots,K_{j}\}$. Both $Y$ and $M$ are assumed to be continuous and there may exist a latent confounder $U$ confounding the relationship between these two variables. In general, we may also consider $M$ as a vector of mediator variables and our results in this article can be straightforwardly generalized from the setting of a single mediator to the setting of multiple mediators.

To formally define causal effects in causal mediation analysis, we make use of the concept of potential outcomes. Potential outcomes present the values of a outcome variable for each individual under varying levels of a treatment variable. We can observe only one of these potential outcomes but can never observe all of them because it is impossible for us to unwind time and go back and manipulate the individual to other treatment levels.

We first make the stable unit treatment value assumption (SUTVA). This assumption requires that the value of the outcome should not be affected by the manner of manipulations providing the same value for the treatment variable, that is, there is only one version of the potential outcomes and there is no interference between individuals Rubin (1980). The SUTVA allows us to uniquely define the potential values for the mediator $M(\mathbf{z})$ and the potential outcome $Y(\mathbf{z})$ if an individual were to receive treatment $\mathbf{Z}=\mathbf{z}$. Let $Y(\mathbf{z},m)$ denote the potential outcome for an individual that would occur if the treatment $\mathbf{Z}$ were set to level $\mathbf{z}$, and if the mediator $M$ were manipulated to level $m$. In contrast, let $Y(\mathbf{z},M(\mathbf{z}^{*}))$ denote the potential outcome for an individual, where we do not specify the actual level of $M$, but set it to what it would have been if treatment had been $\mathbf{Z}=\mathbf{z}^{*}$. To connect the observed random variables with corresponding potential outcomes, we also make the consistency assumption, namely that $M(\mathbf{z})=M$, $Y(\mathbf{z})=Y$ if $\mathbf{Z}=\mathbf{z}$, and $Y(\mathbf{z},m)=Y$ if $\mathbf{Z}=\mathbf{z}$, $M=m$. According to this assumption, we note that the observed outcome $Y$ is just one realization of the potential outcome $Y(\mathbf{z},m)$ with observed treatment level $\mathbf{Z}=\mathbf{z}$ and mediator level $M=m$.

We also assume that the underlying causal model corresponds to a directed acyclic graph (DAG, Pearl (2000)). The DAG is a useful tool for visual representations of qualitative causal relationships between the variables of interest. We show the corresponding DAG of this context in Figure 1. Note that there are no causal relationships between the treatment variables in $\mathbf{Z}$, but they may be associated with each other through some unobserved variable. We use $C$ to denote this unobserved common cause of the treatment variables in $\mathbf{Z}$, and it is assumed to be independent of the other variables conditional on $\mathbf{Z}$. When the treatment variables in $\mathbf{Z}$ are randomly assigned, $C$ is an empty set, and this assumption automatically holds. The symbol ‘ $\vdots$ ’ in the DAG represents other undisplayed treatment variables $Z_{j}$’s, each of which has the directed edges $C\rightarrow Z_{j}$, $Z_{j}\rightarrow M$, and $Z_{j}\rightarrow Y$.

Using the nested potential outcomes notation, we can define the causal parameters of interest in this multiple-treatment model. We first describe the definition of the average causal effect in a single-treatment setting. We use the difference $\textnormal{E}\{Y(z)-Y(z^{*})\}$ to represent the average causal effect of treatment level $z$ versus treatment level $z^{*}$. We now extend this definition to the setting with multiple treatment variables. For notational simplicity, we let $\mathbf{Z}_{-j}=(Z_{1},\ldots,Z_{j-1},Z_{j+1},\ldots,Z_{J})$ and $\mathbf{z}_{-j}=(z_{1},\ldots,z_{j-1},z_{j+1},\ldots,z_{J})$. Then we let $CTE(z_{j},z_{j}^{*}\mid\mathbf{z}_{-j})$ denote the conditional total causal effect of $Z_{j}$ on $Y$ under two differing levels $z_{j}$ and $z_{j}^{*}$ of $Z_{j}$ while conditioning on the interventions $\mathbf{z}_{-j}$ for the other treatments. The formulation is given as follows:

$\displaystyle CTE(z_{j},z_{j}^{*}\mid\mathbf{z}_{-j})=\textnormal{E}\{Y(z_{j},\mathbf{z}_{-j})\}-\textnormal{E}\{Y(z_{j}^{*},\mathbf{z}_{-j})\}.\vspace{-1em}$

In addition, we can also define the average total causal effect, $TE(z_{j},z_{j}^{*})$, which is free of other treatment variables, by taking expectation of $CTE(z_{j},z_{j}^{*}\mid\mathbf{Z}_{-j})$ with respect to $\mathbf{Z}_{-j}$, i.e.,

$\displaystyle TE(z_{j},z_{j}^{*})=\textnormal{E}\{CTE(z_{j},z_{j}^{*}\mid\mathbf{Z}_{-j})\}.$

We next define the conditional natural direct and indirect effects of $Z_{j}$ on $Y$. Let $CNDE(z_{j},z_{j}^{*}\mid\mathbf{z}_{-j})$ denote the conditional natural direct effect of $Z_{j}$ under two different levels $z_{j}$ and $z_{j}^{*}$ while setting the other treatments to $\mathbf{z}_{-j}$ and setting $M$ to the values attained under fixed treatment levels $\mathbf{z}$. We use $CNIE(z_{j},z_{j}^{*}\mid\mathbf{z}_{-j})$ to denote the conditional natural indirect effect which is defined as the difference between two averaged potential outcomes with treatments set to $z_{j}^{*},\mathbf{z}_{-j}$ and the mediator $M$ set to the values attained under differing treatment levels $z_{j}$ and $z_{j}^{*}$ for $Z_{j}$ and $\mathbf{z}_{-j}$ for others. Below we formally give their definitions:

	$\displaystyle CNDE(z_{j},z_{j}^{*}\mid\mathbf{z}_{-j})=$	$\displaystyle\textnormal{E}\{Y(z_{j},\mathbf{z}_{-j},M(\mathbf{z}))\}$
		$\displaystyle~{}-\textnormal{E}\{Y(z_{j}^{*},\mathbf{z}_{-j},M(\mathbf{z}))\},$
	$\displaystyle CNIE(z_{j},z_{j}^{*}\mid\mathbf{z}_{-j})=$	$\displaystyle\textnormal{E}\{Y(z_{j}^{*},\mathbf{z}_{-j},M(z_{j},\mathbf{z}_{-j}))\}$
		$\displaystyle~{}-\textnormal{E}\{Y(z_{j}^{*},\mathbf{z}_{-j},M(z_{j}^{*},\mathbf{z}_{-j}))\}.$

Similarly, we also give the the definitions of the average natural direct and indirect effects of $Z_{j}$ on $Y$ which do not depend on other treatment variables:

	$\displaystyle NDE(z_{j},z_{j}^{*})={}$	$\displaystyle\textnormal{E}\{CNDE(z_{j},z_{j}^{*}\mid\mathbf{Z}_{-j})\},$
	$\displaystyle NIE(z_{j},z_{j}^{*})={}$	$\displaystyle\textnormal{E}\{CNIE(z_{j},z_{j}^{*}\mid\mathbf{Z}_{-j})\}.$

Under the composition assumption (Pearl, 2000) that $Y(\mathbf{z})=Y(\mathbf{z},M(\mathbf{z}))$, we immediately have the following decompositions:

	$\displaystyle CTE(z_{j},z_{j}^{*}\mid\mathbf{z}_{-j})=$	$\displaystyle CNDE(z_{j},z_{j}^{*}\mid\mathbf{z}_{-j})$
		$\displaystyle~{}+CNIE(z_{j},z_{j}^{*}\mid\mathbf{z}_{-j}),$
	$\displaystyle TE(z_{j},z_{j}^{*})=$	$\displaystyle NDE(z_{j},z_{j}^{*})+NIE(z_{j},z_{j}^{*}).$

We use potential outcomes notation to construct causal models that allow us to directly specify the causal effects of interest as functions of parameters in the models. We consider the following causal models including a latent confounder $U$ for potential outcomes, $M(\mathbf{z})$ and $Y(\mathbf{z},m)$, respectively:

$\displaystyle\begin{aligned} M(\mathbf{z})&=g_{M}^{c}(\mathbf{z})+U+\epsilon(\mathbf{z}),\\ Y(\mathbf{z},m)&=g_{Y}^{c}(\mathbf{z})+\beta^{c}m+h(U)+\eta(\mathbf{z},m),\end{aligned}$

(1)

where $\textnormal{E}(U)=\textnormal{E}\{h(U)\}=0$, $g_{M}^{c}(\cdot)$, $g_{Y}^{c}(\cdot)$ and $h(\cdot)$ are unknown functions. Note that the causal models in (1) allow the individual natural direct, indirect and total causal effects vary individual by individual. Similar models were also proposed in Lindquist (2012). Here we use the superscript ‘$c$’ to indicate causal parameters in models (1) for potential outcomes to distinguish the parameters in SEMs for observed variables.

We assume that $\textnormal{E}\{\epsilon(\mathbf{z})\}=0$ and $\textnormal{E}\{\eta(\mathbf{z},m)\}=0$ and that $\epsilon(\mathbf{z})$ and $\eta(\mathbf{z},m)$ are mutually independent and are also independent of $(\mathbf{Z},M,U)$ for all values of $\mathbf{z}$ and the pair $(\mathbf{z},m)$. Then we can obtain $\textnormal{E}\{\eta(\mathbf{z},M(\mathbf{z}^{*}))\}=0$ for all values of $\mathbf{z}$ and $\mathbf{z}^{*}$. With these conditions, we can directly write the average natural direct, indirect and total causal effects of $Z_{j}$ on $Y$ respectively as

	$\displaystyle NDE(z_{j},z_{j}^{*})$	$\displaystyle=\textnormal{E}\{g_{Y}^{c}(z_{j},\mathbf{Z}_{-j})\}-\textnormal{E}\{g_{Y}^{c}(z_{j}^{*},\mathbf{Z}_{-j})\},$		(2)
	$\displaystyle NIE(z_{j},z_{j}^{*})$	$\displaystyle=\beta^{c}\big{[}\textnormal{E}\{g_{M}^{c}(z_{j},\mathbf{Z}_{-j})\}-\textnormal{E}\{g_{M}^{c}(z_{j}^{*},\mathbf{Z}_{-j})\}\big{]},$
	$\displaystyle TE(z_{j},z_{j}^{*})$	$\displaystyle=NDE(z_{j},z_{j}^{*})+NIE(z_{j},z_{j}^{*}).$

The identifiability of these causal effects relies on the identifiability of $g_{M}^{c}(\cdot)$, $g_{Y}^{c}(\cdot)$ and $\beta^{c}$. These parameters encode the causal relationships in models (1) for the potential outcomes. Since the variables in the models are not observed, these parameters cannot be estimated by the ordinary approaches for the models with observed variables. Hence, additional assumptions or (and) more feasible models are required to identify and estimate these parameters.

According to the path diagram in Figure 1, we consider the SEMs for observed variables and a latent confounder $U$ as follows:

	$\displaystyle M(\mathbf{Z})$	$\displaystyle=g_{M}^{s}(\mathbf{Z})+U+\epsilon,$		(3)
	$\displaystyle Y(\mathbf{Z},M)$	$\displaystyle=g_{Y}^{s}(\mathbf{Z})+\beta^{s}M+h(U)+\eta,$

where $\textnormal{E}(\epsilon\mid\mathbf{Z},U)=0$, $\textnormal{E}(\eta\mid\mathbf{Z},M,U)=0$, $g_{M}^{s}(\cdot)$ and $g_{Y}^{s}(\cdot)$ are unknown functions, and the variables $M(\mathbf{Z})$ and $Y(\mathbf{Z},M)$ are observed because of the consistency assumption: $M(\mathbf{Z})=M$ and $Y(\mathbf{Z},M)=Y$. The SEMs are built upon observed variables, which are different from those for the potential outcomes defined in the previous subsection. The superscript ‘$s$’ denotes the parameters occurring in the SEMs, and these parameters are in principle estimable from observed data.

The SEMs in (3) assumes constant individual natural direct, indirect and total causal effects, which is more restrictive than that imposed in the causal models (1). However, we are interested in the average versions of these causal effects, and the average causal effects can be expressed as causal parameters $\{g_{M}^{c}(\cdot)$, $g_{Y}^{c}(\cdot)$, $\beta^{c}\}$ in models (1) according to (2). In general, the parameters $\{g_{M}^{s}(\cdot),g_{Y}^{s}(\cdot),\beta^{s}\}$ in the SEMs (3) are not equal to their counterparts $\{g_{M}^{c}(\cdot)$, $g_{Y}^{c}(\cdot)$, $\beta^{c}\}$ in the causal models. However, under the assumptions encoded by the DAG in Figure 1, we can establish the equality between these two sets of parameters.

Under the models in (3), we can write the parameters $\{g_{M}^{s}(\cdot),g_{Y}^{s}(\cdot),\beta^{s}\}$ as

	$\displaystyle g_{M}^{s}(\mathbf{z})=\textnormal{E}\{M(\mathbf{Z})\mid\mathbf{Z}=\mathbf{z}\},{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}{}$
	$\displaystyle\beta^{s}=\frac{1}{m-m^{*}}\big{[}\textnormal{E}\{Y(\mathbf{Z},M(\mathbf{Z}))\mid M(\mathbf{Z})=m,\mathbf{Z}=\mathbf{z},U\}$
	$\displaystyle-\textnormal{E}\{Y(\mathbf{Z},M(\mathbf{Z}))\mid M(\mathbf{Z})=m^{*},\mathbf{Z}=\mathbf{z},U\}\big{]},$
	$\displaystyle g_{Y}^{s}(\mathbf{z})=\textnormal{E}\{Y(\mathbf{Z},M(\mathbf{Z}))\mid M(\mathbf{Z})=m,\mathbf{Z}=\mathbf{z},U\}$
	$\displaystyle-\beta^{s}m-h(U).$

Noting first from the DAG that $M(\mathbf{z})\mathbin{\mathchoice{\hbox to0.0pt{\hbox{\set@color$\displaystyle\perp$}\hss}\kern 3.46875pt{}\kern 3.46875pt\hbox{\set@color$\displaystyle\perp$}}{\hbox to0.0pt{\hbox{\set@color$\textstyle\perp$}\hss}\kern 3.46875pt{}\kern 3.46875pt\hbox{\set@color$\textstyle\perp$}}{\hbox to0.0pt{\hbox{\set@color$\scriptstyle\perp$}\hss}\kern 2.36812pt{}\kern 2.36812pt\hbox{\set@color$\scriptstyle\perp$}}{\hbox to0.0pt{\hbox{\set@color$\scriptscriptstyle\perp$}\hss}\kern 1.63437pt{}\kern 1.63437pt\hbox{\set@color$\scriptscriptstyle\perp$}}}\mathbf{Z}$ and under the assumption that $\textnormal{E}(U+\epsilon\mid\mathbf{Z}=\mathbf{z})=0$ for any value of $\mathbf{z}$, we immediately have the following result:

$\displaystyle g_{M}^{s}(\mathbf{z})=\textnormal{E}\{M(\mathbf{Z})\mid\mathbf{Z}=\mathbf{z}\}=\textnormal{E}\{M(\mathbf{z})\}=g_{M}^{c}(\mathbf{z}).$

We also note that the potential outcomes for the mediator and outcome are conditionally independent given the latent confounder $U$, i.e., $Y(\mathbf{z},m)\mathbin{\mathchoice{\hbox to0.0pt{\hbox{\set@color$\displaystyle\perp$}\hss}\kern 3.46875pt{}\kern 3.46875pt\hbox{\set@color$\displaystyle\perp$}}{\hbox to0.0pt{\hbox{\set@color$\textstyle\perp$}\hss}\kern 3.46875pt{}\kern 3.46875pt\hbox{\set@color$\textstyle\perp$}}{\hbox to0.0pt{\hbox{\set@color$\scriptstyle\perp$}\hss}\kern 2.36812pt{}\kern 2.36812pt\hbox{\set@color$\scriptstyle\perp$}}{\hbox to0.0pt{\hbox{\set@color$\scriptscriptstyle\perp$}\hss}\kern 1.63437pt{}\kern 1.63437pt\hbox{\set@color$\scriptscriptstyle\perp$}}}M(\mathbf{z})\mid U$ for any value of the pair $(\mathbf{z},m)$. Combining this with the ignorable treatment assignment assumption, we have

		$\displaystyle\textnormal{E}\{Y(\mathbf{Z},M(\mathbf{Z}))\mid M(\mathbf{Z})=m,\mathbf{Z}=\mathbf{z},U\}$
	$\displaystyle={}$	$\displaystyle\textnormal{E}\{Y(\mathbf{z},m)\mid M(\mathbf{z})=m,U\}=\textnormal{E}\{Y(\mathbf{z},m)\mid U\}.$

Consequently,

	$\displaystyle\beta^{s}$	$\displaystyle=\frac{1}{m-m^{*}}\big{[}\textnormal{E}\{Y(\mathbf{z},m)\mid M(\mathbf{z})=m,U\}$
		$\displaystyle~{}~{}~{}~{}~{}~{}~{}-\textnormal{E}\{Y(\mathbf{z},m^{*})\mid M(\mathbf{z})=m^{*},U\}\big{]}$
		$\displaystyle=\frac{1}{m-m^{*}}\big{[}\textnormal{E}\{Y(\mathbf{z},m)\mid U\}-\textnormal{E}\{Y(\mathbf{z},m^{*})\mid U\}\big{]}$
		$\displaystyle=\beta^{c}.$

In addition, it should also be noted that based on the previous results, we can also show the equality of parameter $g_{Y}^{s}(\cdot)$ with parameter $g_{Y}^{c}(\cdot)$ as follows:

$\displaystyle g_{Y}^{s}(\mathbf{z})$

$\displaystyle=\textnormal{E}\{Y(\mathbf{z},m)\mid U\}-\beta^{s}m-h(U)=g_{Y}^{c}(\mathbf{z}).$

Until now, we have shown the equivalence between parameters of the SEMs in (3) with the corresponding parameters of the causal models in (1). Since the SEMs include the unobserved variable $U$, the parameters are still unidentifiable without additional conditions.

In this subsection, we give conditions under which the parameters of the SEMs are identifiable from observed data. Apparently, $g_{M}^{s}(\cdot)$ is identifiable due to the condition that $\textnormal{E}(U+\epsilon\mid\mathbf{Z})=0$ and can be written as follows:

$$g_{M}^{s}(\mathbf{z})=\textnormal{E}(M\mid\mathbf{Z}=\mathbf{z}).$$

In order to guarantee the identifiability of $\beta^{s}$ and $g_{Y}^{s}(\cdot)$, we impose the following condition.

The condition 1 implies that the number of the basis functions of $\mathscr{G}_{M}^{s}$ is finite and a proper subset of these basis functions generates the subspace $\mathscr{G}_{Y}^{s}$. This may be not a stringent condition and can be satisfied in a variety of cases. For example, suppose that $g_{M}^{s}(\mathbf{z})$ is a polynomial function of degree 2 in each component of $\mathbf{z}$. Then, a linear function $g_{Y}^{s}(\mathbf{z})$ of the components satisfies Condition 1.

Based on the discussions about the equality between the parameters of the SEMs and the analogous parameters of the causal models in the previous subsection, we conclude from Theorem 1 that the parameters $g_{M}^{c}(\cdot)$, $g_{Y}^{c}(\cdot)$ and $\beta^{c}$ are also identifiable, which in turn implies the identifiability of the average natural direct, indirect and total causal effects of $Z_{j}$ on $Y$ according to (2).

In this subsection, we provide an approach for estimating the parameters of the SEMs in (3). We also discuss the resampling-based procedures for inference.

Given the known basis functions $\phi_{1}(\mathbf{z}),\ldots,\phi_{L}(\mathbf{z})$ of $\mathscr{G}_{M}^{s}$, we can express $g_{M}^{s}(\mathbf{z})$ and $g_{Y}^{s}(\mathbf{z})$ as linear combinations of them, which have been shown in (4). Denote $\bm{\Phi}(\mathbf{z})=(\phi_{1}(\mathbf{z}),\ldots,\phi_{L}(\mathbf{z}))$, $\bm{\Phi}_{1}(\mathbf{z})=(\phi_{1}(\mathbf{z}),\ldots,\phi_{L_{1}}(\mathbf{z}))$, $\bm{\alpha}=(\alpha_{1},\ldots,\alpha_{L})^{\top}$, $\bm{\gamma}=(\gamma_{1},\ldots,\gamma_{L_{1}})^{\top}$. We then rewrite $g_{M}^{s}(\mathbf{z})$ and $g_{Y}^{s}(\mathbf{z})$ as

$$g_{M}^{s}(\mathbf{z})=\bm{\Phi}(\mathbf{z})\bm{\alpha},~{}~{}\text{and}~{}~{}g_{Y}^{s}(\mathbf{z})=\bm{\Phi}_{1}(\mathbf{z})\bm{\gamma}.$$

By the least-square criterion, the coefficient $\bm{\alpha}$ of the expansion of $g_{M}^{s}(\mathbf{z})$ can be determined by minimizing

$$\hat{\bm{\alpha}}=\arg\!\min_{\bm{\alpha}}\sum_{i=1}^{n}\big{\{}M_{i}-\bm{\Phi}(\mathbf{Z}_{i})\bm{\alpha}\big{\}}^{2}.$$

Then the solution is given by

$$\hat{\bm{\alpha}}=\bigg{\{}\sum_{i=1}^{n}\bm{\Phi}(\mathbf{Z}_{i})^{\top}\bm{\Phi}(\mathbf{Z}_{i})\bigg{\}}^{-1}\bigg{\{}\sum_{i=1}^{n}\bm{\Phi}(\mathbf{Z}_{i})^{\top}M_{i}\bigg{\}},$$

which provides us an estimate of the coefficients of $g_{M}^{s}(\mathbf{z})$.

To derive estimates of $\beta^{s}$ and $\bm{\gamma}$, we utilize the idea of the generalized method of moments (GMM) Hall (2005) and introduce more notations. Let $\bm{\alpha}_{1}=(\alpha_{1},\ldots,\alpha_{L_{1}})^{\top}$, and the subvector being composed of the remaining components in $\bm{\alpha}$ is denoted by $\bm{\alpha}_{2}$. Let $\bm{\Phi}_{2}(\mathbf{z})$ represent the subvector of $\bm{\Phi}(\mathbf{z})$ with components not contained in $\bm{\Phi}_{1}(\mathbf{z})$, i.e., $\bm{\Phi}(\mathbf{z})=\{\bm{\Phi}_{1}(\mathbf{z}),\bm{\Phi}_{2}(\mathbf{z})\}$. Then, we denote $\bm{Y}=(Y_{1},\ldots,Y_{n})^{\top}$, $\bm{\Phi}=\{\bm{\Phi}(\mathbf{Z}_{1})^{\top},\ldots,\bm{\Phi}(\mathbf{Z}_{n})^{\top}\}^{\top}$, $\bm{\Phi}_{1}=\{\bm{\Phi}_{1}(\mathbf{Z}_{1})^{\top},\ldots,\bm{\Phi}_{1}(\mathbf{Z}_{n})^{\top}\}^{\top}$, and $\bm{\Phi}_{2}=\{\bm{\Phi}_{2}(\mathbf{Z}_{1})^{\top},\ldots,\bm{\Phi}_{2}(\mathbf{Z}_{n})^{\top}\}^{\top}$. For notational simplicity, we also let $\bm{\delta}=\bm{\gamma}+\beta^{s}\bm{\alpha_{1}}$. Using these notations, we now rewrite (5) as

$$Y=\bm{\Phi}_{1}(\mathbf{Z})\bm{\delta}+\beta^{s}\bm{\Phi}_{2}(\mathbf{Z})\bm{\alpha}_{2}+\beta^{s}U+h(U)+\beta^{s}\varepsilon+\eta.$$

In view of $\textnormal{E}\{\beta^{s}U+h(U)+\beta^{s}\varepsilon+\eta\mid\mathbf{Z}\}=0$, it follows immediately that

$$\textnormal{E}\big{[}\Phi(\mathbf{Z})^{\top}\{Y-\bm{\Phi}_{1}(\mathbf{Z})\bm{\delta}-\beta^{s}\bm{\Phi}_{2}(\mathbf{Z})\bm{\alpha}_{2}\}\big{]}=\bm{0}.$$

(6)

We then estimate $(\bm{\delta}^{\top},\beta^{s})$ by minimizing the sum of the squares of sample analogues of (6)

		$\displaystyle(\hat{\bm{\delta}}^{\top},\hat{\beta^{s}})=\arg\!\min_{(\bm{\delta}^{\top},\beta^{s})}(\bm{Y}-\bm{\Phi}_{1}\bm{\delta}-\beta^{s}\bm{\Phi}_{2}\bm{\alpha}_{2})^{\top}\bm{\Phi}\bm{\Phi}^{\top}$
		$\displaystyle~{}~{}~{}\qquad\qquad\qquad\qquad~{}~{}(\bm{Y}-\bm{\Phi}_{1}\bm{\delta}-\beta^{s}\bm{\Phi}_{2}\bm{\alpha}_{2}).$

After solving the above optimization problem and substituting the estimator $\hat{\bm{\alpha}}_{2}$ into the solutions, we can easily obtain the estimates of $(\bm{\delta}^{\top},\beta^{s})$ with explicit forms. However, due to their complex expressions, we omit displaying them here. Let $\hat{\bm{\gamma}}=\hat{\bm{\delta}}-\hat{\beta^{s}}\hat{\bm{\alpha}}_{1}$, which gives us an estimate of $\bm{\gamma}$. In addition, all these estimators $\hat{\bm{\alpha}}$, $\hat{\beta^{s}}$ and $\hat{\bm{\gamma}}$ are consistent. Under some regularity conditions, they are also asymptotically normal. Substituting these estimators into (2) and taking averages over samples, we obtain consistent estimators of the average causal effects of interest. Furthermore, by the Delta method, these consistent estimators are also asymptotically normal. Since analytical calculations of the asymptotic variances are difficult, we use a nonparametric bootstrap method to conduct inference.

We consider two different simulation studies where data are generated by using only one and multiple (more than one) available instrumental variables, respectively. Each of the simulations is repeated 1000 times under different sample sizes $n=$ 500, 1000, and 2000. We report the results of estimated causal effects of interest by averaging over the 1000 replications.

In both of the simulation studies, we consider three treatment variables: $Z_{1}$, $Z_{2}$, and $Z_{3}$. Each of them is uniformly generated from $\{1,2,3\}$ with equal probability. The latent confounder $U$ is generated from $N(0,1)$. The mediator $M$ is generated from the following equation:

	$\displaystyle M=$	$\displaystyle Z_{1}+Z_{2}+Z_{3}+Z_{1}*Z_{2}+Z_{1}*Z_{3}$
		$\displaystyle~{}~{}+Z_{2}*Z_{3}+U+\epsilon_{M},$

where $\epsilon_{M}\sim N(0,1)$. The only difference between the two simulation studies is the generation of the outcome $Y$. The details for the generation of $Y$ are described as follows:

In this section, we apply the proposed approach to a real customer loyalty data set from a telecom company in China. The customer loyalty analysis tries to discover key factors affecting loyalty and the affecting process, with which to choose proper actions to maintain customers. In this study, 18553 randomly chosen customers answer questionnaire via personal interview or an online survey. We drop the incompletely observed participants, leaving a total of 9833 participants. In the questionnaire, customers scored their satisfaction about 13 different specific factors, respectively, such as network quality, tariff plan, voice quality, service quality, and so on. The scores are integers ranging from 1 to 10, with 1 denoting “not satisfied at all” and 10 denoting “satisfied”. For the simplicity, we choose two important treatment factors, network quality and tariff plan, as an example to illustrate our approach, which are denoted by $Z_{1}$ and $Z_{2}$, respectively. For each customer, a score indicating the loyalty to the company has also been obtained, and we use the score as the outcome variable $Y$. Investigators also collected the information of each customer’s general satisfaction about the company, and it is used as the mediator variable $M$. We aim to assess whether the treatment factors have significant effects on improving customers’ loyalty. In particular, we wish to figure out whether the effects of the factors on loyalty are mediated by the general satisfaction about the company.

We use the proposed approach for evaluation, and also consider the traditional approach without accounting for possible latent confounders for comparison. For both approaches, we calculate estimators of the average natural direct effect $NDE_{j}(2,1)$, the average natural indirect effect $NIE_{j}(2,1)$, and the average total effect $TE_{j}(2,1)$, and also their 95% confidence intervals. Of note, the average total causal effect $TE_{j}(2,1)$ is estimated by the sum of the estimated average natural direct effect $NDE_{j}(2,1)$ and indirect effect $NIE_{j}(2,1)$. Subscript ‘$j$’ indicates the corresponding causal effects of network quality and tariff plan for $j=1,2$. We report these results in Table 3.

From Table 3, we can see that both the proposed approach and the traditional approach indicate positive and statistically significant causal effects of network quality and tariff plan on loyalty. Since there are no latent confounders between the treatments (network quality and tariff plan) and outcome variable (loyalty) by the randomization of this study, a simple linear regression of the outcome against the treatments can induce good estimates of the average total causal effects. These results are exactly the same as the corresponding ones obtained from the proposed approach, but different from those obtained from the traditional approach. From this point of view, the results for our proposed approach in Table 3 are more reliable, compared with those obtained from the traditional approach which does not consider possible latent confounders between the general satisfaction and loyalty. According to the results, both the estimated average natural direct effects of network quality and tariff plan on loyalty are positive but larger than the corresponding indirect effects through the mediator variable (the general satisfaction about the company). This means that the causal effects of satisfaction about network quality and tariff plan on loyalty are mostly operated in a direct way, but only a small fraction of them are intermediated by the general satisfaction about the company. Note also that the estimated total causal effect of network quality is a bit larger than that of tariff plan, which to some extent implies that the network quality may play a more important role in affecting loyalty than tariff plan.