Model Selection for Causal Modeling in Missing Exposure Problems

Based on the counterfactual causal framework ¹⁷, we denote $(Y_{i}^{1},Y_{i}^{0}),i=1,\dots,n,$ as the potential outcomes if individual $i$ were exposed or unexposed (or equivalently, treated or untreated), respectively. Let $\bm{X}_{i}$ denote the covariates, $A_{i}$ denote the exposure of interest, $Y_{i}$ denote the binary outcome, and $R_{i}=I\{A_{i}\text{ is missing}\}$ denote the indicator of whether the exposure value is missing or not for individual $i=1,\dots,n$. In this study, we focus on the exposure variable being missing, so all other variables are assumed to be completely observed. In addition, we study the model selection problem on a finite set of covariates instead of the high-dimensional setting.

For each individual $i$, $Y_{i}^{1}$ and $Y_{i}^{0}$ cannot be observed at the same time; instead, we only observe $Y_{i}$ and $A_{i}$. The observed dataset is $(\bm{X}_{i},A_{i},Y_{i},R_{i}),i=1,\dots,n$, and the relationship between the observed and potential outcomes can be written as $Y_{i}=A_{i}Y_{i}^{1}+(1-A_{i})Y_{i}^{0}$, for $i=1,\dots,n$. The true propensity score (PS) is defined as $\pi(\bm{X}_{i})=P(A_{i}=1|\bm{X}_{i})$ for individual $i=1,\dots,n$. We denote $\tau_{1}=E(Y^{1})$ and $\tau_{0}=E(Y^{0})$ as the average potential outcomes if treated or untreated, respectively. For the binary outcome in our application studies, we define the causal estimand of interest as the causal risk ratio: $\tau=\tau_{1}/\tau_{0}$. Even though the study focuses on the binary outcome, the theory is developed in general scenarios and the key ideas presented in this article will not be largely changed by the different formats of the causal quantities, so the application can be easily expanded to other cases, such as the continuous outcome or count outcome. To estimate the causal risk ratio $\tau$, three assumptions are required as follows:

1. 1.

   Strongly ignorable treatment assignment assumption (SITA): $(Y^{1},Y^{0})\bot A|\bm{X}$.

2. 2.

   Positivity assumption: $0<P(A=1|\bm{X}=\bm{x})<1,\text{for all possible }\bm{x}$.

3. 3.

   MAR assumption: $R\bot A|\bm{X},Y$.

Here, SITA assumption is also known as the assumption of no unmeasured confounders, which cannot be tested in the application studies. The positivity assumption requires the propensity score to be a positive probability, which can be checked after imputing the missing exposure status correctly ¹⁷. MAR assumption assumes that the missing indicator is conditionally independent of the exposure itself given all observed covariates and outcomes in the dataset, which also cannot be tested on the observed data due to the unknown missing values ³.

The objective of this study is two-fold. Firstly, since different components of $\bm{X}$ can be different as shown in Figure 1, our primary goal is to find the best selection on both imputation and PS models when the exposure is MAR. Secondly, we aim to find a suitable criterion to evaluate the performance of the two models in the application study when a DAG is not known.

To study the different roles of covariates in model selection, we consider a simplified scenario with three covariates, which have different effects on the exposure or the outcome. As shown in Figure 1, $X_{1}$ is the main confounder, $X_{2}$ is the exposure-related covariate, and $X_{3}$ is the covariate related only to the outcome. We denote the $4\times 1$ vector of covariates for subject $i$ as $\bm{X}_{i}=(1,X_{i1},X_{i2},X_{i3})^{T}$. In terms of assumptions described in the general scenario in Section 2.1, we can rewrite three main assumptions based on Figure 1 as follows:

1. 1.

   Strongly ignorable treatment assignment assumption (SITA): $(Y^{1},Y^{0})\bot A|X_{1}$.

2. 2.

   Positivity assumption: $0<P(A=1|\bm{X}=\bm{x})<1,\text{for all possible }\bm{x}$.

3. 3.

   MAR assumption: $R\bot A|X_{1},X_{2}$.

