Incorporation of missing indicator with multiple imputation in propensity score analysis with partially observed covariates: A simulation study

Abstract

One of the primary challenges encountered in propensity score (PS) weighting is the presence of observations with missing covariates. In such cases, several potential solutions based on multiple imputation have been proposed. The most prevalent of these is the MI_te method, which combines treatment effect estimates derived from imputed datasets. A limited number of PS studies have incorporated the MI_te method with the missing indicator method; however, these studies only incorporated the missing indicator into the PS model. The aim of this simulation study is to propose two novel methods that incorporate the missing indicator approach with the MI_te. This incorporation either entails including the missing indicator into the outcome model (MIMI_o) or, alternatively, into both the outcome and PS model (MIMI_pso). The construction of the simulation scenarios was predicated on three elements: the mechanism of missing data, the type of treatment effect, and the presence of unmeasured confounding. In the presence of unmeasured confounding, the MIMI_pso method was the most effective method under the MAR mechanism. In the context of the MNAR mechanism, the method that exhibited the lowest bias was MIMI_o for homogeneous treatment effect and MIMI_pso for heterogeneous treatment effect. The MI_te method exhibited the highest levels of bias and variation. In view of the difficulties involved in identifying the mechanism of missing data, the variability in treatment effects across subgroups and the potential for unmeasured confounding variables in practice, researchers are encouraged to utilize the MIMI_pso method.

Keywords

Propensity score weighting missing data missing indicator multiple imputation unmeasured confounding

1 Introduction

Propensity score (PS) analysis can be employed to obtain unbiased estimates of the treatment effect in nonrandomized studies that contain confounding variables. PS, e(x_i), is defined as the conditional probability of a patient being assigned a treatment given observed covariates, a one-dimensional summary of multidimensional covariates.¹

Inverse probability of treatment weighting (IPTW) is a commonly utilized PS-based balancing strategy, and this paper focuses on PS-IPTW. In IPTW, each patient is typically assigned an unstabilized weight $1 / e (x_{i})$ for the treatment group and $1 / (1 - e (x_{i}))$ for the control group. This may result in the assignment of extreme weights. However, stabilized weights, defined as $Pr (Z = 1) / e (x_{i})$ for patients in the treatment group and $Pr (Z = 0) / (1 - e (x_{i}))$ for the control group (where $Pr (Z = 1)$ and $Pr (Z = 0)$ represent the marginal probabilities of treatment and control in the total sample), have also been proposed as a means of addressing extreme weights. This is particularly relevant for treated patients with low PSs or untreated patients with high PSs.² The use of stabilized weights does not affect the bias of treatment effect estimation; however, it does provide a significant statistical benefit in the form of increased precision in treatment effect estimates, achieved by reducing the variability in the distribution of weights.^3,4

In practice, one of the primary challenges in calculating the PS is the presence of missing values in the covariates. In such cases, the complete case (CC) analysis, the missing indicator method, and multiple imputation (MI) are viable approaches that can be employed. The most straightforward of these methods, CC analysis, is based on estimating the PS only for patients who have observed values for all covariates. CC analysis generally produces biased estimates unless the missing data mechanism is missing completely at random (MCAR), and is inefficient in terms of power under all missing data mechanisms.⁵ Another straightforward method, the missing indicator method (A missing indicator variable, which takes the value 0 if the corresponding value is observed and 1 otherwise, is included in PS model), is seldom employed due to its tendency to introduce bias.^6,7 However, Blake et al. suggest that this method provides unbiased estimates in PS analysis under the following assumptions: no unmeasured confounding; the conditional independence of the covariate values of patients with missing data from treatment and outcome; and the correct specification of the outcome model.⁸ CC analysis and the missing indicator method perform poorly in terms of overall balance, and methods that explicitly account for missing data perform better than these methods, according to Bottigliengo et al.⁹ The final method, MI, is a widely utilized approach for addressing missing data in medical research. The fundamental premise of MI is the iterative estimation of missing values through the sampling of the posterior predictive distribution, conditioned on the observed data. A substantial body of evidence exists demonstrating the superiority of the MI over alternative approaches for handling missing data.^10–12

1.1 Motivation

