Sensitivity Analysis of Causal Treatment Effect Estimation for Clustered Observational Data with Unmeasured Confounding

Consider a binary treatment variable $A$ with value 1 for the treated and $0$ for the control, and a continuous outcome variable $Y$. Denote X as a set of measured confounders, and U as a set of unmeasured confounders. Without loss of generality, we assume one measured confounder and one unmeasured confounder. Let $Y^{A=a}$ (or $Y^{a}$ ) be the counterfactual or potential outcome, which is the outcome that would have been observed if individual had been given treatment $a$. Denote $Y$ as the observed outcome that corresponds to the treatment value actually experienced.

Individual causal effects are defined as the difference between potential outcomes. Because only one of those outcomes is observed for each individual, individual effects cannot be identified. The average treatment effects (ATE) and conditional average treatment effects (CATE) were proposed and can be estimated from the observed data when additional assumptions hold. The causal ATE of $A$ on $Y$ is the difference between the expectation of the potential outcomes in the population while the conditional causal ATE is the ATE in any subpopulation.

	$$\displaystyle\text{ATE}=\mathbb{E}\left(Y^{A=1}\right)-\mathbb{E}\left(Y^{A=0}\right),$$
	$$\displaystyle\text{CATE}=\mathbb{E}\left(Y^{A=1}\mid X=x,U=u\right)-\mathbb{E}\left(Y^{A=0}\mid X=x,U=u\right).$$

To estimate the ATE or CATE, the following causal assumptions are necessary. The stable unit treatment value assumption (SUTVA) indicates that treatment assignment of one subject does not affect the outcome of another subject. The consistency assumption requires that the potential outcome under treatment $A=a$, $Y^{A=a}$, is equal to the observed outcome if treatment $A=a$ is actually received, i.e., $Y^{a}=Y\text{ if }A=a,\text{ for all }a.$ The positivity assumption states that treatment assignment is not deterministic for every set of values of $X$ and $U$, i.e., $0<P(A=a\mid X=x,U=u)<1\text{ for all }a\text{, }x\text{ and }u.$ The exchangeability assumption requires that given $X$ and $U$, treatment assignment is independent of the potential outcomes, i.e., $Y^{a}\perp A\mid X,U$ for all $a$.

Under these assumptions, ATE and CATE can be estimated from the data collected from an observational study:

$$\begin{split}&\mathbb{E}\left(Y^{A=1}\mid X=x,U=u\right)-\mathbb{E}\left(Y^{A=0}\mid X=x,U=u\right)\\ &=\mathbb{E}\left(Y^{A=1}\mid A=1,X=x,U=u\right)-\mathbb{E}\left(Y^{A=0}\mid A=0,X=x,U=u\right)\\ &=\mathbb{E}\left(Y\mid A=1,X=x,U=u\right)-\mathbb{E}\left(Y\mid A=0,X=x,U=u\right).\end{split}$$

In clustered observational study, a common approach to estimate $\mathbb{E}\left(Y\mid A,X=x,U=u\right)$ is to use a generalized mixed effects model or a generalized linear model with GEE, where all potential confounders are included. We will introduce our sensitivity analysis method based on generalized mixed effects model, where random effects will be used to account for within-cluster correlation. Let $i$ $(i=1,\dots,I)$ and $j$ $(j=1,\dots,J)$ be the indicators for the level-1 (e.g., subjects nested in clusters) and the level-2 unit (e.g., clusters), respectively. Consider a continuous outcome $Y_{ij}$ with the marginal mean model:

$$\small\begin{split}\mathbb{E}(Y_{ij}|A_{ij},X_{ij},U_{ij})=\beta_{0}+\beta_{1}A_{ij}+\beta_{2}X_{ij}+\beta_{3}A_{ij}X_{ij}+\theta U_{ij};\end{split}$$

and the corresponding mixed model:

$$\small\begin{split}&Y_{ij}=\beta_{0}+\beta_{1}A_{ij}+\beta_{2}X_{ij}+\beta_{3}A_{ij}X_{ij}+\theta U_{ij}+\zeta_{j}+\epsilon_{ij},\end{split}$$

(1)

where $\boldsymbol{\beta}=(\beta_{0},\beta_{1},\beta_{2},\beta_{3})^{T}$ and $\theta$ are unknown regression parameters; $\zeta_{j}$ and $\epsilon_{ij}$ are the random intercepts and the level-1 error terms, respectively, where $\zeta_{j}\sim N(0,\nu)$ and $\epsilon_{ij}|\zeta_{j}\sim N(0,\phi)$ with unknown variances $\nu$ and $\phi$.

Our goal is to estimate the CATE, denoted as $E_{x}$, with the form

$$\begin{split}E_{x}\coloneqq\mathbb{E}\left(Y\mid A=1,X=x,U=u\right)-\mathbb{E}\left(Y\mid A=0,X=x,U=u\right)=\beta_{1}+x\beta_{3}.\end{split}$$

Since $U$ is unmeasured, we can only fit the following confounded model:

$$Y_{ij}=\beta_{0}^{c}+\beta_{1}^{c}A_{ij}+\beta_{2}^{c}X_{ij}+\beta_{3}^{c}A_{ij}X_{ij}+\zeta_{j}^{c}+\epsilon_{ij}^{c},$$

(2)

