Heckman-type selection models to obtain unbiased estimates with missing measures outcome: theoretical considerations and an application to missing birth weight data

Background In low-income settings, key outcomes such as biomarkers or clinical assessments are often missing for a substantial proportion of the study population. The aim of this study was to assess the extent to which Heckman-type selection models can create unbiased estimates in such settings. Methods We introduce the basic Heckman model in a first stage, and then use simulation models to compare the performance of the model to alternative approaches used in the literature for missing outcome data, including complete case analysis (CCA), multiple imputations by chained equations (MICE) and pattern imputation with delta adjustment (PIDA). Last, we use a large population-representative data set on antenatal supplementation (AS) and birth outcomes from Côte d’Ivoire to illustrate the empirical relevance of this method. Results All models performed well when data were missing at random. When missingness in the outcome data was related to unobserved determinants of the outcome, large and systematic biases were found for CCA and MICE, while Heckman-style selection models yielded unbiased estimates. Using Heckman-type selection models to correct for missingness in our empirical application, we found supplementation effect sizes that were very close to those reported in the most recent systematic review of clinical AS trials. Conclusion Missingness in health outcome can lead to substantial bias. Heckman-selection models can correct for this selection bias and yield unbiased estimates, even when the proportion of missing data is substantial.


Background
A growing literature has highlighted the often substantial differences between evidence based on efficacy trials and empirically observed associations between intervention exposure and health outcomes [1][2][3]. While this gap may to a certain extent reflect differences in programme implementation and differential adherence to treatment protocols in non-clinical settings, biases in observational studies seem also plausible. A substantial body of literature has highlighted the importance of potential confounding variables in observational studies. Slightly less attention has been given to the often substantial degree of missingness in outcome variables [4][5][6]. Missingness in the outcome variable is of particular importance in the context of clinical data in low-income settings, where accurate measures of clinical outcomes often is only available for a relatively small proportion of the population. Even though missing values can in principle be imputed using multiple imputations [4], this approach can lead to biased estimates if unobservable or unmeasured factorssuch as individual health knowledge or attitudesaffect both the outcome of interest and the likelihood of missing data [7].
To illustrate the relevance of such non-random missingness in outcome data as well as the possibility to correct for such biases using Heckman-type selection model, we focus on birth weight (BW) as primary outcome variable in this paper. Low birth weight (birth weight < 2500 g) affects 15.5% of children globally [8], and has been identified as one of the primary causes of the continued high burden of under-5 mortality in low-and middle-income countries [9]. In low income settings, birth weight is only available for women who deliver at a health centres with functioning measurement equipment as well as staff willing and able to record infant weight after birth. Given that institutional deliveries remain scarce in many settings [10], reliable data can often only be attained for a limited proportion of women. Missing outcome data will not cause systematic bias if data are missing at random (MAR). In practice, the MAR assumption will, however, not hold if unobservable traits such as preventive efforts or health knowledge predict both the likelihood to deliver at a facility (the likelihood of having birth weight data available) and the actual health outcome of interest.
The selection model introduced by Heckman [7] provides a potentially useful tool in this situation, since it allows to both test and correct for potential biases created by nonrandom missingness in outcome measures. To illustrate this, we first use Monte-Carlo simulations to assess the relative ability of different models to detect true causal effects. The specific causal effect we investigate is the effect of antenatal supplementation on birth weight. Iron and folic acid supplementation (IFAS) is widely recognized as one of the most effective interventions to address low birth weight (LBW). A meta-analysis of 11 trials revealed a reduction of the risk of LBW by 20% associated with iron supplementation or when iron supplementation was combined with folic acid (relative risk [RR] 0.80, 95% CI 0.71-0.90) [11]. The same patterns have generally not been found in observational studies [12][13][14]. We first assess the extent to which Heckman selection models, namely complete case analysis (CCA), multiple imputations by chained equations (MICE) and pattern imputation with delta adjustment (PIDA), can recover the true causal impact of interest in simulated data in a first step. In a second step, we illustrate these differences using population-representative data on antenatal supplementation (AS) and birth weight from the health and demographic surveillance site (HDSS) in Taabo, Côte d'Ivoire.

