Empirical Likelihood in Nonignorable Covariate-Missing Data Problems

1 Introduction

Missing data is an important practical problem commonly found in medical sciences, social sciences, and many other related disciplines. For example, in sample surveys, some sampled individuals do not complete the questionnaire because of non-contact, refusal to respond, or some other reason. In panel surveys and clinical trials, attrition arises when members of the panel or subjects of the longitudinal study group drop out prior to the end of the study and do not return. In the past decade, there has been a great deal of interest in statistical methods for studying missing data problems, including missing response, missing covariate, or both. There are three major missing-data mechanisms, as discussed in Little and Rubin [1]. The first and simplest case is called missing completely at random (MCAR), i.e., missingness does not depend on any observed or missing quantities. The second case is missing at random (MAR), i.e., missingness depends only on the observed quantities. The third case is called missing not at random (MNAR) or nonignorable missing, i.e., missingness depends on the unobserved quantities; this is the most difficult case to handle. Most of the past works have mainly focused on ignorable missing data problems involving MCAR and MAR missing data mechanisms. For a review on the analyses of ignorable missing data problems, we refer the readers to Little and Rubin [1]. In this paper, we focus on the study of a nonignorable covariate-missing data problem motivated by the alcohol and blood pressure data from the US National Health and Nutrition Examination Survey (NHANES) [2] as described below.

The NHANES is a program of studies designed to assess the health and nutritional status of adults and children in the United States. This survey produces vital and health statistics for the Nation and provides an important basis for public health interventions. For the Demographics Data, the Examination Data, and the Questionnaire Data from the 2003–2004 NHANES, age and gender were fully observed, but the other variables such as education, systolic blood pressure (SBP), diastolic blood pressure (DBP), body mass index (BMI), income, and alcohol consumption had missing values. As argued by Little and Zhang [3], it is plausible and reasonable to assume that the missingness in education, BMI and the two blood pressure measures was missing at random (MAR) due to missed visits, but missingness of household income was thought more likely to be missing not at random (MNAR), since the propensity to respond to the income question in sample surveys is likely to depend on and possibly determined by the underlying value of the income variable, with those with low or high income generally considered less likely to respond than others. In a recent and related work, Bartlett et al. [4] have considered data on alcohol consumption and blood pressure from the 2003–2004 NHANES and have argued that missingness in alcohol consumption is likely to depend largely on alcohol consumption itself and is thus missing not at random (MNAR). To study the effect of socio-economic status (income and education) on blood pressure measures, Little and Zhang [3] have performed the regression analysis of blood pressure on income and education, adjusting for age, gender and body mass index after conditioning on the subsample of cases with household income fully observed. In addition, Bartlett et al. [4] have considered the regression analysis on the dependence of systolic blood pressure on the average number of alcoholic drinks consumed per day, with adjustment for age and body mass index. Both regression analyses can be regarded as nonignorable covariate-missing data problems.

Under the ignorable covariate-missing data mechanism, there are a variety of statistical methods available in the literature on analysis of missing data, including (a) complete case analysis (CCA), (b) inverse probability weighting (IPW) in the spirit of Horvitz and Thompson [5], which weights each of the complete data by the inverse of its inclusion probability, and (c) augmented inverse probability weighted complete-case (AIPWCC) proposed by Robins et al. [6], which incorporates the incomplete data to increase efficiency while reducing the bias. By contrast, under the nonignorable covariate-missing data mechanism, since the probability of nonignorable missingness in a covariate variable is dependent on the value of that variable, it is generally difficult to specify a model for the nonignorable covariate-missing data mechanism. To overcome this difficulty to handle nonignorable covariate-missing data problems, Bartlett et al. [4] have proposed an augmented complete case analysis by modeling the probability of missingness on the fully observed variables under the assumption that missingness in a covariate depends on the value of that covariate, but is conditionally independent of the outcome variable. Under this missing not at random (MNAR) mechanism, the commonly used CCA approach is valid and gives rise to consistent estimates, but is inefficient in that it does not make use of all observed information by disregarding all incomplete cases. Like the AIPWCC approach, the augmented CCA approach utilizes all the observed data by drawing on the information available from both complete and incomplete cases and thus improves upon the efficiency of CCA through specification of an additional model for the probability of missingness given the fully observed variables. The estimating function used in Bartlett et al. [4] can be viewed as the difference of two separate estimating functions, but other linear combinations of these two estimating functions are also possible. Thus, the question arises as how to combine the two sets of estimating functions adopted in the augmented complete case analysis of Bartlett et al. [4].

