Missing data approaches for probability regression models with missing outcomes with applications

In this paper, we investigate several well known approaches for missing data and their relationships for the parametric probability regression model Pβ(Y|X) when outcome of interest Y is subject to missingness. We explore the relationships between the mean score method, the inverse probability weighting (IPW) method and the augmented inverse probability weighted (AIPW) method with some interesting findings. The asymptotic distributions of the IPW and AIPW estimators are derived and their efficiencies are compared. Our analysis details how efficiency may be gained from the AIPW estimator over the IPW estimator through estimation of validation probability and augmentation. We show that the AIPW estimator that is based on augmentation using the full set of observed variables is more efficient than the AIPW estimator that is based on augmentation using a subset of observed variables. The developed approaches are applied to Poisson regression model with missing outcomes based on auxiliary outcomes and a validated sample for true outcomes. We show that, by stratifying based on a set of discrete variables, the proposed statistical procedure can be formulated to analyze automated records that only contain summarized information at categorical levels. The proposed methods are applied to analyze influenza vaccine efficacy for an influenza vaccine study conducted in Temple-Belton, Texas during the 2000-2001 influenza season.


Introduction
Suppose that Y is the outcome of interest and X is a covariate vector. One is often interested in the probability regression model P β (Y|X) that relates Y to X. In many medical and This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/ licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. epidemiological studies, the complete observations on Y may not be available for all study subjects because of time, cost, or ethical concerns. In some situations, an easily measured but less accurate outcome named auxiliary outcome variable, A, is supplemented. The relationship between the true outcome Y and the auxiliary outcome A in the available observations can inform about the missing values of Y. Let V be a subsample of the study subjects, termed the validation sample, for which both true and auxiliary outcomes are available. Thus observations on (X, Y, A) are available for the subjects in V and only (X, A) are observed for those not in V.
It is well known that the complete-case analysis, which uses only subjects who have all variables observed, can be biased and inefficient, cf. Little and Rubin (2002). The issues also rise when substituting auxiliary outcome for true outcome; see Ellenberg and Hamilton (1989), Prentice (1989) and Fleming (1992). Inverse probability weighting (IPW) is a statistical technique developed for surveys by Horvitz and Thompson (1952) to calculate statistics standardized to a population different from that in which the data was collected. This approach has been generalized to many aspects of statistics under various frameworks. In particular, the IPW approach is used to account for missing data through inflating the weight for subjects who are underrepresented due to missingness. The method can potentially reduce the bias of the complete-case estimator when weighting is correctly specified. However, this approach has been shown to be inefficient in several situations, see Clayton et al. (1998) and Scharfstein et al. (1999). Robins et al. (1994) developed an improved augmented inverse probability weighted (AIPW) complete-case estimation procedure. The method is more efficient and possesses double robustness property. The multiple imputation described in Rubin (1987) has been routinely used to handle missing data. Carpenter et al. (2006) compared the multiple imputation with IPW and AIPW, and found AIPW as an attractive alternative in terms of double robustness and efficiency. Using the maximum likelihood estimation (MLE) coupled with the EM-algorithm (Dempster et al. 1977), Pepe et al. (1994) proposed the mean score method for the regression model P β (Y|X) when both X and A are discrete.
In this paper, we investigate several well known approaches for missing data and their relationships for the parametric probability regression model P β (Y|X) when outcome of interest Y is subject to missingness. We explore the relationships between the mean score method, IPW and AIPW with some interesting findings. Our analysis details how efficiency is gained from the AIPW estimator over the IPW estimator through estimation of validation probability and augmentation to the IPW score function. Applying the developed missing data methods, we derive the estimation procedures for Poisson regression model with missing outcomes based on auxiliary outcomes and a validated sample for true outcomes. Further, we show that the proposed statistical procedures can be formulated to analyze automated records that only contain aggregated information at categorical levels, without using observations at individual levels.
The rest of the paper is organized as follows. Section 2 introduces several missing data approaches for the probability regression model P β (Y|X), where the outcome Y may be missing. Section 3 explores the relationships among these estimators. The asymptotic distributions of the IPW and AIPW estimators are derived and their efficiencies are compared. Section 3 investigates efficiency of two AIPW estimators, one is based on the augmentation using a subset of observed variables and the other is based on the augmentation using the full set of observed variables. The procedures for Poisson regression using automated data with missing outcomes are derived in Section 4. The finite-sample performances of the estimators are studied in simulations in Section 5. The proposed method is applied to analyze influenza vaccine efficacy for an influenza vaccine study conducted in Temple-Belton, Texas during the 2000-2001 influenza season. The proofs of the main results are given in the Appendix A, while the proof of a simplified variance formula in Section 4 is placed in the Appendix B.