where $\boldsymbol{\beta}^{c}=(\beta_{0}^{c},\beta_{1}^{c},\beta_{2}^{c},\beta_{3}^{c})^{T}$ are unknown regression parameters; $\zeta_{j}^{c}$ and $\epsilon_{ij}^{c}$ are random intercepts and level-1 error terms, respectively, where $\zeta_{j}^{c}\sim N(0,\nu^{c})$ and $\epsilon_{ij}^{c}|\zeta_{j}^{c}\sim N(0,\phi^{c})$ with unknown variances $\nu^{c}$ and $\phi^{c}$. There is potential violation of assuming a constant variance of random effect and error term, but our methods are robust against this potential violation based on the simulation results.

Denote confounded treatment effect $E_{x}^{c}$ as:

$$E_{x}^{c}\coloneqq\mathbb{E}\left(Y\mid A=1,X=x\right)-\mathbb{E}\left(Y\mid A=0,X=x\right)=\beta_{1}^{c}+x\beta_{3}^{c}.$$

We refer $\boldsymbol{\beta}$ and $\boldsymbol{\beta}^{c}$ as the true and confounded coefficients, respectively. The element in $\boldsymbol{\beta}$ cannot be estimated directly using the observed data while $\boldsymbol{\beta}^{c}$ can. Therefore, our goal is to establish the relationship between $\boldsymbol{\beta}$ and $\boldsymbol{\beta}^{c}$ to infer $\boldsymbol{\beta}$.

Let $F(U|A,X)$ be the cumulative distribution function of $U$ given $A$ and $X$. The conditional expectation of $Y_{ij}$ given $A_{ij}$ and $X_{ij}$ based on model (LABEL:1) becomes:

$$\begin{split}\mathbb{E}(Y_{ij}|A_{ij},X_{ij}=x)&=\mathbb{E_{U}}{\mathbb{E}_{\zeta}\mathbb{E}(Y_{ij}|A_{ij},X_{ij}=x,U_{ij},\zeta_{j})}\\ &=\mathbb{E_{U}}\mathbb{E}(Y_{ij}|A_{ij},X_{ij}=x,U_{ij},\zeta_{j}=0)\\ &=\beta_{0}+\beta_{1}A_{ij}+x\beta_{2}+x\beta_{3}A_{ij}+\theta\int_{-\infty}^{\infty}U_{ij}dF(U_{ij}|A_{ij},X_{ij}=x)\\ &=(\beta_{0}+\theta m_{0x})+\{\beta_{1}+x\beta_{3}+\theta(m_{1x}-m_{0x})\}A_{ij}+x\beta_{2},\end{split}$$

(3)

where $m_{1x}=\mathbb{E}(U_{ij}|A_{ij}=1,X_{ij}=x)$ and $m_{0x}=\mathbb{E}(U_{ij}|A_{ij}=0,X_{ij}=x)$.

Denote the estimated confounded treatment effect as $\widehat{E}_{x}^{c}$ (95% CI, $lb$ to $ub$), where $lb$ and $ub$ are the lower and upper bounds of its 95% confidence interval. Following equation (5), the estimated conditional average treatment effects is $\widehat{E}_{x}^{c}-B_{x}$ (95% CI, $lb-B_{x}$ to $ub-B_{x}$). We say an unmeasured confounder can explain away the estimated confounded treatment effect if the unmeasured confounder is able to shift the confidence interval of the estimated confounded treatment effect to include the null value (no treatment effect). If confidence interval of estimated confounded treatment effect includes the null value, then no confounding is needed to explain away the results or move confidence interval to include the null value.

Once the sensitivity parameters are specified, we can calculate bias factor and the estimated CATE. However, in practice, we do not know the strengths and distributions of unmeasured confounders, thus, the sensitivity parameters are usually unknown. What we can do is to calculate how the estimate is affected under different values of the bias factor. Among all different bias factors that are able to explain away the confounded treatment effect, minimal bias factor is defined as the smallest value of bias factors that is able to explain away the confounded treatment effect, meaning any bias factor larger than the minimal bias factor can also explain away the confounded treatment effect. The robustness of the estimated confounded treatment effect can be measured by the minimal bias factor. A small value of the minimal bias factor indicates the estimated treatment effect is not robust to unmeasured confounders.

Let $\widehat{E}_{x}^{c}>0$ and $\widehat{E}_{x}^{c}<0$ be apparently positive effect and apparently negative effect, respectively. In apparently positive effect case, the confounded treatment effect will be explained away when the lower bound of confidence interval of estimated confounded treatment effect is smaller than 0, i.e., $lb-{B}_{x}\leq 0$. It means any bias factor, $B_{x}$, such that $B_{x}\geq lb$, can explain away the confounded treatment effect. Denote $B_{x}^{*}$ as the minimal bias factor that can explain away the confounded treatment effect. Thus, when $\widehat{E}_{x}^{c}>0$, the minimal bias factor that can explain away the confounded treatment effect is $lb$, i.e., ${B}_{x}^{*}=lb$.