Our objective in this paper is to explore the use of empirical likelihood methods in the nonignorable covariate-missing data problem described in Little and Zhang [3] and Bartlett et al. [4] to effectively incorporate incomplete cases into the analysis of covariate-missing data and thus to improve estimation efficiency when the data are not missing at random. We propose to construct a system of unbiased estimating equations, where there are more equations than unknown parameters of interest. These estimating equations are always unbiased for any given working regression function as long as the working probability model of missingness is correctly specified. One useful feature of these unbiased estimating equations is that they naturally incorporate the incomplete data into the data analysis, making it possible to seek efficient estimation of the parameter of interest even when the working regression function is possibly not specified to be the optimal regression function. We apply the general methodology of empirical likelihood to effectively and optimally combine these unbiased estimating equations when the number of estimating equations is greater than the number of parameters of interest. Furthermore, we propose to study maximum empirical likelihood estimation of the underlying regression parameters based on the aforementioned system of unbiased estimating equations. Moreover, we apply the proposed empirical likelihood estimators to the analysis of the alcohol and blood pressure data from the 2003–2004 NHANES.

As a nonparametric method, empirical likelihood was introduced by [7, 8] for constructing confidence intervals or regions for the mean and other parameters. One advantage of using empirical likelihood is that the shape of the confidence region is determined automatically by the data. Qin and Lawless [9] have demonstrated that the empirical likelihood method can be used to solve estimating equations when the number of estimating equations exceeds the number of parameters. For theoretical developments of the empirical likelihood method, we refer readers to Hall and La Scala [10] and Owen [11]. For the analysis of missing data using the empirical or nonparametric likelihood method, see, for example, Wang and Rao [12, 13], Chen et al. [14], Liang et al. [15], Stute et al. [16], Qin and Zhang [17], and Wang and Dai [18], among others.

This paper is organized as follows. In Section 2, we propose three unbiased estimating functions and study maximum empirical likelihood estimation of the underlying regression parameters. In Section 3, we study efficiency comparison between the proposed and existing estimators. In Section 4, the proposed methodology is illustrated using a data set from the US National Health and Nutrition Examination Survey (NHANES). Simulation results are presented in Section 5 and concluding remarks are given in Section 6. Proofs of the main theoretical results are delegated to an Appendix in Section 7.

2 Methodology

In clinical trials and sample surveys, completely observed covariate information is often not available for every subject. We consider a regression analysis of an outcome variable Y on a vector Z of covariates which are always observed and a vector X of covariates which can be missing for some subjects. Our interest lies in the estimation of the unknown p×1 vector of regression parameters, β, characterized by the conditional mean model

for a specified, possibly nonlinear, link function μ(X,Z,β), where the random error ε satisfies E(ε|X,Z)=0 so that E(Y|X,Z)=μ(X,Z,β), and the joint distribution of the regressors (X,Z) is left completely unspecified. Conditional mean models include familiar linear and logistic regression models. The regression analysis with missing covariate data has received much attention recently in the statistical literature; see, for example, Little and Rubin [1], Little and Schluchter [19], Ibrahim [20], Little [21], Ibrahim and Lipsitz [22], Lipsitz and Ibrahim [23], and Ibrahim et al. [24]. Let D denote a binary non-missing indicator such that D=1 if the covariate X is observed and D=0 if X is missing. When the covariate vector X is missing at random, the probability of X being observed is commonly modeled by a parametric model for the propensity score P(D=1|Y=y,Z=z)=P(D=1|Y=y,X=x,Z=z); the parameters in this parametric propensity score model can be consistently estimated by the fully observed data on (Y,Z,D). The validity of the commonly used IPW and AIPWCC methods is contingent upon the correct specification of a propensity score model for the missingness mechanism P(D=1|Y=y,Z=z) under the MAR assumption.

Let (Y1,X1,Z1,D1),…,(Yn,Xn,Zn,Dn) denote a random sample of (Y,X,Z,D). We are interested in estimating the regression parameter β based on the observed data {(Yi,DiXi,Zi,Di), i=1,…,n} available for analysis. Let U(Y,X,Z,β) denote a p×1 specific full-data unbiased estimating function for β; a common choice for U(Y,X,Z,β) is A(X,Z,β){Y−μ(X,Z,β)} under the conditional mean model (1), where A(X,Z,β) is a vector of p known functions of (X,Z) and β. In the absence of missing data, since E{U(Y,X,Z,β)}=0, a full-data estimator βˆf of β solves the full-data estimating equation: ∑i=1nU(Yi,Xi,Zi,β)=0. With missing covariate data, a complete-case estimator βˆc of β solves the complete-case estimating equation

