Pseudo likelihood and dimension reduction for data with nonignorable nonresponse

ABSTRACT

Tang et al. (2003. Analysis of multivariate missing data with nonignorable nonresponse. Biometrika, 90(4), 747–764) and Zhao & Shao (2015. Semiparametric pseudo-likelihoods in generalized linear models with nonignorable missing data. Journal of the American Statistical Association, 110(512), 1577–1590) proposed a pseudo likelihood approach to estimate unknown parameters in a parametric density of a response Y conditioned on a vector of covariate X, where Y is subjected to nonignorable nonersponse, X is always observed, and the propensity of whether or not Y is observed conditioned on Y and X is completely unspecified. To identify parameters, Zhao & Shao (2015. Semiparametric pseudo-likelihoods in generalized linear models with nonignorable missing data. Journal of the American Statistical Association, 110(512), 1577–1590) assumed that X can be decomposed into U and Z, where Z can be excluded from the propensity but is related with Y even conditioned on U. The pseudo likelihood involves the estimation of the joint density of U and Z. When this density is estimated nonparametrically, in this paper we apply sufficient dimension reduction to reduce the dimension of U for efficient estimation. Consistency and asymptotic normality of the proposed estimators are established. Simulation results are presented to study the finite sample performance of the proposed estimators.

Keywords:

• Dimension reduction
• kernel estimation
• nonignorable nonresponse
• nonresponse instrument
• pseudo likelihood

1. Introduction

Missing data or nonresponse is common in various statistical applications such as sample surveys and biomedical studies. Let Y be a univariate response variable subject to nonresponse, X be a vector of covariates that is always observed, and R be the indicator of whether Y is observed. When the propensity , missing data are ignorable and there is a rich literature on methodology of handling nonresponse (Little & Rubin, 2014). In many applications, however, Y cannot be excluded from the propensity and missing data are nonignorable. With nonignorable nonresponse, it is challenging to estimate unknown characteristics in the conditional distribution of Y given X or the unconditional distribution of Y.

Throughout we use or as a generic notation for the conditional or unconditional probability density with respect to an appropriate measure. When nonresponse is nonignorable and both and are nonparametric, is not identificable (Robins & Ritov, 1997). When both and have parametric forms, maximum likelihood methods have been developed (Baker & Laird, 1988; Greenlees, Reece, & Zieschang, 1982). Since parametric methods are sensitive to model violations, efforts have been made under semiparametric models. Qin, Leung, and Shao (2002) and Wang, Shao, and Kim (2014) imposed a parametric model on but allowed to be nonparametric. Assuming , i.e., the entire covariate vector X can be excluded from the propensity, Tang, Little, and Raghunathan (2003) proposed a pseudo likelihood method in which is nonparametric and is parametric: (1) where θ is an unknown parameter vector and is a conditional density which is known when θ is known. Zhao and Shao (2015) extended the pseudo likelihood method to the case where part of X can be excluded from the propensity, i.e., (2) where and the covariate Z, termed as instrumental variable, cannot be excluded from in (1(1) ). Here is a brief description of what has been done under (1(1) )–(2(2) ). Under (1(1) )–(2(2) ) and the Bayes formula, (3) Let , , be n independent and identically distributed observations from . Based on (3(3) ), if is known, then we can estimate θ by maximising the following likelihood: (4) Usually is unknown. Substituting in (4(4) ) by its estimate results in a pseudo likelihood. For example, Zhao and Shao (2015) assumed a parametric model for and replaced in (4(4) ) by , where is an estimator of η based on , , and is the empirical distribution of , . To avoid model misspecification on , Zhao and Shao (2015) also suggested a nonparametric kernel estimator to replace in (4(4) ).

However, kernel estimation is unstable when the dimension of is not small. The purpose of our work is to propose an alternative way to handle in (4(4) ), which adopts dimension reduction techniques to improve the resulting pseudo likelihood estimator of θ. Our main idea is described in the next section, along with the proposed pseudo likelihood and estimator of θ. Under typical conditions for kernel estimation, our proposed pseudo likelihood estimator is asymptotically normal with convergence rate . We also perform some simulations to examine the finite sample properties of our proposed estimator.

2. Methodology and theory