In the apparently negative case, define $\tilde{B}_{x}=-{B}_{x}$, thus, estimated CATE can be expressed as $\widehat{E}_{x}^{c}+\tilde{B}_{x}$ (95% CI, $lb+\tilde{B}_{x}$ to $ub+\tilde{B}_{x}$). The confounded treatment effect will be explained away when $ub+\tilde{B}_{x}\geq 0$. It means any bias factor, $\tilde{B}_{x}$, such that $\tilde{B}_{x}\geq-ub$, can explain away the confounded treatment effect. Denote $\tilde{B}_{x}^{*}$ as the minimal bias factor that can explain away the confounded treatment effect. Thus, when $\widehat{E}_{x}^{c}<0$, the minimal bias factor that can explain away the confounded treatment effect is $-ub$, i.e., $\tilde{B}_{x}^{*}=-ub$.

As shown above, bias factor is a function of the effects of unmeasured confounders, and the difference of mean of unmeasured confounders between treated and control groups given $X=x$, such that $B_{x}=(m_{1x}-m_{0x})\theta$ and $\tilde{B}_{x}=-(m_{1x}-m_{0x})\theta$. It means that, when $\widehat{E}_{x}^{c}>0$, a set of $\theta$, $m_{1x}$, and $m_{0x}$, such that $(m_{1x}-m_{0x})\theta\geq lb$, can explain away estimated confounded treatment effect. When $\widehat{E}_{x}^{c}<0$, a set of $\theta$, $m_{1x}$, and $m_{0x}$, such that $(m_{0x}-m_{1x})\theta\geq-ub$, can explain away estimated confounded treatment effect. We recommend to provide a table or contour plot which present different values of $\theta$, $m_{1x}$, $m_{0x}$, and bias factor to give researchers a sense of what are combinations of sensitivity parameters that can explain away the confounded treatment effect and therefore a sense of how sensitive the confounded treatment effect is to unmeasured confounders.

For a binary outcome, risk difference (RD) and relative risk (RR) can be used to explore causal treatment effect. In this article, RR will be used to explore the causal relationship. With same notations as in the case of continuous outcome, except that $Y$ denotes a binary variable now, when same assumptions hold:

$$\begin{split}\text{RR}=\frac{\mathbb{E}\left(Y^{A=1}\mid X=x,U=u\right)}{\mathbb{E}\left(Y^{A=0}\mid X=x,U=u\right)}=\frac{\mathbb{E}\left(Y^{A=1}\mid A=1,X=x,U=u\right)}{\mathbb{E}\left(Y^{A=0}\mid A=0,X=x,U=u\right)}=\frac{\mathbb{E}\left(Y\mid A=1,X=x,U=u\right)}{\mathbb{E}\left(Y\mid A=0,X=x,U=u\right)}.\end{split}$$

When the event is rare in both treatment and control group, it can be shown that odds ratio is approximately equal to the relative risk. A logistic regression model is widely used to estimate odds ratio. Let the response variable $Y$ be related to $A$, $X$, and $U$ through following generalized linear mixed model:

$$\text{logit}\left\{\mathbb{P}(Y_{ij}=1|A_{ij},X_{ij},U_{ij},\zeta_{j})\right\}=\beta_{0}+\beta_{1}A_{ij}+\beta_{2}X_{ij}+\beta_{3}A_{ij}X_{ij}+\theta U_{ij}+\zeta_{j},$$

(6)

where $\boldsymbol{\beta}=(\beta_{0},\beta_{1},\beta_{2},\beta_{3})^{T}$ and $\theta$ are unknown regression parameters, and $\zeta_{j}\sim N(0,\nu)$.

When $v$ is small, it can be shown (proof is provided in Appendix A)

$$\text{logit}\{\mathbb{P}(Y_{ij}=1|A_{ij},X_{ij},U_{ij})\}\approx\beta_{0}+\beta_{1}A_{ij}+\beta_{2}X_{ij}+\beta_{3}A_{ij}X_{ij}+\theta U_{ij}.$$

Denote CATE on the log-RR scale as $E_{x}$. Thus, when the event is rare and $v$ is small, $E_{x}$ has the form: $E_{x}\coloneqq\text{log}\frac{\mathbb{E}\left(Y^{A=1}\mid X=x,U=u\right)}{\mathbb{E}\left(Y^{A=0}\mid X=x,U=u\right)}=\text{log}\frac{\mathbb{E}\left(Y\mid A=1,X=x,U=u\right)}{\mathbb{E}\left(Y\mid A=0,X=x,U=u\right)}\approx\beta_{1}+x\beta_{3}.$

Since $U$ is unmeasured, we can only fit the following confounded model:

$$\text{logit}\left\{\mathbb{P}(Y_{ij}=1|A_{ij},X_{ij},\zeta_{j}^{c})\right\}=\beta_{0}^{c}+\beta_{1}^{c}A_{ij}+\beta_{2}^{c}X_{ij}+\beta_{3}^{c}A_{ij}X_{ij}+\zeta_{j}^{c},$$

(7)

where $\zeta_{j}^{c}\sim N(0,\nu^{c})$ and $\boldsymbol{\beta}^{c}=(\beta_{0}^{c},\beta_{1}^{c},\beta_{2}^{c},\beta_{3}^{c})^{T}$ are unknown regression parameters. Similarly, when $\nu^{c}$ is small, it can be shown that for the confounded model:

$$\text{logit}\left\{\mathbb{P}(Y_{ij}=1|A_{ij},X_{ij})\right\}\approx\beta_{0}^{c}+\beta_{1}^{c}A_{ij}+\beta_{2}^{c}X_{ij}+\beta_{3}^{c}A_{ij}X_{ij}.$$

(8)