A number of potential solutions for combining post-MI estimates in PS analysis have been put forth in the literature. These methods entail the combination of treatment effect estimates (referred to as “MI_te”), PSs (referred to as “MI_ps”), or parameter estimates within the PS model (referred to as “MI_par”) with the objective of estimating the treatment effect.^13,14 The majority of authors in the literature, with the exception of Leite et al., concur that the MI_te approach is the optimal method under PS weighting.^13–18 Qu and Lipkovic and Seaman and White employed a methodology termed “multiple imputation missingness pattern” to assess the efficacy of the MI_te approach. This involved incorporating the pattern of missing data into the PS model.^19,20 In contrast, Choi et al. also evaluated the performance of MI_te by including the missing indicator into the PS model and concluded that the optimal method selection is contingent upon the missing data structure and the presence of effect modification.²¹ One area of the literature that requires further development is the refinement of the MI_te method, specifically the addition of the missing indicator to the PS model but not the outcome model. The explicit modelling of missing data, including the missing indicator within the outcome model and both the PS and the outcome model, may prove beneficial for the estimation of treatment effects. In light of these considerations, the aim of this simulation study is to propose two novel methods, namely MI_te missing indicator in outcome model (MIMI_o) and MI_te missing indicator in PS and outcome model (MIMI_pso), which seek to incorporate the missing indicator approach with MI_te by including the missing indicator into the outcome model or, alternatively, into both the outcome and PS models.

2 Design of the simulation study

A data generation mechanism comparable to that described by Choi et al. was adopted along with a series of supplementary scenarios, and previously employed approaches were compared with the novel methods.²¹ All simulations were conducted using the R program, version 4.2.1.²²

2.1 Data generation mechanisms

Two scenarios were considered: no unmeasured confounding and unmeasured confounding. In the simulation scenarios, two continuous covariates ( $X_{1}$ and $X_{2}$ ) were generated for $n = 1000$ subjects i. $X_{1}$ was fully observed whereas $X_{2}$ was partially observed.

2.1.1 Scenario 1: no unmeasured confounding

Covariates: $X_{1}$ follows a normal distribution of mean 0 and standard deviation (SD) of 1. $X_{2}$ depends on $X_{1}$ .

X_{2 i} = 0.5 X_{1 i} + ε_{i}, ε_{i} \sim N (0, 0.75)

In this way the SD of $X_{2}$ is 1 and the correlation between $X_{1}$ and $X_{2}$ is equal to 0.5.

Treatment assignment: The treatment T was generated from the binomial distribution depend on $X_{1}$ and $X_{2}$ according to following model:

logit (P (T_{i} = 1 | X_{1 i}, X_{2 i})) = - 0.8 + 0.5 X_{1 i} + 0.5 X_{2 i}

In this model, 330 subjects, or approximately 33% of the subjects, received treatment.

Outcome: Depending on $X_{1}$ and $X_{2}$ , and different treatment effect conditions, three models were used to generate the continuous outcome $Y$ :

(1)
Homogeneous null treatment effect: $Y_{i} = X_{1 i} + X_{2 i} + ϵ_{i}, ϵ_{i} \sim N (0, 1)$
(2)
Homogeneous non-null treatment effect: $Y_i = X_{1 i} + X_{2 i} + 2 T_{i} + ϵ_{i}, ϵ_{i} \sim N (0, 1)$
(3)
Heterogeneous treatment effect: $Y_i = X_{1 i} + X_{2 i} + T_{i} X_{2 i} + ϵ_{i}, ϵ_{i} \sim N (0, 1)$

The homogeneous treatment effect assumes that the treatment had a constant effect on the outcome. In the case of heterogeneous treatment effect, adding effect modification $(T_{i} X_{2 i})$ to the model ensures that the treatment effect varies across subgroups.
2.1.2 Scenario 2: unmeasured confounding

Covariates: $X_{1}$ follows a normal distribution of mean 0 and SD of 1. Unmeasured confounder U also follows a normal distribution with mean 0 and SD 1, but is independent of $X_{1}$ . $X_{2}$ depends both on $X_{1}$ and U.

X_{2 i} = 0.5 X_{1 i} + 0.5 U_{i} + ε_{i}, ε_{i} \sim N (0, 0.5)

Treatment assignment: The treatment T was generated from the binomial distribution depend on $X_{1}$ , $X_{2}$ and $U$ according to following model:

logit (P (T_{i} = 1 | X_{1 i}, X_{2 i}, U_{i})) = - 0.85 + 0.5 X_{1 i} + 0.5 X_{2 i} + 0.5 U_{i}

In this model, 330 subjects, or approximately 33% of the subjects, received treatment.

