Model Selection for Causal Modeling in Missing Exposure Problems

Abstract

[Abstract]In causal inference, properly selecting the propensity score (PS) model is a popular topic and has been widely investigated in observational studies. In addition, there is a large literature concerning the missing data problem. However, there are very few studies investigating the model selection issue for causal inference when the exposure is missing at random (MAR). In this paper, we discuss how to select both imputation and PS models, which can result in the smallest RMSE of the estimated causal effect. Then, we provide a new criterion, called the “rank score" for evaluating the overall performance of both models. The simulation studies show that the full imputation plus the outcome-related PS models lead to the smallest RMSE and the rank score can also pick the best models. An application study is conducted to study the causal effect of CVD on the mortality of COVID-19 patients.

1 Introduction

In real application studies, missing data can occur in multiple ways ¹. One of the common cases is that the exposure of interest may not be fully observed ². For example, a patient’s exposure status may not be fully recorded when the patient suddenly drops out of the clinical study; or an individual may decline to answer a sensitive question regarding his/her health problem in a survey questionnaire. In addressing the challenge of missing data, researchers have proposed various methodologies for handling cases where the outcome is missing³. However, when the exposure status is missing, few approaches have been provided in literature to deal with this issue ^{4, 5}.

In causal inference, propensity score (PS) analysis is a popular tool to address the confounding issue ⁶. However, when the exposure status is missing at random (MAR), estimation of the causal effect is challenging. The main reason is that MAR on exposure significantly differs from MAR on the outcome in this case. When the outcome is missing, only the outcome model is influenced by the missing data, and the PS model will not be affected by the missing outcome. In contrast, in cases where the exposure is missing, all three models — imputation, PS, and outcome models— are largely affected. Consequently, although PS-based methods have been discussed when the outcome is MAR ^{7, 8}, their direct application to cases where the exposure is missing is not straightforward. Careful adjustments regarding both missing and confounding issues are necessary in such scenarios.

As a motivating example, we consider the case in which researchers aim to determine the causal relationship between CVD (as the exposure) and COVID-19 patients’ mortality (as the outcome) in an observational study where CVD status is not completely observed. In this situation, we need to make adjustments for both missing CVD status and other confounding factors. However, which variables should be selected for the imputation and PS models is not clear because we do not know the true missing values. If we only include some clinical variables related to CVD, such as diabetes and metabolic status without including the outcome and outcome-related variables in the imputation and PS models, both models may not be valid, so the imputed exposure may be deviated from the true missing values and extreme weights may occur, leading to a large bias and variance for the estimated causal effect.

One of the common approaches to dealing with missing data is via imputation. Single imputation (SI) is easy to conduct based on the observed data ^{8, 9}. However, imputing data only once may not be reliable and a large variation may be induced ¹⁰. In contrast, the multiple imputations chained equations (MICE) approach is proposed through random sampling, and applying Rubin’s rules can account for both the sampling variability and the uncertainty in the missing values^{11, 12}. After addressing the missing data issue, we can apply either inverse probability weighting (IPW) or double robust (DR) method based on the imputed dataset ^{13, 14, 15}.

If we rely on parametric methods to fit the imputation and PS models, model selection is an essential step to obtain valid results. In the previous literature, researchers have shown that improper model selection for PS model can result in higher bias and variance ¹⁶. For example, if we include variables that are only related to the exposure in the PS model, the variance of estimated causal effects will become larger. On the other hand, if we only include outcome-related variables in the PS model, the estimated causal effects will become more efficient.

However, little attention has been given to addressing situations when both PS and imputation models need to be selected. That motivates us to conduct studies to investigate which variables should be included in the imputation and PS models so that the smallest bias and variance of the causal estimates can be achieved among candidate models. Furthermore, since we do not know the true causal effect in the application, evaluating those selected models through some model selection criteria is also essential. Therefore, we propose a new rank-based criterion to combine the performance of the imputation and PS models, compared with the traditional criterion.

The rest of this article is organized as follows. In Section 2, we describe the goal and the notation; we design a simulation study to investigate model selection in both imputation and PS models. Section 3 discusses the traditional and newly proposed criteria for conducting model selection. In Section 4, an application study is conducted to identify the causal effect of CVD on mortality for a cohort of COVID-19 patients. The article ends with a discussion in Section 5.

2 Notation and Methods

2.1 Framework and Assumptions

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$ is exposed or unexposed (or equivalently, the treated or untreated), respectively. To study the different roles of covariates in model selection, we consider a simplified scenario with three covariates. 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}$. Let $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. 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.

Figure 1: Causal diagram for simulation studies. The black arrows refer to the causal relationship among confounders, the exposure, and the outcome. The double-sided dash arrows refer to the associational relationship among the missing indicator, the outcome, and the covariates. The red arrow is the causal effect of primary interest. The dashed short line refers to no association between the missing indicator and the exposure given covariates and the outcome due to MAR assumption.

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,2,\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 effect for treated and untreated, respectively. We define the causal estimate of interest $\tau$ as the risk ratio: $\tau=\tau_{1}/\tau_{0}$. For continuous outcomes, we can also define the causal effect as the average causal effect, but the key idea for the model selection process will not be largely changed by the different formats of the causal quantities. To estimate the risk ratio, three assumptions are required as follows:

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

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

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

SITA assumption is also known as the assumption of no unmeasured confounder, which cannot be tested in the application studies. The second 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 covariates in the dataset, which also cannot be tested on the observed data due to the unknown missing values.

To investigate the role of each variable in the imputation model, one of the approaches is to use 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 model. From 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 from Figure 1, $Y$ is directly correlated with the exposure variable, so adding $Y$ can help us predict the missing exposure values.

Next, we aim to study the role of $X_{3}$ in the imputation model. Notice that 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}$ from Figure 1. Therefore, $X_{3}$ should be conditioned on for MAR assumption to hold. In addition, there is a new finding between $X_{3}$ and $A$. 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 can largely improve the predictive ability of the imputation model. In summary, even though MAR assumption only requires us to control $(X_{1},X_{2},X_{3})$, we should include $(X_{1},X_{2},X_{3},Y)$ into the imputation model.

The objective of this study is two-fold. Firstly, since the roles of $X_{1},X_{2},X_{3}$ are quite 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.

2.2 Estimation

Before we discuss the model selection strategy, we first provide the algorithm 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. 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 imputed exposure.

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

3. Apply inverse weighting to adjust for the confounding issue, IPW estimator can be written as ¹⁹:

   $$\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, IPW estimator for $\tau$ will be: $\hat{\tau}^{\text{IPW}}=\hat{\tau}_{1}^{\text{IPW}}/\hat{\tau}_{0}^{\text{IPW}}$.

   or the DR estimator can be written as:

   $$\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}}]$; $\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}}]$. Here, $\hat{\bm{\beta}}$ is estimated parameters from the outcome model, written as $Y\sim A+X$. Then, DR estimator for $\tau$ will be: $\hat{\tau}^{\text{DR}}=\hat{\tau}_{1}^{\text{DR}}/\hat{\tau}_{0}^{\text{DR}}$.