Denote confounded treatment effect on the log-RR scale as $E_{x}^{c}$ which has the form: $E_{x}^{c}\coloneqq\text{log}\frac{\mathbb{E}\left(Y\mid A=1,X=x\right)}{\mathbb{E}\left(Y\mid A=0,X=x\right)}\approx\beta_{1}^{c}+x\beta_{3}^{c}.$

When $U$ is binary and both $\theta$ and $\nu$ are small, we have (proof is provided in Appendix B):

$$\small\begin{split}&\text{logit}\left\{\mathbb{P}(Y=1|A,X=x)\right\}\approx\beta_{0}+\beta_{1}A+\beta_{2}x+\beta_{3}Ax+\text{log}\{P_{Ax}\text{exp}(\theta+\frac{\nu}{2})+(1-P_{Ax})\text{exp}(\frac{\nu}{2})\}\\ &=\beta_{0}+\beta_{2}x+\text{log}\{P_{0x}\text{exp}(\theta+\frac{\nu}{2})+(1-P_{0x})\text{exp}(\frac{\nu}{2})\}+\{\beta_{1}+\beta_{3}x+\text{log}\frac{P_{1x}\text{exp}(\theta)+1-P_{1x}}{P_{0x}\text{exp}(\theta)+1-P_{0x}}\}A,\end{split}$$

(9)

where $P_{AX}=P(U_{ij}=1|A_{ij},X_{ij})$.

Following the confounded model described in equation (8), we have:

$$\text{logit}\left\{\mathbb{P}(Y=1|A,X=x)\right\}\approx\beta_{0}^{c}+\beta_{1}^{c}A_{ij}+\beta_{2}^{c}x+\beta_{3}^{c}Ax=\beta_{0}^{c}+\beta_{2}^{c}x+(\beta_{1}^{c}+\beta_{3}^{c}x)A.$$

(10)

When $\theta$, $\nu$, and $\nu^{c}$ are small, comparing equations (9) and (10), we have the following relationship:

$$E_{x}=E_{x}^{c}-B_{x},$$

(11)

where $E_{x}=\beta_{1}+x\beta_{3}$, $E_{x}^{c}=\beta_{1}^{c}+x\beta_{3}^{c}$, and $B_{x}=\text{log}\frac{P_{1x}\text{exp}(\theta)+1-P_{1x}}{P_{0x}\text{exp}(\theta)+1-P_{0x}}$.

When $U$ is normally distributed such that $U\sim N(\mu_{AX},\sigma_{u})$, i.e., the expectation of $U$ is dependent on both $A$ and $X$, then if $\theta$, $\sigma_{u}$, and $\nu$ are small, we have (proof is provided in Appendix C):

$$\begin{split}\text{logit}\left\{\mathbb{P}(Y=1|A,X=x)\right\}&\approx\beta_{0}+\beta_{1}A+\beta_{2}x+\beta_{3}Ax+\theta\mu_{Ax}+\frac{\theta^{2}\sigma_{u}+\nu}{2}\\ &=\beta_{0}+\beta_{2}x+\mu_{0x}+\frac{\theta^{2}\sigma_{u}+\nu}{2}+(\beta_{1}+x\beta_{3}+\theta\mu_{1x}-\theta\mu_{0x})A.\end{split}$$

(12)

When $\theta$, $\sigma_{u}$, $\nu$, and $\nu^{c}$ are small, comparing equations (10) and (12), we have the following relationship:

$$E_{x}=E_{x}^{c}-B_{x},$$

(13)

where $E_{x}=\beta_{1}+x\beta_{3}$, $E_{x}^{c}=\beta_{1}^{c}+x\beta_{3}^{c}$, and $B_{x}=\theta(\mu_{1x}-\mu_{0x})$.

Intraclass correlation coefficient (ICC) is the proportion of total variance of the outcome that is attributable to the between-subject variation, i.e., how much variation of the outcome can be explained by the random effects. For a logistic model, $\text{ICC}=\frac{\nu}{\nu+\pi/3}.$ It means when $\theta$, ICC, and confounded ICC are small, equations (11) and (13) hold.

With equations (11) and (13), we can calculate the minimal bias factor that can explain away the confounded treatment effect by moving the confidence interval to include the null value. Similarly, if a small minimal bias factor is able to explain away the confounded treatment effect, the confounded treatment effect is not robust to unmeasured confounders.

When multiple observational studies is involved where only study-level estimates not individual-level data can be merged, meta-analysis is the most commonly used analysis technique to estimate the population covariate effects which combines estimates from all studies and takes study heterogeneous into consideration. In this section, we introduce our proposed method to assess how sensitive the estimated population covariate effects from meta-analysis are to the unmeasured confounders.

Suppose we have $K$ studies which are random samples from a population of interest. For study $k$, let $Y_{ijk}$ be the continuous response variable, which is related to treatment assignment $A_{ijk}$, measured confounders $X_{ijk}$, and unmeasured confounder $U_{ijk}$ through the following random intercepts model:

$$Y_{ijk}=\beta_{0k}+\beta_{1k}A_{ijk}+\beta_{2k}X_{ijk}+\beta_{3k}A_{ijk}X_{ijk}+\theta_{k}U_{ijk}+\zeta_{jk}+\epsilon_{ijk},$$