Missing data approaches
Consider the probability regression model P β (Y|X), where Y is the outcome of interest and X is a covariate vector. Let A be the auxiliary outcome for Y and V be the validation set such that observations on (X, Y, A) are available for the subjects in V and only (X, A) are observed for those in V, the complement of V. In practice, the validation sample may be selected based on the characteristics of a subset, Z, of the covariates in X. We write X = (Z, Z c ). For example, Z may include exposure indicator and other discrete covariates and Z c may be the exposure time.
Most statistical methods for missing data require some assumptions on missingness mechanisms. The commonly used ones are missing completely at random (MCAR) and missing at random (MAR). MCAR assumes that the probability of missingness in a variable is independent of any characteristics of the subjects. MAR assumes that the probability that a variable is missing depends only on observed variables. In practice, if missingness is a result by design, it is often convenient to let the missing probability depend on the categorical variables only. There is also simplicity in statistical inference by modeling the missing probability based on the categorical variables. We introduce the following missing at random assumptions.
Since the conditional density f(y, z c |ξ, z, a) = f(z c |ξ, z, a)f(y|z c , ξ, z, a) = f(z c |z, a) f(y|z c , z, a) = f(y, z c |z, a), MAR I implies MAR II. It is also easy to show that MAR II implies MAR.
Let πî be the estimator of the conditional probability π i = P(ξ i = 1|X i , A i ), and the estimator of . Let S β (Y|X) denote the partial derivatives of log P β (Y|X) with respect to β. Let Ê{S β (Y | X i )|X i , A i } be the estimator of the conditional expectation (1) where W i takes one of the following forms: The estimator βÎ 1 obtained by using is an IPW estimator where a subject's validation probability depends only on the category defined by (Z i , A i ). Because , the estimator βÎ 1 is approximately unbiased. The estimator βÎ 2 obtained by using is also an IPW estimator but with the validation probability π i depending on the category defined by (Z i , A i ) and the additional covariate .
The estimator βÊ 1 obtained by using is the mean score estimator where the scores S β (Y i | X i ) for those with missing outcomes are replaced by the estimated conditional expectations given (Z i , A i ). The estimator βÊ 2 obtained by using is the mean score estimator where the scores S β (Y i |X i ) for those with missing outcomes are replaced by the estimated conditional expectations given (X i , A i ). The estimator βÊ 2 is the mean score estimator in The estimator βÂ 1 obtained using is the AIPW estimator augmented with the estimated conditional expectation Ê {S β (Y|X i )|Z i , A i }. The estimator βÂ 2 obtained using is the AIPW estimator augmented with the estimated conditional expectation Ê{S β (Y|X i )|X i , A i }.
The estimator βÂ 3 is obtained using . The differs from in that the estimated validation probability is πî instead of .
Suppose that is an asymptotically unbiased estimator of and that holds, then which entails that the estimator βÂ 1 has the double robust property in the sense that it is a consistent estimator of β if either is a consistent estimator of or Similarly, under MAR I, the estimator βÂ 2 possesses the double robust property in that βÂ 2 is a consistent estimator of β if either is a The estimator βÂ 3 has similar double robust property as βÂ 2 .