Objective and modelling background
The main objective of this paper is to compare Heckman-type selection models to alternative approaches used to deal with missing outcome data in the literature. The Heckman model includes two separate equationsone focusing on selection into the sample (outcome being observedthe sample selection equation), and the main equation linking the covariates of interest to the outcome.
The two Heckman equations for two latent responses y Ã i (the outcome) and s Ã i (the selection propensity variable) can be stated as follows [15]: Where y Ã i and s Ã i are unobserved latent continuous variables, x 0 i and z 0 i are vectors of predictor variables. In general, x is assumed to be a subset of z, which means that all factors predicting the main outcome of interest (y) also predict selection s. μ and ν are normally distributed error term, and β is the primary parameter vector of interest. Outcome variables are observed if the latent selection propensity exceeds zero, i.e.: The main idea of the Heckman model is that it seems theoretically rather likely that unobservable or unmeasured factors may affect both the outcome y and the probability of selection s; these unmeasured factors would be contained in the residuals of both equation (1) and equation (2). Given selection into the main sample, the expected value of the outcome in the main equation is given by: Given that the covariates x and v jointly determine selection into the sample, cov(x, v|s = 1) is non-zero in general, so that beta estimates will be both biased and inconsistent if μ and ν are correlated. This correlation is straightforward to estimate empirically by fitting independent models for y and s, and computing the covariance between the two residual terms. Heckman shows that this bias can be corrected by computing the expected value of v conditional on z and being in the sample, and by including this term in the main empirical model. Consistent estimators can be obtained by maximum likelihood jointly estimating the first stage with a probit model as well as the main equation of interest including the expected value of the selection equation residuals [15].

Study variables
The main outcome variables used were continuous birth weight (BW) as well as binary indicator for LBW (weight < 2500 g).
Additional variables used for the analysis of our demographic surveillance data are socioeconomic status and distance to facility. Socioeconomic status was determined using a household-based asset approach and principal component analysis (PCA) to divide households into wealth quintile (poorest, poor, medium, rich and richest) [16]. Using household and health centres geographical coordinates, we estimated the distance from mother's place of residence to the nearest health facility by means of the Statageodist package [17].

Simulations and statistical analysis
Our empirical analysis is divided into two parts. In the first part, we use Monte-Carlo simulations to illustrate the empirical performance of CCA, MICE, PIDA and Heckman with missing outcome data. Based on the empirical data used in the second part of the analysis, we assume a sample size of 10,000 births, and normally distributed birth weight with mean 3000 g, and standard deviation of 500 g. Based on the most recent systematic review, we assume that supplementation linearly increases birth weight by 50 g. We first assume that 40% of the outcome data are missing at random, and plot the estimated coefficients on supplementation based on 1000 randomly created data sets. In a second step, we assume that missing outcome data is a function of unobserved health knowledge, and that unobserved health knowledge is also predictive of birth weight. For the data generating process, we assume that health knowledge follows a standard normal distribution, and that each standard deviation (SD) increase in health knowledge increases birth weight by 100 g. We also assume the probability of delivering increases with the unobserved health knowledge variable, and decreases with household distance from the facility. We then test the various modelling approaches under this "endogenous selection" (as Heckman refers to it) scenario.
To illustrate the empirical relevance of this approach, we use a large population-representative data set on antenatal iron and folic acid supplementation (IFAS) and birth outcomes from the Taabo HDSS in Côte d'Ivoire. We first use the Heckman model to directly test for endogenous sample selection, and then compare Heckman-corrected estimates to complete case analysis. In a second step, we also explore MICE and PIDA model to compare the relative performance of these tools in the setting studied. MICE was done with a number of 150 imputations using Stata multiple imputations (mi) package [18]. For the multiple imputations, we created a prediction model for BW with missing values from all other variables. All variables included in Appendix Table 4 were included in the imputation models.
All statistical analyses were performed in Stata version 12.0 (StataCorp; College Station, TX, USA).

