 Research article
 Open Access
 Published:
Heckman imputation models for binary or continuous MNAR outcomes and MAR predictors
BMC Medical Research Methodology volume 18, Article number: 90 (2018)
Abstract
Background
Multiple imputation by chained equations (MICE) requires specifying a suitable conditional imputation model for each incomplete variable and then iteratively imputes the missing values. In the presence of missing not at random (MNAR) outcomes, valid statistical inference often requires joint models for missing observations and their indicators of missingness. In this study, we derived an imputation model for missing binary data with MNAR mechanism from Heckman’s model using a onestep maximum likelihood estimator. We applied this approach to improve a previously developed approach for MNAR continuous outcomes using Heckman’s model and a twostep estimator. These models allow us to use a MICE process and can thus also handle missing at random (MAR) predictors in the same MICE process.
Methods
We simulated 1000 datasets of 500 cases. We generated the following missing data mechanisms on 30% of the outcomes: MAR mechanism, weak MNAR mechanism, and strong MNAR mechanism. We then resimulated the first three cases and added an additional 30% of MAR data on a predictor, resulting in 50% of complete cases. We evaluated and compared the performance of the developed approach to that of a complete case approach and classical Heckman’s model estimates.
Results
With MNAR outcomes, only methods using Heckman’s model were unbiased, and with a MAR predictor, the developed imputation approach outperformed all the other approaches.
Conclusions
In the presence of MAR predictors, we proposed a simple approach to address MNAR binary or continuous outcomes under a Heckman assumption in a MICE procedure.
Background
In clinical epidemiology, missing data are generally classified as (i) missing completely at random (MCAR); (ii) missing at random (MAR) when, conditional on the observed data, the probability of data being missing does not depend on unobserved data; or (iii) missing not at random (MNAR) when, conditional on the observed data, the probability of data being missing still depends on unobserved data, i.e., neither MCAR nor MAR [1, 2]. Unfortunately, the missing data mechanisms of MNAR, MAR and MCAR are generally not testable unless there are direct modelisations of the missing data mechanisms. Although methods for handling MCAR or MAR data in clinical epidemiology have been widely described and studied, methods adapted for MNAR mechanisms are less studied.
In the presence of MNAR missing outcomes, valid statistical inference implies describing the missing data mechanism [1, 3]. Hence, it often requires joint models for missing outcomes and their indicators of missingness [4]. Two principal factorisations of these joint models have been proposed: patternmixture models and selection models [1, 5–7]. The first consists of using different distributions to model individuals with and without missing observations [8, 9]. The second directly models the relationship between the risk of a variable being missing and its unseen value. It involves defining an analysis model for the outcome and a selection model (i.e. the missing data mechanism). It generally relies on a bivariate distribution to model the outcome and its missing binary indicator simultaneously [10]. This approach, called sample selection model, Tobit type2 model [11] or Heckman’s model, was first introduced by Heckman for continuous outcomes [12, 13]. For continuous outcomes, two approaches have been proposed to estimate the model parameters: a onestep process that directly estimates all parameters of the joint model using the maximum likelihood estimator [11] and a twostep process [12, 13]. The first step of the latter consists of estimating the parameters of the selection model. The second step consists of fitting the outcome model adjusted on a correction term named “inverse Mills ratio” (IMR), which is obtained via the first step. IMR corresponds to the mean of the conditional distribution of the outcome within the bivariate normal distribution knowing that the outcome has been observed [14]. This allows unbiased estimates of the parameters of the outcome model to be calculated.
For binary outcomes, sample selection methods rely on a different model. This model is not simply an adaptation of the continuous case and notably is not simply an adaptation of the twostep estimator with a different outcome model as a generalised linear model. In the setting of binary outcomes, the use of a bivariate probit model and a onestep maximum likelihood estimator is mandatory [10]. Indeed, the use of a Heckman’s model implies linking the outcome model and the selection model by their error terms. Some authors, through analogy with Heckman’s twostep estimator, proposed modelling binary outcomes using a probit model adjusted on the IMR [15]. Despite the misuse of such approaches, it has been specifically demonstrated that the use of a twostep approach including the IMR in a probit model for binary outcomes is not valid [10, 16]. More generally, Heckman’s twostep estimator could not be extended straightforwardly to general linear outcome models by plugging IMR into the linear predictor. It relies on the fact that outcome expectation in nonlinear models subject to selection does not involve a simple corrector term in the linear predictor [16].
If Heckman’s model handles MNAR missing binary outcomes well using a bivariate probit model, then in the presence of additional missing data on predictors, there is no process that can address all the missing data simultaneously. In this setting, missing data on predictors are typically treated using a nonsatisfactory completepredictors approach, i.e., cases with at least one missing predictor are removed from the analysis. In the presence of missing data on more than one variable (including the outcome), multiple imputation (MI) appears to be one of the most flexible and easiest method to apply due to the numerous types of variables handled and the extensive development of statistical packages dedicated to its implementation [17]. Galimard et al. [18] previously developed an approach based on a conditional imputation model for an MNAR mechanism using a Heckman’s model and a twostep estimator to impute MNAR missing continuous outcomes. This approach allows imputing MAR missing covariates and MNAR missing outcomes within a multiple imputation by chained equations (MICE) procedure [18]. MICE specifies a suitable conditional imputation model for each incomplete variable and iteratively imputes the missing values until convergence. The key concept of MI procedures is to use the distribution of the observed data to draw a set of plausible values for the missing data. Thus, imputing missing MNAR binary outcomes implies developing valid methods to obtain a valid distribution of missing binary outcomes. As mentioned above, the direct extension of the work of Galimard et al. [18] on continuous outcomes cannot be considered because it involves a twostep estimator which is not compatible with Heckman’s model with binary outcomes.
Aims of this work
The first aim of this work is to propose an approach to handle MNAR binary outcomes. To our knowledge, the use of sample selection models as imputation models has never been proposed for missing binary outcomes, which is a current framework in clinical research. Thus, we propose developing an imputation method for binary outcomes based on a bivariate probit model associated with a onestep maximum likelihood estimator.
The second aim is to extend this approach for continuous outcomes proposing a new approach for the issue raised by Galimard et al. [18]. Indeed, for continuous outcomes, one of the main drawbacks of Heckman’s twostep estimator is that the uncertainties of the first step estimates are not taken into account in the second step. Indeed, IMR is considered as known observed values in the second step, whereas they have been estimated in the first step. Thus, the uncertainties around the final estimates are not fully assessed using a twostep estimator [19]. This point could impact the quality of the imputation. This is the reason why we hypothesised that the use of a onestep estimator could also improve the performance of Heckman’s model as an imputation model for continuous outcomes. Therefore, we also proposed a new approach for continuous missing outcomes.
The final aim is to integrate the current developed MNAR model into a MICE procedure. It will handle both MNAR outcomes and MAR predictors in the same process.
In what follows, we introduce the study that motivated this work. Then, the “Methods” section section develops our proposed imputation model using onestep ML estimation for binary and continuous outcomes. The “Results” section section presents the evaluation of its performance using a simulation study and an illustrative example using data from our motivating example. Finally, a discussion and some conclusions are provided.
Motivating example: the BIVIR study
The BIVIR study was a threearm, parallel, randomised clinical trial that aimed to assess the efficacy of the OseltamivirZanamivir combination relative to each monotherapy in patients with seasonal influenza. This study was conducted by 145 general practitioners throughout France during the 20082009 seasonal influenza epidemic and included 541 patients. Primary analyses of the trial showed that the OseltamivirZanamivir combination is less effective than Oseltamivir monotherapy and not significantly more effective than Zanamivir monotherapy based on the proportion of patients with nasal influenza reverse transcription (RT)PCR below 200 copies genome equivalent (cgeq)/ μl at day 2 after randomisation [20]. We focused our work on evaluating the impact of the treatment group on adherence adjusted on the first day severity score of flu symptoms. Adherence was defined as completing the full treatment between day 1 and day 5 and was selfreported by the patient. Unfortunately, adherence was missing for 115 (21%) patients. It was reasonable to suspect that patients who decided to stop treatment might be more likely to not record data on their adherence, resulting in an MNAR mechanism. The severity score corresponding to flu symptoms was measured as a weighted sum (ranging from 0 to 78) of 13 intensity symptoms [21]. The score was missing for 114 (21%) patients, and a MAR mechanism was suspected.
Methods
Heckman’s model
Let Y_{i} be a binary outcome and X_{i} be a pvector of covariates for individual i=1,...,n. Adopt the following probit regression model as the outcome model:
where Φ is the standard normal cumulative distribution function and β is a pvector of fixed effects. Assuming an underlying MNAR mechanism for Y, introduce a selection model that represents the nonrandom sampling of the missingness process:
where R_{yi} is an indicator of Y_{i} missingness (equal to 1 if Y_{i} is observed and 0 if Y_{i} is missing), \(X_{i}^{s}\) is a qvector of observed covariates potentially associated with the missingness mechanism, and β^{s} is an unknown qvector of coefficients.
According to the bivariate probit model, define Y^{′} and Ry′ as two latent normally distributed variables associated with Y and R_{y}, respectively, such that for individual i, Y_{i}=1 if Yi′>0 and Y_{i}=0 otherwise and R_{yi}=1 if Ryi′>0 and R_{yi}=0 otherwise. Heckman’s model considers that the two latent formulations of the selection and outcome models are linked through their error terms, which follow a bivariate normal distribution. The joint model of the outcome and selection models is defined as:
where ρ corresponds to the correlation coefficient between the error terms of the selection model \(\left (\varepsilon ^{s}_{i}\right)\) and outcome model (ε_{i}). When ρ equals 0, the selection and outcome models are independent, \(E\left (R_{yi}'Y_{i},X_{i},X_{i}^{s}\right)\) does not depend on Y_{i}, and the mechanism is MAR. When ρ is not equal to 0, \(E\left (R_{yi}'Y_{i},X_{i},X_{i}^{s}\right)\) depends on Y_{i}, and the mechanism is MNAR. The larger ρ is, the stronger the MNAR mechanism is.
For a continuous outcome, Heckman’s model given in Eq. (3) is simplified as Y_{i}, the nonlatent outcome instead of Yi′, is directly inserted in the joined model. The joint model for continuous outcomes is presented below:
where σ_{ε} is the variance of error terms (ε_{i}).
Model estimation
Maximum likelihood estimator
The parameters of the Heckman’s model (β,β^{S},ρ) are directly obtained by maximising the following loglikelihood of the joint bivariate probit model [10, 15, 19]:
where Φ_{2} corresponds to the binormal cumulative density function.
For a continuous outcome, the onestep estimator consists of estimating the parameters of the joint model (β,β^{S},ρ,σ_{ε}) via the following loglikelihood [14]:
Twostep estimator
For a continuous outcome, Heckman proposed a twostep approach to estimate the parameters of the joint model given in Eq. (4). His development comes from the expression of the following conditional expectation of the outcome [10]:
where \(\lambda _{i}=\phi \left (X^{s}_{i}\beta ^{s}\right)/\Phi \left (X^{s}_{i}\beta ^{s}\right)\) is called the “inverse Mills ratio” (IMR); ϕ corresponds to the probability density function of the normal distribution. As the IMR of each individual corresponds to an error term resulting from the probit selection model [22], Heckman proposed the following twostep procedure:

1
Estimate selection model parameters \(\left (\widehat {\beta ^{S}}\right)\) by maximum likelihood

2
For each observed i, compute \(\widehat {\lambda _{i}}\) using \(\widehat {\beta ^{S}}\)

3
Estimate \(\widehat {\beta }\) from Eq. (5)
Exclusionrestriction rule
In practice, Heckman’s model must avoid collinearity between the two linear predictors of the outcome model and the selection model. Indeed, if the variables included in the selection and outcome models are exactly the same, then E[ Y_{i}X_{i},R_{yi}=1]=X_{i}β+ρσ_{ε}λ_{i} is only identified through the IMR (λ) producing collinearity issues and possibly erroneous estimation. To avoid this concern, it has been recommended to include at least a supplementary variable in the selection equation [14, 22, 23]. Ideally, this supplementary variable should be linked to the indicator of missingness and linked to the outcome [24].
Imputation model using Heckman’s model
Under the MAR mechanism, imputation approaches use the conditional distribution of observed Y given the other covariates to impute the missing Y. However, in Heckman’s model, the conditional expectations of the observed and missing Y are different. For a binary outcome ([10], p. 921):
We propose using Eq. (6) to define the imputation model for binary outcomes.
Imputation algorithm
For a binary outcome, consider Heckman’s model parameters θ=(β,β^{s},ρ). The imputation algorithm consists of the following steps:

1
Use the onestep estimator to obtain Heckman’s model parameters \(\left (\widehat {\theta },\widehat {\Psi }\right)\) where \(\widehat {\Psi }\) is the variancecovariance matrix of \(\widehat {\theta }\)

2
Draw θ^{∗} from \(N\left (\widehat {\theta }\text {,} \widehat {\Psi }\right)\)

3
Draw \(Y^{*}_{i}\) from a Bernoulli distribution with parameter \(p^{*}_{i}\) from:
$$p^{*}_{i}=\frac{\Phi_{2}\left(X_{i}\beta^{*},X^{s}_{i} \beta^{s*},\rho^{*}\right)}{\Phi\left( X^{s}_{i} \beta^{s*}\right)}$$
For a continuous outcome, Eq. (6) becomes ([10], p. 913):
With model parameters θ=(β,β^{s},σ_{ε},ρ), in the third step of the imputation algorithm:

3
Draw Y^{∗} from:
$$Y^{*}_{i}=X_{i} \beta^{*} + \rho^{*} \sigma_{\varepsilon}^{*}\frac{\phi\left(X^{s}_{i} \beta^{s*}\right)}{\Phi\left(X^{s}_{i}\beta^{s*} \right)} + \varepsilon^{*} ~~\text{with}~~\varepsilon^{*} \sim N\left(0,~\sigma_{\varepsilon}^{*2}\right) $$
Multiple imputation by chained equations using Heckman’s imputation model
The final aim of this work is to provide a global framework to impute MNAR outcomes and MAR predictors through a MICE procedure. This procedure requires specifying conditional imputation models for each variable with missing data. The global procedure starts with an initial fill of all missing data using random draws from observed values. The posterior predictive distribution of the first incomplete variable is obtained using all observed values. Then, for a given observation with a missing value of the first variable, imputations are generated given all the other variables. Following variables with missing values are similarly repeatedly imputed in an iterated sequence. The key point of chained equation is that consecutive iterations use imputed values of the previous. Then missing value are iteratively imputed until convergence (at least 10 cycles) [17]. The theoretical properties of MICE are not well understood: except in simple cases, conditional imputation models do not correspond to any joint model [25, 26]. However, it performs well in practice [27, 28]. This procedure is realised in parallel to obtain several imputed datasets. Analyses and Rubin’s rules are then applied to obtain the final estimations of the parameters of interest.
We propose using Heckman’s imputation model for MNAR outcomes and standard imputation regression models for missing predictors, such as linear models for continuous covariates and logistic models for binary covariates. In this framework, Galimard et al. [18] proved that the missing data indicator of MNAR outcomes should be included in the imputation models of all other variables. The MICE algorithm involves defining conditional imputation models. In our case the definition of such imputation models will depend on the type of the missing mechanism:

Heckman’s imputation model for MNAR outcome, specifying outcome and selection models