Method comparisons and asymptotic results
Let V(X i , A i ) denote the subjects in V with values of (X, A) equal to (X i , A i ), n V (X i , A i ) the number of subjects in V(X i , A i ), and n(X i , A i ) the number of subjects with values of (X, A) equal to (X i , A i ). When X and A are discrete and their dimensionality is reasonably small, the probability π i = P(ξ i = 1|X i , A i ) can be estimated by πî = n V (X i , A i )/n(X i , A i ). The conditional expectation E {S β (Y|X i )|X i , A i } can be nonparametrically estimated based on the validation sample, Proposition 1. Suppose that X = (Z, Z c ) and A are discrete and their dimensionality is reasonably small. Under the nonparametric estimators , πî = n V (X i , A i )/n(X i , A i ) and the estimators for the conditional expectation defined in (9) and (10), the estimators βÎ 1 , βÊ 1 and βÂ 1 are equivalent, and the estimators βÎ 2 , βÊ 2 , βÂ 2 and βÂ 3 are equivalent. However, the estimator βÂ 2 is different from βÂ 1 unless is linearly related to Z i in which case β is not identifiable.
The results of Proposition 1 are very intriguing since research has shown that the AIPW and the mean score methods are more efficient than the IPW method. It is also intriguing that the AIPW estimators βÂ 2 and βÂ 3 are actually the same estimators, not affected by the validation probability. To further understand these approaches, we investigate the asymptotic properties of these methods where (X, A) are not necessarily discrete. Through the asymptotic analysis, we gain insights about what matters to the efficiency in terms of the selections of the validation sample and the augmentation function.
be the parametric model for the validation probability π i , where ψ is a qdimensional parameter. We show in Corollary 2 that the nonparametric estimator of π(X i , A i , ψ) can also be expressed in the parametric form when (X i , A i ) are discrete. Let ψ 0 be the true value of ψ. Under MAR I, the MLE ψ̂ = (ψ̂1, …, ψq) of ψ = (ψ 1 , …, ψ q ) is obtained by maximizing the observed data likelihood, The validation probability π i is estimated by πĩ = π (X i , A i , ψ). Then by the standard likelihood based analysis, we have the approximation (11) where and I ѱ are the score vector and information matrix for ѱ̂ defined by (12) where a ⊗2 = aa′.
Consider the IPW estimator βÎ obtained by solving the estimating equation  (13) and the AIPW estimator βÂ based on solving the estimating equation (14) Theorem 1. Assume that P β (Y|X) and π(X, A, ѱ) have bounded third-order derivatives in a neighborhood of the true parameters and are bounded away from 0 almost surely, both −E {(∂ 2 /∂β 2 ) (log P β (Y|X))} and I ѱ are positive definite at the true parameters. Then, under MAR I, and .
Both n 1/2 (β̂I − β) and n 1/2 (β̂A − β) have asymptotically normal distributions with mean zero and covariances equal to and , respectively. Further, (15) and (16) where Suppose that the validation probability Suppose that πĩ is the MLE of under the parametric family ѱ = (Z i , A i , ѱ). Let βÂ 1 be the estimator obtained by solving (14) where the augmented term, Let βÂ 2 be the estimator obtained by solving (14) where The following corollary presents the asymptotic results for two AIPW estimators of β, one that corresponds to the augmentation based on a subset, (Z, A), of observed variables and the other that corresponds to the augmentation based on the full set, (X, A), of the observed variables.
Corollary 1. Suppose that the validation probability and (18) where and . The asymptotic variance of βÂ 2 is smaller than the asymptotic variance of βÂ 1 if the covariates Z i are a proper subset of X i .
Suppose that (Z, A) are discrete taking values (z, a) in a set Ƶ of finite number of values. If the number of parameters in ѱ equals the number of values ѱ z, a = P(ξ i = 1 |Z i =z, A i = a) for all distinct pairs (z, a), then ѱ = {ѱ z, a } and n(z, a, ѱ) = ѱ z, a . Further, can be viewed as a column vector with 1 in the position for ѱ z, a and 0 elsewhere. The information matrix I ѱ defined in (12) has the expression, where ρ(z, a) = P(Z i = z, A i = a). It follows that I ѱ is a diagonal matrix and its inverse matrix is also diagonal. The MLE ѱẑ , a = n V (z, a)/n(z, a) is in fact the nonparametric estimator for ѱ z, a based on the proportion of validated samples in the category specified by (z, a). The equation (11) can be expressed as for (z, a) ∈ Ƶ. By Threom 1, the possible efficiency gain of the AIPW estimator over the IPW estimator is shown through the equation (15). The AIPW estimator is more efficient unless Var(B i + O i ) = 0. In particular, from the proof of Theorem 1, we have (19) (20) where B i and O i are defined following (16). The following corollary presents the analysis of the term when (Z i , A i ) are discrete to understand how efficiency may be gained from the AIPW estimator over the IPW estimator.
By Corollary 2, (19) and (20), βÂ is more efficient than βÎ unless Var{S β (Y j |X j )|Z j = z, A j = a} = 0 for all (z, a) for which P(Z i = z, A i = a) ≠ 0. If X = Z and the validation probability π i = P(ξ i = 1 |X i , A i ) is nonparametrically estimated with the cell frequencies ѱẑ , a = n V (z, A)/n(z, a), then βâ and βÎ are asymptotically equivalent.
Remark Consider the estimators of β obtained based on the estimating equation (1) corresponding to different choices of W i given in (2) to (8). If (Z, A) are discrete and the validation probability is estimated nonparametrically by the cell frequency, then by Theorem 1 and Corollary 2, βÂ 1 and βÎ 1 have same asymptotic normal distributions as long as are estimated nonparametrically or based on some parametric models. In addition, by Theorem 1, Corollary 1 and 2, βÂ 3 and βÎ 2 have the same asymptotic normal distributions as long as

