Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
88,113 characters · 11 sections · 27 citation commands
Causal Simulation Experiments: Lessons from Bias Amplification
\runninghead{Stokes et al}
\corrauth{Russell Steele, McGill University Department of Mathematics and Stistics, Burnside Hall, Room 1005 805 Sherbrooke Street West Montreal, Quebec Canada H3A 0B9}
\email{[email removed]}
\setcounter{page}{1}
Causal identification strategies aim to condition on a sufficient set of observables such that the potential outcomes are conditionally independent of the treatment of interest rubin1974estimating, rosenbaum1983central,wooldridge2010econometric. Causal variable selection procedures often assume that at least one subset of the observed variables forms such a sufficient set witte2018covariate. The object in causal variable selection then becomes how to separate variables which are necessary for identification of the causal effect from those variables which are extraneous witte2018covariate,hernan2002 in the interest of reducing estimate variance or covariate dimensionality greenland2015statistical,witte2018covariate.\\ In non-experimental observational studies, we do not have full access to a sufficient set in many realistic settings, and important confounding pathways remain unblocked VanderWeelePeng_Sensitivity2017,hill_sensitivity. This is referred to as unmeasured confounding or endogeneity in the statistics and econometrics literatures respectively. However, applied researchers currently rely on variable selection techniques such as lasso, step-wise, change-in-estimator selection, and outcome and/or treatment oriented approaches Talbot2019 despite violating their underlying assumptions.\\ The use of variable selection techniques is to avoid conditioning on negligible confounding pathways without introducing meaningful bias to the estimator. In this paper we explore how this intuition can break down under even mild violations of the underlying assumptions, particularly under the threat of bias amplification.\\ First consider data generated from the following directed acyclic graph (DAG) (Figure (ref)) and set of structural equations:
where $\boldsymbol{Y}$ is the outcome, $\boldsymbol{A}$ is the treatment of interest, $\boldsymbol{U}$ is an unmeasured variable and $\boldsymbol{BAV}$ refers to 10 different potential bias amplifying variables that are measured and affect both $\boldsymbol{A}$ and $\boldsymbol{U}$ but have no direct effect on $\boldsymbol{Y}$. This model contains one confounding path that cannot be blocked ($\boldsymbol{A}\leftarrow \boldsymbol{U} \rightarrow \boldsymbol{Y}$) and 10 confounding paths ($\boldsymbol{A}\leftarrow \boldsymbol{BAV_1} \rightarrow \boldsymbol{U} \rightarrow \boldsymbol{Y}$; \dots ; $\boldsymbol{A}\leftarrow \boldsymbol{BAV_{10}} \rightarrow \boldsymbol{U} \rightarrow \boldsymbol{Y}$) that can be blocked by including the measured $\boldsymbol{BAV_i}$ variables in the model. However, including any of these $\boldsymbol{BAV_i}$ might also increase bias (potential bias amplifying variables). Our goal is to find the least biased estimator of the average causal effect of treatment ($\beta_a$).\\ By including more $\boldsymbol{BAV}s$, intuition suggests the remaining unmeasured confounding bias should decrease. However as demonstrated in the bias amplification literature pearl2012class,pearl2011invited,middleton2016bias conditioning on confounders may still increase bias. For example, suppose further the 10 observable variables account for 90% of the variance in the variable $\boldsymbol{U}$ responsible for unmeasured confounding. The blue violin plot in Figure (ref) represents the density of the estimates from the true outcome model with the treatment and both measured/unmeasured confounding variables included as regressors. As expected the estimates are approximately normally distributed around the true value $\beta_a =0.7$. The green violin plot represents the biased estimates from the naive model, the simple regression of the outcome, $\boldsymbol{Y}$, on the treatment $\boldsymbol{A}$, which does not include any of the confounders (measured or unmeasured). The red violin plot represents the linear model adjusted for all 10 measured confounders which account for $90\%$ of the unmeasured confounding. The adjusted model performs much worse than the naive model both in terms of bias (0.73 compared to 0.43, interpretable as standard deviations) and variance (standard deviation of 0.1 compared to 0.02). In fact, in 4990 of 5000 simulations the adjusted estimate was farther from the truth than the naive estimate and nearly 65% of the adjusted estimates had the incorrect effect sign.
The purpose of this paper is to explain why model selection intuition fails us in this case and how we can use a combination of data and simulation approaches to improve model selection. We build upon an emerging theoretical literature exploring a class of variables which can amplify existing unmeasured confounding bias WOOLDRIDGE2016232,pearl2012class,ding2017instrumental. This class of variables is potentially very large and common in practical applications. Finally, we discuss possible model selection strategies to minimise bias and variance when unmeasured confounding is believed to be present.\\ We adopt a matrix notation framework to characterize this problem because 1) we can easily generalize to a much larger class of directed acyclic graphs and structural equations than previously studied, 2) it offers a unifying geometric explanation in the context of least squares estimation and 3) it offers a solid foundation for how to build data informed model selection procedures. Finally, we develop a procedure for simulating from a more complete parameter space in a way that respects the underlying amplification process. In addition to lending itself better to articulating and answering causal simulation questions this procedure helps explain why some previous studies have incorrectly concluded that applied investigators need not worry about amplification in practicemeyers2011. We evaluate the challenges of implementing this approach with a real clinical example with binary treatment.
Figure (ref) shows a directed acyclic graph (DAG) for a simpler model containing both measured and unmeasured confounding. Let $\boldsymbol{Y}$ represent the outcome and $\boldsymbol{A}$ is the treatment or variable of interest. Let $\boldsymbol{U}$ be an unmeasured confounding variable that we cannot include in a regression model, but which has a functional relationship with both $\boldsymbol{Y}$ and $\boldsymbol{A}$. The bias amplifying variable ($\boldsymbol{BAV}$) in this DAG is analagous to $\boldsymbol{U}$ in that it is a cause of $\boldsymbol{Y}$ and a cause of $\boldsymbol{A}$, however we are able to measure it. It could naturally be included in any reasonable regression modeling scheme. Intuition from causal variable selection techniques would tell us to include $\boldsymbol{BAV}$ in the regression to reduce bias because it forms a confounding path ($\boldsymbol{A} \leftarrow \boldsymbol{BAV} \rightarrow \boldsymbol{Y}$). However, as has been demonstrated pearl2012class,pearl2011invited,ding2017instrumental,middleton2016bias blocking this confounding path can actually increase or amplify the bias relative to the naive estimate only including $\boldsymbol{A}$.\\ Here we will consider a special case with a linear system of equations. The target estimand is the average causal effect (ACE), which is simply $\beta_a$ in the linear model case (See appendix section (ref)). $\boldsymbol{U}$ is unmeasured and thus we cannot identify the ACE from the observed data, but we are interested in estimating the quantity with as little bias as possible.\\ The true model representing Figure (ref) under the linear association assumption are:
where $\alpha_y$ and $\alpha_a$ are the intercept terms for Y and A respectively. We use the form $\beta_x$ throughout this paper to denote true linear regression coefficients for some variable $X$ on the outcome Y. For example, the true regression coefficient for U on Y is $\beta_u$. Analogously, the true regression parameter for some variable $\boldsymbol{X}$ on the treatment A is represented by $\gamma_x$. The estimates of these parameters by OLS are denoted by $\hat{\beta_x}$, $\hat{\gamma_x}$ with additional superscripts to clarify which set of estimating equations the estimator is derived from. By assumption $\boldsymbol{\epsilon_1}$ and $\boldsymbol{\epsilon_2}$ are error terms independent of each other and all other variables represented in the DAG. We assume that $\boldsymbol{\epsilon_1}$ and $\boldsymbol{\epsilon_2}$ have mean 0, and some variance $\sigma^2_{\epsilon_{1,2}}$. In simulation experiments, we additionally assume that the error terms are normally distributed, but this is more than what is necessary for the theoretical results to hold.
To tackle the question of model selection we must derive properties of the feasible $\hat{\beta_a}$ estimators. To this aim we propose expressing OLS estimates using matrix notation and ideas borrowed from the partial regression literature. Further, we propose considering also the probability limits of the estimators to extend our results to more general and realistic cases of bias amplification (See appendix (ref)). For the naive estimator we are estimating the following simple regression:
Notice that $\boldsymbol{\upsilon_1}$ represents the error term in the estimating equation as opposed to $\boldsymbol{\epsilon_1}$ in the true underlying model. We can write $\boldsymbol{\upsilon_1}$ as , $\boldsymbol{\upsilon_1} = \boldsymbol{U}\beta_u +\boldsymbol{BAV}\boldsymbol{\beta_{bav}} + \boldsymbol{\epsilon_1}$. Unbiased estimation of the the naive model by OLS requires the assumption that $E[\boldsymbol{\upsilon_1}|\boldsymbol{A}] = 0$, but this of course is not true. The bias is a result of this erroneous assumption. The naive estimator bias is a special case of the classic omitted variables problem, where we have two omitted variables which are related to both the treatment and the exposure, $\boldsymbol{U}$ and $\boldsymbol{BAV}$.\\ Let $\hat{\beta_a}^{naive}$ be the estimate of $\beta_a$ from the naive model (ref). Throughout this paper we will consider the matrix $\boldsymbol{Z}$ to be a matrix of all the variables that we include in a regression that are not the variable of interest $\boldsymbol{A}$, in other words control variables in a selection on observables approach. In the naive model, $\boldsymbol{Z} = \boldsymbol{1}$ where throughout $\boldsymbol{1}$ will denote an $n\times 1$ vector of 1s. In matrix notation, applying the Frisch-Waugh-Lovell (FWL) theorem (see appendix (ref)), we can write $\hat{\beta_a}^{naive}$ as:
where $\boldsymbol{M_{\boldsymbol{1}}}$ is a centering projection matrix, defined and described in detail in appendix section (ref). In the case of linear relationships between all the variables, following pearl2012class, this estimator has the following expectation:
The absolute bias for the ACE then clearly is $|\beta_u\frac{\gamma_u\sigma_u^2}{\sigma_a^2} + \beta_{bav}\frac{\gamma_{bav}\sigma_{bav}^2}{\sigma_a^2}|$. Now consider the estimates resulting from further conditioning on the observable $\boldsymbol{BAV}$ variable, i.e fitting the following model:
We will denote the resulting estimator $\hat{\beta_a^{|bav}}$ which can be written as follows by again applying the FWL theorem:
where $\boldsymbol{Z} = \boldsymbol{[\operatorname*{\boldsymbol{1}},BAV]}$ and $\boldsymbol{M_z}$ is the annihilator projection matrix of the matrix Z (see appendix (ref) for details and properties). Again following Pearl pearl2012class, the expectation of $\hat{\beta_a^{|bav}}$ is:
In the appendix (see (ref)) we explicitly show Pearl's derivation and how it relies on the conditional expectation $E[\boldsymbol{U}|\boldsymbol{A}, \boldsymbol{BAV}]$ being linear in both $\boldsymbol{A}$ and $\boldsymbol{BAV}$. Pearl's derivation is limited in that it is cumbersome and does not generalize well to a broad class of DAGs and functional forms. A simple example where we are unable to use Pearl's method is the case of an interaction term in the exposure structural equation between $\boldsymbol{U}$ and $\boldsymbol{BAV}$. Suppose we replace equation (ref) with:
The above equations show that $E[\boldsymbol{U}|\boldsymbol{A},\boldsymbol{BAV}]$ is nonlinear in $\boldsymbol{A}$ in $\boldsymbol{BAV}$, and cannot be represented by an unbiased least squares projection of the form $\boldsymbol{U} = \alpha_u + \boldsymbol{A}\zeta_a + \boldsymbol{BAV}\zeta_{bav} + \boldsymbol{\epsilon_3}$ as required by Pearl's derivation method (see Appendix (ref) for details), where $\zeta_i$ represents the true regression coefficient for variable $i$. If we impose further strict distributional assumptions over all the variables, we may still be able to directly solve the conditional expectation and find an expression for bias in terms of the underlying parameters. In many applied cases, these distributional assumptions will not be justified, particularly assuming a distribution for the unmeasured confounding which will always be untestable.\\ In contrast, if we consider the probability limits, we do not need to assume that $E[\boldsymbol{U}|\boldsymbol{A}, \boldsymbol{BAV}]$ is linear, nor do we have to make any additional distributional assumptions to find meaningful limiting expressions for our estimators in a broad class of clinically relevant circumstances. In addition to giving rise to a meaningful interpretation, the closed form asymptotics we derive allow us to more easily harness domain knowledge about the underlying causal process for the purpose of model selection.\\ Since we are still interested in the finite sample expectation of the estimators and the bias directly, we report the expectations when appropriate and feasible. The probability limit facilitates insight under weaker assumptions than those necessary to derive exact forms of the expectations. Additionally, in some cases, like the linear model of Pearl pearl2012class, the probability limits for $\hat{\beta_a}^{naive}$ and $\hat{\beta_a^{|bav}}$ are precisely equal to their expectations (see appendix (ref)).\\
Pearl pearl2012class presented bias amplification results under the assumption of standard normal variables. Here we do not make any assumptions about distribution, mean, or variances for two reasons. First, these assumptions are not strictly necessary to the result. Second, avoiding these assumptions helps clarify some of the mechanics and the intuition behind the phenomenon of bias amplification. As the amplifying term in the denominator of equation (ref), $\sigma_a^2 - \gamma_{bav}^2\sigma_{bav}^2$, gets smaller the bias due to the unmeasured confounding path ($\boldsymbol{A} \leftarrow \boldsymbol{BAV} \rightarrow \boldsymbol{Y}$), $\beta_u\times\gamma_u\sigma_u^2$, increases. This is because when we specify the functional form of a system of random variables and conditional independence assumptions, we are also determining a formula for its variance. Under 1) the structural equation we specified for the exposure (equation (ref)), and 2) the independence assumption between the unmeasured confounding and the bias amplifying variables, the variance is equivalent to:
Rearranging equation (15) to $\sigma_a^2 - \gamma_{bav}^2\sigma_{bav} = \gamma_u^2\sigma_u^2 + \sigma_{\epsilon_2}^2$, it becomes clear the residual variance in $\boldsymbol{A}$ (i.e. not due to the $\boldsymbol{BAV}$ variable) is equal to the sum of the variance due to $\boldsymbol{U}$ and the independent variation $\sigma_{\epsilon_2}$. Therefore, the amplification of the bias in the general case depends not only on the magnitude of $\gamma_{bav}$ (i.e the strength of association between $\boldsymbol{BAV}$ and $\boldsymbol{A}$), but how much of the treatment variance the $\boldsymbol{BAV}$ variable linearly accounts for. When we assume that all the variables are standard normal, the amplification becomes $1 -\gamma_{bav}^2$ as presented in Pearlpearl2012class,because the variance of standard normal variables is equal to 1 ($\sigma_a^2,\sigma_{bav}^2 = 1$).\\ In order to visualize this phenomenon, we use ideas from partial regression plots velleman1981efficient. By the FWL theorem, we can always pre-multiply an estimating equation by the residual-making variables of a set of regressors and get the same estimates (see appendix (ref) for further details). For example, the following two regression equations produce the same numerical estimates of $\hat{\beta_a}^{naive}$:
Equation (ref) is the model for a simple linear regression of a modified outcome, $\boldsymbol{M_{\boldsymbol{1}}Y}$ on a modified treatment, $\boldsymbol{M_{\boldsymbol{1}}A}$ (See appendix (ref)). There is no intercept term as the mean of the modified treatment must be equal to zero.
Similarly equations (ref) and (ref) produce equivalent estimates of $\hat{\beta}_a$, where $\boldsymbol{Z} = [\boldsymbol{1} \quad \boldsymbol{BAV}]$ is a column of 1s and the $\boldsymbol{BAV}$ variable. Equation (ref) is a single variable regression on a transformed set of variables. The modified $\boldsymbol{Y}$ is produced by taking the residuals from regressing $\boldsymbol{Y}$ on a column of 1s and $\boldsymbol{BAV}$, in other words the dependent variable is the remaining variation in $\boldsymbol{Y}$ which is not linearly associated with an intercept and $\boldsymbol{BAV}$. The independent variable is the remaining variation in $\boldsymbol{A}$ not linearly associated with a column of 1s and $\boldsymbol{BAV}$. Since we have reduced the multi-variable regression to a simple linear regression we can easily visualize the amplification process via a partial regression plot.
In Figure (ref) the left visualizes the naive regression equation (ref), whereas the blue graph on the right visualizes the regression equation (ref) that includes BAV. The data was simulated from a special case of equations (ref) and (ref), with $n= 1000$. Details can be found in the appendix (Section (ref)).\\ The unbiased ACE is the slope of the black line ($\beta_a = 0.2$) in these plots. The slope of the blue line (equal to the OLS estimator from the amplifying model) is clearly farther away from the true slope (in black) compared to the slope of the red line from the naive model, and thus the conditional estimator is more biased.\\ Note first that including $\boldsymbol{BAV}$ in the model reduces the variance in the adjusted treatment, which can be seen by comparing the relative sparsity of points along the x-axis in red compared to the relative density of points along the x-axis in blue. However, if we inspect the spread of points vertically along the y-axis, we can see that the red and blue samples are similarly dispersed in this dimension because conditional on the treatment, linear combinations of $\boldsymbol{BAV}$ explain very little of the variance in the outcome. Most importantly, including $\boldsymbol{BAV}$ does not change the variance in $\boldsymbol{Y}$ due to $\boldsymbol{U}$, the unmeasured confounder. As a result, the line of best fit of the adjusted model must be steeper in absolute terms in order to maintain the association between the adjusted response and adjusted treatment over the narrower variation of the adjusted treatment variable. When we add the $\boldsymbol{BAV}$ to the regression model, the bias is 0.14 larger in absolute terms (or approximately 65% greater in relative terms) than the naive estimate, even though it blocks a confounding path between the treatment $\boldsymbol{A}$ and the outcome $\boldsymbol{Y}$. More simply, trying to block a confounding path with weak response association can amplify bias in causal effect estimation because it increases the proportion of treatment association due to unmeasured confounding on unblocked paths.\\ The magnitude of bias amplification can be potentially very large. The absolute bias of the $\boldsymbol{BAV}$ estimator will be larger than the absolute bias of the naive estimator whenever $|\beta_u\frac{\gamma_u\sigma_u^2}{\sigma_a^2} + \beta_{bav}\frac{\gamma_{bav}\sigma_{bav}^2}{\sigma_a^2}| < |\beta_u\frac{\gamma_u\sigma_u^2}{\sigma_a^2 - \gamma_{bav}^2\sigma_{bav}^{bav}}|$ if the relationships are linear. In particular, the bias is greater if $\boldsymbol{BAV}$ is not strongly associated with the outcome (i.e small values of $\beta_{bav}\sigma_{bav}^2$)). The special case $\beta_{bav} = 0$ implies that $\boldsymbol{BAV}$ is a true instrumental variable. Instrumental variables were in fact the leading case for the discovery of this class of bias amplifiersWOOLDRIDGE2016232,pearl2012class. If there are no interaction terms (i.e a model that is linear in the original variables) adding an instrumental variable always weakly increases absolute bias in OLS relative to the naive model, with equality only when there is no unmeasured confounding WOOLDRIDGE2016232,pearl2012class.\\ In summary, variable selection approaches which aggressively target confounding paths with strong associations with treatment and weak associations with outcome are at grave risk of bias amplification as they are much more sensitive to the assumption that a full sufficient set is measurable. Adding controlling variables in proportion to their ability to predict the treatment in linear models only becomes a bias reducing approach if the resulting variable set satisfies ignorability assumptions. This is often not possible or extremely unlikely in many non-experimental settings.
The danger of including variables that are strongly associated to the treatment is that we cannot identify unmeasured confounding. Consider the probability limit of the estimator in equation ((ref)):
where $\sum_{i=1}^{n}\hat{\xi_i}^2$ is the estimated sum of squared residuals from the regression of the treatment on $\boldsymbol{BAV}$ and an intercept term (note that $\boldsymbol{\xi} = \boldsymbol{U}\gamma_u + \boldsymbol{\epsilon_2}$ from equation (ref)). Since all the variables necessary to estimate $\hat{\xi}$ are observable, we can identify the denominator, or the amplifying term. Equation (ref) shows that the amplifying term is numerically equivalent to the sum of squared residuals of $\boldsymbol{A}$ on an intercept column and $\boldsymbol{BAV}$. Note that this result does not even require a limit or expectation to hold. However, it is typically more useful to think about the probability limit in an applied application since the numerator simplifies to a covariance term in the case that the added variables are independent of the unmeasured confounding.\\ In the case that individual treatment assignments are independent and have variance 1, the probability limit of the average of the squared residuals is equal to one minus the proportion of variance of $\boldsymbol{A}$ explained by $\boldsymbol{Z}$ middleton2016bias, i.e.
Thus far we have only considered the DAG in Figure (ref) under the restrictive assumption of linear associations amongst variables. The identifiability of the bias amplification term of the preceding section can be extended in two important ways. First, the results extend to the addition of any $p$ bias-amplifying variables by simply increasing the number of columns of $\boldsymbol{Z}$ to include any number of bias-amplifying variables, as the FWL theorem allows for arbitrary numbers of columns as long as $\boldsymbol{Z}$ is of full rank. We include examples of multiple bias-amplifying variables in section (ref). Second, in the following subsection, we provide the details of how to extend the result to non-linear associations.
In the previous sections we assume that the data generating processes governing the treatment and the outcome are linear. However, in order to identify the amplification term with observable variables this is not strictly necessary. First, we relax the assumption that the data generating process of $\boldsymbol{A}$ is linear, allowing it to be any arbitrary function $\boldsymbol{A} = f(\boldsymbol{U},\boldsymbol{BAV}, \boldsymbol{\epsilon_2})$ but let the model for $\boldsymbol{Y}$ remain unchanged such that equation (ref) holds. In the appendix (section (ref)) we show that the numerical form of the conditional estimator to be $\operatorname*{\hat{\beta}_a^{|bav}} = \beta_a + \beta_u\boldsymbol{\frac{A^TM_zU}{A^TM_z}} + \boldsymbol{\frac{A^TM_z\epsilon_2}{A^TM_zA}}$, where $\boldsymbol{M_zA}$ is by definition the vector of residuals from the regression of $\boldsymbol{A}$ on the columns of $\boldsymbol{Z}$, which we specified to mean $\boldsymbol{BAV}$ and a column of 1s. As in the fully linear case, $\boldsymbol{A^TM_zA}$ is the sum of squared residuals. The residuals will be the treatment $\boldsymbol{A}$, removed of linear components of $\boldsymbol{Z}$ which do not directly depend on the underlying function governing the relationship between $\boldsymbol{A}$ and $\boldsymbol{BAV}$. Note that $\boldsymbol{A^TM_zA} \leq \boldsymbol{A^TM_{\boldsymbol{1}}A}$ since $\boldsymbol{Z}$ contains the column of 1's, and so amplification will occur as long as the treatment is some function of $\boldsymbol{BAV}$ producing a positive correlation between the treatment and bias amplifying variable. As shown in section (ref), the extent of the amplification will be determined by the linear correlation between $\boldsymbol{A}$ and $f(\boldsymbol{U},\boldsymbol{BAV},\boldsymbol{\epsilon_2})$.\\ When we allow for non-linear associations in the outcome, an important point to clarify is that adjusting for $\boldsymbol{BAV}$ in OLS will not necessarily be sufficient to block the confounding path that $\boldsymbol{BAV}$ forms. Consider a very simple extension to the outcome model as follows:
where for simplicity $f(\boldsymbol{BAV})$ is a high-ordered polynomial term that is correlated with $\boldsymbol{BAV^2}$ after adjusting for the linear term. If we adjust for $\boldsymbol{BAV}$ and not the squared term, the causal estimates will clearly suffer from omitted variable bias since the squared term remains correlated with both the treatment and outcome. However, we can still identify the amplification. The proper ACE under equation ((ref)) are $\frac{\partial E[\boldsymbol{Y}|\boldsymbol{A},\boldsymbol{U},\boldsymbol{BAV}]}{\partial \boldsymbol{A}} = f_1^{\prime}(\boldsymbol{A})\beta_a$. Bias needs to be evaluated as deviations from the true causal effects ($f_1^{\prime}(\boldsymbol{A})\beta_a$) and not the parameter $\beta_a$.
where $\hat{\xi_i}^2$ are the squared residuals from the regression of $\boldsymbol{A}$ on $\boldsymbol{BAV}$ and a constant. Although the amplified bias is different from the simple case, the factor by which the bias is amplified is still identifiable using a regression depending only on observables. The amplification results from misspecification of the relationship between $\boldsymbol{BAV}$ and the outcome. Similar to before there will be some cases where the amplification of the U-bias is outweighed by the reduction in omitted variable bias due to $\boldsymbol{BAV}$ and other cases in which this would not be the case.\\ Now consider a fully non-linear, but still additive, model specification:
If we estimate the linear naive and linear adjusted models as before we get the following estimates:
Looking at equation (ref), the unmeasured confounding pathway remains amplified and we can estimate the residuals which cause the amplification from observable quantities. However, the direction of the shift in overall bias is unclear when we condition on $\boldsymbol{BAV}$. Note that the middle term is unambiguaously larger for the $\boldsymbol{BAV}$ model than the naive model. The third term in the $\boldsymbol{BAV}$ model will be smaller than the naive model in the numerator, but larger in the denominator. Most troubling, by allowing $\boldsymbol{Y}$ to be a nonlinear function of $\boldsymbol{A}$, we can no longer predict whether the first term is getting closer or farther from the truth. However, if $f_1(\boldsymbol{A})$ and $f_3(\boldsymbol{BAV})$ are known or can be well approximated, we can estimate $\frac{\operatorname*{plim\,} \frac{1}{n}\boldsymbol{A^TM_z} f_1(\boldsymbol{A})}{\operatorname*{plim\,} \frac{1}{n} \sum_{i=1}^n \hat{\xi_i}^2}$ and $\frac{\operatorname*{plim\,} \frac{1}{n}\boldsymbol{A^TM_z} f_3(\boldsymbol{BAV})}{\operatorname*{plim\,} \frac{1}{n} \sum_{i=1}^n \hat{\xi_i}^2}$, since they are functions entirely of observables. To estimate the numerators, three regressions should be run: the exposure on $\boldsymbol{Z}$, $f_1(\boldsymbol{A})$ on $\boldsymbol{Z}$, and $f_3(\boldsymbol{BAV})$ on $\boldsymbol{Z}$. By storing the residuals and combining them appropriately, the first and third term in equation (ref) can be estimated up to $\beta_a$ and $\beta_{bav}$ (For more details see appendix (ref)).\\ If we make no assumptions about the functional form, and allow for non-linearities, and interactions between all variables, we can show that the OLS adjusted estimator is always the expression below (see appendix (ref)):
When $\boldsymbol{Z}$ includes an intercept column, both $\boldsymbol{(M_zA)^T}$ and $\boldsymbol{M_zY}$ will have mean zero and thus we can think of the numerator as an empirical estimate of the covariance between the residuals from the regression of $\boldsymbol{A}$ on $\boldsymbol{Z}$ and the residuals from the regression of $\boldsymbol{Y}$ on $\boldsymbol{Z}$. Unmeasured confounding bias in OLS occurs when after projecting out linear combinations of the controlling variables, $\boldsymbol{Z}$, there remain linear associations between the outcome and the treatment due to unobserved variables. The part of the bias due to unmeasured confounding is amplified whenever the control variables explain variance in the treatment. Holding all else constant, as the residuals from the regression of the treatment on $\boldsymbol{Z}$ decrease in magnitude, the absolute value of the estimator $\hat{\beta_a^{|z}}$ will increase in magnitude. This is a general form of the result we showed in the previous section which is extremely powerful in that it encaptures a very large class of structural equations and DAGs. However, the cost of this generality is that without making more specific assumptions about the particular form of the model, and in particular the outcome model, it becomes more difficult to incorporate the knowledge of the amplification factor into our model selection and thus apriori know which of the two estimators, $\operatorname*{\hat{\beta}_a^{naive}}$ or $\operatorname*{\hat{\beta}_a^{|bav}}$, will be less biased. Interaction terms, for example, are an additional difficulty. Pearl pearl2012class, for example, showed that under a simple interaction effect between the unmeasured confounding and some function of a pure instrument, the adjusted estimator can be less biased than the naive case. To properly evaluate estimators in the context of bias amplification requires appropriate simulations. In the next section, we describe how to avoid the pitfalls of previous simulation work meyers2011.\\
In our experience, simulating bias amplification is challenging in a number of subtle, but important ways. Our context of interest is assessing the potential for bias amplification in an analysis of an observational study in which we have measured several independent variables and the outcome but there might be an unmeasured confounder. We are interested in evaluating the feasible estimators we have developed in the previous sections, $\operatorname*{\hat{\beta}_a^{naive}}$ and $\operatorname*{\hat{\beta}_a^{|bav}}$ for example, with respect to possible data sets generated by a class of DAGs and structural equations. In this section we show that if we constrain certain aspects of the simulated data (in particular, the marginal variances of observed quantities), we are better able to articulate and answer causal questions about the effect of bias amplification on proposed estimators. While we discuss the example of bias amplification simulations specifically, this section has implications for simulating data to test causal estimators more broadly.
Now, consider the challenge of determining the effect of increasing unmeasured confounding on bias amplification in Figure (ref). We might, for example, be interested in how large an unmeasured confounder must be, with fixed amplifying variables, to cross some threshold of bias in the adjusted model as part of a sensitivity analysis. To answer such a question we must define clearly what is meant by the strength of an unmeasured confounder. In Figure (ref), there are two edges which determine the overall bias due to the unmeasured confounding path through U: the edge from $\boldsymbol{U}$ to $\boldsymbol{A}$ and the edge from $\boldsymbol{U}$ to $\boldsymbol{Y}$. The bias due to the unmeasured confounding path through U in the naive model is simply the product of the weight of these two edges, scaled by the variance of the treatment as shown in equation (ref). The extent to which bias can become amplified, however, is not symmetric with respect to the weight of the edges $\boldsymbol{U} \rightarrow \boldsymbol{A}$ and $\boldsymbol{U} \rightarrow \boldsymbol{Y}$, since amplification is the result of variance explained in the treatment as discussed in section (ref). There is more potential for amplification of a strong unmeasured confounder (in the sense the product of the confounding edges is large) when the strength is due to $\boldsymbol{U}$ being a strong cause of $\boldsymbol{Y}$ compared to a strong cause of $\boldsymbol{A}$. This is because when $\boldsymbol{U}$ is a strong cause of $\boldsymbol{A}$, the $\boldsymbol{BAV}$ can only explain a small amount of the variance of $\boldsymbol{A}$, limiting the possible amount of bias amplification. Thus to answer a causal question about the effect of increased unmeasured confounding on bias amplification we should only vary one of the confounding edges, holding all other edges fixed.\\ As an example, suppose we are interested in the change in bias amplification when we increase the strength of the edge from $\boldsymbol{U}$ to $\boldsymbol{A}$, holding all else constant. This notion of intervening on a single edge of our DAG while holding the others fixed should be familiar to causal inference practitioners since it is the principle behind counterfactual analysis more broadly. Here we want to ensure that our results from varying a single edge are not confounded by variations in other edges as the result of unintended consequences or induced associations.\\ Because the goal is to increase the strength of a single edge, holding all else constant, we must specify a metric by which we measure the strength of the edge. In a fully linear system, we might consider the strength of the edge as the regression coefficient itself, $\gamma_u$, or the proportion of variance explained by $\boldsymbol{U}$, $\frac{\gamma_u^2\sigma_u^2}{\sigma_a^2}$ and the sign of $\gamma_u$. It is tempting to see the two measures as equivalent with different scalings, but this is only true in the context of simulating a single equation. In the context of a system of linear equations, especially with the potential for bias amplification, we argue the relevant quantity is the proportion of variance explained by each child node of the parent variable. This can be seen most easily by examining the bias formula in equation (ref), where the amplifying term is the remaining variation in $\boldsymbol{A}$ unexplained by the potential bias amplifying variables.\\ Consider the implications of treating the coefficients themselves as the relevant measure of edge strength in a simulation trying to determine the effect of increasing the causal association along the path from $\boldsymbol{U}$ to $\boldsymbol{A}$. If we want to increase $\gamma_u$ to $\gamma_u^\prime > \gamma_u$ without changing any other parameters, we must also increase the total variance in the treatment, $\boldsymbol{A}$, since $\sigma_a^2 = \gamma_u^2\sigma_u^2 + \gamma_{bav}^2\sigma_{bav}^2 + \sigma_{\epsilon_2}^2$. A treatment with a larger variance is in some sense a different intervention, and thus this simulation is not compatible with the class of experiments which generated the original data with parameter $\gamma_u$. Further, from the previous sections we know this implies the total amount of variance explained from the bias amplifier $\boldsymbol{BAV}$ is reduced, since $\frac{\gamma_{bav}^2\sigma^2_{bav}}{(\sigma_a^2)^\prime}$ has been reduced. Although we have not changed the parameter $\gamma_{bav}$ we have decreased the extent to which $\boldsymbol{BAV}$ amplifies the bias as seen by examining equation (ref). The increased variance in $\boldsymbol{A}$ in turn modifies the total variance of $\boldsymbol{Y}$. Therefore, the relative proportion of variance of $\boldsymbol{Y}$ that is explained by $\boldsymbol{BAV}$ is modified by changing the causal effect of $\boldsymbol{U}\rightarrow \boldsymbol{A}$, as are the measured proportion of variance of $\boldsymbol{BAV}\rightarrow \boldsymbol{A}$, $\boldsymbol{A} \rightarrow \boldsymbol{Y}$, $\boldsymbol{U}\rightarrow \boldsymbol{Y}$, and $\boldsymbol{BAV}\rightarrow \boldsymbol{Y}$ and their associated covariance terms.\\ We can see in Figure (ref) that by modifying a single coefficient and leaving all other coefficients unchanged we have inadvertently modified the relative proportion of variance explained by the 4 other edges ($\boldsymbol{BAV} \to \boldsymbol{A}$, $\boldsymbol{BAV} \to \boldsymbol{Y}$, $\boldsymbol{U} \to \boldsymbol{Y}$, and $\boldsymbol{A} \to \boldsymbol{Y}$) represented by the wavy arrows. Data generated by the second set of structural equations are not compatible with the constraints of the experiment which generated the first data and by intervening on a single edge we have modified all of the competing effects of interest. Comparing the distribution of estimates produced under $\gamma_u$ and $\gamma^{\prime}_u$ gives us a confounded and thus biased estimate of the impact of increasing the unmeasured confounding through its causal pathway to the treatment on the estimators or functions thereof. We will show that this bias can result in under-estimating the impact of bias amplifying variables.\\ In general, when we vary one of the regression coefficients along a causal pathway, this has upstream and downstream effects on the proportion of variance explained by all variables going into or out of the varied node. In order to keep the proportional effects of the other edges constant, we need to use the error terms of the structural equations ($\boldsymbol{\epsilon_1}$ and $\boldsymbol{\epsilon_2}$) to absorb the shocks to the marginal variances.\\ In Figure (ref), if we change $\gamma_u$ and simply adjust the structural error term $\epsilon_a$ such that the total variance in $\boldsymbol{A}$ remains constant, we can isolate the effect of modifying $\boldsymbol{U}\to \boldsymbol{A}$. Below in Figure (ref) we visualize the consequences of failing to hold the variance of the treatment when we modify $\gamma_u$.
In red, for Figure (ref), we simulate bias amplification where $\gamma_u = 0.3$. In green, we simulate bias amplfication where $\gamma_u$ is increased to $0.55$ holding all other parameters constant, thus allowing the total variance of the treatment to grow from $1$ to $1.21$. This has the downstream effect of also increasing the variance of the outcome from $1$ to $1.02$. This also then impacts the relative proportions of variance explained of the treatment and the outcome that are explained by $\boldsymbol{U}$ and $\boldsymbol{BAV}$ respectively. Notice that the bias increases from $20\%$ ($\frac{0.36-0.3}{0.3})$ to $43\%$ $(\frac{0.43-0.3}{0.3})$. In blue, we increase $\gamma_u$ from $0.3$ to $0.55$, but re-normalize the variance in the treatment to remain constant at $1$. The bias now increases further to $67\%$ ($\frac{0.5-0.3}{0.5}$) with respect to the original simulation in red. We do this by decreasing the variance of the independent noise term, $\boldsymbol{\epsilon_a}$ to $\boldsymbol{\epsilon^\prime_a}$, allowing it to absorb the increase in variation from $\boldsymbol{U}$. When we do not fix the variance, we underestimate the impact of the amplifier on both the bias and the variance {\em because the unfixed variance case simulates a different kind of intervention due to the change in variance of the treatment variable}. In the simulation above, by not keeping the variance fixed in $\boldsymbol{A}$ we implicitly reduced the amount of variance that $\boldsymbol{BAV}$ accounts for in the treatment from $36\%$ to $30\%$. In effect, we were comparing the distribution of
when a more fair causal counterfactual would be to compare the distribution of
Therefore, our simulation experiment results in green are distorted because when we increased the unmeasured confounding through $\boldsymbol{U} \to \boldsymbol{A}$, we also decreased the strength of the bias amplifying variable through the pathway $\boldsymbol{BAV} \to \boldsymbol{A}$. Notice that this bias will impact decisions and conclusions we might make about the merits of different estimators in this context. For example, below we compare the conditional estimator, $\operatorname*{\hat{\beta}_a^{|bav}}$, to the naive estimator, $\operatorname*{\hat{\beta}_a^{naive}}$ with respect to their bias in the same three simulation set ups.\\ In Figure (ref) we show the direct comparison of the bias for the conditional and the naive estimators. When we increase the unmeasured confounding through $\gamma_u$ but fail to renormalize the treatment variance, we do not capture the full extent to which the conditional estimator amplifies the bias. If we compared the green and red plot it would seem that nearly doubling the unmeasured confounding coefficient only has a small impact on the relative bias of the naive and conditional estimator, since the relative bias only increased from 0.07 to 0.09 ($29\%$). By comparing the green density plot to the blue, we see the relative bias doubles (from 0.07 to 0.14). Therefore, the decision to use the naive or conditional estimator is in fact much more sensitive to the amount of unmeasured confounding than it would appear under the improper simulation with floating variance. It is extremely important to do these kinds of simulations properly particularly in the context of sensitivity analysis where we are testing the performance of estimators with respect to untestable assumptions such as unmeasured confounding.\\ To properly simulate bias amplification and answer questions of clinical concern with respect to the merits of potential estimators, we must think of the structural equations as an interconnected system. While we typically specify such equations from the perspective of determining their conditional means, the structural equations along with our independence assumptions determine the variances of the variables in the system. Above, this necessitates increasing the strength of the edge $\boldsymbol{U}\rightarrow \boldsymbol{A}$ while holding all other edges constant, which requires us to re-normalize the variances to maintain the strength of the edge $\boldsymbol{BAV} \rightarrow \boldsymbol{A}$.\\ In appendix section (ref), we consider the properties of a simulation experiment aiming to vary the strength of the edge $\boldsymbol{BAV} \to \boldsymbol{A}$. We show that in the case that we fail to fix the variance of the treatment that the bias of the conditional estimator $\operatorname*{\hat{\beta}_a^{|bav}}$ is invariant to $\gamma_{bav}, \forall \gamma_{bav} \in (-\infty, \infty)$, but that the naive estimator is strictly increasing in $\gamma_{bav}$. It is clear from the theory we developed in section (ref) that if we increase the edge from $\boldsymbol{BAV} \to \boldsymbol{A}$ that amplification should strictly increase, but if we allow the variance in $\boldsymbol{A}$ to increase as the parameter increases, the amplification effect is precisely cancelled out.\\ In general terms, simulating linear systems of location-scale family random variables requires first fixing the variances of the variables in the DAG. The relevant quantity determining the strength of the various edges are ratios of variances and covariances of the upstream parent nodes to the variance of the child node in determining the edge's strength. Since the effects are relative, in a simulation context we can normalize the variances to 1 or set them to the expected/observed variances of the data in a particular context. For simplicity we will demonstrate the normalized approach. In Figure (ref), this means that $\sigma_u^2 = \sigma_a^2 = \sigma_{bav}^2 = \sigma_y^2 = 1$.\\ The second step is to be explicit about independence and conditional independence assumptions. Given the independence assumptions, we can specify the covariance matrix of each child variable $\boldsymbol{Y_{child}}$ in terms of the matrix of $k$ arbitrary parent variables which form the edges going into the child variable, $\boldsymbol{{Y_{child}}}$.
The diagonal of all the parent covariance matrices is 1 since we have normalized all variables pictured in the DAG. The covariances themselves will be determined by the independence assumptions, the edges connecting the child nodes, and their structural equations. Essentially we are choosing the proportion of the child variation that the variances and the covariances of the parent variances explain. The error terms, $\epsilon$'s are the only non-normalized variances, and they absorb the shocks when we increase and decrease the strength of the edges of the non-error variables. This maintains the strength of all other relations visualized on the DAG.\\ Since all variance terms must be non-zero (or equivalently that $\boldsymbol{\beta_{parent}^TVar(X_{parent})\beta_{parent}} \leq 1 = \sigma_{child}^2$), the variance equations define bounds on the simulation parameter space. In the above example, conditional on holding the strength of the edges $\boldsymbol{U} \rightarrow \boldsymbol{Y}, \boldsymbol{BAV} \rightarrow \boldsymbol{A}, \boldsymbol{BAV} \rightarrow \boldsymbol{Y}, \boldsymbol{A} \rightarrow \boldsymbol{Y}$, $\gamma_u \in (-0.893,0.893)$ defines the feasible range. That is, the edge $U\rightarrow A$ can explain up to 79.75% of the variation in $\boldsymbol{A}$ ($\frac{\gamma_u^2\sigma_u^2}{\sigma_a^2} = \frac{\gamma_u^2}{1}$) since the edge $\boldsymbol{BAV}\rightarrow \boldsymbol{A}$ explains 20.25% of the variation already. In general, the extent to which an edge can explain variation in the child node is constrained by the other child nodes and the covariance structure between those variables. A parameter however, such as $\gamma_u$ may be constrained by more than one set of inequalities. In this particular case $\gamma_u$ has to satisfy the following inequalities:
where conditional on the strength of the particular edges ($\boldsymbol{U}\rightarrow \boldsymbol{Y}$, $\boldsymbol{BAV} \rightarrow \boldsymbol{A}$, $\boldsymbol{BAV} \rightarrow \boldsymbol{Y}$, $\boldsymbol{A}\rightarrow \boldsymbol{Y}$) in the above simulation, only the first inequality was binding.\\ The nuance here is that the extent to which we can simulate unmeasured confounding depends upon not only how much amplifying we have simulated, but also on the true effect of the treatment on the outcome $\boldsymbol{A} \rightarrow \boldsymbol{Y}$. Since this is an interdependent system of equations, all of the parameters are competing for shares of fixed variances. If the treatment, independent of $\boldsymbol{U}$ and $\boldsymbol{BAV}$, explains the large majority of the outcome variance (i.e, the edge $\boldsymbol{A}\rightarrow \boldsymbol{Y}$), it means the weight of the edge $\boldsymbol{U} \rightarrow \boldsymbol{Y}$ must be relatively small, opposite signed, or the structural equations contain an effect modifier. This in turn constrains $\gamma_u$.\\ Consider again the above simulation experiment where we are interested in varying the strength of $\boldsymbol{U}\rightarrow \boldsymbol{A}$ conditional on all other pathways. Suppose that the pathway $\boldsymbol{A} \rightarrow \boldsymbol{Y}$ explains $64\%$ of the variance in $\boldsymbol{Y}$, i.e that $\beta_a = 0.8$. Now both constraints on $\gamma_u$ are binding and the simulation parameter space is $\gamma_u \in (-0.893, 0.3916)$.\\ In summary, when simulating linear location-scale family systems of equations we start by identifying the DAG and the independence assumptions between variables. Second, our simulation experiment should attempt to answer a causal question about how a proposed estimator behaves in response to an intervention on the weights of causal DAG. Just like experimental design, properly estimating the relevant counterfactural requires that the difference in distributions between our intervention(s) and the control is the effect of the intervention(s) themselves. As demonstrated in this section, simluating linear systems of equations requires varying one of the edges of the DAG holding all else constant, and matching the means and variances of the simulated variables with that of the target observational study we are trying to mimic. This allows us to generate simulations whose distributions are proper counterfactuals. Third, conditional on the other edges, the covariance matrices impose bounds for the parameter space that we can simulate and thus the extent to which we can vary the edge of interest. For a specific realization of the experiment and accompanying valid parameters, the variables are constructed in the downstream direction, that is from parent nodes to child.\\ In the example of simulating the proper intervention in Figure (ref), we first simulate $\boldsymbol{U}$ and $\boldsymbol{BAV}$ independently with variance 1 respectively. Given $\gamma_u$ and $\gamma_{bav}$, the variance of the error term $\boldsymbol{\epsilon_2}$ from equation (ref) is implied and can be simulated. Having $\boldsymbol{U}$, $\boldsymbol{BAV}$ and $\boldsymbol{\epsilon_2}$ allows us to simulate the treatment $\boldsymbol{A}$. Conditional on the already simulated variables, their associated parameters, and $\beta_a$, $\beta_{bav}$, and $\beta_u$, the variance of the error term $\boldsymbol{\epsilon_1}$ is implied and can be simulated. Finally, since all of the child variables for the outcome have been simulated, we can simulate the outcome. To be clear, we can fix proportions of variance explained by each edge in any order we'd like as long as we respect the underlying constraints. However, given an admissible set of weights of the edges we must proceed from parent to child nodes to conduct the simulation.\\ While this method requires us to calculate inequalities and make explicit the implied variance formulas for our variables, the benefits are that we can view our simulation as a well-defined causal experiment matching the constraints of our target study and we get sets of parameter bounds. When we do not keep the variance fixed, there are no defined bounds beyond heuristics, and more importantly, we are no longer matching the data to our target observational study. In many small systems, such as the one in Figure (ref), it is often computationally inexpensive to simulate a discretized approximation to all possible parameter configurations. In extremely large systems we can use domain knowledge to make refinements on these bounds and simulate a reasonable subset of the parameter space. This method allows us to make refinements over edges with strong priors while simulating the entirety of edges with greater uncertainty.
Here we conduct a data simulation for an observational study. We want to consider a medical example with realistic amounts of variance in the treatment and the outcome. Further, we specifically consider the case of a binary treatment which is common in medical applications, biostatistics, and epidemiology. The difficulty, in general, when simulating with real data is that you do not know the true underlying parameter values. In this section, we start with a randomized controlled trial (RCT) and modify it appropriately, so that we can take the intention to treat (ITT) estimate as the true underlying effect for the foundation of our simulations.\\ In our simulation experiment, we keep the treatment data unchanged (thus fixing their variance), and then simulate unmeasured confounding ($\boldsymbol{U}$) and bias amplifiers ($\boldsymbol{BAV}$) in order to modify selected covariates ($\boldsymbol{X}$) and the outcome ($\boldsymbol{Y}$) to produce a synthetic observational experiment. In order to precisely control the relationships between the simulated variables and the real variables we treat the binary treatment, $\boldsymbol{A}$, as though it comes from a latent probit model.
where $\boldsymbol{\tilde{X}} = \frac{\boldsymbol{X}}{\sigma^\prime} + \boldsymbol{BAV}$, and $\sigma^\prime$ is a scaling variable such that $\boldsymbol{X}$ and $\tilde{\boldsymbol{X}}$ have the same population variance. All of the latent variables ($\boldsymbol{U}$, $\boldsymbol{BAV_{n\times k}}$, $\boldsymbol{\epsilon_2}$, and hence $\boldsymbol{A^\star} = \alpha_a + \boldsymbol{U}\gamma_u + \boldsymbol{BAV}\gamma_{bav} + \boldsymbol{\epsilon_2} > \boldsymbol{0}$) are set to come from normal distributions. The details of the how the simulation is performed are in appendix section (ref).\\ For this paper, we use data from R2KJHK_2019, a published RCT with 294 participants and relatively balanced distribution of covariates. While the reseachers examined many outcomes we will focus on the effects of an e-Health intervention in infants on child eating behaviours. The researchers gave the parents in the treatment group access to a "monthly age-appropriate video addressing infant feeding topics together with corresponding cooking films/recipes", and the outcome was eating habits of the child at a later point in time. In the observational study that we want to create, (target observational study) we want to estimate the effect of the treatment on emotional overeating as measured by the Child Eating Behavior Questionnaire (CEBQ).
Our foundation is the unbiased ITT effect from the RCT data regressing the treatment on the outcome ($Y \sim A$) shown in the first column of Table (ref).
In column 1 of table (ref), we see that the ITT estimate is $0.12$. As this is an RCT, we do not expect baseline covariates [Child Food Neophobia Score ($\boldsymbol{CFNS}$), Child Feeding Questionaire ($\boldsymbol{CFQ}$) subscale pressure, and Age of mother ($\boldsymbol{Age_{mother}}$)] to be associated with exposure. We thus assume that the experimental data are generated from the causal DAG in Figure (ref), where $\boldsymbol{X}$ represents the matrix of all three covariates ($\boldsymbol{CFNS}$,$\boldsymbol{CFQ}$, and $\boldsymbol{Age_{mother}}$) after they have been individually standardized to have mean 0 and variance 1. To verify that these variables are not bias amplifiers, that is explain only a negligible proportion of the treatment variance, we also present the results of the regression of the treatment on the 3 covariates in column 3 in table (ref). We can see that jointly and individually the three covariates explain very little of the variance in the treatment, $\mathcal{R^2} = 0.009$. This should be expected in a truly randomized experiment set-up since proper randomization breaks the causal association from the covariates to the treatment.
Since these covariates do not cause $\boldsymbol{A}$ and we have assumed that the ITT estimator is unbiased, when we estimate $\boldsymbol{Y} \sim \boldsymbol{A} + \boldsymbol{CFNS_{score}}+ \boldsymbol{CFQ_{pressure}}+ \boldsymbol{Age_{mother}}$ the expectation and probability limit of $\hat{\beta_a}$ remains unchanged regardless of the strength of association between the covariates and the outcome. However, actual results may vary due to final sample variation. In our RCT data, the unadjusted model estimates a treatment effect of $0.122$ and the adjusted model estimates $0.137$. Since simulation experiments performed in section (ref) all condition on covariates, we consider the covariate adjusted results from the RCT as the gold standard for determining bias due to unmeasured confounding in our simulated data.
Our objective is to simulate data according to the DAG in Figure (ref). To produce the simulations, we took 10000 bootstrap replications of the original outcome, treatment and covariates. From each bootstrap sample of the treatment, $\boldsymbol{A_{bootstrap}}$, of size $n=294$ we simulated the latent variable $\boldsymbol{A^\star}$ using the procedure outlined in the appendix (section (ref)). Next, conditional on the drawn latent samples of $\boldsymbol{A^\star}$ and the bootstrapped covariates, we drew samples for the unmeasured confounding, $\boldsymbol{U}$, and bias amplifying variable, $\boldsymbol{BAV}$. The modified random control variables, $\boldsymbol{\Tilde{X}} = \frac{\boldsymbol{X}}{\sigma^\prime} + \boldsymbol{BAV}$, were produced by adding the bias amplifying variables to a scaled version of the original control variables. Linear combinations of the unmeasured confounding and modified covariates were then added with reasonable values to the outcome such that the following DAG and equations hold (see simulation results).
In section (ref) we showed that the true treatment effect was $0.137$ conditional on the covariates $\boldsymbol{X}$. In the boostrap simulation pictured below, the unbiased model conditional on both the modified covariates, $\boldsymbol{\Tilde{X}}$, and the unmeasured confounding $\boldsymbol{U}$ is $0.136$ as expected. The naive model estimator had an average estimate of $0.234$ in the simulations and thus an absolute estimated bias of $0.097$, or a relative bias of 1.8 standard deviations ($\frac{0.234 - 0.137}{0.053}$) with respect to the unbiased estimate in section (ref).\\ When we further condition on the modified covariates, the absolute bias ($E[|\hat{\beta}_a^{|\tilde{x}}- 0.137|]$) more than doubles to 0.225, and the relative bias increases to 4.3 standard deviations ($\frac{0.36 - 0.137}{0.053}$) with respect to the unbiased estimate in section (ref).
The simulations confirm that bias amplification can be significant even when constrained to problems of realistic variance. Further, we see that bias amplification is potentially a problem for binary outcomes. This underscores the theoretical points made in sections (ref) and (ref) where we showed that the phenomenon behind bias amplification does not require specific distributional assumptions of the variables in the model.\\ More importantly, by combining the methodology outline in the appendix (See (ref)) to simulate measured confounding using real data and the principles for simulating systems of equations in section (ref), we can produce realistic and complete simulations of parameter spaces which match the underlying characteristics of the data. Investigators who choose covariates based on the assumption of no unmeasured confounding can now evaluate the amount of bias amplification that would occur if this assumption does not hold.\\ Finally, in the appendix (section (ref)) we consider an example of a causal simulation experiment with a binary treatment variable under the DAG in Figure (ref) and structural equations (ref), (ref), and (ref). The experiment involves modifying the strength of the edge $\boldsymbol{\tilde{X_1}} \rightarrow \boldsymbol{A}$ and evaluating the impact on the naive and conditional estimators. With binary treatment ($\boldsymbol{A}$), we show that if we fail to hold the variance of the latent treatment ($\boldsymbol{A^\star}$) constant and increase $\gamma_{\tilde{x_1}}$, then it is possible to decrease the amount of observed treatment variance ($\sigma_{a}^2$) explained by $\boldsymbol{\tilde{X_1}}$. Further, the increased treatment variance also decreases the strength of the edge $\boldsymbol{U} \rightarrow \boldsymbol{A}$. As a result of performing the causal simulation experiment improperly, it appears as though that varying the strength of the potential amplifiers has a negligible or negative impact on the resulting bias amplification. The improper and proper approaches to intervention are shown in Figure (ref) and Figure (ref) respectively and the results from these simulations are visualized in the appendix in Figure (ref) in Appendix section (ref). This of course leads to improper inferences regarding the relative merits of the naive and conditional estimators as well. This highlights once again the importance of comparing simulations with comparable properties and ensuring that when we intervene on the edges of our causal diagram that we are not inadvertently varying the edges we mean to keep fixed. Just as in the experimental context, our simulation results become muddled or meaningless if we are not evaluating well-articulated counterfactuals.
Causal model selection techniques have largely been developed under the assumption that a sufficient set of variables is available to create ignorability. When a sufficient set is not available or when a causal variable selection technique does not correctly identify the sufficient set, we are at risk of bias amplification. In the first simulation in section (ref), we showed that even under mild perturbations of the usual assumptions, conditioning on a set of jointly strong proxy variables for $\boldsymbol{A}$ in OLS led to a very biased estimator (0.73 standard deviations on average). Further, most current causal variable selection techniques are likely to include this set of variables since they are significant predictors of the outcome and the treatment as well as variables which cause large changes in estimates when included sequentially.\\ Under threat of bias amplification, treatment-oriented selection techniques for regression analyses using continuous exposure regimes should be used cautiously unless one has strong priors that a sufficient set is available and likely to be identified. We showed in section (ref) that it is precisely the amount of variance in the treatment explained by the observables in our model which is responsible for bias amplification. Similarly, we can see that a significant change in estimate is not sufficient to suggest that overall bias is decreasing since this could be the result of further bias amplification.\\ These results call for new techniques to be developed for observational studies which can accommodate unmeasured confounding to help researchers choose reasonable and least-biased methods. We suggest to first identify the most plausible causal DAG. From the DAG and basic structural equation assumptions, an expression for asymptotic bias can often be derived. Further, we suggest to estimate the always-identifiable amplification term in observational settings and to assess the risk of bias amplification. With a measure for amplification and a limiting bias expression, a sensitivity analyses can be performed. One reasonable sensitivity analysis approach would be to estimate the amount of unmeasured confounding required in the spirit of E-values VanderWeelePeng_Sensitivity2017 to determine the strength of confounding associations required to "explain away the treatment effect" VanderWeelePeng_Sensitivity2017 and to make principled inferences from the data. This would require, as we have shown, properly simulating the unmeasured confounded as to respect the properties of the original data and such that the other competiting effects, i.e edges of the DAG, are not inadvertently altered. In such a set-up, large effects and relatively small amplifying terms lend credibility to results as being robust to unmeasured confounding, particularly in cases when suitable priors can be placed on the variables along the unmeasured confounding pathway. Alternatively, one could follow the approach of hill_sensitivity and use the underlying structural equations and the data to generate candidate values of the unmeasured confounding. As we showed in section (ref) it is important that any such simulation method take into account the asymmetry of bias amplification with respect to the weight of the edge $\boldsymbol{U} \to \boldsymbol{A}$ and $\boldsymbol{U} \to \boldsymbol{Y}$.\\ Ultimately, simulation experiments must aim to produce data from which we can draw causal conclusions to questions about estimators or functions. This means having well-defined interventions on the edges of the causal graphs and holding the other edges constant. In linear systems of equations, this requires keeping the moments of the variables, in particular variance, fixed when modifying the weight of the DAG's edges. If we allow the treatment variance to vary incidentally as we increase confounding effects, the intervention arm of our simulations will no longer match the target observational study in the control arm. As a further consequence, the additional variance in the exposure may absorb much of the amplifying effect. This leads to systematic underestimation of bias amplification and may be an explanation for why the threat of bias amplification has not been appreciated as a concern for applied researchersmeyers2011. Fixing the variance of the variables has the additional benefit of defining the feasible parameter space. By constraining the underlying parameters by the implied variance equations (e.g equation ((ref))), it is computationally and conceptually easier to simulate the entire range of plausible treatment effects and biases. This leads to more representative simulations and more principled inferences.
\nocite{*}