Study area
The empirical data used in this study were collected through the Taabo HDSS [19,20]. The Taabo HDSS is located in the Agnéby-Tiassa region in south-central Côte d'Ivoire. It covers a surface area of approximately 980 km 2 located between latitude 6°0′ and 6°20′ N and between longitude 4°55′ and 5°15′ W.
The area is predominantly rural, with 13 main villages and more than 100 small hamlets. Within the study zone there are 11 health facilities, including seven health centres and four dispensaries in the rural area, and a 12-bed hospital located in Taabo-Cité considered as semi-urban ( Fig. 1).

Data collection
All women of reproductive age (15-49 years) from the Taabo HDSS whose pregnancy started and ended between January 1, 2012 and December 31, 2017 were included. Each household of the Taabo HDSS was visited at least three times a year during this period for detailed surveillance of vital events (i.e. birth, death, in-migration, out-migration and pregnancy). During each surveillance round, new pregnancies were systematically listed and followed-up longitudinally. When a pregnancy was completed (independent of the outcome), a standardized questionnaire on pregnancy-related morbidity was administered by field-enumerators through a personal interview with mothers [21]. This questionnaire included information on pregnancy outcome and morbidity, iron and folic acid supplementation (IFAS), birth weight, place of delivery, and birth assistance. All data were double-entered, cross-checked, and managed using a household registration system implemented in Windev version 12.0 (PC Soft, Montpellier, France) [22]. Figure 2 summarizes the main results from the Monte-Carlo simulations. One thousand random data sets with 10, 000 observations in each random draw were created and analyzed. Without missingness, estimated ordinary least squares (OLS) coefficients were normally distributed around the true causal effect of 50 as expected (Fig. 2, panel  1). With 40% missing at random (Fig. 2, panel 2), OLS is still unbiased with slightly decreased efficiency. In panels 3-5, we present results under the assumption that missingness is correlated with an unobserved determinant of birth weight. As seen in panel 3, OLS estimates are severely biased towards zero in this scenario. MICE (panel 4) changes these results only marginally. As shown in panel 5, the Heckman model is able to remove this bias completely and recovers unbiased estimates, even though the variation   Empirical application: antenatal supplementation and birth weight in the Taabo HDSS

Description of study population
Between 2012 and 2017, a total of 7619 pregnancies were reported and 7602 pregnancies were followed up after delivery (Fig. 3). Seventeen pregnancies were lost to follow-up due to out-migration of the women. Twenty records were dropped due to missing information on pregnancy-related morbidity. Overall 7542 monitored pregnancies had complete data records, and hence, were considered as final study sample. Within these fully monitored pregnancies, 7325 resulted in live births, 185 were still births, and 73 were miscarriages. Birth weight was observed for 4510 births, and unobserved for 2815 births.
Appendix Table 4 shows characteristics of women in the sample overall as well as for women who benefitted from IFAS. Over half of the women in the study had no educational attainment (54.9%) and could not write and read (73.6%). IFAS was received by 3260 (44.5%) of pregnant women. Table 1 shows the main estimation results for continuous birth weight. In fully adjusted OLS models, IFAS was associated with a non-significant 22.5 g increase (95% confidence interval (95% CI) = − 13.7, 58.7; pvalue = 0.224) in BW using OLS. In the Heckman model, the estimated increase in weight was 53.2 g (95% CI: 12.7, 93.6; p-value = 0.010). Appendix Table 5 and Appendix Table 6 compare the predicted effects of IFAS on LBW, as well as the estimated association between LBW and the other variables, from alternative multilevel logistic, MICE and Heckman probit models, respectively. In the complete cases logistic model with controls for village of residence, IFAS was not associated with higher odds of having a LBW (p-value = 0.626). Using Heckman's model to correct for endogenous selection IFAS was associated with a 10.4 percentage point reduction in the probability of LBW (95% CI: 0.169; − 0.039; p-value =0.002).