Since $X_{1}$ is the “sufficient set” for confounding adjustment in this specific case ¹⁸, SITA assumption only requires conditioning on the main confounder $X_{1}$ instead of all covariates. In addition, we do not include the association between $R$ and $Y$ as shown in Figure 1 since we want to investigate whether adding $X_{3}$ and $Y$ will increase the accuracy of the imputation model for the missing exposure. In that way, we only need to condition on $X_{1},X_{2}$ for MAR assumption based on Figure 1.

To investigate the role of each variable in the imputation model, one of the approaches is to rely on the theories of directed acyclic graph (DAG), which is a common tool using “d-separation” ¹⁹ to check conditional independence among variables without specifying the form of the models. In our specific example in Figure 1, we know that $X_{1},X_{2}$ should be included in the imputation model because they affect $A$ and $R$ directly. The main question is whether $Y$ should be included in the imputation model or not. Even though we do not include the outcome in the true missingness model in Figure 1, $Y$ is directly correlated with the exposure variable, so adding $Y$ may help predict the missing exposure values.

Next, we aim to study the role of $X_{3}$ in the imputation model. Notice that $X_{3}\bot A$, but if given $(X_{1},X_{2})$, $Y$ is a function of $X_{3}$ and $R$ is also correlated with $X_{3}$, as shown in Figure 1. In addition, if we condition on $Y$, $Y$ will be the collider in the path of $A\rightarrow Y\leftarrow X_{3}$, so $X_{3}\not\!\perp A|(Y,X_{1},X_{2})$. In other words, including the outcome and outcome-related variables in the imputation model is expected to improve the predictive ability of the imputation model. In summary, we should include all $(X_{1},X_{2},X_{3},Y)$ into the imputation model, which will be verified later in simulation studies.

Before we discuss the model selection strategy, we first discuss how to estimate $\tau$. As we have previously mentioned, estimating $\tau$ requires us to account for both missing and confounding issues. For the missing data problem, we will impute the exposure status via MICE based on MAR assumption ⁷, ⁸. After the missing exposures are imputed, $\tau$ can be estimated using inverse-probability weighting (IPW) or double robust (DR) method to account for the confounding issues.

In summary, we present an algorithm with three key steps as follows:

1. 1.

   Fit imputation (Imp) model based on MAR assumption using MICE: fit logistic regression model for $A\sim\bm{X}+Y$ to obtain $A_{i}^{\text{imp}}$ as the exposure on the imputed dataset, which includes the original exposure for individuals without missing values and imputed exposure for individuals with missing values.

2. 2.

   Fit PS model using the logistic model on the imputed dataset: $A^{\text{imp}}\sim\bm{X}$ to obtain fitted PS values, denoted as $\hat{\pi}_{i}(\bm{X})$.

3. 3.

   Apply inverse weighting to adjust for the confounding issue, and the IPW estimator can be written as ²⁰:

   $\displaystyle\hat{\tau}_{1}^{\text{IPW}}=\frac{1}{n}\sum_{i=1}^{n}\frac{A_{i}^{\text{imp}}Y_{i}}{\hat{\pi_{i}}(\bm{X}_{i})},\ \ \hat{\tau}_{0}^{\text{IPW}}=\frac{1}{n}\sum_{i=1}^{n}\frac{(1-A_{i}^{\text{imp}})Y_{i}}{1-\hat{\pi_{i}}(\bm{X}_{i})}$

   (1)

   Then, an IPW estimator for $\tau$ is: $\hat{\tau}^{\text{IPW}}=\hat{\tau}_{1}^{\text{IPW}}/\hat{\tau}_{0}^{\text{IPW}}$ or the DR estimator can be written as:

   $\displaystyle\hat{\tau}_{1}^{\text{DR}}=\frac{1}{n}\sum_{i=1}^{n}\left[\frac{A_{i}^{\text{imp}}Y_{i}}{\hat{\pi_{i}}(\bm{X}_{i})}-\frac{\hat{\pi_{i}}(X)-A_{i}^{\text{imp}}}{\hat{\pi_{i}}(\bm{X}_{i})}\hat{m_{1}}(A_{i}^{\text{imp}},\bm{X}_{i})\right],\ \ \hat{\tau}_{0}^{\text{DR}}=\frac{1}{n}\sum_{i=1}^{n}\left[\frac{(1-A_{i}^{\text{imp}})Y_{i}}{1-\hat{\pi_{i}}(\bm{X}_{i})}-\frac{A_{i}^{\text{imp}}-\hat{\pi_{i}}(X)}{1-\hat{\pi_{i}}(\bm{X}_{i})}\hat{m_{0}}(A_{i}^{\text{imp}},\bm{X}_{i})\right]$

   (2)

   where $\hat{m_{1}}(A_{i}^{\text{imp}},\bm{X}_{i})$ is the fitted response for the treatment group, i.e., $\hat{m_{1}}(A_{i}^{\text{imp}},\bm{X}_{i})=E[Y_{i}|A_{i}^{\text{imp}}=1,\bm{X}_{i};\hat{\bm{\beta}}_{1}]$; $\hat{m_{0}}(A_{i}^{\text{imp}},\bm{X}_{i})$ is the fitted response for the control group, i.e., $\hat{m_{0}}(A_{i}^{\text{imp}},\bm{X}_{i})=E[Y_{i}|A_{i}^{\text{imp}}=0,\bm{X}_{i};\hat{\bm{\beta}}_{0}]$. Here, $(\hat{\bm{\beta}_{1}},\hat{\bm{\beta}_{0}})$ are estimated parameters for the treatment or control group from the outcome model, written as $Y\sim A+X$. Then, a DR estimator for $\tau$ is: $\hat{\tau}^{\text{DR}}=\hat{\tau}_{1}^{\text{DR}}/\hat{\tau}_{0}^{\text{DR}}$.