As pointed out in the previous section, the main contribution of this paper is to estimate in the denominator of (4(4) ) in a more efficient way, especially when the covariate U is of high dimension. Often times, the dimension of the instrument variable Z is small, while the covariate U contains a lot of variables (demographic variables for instance), and its dimension p is not small. A straightforward idea is to split as , instead of as in Zhao and Shao (2015). Since does not involve θ, the likelihood in (4(4) ) is equivalent to (5) and our main task is to estimate the denominator in (5(5) ), i.e., the integral (6) for arbitrarily given θ, y and u. Since the real form of is not our main concern, for robustness, we adopt a nonparametric kernel regression, called Nadaraya-Watson (NW) estimation (Nadaraya, 1964; Watson, 1964), to estimate the conditional expectation in (6(6) ). The way we split , along with the nonparametric estimation, actually frees ourselves from parameterising or modelling the relationship between Z and U.

Before building our estimators, we introduce a generic notation for a kernel with an appropriate dimension and bandwidth h, i.e., appeared in different places may be different. In what follows is chosen to be a product kernel of dimension s and order in the sense that , where is the jth component of the s-dimensional x and is a bounded and Lipschitz continuous univariate kernel having a compact support and satisfying , is finite and nonzero, and for all 0<l<m. The NW estimator of in (6(6) ) is (7) Substituting this estimator into (5(5) ) leads to a maximum pseudo likelihood estimator of θ, (8) Although p, the dimension of U, does not exert a direct influence on the convergence rate of as shown in Theorem 2.1, it is well known in the literature that kernel estimators do not perform well when p is very large. For example, their convergence requires a very large sample size n.

Fortunately, there is a well-developed nonparametric method called Sufficient Dimension Reduction (SDR) (Cook & Weisberg, 1991; Li & Wang, 2007; Li, 1991; Ma & Zhu, 2012; Xia, Tong, Li, & Zhu, 2002), which helps to reduce the dimension of predictors by finding a matrix B with the smallest possible such that , meaning that Z is only related to U via , where is the transpose of B. It is common that d<p and we are able to improve the estimation in (7(7) ) and (8(8) ) by applying NW estimation directly to Z and with being an SDR estimator of B. Starting from sliced inverse regression (SIR) (Li, 1991), research has been done in the literature to develop SDR estimators of B, including sliced average variance estimation (SAVE) (Cook & Weisberg, 1991), directional regression (DR) (Li & Wang, 2007), (conditional) minimum average variance estimation (MAVE ) (Xia et al., 2002), semiparametric approach to dimension reduction (Ma & Zhu, 2012) and etc. We adopt SIR method developed by Li (1991) to estimate B in our simulation studies in Section 3 because it is easy to implement and works out well in practice.

When , in (6(6) ) is equal to , which can be estimated by a new kernel estimator Then, a new maximum pseudo likelihood estimator of θ, (9) In the rest of this section we establish asymptotic consistency and normality of our proposed estimator in (8(8) ) and in (9(9) ). We first introduce some notation. Let and be the first derivative with respect to θ. For the estimator without SDR, denote (10) and the kernel estimator of γ is denoted by (11) Therefore, maximising the pseudo likelihood in (8(8) ) is the same as maximising the pseudo log-likelihood (12) Differentiating with respect to θ, we obtain the score function (13) Thus, can also be obtained by solving .

We use the same notation γ as in (10(10) ) for the estimator with SDR, (14) Its kernel estimator is (15) The pseudo log-likelihood and score function of the can be obtained by replacing in (12(12) ) and (13(13) ) by .

As the structure of two etimators are alike, we only provide the sufficient conditions for in Assumptions 2.1–2.2 where γ is defined as in (10(10) ). For , however, Assumptions 2.1–2.2 need to be modified by replacing (10(10) ) by (14(14) ), (11(11) ) by (15(15) ) and p by d. In addition, in order for to converge with rate , should also possess some asymptotic properties as indicated in Assumption 2.3 with γ defined in (14(14) ). In fact, SDR estimator and kernel estimators and all have good asymptotic performances (Bierens, 1987; Li, 1991; Nadaraya, 1964; Watson, 1964).

Throughout, denotes the true but unknown value of θ.

Assumption 2.1

There exists a constant c>0 such that . The function has bounded mth derivative with respect to u. There exists a q>1 such that as and is bounded.

Let be the sup-norm. Under Assumption 2.1, and , which means that the estimator converges uniformly to γ at a certain rate, i.e., (Hansen, 2008; Newey & McFadden, 1994).

Assumption 2.2

1. There exists such that where is defined in (12(12) ).

2. for with small enough and and are bounded as functions of y and u.

3. For small enough is continuously differentiable in θ on a neighbourhood of where is defined in (13(13) ). There exists a with such that for an . exists and is nonsingular.

4. and as .

Assumptions 2.1–2.2 together guarantee that is consistent for . As to , there is an extra step of SDR estimation requiring Assumption 2.3.

Assumption 2.3

Let for some where is the dimension of z, and γ is as defined in (14(14) ).