If the imputation model is correctly specified, based on MAR assumption, we know $P(A^{\text{imp}}=1|\bm{X},Y)=P(A=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 ^{14, 20}. 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$ ¹⁵.

2.3 A Simulation Study

To investigate which covariates should be included in the imputation and PS models to achieve the smallest RMSE of $\hat{\tau}$, a simulation study is conducted. 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:

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

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

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

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

When we conduct imputation using MICE, multiple imputation times $m=20$ are chosen for imputing the missing data. Then, IPW and DR estimators for $\tau$ can be constructed on each imputed dataset and Rubin’s rule can be applied to obtain final estimated values based on the Section 2.2 ¹².

For each simulated data set, we investigate five different imputation models and four different PS models. In total, twenty possible combinations of the imputation and PS models are investigated. We name these models in Tables 1 and 2:

Table 1: Five imputation models in the simulation studies

Model Name (shortened form)	Model Form
exposure-related model (Exp)	$A\sim X_{1}+X_{2}$
outcome-related model (Out)	$A\sim X_{1}+X_{3}$
all covariates-included model (Covs)	$A\sim X_{1}+X_{2}+X_{3}$
response-included model (Res)	$A\sim X_{1}+X_{2}+Y$
full model (Full)	$A\sim X_{1}+X_{2}+X_{3}+Y$

We provide some rationales 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.

In contrast, the goal of fitting PS model is not to increase the predictive ability, but to adjust for the confounding issue. Therefore, we start with a naive PS model with only $X_{1}$ as the main confounder. Then, we gradually add exposure-related or outcome-related variables until PS model includes all covariates.

Table 2: Four PS models in the simulation studies

Model Name (shortened form)	Model Form
naive model (Naive)	$A^{\text{imp}}\sim X_{1}$
exposure-relate model (Exp)	$A^{\text{imp}}\sim X_{1}+X_{2}$
outcome-relate model (Out)	$A^{\text{imp}}\sim X_{1}+X_{3}$
all covariates-included model (Covs)	$A^{\text{imp}}\sim X_{1}+X_{2}+X_{3}$

2.4 Simulation Results

In this subsection, we investigate which combination of imputation and PS models can lead to the optimal RMSE of $\hat{\tau}$. Table 3 shows simulation results for IPW estimator when $n=500$ and Figure 2 is plotted to visualize the results. We also provide the simulation results for risk ratio and average causal effects in Appendix A.6. To evaluate the performance of different models, we present the bias rate as the main interest: $(\bar{{\tau}}-\tau)/\tau$, where $\bar{{\tau}}=\frac{1}{N}\sum_{j=1}^{N}\hat{{\tau}_{j}}$, and $\hat{{\tau}}_{j}$ is the point estimate in the $j^{th}$ replication, $j=1,2,\dots,N$. The empirical standard error (ESE), which is also called the sampling standard error, is written as $\text{ESE}=\sqrt{\frac{1}{N-1}\sum_{j=1}^{N}(\hat{{\tau}_{j}}-\bar{{\tau}})^{2}}$. The square root of mean squared error (RMSE) is defined as $\text{RMSE}=\sqrt{\frac{1}{N-1}\sum_{j=1}^{N}(\hat{{\tau}}_{j}-\tau)^{2}}$.

Figure 2: Boxplot of estimated values for IPW estimators in the simulation studies when $n=500$. The x-axis refers to which imputation and PS models are used. Here, Exp, Out, Res, and Covs, refer to the exposure-related, outcome-related, response-included, and covariates-included models, respectively. The red line is a true causal effect.

Table 3: Performance of different models using IPW method when $n=500$

Imp, PS Models	Bias	Bias Rate	ESE	RMSE
Exp, naive	-0.186	-12.187	0.085	0.204
Exp, exp	-0.217	-14.283	0.092	0.236
Exp, out	-0.188	-12.324	0.067	0.199
Exp, covs	-0.219	-14.384	0.080	0.233
Out, naive	-0.290	-19.069	0.109	0.310
Out, exp	-0.333	-21.858	0.103	0.348
Out, out	-0.230	-15.074	0.059	0.237
Out, covs	-0.263	-17.301	0.059	0.270
Covs, naive	-0.181	-11.920	0.123	0.219
Covs, exp	-0.211	-13.883	0.138	0.252
Covs, out	-0.188	-12.335	0.067	0.199
Covs, covs	-0.218	-14.295	0.080	0.232
Res, naive	0.102	6.679	0.176	0.204
Res, exp	0.137	9.026	0.217	0.257
Res, out	-0.046	-2.996	0.107	0.116
Res, covs	-0.038	-2.523	0.137	0.143
Full, naive	-0.017	-1.096	0.160	0.161
Full, exp	-0.011	-0.710	0.194	0.194
Full, out	0.010	0.640	0.110	0.111
Full, covs	0.023	1.498	0.137	0.139

• Here, Exp, Out, Res, and Covs, refer to the exposure-related, outcome-related, response-included, and covariates-included models, respectively.

2.4.1 Comparison of Imputation Models

Based on Table 3 and Figure 2, if we focus on the same PS model and compare different imputation models, we find that the full imputation model has resulted in the smallest bias of $\hat{\tau}$ compared to other imputation models. In contrast, the exposure-related, outcome-related, and all covariates-included imputation models generate huge bias rates.

Since $(X_{1},X_{2})$ are exposure-related variables from Figure 1, some researchers may argue that the preferred imputation model only needs to include $(X_{1},X_{2})$ without $(X_{3},Y)$ . However, our simulation results indicate that adding $(X_{3},Y)$ into the imputation model is necessary to reduce the bias of the estimated causal effect. For example, comparing the exposure-related model with the response-included imputation model, Table 3 shows that adding $Y$ into the imputation model reduces about 10% of bias. To study the effect of $X_{3}$, if we compare the exposure-related model versus the covariate-included imputation model, adding $X_{3}$ into the imputation model does not reduce the large bias if we do not condition on $Y$. In contrast, when we condition on $Y$, compared with the response-included imputation versus the full imputation model, Figure 2 shows that adding $X_{3}$ leads to a reduction of bias. These findings are consistent with the discussion of DAG in Section 2, which indicate that even though the true missing model does not include $Y$ and $A$ is not directly correlated with $X_{3}$, adding both $X_{3}$ and $Y$ into the imputation model can improve the predictive ability.

In summary, even though the true missingness model does not include $Y$, all exposure-related, outcome-related covariates and the outcome itself should be included in the imputation model, so the full imputation model results in the smallest bias among others.

2.4.2 Comparison of PS Models

Based on the above discussion, we next focus on the full imputation model and compare the performance of different PS models in Table 3 and Figure 2. We find that all PS models that incorporate the main confounder $X_{1}$ can lead to very small bias, as $X_{1}$ is the "sufficient set" for confounding adjustment, in this case, ²¹. Even though $X_{2}$ directly causes $A$, adding $X_{2}$ does not reduce the bias because it only causes the change of the exposure and has no direct impact on the outcome $Y$, so $X_{2}$ is not considered as the main confounder.

In terms of standard errors, the outcome-related PS model with $X_{3}$ can generate the smallest ESE when we fit the full imputation model. In contrast, both exposure-related and all-covariates PS models lead to an increase in ESE because adding exposure-related variable $X_{2}$ may generate more extreme PS values to enlarge the variation in finite samples ¹⁶. In summary, we recommend choosing the outcome-related PS model over the exposure-related PS model.

The results using DR estimator are presented in Table 4 and visualized in Figure 3. Overall, we observe a similar trend of model selection using DR estimator as IPW estimator, suggesting that full imputation plus outcome-based PS models can effectively reduce RMSE of $\hat{\tau}$ compared to other candidate models.

Overall, the best model selection is to choose the full imputation plus outcome-related PS models for both IPW or DR estimator, which results in the smallest RMSE of $\hat{\tau}$. In contrast, using outcome-related imputation plus exposure-related PS models leads to the largest RMSE, which is not recommended.

Figure 3: Boxplot of estimated values for DR estimators in the simulation studies when $n=500$. The x-axis refers to which imputation and PS models are used. Here, Exp, Out, Res, and Covs, refer to the exposure-related, outcome-related, response-included, and covariates-included models, respectively. The red line is the true causal effect.

Table 4: Performance of different models using DR method when $n=500$

Imp, PS Models	Bias	Bias Rate	ESE	RMSE
Exp, naive	-0.186	-12.198	0.084	0.204
Exp, exp	-0.222	-14.553	0.085	0.237
Exp, out	-0.188	-12.374	0.065	0.199
Exp, covs	-0.225	-14.804	0.066	0.235
Out, naive	-0.291	-19.087	0.109	0.310
Out, exp	-0.321	-21.056	0.103	0.337
Out, out	-0.226	-14.859	0.058	0.233
Out, covs	-0.264	-17.354	0.057	0.270
Covs, naive	-0.182	-11.933	0.123	0.219
Covs, exp	-0.215	-14.105	0.134	0.253
Covs, out	-0.189	-12.392	0.065	0.200
Covs, covs	-0.225	-14.786	0.067	0.235
Res, naive	0.102	6.679	0.176	0.204
Res, exp	0.133	8.763	0.211	0.250
Res, out	-0.056	-3.654	0.100	0.114
Res, covs	-0.066	-4.350	0.115	0.132
Full, naive	-0.017	-1.102	0.160	0.161
Full, exp	-0.015	-0.977	0.189	0.189
Full, out	-0.010	-0.660	0.104	0.104
Full, covs	-0.011	-0.691	0.120	0.120

• Here, Exp, Out, Res, and Covs, refer to the exposure-related, outcome-related, response-included, and covariates-included models, respectively.

3 Model Selection Criteria

In application studies, the true causal effect is unknown, so we are not able to find which combination of models leads to the best performance. In such a case, we need some model selection criteria to choose the best 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.

3.1 Weighted Accuracy ($\text{Accuracy}^{(w)})$

We first discuss a theoretical way to evaluate the performance of the imputation model via weighted accuracy based on the observed data. Due to the missing exposure, even though we can impute the missing data, we do not know the true missing values, which makes it challenging to estimate the accuracy only based on the observed data. A naive way is to fit the imputation model based on the observed data and compare the imputed exposure with those known values. However, since the accuracy is a function of $(\bm{X},Y)$ and the distribution of $(\bm{X},Y)$ on observed data has been changed on the whole data, the way to estimate accuracy from the observed data is no longer valid.

One of the approaches to deal with this issue is to apply inverse probability weighting on the accuracy when we consider $w_{i}=P(R_{i}=1|\bm{X_{i}},Y_{i})$ as the propensity for the missingness. Obtaining the weighted accuracy can be described in the following four steps:

1. First, we randomly split the individuals into either the training or the testing data according to a ratio $q$, where $q=\frac{n_{\text{test}}}{n_{\text{obs}}}$. Here, $n_{\text{obs}},n_{\text{test}}$ is the sample size for the whole observed data and for the testing data, respectively. Certainly, $q$ can be arbitrarily chosen by the user, but the simplest way is to set the ratio equal to the missing rate in the original dataset.

2. Next, we fit the imputation model on the training data and impute the exposure on the testing data to compare with the known exposure status.

3. After that, we fit the full missingness model on the whole dataset, written as $\text{logit}(w_{i})=\bm{X_{i}^{T}\gamma}$, where $\bm{\gamma}$ is the vector of coefficients including the intercept to be estimated so that we can obtain $\hat{w_{i}}$ as the estimated propensity of missingness for individual $i=1,2,\dots,n$.

4. Finally, we calculate the weighted accuracy, called “$\text{Accuracy}^{(w)}$":

   $$\text{Accuracy}^{(w)}(A^{\text{imp}})=\frac{\sum_{i=1}^{n}\frac{1-R_{i}}{1-\hat{w}_{i}}\delta_{i}\mathbbm{1}(A_{i}=A^{\text{imp}})}{n\times q},$$ (3)

where $\delta_{i}=\mathbbm{1}(\text{individual $i$ is chosen in the test set})$ and $\mathbbm{1}(A_{i}=A^{\text{imp}})$ is an indicator for whether the imputed value equals to the observed value. Here, $\hat{w}_{i}$ is estimated from the full missingness model including all covariates and the outcome.

The weighted accuracy is consistent to the true accuracy when the full missingness model is correctly specified. The detailed proof is attached in Appendix A.7.1. We also verify this conclusion in the simulation study, in which we find that the difference between the true accuracy and the weighted accuracy is close to 0.01 (the results are omitted).

3.2 ASMD, KS, and BIC

To evaluate the PS model, the traditional approaches are based on the balance statistic in the confounders. For example, one can calculate absolute standardized mean difference (ASMD) and Kolmogorov-Smirnov (KS) statistics in the covariates after propensity score adjustment.

Another approach to evaluate the PS model is using BIC criterion to select the best model adjusting for the sample size and number of variables in the model. Since in PS model selection, we aim to select the outcome-related variables instead of the exposure-related variables, we suggest using BIC of the outcome model $Y\sim A+\bm{X}$ as the selection criterion. From Figure 1, since $X_{2}\bot Y|A=a$, a smaller BIC will indicate that we have included outcome-related variables and excluded the exposure-related variables. For example, based on the simulation results in Table 5, the naive and exposure-related PS models generally have larger BIC than outcome-related and covariates-included PS models. Notice that we do not recommend using the c-statistics or accuracy of PS model to select variables because in PS analysis, we aim to balance the confounder across the exposure and non-exposure group instead of maximizing the predictability of the PS model ^{22, 23}.

We should also notice that only using one traditional criterion may not be appropriate to select both imputation and PS models, such as either $\text{Accuracy}^{(w)}$, ASMD, or KS statistics, because they just focus on the performance of a single model. For example in Table 5, if we just adopt $\text{Accuracy}^{(w)}$ to select imputation and PS models, we cannot distinct between different PS models because the same imputation model will lead to the same $\text{Accuracy}^{(w)}$ value. In contrast, if we only choose BIC to select two models and we choose the same exposure-related PS model, we cannot easily distinguish the performance of the outcome-related vs covariates-included imputation models. That motivated us to combine the evaluation of $\text{Accuracy}^{(w)}$ and BIC into an integrated criterion, called rank score proposed in the next section.

3.3 Rank Score

The idea of the new criterion is to combine the performance of both imputation and PS models, so the “rank score" is proposed to average over the ranks of $\text{Accuracy}^{(w)}$ and BIC among all candidates of models. We can also think of the rank score similar as the “best subset" strategy, which considers all possible combination of the candidate models and select the best choice in contrast to the “stepwise" strategy. For a given model, we can calculate its rank score value by:

$$\text{ Rank Score}=\frac{\text{Rank(1-$\text{Accuracy}^{(w)}$)+Rank(BIC)}}{2},$$ (4)

where “Rank($1-\text{Accuracy}^{(w)}$)" is the rank of the model based on the value of $1-\text{Accuracy}^{(w)}$ if we order $1-\text{Accuracy}^{(w)}$ from the largest to the smallest. “Rank(BIC)" is the rank for a given model based on the value of BIC for regressing $Y$ on $A$ and the covariates selected in the given model. We can calculate the rank score for every possible combination of the imputation and PS models. Then, we select the smallest rank score, which leads to the highest $\text{Accuracy}^{(w)}$ and the smallest BIC, so the best imputation and PS models can be successfully chosen. From the simulation studies in Table 5, the full imputation plus outcome-related PS models will result in the smallest rank score as the best choice, consistent with the smallest RMSE. In comparison, if we choose the outcome imputation plus exposure-related model, the rank score has the largest value, which is also consistent with the highest RMSE.

In summary, the main advantage of the rank score is to combine the performance of both imputation and PS models and directly find the best model based on this rank-based criterion. In addition, the rank score is a unit-free score, which is not affected by the different magnitudes of the accuracy and BIC. The evaluation of those criteria is shown in Table 7.

Table 5: Criterion values for different models using IPW method when $n=500$

Imp, PS Models	$\text{Accuracy}^{(w)}$	Out BIC	ASMD	KS	ABIC	Rank Score	Rank(RMSE)
Exp, naive	0.586	677.561	0.161	0.130	1.2155	11.789	10
Exp, exp	0.586	720.437	0.022	0.130	1.3491	14.272	14
Exp, out	0.586	430.074	0.149	0.130	0.4498	7.849	7
Exp, covs	0.586	433.966	0.014	0.130	0.4619	8.959	13
Out, naive	0.507	683.862	0.103	0.148	1.8608	15.947	19
Out, exp	0.507	724.226	0.031	0.148	1.9864	18.306	20
Out, out	0.507	434.737	0.075	0.148	1.0901	12.460	15
Out, covs	0.507	435.659	0.008	0.148	1.0929	12.918	18
Covs, naive	0.585	676.668	0.170	0.136	1.2244	11.875	11
Covs, exp	0.585	719.239	0.033	0.136	1.3570	14.354	16
Covs, out	0.585	430.103	0.148	0.136	0.4616	8.043	8
Covs, covs	0.585	433.958	0.014	0.136	0.4736	9.136	12
Res, naive	0.609	652.471	0.182	0.124	0.9677	8.103	9
Res, exp	0.609	694.546	0.040	0.124	1.0989	10.450	17
Res, out	0.609	411.095	0.154	0.124	0.2209	3.993	2
Res, covs	0.609	415.906	0.015	0.124	0.2358	4.580	4
Full, naive	0.622	662.369	0.172	0.138	0.8601	7.463	5
Full, exp	0.622	705.051	0.034	0.138	0.9931	9.941	6
Full, out	0.622	401.525	0.149	0.138	0.0536	1.989	1
Full, covs	0.622	406.367	0.014	0.138	0.0687	2.574	3

• Here, $\text{Accuracy}^{(w)}$: the weighted accuracy on the observed data. Out BIC is BIC of the outcome model. ASMD: absolute standardized mean difference. KS: Kolmogorov-Smirnov distance for $(X_{1},X_{2},X_{3})$. ABIC is the rescaled form of Accuracy and BIC. Rank Score: average ranks of $1-\text{Accuracy}^{(w)}$ and BIC.

Table 6: Criterion values for different models using DR method when $n=500$

Imp, PS Models	$\text{Accuracy}^{(w)}$	Out BIC	ASMD	KS	ABIC	Rank Score	Rank(RMSE)
Exp, naive	0.586	677.561	0.161	0.130	1.2155	11.789	10
Exp, exp	0.586	720.437	0.022	0.130	1.3491	14.272	15
Exp, out	0.586	430.074	0.149	0.130	0.4498	7.849	7
Exp, covs	0.586	433.966	0.014	0.130	0.4619	8.959	14
Out, naive	0.507	683.862	0.103	0.148	1.8608	15.947	19
Out, exp	0.507	724.226	0.031	0.148	1.9864	18.306	20
Out, out	0.507	434.737	0.075	0.148	1.0901	12.460	12
Out, covs	0.507	435.659	0.008	0.148	1.0929	12.918	18
Covs, naive	0.585	676.668	0.170	0.136	1.2244	11.875	11
Covs, exp	0.585	719.239	0.033	0.136	1.3570	14.354	17
Covs, out	0.585	430.103	0.148	0.136	0.4616	8.043	8
Covs, covs	0.585	433.958	0.014	0.136	0.4736	9.136	13
Res, naive	0.609	652.471	0.182	0.124	0.9677	8.103	9
Res, exp	0.609	694.546	0.040	0.124	1.0989	10.450	16
Res, out	0.609	411.095	0.154	0.124	0.2209	3.993	2
Res, covs	0.609	415.906	0.015	0.124	0.2358	4.580	4
Full, naive	0.622	662.369	0.172	0.138	0.8601	7.463	5
Full, exp	0.622	705.051	0.034	0.138	0.9931	9.941	6
Full, out	0.622	401.525	0.149	0.138	0.0536	1.989	1
Full, covs	0.622	406.367	0.014	0.138	0.0687	2.574	3

• Here, $\text{Accuracy}^{(w)}$: the weighted accuracy on the observed data. Out BIC is BIC of the outcome model. ASMD: absolute standardized mean difference. KS: Kolmogorov-Smirnov distance for $(X_{1},X_{2},X_{3})$. ABIC is the rescaled form of Accuracy and BIC. Rank Score: average ranks of $1-\text{Accuracy}^{(w)}$ and BIC.

3.4 Other Criterion

Another possible criterion to combine accuracy and BIC is to directly conduct a product of two values, written as:

$$(1-\text{Accuracy}^{(w)})\times\text{BIC}.$$ (5)

However, since the magnitude of $\text{Accuracy}^{(w)}$ and BIC are very different, the product values are easily affected by the extreme value, so this criterion is not very helpful in selecting both models.

Even though BIC and $\text{Accuracy}^{(w)}$ do not have the same magnitude, we try to rescale both terms to the range of $[0,1]$. We call the following criterion as “ABIC":

$$\text{ABIC}=\left(1-\frac{\text{$\text{Accuracy}^{(w)}$}-\text{min}(\text{$\text{Accuracy}^{(w)}$})}{\text{max(\text{$\text{Accuracy}^{(w)}$})}-\text{min(\text{$\text{Accuracy}^{(w)}$})}}\right)+\frac{\text{BIC}-\text{min}(\text{BIC})}{\text{max(\text{BIC})}-\text{min(\text{BIC})}},$$ (6)

where “min" means the smallest value and “max" means the largest value among all candidates of models. Even though ABIC solves the different magnitude issues, the min and max values of ABIC can still be affected by the extreme value of BIC or $\text{Accuracy}^{(w)}$. Then, the distinction of ABIC values among the candidate models may be largely shrunk to a quite small value, which makes it hard for us to distinguish the model performance, shown in Table 5-6.

3.5 Results of Selection Criteria

To evaluate the performance of different criteria, we also provide the correlation between RMSE and each criterion in Table 7. More specifically, after $N=1000$ replications, we have obtained RMSE for all candidate models. Within each replication of simulation studies, we can evaluate each candidate model with different criteria. Then, in the $j^{\text{th}}$ replication, we can calculate the Spearman correlation between the proposed criterion and RMSE of $\hat{\tau}$. Table 7 shows the average of the Spearman correlations over 1000 replications.

From Table 7, compared with other criteria, the rank score shows a stronger correlation (larger than 0.8) with the RMSE of $\hat{\tau}$ using either IPW or DR method, which indicates that a smaller value of the rank score is more likely to result in a smaller RMSE. In other words, even though we do not know the true causal effect and RMSE in the application, the rank score is still valid for researchers to evaluate the performance of imputation and PS models together. From Table 5 and 6, we can confirm that full imputation plus outcome-related PS models have a higher chance to be selected via the “rank score" criterion.

Table 7: Spearman correlation between criteria and RMSE when $n=500$

Criterion	IPW Method	DR Method
1-$\text{Accuracy}^{(w)}$	0.678	0.655
ASMD	-0.306	-0.337
KS	0.212	0.191
Out BIC	0.682	0.702
ABIC	0.793	0.785
Rank Score	0.841	0.837

• Here, Out BIC is BIC of the outcome model. ASMD: absolute standardized mean difference for $(X_{1},X_{2},X_{3})$. KS: Kolmogorov-Smirnov distance for $(X_{1},X_{2},X_{3})$. ABIC is the rescaled form of Accuracy and BIC. Rank Score: averaged ranks of $1-\text{Accuracy}^{(w)}$ and BIC.

In summary, our proposed rank score combines the performance of both imputation and PS models. As a rank-based criterion, the rank score shows a stronger correlation with the RMSE of $\hat{\tau}$ compared with other traditional criteria. Based on our simulation results, the rank score can also successfully select the full imputation plus the outcome-related PS models with the smallest RMSE as the best choice.

4 Application Study

4.1 Background

COVID-19 is a highly contagious disease that can lead to severe clinical manifestations and long-term complications for patients. Some researchers have shown that the presence of cardiovascular disease (CVD) is associated with significantly worse mortality in patients ²⁴. However, the associational analysis cannot answer the question of whether CVD directly causes higher mortality or whether other risk factors that lead to CVD are also a common cause of COVID-19 death. In this section, we present an application study on a COVID-19 dataset to investigate the causal effect of CVD (as the exposure) on the mortality of COVID-19 patients (as the outcome) when CVD status is not fully observed. We aim to select both imputation and PS models using the proposed methods and compare their performance.

The data were collected during the second wave of the pandemic in Brazil from June 25th to December 31st, 2020. The dataset contains 2878 COVID-19 patients, of whom 1253 died from COVID-19. Although all patients were exposed to COVID-19, 548 patients had missing information on CVD, resulting in a missing rate of 19%. Demographic and clinical variables were also collected including age, sex, and diabetes without missing data. A summary of the dataset is provided in Table 8-9, stratified by the mortality or the missingness. The independence between mortality and covariates was tested using Fisher test for categorical variables and the Mann-Whitney test for continuous variables. From Table 8, we find a significant association between age and mortality, which is a well-known risk factor. There are also significant associations between CVD and mortality, which may be influenced by both confounders and missingness in the exposure. However, no significant association is identified between diabetes, sex, and mortality.

Table 8: Summary table of COVID-19 data stratified by mortality

Covariate	Full Sample (n=2878)	Alive (n=1625)	Died (n=1253)	p-value
CVD				0.022
No	802 (34)	478 (36)	324 (32)
Yes	1528 (66)	835 (64)	693 (68)
Missing	548	312	236
age				$<$0.001
Mean (sd)	68.9 (16.3)	65.4 (16.6)	73.3 (14.7)
Median (Min,Max)	72 (0,107)	67 (0,104)	76 (0,107)
sex				0.054
Male	1551 (54)	850 (52)	701 (56)
Female	1327 (46)	775 (48)	552 (44)
diabetes				0.3
No	1117 (39)	617 (38)	500 (40)
Yes	1761 (61)	1008 (62)	753 (60)

Table 9: Summary table of COVID-19 data stratified by the missing indicator

Covariate	Full Sample (n=2878)	Missing (n=2330)	Observed (n=548)	p-value
age				0.099
Mean (sd)	68.9 (16.3)	69 (16.6)	68.3 (14.8)
Median (Min,Max)	72 (0,107)	72 (0,104)	71 (1,107)
sex				0.37
Male	1551 (54)	1246 (53)	305 (56)
Female	1327 (46)	1084 (47)	243 (44)
diabetes				$<$0.001
No	1117 (39)	1111 (48)	6 (1)
Yes	1761 (61)	1219 (52)	542 (99)

In this application, we aim to estimate the causal risk ratio of mortality for patients with CVD versus those without CVD as our main interest in the application study. However, in such observational studies, we cannot simply randomize the exposure status for patients, making it challenging to estimate the causal effect of CVD on COVID-19 mortality in the presence of confounding issues. Meanwhile, since CVD status is not fully observed, we also have to deal with the missing exposure problem, so we need to make assumptions about the missing exposure and confounding issues. From Table 9, we find a significant association between missingness and diabetes. Older COVID-19 patients with chronic diseases, such as diabetes, may be more likely to report whether they have CVD or not, so assuming MCAR in this study is not reasonable. Instead, we adopt the MAR assumption, which assumes that the missing indicator for CVD is conditionally independent of CVD status itself, given other clinical variables. In addition, we assume SITA assumption holds, i.e. the potential outcomes are conditionally independent of the CVD given all pre-measured covariates. After that, we intend to apply MICE to impute the missing values and IPW-based estimators to account for confounding issues.

4.2 Model Selection

Since the effect of each clinical variable on CVD status or the outcome may be different, properly selecting imputation and PS models is essential to obtain a valid estimation of the causal effect. However, a DAG in this case is not readily available if we lack sufficient background information in the application. Previous clinical papers show that age and sex are important risk factors for death ^{25, 26}. In addition, many clinical studies find that diabetes and age can affect the risk of CVD ²⁷. Based on the prior knowledge, we can draw a hypothetical causal diagram in Figure 4.

Figure 4: A hypothetical causal diagram for application study. The black arrows refer to the causal relationship among confounders, the exposure, and the outcome. The double-sided dash arrows refer to the associational relationship among the missing indicator, the outcome, and the covariates. The red arrow is the effect of primary interest. The dashed short line refers to no association between the missing indicator and the exposure given covariates and the outcome due to MAR assumption.

For simplicity, we rename the age as $\bm{X_{1}}$, diabetes as $X_{2}$, sex as $X_{3}$, CVD as the exposure $A$, and death as the outcome $Y$. Then, following the same model setup in the simulation study described in Section 2.3, we consider five different imputations and four different PS models as candidate models.

Table 10: Results of model selection for COVID-19 data using IPW and DR methods

Imp Model	PS Model	IPW Estimate	DR Estimate	$\text{Accuracy}^{(w)}$	Out BIC	ASMD	KS	Rank Score
$A\sim\bm{X_{1}}+X_{2}$	$A^{\text{imp}}\sim X_{1}$	0.945	0.986	0.607	3790.876	0.050	0.292	7.500
$A\sim\bm{X_{1}}+X_{2}$	$A^{\text{imp}}\sim\bm{X_{1}}+X_{2}$	0.961	0.988	0.607	3796.287	0.006	0.290	12.000
$A\sim\bm{X_{1}}+X_{2}$	$A^{\text{imp}}\sim\bm{X_{1}}+X_{3}$	0.950	0.991	0.607	3789.885	0.048	0.292	5.000
$A\sim\bm{X_{1}}+X_{2}$	$A^{\text{imp}}\sim\bm{X_{1}}+X_{2}+X_{3}$	0.963	0.989	0.607	3795.815	0.008	0.290	9.500
$A\sim\bm{X_{1}}+X_{3}$	$A^{\text{imp}}\sim X_{1}$	0.946	0.983	0.581	3790.743	0.042	0.267	12.000
$A\sim\bm{X_{1}}+X_{3}$	$A^{\text{imp}}\sim\bm{X_{1}}+X_{2}$	0.960	0.984	0.581	3796.138	0.006	0.264	16.500
$A\sim\bm{X_{1}}+X_{3}$	$A^{\text{imp}}\sim\bm{X_{1}}+X_{3}$	0.951	0.988	0.581	3789.763	0.040	0.267	9.500
$A\sim\bm{X_{1}}+X_{3}$	$A^{\text{imp}}\sim\bm{X_{1}}+X_{2}+X_{3}$	0.962	0.986	0.581	3795.679	0.007	0.264	14.500
$A\sim\bm{X_{1}}+X_{2}+X_{3}$	$A^{\text{imp}}\sim X_{1}$	0.942	0.984	0.606	3790.866	0.050	0.294	11.000
$A\sim\bm{X_{1}}+X_{2}+X_{3}$	$A^{\text{imp}}\sim\bm{X_{1}}+X_{2}$	0.958	0.985	0.606	3796.309	0.006	0.291	16.500
$A\sim\bm{X_{1}}+X_{2}+X_{3}$	$A^{\text{imp}}\sim\bm{X_{1}}+X_{3}$	0.947	0.990	0.606	3789.882	0.049	0.294	8.500
$A\sim\bm{X_{1}}+X_{2}+X_{3}$	$A^{\text{imp}}\sim\bm{X_{1}}+X_{2}+X_{3}$	0.960	0.986	0.606	3795.837	0.008	0.291	14.000
$A\sim\bm{X_{1}}+X_{2}+Y$	$A^{\text{imp}}\sim X_{1}$	0.942	0.985	0.606	3790.846	0.050	0.293	8.500
$A\sim\bm{X_{1}}+X_{2}+Y$	$A^{\text{imp}}\sim\bm{X_{1}}+X_{2}$	0.959	0.986	0.606	3796.272	0.006	0.291	13.500
$A\sim\bm{X_{1}}+X_{2}+Y$	$A^{\text{imp}}\sim\bm{X_{1}}+X_{3}$	0.947	0.990	0.606	3789.852	0.049	0.293	6.000
$A\sim\bm{X_{1}}+X_{2}+Y$	$A^{\text{imp}}\sim\bm{X_{1}}+X_{2}+X_{3}$	0.961	0.987	0.606	3795.797	0.008	0.291	11.000
$A\sim\bm{X_{1}}+X_{2}+X_{3}+Y$	$A^{\text{imp}}\sim X_{1}$	0.946	0.986	0.608	3790.742	0.048	0.289	3.500
$A\sim\bm{X_{1}}+X_{2}+X_{3}+Y$	$A^{\text{imp}}\sim\bm{X_{1}}+X_{2}$	0.960	0.987	0.608	3796.144	0.006	0.287	9.000
$A\sim\bm{X_{1}}+X_{2}+X_{3}+Y$	$A^{\text{imp}}\sim\bm{X_{1}}+X_{3}$	0.950	0.991	0.608	3789.735	0.048	0.289	1.000
$A\sim\bm{X_{1}}+X_{2}+X_{3}+Y$	$A^{\text{imp}}\sim\bm{X_{1}}+X_{2}+X_{3}$	0.962	0.988	0.608	3795.661	0.008	0.287	6.000

Without assuming the hypothetical DAG is correct, we apply the rank score to select the best imputation and PS models among the candidates. Table 10 reports the values for each criterion and based on rank score, the best imputation model includes “age", “sex", “diabetes", and “death" as the predictors and the best PS model selects “age", “sex", and “diabetes" as the predictors. We find that the above selection result leads to the smallest rank score, which is consistent with the simulation results. In addition, if we do not know DAG in Figure 4 is correct or not, using rank score can help us select the full imputation model plus the outcome-related PS model as the best choice.

4.3 Estimation Results

After the model selection, we obtain the estimated causal effect of CVD on mortality using IPW and DR estimators. We also provide confidence intervals (CI) after running $B=2000$ bootstrap resamples, as shown in Table 11. We compare two different types of CI including a CI based on bootstrap standard error (BSE) with a normal approximation, called “CI Normal", and another CI based on bootstrap percentiles (between 2.5% and 97.5% percentiles of bootstrap point estimates), called “CI Percentile" ²⁸.

The estimated causal risk ratio of CVD on the mortality for COVID-19 patients is 0.95 using IPW method with a 95% CI of [0.867, 1.052] CI Percentile. Using DR method, we also obtain a similar result as IPW method.

Although some researchers may think that patients with CVD status increase the risk of death from the common sense, our analysis does not support this assumption. Since CI covers the value of 1, we would expect that the statistical test for whether the ratio of patients with CVD is higher than those without CVD is not significant after adjusting for missing data and confounding factors. One possible explanation for this finding is due to the confounding issue, i.e. age is strongly associated with both CVD status and death. Additionally, researchers may have ignored the effect of missing exposure on the estimated results. These results align with recent clinical research ^{29, 30}, which suggests that those CVD risk factors, such as age, rather than CVD status itself, are the primary contributors to mortality.

One of the limitations of the application study is that we only include four clinical variables. Since we lack enough observations of other variables in the real dataset, we may miss some confounders in the analysis. In addition, the model selection and estimated results are only based on the samples in Brazil instead of the whole population in the world. These issues can be addressed if we have individual-level health data with better quality.

Table 11: Estimated results using IPW and DR estimators for COVID-19 data

Method	Selected Imputation, PS Models	Estimate (risk ratio)	BSE	CI Normal	CI Percentile
IPW	$A\sim\bm{X_{1}}+X_{2}+X_{3}+Y$, $A^{\text{imp}}\sim\bm{X_{1}}+X_{3}$	0.950	0.048	(0.856,1.044)	(0.867,1.052)
DR	$A\sim\bm{X_{1}}+X_{2}+X_{3}+Y$, $A^{\text{imp}}\sim\bm{X_{1}}+X_{3}$	0.991	0.050	(0.897,1.085)	(0.905,1.098)

• Estimate: causal effect (as the risk ratio) of CVD on mortality. BSE: bootstrap standard errors. CI-Normal: 95% bootstrap normal confidence interval using BSE. CI-Percentile: 95% bootstrap percentile confidence interval.

5 Discussion

Model selection on both imputation and PS models when the exposure is missing at random (MAR) is a challenging problem, and few studies discuss the appropriate criterion for selecting these models in practice. In this paper, our simulation studies investigate the effect of different imputation and PS models on the estimated causal effect. Unexpectedly, we find that even though in the simulation setup the missing mechanism does not depend on the outcome, selecting the outcome and outcome-related variables for both the imputation and PS models can improve the accuracy of the imputed exposure and reduce the RMSE of the estimated causal effect. In addition, we propose “rank score" as a rank-based and unit-free criterion to evaluate the performance of both imputation and PS models, which shows a higher correlation with RMSE in simulation studies.

The limitation of the rank score is that it requires ranking both $\text{Accuracy}^{(w)}$ and BIC among all candidate models, which may be affected by the pool of candidate models. If the number of covariates in the dataset is small, we can recommend fitting all possible combinations of models. On the other hand, when the number of candidate models is very large, ranking all candidate models becomes time-consuming. To address this issue, researchers can perform model selection sequentially, called “stepwise" selection, i.e., we use the $\text{Accuracy}^{(w)}$ to select the best imputation model, and then we select PS model through BIC given that imputation model. However, such a sequential strategy is not ideal in some cases because the selection of PS model is only conditioned on the selected imputation model from the first step, which may not lead to the overall best models with the smallest RMSE. On the other side, the rank score is similar to the “best subset" strategy because both accuracy and BIC contribute equally to the rank score, ensuring a balance between the rankings of imputation and PS models.

Further research can expand the rank score into multiple areas such as MAR on both exposure and outcome. Besides that, one can also consider other missing mechanisms such as MNAR. It may be interesting for researchers to incorporate missing exposure into other complex settings, such as longitudinal or survival data in future studies.

6 Acknowledgements

We appreciate IntegraSUS for publicly sharing the COVID-19 data in Brazil. We also appreciate Wei Liang for his great suggestions on the Appendix.

7 Appendix

A.7.1 Proof of Weighted Accuracy

The key idea of weighted accuracy is to create a pseudo sample after weighting on those observed individuals so that the weighted accuracy on the observed dataset can represent the original accuracy on the whole dataset. Based on MAR assumption, if the missingness model is correct, under the true value of $w_{i}$, we can prove that,

$$\text{$\text{Accuracy}^{(w)}$}(A^{\text{imp}})\xrightarrow[]{\text{p}} E\left[\frac{1-R_{i}}{1-w_{i}}\mathbbm{1}(A_{i}=A^{\text{imp}})\right]\frac{E(\delta_{i})}{E(\delta_{i})}$$ (A.7)

$$=E\left\{\frac{E[1-R_{i}|\bm{X},Y]}{1-w_{i}}E[\mathbbm{1}(A_{i}=A^{\text{imp}})|\bm{X},Y]\right\}$$

$$=E\left\{E[\mathbbm{1}(A_{i}=A^{\text{imp}})|\bm{X},Y]\right\}=P(A_{i}=A^{\text{imp}}).$$

Here the first equality holds because $\delta_{i}$ is an indicator of whether the individual is randomly selected into the testing data or not, which is free of $(\bm{X},Y)$. Its expectation is an arbitrarily chosen ratio, which can also be canceled. The second equality holds because after taking double expectation with the condition of $(\bm{X},Y)$, we know $A_{i}\bot R_{i}|\bm{X},Y$ based on the MAR assumption. The third equality holds because the full missingness model is correctly specified, so we know $E(R_{i}|\bm{X},Y)=w_{i}$. After taking the double expectation again, we know $\text{Accuracy}^{(w)}$ will be converged to the true accuracy on the original data as $n\to\infty$.

Notice that Equation A.7 also describes why we cannot only use the naive way to estimate the accuracy on the observed data, written as: $\text{Naive-Accuracy}=\sum_{i=1}^{n}(1-R_{i})\mathbbm{1}(A_{i}=A^{\text{imp}})$. In such a way, the missing indicator will largely affect the naive estimator, which cannot be canceled without the weights $w_{i}$. Therefore, to properly evaluate the performance of the imputation model, we need to adopt the inverse weights of the missingness on the estimated accuracy of the observed data.

In the simulation studies, since we know the true missing values, we can also impute all missing values on the whole data via MICE and calculate the standard accuracy, called the “Benchmark-Accuracy" on the original data. In such a way, the bias between $\text{Accuracy}^{(w)}$ and Benchmark-Accuracy is relatively small, which supports the statement that the weighted accuracy will converge to the true accuracy as $n\to\infty$. However, we will not know the true accuracy in the application due to the unknown missing data, so $\text{Accuracy}^{(w)}$ will be a useful tool for researchers to evaluate the selected imputation model.

Additionally, IPW method for accuracy can also be easily extended into other criteria, such as sensitivity or specificity. We only need to modify the indicator function in Equation A.7 to represent the corresponding criterion of interest. Then, following the same steps, we can obtain the weighted criterion, which can evaluate the performance of the imputation model only based on the observed data.

A.7.2 Estimate ACE for Continuous Outcome

In this subsection, we provide the selection and estimated results when we intend to estimate the average causal effect (ACE), defined as $\tau=\tau_{1}-\tau_{0}$ for continuous outcome $Y$. The other setup for the simulation studies are same as before.

Similarly, Figure A.5-A.6 show that the full imputation model plus the outcome-related PS model leads to the smallest RMSE among other candidate models. In contrast, the outcome-related imputation model plus the exposure-related PS model leads to the largest RMSE, which is the worst model. Additionally, the correlation of the rank score with RMSE is higher than 0.8 and the rank score is the smallest value for the full imputation plus the outcome-related PS model. Even though the correlation between weighted accuracy and RMSE is also higher than 0.8 in this case, we cannot select the best PS model if we only use weighted accuracy. Compared with other traditional criteria, the rank score combines the performance of both models, so it can successfully pick up the best models together.

Figure A.5: Boxplot of estimated values for IPW estimators in the simulation studies for the continuous outcome to estimate ACE when $n=500$. The x-axis refers to which imputation and PS models are used. Here, Exp, Out, Res, and Covs, refer to the exposure-related, outcome-related, response-included, and covariates-included models, respectively. The red line is a true causal effect.

Figure A.6: Boxplot of estimated values for DR estimators in the simulation studies for the continuous outcome to estimate ACE when $n=500$. The x-axis refers to which imputation and PS models are used. Here, Exp, Out, Res, and Covs, refer to the exposure-related, outcome-related, response-included, and covariates-included models, respectively. The red line is the true causal effect.

Table A.12: Performance of different models using IPW method when $n=500$ for continuous outcome to estimate ACE

Imp, PS Models	Bias	Bias Rate	ESE	RMSE
Exp, naive	-0.754	-37.708	0.186	0.777
Exp, exp	-0.895	-44.766	0.203	0.918
Exp, out	-0.757	-37.865	0.088	0.762
Exp, covs	-0.898	-44.885	0.122	0.906
Out, naive	-1.163	-58.158	0.366	1.219
Out, exp	-1.277	-63.846	0.366	1.328
Out, out	-0.910	-45.479	0.082	0.913
Out, covs	-1.005	-50.270	0.090	1.009
Covs, naive	-0.752	-37.589	0.358	0.833
Covs, exp	-0.889	-44.452	0.419	0.983
Covs, out	-0.758	-37.915	0.088	0.763
Covs, covs	-0.897	-44.874	0.123	0.906
Res, naive	0.428	21.390	0.299	0.522
Res, exp	0.548	27.392	0.362	0.656
Res, out	-0.466	-23.292	0.113	0.479
Res, covs	-0.538	-26.893	0.195	0.572
Full, naive	-0.001	-0.033	0.298	0.298
Full, exp	0.000	0.012	0.361	0.361
Full, out	-0.018	-0.922	0.101	0.103
Full, covs	-0.019	-0.935	0.141	0.143

• Here, Exp, Out, Res, and Covs, refer to the exposure-related, outcome-related, response-included, and covariates-included models, respectively.

Table A.13: Performance of different models using DR method when $n=500$ for continuous outcome to estimate ACE

Imp, PS Models	Bias	Bias Rate	ESE	RMSE
Exp, naive	-0.754	-37.718	0.186	0.777
Exp, exp	-0.897	-44.845	0.204	0.920
Exp, out	-0.757	-37.831	0.088	0.762
Exp, covs	-0.901	-45.051	0.097	0.906
Out, naive	-1.163	-58.174	0.366	1.220
Out, exp	-1.271	-63.529	0.367	1.323
Out, out	-0.903	-45.130	0.082	0.906
Out, covs	-1.043	-52.127	0.087	1.046
Covs, naive	-0.752	-37.599	0.358	0.833
Covs, exp	-0.890	-44.476	0.422	0.984
Covs, out	-0.756	-37.817	0.088	0.761
Covs, covs	-0.901	-45.043	0.098	0.906
Res, naive	0.428	21.389	0.299	0.522
Res, exp	0.551	27.542	0.363	0.660
Res, out	-0.492	-24.621	0.103	0.503
Res, covs	-0.569	-28.453	0.133	0.584
Full, naive	-0.001	-0.036	0.297	0.297
Full, exp	0.003	0.130	0.363	0.363
Full, out	-0.018	-0.878	0.101	0.102
Full, covs	-0.019	-0.932	0.117	0.119

• Here, Exp, Out, Res, and Covs, refer to the exposure-related, outcome-related, response-included, and covariates-included models, respectively.

Table A.14: Criterion values for different models using IPW method when $n=500$ for continuous outcome to estimate ACE

Imp, PS Models	$\text{Accuracy}^{(w)}$	Out BIC	ASMD	KS	ABIC	Rank Score	Rank(RMSE)
Exp, naive	0.592	2470.074	0.161	0.130	1.7032	13.364	11
Exp, exp	0.592	2473.226	0.022	0.130	1.7062	13.946	16
Exp, out	0.592	1680.469	0.149	0.130	0.9350	8.879	9
Exp, covs	0.592	1676.146	0.014	0.130	0.9308	8.284	14
Out, naive	0.512	2481.170	0.102	0.148	1.9921	18.033	19
Out, exp	0.512	2479.270	0.031	0.148	1.9903	17.950	20
Out, out	0.512	1707.195	0.075	0.148	1.2392	13.368	15
Out, covs	0.512	1684.670	0.008	0.148	1.2172	12.552	18
Covs, naive	0.591	2467.830	0.170	0.136	1.7068	13.376	12
Covs, exp	0.591	2470.496	0.033	0.136	1.7094	13.835	17
Covs, out	0.591	1680.727	0.148	0.136	0.9410	9.018	10
Covs, covs	0.591	1676.198	0.014	0.136	0.9366	8.403	13
Res, naive	0.639	2396.528	0.219	0.153	1.4664	8.541	6
Res, exp	0.639	2400.459	0.077	0.153	1.4703	8.949	8
Res, out	0.639	1627.598	0.156	0.152	0.7184	4.648	5
Res, covs	0.639	1628.720	0.018	0.152	0.7195	4.836	7
Full, naive	0.794	2429.053	0.167	0.130	0.9480	7.019	3
Full, exp	0.794	2433.731	0.028	0.130	0.9526	7.489	4
Full, out	0.794	1454.705	0.149	0.130	0.0004	1.024	1
Full, covs	0.794	1459.247	0.014	0.130	0.0049	1.488	2

• $\text{Accuracy}^{(w)}$: the weighted accuracy on the observed data. Out BIC is BIC of the outcome model. ASMD: absolute standardized mean difference. KS: Kolmogorov-Smirnov distance for $(X_{1},X_{2},X_{3})$. ABIC is the rescaled form of Accuracy and BIC. Rank Score: average ranks of $1-\text{Accuracy}^{(w)}$ and BIC.

Table A.15: Criterion values for different models using DR method when $n=500$ for continuous outcome to estimate ACE

Imp, PS Models	$\text{Accuracy}^{(w)}$	Out BIC	ASMD	KS	ABIC	Rank Score	Rank(RMSE)
Exp, naive	0.592	2470.074	0.161	0.130	1.7032	13.364	11
Exp, exp	0.592	2473.226	0.022	0.130	1.7062	13.946	16
Exp, out	0.592	1680.469	0.149	0.130	0.9350	8.879	10
Exp, covs	0.592	1676.146	0.014	0.130	0.9308	8.284	14
Out, naive	0.512	2481.170	0.102	0.148	1.9921	18.033	19
Out, exp	0.512	2479.270	0.031	0.148	1.9903	17.950	20
Out, out	0.512	1707.195	0.075	0.148	1.2392	13.368	15
Out, covs	0.512	1684.670	0.008	0.148	1.2172	12.552	18
Covs, naive	0.591	2467.830	0.170	0.136	1.7068	13.376	12
Covs, exp	0.591	2470.496	0.033	0.136	1.7094	13.835	17
Covs, out	0.591	1680.727	0.148	0.136	0.9410	9.018	9
Covs, covs	0.591	1676.198	0.014	0.136	0.9366	8.403	13
Res, naive	0.639	2396.528	0.219	0.153	1.4664	8.541	6
Res, exp	0.639	2400.459	0.077	0.153	1.4703	8.949	8
Res, out	0.639	1627.598	0.156	0.152	0.7184	4.648	5
Res, covs	0.639	1628.720	0.018	0.152	0.7195	4.836	7
Full, naive	0.794	2429.053	0.167	0.130	0.9480	7.019	3
Full, exp	0.794	2433.731	0.028	0.130	0.9526	7.489	4
Full, out	0.794	1454.705	0.149	0.130	0.0004	1.024	1
Full, covs	0.794	1459.247	0.014	0.130	0.0049	1.488	2

• $\text{Accuracy}^{(w)}$: the weighted accuracy on the observed data. Out BIC is BIC of the outcome model. ASMD: absolute standardized mean difference. KS: Kolmogorov-Smirnov distance for $(X_{1},X_{2},X_{3})$. ABIC is the rescaled form of Accuracy and BIC. Rank Score: average ranks of $1-\text{Accuracy}^{(w)}$ and BIC.

Table A.16: Spearman correlation between criteria and RMSE when $n=500$ for continuous outcome to estimate ACE

Criterion	IPW Method	DR Method
1-$\text{Accuracy}^{(w)}$	0.867	0.867
ASMD	-0.339	-0.339
KS	0.203	0.201
Out BIC	0.583	0.583
ABIC	0.730	0.730
Rank Score	0.854	0.854

• Here, ASMD: absolute standardized mean difference for $(X_{1},X_{2},X_{3})$. KS: Kolmogorov-Smirnov distance for $(X_{1},X_{2},X_{3})$. ABIC is the rescaled form of Accuracy and BIC. Rank Score: averaged ranks of $1-\text{Accuracy}^{(w)}$ and BIC.