Poisson regression using the automated data with missing outcomes
Many medical and public health data are available only in aggregated format, where the variables of interest are aggregated counts without being available at individual levels. Many existing statistical methods for missing data require observations at individual levels. Applying the missing data methods presented in Section 3, we derive some estimation procedures for the Poisson regression model with missing outcomes based on auxiliary outcomes and a validated sample for true outcomes. Further, we show that, by stratifying based on a set of discrete variables, the proposed statistical procedure can be formulated so that it can be used to analyze automated records which do not contain observations at individual levels, only summarized information at categorical levels.
Let Y represent the number of events occurring in the time-exposure interval [0, T] and Z the covariates. We consider the Poisson regression model, (22) where Z is a vector of k + 1 covariates and β a vector of k + 1 regression coefficients. In practice, the exact number of true events may not be available for all subjects. We may instead have the number of possible events, namely, the auxiliary events. For example, in the study of vaccine adverse events associated with childhood immunizations, the number of auxiliary events A for MAARI is collected based on ICD-9 codes through hospital records. Further diagnosis may indicate that some of these events are false events. The number of true vaccine adverse events, Y, can only be confirmed for the subjects in the validation set V. Suppose that Z is the vaccination status, 1 for the vaccinated and 0 for the unvaccinated. Then, under Poisson regression, exp(β) is the relative rate of event occurrence per unit time of the exposed versus unexposed. We assume that the number of automated events A can be expressed as A = Y + W, where W is the number of false events independent of Y conditional on (Z, T). Suppose that W follows the Poisson regression model (23) where γ′ = (a 0 , a 1 , γ 1 , ⋯, γ k−1 ).
We apply the missing data methods introduced in Section 3 on model (22). The variables While the covariate Z can be considered as categorical, it is natural to consider the exposure time T as a continuous variable. We assume that the validation probability depends only on the stratification of (Z, A). That is, the validation sample is a stratified random sample by the categories defined by (Z, A). Of those estimators discussed in Section 2, there are only two different estimators, βÎ 1 and βÂ 2 . We show in Section 4.3 that the proposed method can be formulated so that it can be used to analyze the automated records with missing outcomes. First we derive the explicit estimation procedures for βÎ 1 and βÂ 2 and their variance estimators under model (22).