1. Uniformly in Ω, The mth derivatives of and on are locally Lipschitz-continuous as functions of . and each entry in the matrices are locally Lipschitz-continuous and bounded from above as a functions of .

2. in probability as and where .

3. The bandwidth for .

Note that B is a matrix except when d=1, and when it comes to the calculation of derivatives or norms, we are treating B as , which denotes the vector formed by concatenating the columns of B. For simplicity, we are still using B to represent . Similarly, we denote as .

Theorem 2.1

1. If Assumptions 2.1–2.2 hold. Then, as in probability and where means convergence in distribution, and

2. If Assumptions 2.1–2.2 hold with (10(10) ), (11(11) ) and p replaced by (14(14) ), (15(15) ) and d, respectively, and if Assumption 2.3 also holds. Then, as in probability and where

It follows from Theorem 2.1 that the asymptotic variance of may or may not be smaller than that of , partly due to the variability of SDR estimation. Owing to SDR, the dimension of kernel estimation is reduced from p to d, which is the main advantage of over , although they have the same convergence rate. A problem with a not so small p is the selection of bandwidth h for kernel estimation. Specifically, if with a , then from Assumptions 2.1 and 2.2(iv), τ must be between and for , where m is the order of kernel. When p is not so small, the upper bound is pretty small, and we may need to increase the kernel order m in order to find such τ satisfying the constraints. If m>2, the kernel takes negative values, which may reduce the stability of kernel estimation. For example, if q=2 in Assumption 2.1, then we must have and we cannot use m=2 when p>3; if p=8, then we must use a kernel of order m=5. On the other hand, if we reduce U to with dimension d, the upper bound increases to , allowing enough flexibility in choosing τ to be larger than the lower bound , which results in a low-order kernel. If q=2, then we may use m=2 when.

Although both and have convergence rate , may still have an edge over in finite sample performance, because of smaller dimension and order used in the kernel estimation. This is supported by the simulation results in Section 3.

3. Simulation studies

We study the finite-sample performance of the proposed pseudo likelihood estimators in two simulation studies. Four estimators are compared: the estimator based on full data assuming no missing data (), the estimator based on complete case analysis (), the maximum pseudo likelihood estimator without dimensional reduction ( defined in (8(8) )), and the maximum pseudo likelihood estimator with dimensional reduction ( defined in (9(9) )). The NW estimations are computed using a Gaussian kernel . The bandwidths are selected by rule of thumb proposed by Silverman (1986), i.e., , where is the estimated standard deviation of U. For SDR, we apply SIR in Li (1991) with the number of slices . All results are based on 1000 simulation replications and sample size n=400.

3.1. Experiment 1

In the first experiment, , where and . Hence, the dimension p equals 6. To evaluate the performance of SDR, two models of are studied: and , where d equals 1 and 2, respectively. The response Y is generated as , where . The propensity , resulting in unconditional response rates 73% and 74% in two cases, respectively. Table 1 reports the estimators and their standard deviations.

Table 1. Simulation results of experiment 1.

		Estimator				Standard Deviation
d=1		0.9935	1.8226	1.5812	0.9523	0.1526	0.2002	0.2429	0.2567
		1.0017	0.8448	0.8462	1.0133	0.0842	0.0950	0.1072	0.1042
		1.0005	0.8417	0.8285	1.0031	0.0819	0.0946	0.1094	0.1066
		1.0046	0.8449	0.8331	1.0065	0.0850	0.0957	0.1062	0.1086
		1.0016	0.8448	0.8403	1.0161	0.0818	0.0942	0.1031	0.1060
		−1.0008	−0.9620	−0.9230	−0.9979	0.0416	0.0469	0.0612	0.0530
		1.9725	1.7240	1.9107	1.9991	0.1415	0.1442	0.1854	0.2265
d=2		0.9995	1.8152	1.3759	0.9347	0.1654	0.2637	0.2453	0.2402
		0.9959	0.7565	0.9055	1.0319	0.0910	0.1231	0.1164	0.1226
		0.9960	0.7596	0.9031	1.0323	0.0906	0.1240	0.1186	0.1271
		1.0001	0.7629	0.9090	1.0340	0.0907	0.1240	0.1189	0.1283
		1.0020	0.7630	0.9079	1.0362	0.0914	0.1201	0.1163	0.1206
		−0.9983	−0.8816	−0.9408	−1.0077	0.0436	0.0588	0.0623	0.0647
		1.9750	1.7556	1.9051	2.0024	0.1429	0.1479	0.1884	0.2170