it is well known that βˆc can be biased when X is not missing completely at random (MCAR). Under the nonignorable missing-data mechanism in which the missingness of X is dependent on the value of X itself, we cannot directly estimate the underlying missingness mechanism P(D=1|Y=y,X=x,Z=z) on the basis of the observed data {(Yi,DiXi,Zi,Di), i=1,…,n} due to missingness of Xi whenever Di=0. As a result, the propensity score approach based on modeling P(D=1|Y=y,X=x,Z=z) is not applicable to the estimation of β under model (1) with nonignorable missingness of X; in particular, the MAR-based IPW and AIPWCC methods cannot be applied to estimate β in this context.

Under the assumption that the missingness of X is independent of the outcome variable Y conditional on covariates X and Z, namely, D and Y are conditionally independent given X and Z or D⊥Y|X,Z, the complete-case estimator βˆc is a consistent estimator of the regression parameter β. Although the complete-case analysis using the estimating eq. (2) provides valid inferences for β under the conditional independence assumption of D⊥Y|X,Z, the resulting complete-case estimator βˆc is inefficient in that it fails to draw on the observed information contained in the incomplete cases. To attempt to improve upon the efficiency of the complete-case analysis, Bartlett et al. [4] have proposed an additional model for the probability of missingness given the fully observed outcome variable Y and the fully observed covariate Z by making the same conditional independence assumption for missingness as in the complete-case analysis. Specifically, under the conditional independence assumption of D⊥Y|X,Z, the probability of X being observed is described by the probability model π0(y,z)=P(D=1|Y=y,Z=z). Let π(y,z;γ) represent a working probability model for π0(y,z), where π(y,z;γ) is a strictly positive function of (y,z,γ) and γ is a q×1 vector parameter. Let γ0 denote the truth of γ; if the working probability model is correctly specified, then π(y,z;γ0)=π0(y,z)=P(D=1|Y=y,Z=z). Since the model π(y,z;γ) for P(D=1|Y=y,Z=z) only involves fully observed variables Y and Z, we can estimate γ by the method of maximum likelihood based on the fully observed data (Y1,Z1,D1),…,(Yn,Zn,Dn). Indeed, we can estimate γ by the maximum likelihood estimator γˆn, which maximizes the binomial likelihood:

and is a solution to the following system of score equations:

where πγ(y,z;γ)=∂π(y,z;γ)/∂γ. Although the additional model π(y,z;γ) for the probability of missingness is not a model for the underlying missingness mechanism P(D=1|Y=y,X=x,Z=z) under the conditional independence assumption, it allows us to develop an alternative approach to estimation of β. Under the conditional independence assumption of D⊥Y|X,Z and based on the working probability model π(y,z;γ) for the missingness of X, Bartlett et al. [4] have proposed an augmented CCA estimation method for estimating β. To describe their method, we posit a working conditional score model for m0(Y,Z;β,ξ0)=E{U(Y,X,Z,β)|Y,Z,D=1} in terms of a parametric model m(Y,Z;β,ξ) with the same dimension as β, where m(Y,Z;β,ξ) is an arbitrary, user specified, and possibly data dependent working regression function of (Y,Z;β,ξ) and ξ is an additional finite-dimensional parameter with its true value denoted as ξ0. For each fixed (γ,ξ), let βˆacc(γ,ξ) be defined as a solution to the system of augmented complete-case estimating equations:

Then the augmented complete case (ACC) estimators of β are given, respectively, by βˆacc0=βˆacc(γ0,ξˆn) when γ0 is known and βˆacc=βˆacc(γˆn,ξˆn) when γ0 is unknown, where ξˆn is an estimator of ξ based on the complete data {(Yi,Xi,Zi):Di=1,i=1,…,n}; for example, ξˆn may be chosen to be the possibly nonlinear least squares estimator of ξ from the regression of U(Yi,Xi,Zi;βˆc) on (Yi,Zi) among the complete cases {i:Di=1}. We assume throughout that ξˆn is n-consistent for some constant ξ∗ or ξˆn−ξ∗=Op(n−1/2) under suitable regularity conditions. When the working probability model π(y,z;γ) is correctly specified for P(D=1|Y=y,Z=z), Bartlett et al. [4] have shown under suitable regularity conditions that the augmented complete case (ACC) estimator βˆacc is consistent and asymptotically normal for any working regression function m(Y,Z;β,ξ). They have also shown that for a given choice of U(Y,X,Z,β), the optimal choice for the working regression function m(Y,Z;β,ξ) is given by