Inverse probability weighting estimation
We adopt all notations introduced in Section 3. In particular, let and . Let X = (Z, T) and X i = (Z i , T i ) to be consistent with earlier notations. The score function for subject i under model (22) is . The estimator βÎ 1 is obtained by solving , where . By Corollary 1, √n(β̂I 1 − β) converges in distribution to a normal distribution with mean zero and the variance The information matrix can be estimated by Î(β) which is obtained by replacing β with βÎ 1 , P(Z i = z) by the sample proportion of the event {Z i = z}, and E(T i |Z i = z) with the sample average exposure time for those with covariates Z i = z. The matrix Σ A1 (β) can be estimated by (24) where Since A is observed for all subjects, W can be determined if Y is known, and undetermined otherwise. The IPW estimator, γÎ 1 , of γ can be estimated by solving the equation

Augmented inverse probability weighted estimation
Under the assumption that W follows the Poisson regression model (23) and is independent of Y conditional on (Z, T), for a given β by substituting γ by its estimator γÎ 1 of Section 4.1. Then the estimator βÂ 2 is obtained by solving (25) By Corollary 1, √n(β̂A 2 − β) converges in distribution to a normal distribution with mean zero and the variance matrix where I −1 (β) +I −1 (β) Σ A2 (β)I −1 (β), where . The information matrix I(β) can be estimated by Î(β) given in Section 4.1. The conditional variance Var {S β (Y |X)|Z = z, T, A = a} = ap z (1 − p z )z ⊗2 can be estimated by ap̂z (1 − p̂z), where p̂z = exp(β′ z )/(exp(β′z) + exp(γ′ z)). It follows that Σ A2 (β) can be consistently estimated by where ρ(a, z) is the estimator of P{A i = a, Z i = z} and ρ̂v(a, z) is the estimator of P{i ∈ V|A i = a, Z i = z}.

Estimation using the automated data
This section formulates the missing data estimation procedure for (22) based on the automated (summarized) information at categorical levels defined by relevant covariates of the model. In particular, we show that βÎ 1 and βÂ 2 and their variance estimators can be formulated using the automated data at categorical levels.
The following notations are used to show that the estimators of β and their variance estimators can be calculated using the automated information at the categorical levels. Let  V(a, l, m) V(a, l, m), for the number of subjects in , y alm for the number of events for subjects in V(a, l, m), y lm for the number of events for subjects in V(l, m), t alm for the total exposure time for subjects with (A = a, Z ⋍ (l, m)), t 2, alm for the total squared exposure time for subjects with (A = a, Z ⋍ (l, m)), t lm for the total exposure time for subjects with Z ⋍ (l, m), α lm for the number of automated events for subjects with Z ⋍ (l, m).
Estimation with βÎ 1 using the automated data-The validation probability can be estimated by 1/λ alm when A i = a, Z i ⋍ (l, m). It can be shown that the estimating equation for βÎ 1 is equivalent to the following nonlinear equations for {b lm , for l = 0, 1, m = 1, ⋯, k}, for l = 0, 1 and m = 1, …, k − 1. When k > 1, the equations have no explicit solutions.
In the following, we show that the asymptotic variance of βÎ 1 can be consistently estimated by only using the automated information at categorical levels. The information matrix is a (k + 1) × (k + 1) symmetric matrix given by where r = k − 1 and q lm = E(T i e blm I{individual i in category (l, m)}). The consistent estimator, Î(β), of I(β) is thus obtained by replacing q lm with exp(b̂l m )t lm /n. Under model (23), the expected number of false events for a subject in category (l, m) with the time-exposure interval [0, T] is T exp(d lm ), for l = 0, 1 and m = 1, ⋯, k, where d 1k = a 0 + a 1 , d 0k = a 0 , d 1m = a 0 + a 1 + γ m and d 0m = a 0 + γ m for 1 ≤ m ≤ k − 1. The conditional distribution of Y given A = a, T, and Z ⋍ (l, m) is Binomial (a, p lm ), where p lm = exp(b lm )/ (exp(b lm ) + exp(d lm )) for a ≥ 1. Then Var(Y|A = a, Z ⋍ (l, m)) can be estimated by ap̂l m (1 −p lm ), where p lm = e bl m /(e bl m + e dl m ), and Var(T|A = a, Z ⋍ (l, m)) can be estimated by ν a, l, m = t 2, a, l, m /n alm − (t alm /n alm ) 2 . By (24) and the discussion that follows, Σ A1 (β) can be estimated by (26)  where and G lm be the value of when subject i belongs to the category (l, m). Hence the covariance matrix of βÎ 1 can be estimated by Î −1 (β) + Î −1 (β) Σ^A 1 (β)I −1 (β) using the automated data.