Outcome: Depending on $X_{1}$ , $X_{2}$ and U, different treatment effect conditions, three models were used to generate the continuous outcome $Y$ :

(1)
Homogeneous null treatment effect: $Y_{i} = X_{1 i} + X_{2 i} + U_{i} + ϵ_{i}, ϵ_{i} \sim N (0, 1)$
(2)
Homogeneous non-null treatment effect: $Y_i = X_{1 i} + X_{2 i} + U_{i} + 2 T_{i} + ϵ_{i}, ϵ_{i} \sim N (0, 1)$
(3)
Heterogeneous treatment effect: $Y_i = X_{1 i} + X_{2 i} + U_{i} + T_{i} X_{2 i} + ϵ_{i}, ϵ_{i} \sim N (0, 1)$

In both scenarios, missing values were generated for 50% of the $X_{2}$ observations depending on the missing data mechanisms: −
MCAR: A random 50% of the observations in $X_{2}$ are set to be missing.
−
MAR: Let R be the binary missing indicator for $X_{2}$ .The probability of $X_{2}$ being missing is defined as $logit (P (R_{i} = 1)) = X_{1 i}$ . In this case, probability that $X_{2}$ is missing increases with the value of $X_{1}$ .
−
MNAR: Scenario 1: The probability of $X_{2}$ being missing is defined as
$logit (P (R_{i} = 1)) = X_{2 i}$

Scenario 2: The missingness in $X_{2}$ is generated if U > 0.
2.2. Analysis of the simulated datasets

Firstly, it should be noted that in all simulated datasets, the covariates were heterogeneous across treatment groups. However, further details on this are not presented in this study, as it was beyond the scope of the research. The PSs were obtained through logistic regression, after which the stabilized weights were calculated using the PSs. The average treatment effect was estimated as the coefficient of the treatment in the weighted regression model, where the outcome was predicted solely by the treatment. Following the generation of missing values in $X_{2}$ , the MI was implemented using a fully conditional specification procedure with the R package mice, resulting in 10 imputed datasets.²³ The outcome was included in the imputation model since de Vries and Groenewold highlighted its importance.²⁴ Furthermore, given that the imputation model that excluded interactions yielded biased estimates and low coverage, two distinct approaches were employed in constructing the imputation model: Only main effects $(X_{1}, T, Y)$ , Main effects and interaction terms $(X_{1} T, X_{1} Y, T Y)$ .²⁵

2.2.1 Methods applied after MI

The following four methods, which are to be employed after imputing missing data, are summarized below. In order to apply the methods outlined below, a missing indicator variable, denoted as R, was defined with a value of 1 for each subject with missing data, and 0 otherwise.

(1)
MI_te method: For each imputed dataset, (a) PSs were calculated by logistic regression, (b) Stabilized weights were obtained, (c) Treatment effects were estimated by weighted linear regression. Finally, using Rubin's rule, treatment effect estimates were combined into a single treatment effect across the estimated data sets.¹⁴
(2)
MI_te missing indicator-PS model (MIMI_ps): The first step of the MI_te method was modified by adding R in the PS model along with the covariates. The remainder of the procedure remains unchanged.
(3)
MI_te missing indicator-outcome model (MIMI_o): The first two steps of the MI_te method were applied. In the next third step, the outcome model, which estimates the treatment effect, was modified by adding R.
(4)
MI_te missing indicator-PS and outcome models (MIMI_pso): The PS and outcome models in steps 1 and 3 of the MI_te method were modified by adding R.

2.3. Measures of performance

In all scenarios, 1000 datasets were generated through simulation. In order to facilitate comparisons between the methods, the bias and SD of the treatment effect were taken into account as performance measures.

3 Results

The performance measures of the MI_te, MIMI_ps, MIMI_o and MIMI_pso methods for simulation scenarios 1 and 2 were presented in Tables 1, 2 and 3 and Figures 1 and 2, respectively, for the homogeneous null, homogeneous non-null and heterogeneous treatment effects under different missing data mechanisms.

Table 1.
Treatment effect bias and SD in the homogeneous (Null) treatment effect.

Scenario 1: No unmeasured confounding Scenario 2: Unmeasured confounding

Only main effects Main effects and interaction terms Only main effects Main effects and interaction terms

Imputation model type Bias SD Bias SD Bias SD Bias SD

Baseline 1.233 0.127 1.619 0.160 1.930 0.144 1.936 0.146