where $i$ $(i=1,\dots,I)$ and $j$ $(j=1,\dots,J)$ are the indicators for the level-1 (e.g.,subjects in clusters) and the level-2 unit (e.g., clusters), respectively; and $\beta_{k}$’s and $\theta_{k}$ are the regression parameters adjusting for both measured and unmeasured confounders. Similarly, when causal assumptions hold, the model can be used to estimate the causal ATE. Denote $E_{k|x}$ as the CATE for study $k$ given $X=x$, such that $E_{k|x}=\beta_{1k}+x\beta_{3k}$.

Since $U_{ijk}$ is unmeasured, we can only fit the following confounded model for each of the $K$ studies:

$$Y_{ijk}=\beta_{0k}^{c}+\beta_{1k}^{c}A_{ijk}+\beta_{2k}^{c}X_{ijk}+\beta_{3k}^{c}A_{ijk}X_{ijk}+\zeta_{jk}^{c}+\epsilon_{ijk}^{c},$$

where $\beta_{k}^{c}$’s are confounded regression parameters only adjusting for measured confounders. Denote $E_{k|x}^{c}$ as the confounded treatment effect for study $k$ given $X=x$, such that $E_{k|x}^{c}=\beta_{1k}^{c}+x\beta_{3k}^{c}$.

In a standard setup of parametric random effects meta-analysis (Borenstein et al., 2007), it assumes that each of the $K$ studies measures a potentially unique confounded effect size $E_{k|x}^{c}$, such that $E_{k|x}^{c}\sim N(\mu_{E_{x}^{c}},V_{E_{x}^{c}})$, for $k=1,\dots,K$. Let $\widehat{E}_{k|x}^{c}$ be the point estimate of the treatment effect and $\nu_{k}$ be the known within-study variance for study $k$. In random effects meta-analysis setup, $\widehat{E}_{k|x}^{c}$ is determined by the true confounded effect size $E_{k|x}^{c}$ plus the known within-study variance $\nu_{k}$, while $E_{k|x}^{c}$ is determined by the grand mean, $\mu_{E_{x}^{c}}$ and the between-study variance, $V_{E_{x}^{c}}$. More generally, for any observed effect $\widehat{E}_{k|x}^{c}$,

$$\widehat{E}_{k|x}^{c}=\mu_{E_{x}^{c}}+\varphi_{k}+\kappa_{k},$$

where $\varphi_{k}\sim N(0,V_{E_{x}^{c}})$ and $\kappa_{k}\sim N(0,\nu_{k})$. Fixed effects meta-analysis can be viewed as a special case where $V_{E_{x}^{c}}=0$.

The relationship between the CATE and confounded treatment effects in equation (5) still hold for each of the $K$ studies as we illustrated in the single study case,

$$\displaystyle E_{k|x}=E_{k|x}^{c}-B_{k|x};\hskip 36.135ptB_{k|x}=(m_{k|1x}-m_{k|0x})\theta_{k},$$

where $m_{k|1x}=\mathbb{E}(U_{ijk}|A_{ijk}=1,X_{ijk}=x)$ and $m_{k|0x}=\mathbb{E}(U_{ijk}|A_{ijk}=0,X_{ijk}=x)$, $B_{k|x}\sim N(\mu_{B_{x}},V_{B_{x}})$ and $E_{k|x}^{c}\sim N(\mu_{E_{x}^{c}},V_{E_{x}^{c}})$. It means that we allow the bias factor $B_{k|x}$ to vary across studies. Assume that $B_{k|x}$ is independent of $E_{k|x}$. $B_{k|x}$ is independent of $E_{k|x}$ but not $E_{k|x}^{c}$ because studies with a larger bias factor may observe a larger confounded treatment effect. Thus, the distribution of the CATE can be written as $E_{k|x}\sim N(\mu_{E_{x}},V_{E_{x}})$, where

$$\mu_{E_{x}}=\mu_{E_{x}^{c}}-\mu_{B_{x}};\hskip 36.135ptV_{E_{x}}=V_{E_{x}^{c}}-V_{B_{x}}.$$

Although we use continuous outcome as an example to demonstrate the method in multiple studies, similar derivation can be easily extended to the binary outcome case by replacing the error distribution and the link function in the generalized linear mixed model for binary outcomes. Given equations (11) and (13) hold for each of the $K$ studies, we can similarly derive the distribution of CATE under the case of binary outcomes.

Following the newly designed metrics for meta-analyses of heterogeneous effects in Mathur and VanderWeele (2019a) and Mathur and VanderWeele (2019b), we define that there is a treatment effect if the probability of observing a CATE larger than a scientifically meaningful effect size $q$, is larger than a threshold $r$. This means that a bias factor will explain away the confounded treatment effect in meta-analysis if the bias factor can reduce the probability of observing a conditional average treatment effects larger than $q$ to less than $r$. Note, for continuous outcome, $q$ is on the difference of the continuous outcomes scale, while for binary outcome, $q$ is on the log-RR scale.

Denote $\widehat{\mu}_{E_{x}^{c}}$ and $\widehat{V}_{E_{x}^{c}}$ as the estimators of $\mu_{E_{x}^{c}}$ and $V_{E_{x}^{c}}$. Let $\widehat{\mu}_{E_{x}^{c}}>0$ and $\widehat{\mu}_{E_{x}^{c}}<0$ be the apparently positive effect and apparently negative effect, respectively, in meta-analysis.