Remark
In the special case where ρ(α, l, m) ≈ 0 for α ≥ 2, a much simpler formula for the variance estimator of the log relative risk can be derived. For example in the vaccine safety study, the adverse-event rate is very small. Let Then an estimate of variance of b̂1 is given by (29) which is the weighted sum of the estimated variances for the estimated log relative rate of the exposed versus the unexposed over k groups. The details of deviation are given in the Appendix B.

Qi and Sun
Page 14

A simulation study
We conduct a simulation study to examine the finite sample performance of the IPW estimator βÎ 1 and the AIPW estimator βÂ 2 . We consider the Poisson regression model (22).
The covariates Z 1 and Z 2 are generated from the Bernoulli distributions with the probability of success equals to 0.4 and 0.5 respectively. The exposure time T is generated from a uniform distribution on [0, 10]. Given Z = (Z 1 , Z 2 ) and T, the outcome variable Y follows a Poisson distribution with mean T exp (b 0 + b 1 Z 1 + θZ 2 ) where b 0 = −0.5, b 1 = −0.8 and θ = −0.6, and W follows a Poisson distribution with mean T exp (a 0 + a 1 Z 1 + γZ 2 ) where a 0 = −1.3, a 1 = −1.1, γ = −1. We set A = Y + W.
Four models for the validation sample are considered. Under Model 1, the validation sample is a simple random sample with probability π i = 0.4. Model 2 considers π i = 0.6. In Model 3, the validation probability only depends on A through the logistic regression model logit{π i (X, A)} = A − 0.5 where X = (Z, T). In Model 4, the validation probability depends on A and Z 1 through the logistic regression model logit{π i (X, A)} = A − Z 1 − 0.5.
Tables 1 and 2 present the simulation results for n = 50, 100, 300, 500 and 800. Each entry of the tables is based on 1000 simulation runs. Tables 1 and 2 summarize the bias (Bias), the empirical standard error (SSE), the average of the estimated standard error (ESE), and the empirical coverage probability (CP) of 95% confidence intervals of βÎ 1 and βÂ 2 for β = (b 0 , b 1 , θ). We also compare the performance of the estimators βÎ 1 and βÂ 2 with the completecase (CC) estimator βĈ obtained by simply deleting subjects with missing values of Y i . As a gold standard, we present the estimation results for the full data where all the values of Y i are fully observed. Table 1 presents the results under Model 1 and 2, and Table 2 shows the results under Model 3 and 4. Table 1 shows that under Model 1 and Model 2, the bias of all estimators is very small at a level comparable with that of the full data estimator. The bias decreases with increased sample size and the increased level of the validation probability. The empirical standard errors are in good agreement with the corresponding estimated standard errors, except for the IPW estimator when n ≤ 100 and π ≤ 0.6. Among them, AIPW has the smallest standard errors for all parameters and sample sizes concerned. The coverage probabilities of the confidence intervals for b 0 , b 1 and θ are close to the nominal level 95%. When the sample size and the validation probability are both small, for example, n = 50 and π = 0.4, the IPW has large bias and is unstable but the AIPW still performs well.  Halloran et al. (2003).
Any children representing with history of fever and any respiratory illness were eligible to have a throat swab for influenza virus culture. The decision to obtain specimens was made irrespective of whether a patient had received CAIV-T. The specific case definition was culture-confirmed influenza. Table 3 taken from Halloran et al. (2003) contains information on the number of children in three age groups, the number of children who are vaccinated versus unvaccinated, the number of nonspecific MAARI cases, the number of cultures performed, and the number of cultures positive for each group.
With the method developed in Section 4 for Poisson regression, we compare the risk of developing MAARI for children who received CAIV-T to the risk for children who had never received CAIV-T using the automated information provided in Table 3. The number of nonspecific MAARI cases extracted using the ICD-9 codes is the auxiliary outcome A, whereas the actual number of influenza cases Y is the outcome of interest. Let Z 1 be the treatment indicator (1=vaccine and 0=placebo). Let Z 2 = (η 1 , η 2 ) be the dummy variables indicating three age groups, where η 1 = 1 if the age is in the range 1.5-4, η 1 = 0, otherwise, and η 2 = 1 if the age is in the range 5-9, η 2 = 0, otherwise. The reference group is the age 10-18. The exposure time for all children is taken as T = 1 year.
This data was analyzed by Halloran et al. (2003) and Chu and Halloran (2004). Assuming the binary probability model for P β (Y|X) where X includes the vaccination status and age group indicators, and using the mean score method, Halloran et al. (2003) found that the estimated VE based on the nonspecific MAARI cases alone was 0.18 with 95% confidence interval of (0.11, 0.24). The estimated VE by incorporating the surveillance cultures was 0.79 with 95% confidence interval of (0.51, 0.91). Halloran et al. also reported sample-sizeweighted VE= 0.77 with 95% confidence interval of (0.48, 0.90). Chu and Halloran (2004) have developed a Bayesian method to estimate vaccine efficacy. By Chu and Halloran (2004), the estimated VE was 0.74 with 95% confidence interval (0.50, 0.88) and estimated VE by the multiple imputation method was 0.71 with 95% confidence interval (0.42, 0.86).
Our estimates of the vaccine efficacy are in line with the existing methods. The estimator βÂ 2 yields smaller standard errors and therefore confidence intervals are more precise than the existing methods of Halloran et al. (2003) and Chu and Halloran (2004). Compared to the binary regression, Poisson regression model allows multiple recurrent MAARI cases for each child. Although for this particular application the exposure time is fixed at one year time interval, the proposed method is applicable to the situation where the length of exposure time may be different for different children.