After PS 0.085 0.069 0.107 0.087 0.267 0.130 0.264 0.137

MCAR

MI_te 0.096 0.065 0.096 0.062 0.278 0.130 0.273 0.138

MIMI_ps 0.095 0.062 0.097 0.062 0.287 0.130 0.272 0.131

MIMI_o 0.103 0.069 0.109 0.069 0.285 0.142 0.291 0.143

MIMI_pso 0.098 0.063 0.095 0.061 0.279 0.131 0.283 0.133

MAR

MI_te 0.102 0.069 0.111 0.071 0.275 0.135 0.284 0.141

MIMI_ps 0.099 0.065 0.108 0.066 0.281 0.138 0.282 0.141

MIMI_o 0.104 0.068 0.105 0.073 0.306 0.146 0.304 0.145

MIMI_pso 0.096 0.063 0.105 0.064 0.275 0.130 0.283 0.138

MNAR

MI_te 0.103 0.069 0.111 0.070 0.294 0.137 0.297 0.145

MIMI_ps 0.096 0.063 0.104 0.066 0.283 0.134 0.279 0.136

MIMI_o 0.100 0.069 0.102 0.073 0.250 0.132 0.252 0.131

MIMI_pso 0.095 0.065 0.101 0.064 0.278 0.130 0.279 0.132

	Scenario 1: No unmeasured confounding	Scenario 2: Unmeasured confounding
Baseline	1.233	0.127	1.619	0.160	1.930	0.144	1.936	0.146
After PS	0.085	0.069	0.107	0.087	0.267	0.130	0.264	0.137
MCAR
MI_te	0.096	0.065	0.096	0.062	0.278	0.130	0.273	0.138
MIMI_ps	0.095	0.062	0.097	0.062	0.287	0.130	0.272	0.131
MIMI_o	0.103	0.069	0.109	0.069	0.285	0.142	0.291	0.143
MIMI_pso	0.098	0.063	0.095	0.061	0.279	0.131	0.283	0.133
MAR
MI_te	0.102	0.069	0.111	0.071	0.275	0.135	0.284	0.141
MIMI_ps	0.099	0.065	0.108	0.066	0.281	0.138	0.282	0.141
MIMI_o	0.104	0.068	0.105	0.073	0.306	0.146	0.304	0.145
MIMI_pso	0.096	0.063	0.105	0.064	0.275	0.130	0.283	0.138
MNAR
MI_te	0.103	0.069	0.111	0.070	0.294	0.137	0.297	0.145
MIMI_ps	0.096	0.063	0.104	0.066	0.283	0.134	0.279	0.136
MIMI_o	0.100	0.069	0.102	0.073	0.250	0.132	0.252	0.131
MIMI_pso	0.095	0.065	0.101	0.064	0.278	0.130	0.279	0.132

SD: standard deviation, MIMI_ps: MI_te missing indicator-PS, MIMI_o: MI_te missing indicator-outcome, MIMI_pso: MI_te missing indicator-PS and outcome.

Table 2.

Treatment effect bias and SD in the homogeneous (Non-Null) treatment effect.

	Scenario 1: No unmeasured confounding				Scenario 2: Unmeasured confounding
	Only main effects		Main effects and interaction terms		Only main effects		Main effects and interaction terms
Imputation model type	Bias	SD	Bias	SD	Bias	SD	Bias	SD
Baseline	1.234	0.119	1.237	0.122	1.932	0.145	1.933	0.144
After PS	0.084	0.105	0.084	0.109	0.265	0.165	0.264	0.154
MCAR
MI_te	0.089	0.112	0.092	0.116	0.273	0.169	0.267	0.166
MIMI_ps	0.088	0.112	0.095	0.121	0.271	0.158	0.271	0.168
MIMI_o	0.095	0.122	0.101	0.129	0.268	0.166	0.279	0.188
MIMI_pso	0.091	0.116	0.090	0.114	0.274	0.157	0.274	0.179
MAR
MI_te	0.089	0.112	0.092	0.116	0.276	0.169	0.270	0.165
MIMI_ps	0.093	0.119	0.099	0.128	0.272	0.165	0.283	0.182
MIMI_o	0.093	0.116	0.101	0.129	0.287	0.180	0.295	0.183
MIMI_pso	0.089	0.115	0.094	0.121	0.267	0.156	0.275	0.186
MNAR
MI_te	0.090	0.113	0.102	0.128	0.275	0.188	0.276	0.191
MIMI_ps	0.090	0.115	0.099	0.125	0.277	0.167	0.278	0.178
MIMI_o	0.089	0.113	0.099	0.126	0.234	0.172	0.239	0.176
MIMI_pso	0.089	0.113	0.090	0.113	0.269	0.157	0.273	0.177