The simulation results can be summarised as follows. First, is nearly unbiased for θ no matter d equals 1 or 2. Second, the bias of may not be negligible and, in some cases, the bias is comparable to that is biased in theory. As we discussed in Section 2, compared with , the assumptions are harder to satisfy for . Third, the standard deviations of and are quite close and, in some cases, has slightly smaller standard deviation than that of .

3.2. Experiment 2

In the second experiment, we consider a binary outcome Y . Covariate U is generated the same as that in the first experiment, , or . The binary Y is generated according to , where . The propensity , resulting unconditional response rates 71% or 72%. The simulation results are reported in Table 2. The conclusions are quite similar to those for experiment 1, but the bias of is more serious: the bias of are even larger than the bias of . Meanwhile, is still nearly unbiased.

Table 2. Simulation results of experiment 2.

		Estimator				Standard Deviation
d=1		−1.0369	−0.9579	−3.7354	−0.9329	0.3671	0.5210	0.5141	0.5301
		1.0307	1.1445	1.8894	0.9512	0.1858	0.2649	0.2733	0.2544
		1.0245	1.1253	1.8723	0.9524	0.1911	0.2626	0.2849	0.2657
		1.0236	1.1281	1.8988	0.9612	0.1883	0.2711	0.2782	0.2648
		1.0283	1.1376	1.9018	0.9577	0.1854	0.2623	0.2735	0.2613
		−1.0265	−1.0499	−1.1685	−0.9859	0.1498	0.2181	0.2105	0.2120
d=2		−1.0307	−0.9436	−4.6084	−1.0587	0.4623	0.6637	0.6068	0.6450
		1.0271	1.1336	1.9316	1.0554	0.2154	0.3046	0.3262	0.3091
		1.0318	1.1401	1.9586	1.0627	0.2085	0.3064	0.3353	0.3147
		1.0227	1.1270	2.0029	1.0553	0.2043	0.2973	0.3217	0.2928
		1.0268	1.1367	1.9532	1.0458	0.2140	0.3070	0.3324	0.3011
		−1.0275	−1.0568	−1.2408	−1.0294	0.1495	0.2176	0.2011	0.2104

4. Proofs

Proof of Theorem 2.1

We first prove part (i). Here, γ is defined as in (10(10) ) and is defined as in (11(11) ). When , denote as . Note that, in Assumption 2.1, we assume that there exists a constant c>0 such that . Moreover, the univariate kernel is bounded, hence, there exists a constant b>0 such that .

For the consistency of , we only need to verify the assumptions of Theorem 2(a) in Zhao and Shao (2015). Note that Assumptions 2.1 and 2.2(i), (iv) guarantee that is consistent. Then, Hence, Since , in probability as . Therefore, Following Theorem 2 in Zhao and Shao (2015), with Assumption 2.2(i), in probability as , where is the true value of θ.

For the asymptotic normality of , we only need to verify the assumptions of Theorem 8.11 in Newey and McFadden (1994).

Step 1. Let , and We would like to prove Note that (16) which follows from the Assumptions 2.1 and 2.2(ii).

Step 2. Let We would like to prove and . For , V-statistics method is applied. We use the result of Lemma 8.4 by Newey and McFadden (1994) directly. In our case, . Then Let , then where and Then, following Assumption 2.2(ii) and (iv), we obtain that Hence, .

As to , by Chebychev's Inequality, since , which follows the Assumptions 2.1 and 2.2(ii).

Step 3. Note that Let , and Then, We would like to prove By Chebychev's Inequality, we only need to prove So we only need to prove and . Let Then, Therefore, by Assumption 2.2(ii) and (iv), . And Note that, for some constant C, and Since when n is large enought, and is bounded, Then by conditional dominated convergence theorem, By Dominated Convergence Theorem, .

Then, we combine Step 1–3, based on Theorems 8.11 and 8.12 in Newey and McFadden (1994), with Assumption 2.2(iii), and .

Second, we prove part (ii). Here, γ is defined as in (10(10) ) and is defined as in (15(15) ). We consider the normality of . When is replaced by the true SDR direction B, the proof is exactly the same as in part (i). Then we could easily get that . Now we only need to consider Based on Lemma 3 in Ma and Zhu (2012), with Assumption 2.3(i) and (iii), And, with Assumption 2.3(ii), we can easily prove the normality of Then we finish the proof of part (ii).

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by Division of Mathematical Sciences [1612873] and the Chinese Ministry of Education 111 Project [B14019].

Notes on contributors

Ji Chen

Ji Chen is a PhD candidate in East China Normal University.

Bingying Xie

Bingying Xie is a statistician at Roche in Shanghai, China.

Jun Shao

Jun Shao is a professor in University of Wisconsin-Madison.