Conclusions
In this paper, we investigated the mean score method, the IPW method and the AIPW method for the parametric probability regression model P β (Y|X) when outcome of interest Y is subject to missingness. The asymptotic distributions are derived for the IPW estimator and the AIPW estimator. The selection probability often needs to be estimated for the IPW estimator, and both the selection probability and the conditional expectation of the score function needs to be estimated for the AIPW estimator. We investigated the properties of the IPW estimator and the AIPW estimator when the selection probability and the conditional expectation are implemented differently.
An AIPW estimator is said to be fully augmented if the selection probability and the conditional expectation are estimated using the full set of observed variables; it is partially augmented if the selection probability and the conditional expectation are estimated using a subset of observed variables. Corollary 1 shows that the fully augmented AIPW estimator is more efficient than the partially augmented AIPW estimator. Corollary 2 shows that the AIPW estimator is more efficient than the IPW estimator. However, when the selection probability depends only on a set of discrete random variables, the IPW estimator obtained by estimating the selection probability nonparametrically with the cell frequencies is asymptotically equivalent to the AIPW estimator augmented using the same set of discrete random variables. Proposition 1 shows that the IPW estimator, the AIPW estimator and the mean score estimator are equivalent if the selection probability and the conditional expectation are estimated using same set of discrete random variables.
Applying the developed missing data methods, we derived the estimation procedures for Poisson regression model with missing outcomes based on auxiliary outcomes and a validated sample for true outcomes. By assuming the selection probability depending only on the observed discrete exposure variables, not on the continuous exposure time, we show that the IPW estimator and the AIPW estimator can be formulated to analyze data when only aggregated/summarized information are available. The simulation study shows that for a moderate sample size and selection probability, the IPW estimator and AIPW estimator perform better than the complete-case estimator. The AIPW estimator is more efficient and more stable than the IPW estimator. The proposed methods are applied to analyze a data set from for an influenza vaccine study conducted in Temple-Belton, Texas during the 2000-2001 influenza season. The data set presented in Table 3 only contains summarized information at categorical levels defined by the three age groups and vaccination status. The actual number of influenza cases (the number of positive cultures) out of the number of MAARI cases cultured, along with the number of MAARI cases, are available for each category. Our analysis using the AIPW approach shows that the age-adjusted relative rate in the vaccinated group compared to the unvaccinated group equals 0.1641, which represents about 84% reduction in the risk of developing MAARI for the vaccinated group compared to the unvaccinated group.
we have . Thus the AIPW estimator βÂ 1 , the IPW estimator βÎ 1 and the mean score estimator βÊ 1 are equivalent to each other.
Note that (A.2) which is not zero unless is linearly related to Z i and in this case β is not identifiable. Hence the AIPW estimator βÂ 2 is different from the AIPW estimator βÂ 1 .