SD: standard deviation, MIMI_ps: MI_te missing indicator-PS, MIMI_o: MI_te missing indicator-outcome, MIMI_pso: MI_te missing indicator-PS and outcome.

Table 3.

Treatment effect bias and SD in the heterogeneous treatment effect.

	Scenario 1: No unmeasured confounding				Scenario 2: Unmeasured confounding
	Only main effects		Main effects and interaction terms		Only main effects		Main effects and interaction terms
Imputation model type	Bias	SD	Bias	SD	Bias	SD	Bias	SD
Baseline	1.622	0.146	1.614	0.154	2.409	0.179	2.423	0.175
After PS	0.104	0.085	0.109	0.087	0.292	0.153	0.291	0.158
MCAR
MI_te	0.130	0.084	0.135	0.088	0.306	0.156	0.307	0.154
MIMI_ps	0.135	0.090	0.138	0.092	0.269	0.136	0.309	0.156
MIMI_o	0.152	0.110	0.150	0.089	0.288	0.159	0.332	0.169
MIMI_pso	0.130	0.083	0.130	0.079	0.273	0.144	0.312	0.155
MAR
MI_te	0.129	0.082	0.135	0.088	0.309	0.156	0.309	0.154
MIMI_ps	0.155	0.113	0.153	0.095	0.234	0.124	0.328	0.170
MIMI_o	0.193	0.138	0.186	0.132	0.336	0.183	0.294	0.195
MIMI_pso	0.148	0.100	0.140	0.082	0.229	0.134	0.328	0.160
MNAR
MI_te	0.153	0.106	0.148	0.088	0.300	0.149	0.336	0.168
MIMI_ps	0.203	0.149	0.192	0.135	0.187	0.119	0.256	0.181
MIMI_o	0.209	0.141	0.189	0.131	0.202	0.139	0.243	0.200
MIMI_pso	0.191	0.128	0.173	0.125	0.182	0.138	0.219	0.163

SD: standard deviation, MIMI_ps: MI_te missing indicator-PS, MIMI_o: MI_te missing indicator-outcome, MIMI_pso: MI_te missing indicator-PS and outcome.

3.1 Scenario 1: no unmeasured confounding

3.1.1 Homogeneous treatment effect

Among the models for all methods, the SD was observed to be lower when the treatment effect was null, and the bias was found to be slightly lower when the treatment effect was non-null. In the context of a null treatment effect, the MIMI_ps method in MCAR, the MIMI_pso method in MAR, and the MIMI_ps and MIMI_pso methods in MNAR were identified as the methods with the smallest bias and variation. In the context of a non-null treatment effect, the methods exhibiting the smallest bias and variation were MI_te and MIMI_ps in the MCAR scenario. In contrast, the performance of the four methods in the MAR and MNAR scenarios was found to be almost similar. The efficacy of the methods was diminished when interaction terms were included into the imputation model (Tables 1 and 2 and Figure 1).

Figure 1.

Bias of methods applied after MI in SCENARIO 1: no unmeasured confounding.

Figure 2.

Bias of methods applied after MI in SCENARIO 2: unmeasured confounding.

3.1.2 Heterogeneous treatment effect

The most optimal results were achieved through the use of the MI_te and MIMI_pso methods in the context of MCAR, the MI_te method in the context of MAR, and the MI_te method that includes an imputation model with interaction terms in the context of MNAR (Table 3, Figure 1).

3.2 Scenario 2: unmeasured confounding

3.2.1 Homogeneous treatment effect

The methods with the lowest bias and SD when the treatment effect was null were the MI_te and MIMI_pso in MCAR and MAR, and the MIMI_o method in MNAR. When the treatment effect was non-null, the method with the lowest bias was MIMI_o in MCAR and MNAR and MIMI_pso in MAR (Tables 1 and 2 and Figure 2).

3.2.2 Heterogeneous treatment effect