Associations between IFAS, LBW and birth weight
In both the binary dependent variable model (Appendix Table 4) and the continuous variable model (Table 1), the null hypothesis of independent residuals (cov(u,v) = 0) was rejected with p-value < 0.01. Table 2 shows the results from the selection equation. As expected, data availability was strongly correlated with socioeconomic variables as well as supplementation. Compared to women without schooling, women with secondary or higher education had an 8.0 percentage points (95% CI: 0.044, 0.117; p-value < 0.00) higher propensity to have data available. Similarly, compared to the poorest households, women from the top two wealth quintiles of households had 10.0 (95% CI: 0.068, 0.131; pvalue = 0.00) and 10.5 percentage points (95% CI: 0.066, 0.145; p-value < 0.00) higher probability of having data available. IFAS increased the probability by 24.3 percentage points (95% CI: 0.221, 0.265; p-value < 0.00).  Table 3 shows results of mean imputation, MICE and three potential PIDA scenarios. Using mean imputation and MICE, a non-significant association was found between BW and IFAS. While the estimates from the MICE model were almost identical to those found in the CCA (Table 1), mean imputation lowered the estimated association to a non-significant 7.02 g (95% CI: − 13.97; 28.01). The right hand side of Table 3 shows the PIDA results, and strongly highlights the sensitivity of the empirical model to the assumed patterns in the missing data. In PIDA scenario 1 (where missing BW data were replaced with group means) and scenario 3 (when missing BW data were replaced with BW half a SD above the mean) IFAS was associated with significant 18.7 g (95% CI: − 2.27, 39.69) and 76.8 g (95% CI: 54.10, 98.57) increase in BW. When missing values were replaced with values half a standard deviation below the mean (scenario 2), IFAS was associated with a 51.8 g (95% CI: − 73.76; − 29.80) decrease in BW. Figure 4 summarizes the estimated coefficients of all models considered and shows them relative to the latest systematic review.  Missing values were replaced with a birth weight half a standard deviation lower than the observed mean.

Estimated selection probabilities: birth weight availability
c Missing values were replaced with a birth weight half a standard deviation above the observed mean. *** p < 0.01. ** p < 0.05.* p < 0.1; (+) Reference category; Adjusted model controls village fixed effects; CI: confidence interval

Discussion
In this study, we have shown that Heckman-type selection models can be used to assess and correct potential nonrandom missingness of outcome data in the context of BW and micronutrient supplementation in low-income setting. Using simulated data, we show that bias will always emerge in standard empirical models if unobserved determinants of the outcome also predict the availability of outcome measures. Using recent data from a HDSS in Côte d'Ivoire, we then show that missingness in BW does indeed seem to correlate with unobserved maternal traits that jointly predict availability and health outcomes. This correlation between unobserved selection determinants and health outcomes leads to substantial biases in traditional regression models that cannot be removed by alternative imputation models, but generally appears to be well accounted for in Heckman models. In terms of alternative approaches, we also show that PIDA can in principle recover unbiased estimates. The main challenge with this approach is that identifying the most realistic scenario is not obvious. Given that the range of potential assumptions is rather large, PIDA methods seem most useful for illustrating the sensitivity of regression results with respect to missing data assumptions. The study presented here has several limitations. First, given the observational nature of the data, we do not know the true causal effect of IFAS in our empirical application; while we can use the latest systematic review on this intervention as reference benchmark; this benchmark does not need to necessarily hold in our setting so that we cannot directly assert the unbiasedness of the Heckman estimation. Our simulation model also assumed normal residuals, which may not always be the case. Several recent papers suggest that nonnormal residual distributions can relatively easily be incorporated in this model [23,24]. Second, it is also important to highlight that the rate of missing BW data is rather high in our study setting, so that the differences we found across models would likely be smaller in settings with better data coverage. From a health perspective, the data used in the last part of the study is relatively coarse and did not allow to separate the effects of iron and folic acid supplementation. Similarly, we were also not able to test for frequency of dosage effects of these supplements, which have been shown to be important in previous studies [25,26].

Conclusion
The results presented in this study suggest that missing outcome data can lead to substantial biases in observational studies assessing the cross-sectional associations between programme coverage and health outcomes. Heckman selection models appear to be well suited to address this potential bias and should be more widely used to address nonrandom missingness in outcome data.  For PIDA (1), missing BW data were replaced with group means. For PIDA (2) missing BW data were replaced with BW half a standard deviation below the mean. For PIDA (3), missing BW data were replaced with BWt half a standard deviation above the mean. Effect sizes (ES) represent grams, with 95% confidence intervals in parentheses