For the apparently positive effect, we define the probability of observing a CATE larger than $q$ as $p(q)_{x}=P\left(E_{k|x}>q\right)$ where $E_{k|x}\sim N(\mu_{E_{x}},V_{E_{x}})$. $p(q)_{x}$ can be expressed using $\mu_{B_{x}}$, $V_{B_{x}}$, $\mu_{E_{x}^{c}}$, and $V_{E_{x}^{c}}$:

$${p}(q)_{x}=1-\Phi\left(\frac{q+\mu_{B_{x}}-\mu_{E_{x}^{c}}}{\sqrt{V_{E_{x}^{c}}-V_{B_{x}}}}\right),\quad V_{E_{x}^{c}}>V_{B_{x}}$$

(14)

where $\Phi$ is the cumulative function of standard normal distribution. $\mu_{B_{x}}$ and $V_{B_{x}}$ such that $1-\Phi\left(\frac{q+\mu_{B_{x}}-\widehat{\mu}_{E_{x}^{c}}}{\sqrt{\widehat{V}_{E_{x}^{c}}-V_{B_{x}}}}\right)<r$ will explain away the estimated confounded treatment effect in meta-analysis. When $0<r<0.5$, a constant bias factor among studies, i.e., $V_{B_{x}}=0$, will lead to an upper bound of $\widehat{p}(q)_{x}$, i.e.,

$$1-\Phi\left(\frac{q+\mu_{B_{x}}-\widehat{\mu}_{E_{x}^{c}}}{\sqrt{\widehat{V}_{E_{x}^{c}}-V_{B_{x}}}}\right)\leq 1-\Phi\left(\frac{q+\mu_{B_{x}}-\widehat{\mu}_{E_{x}^{c}}}{\sqrt{\widehat{V}_{E_{x}^{c}}}}\right),\ \ 0<r<0.5.$$

It means a constant bias factor, such that $B_{X}\geq B_{x}^{*}$ where

$$1-\Phi\left(\frac{q+B_{x}^{*}-\widehat{\mu}_{E_{x}^{c}}}{\sqrt{\widehat{V}_{E_{x}^{c}}}}\right)\equiv r\rightarrow B_{x}^{*}=\Phi^{-1}(1-r)\widehat{V}_{E_{x}^{c}}^{\frac{1}{2}}-q+\widehat{\mu}_{E_{x}^{c}},$$

(15)

and any other non-constant bias factors with mean larger or equal to $B_{x}^{*}$ can explain away the confounded treatment effect in meta-analysis. Thus, $B_{x}^{*}$ can be viewed as the minimal common constant bias factor that reduce the probability of observing a CATE larger than a scientifically meaningful effect size $q$, to a threshold $r$.

For the apparently negative case, define $\tilde{B}_{x}=-B_{x}$; therefore $E_{k|x}=E_{k|x}^{c}+\tilde{B}_{k|x}$ and $E_{k|x}\sim N(\mu_{E_{x}},V_{E_{x}})$ where $\mu_{E_{x}}=\mu_{E_{x}^{c}}+\mu_{\tilde{B}_{x}}$ and $V_{E_{x}}=V_{E_{x}^{c}}-V_{\tilde{B}_{x}}$. We define the probability of observing a CATE larger than $q$ as $p(q)_{x}=P\left(E_{k|x}<q\right)$. $p(q)_{x}$ can be expressed using $\mu_{\tilde{B}_{x}}$, $V_{\tilde{B}_{x}}$, $\mu_{E_{x}^{c}}$, and $V_{E_{x}^{c}}$:

$${p}(q)_{x}=\Phi\left(\frac{q-\mu_{\tilde{B}_{x}}-\mu_{E_{x}^{c}}}{\sqrt{V_{E_{x}^{c}}-V_{\tilde{B}_{x}}}}\right),\quad V_{E_{x}^{c}}>V_{\tilde{B}_{x}},$$

Similarly, when $0<r<0.5$, a constant bias factor among studies will lead to an upper bound of $\widehat{p}(q)_{x}$. Let $\tilde{B}_{x}^{*}$ be the minimal common constant bias factor that can reduce the probability of observing a conditional average treatment effects larger than a scientifically meaningful effect size $q$, to a threshold $r$ such that

$$\Phi\left(\frac{q-\tilde{B}_{x}^{*}-\widehat{\mu}_{E_{x}^{c}}}{\sqrt{\widehat{V}_{E_{x}^{c}}}}\right)\equiv r\rightarrow\tilde{B}_{x}^{*}=q-\Phi^{-1}(r)\widehat{V}_{E_{x}^{c}}^{\frac{1}{2}}-\widehat{\mu}_{E_{x}^{c}}.$$

(16)

Consider fixed $m_{k|1x}$, $m_{k|0x}$, $\theta_{k}$ and therefore a fixed bias factor across studies, meaning the unmeasured confounder has the same distribution and effect across all $K$ studies. In the apparently positive effect case, this means that if each of the $K$ studies suffers from unmeasured confounder with the same strength such that $(m_{k|1x}-m_{k|0x})\theta_{k}\geq\Phi^{-1}(1-r)\widehat{V}_{E_{x}^{c}}^{\frac{1}{2}}-q+\widehat{\mu}_{E_{x}^{c}}$, the confounded treatment effect will be explained away. In the apparently negative effect case, it means that if $(m_{k|0x}-m_{k|1x})\theta_{k}\geq q-\Phi^{-1}(r)\widehat{V}_{E_{x}^{c}}^{\frac{1}{2}}-\widehat{\mu}_{E_{x}^{c}}$ in each of the $K$ studies, the confounded treatment effect will be explained away.