In MCAR, the MIMI_ps method exhibited the lowest level of bias. In contrast, in MAR and MNAR, the MIMI_pso method demonstrated the lowest level of bias. The inclusion of interaction terms into the imputation model yielded no discernible impact, akin to the MNAR scenario in Scenario 1 (Table 3, Figure 2).

In this scenario, the method with the highest degree of bias and variation in MNAR mechanism was the MI_te method.

The incorporation of a missing indicator into the PS and outcome model within the MAR mechanism resulted in a favorable and positive impact on bias reduction, irrespective of whether the treatment effect exhibited homogeneity or heterogeneity. This impact is also applicable in the case of MNAR, where the missing indicator was included into the outcome model for homogeneous treatment effects and into the PS and outcome models for heterogeneous treatment effects.

The impact of a 25% missing data is also assessed in this study, with the relevant outcomes presented in the Supplementary material. It was observed that the most efficacious methods for 25% and 50% missing data were not radically different and, in fact, were largely parallel.

4. Discussion

This study investigates the impact of incorporating the missing indicator into the MI_te method in different ways on the treatment effect estimates in the presence of missing data in the PS-IPTW analysis. In this context, the research question guiding the study was whether the missing indicator should be incorporated into the MI_te method with the PS model, the outcome model, or both the PS and the outcome model. This evaluation was conducted through a simulation scenario. To ensure comparability of results, we employed the data generation scenario proposed by Choi et al., but pursued a distinct analytical trajectory in examining the simulated datasets.²¹ The initial discrepancy pertained to the utilization of standardized weights, a methodology unfamiliar to us in the context of PS-IPTW analysis when dealing with missing data. Secondly, in all scenarios, interaction terms were included into the imputation model, irrespective of whether heterogeneous treatment effects were also considered. In this context, as in the study of Choi et al., the MI_te and MIMI_ps methods were used in PS-IPTW analysis to handle missing data. In light of these considerations, the scenarios in which our study demonstrated superior performance were identified through a comparative analysis of the two studies. These findings are summarized below: In the context of a homogeneous null treatment effect, the SD was observed to be lower. Similarly, in the presence of unmeasured confounding and a homogeneous null treatment effect, the overall performance in MCAR and the SD in MNAR were superior. Furthermore, in the heterogeneous treatment effect, the SD and also the bias in MNAR were lower. It is noteworthy that the SDs of the treatment effect estimates in our study were frequently lower than those reported by Choi et al. (2019). It can be posited that the primary reason for this discrepancy is the utilization of stabilized weights, which differs from the approach employed by Choi et al. (2019). The application of stabilized weights serves to reduce the variability inherent in the weight distribution. However, even with disparate weightings, it was observed that the confidence intervals for the estimated treatment effects exhibited overlap, indicating that the fluctuations in bias and SD are mutually compensating.

Choi et al. evaluated the CC analysis and missing indicator methods, as well as the MI_te and MIMI_ps methods.²¹ They concluded that the optimal approach to handling missing data depends on the specific missing data mechanism and the presence of effect modification. The findings of their study indicate that the missing indicator method and CC analysis are not generally recommended for handling missing values. In light of the extensive criticism directed at these methods in the literature, we chose not to address the CC and missing indicator approaches in this study.^11,26,27 When these two methods are disregarded, irrespective of the missing data mechanism, MIMI_ps for the homogeneous treatment effect and MI_te, which includes interaction terms into the imputation model for the heterogeneous treatment effect, were the most effective approaches in Choi et al. However, despite the low probability of reaching a general conclusion, regardless of the missing data mechanism in this study, the overall results showed that the MI_te method was less effective than alternative approaches in a large number of scenarios. In particular, when an unmeasured confounder is present, the MI_te method should be applied by incorporating the missing indicator into either the outcome model or the PS and outcome models. This suggests that, in the presence of unmeasured confounding, the impact of missing values in the covariate should be taken into account following imputation. Furthermore, this is made evident when missing values are taken into account in both the PS and outcome models, as opposed to solely in the PS model.