General linear imputation models for MAR predictors as described by van Buuren et al. [2] adding R_{y} and the outcome to other variables in the linear predictors
Simulation study
Datagenerating process
We generated three normally independent and identically distributed variables, X_{1}, X_{2} and X_{3}, with X_{j}∼N(0,σ^{2}). Two error terms, ε and ε^{s}, were generated using ρ fixed at 0, 0.3 and 0.6 to simulate MAR, light MNAR and heavy MNAR settings from a bivariate normal distribution according to the model given in Eq. (3).
For binary outcomes, Y was generated as follows: if β_{0}+β_{1}X_{1}+β_{2}X_{2}+ε>0, then Y=1; otherwise, Y=0. The missing indicator R_{y} of Y was generated according to the following algorithm: if \(\beta _{0}^{s}+\beta _{1}^{s} X_{1}+\beta _{2}^{s} X_{2}+\beta _{3}^{s} X_{3}+\varepsilon ^{s}>0\), then R_{y}=1; otherwise, R_{y}=0.
For continuous outcomes, Y was generated according to Y=β_{0}+β_{1}X_{1}+β_{2}X_{2}+ε. Note that in that case and according to the model given in Eq. (4), σ_{ε}=1.
We fixed σ^{2} to 0.5 and (β_{0},β_{1},β_{2}) to (0,1,1). \(\left (\beta _{0}^{s},\beta _{1}^{s},\beta _{2}^{s},\beta _{3}^{s}\right)\) were fixed to (0.75,1,0.5,1), which resulted in approximately 30% missing data for the outcome.
To evaluate the robustness of our approach, we also generated a nonHeckman MNAR mechanism by directly including Y in the following selection equation: \(P(R_{y}=1)= logit\left (\beta _{0}^{sl}+X_{1}0.5 \times X_{2}+X_{3}+\beta _{Y}^{sl} Y\right)\). Two sets of parameters were considered. To obtain approximately 30% missing data on Y, we fixed \(\beta _{0}^{sl}\) to 0.60 and 0.20 for binary outcomes and to 1.31 and 1.86 for continuous outcomes, with \(\beta _{Y}^{sl}\) equal to 0, 1 and 2.
We first simulated scenarios with only missing outcomes to validate our approach in a simple setting. Then, to evaluate the performance of the MICE process, we generated missing data on X_{2} using two MAR mechanisms depending on either (X_{1},Y) or (X_{1},X_{3}). Thus, R_{2}, the indicator of X_{2} missingness, was defined by either:

P(R_{2}=1X_{1},X_{3})=Φ(0.25+X_{1}+X_{3})

\(P(R_{2}=1X_{1},Y)=\Phi \left (\beta _{0}^{R_{2}}+X_{1}+Y\right)\)
\(\beta _{0}^{R_{2}}\) was fixed to 1.10 and 0.25 for binary and continuous outcomes, respectively. We obtained approximately 30% missing data for X_{2}.
A total of N=1000 independent datasets of size 500 were generated for each setting. The sample size was chosen to be similar to our motivating example.
Analysis methods
The analysis models were probit models and linear models for binary and continuous outcomes, respectively, including X_{1} and X_{2} as predictors. The simulated data were first analysed prior to data deletion as a benchmark. The incomplete data were then analysed using the following methods:

Complete case analysis (CCA).

Heckman’s model (HEml) consisting of onestep ML estimation, as described in the “Methods” section for binary and continuous outcomes.

Multiple imputation using Heckman’s onestep ML estimation (MIHEml), as described in the “Methods” section.
For continuous outcomes exclusively, twostep approaches have also been performed.

Heckman’s twostep estimation (HE2steps) consisting of Heckman’s twostep estimator for continuous outcomes as described in the “Methods” section for continuous outcomes.