When the working regression function m(y,x;β,ξ) is correctly specified to be the optimal regression function in eq. (5), we have ξ∗=ξ0. In this case, the augmented complete case estimator βˆacc improves upon the efficiency of the complete case estimator βˆc by minimizing the variance of βˆacc among all the choices of the working regression function m(Y,Z;β,ξ) for any given choice of U(Y,X,Z,β). In addition, Bartlett et al. [4] have proposed a modified estimator of β, which is guaranteed to be at least as efficient as the complete case estimator βˆc for any choice of the working regression function m(Y,Z;β,ξ).

To explore the empirical likelihood approach to the estimation of the regression parameter β under model (1) with covariate X subject to nonignorable missing, we are motivated by eqs (3) and (4) to define the following estimating functions on the basis of the working probability model π(y,z;γ) and working regression model m(y,z;β,ξ):

The first estimating function g1(Y,Z,D,X;β) is identical to the complete case estimating function in eq. (2) and has mean zero or EF{g1(Y,Z,D,X;β0)}=0 when the conditional independence assumption of D⊥Y|X,Z holds, where (Y,Z,D,X)∼F. The second estimating function g2(Y,Z,D;β,γ,ξ) involves both subjects with X observed and those with X missing and has mean zero or EF{g2(Y,Z,D;β,γ0,ξ)}=0 for any working regression model m(y,z;β,ξ) provided that the probability model π(y,z;γ) for P(D=1|Y=y,Z=z) is correctly specified, since then E{D−π(Y,Z;γ0)|Y,Z}=0. The third estimating function g3(Y,Z,D;γ) is the score function of γ and is optional when the true probability model π(y,z;γ0) for P(D=1|Y=y,Z=z) is completely known. Furthermore, when the working propensity score π(y,z;γ) is correctly specified for P(D=1|Y=y,Z=z), we have EF{g3(Y,Z,D;γ0)}=0.

3 Efficiency comparison

To compare the proposed maximum empirical likelihood estimators βˆel1, βˆel2, and βˆel3 with other existing estimators, we employ the notion of influence functions. For an exposition on influence functions, see, for example, Tsiatis [25] Chapter 3. For two matrices M1 and M2, the notations M1≤M2 and M2≥M1 mean that M2−M1 is nonnegative definite. Recall that βˆc is the complete-case estimator defined in eq. (2), and βˆacc0 and βˆacc are the augmented complete-case estimators defined in eq. (4). In addition to βˆacc, Bartlett et al. [4] have also proposed a modified estimator βˆacc2 which, for a given choice of working regression model m(y,z;β,ξ), solves

where Cˆ7=(Cˆ4−Cˆ6Sˆγ−1Cˆ2τ)(Cˆ5−Cˆ2Sˆγ−1Cˆ2τ)−1 is a consistent estimator of C7 defined in eq. (19) and

Throughout En(⋅) represents the empirical mean operator. Bartlett et al. [4] have shown that βˆacc2 is at least as efficient as βˆc.

We first compare βˆel1 with the complete-case estimator βˆc and the augmented complete-case estimator βˆacc0. According to Bartlett et al. [4], the influence function of βˆc is given by

Furthermore, it can be shown that the influence function of βˆacc0 is equal to

Moreover, the asymptotic expansion of βˆel1 in Theorem 2.1 implies that the influence function of βˆel1 is given by

To compare the asymptotic performances of βˆel1, βˆc, and βˆacc0, let

denote the linear space spanned by g2(Y,Z,D;β0,γ0,ξ∗) and let Π{g1(Y,Z,D,X;β0)|Λ2} stand for the projection of g1(Y,Z,D,X;β0) onto the space Λ2. It can be shown after some algebra that

Thus, the influence function of βˆel1 can be alternatively written as

This implies that ψel1(Y,Z,D,X;β0,γ0,ξ∗) is the residual from the projection of ψc(Y,Z,D,X;β0) onto the space Λ2. It now follows from the theory of influence functions that

In addition, since the second term in ψacc0(Y,Z,D,X;β0,γ0) belongs to the space Λ2, we have