The simulation study conducted by Noghrehchi et al. demonstrated that imputation models that are misspecified or overfitted can lead to bias in parameter estimates.²⁸ As demonstrated by Choi et al., the inclusion of interaction terms among covariates, treatment, and outcome in the imputation model serves to reduce bias in instances where a heterogeneous treatment effect is present. This phenomenon was limited to the MNAR mechanism in the present study, resulting in a reduction in standard deviation and the introduction of bias. The same mechanism yielded disparate results in the presence of an unmeasured confounder, namely that the efficacy of the methods was diminished by the inclusion of interaction terms into the imputation model. In other words, when the mechanism is MNAR and there are no unmeasured confounders, the addition of interaction terms to the imputation model improves the estimates of the treatment effect if it varies among subjects. The addition of interaction terms to the imputation model was found to result in an increase in bias and standard deviation, rather than providing any benefit, when the treatment effect was homogeneous. Therefore, in the case of a homogeneous treatment effect and a heterogeneous treatment effect with no unmeasured confounding, the appropriate imputation model may be one that includes only the main effects. Conversely, in the case of a heterogeneous treatment effect with unmeasured confounding, the model should include both the main effects and interaction terms. Further investigation is required to ascertain the optimal imputation model for missing data in PS-IPTW.

In response to the criticisms raised by Choi and Gelfond regarding the exclusive consideration of the homogeneous null treatment effect, Choi et al. conducted a replication of their simulations by setting the treatment effect to 1 and reported that the results remained consistent when the treatment effect was 0 and 1.^29,30 In contrast with the findings of Choi et al., our results did not yield comparable outcomes in both scenarios. This leads us to conclude that the magnitude of the treatment effect may have an impact on the selection of the appropriate methodology. Furthermore, a reduction in the variability of the estimated treatment effect was observed when the treatment effect was null. This indicates that treatment effect estimates are less consistent when there is a homogeneous non-null treatment effect. Further research is required to facilitate more detailed inferences.

In the context of routinely collected healthcare data, it is challenging to ascertain the specific missing data mechanism that is responsible for the observed missing data. A review of the literature reveals that, to the best of our knowledge, only Little's MCAR test is available for evaluating whether missing data are MCAR.³¹ Once it has been established that the missing data are not MCAR, it is not possible to determine whether the missing data are MAR or MNAR. As demonstrated in this study, the most suitable method for incorporating the missing indicator approach with MI_te is contingent upon the nature of the treatment effect, mechanism for missing data, and the existence of unmeasured confounding. In other words, irrespective of the missing data mechanism, it is evident that it is not possible to discuss the superiority of a method in terms of homogeneous or heterogeneous treatment effects. In light of these and the fact that the treatment effect often varies across subgroups in practice, the researcher may choose to commence a study with the presumption that the treatment effect is heterogeneous, that missing data are MNAR, and that there are unmeasured confounding variables. In such a framework, the most optimal method for obtaining treatment effect estimates with PS-IPTW is to utilize the MI_te approach, in which the indicator variable is included into the PS and outcome model, referred to as MIMI_pso in this study. Consequently, rather than employing the MI_te approach as the default, researchers should incorporate it with the missing indicator method.

Supplemental Material

sj-docx-1-smm-10.1177_09622802251338365 - Supplemental material for Incorporation of missing indicator with multiple imputation in propensity score analysis with partially observed covariates: A simulation study

Supplemental material, sj-docx-1-smm-10.1177_09622802251338365 for Incorporation of missing indicator with multiple imputation in propensity score analysis with partially observed covariates: A simulation study by Sevinc Puren Yucel Karakaya and Ilker Unal in Medical Research

Footnotes

Authors’ contributions

All authors contributed to the study conception and design. Data generation and analysis were performed by Sevinc Puren Yucel Karakaya and Ilker Unal. The first draft of the manuscript was written by Sevinc Puren Yucel Karakaya and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Ethics approval

The study does not require ethical approval.

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Sevinc Puren Yucel Karakaya

Supplemental material

Supplemental material for this article is available online.

References

Rosenbaum

Rubin

. The central role of the propensity score in observational studies for causal effects. Biometrika 1983; 70: 41–55.

Austin

Stuart

. Moving towards best practice when using inverse probability of treatment weighting (IPTW) using the propensity score to estimate causal treatment effects in observational studies. Stat Med 2015; 34: 3661–3679.

Brookhart

Wyss

Layton

, et al. Propensity score methods for confounding control in nonexperimental research. Circ Cardiovasc Qual Outcomes 2013; 6: 604–611.

Cole

Hernán

. Constructing inverse probability weights for marginal structural models. Am J Epidemiol 2008; 168: 656–664.

Van der Heijden

Donders

ART

Stijnen

, et al. Imputation of missing values is superior to complete case analysis and the missing-indicator method in multivariable diagnostic research: a clinical example. J Clin Epidemiol 2006; 59: 1102–1109.