Multiple imputation using Heckman’s twostep model estimation (MIHE2steps) for continuous outcomes, as described in Galimard et al. [18].
For HEml, MIHEml, HE2steps, and MIHE2steps, the selection equation included X_{1}, X_{2} and X_{3}. For MIHEml and MIHE2steps, the incomplete data were imputed m=50 times, and final estimates were obtained by applying Rubin’s rules for small samples [29].
For scenarios with missing X_{2}: (1) for the HEml and HE2steps approaches, observations with missing X_{2} were deleted from the analyses as previously described in the completepredictors approach; (2) for MIHEml and MIHE2steps, a MICE procedure was applied. X_{2} was imputed using a linear regression model and an approximate proper imputation algorithm [2]. As recommended, we included R_{y} and Y in its imputation model [2, 18]. Twenty iterations of the chained equation process were applied.
In each datagenerating scenario, the performance of each method was assessed by computing the percent relative bias (%Rbias), the root mean square of the estimated standard error (SE_{cal}), the empirical Monte Carlo standard error (SE_{emp}), the root mean square error (RMSE) and the percent of the coverage of nominal 95% confidence intervals (Cover) of β_{1} and β_{2}.
Computational settings
Simulations and analyses were performed using R statistical software, version 3.3.0 [30]. We computed the imputation procedure within the mice R package version 2.25 [31]. Heckman’s Onestep model estimator was supplied by functions semiParBIV() and copulaSampleSel() of the GJRM R package version 0.11, for binary and continuous cases respectively [19, 32]. Our code is available in the supplementary materials (S1 for binary outcomes and S2 for continuous outcomes). Heckman’s twostep model estimator was performed using the function heckit() of package sampleSelection version 1.04 [14].
Results
In this section, only the results of β_{1} estimations are presented. β_{2} estimations are presented in Additional file 1.
Only missing data on outcome Y
Table 1 (Fig. 1) presents the results of the simulation study based on a scenario with missing binary outcome Y and complete predictors X. When Y is missing due to a MAR mechanism (ρ=0), all methods provide unbiased estimates of β_{1} (relative biases less than 2%). The standard errors of the approaches using Heckman’s model are greater than those of CCA. Nevertheless, all coverages are close to their nominal values. In the presence of an MNAR mechanism, CCA is biased 6.1% with ρ=0.3 and 11.9% with ρ=0.6. HEml and MIHEml are unbiased. The results for β_{2} are similar (Additional file 1: Table S8).
The results of the simulations that considered missing data on a continuous outcome are presented in Table 2 (Fig. 2). Compared to a binary outcome, similar results are observed. HEml, HE2steps, MIHEml and MIHE2steps presented similar results concerning biases; nevertheless, the standard errors obtained with HEml and MIHEml with ρ≠0 are smaller than those observed with HE2steps and MIHE2steps, while the confidence intervals remain near 95%. The results for β_{2} are similar (Additional file 1: Table S9).
The results of the simulations with data created using a logit selection model including Y as a covariate (i.e., a nonHeckman MNAR mechanism) are presented in Table 3 (Fig. 1) for binary outcomes and in Table 4 (Fig. 2) for continuous outcomes. CCA is not biased for \(\beta ^{sl}_{Y}=0\) and is biased for \(\beta ^{sl}_{Y} \neq 0\). The biases increase with the effect of Y. For MNAR binary outcomes, HEml and MIHEml are biased from 2.5 to 4.2% but are less biased than CCA. For continuous outcomes, HEml, HE2steps, MIHEml and MIHE2steps are slightly biased for \(\beta ^{sl}_{Y}\neq 0\), and lower standard errors are obtained using HEml and MIHEml, while biases appear to be very slightly greater.
Missing data on outcome Y and covariate X _{2}
The results of the simulations that considered missing data on a binary outcome Y and on X_{2} depending on X_{1} and X_{3} are presented in Table 5 (Fig. 1). Approximately 50% of the cases were analysed with CCA, while 70% were analysed with HEml and the entire dataset with MIHEml. Under a MAR mechanism for the missing outcome (ρ=0), the biases for CCA, HEml and MIHEml range from 1.0 to 2.1%. The smallest standard error is obtained using CCA. If the missing mechanism is MNAR, then CCA is biased from 3.8 to 8.4%, whereas the biases of HEml and MIHEml remain less than 2.5%. MIHEml provides lower standard errors than HEml notably because HEml uses only approximately 70% of the observations. The results for β_{2} are similar (Additional file 1: Table S10).
The results of the simulations that considered missing data on binary outcome Y and on X_{2} depending on X_{1} and Y are presented in Table 5 (Fig. 1). Regardless of ρ, CCA and HEml are biased from 20% to more than 33%. Regardless of ρ, MIHEml is almost unbiased (relative bias of less than 2.5%). The results for β_{2} are similar excepted for unbiased HEml (Additional file 1: Table S10).
The results of the simulation studies with missing continuous outcomes Y and missing X_{2} depending on X_{1} and X_{3} are presented for β_{1} in Table 6 (Fig. 2). When ρ=0, all methods are unbiased (relative biases of less than 1%). The smallest standard error is obtained with CCA. When ρ≠0, CCA is biased from 6.3 to 13.3%. The other methods are almost unbiased (relative biases of less than 2.2%). The results for β_{2} are similar (Additional file 1: Table S11).
The results of the simulations with missing continuous outcomes Y and missingness of X_{2} depending on X_{1} and Y are presented for β_{1} in Table 6 (Fig. 2). Regardless of ρ, CCA, HEml and HE2steps are biased from 27.7% to more than 37.7%. CCA presents the smallest standard error. Regardless of ρ, MIHEml and MIHE2steps are unbiased (relative biases of less than 2%). The standard errors observed for MIHEml are smaller than those observed for MIHE2steps, while the coverage remains close to 95%. The results for β_{2} are similar (Additional file 1: Table S11). However, when ρ=0.6, MIHEml and MIHE2steps are slightly biased for β_{2} (3.4% and 4.5%, respectively).
Similar results are observed when the sample size decreased down to 200, although biases ans SEs slightly increased (Additional file 1: Tables S12, S13, S14 and S15).
Application to illustrative examples
The impact of treatment group on adherence has been assessed using a probit model adjusted on severity score. Adherence presented 115 (21%) missing data. There were 51 and 375 nonadherent and adherent patients, respectively. The missing data mechanism of adherence was strongly suspected to be MNAR. The severity score was missing for 114 (21%) patients, and its missing data mechanism was suspected to be MAR. Four methods were applied: CCA, HEml, MIHEml and MI. A standard MI approach was added using a MICE procedure with a linear imputation model for severity score and a probit imputation model for adherence. The aim of the latter model was to assess the performance of an available misspecified but widely used approach. The missing data mechanisms assumed by each method are presented in Table 7. The HEml and MIHEml selection equations for adherence included treatment group, severity score and antibiotic treatment. The latter binary variable was chosen to fulfill the exclusionrestriction criterion. The MAR variables were imputed using linear and probit regression models for continuous and binary variables, respectively. Using MIHEml, the indicator of adherence missingness was included in the severity score imputation model. The MICE procedure was applied for 20 iterations, and m=100 datasets were generated. Finally, Rubin’s rules for small samples were applied.
The results are presented in Table 7. The reference group for treatment is the combination group. The Severity score coefficient corresponds to an increase of 20 units. CCA includes only 359 cases, i.e., 66% of the entire dataset. Observations with missing predictors are ignored in the HEml analyses, i.e., only 427 (79%) cases are retained. MI and MIHEml consider all observations. As expected, MI and MIHEml have lower standard errors than those of CCA and HEml. The coefficients estimated for OseltamivirPlacebo with MI and MIHEml are similar and higher than those obtained with CCA or HEml. The effect of OseltamivirPlacebo reached significance with MIHEml, thus enabling the assessment of the impact of OseltamivirPlacebo on adherence. The estimated coefficients of ZanamivirPlacebo and severity score are similar for CCA and MI, slightly higher for HEml and higher for MIHEml. Not surprisingly, the proportion of imputed values corresponding to the nonadherent outcome are 13% and 47% for MI and MIHEml, respectively, indicating that missing values on selfreported adherence are more likely to correspond to nonadherent patients.
We also challenged the MAR assumption concerning the missing mechanism associated with the severity score. Thus, we performed a new MICE procedure encoding two Heckman’s imputation models for adherence and severity score. It involves defining selection and outcome models for severity score. The results for the effects were similar: 0.376 (0.186) and 0.096 (0.179) for OseltamivirPlacebo and ZanamivirPlacebo, respectively. These results suggest a weak impact of the MNAR mechanism for severity score.
Discussion
The first aim of this work was to propose a unique approach to address binary outcomes according to an MNAR mechanism and missing predictors with a MAR mechanism. According to our simulation results, for MNAR outcomes, only MIHEml and HEml were unbiased. Our simulation studies were generated using a real Heckman’s model. Thus, we generated MNAR outcomes using a logistic selection model, directly including Y as a predictor, i.e. an MNAR mechanism that is noncompatible with Heckman’s model. Although our results remain biased, the use of MIHEml reduced the biases compared to CCA. Because it is not possible to confirm the validity of Heckman’s model from the observed data alone [17, 33], the developed approach appears to at least reduce the biases under an MNAR mechanism if the Heckman’s hypotheses do not hold.
To thoroughly evaluate our approach in a MICE procedure, we simulated missing data on predictors following two scenarios: one where the MAR mechanism for X_{2} depended on the fully observed X_{1} and X_{3}, and one where the mechanism depended on X_{1} and Y. For these two scenarios, Heckman’s model (HEml) used only cases with complete predictors to estimate the model parameters, i.e, did not use all available information. This loss of information produced larger standard errors, particularly for β_{1} and only slightly for β_{2}. This result is not surprising because the information lost, resulting from ignoring patients with missing X_{2}, primarily affected X_{1}. In terms of bias, the first scenario presented similar results to those obtained without missing X_{2} data. In the second scenario, where the missing mechanism for X_{2} also depended on Y, MIHEml outperformed all the other methods. The second aim was to validate the proposition of Galimard et al. [18] using a onestep ML estimator for continuous outcome. Our simulations showed that MIHEml performs slightly better than MIHE2steps in terms of standard errors for the missing MNAR outcomes.
Although our method performs well in the presence of a MAR mechanism, i.e., when ρ=0, it is preferable to determine whether the missing data mechanism is most likely to be MNAR or MAR to avoid modelling a selection equation. Indeed, the standard errors are greater than those of the standard approaches for ρ=0. Unfortunately, it is not possible to distinguish between MAR and MNAR from the observed data alone [17, 33]. Hence, sensitivity analyses are often performed to evaluate departures from MAR. Some authors have proposed a pattern mixture model using δ adjustment, i.e., systematically adding a certain increment δ to the linear predictors of the imputed values. Despite its simplicity, van Buuren considered this method to be a powerful approach for evaluating the MAR mechanism by varying δ [2, 8, 17]. This method identifies two patterns: one for the observed data and one for the unobserved data. Missing values are imputed conditionally on the observed data with an additional shift parameter δ, which is the magnitude of departure from MAR. Then, the model for the observed data is different from the model for the missing data. Similarly, MIHEml can be viewed as a method that applies a shift term or a correction term for the selection bias in the imputation model specific to each observation i. Precisely, as \(E(Y_{i}R_{yi}=0)=X_{i}\beta +\rho \sigma _{\varepsilon }\left (\phi \left (X^{s}_{i}\beta ^{s}\right)\right)/\Phi \left (X^{s}_{i}\beta ^{s}\right)\), MIHEml uses a selection correction term that can be considered as an individual δ_{i} for each patient (adjusted on the parameters of the selection equation). In this sense, we obtained a more precise δadjustment approach.
The construction of the selection model follows strict rules [14, 23]. In our experience, respect of the exclusionrestriction criterion should be strict. Indeed, Heckman’s model can inflate standard errors due to the collinearity between the regressors and IMR, and this problem is exacerbated when the exclusionrestriction criterion does not hold [34]. Moreover, MICE (or full conditional specification) follows certain rules. Each variable with missing data requires a specific conditional imputation model that is generally defined by a link function and a linear predictor with its set of predictors. Theoretically, imputation models should be derived from the global joint distribution of the variables, including the outcome [2, 35], and misspecification may result in biased parameter estimates [36]. Despite recent work in simple cases, the theoretical properties of MICE are not fully understood [25, 26, 28, 37]. Nevertheless, it performs well in practice, particularly when the conditional imputation models are well accommodated to the substantive model. The efficiency of the MICE approach is generally validated by simulation studies, and the results appear robust even when the compatibility between the full conditional distribution and the global joint distribution is not proven [2]. Although simulation is never sufficiently complete, these simulations suggest that our approach of multiple imputation using Heckman’s model and its use in a MICE process are valid and could be useful when the MNAR mechanism on the outcome is compatible with Heckman’s model. To avoid the bivariate normality assumption of Heckman’s model, Marchenko and Genton [38] proposed a Heckman’s model with a bivariate Student distribution for error terms. Ogundimu and Collins [39] developed an imputation model using this selectiont model. Unfortunately, their imputation model is only available for continuous outcomes. We compare the proposition in the current paper for continuous outcome to the propositions of Ogundimu and Collins [39] and Galimard et al. [18] in Additional file 2. Not surprisingly, the results were similar. Indeed, tdistributions are very close to a normal distribution for high degrees of freedom. In this paper, we focused on frequentist sample selection approaches within a MICE procedure. Nevertheless, Bayesian posterior distribution of sample selection models can be obtained using Gibbs sampling and data augmentation [40, 41]. Such a fully Bayesian framework could improve the imputation when based on small samples; this could be evaluated in further research.
Finally, our simulation study does not explore MNAR mechanisms on covariates and outcomes. Such a situation requires specifying a Heckman’s imputation model for each MNAR variable (i.e. selection and outcome models). Nevertheless, we used this type of approach in our example analysis to evaluate the departure from MAR for the missing predictors.
Conclusion
In the presence of MAR predictors, we proposed a simple approach to address MNAR binary or continuous missing outcomes under a Heckman assumption in a MICE procedure. This approach could be either directly used to handle such a framework (MNAR outcomes and MAR predictors) or to challenge the robustness of a suspected MAR mechanism for missing outcomes, such as in a sensitivity analysis. Finally, a R package, named “miceMNAR”, dedicated to the proposed approaches has been implemented and is available on the CRAN (https://CRAN.Rproject.org/package=miceMNAR).
Abbreviations
 CCA:

Complete case analysis
 Cover:

Coverage of nominal 95% confidence intervals
 HEml:

Heckman’s onestep ML estimation
 HE2steps:

Heckman’s twostep estimation
 IMR:

Inverse Mills ratio
 MAR:

Missing at random
 MCAR:

Missing completely at random
 MNAR:

Missing not at random
 MI:

Multiple imputation
 MICE:

Multiple imputation by chained equations
 MIHEml:

Multiple imputation using Heckman’s onestep ML estimation
 MIHE2steps:

Multiple imputation using Heckman’s twostep estimation
 Rbias:

Relative bias
 RMSE:

Root mean square error
 S E _{ cal } :

Root mean square of estimated standard errors
 S E _{ emp } :

Empirical Monte Carlo standard errors
References
 1
Little RJ, Rubin DB. Statistical Analysis with Missing Data. New York: Wiley; 2002.
 2
van Buuren S. Flexible Imputation of Missing Data. Boca Raton: CRC press; 2012.
 3
Thijs H, Molenberghs G, Michiels B, Verbeke G, Curran D. Strategies to fit patternmixture models. Biostatistics. 2002; 3(2):245–65.
 4
Fitzmaurice GM, Kenward MG, Molenberghs G, Verbeke G, Tsiatis AA. Missing data: Introduction and statistical preliminaries. In: Handbook of Missing Data Methodology. Boca Raton: Chapman and Hall/CRC Press: 2014. p. 3–22.
 5
Little RJ. Patternmixture models for multivariate incomplete data. J Am Stat Assoc. 1993; 88(421):125–34.
 6
Rubin DB. Formalizing subjective notions about the effect of nonrespondents in sample surveys. J Am Stat Assoc. 1977; 72(359):538–43.
 7
Glynn RJ, Laird NM, Rubin DB. Selection modeling versus mixture modeling with nonignorable nonresponse. In: Drawing Inferences from Selfselected Samples. New York: Springer: 1986. p. 115–42.
 8
van Buuren S, Boshuizen HC, Knook DL. Multiple imputation of missing blood pressure covariates in survival analysis. Stat Med. 1999; 18(6):681–94.
 9
Ratitch B, O’Kelly M, Tosiello R. Missing data in clinical trials: From clinical assumptions to statistical analysis using pattern mixture models. Pharm Stat. 2013; 12(6):337–47.
 10
Greene WH. Econometric Analysis: International Edition (7th Ed.)Edinburgh: Pearson; 2011.
 11
Amemiya T. Tobit models: A survey. J Econom. 1984; 24(1):3–61.
 12
Heckman JJ. The common structure of statistical models of truncation, sample selection and limited dependent variables and a simple estimator for such models. Ann Econ Soc Meas. 1976; 5(4):475–92.
 13
Heckman JJ. Sample selection bias as a specification error. Econometrica. 1979; 47(1):153–61.
 14
Toomet O, Henningsen A. Sample selection models in R: Package sampleSelection. J Stat Softw. 2008; 27(7):1–23.
 15
Van de Ven WPMM, Van Praag BMS. The demand for deductibles in private health insurance: A probit model with sample selection. J Econom. 1981; 17(2):229–52.
 16
Greene W. A stochastic frontier model with correction for sample selection. J Prod Anal. 2010; 34(1):15–24.
 17
White IR, Royston P, Wood AM. Multiple imputation using chained equations: Issues and guidance for practice. Stat Med. 2011; 30(4):377–99.
 18
Galimard JE, Chevret S, Protopopescu C, RescheRigon M. A multiple imputation approach for MNAR mechanisms compatible with Heckman’s model. Stat Med. 2016; 35(17):2907–20.
 19
Marra G, Radice R. A penalized likelihood estimation approach to semiparametric sample selection binary response modeling. Electron J Stat. 2013; 7:1432–55.
 20
Duval X, van der Werf S, Blanchon T, Mosnier A, BouscambertDuchamp M, Tibi A, Enouf V, CharloisOu C, Vincent C, Andreoletti L, Tubach F, Lina B, Mentré F, Leport C, and the Bivir Study Group. Efficacy of oseltamivirzanamivir combination compared to each monotherapy for seasonal influenza: A randomized placebocontrolled trial. PLoS Med. 2010; 7(11):1000362.
 21
Treanor JJ, Hayden FG, Vrooman PS, Barbarash R, Bettis R, Riff D, Singh S, Kinnersley N, Ward P, Mills RG, et al. Efficacy and safety of the oral neuraminidase inhibitor oseltamivir in treating acute influenza: a randomized controlled trial. JAMA. 2000; 283(8):1016–24.
 22
Vella F. Estimating models with sample selection bias: A survey. J Hum Resour. 1998; 33(1):127–69.
 23
Puhani P. The Heckman correction for sample selection and its critique. J Econ Surveys. 2000; 14(1):53–68.
 24
Marra G, Radice R, Bärnighausen T, Wood SN, McGovern ME. A simultaneous equation approach to estimating hiv prevalence with nonignorable missing responses. J Am Stat Assoc. 2017; 112(518):484–96.
 25
Chen HY. Compatibility of conditionally specified models. Stat Probab Lett. 2010; 80(7):670–7.
 26
Hughes RA, White IR, Seaman SR, Carpenter JR, Tilling K, Sterne JA. Joint modelling rationale for chained equations. BMC Med Res Methodol. 2014; 14(1):28.
 27
van Buuren S. Multiple imputation of discrete and continuous data by fully conditional specification. Stat Meth Med Res. 2007; 16(3):219–42.
 28
van Buuren S, Brand JP, GroothuisOudshoorn C, Rubin DB. Fully conditional specification in multivariate imputation. J Stat Comput Simul. 2006; 76(12):1049–64.
 29
Rubin DB. Multiple Imputation for Nonresponse in Surveys. NewYork: Wiley; 1987.
 30
R Core Team. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing; 2016. R Foundation for Statistical Computing. http://www.Rproject.org/.
 31
van Buuren S, GroothuisOudshoorn K. mice: Multivariate imputation by chained equations in R. J Stat Softw. 2011; 45(3):1–67.
 32
Marra G, Radice R. Estimation of a regression spline sample selection model. Comput Stat Data Anal. 2013; 61:158–73.
 33
Kaambwa B, Bryan S, Billingham L. Do the methods used to analyse missing data really matter? An examination of data from an observational study of intermediate care patients. BMC Res Notes. 2012; 5(1):330.
 34
Bushway S, Johnson BD, Slocum LA. Is the magic still there? the use of the Heckman twostep correction for selection bias in criminology. J Quant Criminol. 2007; 23(2):151–78.
 35
Gilks WR, Richardson S, Spiegelhalter DJ. Introducing markov chain monte carlo. In: Markov Chain Monte Carlo in Practice. Boca Raton: CRC Press: 1996. p. 75–88.
 36
Meng XL. Multipleimputation inferences with uncongenial sources of input. Stat Sci. 1994; 9:538–58.
 37
Liu J, Gelman A, Hill J, Su YS, Kropko J. On the stationary distribution of iterative imputations. Biometrika. 2014; 101(1):155–73.
 38
Marchenko YV, Genton MG. A heckman selectiont model. J Am Stat Assoc. 2012; 107(497):304–17.
 39
Ogundimu EO, Collins GS, A robust imputation method for missing responses and covariates in sample selection models. Stat Meth Med Res. 2017;0(0). https://doi.org/10.1177/0962280217715663.
 40
Kai L. Bayesian inference in a simultaneous equation model with limited dependent variables. J Econom. 1998; 85(2):387–400.
 41
van Hasselt M. Bayesian inference in a sample selection model. J Econom. 2011; 165(2):221–32.
Acknowledgements
We thank the scientific committee of the BIVIR study group for the permission to use their data (see Additional file 3).
Funding
MRR, SC and JEG are funded by Paris Diderot University, Paris, France. EC is granted by Paris Descartes University, Paris, France. MRR, SC and EC are funded by APHP (Assistance publique  Hôpitaux de Paris), Paris, France. The funding sources had no role in the study design, data collection, data analysis, data interpretation, or manuscript writing.
Availability of data and materials
The R codes for imputation models using Heckman’s model are available in additional files (Additional file 4 for binary outcome and Additional file 5 for continuous outcome) and can easily be used with the MICE package [31]. The R code corresponding to the datagenerating process can be obtained on request to JacquesEmmanuel Galimard (jacquesemmanuel.galimard@inserm.fr).
The real dataset supporting the findings (BIVIR Study) can be obtained on request to the scientific committee of the BIVIR study group by contacting Professor Catherine Leport (catherine.leport@bch.aphp.fr).
Author information
Affiliations
Contributions
JEG, SC, EC and MRR contributed to the design of the paper and the writing and revision of the manuscript. JEG performed the simulations, prepared and analysed the data. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
All the data have already been published in: “Efficacy of oseltamivirzanamivir combination compared to each monotherapy for seasonal influenza: a randomized placebocontrolled trial.” (http://dx.doi.org/10.1371/journal.pmed.1000362). This study was approved on July 18, 2008 by the Ethics Committee of Ile de France 1 (“CPP Ile de France 1”) and the French drug administration (AFSSAPS). We used already analysed data and a prespecified secondary outcome on compliance to antiviral treatment (Trial registration: http://www.ClinicalTrials.govNCT00799760).
Consent for publication
All the data have already been published in: “Efficacy of oseltamivirzanamivir combination compared to each monotherapy for seasonal influenza: a randomized placebocontrolled trial.” (http://dx.doi.org/10.1371/journal.pmed.1000362). This study was approved on July 18, 2008 by the Ethics Committee of Ile de France 1 (“CPP Ile de France 1”) and the French drug administration (AFSSAPS). We used already analysed data and a prespecified secondary outcome on compliance to antiviral treatment (Trial registration: http://www.ClinicalTrials.govNCT00799760).
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional files
Additional file 1
Additional tables. (PDF 108 kb)
Additional file 2
Comparison to Ogundimu and Collins. (PDF 129 kb)
Additional file 3
BIVIR study group. (PDF 78 kb)
Additional file 4
R code to impute binary outcome. (R 1 kb)
Additional file 5
R code to impute continuous outcome. (R 1 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Galimard, J., Chevret, S., Curis, E. et al. Heckman imputation models for binary or continuous MNAR outcomes and MAR predictors. BMC Med Res Methodol 18, 90 (2018). https://doi.org/10.1186/s1287401805471
Received:
Accepted:
Published:
Keywords
 Heckman’s model
 Missing data
 Missing not at random (MNAR)
 Multiple imputation by chained equation (MICE)
 Sample selection method