Consider variant $m_{k|1x}$, $m_{k|0x}$, $\theta_{k}$ and therefore a variant bias factor across studies. In apparently positive effect case, any random bias factor with mean equal or larger than $\Phi^{-1}(1-r)\widehat{V}_{E_{x}^{c}}^{\frac{1}{2}}-q+\widehat{\mu}_{E_{x}^{c}}$ can explain away the confounded treatment effect. In apparently negative effect case, any random bias factor with mean equal or larger than $q-\Phi^{-1}(r)\widehat{V}_{E_{x}^{c}}^{\frac{1}{2}}-\widehat{\mu}_{E_{x}^{c}}$ can explain away the confounded treatment effect.

Similar to the case of the single study, the robustness of the estimated confounded treatment effect can be measured by the minimal common constant bias factor. A small minimal common constant bias factor indicates that the estimated confounded treatment effect is not robust to unmeasured confounders.

Table 1 summarizes how to calculate the minimum (common constant) bias factor for a single study and for multiple studies (meta-analysis) under different scenarios. Once $B_{x}^{*}$ or $\tilde{B}_{x}^{*}$ is calculated, one can determine the minimum effect of unmeasured confounders on response and difference between conditional mean of unmeasured confounders, that are sufficient to explain away the confounded treatment effect.

We take a prospective clustered study, which focuses on identifying potential causal effect of breastfeeding on children’s intelligence quotient (IQ), as an example to explain how to apply our method to assess the robustness of results from an observational study. Luby et al. (2016) examined relationship between breastfeeding and children’s IQ. After adjusting for children’s age, sex, history of depression, caregiver’s highest level of education, and familial income-to-needs ratio, the authors found that the IQ of breastfed child was $5.49$ (95% CI: $0.75$ to $10.23$) higher than the IQ of non-breastfed child. The investigators did not control for maternal smoking during pregnancy while maternal smoking may be a confounder as it could associate with less breastfeeding as well as lower offspring IQ. We are worried this result is confounded by maternal smoking status, so we conduct a sensitivity analysis to assess how robust the estimate is to the unmeasured confounding: maternal smoking status.

As the study indicated an apparently positive effect of breastfeeding on children’s IQ, the minimal bias factor, $B_{x}^{*}$, that is able to explain away the confounded effect equals the lower bound of confidence interval, i.e., $B_{x}^{*}=0.75$, which means that any bias factor larger than $0.75$ can also explain away the confounded effect. The minimum bias factor $0.75$ can be decomposed into the product of two factors: the effect of maternal smoking on children’s IQ adjusting for other covariates and the difference between $\mathbb{E}(U|A=1,X=x)$ and $\mathbb{E}(U|A=0,X=x)$. Our methods allow investigators to make statements of the following form: the confounded effect of breastfeeding on children’s IQ could be explained away if the difference between the prevalence of maternal smoking in breastfed children and non-breastfed children is larger than $0.25$ and the effect of smoking on IQ is greater than 3. There are also other decompositions that could explain away the confounded effect of breastfeeding on children’s IQ. A contour plot depicted in Figure 1 can be used to visualize the relationship between bias factors, difference between $\mathbb{E}(U|A=1,X=x)$ and $\mathbb{E}(U|A=0,X=x)$ (i.e., $m_{1x}-m_{0x}$), and the effect of unmeasured confounder $\theta$. Any combination of $m_{1x}-m_{0x}$ and $\theta$ in the contour plot that corresponds to a bias factor equal or larger than $0.75$ are able to explain away the confounded effect. Using both the information of the minimal bias factor and different combinations of sensitivity parameters that are able to explain away the confounded treatment effect, investigators can get a sense of how sensitive the conclusions are to potential unmeasured confounders.

Applications of our methods for meta-analysis and single study are similar. Investigators can calculate the minimal common constant bias factor based on equations (15) or (16) by plugging in $\widehat{V}_{E_{x}^{c}}$, $\widehat{\mu}_{E_{x}^{c}}$, $r$, and $q$. If each of the studies suffers from an unmeasured confounder with the same distribution and effect on outcome that lead to a bias factor equal or greater than the minimal common constant bias factor, the confounded treatment effect from the meta-analysis can be explained away. Further, any non-fixed unmeasured confounder across studies with expectation equal or larger than minimal common constant bias factor can also explain away the confounded effect. Similar to single study, investigators can use our method to assess whether a set of sensitivity parameters can explain away the confounded effect and get a sense of how sensitive the conclusions are to potential unmeasured confounders.