By (A.1) and (A.2), we have
Following the same arguments leading to (A.1), we also have . Hence, the estimators βÎ 2 , βÊ 2 and βÂ 2 are equivalent. By following the steps in (A.2), we also have . Hence, βÂ 3 is the same as βÎ 2 .
Therefore, these are essentially two different estimators.
From ( The second term of (A.3) is (A.4) By (11) Now consider the AIPW estimator βÂ based on solving the estimating equation (14). For Suppose that πĩ and Ẽ i are the estimates of π i and E i based on some parametric or nonparametric models. Then it can be shown using Taylor expansion and standard probability arguments that the second term is at the order of o p (1) under MAR I. Hence It can be shown that under MAR I, and . By routine derivations, we have By the central limit theorem, both n 1/2 (β̂I − β) and n 1/2 (β̂A − β) have asymptotically normal distributions with mean zero and covariances equal to and , respectively.
Next, we examine the covariance matrices and to understand the efficiency gain of βÂ over βÎ. Note that and .
Denote The second term of equals Σ A1 (β) and the second term of equals Σ A2 (β). Then it follows from the main results in Theorem 1 that (17) and (18) hold.
Also by Theorem 1, the difference in the variances of and contributes to the difference in the asymptotic variances of βÂ 1 and βÂ 2 . Since under MAR I, which is less than Σ A1 (β) if the covariates Z i is a proper subset of X i .

Proof of Corollary 2
Consider the definitions of B i and O i given following (16) From the discussions preceding Corollary 2, ѱ = {ѱ z, a } and π(z, a, ѱ) = ѱ z, a , where ѱ z, a = P(ξ i = 1|Z i = z, A i = a) for all distinct pairs (z, a). Hence, is a column vector with 1 in the position for ѱ z, a and 0 elsewhere. And I ѱ is a diagonal matrix and its inverse matrix is also diagonal.