Notice that when we conduct imputation using MICE, we choose $m=20$ as the number of imputations. Then, IPW and DR estimators for $\tau$ can be constructed on each imputed dataset and Rubin’s Rules is applied to obtain final estimated values based on the algorithm described above 2.3 ⁹. If the imputation model is correctly specified, based on MAR assumption, we can approximate $P(A=1|\bm{X},Y)$ by $P(A^{\text{imp}}=1|\bm{X},Y)$. Furthermore, if PS model is also correct, due to SITA assumption, $\hat{\tau}^{\text{IPW}}\xrightarrow[]{\text{p}}\tau$ and its asymptotic normality holds as $n\to\infty$ ²⁰. For DR estimator, it can protect against misspecification of either PS or the outcome model ^{12, 11}. In other words, if the imputation model is correct and either PS or the outcome model is correct, we know $\hat{\tau}^{\text{DR}}\xrightarrow[]{\text{p}}\tau$ and its asymptotic normality holds as $n\to\infty$ ².

In application studies, since the true causal effect is unknown, we are not able to find which combination of imputation and PS models leads to the best performance in terms of some performance metrics, such as RMSE. In such a case, we need some model selection criteria to choose the best combination of imputation and PS models. In this section, we discuss some traditional criteria for model selection and propose a new criterion that takes into consideration both models.

To investigate which covariates should be included in the imputation and PS models, a simulation study is conducted based on the causal diagram in Figure 1. The number of replications is $N=1000$, and the sample size is $n=500$ in each data generation. Three covariates $X_{1},X_{2},X_{3}$ are generated from $N(0,1)$ independently. $Y,A$, and the missing indicator $R$ are generated by the following three models:

• •

  missingness model: $\text{logit}\{P(R=1|\bm{X})\}=-0.3+0.4X_{1}+0.6X_{2}+1.8X_{3}$;

• •

  treatment model: $\text{logit}\{P(A=1|\bm{X})\}=-0.2+0.3X_{1}+1X_{2}$;

• •

  outcome model: $\text{logit}\{P(Y=1|A,\bm{X})\}=-0.2+2A-0.3X_{1}+2.5X_{3};$

The missing rate is about 48% in this scenario and the true causal effect $\tau=1.523$. Notice that the true missingness model does not include the outcome $Y$, but one of our aims is to investigate whether adding the outcome into the imputation model will improve the performance of the estimated causal effect in the simulation studies.

For each simulated data set, we investigate five different imputation models and four different PS models, which we have specified and named in Tables 3.1, respectively. In total, twenty possible combinations of the imputation and PS models are investigated. We provide some rationale behind the specification of these imputation and PS models. When selecting the imputation model, we start with an exposure-related model, which includes only $X_{1}$ and $X_{2}$ as they directly affect $A$. To increase the predictive ability, we gradually modify the imputation model by adding $X_{3}$ and $Y$ in a step-by-step process until the full imputation model is obtained.