Malla

Perera-Salazar

McFadden

, et al. Handling missing data in propensity score estimation in comparative effectiveness evaluations: a systematic review. J Comp Eff Res 2018; 7: 271–279.

Groenwold

White

Donders

ART

, et al. Missing covariate data in clinical research: when and when not to use the missing-indicator method for analysis. CMAJ 2012; 184: 1265–1269.

Blake

Leyrat

Mansfield

, et al. Estimating treatment effects with partially observed covariates using outcome regression with missing indicators. Biom J 2020; 62: 428–443.

Bottigliengo

Lorenzoni

Ocagli

, et al. Propensity score analysis with partially observed baseline covariates: a practical comparison of methods for handling missing data. Int J Environ Res Public Health 2021; 18: 66–94.

10.

Lee

Simpson

. Introduction to multiple imputation for dealing with missing data. Respirology 2014; 19: 162–167.

11.

Sterne

White

Carlin

, et al. Multiple imputation for missing data in epidemiological and clinical research: potential and pitfalls. Br Med J 2009; 338: b2393.

12.

Mackinnon

. The use and reporting of multiple imputation in medical research—a review. J Intern Med 2010; 268: 586–593.

13.

Mitra

Reiter

. A comparison of two methods of estimating propensity scores after multiple imputation. Stat Methods Med Res 2016; 25: 188–204.

14.

Leyrat

Seaman

White

, et al.

Propensity score analysis with partially observed covariates: how should multiple imputation be used?

Stat Methods Med Res 2019; 28: 3–19.

15.

de Vries

BBP

Groenwold

. A comparison of two approaches to implementing propensity score methods following multiple imputation. Epidemiol Biostat Public Health 2017; 14: e12630-1–e12630-21.

16.

Leite

Aydin

Cetin-Berber

. Imputation of missing covariate data prior to propensity score analysis: a tutorial and evaluation of the robustness of practical approaches. Eval Rev 2021; 45: 34–69.

17.

Nguyen

Stuart

. Multiple imputation for propensity score analysis with covariates missing at random: some clarity on within and across methods. Am J Epidemiol 2024; 193: kwae105.

18.

Hill

. Reducing bias in treatment effect estimation in observational studies suffering from missing data. ISERP Working Papers, Institute for Social and Economic Research and Policy, Columbia University, 2004.

19.

Lipkovich

. Propensity score estimation with missing values using a multiple imputation missingness pattern (MIMP) approach. Stat Med 2009; 28: 1402–1414.

20.

Seaman

White

. Inverse probability weighting with missing predictors of treatment assignment or missingness. Commun Stat Theory Methods 2014; 43: 3499–3515.

21.

Choi

Dekkers

le Cessie

. A comparison of different methods to handle missing data in the context of propensity score analysis. Eur J Epidemiol 2019; 34: 23–36.

22.

R Core Team. R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing, 2021.

23.

Buuren

Groothuis-Oudshoorn

. Mice: multivariate imputation by chained equations in R. J Stat Softw 2010; 45: 1–67.

24.

Penning de Vries

BBL

Groenwold

. Comments on propensity score matching following multiple imputation. Stat Methods Med Res 2016; 25: 3066–3068.

25.

Tilling

Williamson

Spratt

, et al. Appropriate inclusion of interactions was needed to avoid bias in multiple imputation. J Clin Epidemiol 2016; 80: 107–115.

26.

Crowe

Lipkovich

Wang

. Comparison of several imputation methods for missing baseline data in propensity scores analysis of binary outcome. Pharm Stat 2010; 9: 269–279.

27.

White

Carlin

. Bias and efficiency of multiple imputation compared with complete-case analysis for missing covariate values. Stat Med 2010; 29: 2920–2931.

28.

Noghrehchi

Stoklosa

Penev

, et al. Selecting the model for multiple imputation of missing data: just use an IC!. Stat Med 2021; 40: 2467–2497.

29.

Choi

Gelfond

. The validity of propensity score analysis using complete cases with partially observed covariates. Eur J Epidemiol 2020; 35: 87–88.

30.

Choi

Dekkers

le Cessie

. Authors’ reply: a comparison of different methods to handle missing data in the context of propensity score analysis. Eur J Epidemiol 2020; 35: 89–91.

31.

Little

. A test of missing completely at random for multivariate data with missing values. J Am Stat Assoc 1988; 83: 1198–1202.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.03 MB