We take a meta-analysis, which focuses on identifying potential causal effect of body mass index (BMI) in midlife on dementia, as an example to explain how to apply our method to assess the robustness of results from meta-analysis. Albanese et al. (2017) examined relationship between BMI in midlife and dementia by combining results, on relative risk (RR) scale, from 19 longitudinal studies on 589,649 participants (2,040 incident dementia cases) followed up for up to 42 years. The authors found that midlife (age 35 to 65 years) obesity (BMI $\geq$ 30) was associated with larger risk of dementia in late life (RR, 1.33; 95% CI, 1.08 to 1.63). The study reported $\widehat{\mu}_{E_{x}^{c}}=\text{log}1.33$ ($\widehat{\mu}_{E_{x}^{c}}$ is on log RR scale), while the study level summary measures reported from Albanese et al. (2017) can be used to obtain $\widehat{V}_{E_{x}^{c}}$. We estimated $\widehat{V}_{E_{x}^{c}}=0.08$ via methods illustrated in Borenstein et al. (2007). Note that the choice of $q$ and $r$ should be adapted based on the context and specific research questions, which means that different analysis can have difference choice on $q$ and $r$. In this example, let us consider a RR of 1.2 as a scientifically meaningful effect size, meaning $q=\text{log}1.2$, and consider having at least 40% chance to observe an effect size above $q$ for results to be of scientific interest, meaning $r=0.4$.

As the meta-analysis indicated an apparently positive effect of BMI on dementia, we have the minimal common constant bias factor, which can reduce the probability of observing a conditional average treatment effects larger than a scientifically meaningful effect size $q$ by a threshold $r$, such that

$$B_{x}^{*}=\Phi^{-1}(1-r)\widehat{V}_{E_{x}^{c}}^{\frac{1}{2}}-q+\widehat{\mu}_{E_{x}^{c}}=\Phi^{-1}(1-0.4)*\sqrt{0.08}-\text{log}1.2+\text{log}1.33=0.17.$$

It means that any constant bias factor across studies such that $B_{x}\geq 0.17$ and any other non-constant bias factors across studies with mean larger or equal to $0.17$ can explain away the confounded treatment effect in meta-analysis. The minimum common constant bias factor $0.17$ can be decomposed into the product of two factors: the effect of unmeasured confounder on dementia and the difference between the mean of unmeasured confounder in different BMI groups. Our methods allow investigators to make statement of the following form: the confounded effect of BMI on dementia could be explained away if an unmeasured confounder exists such that its mean difference between BMI groups is larger than 0.17 and its effect on dementia is greater than 1 in each of the studies. We can similarly refer to the contour plot in Figure 1 to visualize combinations of sensitivity parameters that corresponds to a bias factor equal or larger than 0.17 to explain away the confounded effect. Using both the information of the minimal common constant bias factor and different combination of sensitivity parameters that are able to explain away the confounded treatment effect, investigators can get a sense of how sensitive the conclusions are to potential unmeasured confounders.

First, we assessed the performance of the estimator of $E_{x}$, $\widehat{E}_{x}=\widehat{\beta}_{1}^{c}+x\widehat{\beta}_{3}^{c}-B_{x}$. We assumed that $A$, $X$, and $U$ are dependent, which is a frequently encountered condition in observational studies.

Let the number of subjects be $J$ and the number of repeated measures be $I$ per subject. We simulated a normally distributed unmeasured confounder $U_{ij}\sim N(0,\sigma_{u}^{2})$, one measured binary confounder $X_{ij}\sim Bin(0.5)$ if $U_{ij}<1$ and $X_{ij}\sim Bin(0.4)$ if $U_{ij}\geq 1$, and treatment assignment $A_{ij}\sim Bin(0.4)$ if $U_{ij}+X_{ij}<2$ and $A_{ij}\sim Bin(0.5)$ if $U_{ij}+X_{ij}\geq 2$. We simulated random intercept $\zeta_{j}$ and the level-1 residuals $\epsilon_{ij}$ from the distributions $N(0,4)$ and $N(0,1)$, respectively. We proceeded to simulate the response variable $Y$ using the underlying true model described in (LABEL:1) setting $\beta_{0}=1$ and $\beta_{2}=3$ while varying the values of $\beta_{1}$, $\beta_{3}$, and $\theta$. $\widehat{E}_{x}$ can be calculated after we estimated the confounded coefficient $\beta_{1}^{c}$ and $\beta_{3}^{c}$ in (2) using R package lme4 (Bates et al., 2015). We simulated 1,000 such studies and inferred the estimated conditional average treatment effects described in equation (5) under $X=0$ and $X=1$. We repeated the simulations with different values of $I$, $J$, $\beta_{1}$, $\beta_{3}$, $\theta$, and $\sigma_{u}^{2}$.

To measure the performance of our methods in estimating conditional causal treatment effect, we assessed the bias, standard error (SE) and coverage probability (CP) of 95% confidence intervals of $\widehat{E}_{x}$. Results in Table 2 show that the estimated conditional average treatment effects, $\widehat{E}_{x}$, is unbiased and the standard error is relatively small. The coverage of 95% confidence interval is nominal. The accuracy of the point estimate appears to be independent of $X$.

Then, we assessed the performance of the estimator of $E_{x}$ in equation (13) for a single study setting with binary outcome and continuous unmeasured confounder. Assuming 200 subjects and the 4 repeated measures per subject, $A$, $X$, and $U$ were simulated similarly as for continuous outcome. We simulated random intercept $\zeta_{j}$ and the level-1 residuals $\epsilon_{ij}$ from the distributions $N(0,\nu)$ and $N(0,1)$, respectively. We proceeded to simulate the response variable $Y$ using the underlying true model described in (6) setting $\beta_{0}=-4.5$, $\beta_{1}=1$, $\beta_{2}=3$, and $\beta_{3}=-0.5$ while varying the value of $\theta$.