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

Abstract

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. Mathematics Subject Classification Primary 62J02; Secondary 62F12

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 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 $\bar{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. Let (Z _i ,X _i ,Y _i ,A _i ), i=1,…,n, be independent identically distributed (iid) copies of (Z,X,Y,A). Let ξ _i =I(i∈V) be the selection indicator.

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.

MAR I:ξ _i is independent of Y _i conditional on (X _i ,A _i ) and ξ _i is independent of $Z_i^c$ conditional on (Z _i ,A _i ).

MAR II:ξ _i is independent of $(Y_i,Z_i^c)$ conditional on (Z _i ,A _i ).

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 ${\widehat{\pi }}_i$ be the estimator of the conditional probability π _i =P(ξ _i =1|X _i ,A _i ), and ${\widehat{\pi }}_i^z$ the estimator of ${\pi }_i^z=P({\xi }_i=1|Z_i,A_i)$. Let S _β (Y|X) denote the partial derivatives of logP _β (Y|X) with respect to β. Let $Ê\left\{S_{\beta }\left(Y\right|X_i\left)\right|X_i,A_i\right\}$ be the estimator of the conditional expectation E{S _β (Y|X _i )|X _i ,A _i }, and $Ê\left\{S_{\beta }\left(Y\right|X_i\left)\right|Z_i,A_i\right\}$ the estimator of E{S _β (Y|X _i )|Z _i ,A _i }. We investigate several estimators of β based on the following estimating equations with different choices of W _i :

$$\sum _{i=1}^n{W}_i=0,$$

(1)

where W _i takes one of the following forms:

$$W_i^{I1}=\frac{{\xi }_i}{{\widehat{\pi }}_i^z}{S}_{\beta }\left(Y_i\right|X_i)$$

(2)

$$W_i^{E1}={\xi }_i{S}_{\beta }\left(Y_i\right|X_i)+(1-{\xi }_i)Ê\left\{S_{\beta }\left(Y\right|X_i\left)\right|Z_i,A_i\right\}$$

(3)

$$W_i^{A1}=\frac{{\xi }_i}{{\widehat{\pi }}_i^z}{S}_{\beta }\left(Y_i\right|X_i)+\left(1-\frac{{\xi }_i}{{\widehat{\pi }}_i^z}\right)Ê\left\{S_{\beta }\left(Y\right|X_i\left)\right|Z_i,A_i\right\}$$

(4)

$$W_i^{I2}=\frac{{\xi }_i}{{\widehat{\pi }}_i}{S}_{\beta }\left(Y_i\right|X_i)$$

(5)

$$W_i^{E2}={\xi }_i{S}_{\beta }\left(Y_i\right|X_i)+(1-{\xi }_i)Ê\left\{S_{\beta }\left(Y\right|X_i\left)\right|X_i,A_i\right\}$$

(6)

$$W_i^{A2}=\frac{{\xi }_i}{{\widehat{\pi }}_i^z}{S}_{\beta }\left(Y_i\right|X_i)+\left(1-\frac{{\xi }_i}{{\widehat{\pi }}_i^z}\right)Ê\left\{S_{\beta }\left(Y\right|X_i\left)\right|X_i,A_i\right\}.$$

(7)

$$W_i^{A3}=\frac{{\xi }_i}{{\widehat{\pi }}_i}{S}_{\beta }\left(Y_i\right|X_i)+\left(1-\frac{{\xi }_i}{{\widehat{\pi }}_i}\right)Ê\left\{S_{\beta }\left(Y\right|X_i\left)\right|X_i,A_i\right\}.$$

(8)

The estimator ${\widehat{\beta }}_{I1}$ obtained by using $W_i^{I1}$ is an IPW estimator where a subject’s validation probability ${\pi }_i^z$ depends only on the category defined by (Z _i ,A _i ). Because $E\left\{{\left({\pi }_i^z\right)}^{-1}{\xi }_i{S}_{\beta }\left(Y_i\right|X_i)\right\}=E\left\{S_{\beta }\left(Y_i\right|X_i)\right\}=0$, the estimator ${\widehat{\beta }}_{I1}$ is approximately unbiased. The estimator ${\widehat{\beta }}_{I2}$ obtained by using $W_i^{I2}$ 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 $Z_i^c$.

The estimator ${\widehat{\beta }}_{E1}$ obtained by using $W_i^{E1}$ 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 ${\widehat{\beta }}_{E2}$ obtained by using $W_i^{E2}$ 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 ${\widehat{\beta }}_{E2}$ is the mean score estimator in Pepe et al. ([1994]). The mean score estimator is the MLE estimator employing the EM-algorithm ([Dempster et al. 1977]) under the assumption that the auxiliary outcome is noninformative in the sense that the probability model P _θ (A|Y,X) is unrelated to β.

The estimator ${\widehat{\beta }}_{A1}$ obtained using $W_i^{A1}$ is the AIPW estimator augmented with the estimated conditional expectation $Ê\left\{S_{\beta }\left(Y\right|X_i\left)\right|Z_i,A_i\right\}$. The estimator ${\widehat{\beta }}_{A2}$ obtained using $W_i^{A2}$ is the AIPW estimator augmented with the estimated conditional expectation $Ê\left\{S_{\beta }\left(Y\right|X_i\left)\right|X_i,A_i\right\}$. The estimator ${\widehat{\beta }}_{A3}$ is obtained using $W_i^{A3}$. The $W_i^{A3}$ differs from $W_i^{A2}$ in that the estimated validation probability is ${\widehat{\pi }}_i$ instead of ${\widehat{\pi }}_i^z$.

Suppose that ${\widehat{\pi }}_i^z$ is an asymptotically unbiased estimator of ${\bar{\pi }}_i^z$ and that $Ê\left\{S_{\beta }\left(Y\right|X_i\left)\right|Z_i,A_i\right\}$ is asymptotically unbiased of $Ē\left\{S_{\beta }\left(Y\right|X_i\left)\right|Z_i,A_i\right\}$, where both ${\bar{\pi }}_i^z$ and $Ē\left\{S_{\beta }\left(Y\right|X_i\left)\right|Z_i,A_i\right\}$ are functions of (Z _i ,A _i ). Under MAR II, if one of the equalities, ${\bar{\pi }}_i^z={\pi }_i^z$ and $\left.Ē\left[S_{\beta }\left(Y\right|X_i\left)\right|Z_i,A_i\right]\right\}\left.=E\left[S_{\beta }\left(Y\right|X_i\left)\right|Z_i,A_i\right]\right\}$, holds, then

$$E\left\{{\left({\bar{\pi }}_i^z\right)}^{-1}{\xi }_i{S}_{\beta }\left(Y_i|X_i\right)\right\}+E\left\{\left(1-{\left({\pi }_i^z\right)}^{-1}{\xi }_i\right)Ē\left[S_{\beta }\left(Y|X_i\right)|Z_i,A_i\right]\right\}=E\left\{S_{\beta }\left(Y_i\right|X_i)\right\}=0,$$

which entails that the estimator ${\widehat{\beta }}_{A1}$ has the double robust property in the sense that it is a consistent estimator of β if either ${\widehat{\pi }}_i^z$ is a consistent estimator of ${\pi }_i^z$ or $Ê\left[S_{\beta }\left(Y\right|X_i\left)\right|Z_i,A_i\right]$ is a consistent estimator of E[S _β (Y|X _i )|Z _i ,A _i ]. Similarly, under MAR I, the estimator ${\widehat{\beta }}_{A2}$ possesses the double robust property in that ${\widehat{\beta }}_{A2}$ is a consistent estimator of β if either ${\widehat{\pi }}_i^z$ is a consistent estimator of ${\pi }_i^z$ or $Ê\left[S_{\beta }\left(Y\right|X_i\left)\right|X_i,A_i\right]$ is a consistent estimator of E[S _β (Y|X _i )|X _i ,A _i ]. The estimator ${\widehat{\beta }}_{A3}$ has similar double robust property as ${\widehat{\beta }}_{A2}$.

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 ${\widehat{\pi }}_i=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,

$$Ê\left\{S_{\beta }\left(Y\right|X_i\left)\right|X_i,A_i\right\}=\sum _{j\in V(X_i,A_i)}{S}_{\beta }\left(Y_j|X_j\right)/n^V(X_i,A_i),$$

(9)

Under MAR I, $Ê\left\{S_{\beta }\left(Y\right|X_i\left)\right|X_i,A_i\right\}$ is an unbiased estimator of E{S _β (Y|X _i )|X _i ,A _i }. Now we let V(Z _i ,A _i ) denote the subjects in V with values of (Z,A) equal to (Z _i ,A _i ), n^V(Z _i ,A _i ) the number of subjects in V(Z _i ,A _i ), and n(Z _i ,A _i ) the number of subjects in the sample with values of (Z,A) equal to (Z _i ,A _i ). A nonparametric estimator of ${\pi }_i^z=P({\xi }_i=1|Z_i,A_i)$ is given by ${\widehat{\pi }}_i^z=n^V(Z_i,A_i)/n(Z_i,A_i)$. A nonparametric estimator of E{S _β (Y|X _i )|Z _i ,A _i } is given by

$$Ê\left\{S_{\beta }\left(Y\right|X_i\left)\right|Z_i,A_i\right\}=\sum _{j\in V(Z_i,A_i)}{S}_{\beta }\left(Y_j|X_j\right)/n^V(Z_i,A_i).$$

(10)

Under MAR II, (Y _i ,X _i ) is independent of ξ _i conditional on (Z _i ,A _i ), then $Ê\left\{S_{\beta }\left(Y\right|X_i\left)\right|Z_i,A_i\right\}$ is an unbiased estimator of E{S _β (Y|X _i )|Z _i ,A _i }.

Proposition 1.

Suppose that X=(Z,Z^c) and A are discrete and their dimensionality is reasonably small. Under the nonparametric estimators ${\widehat{\pi }}_i^z=n^V(Z_i,A_i)/n(Z_i,A_i)$, ${\widehat{\pi }}_i=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 ${\widehat{\beta }}_{I1}$, ${\widehat{\beta }}_{E1}$ and ${\widehat{\beta }}_{A1}$ are equivalent, and the estimators ${\widehat{\beta }}_{I2}$, ${\widehat{\beta }}_{E2}$, ${\widehat{\beta }}_{A2}$ and ${\widehat{\beta }}_{A3}$ are equivalent. However, the estimator ${\widehat{\beta }}_{A2}$ is different from ${\widehat{\beta }}_{A1}$ unless $Z_i^c$ 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 ${\widehat{\beta }}_{A2}$ and ${\widehat{\beta }}_{A3}$ 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.

Suppose that $\tilde{E}\left\{S_{\beta }\left(Y\right|X_i\left)\right|X_i,A_i\right\}$ is a consistent parametric/nonparametric estimator of E _a {S _β (Y|X _i )|X _i ,A _i }, where E _a {S _β (Y|X _i )|X _i ,A _i } is E{S _β (Y|X _i )|X _i ,A _i } or E{S _β (Y|X _i )|Z _i ,A _i }. Let π(X _i ,A _i ,ψ) be the parametric model for the validation probability π _i , where ψ is a q-dimensional 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 ψ₀ be the true value of ψ. Under MAR I, the MLE $\widehat{\psi }=\left({\widehat{\psi }}_1,\dots ,{\widehat{\psi }}_q\right)$ of ψ=(ψ₁,…,ψ _q ) is obtained by maximizing the observed data likelihood,

$$\prod _{i=1}^n{\left\{\pi (X_i,A_i,\psi )\right\}}^{{\xi }_i}{\left\{1-\pi (X_i,A_i,\psi )\right\}}^{1-{\xi }_i}.$$

The validation probability π _i is estimated by ${\tilde{\pi }}_i=\pi \left(X_i,A_i,\widehat{\psi }\right)$. Then by the standard likelihood based analysis, we have the approximation

$$\widehat{\psi }-{\psi }_0=n^{-1}\sum _{i=1}^n{\left(I^{\psi }\right)}^{-1}{S}_i^{\psi }+o_p\left(n^{-1/2}\right),$$

(11)

where $S_i^{\psi }$ and I^ψ are the score vector and information matrix for $\widehat{\psi }$ defined by

$$\begin{array}{lcr}{S}_i^{\psi }& =& \frac{({\xi }_i-\pi (X_i,A_i,{\psi }_0\left)\right)}{\pi (X_i,A_i,{\psi }_0)(1-\pi (X_i,A_i,{\psi }_0\left)\right)}\frac{\mathrm{\partial \pi }(X_i,A_i,{\psi }_0)}{\mathrm{\partial \psi }},\\ I^{\psi }& =& E\left\{\frac{1}{\pi (X_i,A_i,{\psi }_0)(1-\pi (X_i,A_i,{\psi }_0\left)\right)}{\left(\frac{\mathrm{\partial \pi }(X_i,A_i,{\psi }_0)}{\mathrm{\partial \psi }}\right)}^{\otimes 2}\right\},\end{array}$$

(12)

where a^⊗2=a a^′.

Consider the IPW estimator ${\widehat{\beta }}_I$ obtained by solving the estimating equation

$$U_I=\sum _{i=1}^n\frac{{\xi }_i}{{\tilde{\pi }}_i}{S}_{\beta }\left(Y_i\right|X_i)$$

(13)

and the AIPW estimator ${\widehat{\beta }}_A$ based on solving the estimating equation

$$U_A=\sum _{i=1}^n\left[\frac{{\xi }_i}{{\tilde{\pi }}_i}{S}_{\beta }\left(Y_i|X_i\right)+\left(1-\frac{{\xi }_i}{{\tilde{\pi }}_i}\right)\tilde{E}\left\{S_{\beta }\left(Y\right|X_i\left)\right|X_i,A_i\right\}\right].$$

(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{(∂²/∂ β²)(logP _β (Y|X))} and I^ψ are positive definite at the true parameters. Then, under MAR I,

$$\begin{array}{l}{n}^{1/2}\left({\widehat{\beta }}_I-\beta \right)=I^{-1}\left(\beta \right)n^{-1/2}\sum _{i=1}^n{Q}_i^I+o_p\left(1\right),\\ n^{1/2}\left({\widehat{\beta }}_A-\beta \right)=I^{-1}\left(\beta \right)n^{-1/2}\sum _{i=1}^n{Q}_i^A+o_p\left(1\right),\end{array}$$

where I(β)=E{−(∂²/∂ β²)(logP _β (Y|X))}=Var(S _β (Y _i |X _i )),

$$Q_i^I={\xi }_i/{\pi }_i{S}_{\beta }\left(Y_i\right|X_i)-E\left\{{{\pi }_i}^{-2}{\xi }_i{S}_{\beta }\left(Y_i\right|X_i){\left(\mathrm{\partial \pi }(X_i,A_i,{\psi }_0)/\mathrm{\partial \psi }\right)}^{\prime }\right\}{\left(I^{\psi }\right)}^{-1}{S}_i^{\psi }$$

and $Q_i^A={\xi }_i/{\pi }_i{S}_{\beta }\left(Y_i\right|X_i)+\left(1-{\xi }_i/{\pi }_i\right)E_a\left\{S_{\beta }\left(Y\right|X_i\left)\right|X_i,A_i\right\}$.

Both $n^{1/2}\left({\widehat{\beta }}_I-\beta \right)$ and $n^{1/2}\left({\widehat{\beta }}_A-\beta \right)$ have asymptotically normal distributions with mean zero and covariances equal to $I^{-1}\left(\beta \right)\text{Var}\left(Q_i^I\right)I^{-1}\left(\beta \right)$ and $I^{-1}\left(\beta \right)\text{Var}\left(Q_i^A\right)I^{-1}\left(\beta \right)$, respectively. Further,

$$\begin{array}{l}\text{Var}\left(Q_i^I\right)=\text{Var}\left(Q_i^A\right)+\text{Var}(B_i+O_i)\end{array}$$

(15)

and

$$\begin{array}{l}\text{Var}\left(Q_i^A\right)=I\left(\beta \right)+\text{Var}\left(\left(1-\frac{{\xi }_i}{{\pi }_i}\right)\left\{S_{\beta }\left(Y_i\right|X_i)-E_a\left\{S_{\beta }\left(Y\right|X_i\left)\right|X_i,A_i\right\}\right\}\right),\end{array}$$

(16)

where $O_i=E\left\{{{\pi }_i}^{-2}{\xi }_i{S}_{\beta }\left(Y_i\right|X_i){\left(\mathrm{\partial \pi }(X_i,A_i,{\psi }_0)/\mathrm{\partial \psi }\right)}^{\prime }\right\}{\left(I^{\psi }\right)}^{-1}{S}_i^{\psi }$ and B _i =(1−ξ _i /π _i )E _a {S _β (Y|X _i )|X _i ,A _i }.

Suppose that the validation probability π _i = P(ξ _i = 1|X _i ,A _i ) depends only on (Z _i ,A _i ). That is, ${\pi }_i={\pi }_i^z=P\left({\xi }_i=1|Z_i,A_i\right)$. Suppose that ${\tilde{\pi }}_i$ is the MLE of ${\pi }_i^z$ under the parametric family ψ(Z _i ,A _i ,ψ). Let ${\widehat{\beta }}_{A1}$ be the estimator obtained by solving (14) where the augmented term, $\tilde{E}\left\{S_{\beta }\left(Y\right|X_i\left)\right|X_i,A_i\right\}$, is a consistent parametric/nonparametric estimator of E{S _β (Y|X _i )|Z _i ,A _i }. Let ${\widehat{\beta }}_{A2}$ be the estimator obtained by solving (14) where $\tilde{E}\left\{S_{\beta }\left(Y\right|X_i\left)\right|X_i,A_i\right\}$ is a consistent parametric/nonparametric estimator of E{S _β (Y|X _i )|X _i ,A _i }. 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 π _i =P(ξ _i =1|X _i ,A _i ) depends only on (Z _i ,A _i ). Under the conditions of Theorem 1,

$$n^{1/2}\left({\widehat{\beta }}_{A1}-\beta \right)\stackrel{\mathcal{D}}{\to }N\left(0,I^{-1}\left(\beta \right)+I^{-1}\left(\beta \right){\Sigma }_{A1}\left(\beta \right)I^{-1}\left(\beta \right)\right),$$

(17)

and

$$n^{1/2}\left({\widehat{\beta }}_{A2}-\beta \right)\stackrel{\mathcal{D}}{\to }N\left(0,I^{-1}\left(\beta \right)+I^{-1}\left(\beta \right){\Sigma }_{A2}\left(\beta \right)I^{-1}\left(\beta \right)\right),$$

(18)

where ${\Sigma }_{A1}\left(\beta \right)=E\left[\left(\right(1-{\pi }_i^z)/{\pi }_i^z)\text{Var}\left\{S_{\beta }\right(Y_i\left|X_i\right)|Z_i,A_i\}\right]$ and ${\Sigma }_{A2}\left(\beta \right)=E\left[\left(\right(1-{\pi }_i^z)/{\pi }_i^z)\text{Var}\left\{S_{\beta }\left(Y_i\right|X_i\left)\right|X_i,A_i\right\}\right]$. The asymptotic variance of ${\widehat{\beta }}_{A2}$ is smaller than the asymptotic variance of ${\widehat{\beta }}_{A1}$ 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 π(z,a,ψ)=ψ_z,a. Further, $\frac{\mathrm{\partial \pi }(z,a,{\psi }_0)}{\mathrm{\partial \psi }}$ 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,

$$\begin{array}{l}{I}^{\psi }=\sum _{z,a}\rho (z,a)\left\{\frac{1}{\pi (z,a,{\psi }_0)(1-\pi (z,a,{\psi }_0\left)\right)}\frac{\mathrm{\partial \pi }(z,a,{\psi }_0)}{\mathrm{\partial \psi }}{\left(\frac{\mathrm{\partial \pi }(z,a,{\psi }_0)}{\mathrm{\partial \psi }}\right)}^{\prime }\right\},\end{array}$$

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 ${\widehat{\psi }}_{z,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

$${\widehat{\psi }}_{z,a}-\pi (z,a,{\psi }_0)=n^{-1}\frac{n^V(z,a)-n(z,a)\pi (z,a,{\psi }_0)}{\rho (z,a)}+o_p\left(n^{-1/2}\right),$$

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

$$\begin{array}{lcr}{n}^{-1/2}{U}_A& =& n^{-1/2}\sum _{i=1}^n{Q}_i^A+o_p\left(1\right)\end{array}$$

(19)

$$\begin{array}{lcr}{n}^{-1/2}{U}_I& =& n^{-1/2}\sum _{i=1}^n{Q}_i^A-n^{-1/2}\sum _{i=1}^n(B_i+O_i)+o_p\left(1\right),\end{array}$$

(20)

where B _i and O _i are defined following (16). The following corollary presents the analysis of the term $n^{-1/2}\sum _{i=1}^n(B_i+O_i)$ when (Z _i ,A _i ) are discrete to understand how efficiency may be gained from the AIPW estimator over the IPW estimator.

Corollary 2.

Under the conditions of Theorem 1, suppose that X=(Z,Z^c) and (Z,A) are discrete taking values (z,a) in a set of finite number of values. Suppose that the validation probability π _i =P(ξ _i =1|X _i ,A _i ) only depends on (Z _i ,A _i ) and ψ_z,a=P(ξ _i =1|Z _i =z,A _i =a) is estimated nonparametrically by ${\widehat{\psi }}_{z,a}=n^V(z,a)/n(z,a)$. Then

$$\begin{array}{l}{n}^{-1/2}\sum _{i=1}^n(O_i+B_i)\\ =-\sum _{z,a}{n}^{-1/2}\sum _{j=1}^n\frac{{\xi }_j-\pi (z,a,{\psi }_0)}{\pi (z,a,{\psi }_0)}I(Z_j=z,A_j=a)\\ \left[E\left\{S_{\beta }\right(Y_j\left|X_j\right)|Z_j=z,A_j=a\}-E\left\{S_{\beta }\left(Y_j\right|X_j\left)\right|Z_j=z,Z_j^c,A_j=a\right\}\right].\end{array}$$

(21)

By Corollary 2, (19) and (20), ${\widehat{\beta }}_A$ is more efficient than ${\widehat{\beta }}_I$ unless V a r{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 ${\widehat{\psi }}_{z,a}=n^V(z,a)/n(z,a)$, then ${\widehat{\beta }}_A$ and ${\widehat{\beta }}_I$ 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 ${\pi }_i^z=P({\xi }_i=1|Z_i,A_i)$ is estimated nonparametrically by the cell frequency, then by Theorem 1 and Corollary 2, ${\widehat{\beta }}_{A1}$ and ${\widehat{\beta }}_{I1}$ have same asymptotic normal distributions as long as $Ê\left[S_{\beta }\right(Y\left|X_i\right)|Z_i,A_i]$ is a consistent estimator of E[S _β (Y|X _i )|Z _i ,A _i ]. But ${\widehat{\beta }}_{A2}$ is more efficient than ${\widehat{\beta }}_{I1}$ as long as $Ê\left[S_{\beta }\right(Y\left|X_i\right)|X_i,A_i]$ is a consistent estimator of E[S _β (Y|X _i )|X _i ,A _i ] since Var(B _i +O _i ) is not zero by (21). These results are not affected by whether E[S _β (Y|X _i )|Z _i ,A _i ] and E[S _β (Y|X _i )|X _i ,A _i ] are estimated nonparametrically or based on some parametric models. In addition, by Theorem 1, Corollary 1 and 2, ${\widehat{\beta }}_{A3}$ and ${\widehat{\beta }}_{I2}$ have the same asymptotic normal distributions as long as $Ê\left[S_{\beta }\right(Y\left|X_i\right)|X_i,A_i]$ is a consistent estimator of E[S _β (Y|X _i )|X _i ,A _i ].

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,

$$P(Y=y|Z,T)=exp\left\{-Texp\left({\beta }^{\prime }Z\right)\right\}{\left\{Texp\left({\beta }^{\prime }Z\right)\right\}}^y/y!,$$

(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

$$P(W=w|Z,T)=exp\left\{-Texp\left({\gamma }^{\prime }Z\right)\right\}{\left\{Texp\left({\gamma }^{\prime }Z\right)\right\}}^w/w!,$$

(23)

where γ^′=(a₀,a₁,γ₁,⋯,γ_k−1).

We apply the missing data methods introduced in Section 3 on model (22). The variables (Z _i ,T _i ,Y _i ,A _i ) are observed for the validation sample V and (Z _i ,T _i ,A _i ) are observed for the nonvalidation sample $\bar{V}$. 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, ${\widehat{\beta }}_{I1}$ and ${\widehat{\beta }}_{A2}$. 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 ${\widehat{\beta }}_{I1}$ and ${\widehat{\beta }}_{A2}$ and their variance estimators under model (22).

4.1 Inverse probability weighting estimation

We adopt all notations introduced in Section 3. In particular, let ${\pi }_i^z=P({\xi }_i=1|Z_i,A_i)$ and ${\widehat{\pi }}_i^z=n^V(Z_i,A_i)/n(Z_i,A_i)$. 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 $S_{\beta }\left(Y_i\right|X_i)=Z_i^{\prime }(Y_i-T_{i}exp\left({\beta }^{\prime }{Z}_i\right)).$ The estimator ${\widehat{\beta }}_{I1}$ is obtained by solving $\sum _{i=1}^n({\xi }_i/{\widehat{\pi }}_i^z)S_{\beta }\left(Y_i\right|X_i)=0,$ where $S_{\beta }\left(Y_i\right|X_i)=Z_i^{\prime }(Y_i-T_{i}exp\left({\gamma }^{\prime }{Z}_i\right))$. By Corollary 1, $\sqrt{n}({\widehat{\beta }}_{I1}-\beta )$ converges in distribution to a normal distribution with mean zero and the variance matrix I⁻¹(β)+I⁻¹(β)Σ_{A 1}(β)I⁻¹(β), where ${\Sigma }_{A1}\left(\beta \right)=E\left[\left(\right(1-{\pi }_i^z)/{\pi }_i^z)\text{Var}\left\{S_{\beta }\left(Y_i\right|X_i\left)\right|Z_i,A_i\right\}\right]$.

The information matrix $I\left(\beta \right)=E(Z_i{Z}_i^{\prime }{T}_{i}exp({\beta }^{\prime }{Z}_i\left)\right)=\sum _{z}P(Z_i=z)z{z}^{\prime }exp\left({\beta }^{\prime }z\right)E\left(T_i\right|Z_i=z)$ can be estimated by $Î\left(\widehat{\beta }\right)$ which is obtained by replacing β with ${\widehat{\beta }}_{I1}$, 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 Σ_{A 1}(β) can be estimated by

$${\widehat{\Sigma }}_{A1}\left(\widehat{\beta }\right)=\sum _{a,z}\widehat{\rho }(a,z)\frac{1-{\widehat{\rho }}^v(a,z)}{{\widehat{\rho }}^v(a,z)}\widehat{\text{Var}}\left\{S_{\beta }\right(Y\left|X\right)|A=a,Z=z\},$$

(24)

where $\widehat{\rho }(a,z)$ is the estimator of P{A _i =a,Z _i =z}, ${\widehat{\rho }}^v(a,z)$ is the estimator of P{i∈V|A _i =a,Z _i =z}, and $\widehat{\text{Var}}\left\{S_{\beta }\left(Y\right|X\left)\right|A=a,Z=z\right\}$ is an estimator of Var{S _β (Y _i |X _i )|Z _i ,A _i }] which is derived in the following.

Since A is observed for all subjects, W can be determined if Y is known, and undetermined otherwise. The IPW estimator, ${\widehat{\gamma }}_{I1}$, of γ can be estimated by solving the equation $\sum _{i=1}^n({\xi }_i/{\widehat{\pi }}_i^z)S_{\gamma }\left(W_i\right|X_i)=0,$ where $S_{\gamma }\left(W_i\right|X_i)=Z_i^{\prime }(W_i-T_{i}exp\left({\gamma }^{\prime }{Z}_i\right))$. The conditional distribution of Y given A=a, T, and Z=z is Binomial (a,p _z ), where p _z = exp(β^′z)/(exp(β^′z)+ exp(γ^′z)). Since this conditional distribution does not depend on T, the outcome Y and T are conditionally independent given (A,Z). Therefore, Var{S _β (Y|X)|A,Z}=Z Z^′{Var(Y|A,Z)+ exp(2β^′Z)Var(T|A,Z)}. The variance Var(Y|A=a,Z=z) can be estimated by $a{\widehat{p}}_z(1-{\widehat{p}}_z)$, where ${\widehat{p}}_z=exp\left({\widehat{\beta }}^{\prime }z\right)/\left\{exp\left({\widehat{\beta }}^{\prime }z\right)+exp\left({\widehat{\gamma }}^{\prime }z\right)\right\}$, and Var(T|A=a,Z=z)=E(T²|A=a,Z=z)−{E(T|A=a,Z=z)}² can be estimated nonparametrically using the first and the second sample moments conditional on each category with A=a and Z=z.

4.2 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), $E\left\{S_{\beta }\left(Y\right|X\left)\right|Z,T,A\right\}=A{Z}^{\prime }\frac{exp\left({\beta }^{\prime }Z\right)}{exp\left({\beta }^{\prime }Z\right)+exp\left({\gamma }^{\prime }Z\right)}-T{Z}^{\prime }exp\left({\beta }^{\prime }Z\right).$ Let $Ê\left\{S_{\beta }\left(Y\right|X_i\left)\right|X_i,A_i\right\}$ be the estimator of E{S _β (Y|X _i )|X _i ,A _i } for a given β by substituting γ by its estimator ${\widehat{\gamma }}_{I1}$ of Section 4.1. Then the estimator ${\widehat{\beta }}_{A2}$ is obtained by solving

$$\sum _{i=1}^n\frac{{\xi }_i}{{\widehat{\pi }}_i^z}{S}_{\beta }\left(Y_i|X_i\right)+\left(1-\frac{{\xi }_i}{{\widehat{\pi }}_i^z}\right)Ê\left\{S_{\beta }\left(Y\right|X_i\left)\right|X_i,A_i\right\}=0.$$

(25)

By Corollary 1, $\sqrt{n}({\widehat{\beta }}_{A2}-\beta )$ converges in distribution to a normal distribution with mean zero and the variance matrix where I⁻¹(β)+I⁻¹(β)Σ_{A 2}(β)I⁻¹(β), where ${\Sigma }_{A2}\left(\beta \right)=E\left[\left(\right(1-{\pi }_i^z)/{\pi }_i^z)\text{Var}\left\{S_{\beta }\right(Y_i\left|X_i\right)|X_i,A_i\}\right]$. The information matrix I(β) can be estimated by $Î\left(\widehat{\beta }\right)$ given in Section 4.1. The conditional variance Var{S _β (Y|X)|Z=z,T,A=a}=a p _z (1−p _z )z^⊗2 can be estimated by $a{\widehat{p}}_z(1-{\widehat{p}}_z)$, where ${\widehat{p}}_z=exp\left({\widehat{\beta }}^{\prime }z\right)/(exp({\widehat{\beta }}^{\prime }z)+exp({\widehat{\gamma }}^{\prime }z\left)\right)$. It follows that Σ_{A 2}(β) can be consistently estimated by

$${\widehat{\Sigma }}_{A2}\left(\beta \right)=\sum _{a,z}\widehat{\rho }(a,z)\frac{1-{\widehat{\rho }}^v(a,z)}{{\widehat{\rho }}^v(a,z)}a{\widehat{p}}_z(1-{\widehat{p}}_z)z^{\otimes 2},$$

where $\widehat{\rho }(a,z)$ is the estimator of P{A _i =a,Z _i =z} and ${\widehat{\rho }}^v(a,z)$ is the estimator of P{i∈V|A _i =a,Z _i =z}.

4.3 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 ${\widehat{\beta }}_{I1}$ and ${\widehat{\beta }}_{A2}$ and their variance estimators can be formulated using the automated data at categorical levels.

In many applications it is convenient to write Z=(1,Z₍₁₎,Z₍₂₎) and β=(b₀,b₁,θ^′)^′, where Z₍₁₎ is the treatment indicator (Z₍₁₎=1 for the exposed group and Z₍₁₎=0 for the unexposed group) and Z₍₂₎=(η₁,⋯,η_k−1)^′ as the other covariates, and θ=(θ₁,⋯,θ_k−1)^′. For the applications involving the automated data records, we let η₁,⋯,η_k−1 be k−1 dummy variables representing k groups. Without loss of generality, we choose the k th group as the reference group, η₁=1, η₂=0, ⋯, η_k−1=0 for group 1, η₁=0, η₂=1, ⋯, η_k−1=0 for group 2, so on and η₁=0, η₂=0, ⋯, η_k−1=0 for group k. Thus each value of Z denotes a category which can be represented by (l,m) for l=0,1 and m=1,⋯,k. This correspondence is denoted by Z≃(l,m) for convenience. For l=0,1 and m=1,⋯,k−1, category (l,m) is defined by Z with Z₍₁₎=l, η _m =1 and η _j =0 for j≠m, j=1,…,k, and category (l,k) is defined by Z₍₁₎=l and η _j =0 for j=1,…,k−1. Under model (22), the expected number of events for a subject in category (l,m) with the time-exposure interval [ 0,T] is T exp(b _lm ), for l=0,1 and m=1,⋯,k, where b_1k=b₀+b₁, b_0k=b₀, b_1m=b₀+b₁+θ _m and b_0m=b₀+θ _m for 1≤m≤k−1. The parameter b₁ represents the log-relative rate of the exposed versus the unexposed adjusted for other factors.

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) denote the set of subjects in V with (A=a, Z≃(l,m)), V(l,m) for the set of subjects in V with (Z≃(l,m)), n _alm for the number of subjects with (A=a, Z≃(l,m)), $n_{\mathit{\text{alm}}}^v$ for the number of subjects in V(a,l,m), $n_{\mathit{\text{lm}}}^v$ for the number of subjects in V(l,m), ${\lambda }_{\mathit{\text{alm}}}=n_{\mathit{\text{alm}}}/n_{\mathit{\text{alm}}}^v$, 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,a l m} 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${\widehat{\beta }}_{I1}$using the automated data. The validation probability ${\pi }_i^z$ can be estimated by 1/λ _alm when A _i =a, Z _i ≃(l,m). It can be shown that the estimating equation for ${\widehat{\beta }}_{I1}$ is equivalent to the following nonlinear equations for {b _lm , for l=0,1, m=1,⋯,k},

$$\begin{array}{l}\sum _{m=1}^k\left(\sum _{a\in A}{y}_{\mathit{\text{alm}}}{\lambda }_{\mathit{\text{alm}}}-e^{b_{\mathit{\text{lm}}}}\sum _{a\in A}{t}_{\mathit{\text{alm}}}{\lambda }_{\mathit{\text{alm}}}\right)=0,\\ \sum _{l=0,1}\left(\sum _{a\in A}{y}_{\mathit{\text{alm}}}{\lambda }_{\mathit{\text{alm}}}-e^{b_{\mathit{\text{lm}}}}\sum _{a\in A}{t}_{\mathit{\text{alm}}}{\lambda }_{\mathit{\text{alm}}}\right)=0,\end{array}$$

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 ${\widehat{\beta }}_{I1}$ 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

$$\begin{array}{l}I\left(\beta \right)=E(Z_i{Z}_i^{\prime }{T}_{i}exp({\beta }^{\prime }{Z}_i\left)\right)\\ =\left[\begin{array}{ccccc}\sum _{l,m}{q}_{\mathit{\text{lm}}}& \sum _{m=1}^k{q}_{1m}& q_{11}+q_{01}& \cdots & q_{1r}+q_{0r}\\ \sum _{m=1}^k{q}_{1m}& \sum _{m=1}^k{q}_{1m}& q_{11}& \cdots & q_{1r}\\ q_{11}+q_{01}& q_{11}& q_{11}+q_{01}& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ q_{1r}+q_{0r}& q_{1r}& 0& \cdots & q_{1r}+q_{0r}\end{array}\right],\end{array}$$

where r=k−1 and $q_{\mathit{\text{lm}}}=E\left(T_i{e}^{b_{\mathit{\text{lm}}}}I\right\{\text{individual}i\text{in category}(l,m)\left\}\right)$. The consistent estimator, $Î\left(\widehat{\beta }\right)$, of I(β) is thus obtained by replacing q _lm with $exp\left({\widehat{b}}_{\mathit{\text{lm}}}\right)t_{\mathit{\text{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₀+a₁, d_0k=a₀, d_1m=a₀+a₁+γ _m and d_0m=a₀+γ _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 $a{\widehat{p}}_{\mathit{\text{lm}}}(1-{\widehat{p}}_{\mathit{\text{lm}}})$, where ${\widehat{p}}_{\mathit{\text{lm}}}=e^{{\widehat{b}}_{\mathit{\text{lm}}}}/(e^{{\widehat{b}}_{\mathit{\text{lm}}}}+e^{{\widehat{d}}_{\mathit{\text{lm}}}})$, 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 )². By (24) and the discussion that follows, Σ_{A 1}(β) can be estimated by

$${\widehat{\Sigma }}_{A1}\left(\widehat{\beta }\right)=\sum _{a,l,m}\widehat{\rho }(a,l,m)\frac{1-{\widehat{\rho }}^v(a,l,m)}{{\widehat{\rho }}^v(a,l,m)}{G}_{\mathit{\text{lm}}}\left\{a{\widehat{p}}_{\mathit{\text{lm}}}\right(1-{\widehat{p}}_{\mathit{\text{lm}}})+{\nu }_{a,l,m}\},$$

(26)

where $\widehat{\rho }(a,l,m)=n_{\mathit{\text{alm}}}/n$, ${\widehat{\rho }}^v(a,l,m)=n_{\mathit{\text{alm}}}^v/n_{\mathit{\text{alm}}}$ and G _lm be the value of $G_i=z_i^{\otimes 2}$ when subject i belongs to the category (l,m). Hence the covariance matrix of ${\widehat{\beta }}_{I1}$ can be estimated by ${Î}^{-1}\left(\widehat{\beta }\right)+{Î}^{-1}\left(\widehat{\beta }\right){\widehat{\Sigma }}_{A1}\left(\widehat{\beta }\right){Î}^{-1}\left(\widehat{\beta }\right)$ using the automated data.

Estimation with${\widehat{\beta }}_{A2}$using the automated data. The estimating equations (25) are equivalent to the following nonlinear equations for {b _lm , for l=0,1, m=1,⋯,k},

$$\begin{array}{l}\sum _{m=1}^k\left\{\sum _{a\in A}{y}_{\mathit{\text{alm}}}{\lambda }_{\mathit{\text{alm}}}-e^{b_{\mathit{\text{lm}}}}{t}_{\mathit{\text{lm}}}+\frac{e^{b_{\mathit{\text{lm}}}}}{e^{b_{\mathit{\text{lm}}}}+e^{{\widehat{d}}_{1m}}}\left({\alpha }_{\mathit{\text{lm}}}-\sum _{a\in A}a{n}_{\mathit{\text{alm}}}{\lambda }_{\mathit{\text{alm}}}\right)\right\}=0,\\ \sum _{l=0,1}\left\{\sum _{a\in A}{y}_{\mathit{\text{alm}}}{\lambda }_{\mathit{\text{alm}}}-e^{b_{\mathit{\text{lm}}}}{t}_{\mathit{\text{lm}}}+\frac{e^{b_{\mathit{\text{lm}}}}}{e^{b_{\mathit{\text{lm}}}}+e^{{\widehat{d}}_{\mathit{\text{lm}}}}}\left({\alpha }_{\mathit{\text{lm}}}-\sum _{a\in A}a{n}_{\mathit{\text{alm}}}{\lambda }_{\mathit{\text{alm}}}\right)\right\}=0,\end{array}$$

for l=0,1 and m=1,…,k−1.

Since Var{S _β (Y|X)|Z≃(l,m),T,A=a}=a p _lm (1−p _lm )G _lm , Σ_{A 2}(β) can be consistently estimated by

$${\widehat{\Sigma }}_{A2}\left(\widehat{\beta }\right)=\sum _{a,l,m}\widehat{\rho }(a,l,m)\frac{1-{\widehat{\rho }}^v(a,l,m)}{{\widehat{\rho }}^v(a,l,m)}a{\widehat{p}}_{\mathit{\text{lm}}}(1-{\widehat{p}}_{\mathit{\text{lm}}})G_{\mathit{\text{lm}}}.$$

Hence the covariance matrix of ${\widehat{\beta }}_{A2}$ can be estimated by ${Î}^{-1}\left(\widehat{\beta }\right)+{Î}^{-1}\left(\widehat{\beta }\right){\widehat{\Sigma }}_{A2}\left(\widehat{\beta }\right){Î}^{-1}\left(\widehat{\beta }\right)$ 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

$$w_m=\frac{{\alpha }_{0m}{\alpha }_{1m}{y}_{0m}{y}_{1m}}{{\alpha }_{0m}{y}_{0m}{n}_{1m}^v+{\alpha }_{1m}{y}_{1m}{n}_{0m}^v}/\sum _{m=1}^k\frac{{\alpha }_{0m}{\alpha }_{1m}{y}_{0m}{n}_{1m}^v}{{\alpha }_{0m}{y}_{0m}{n}_{1m}^v+{\alpha }_{1m}{y}_{1m}{n}_{0m}^v}.$$

Then an estimate of variance of ${\widehat{b}}_1$ is given by

$$\widehat{\text{Var}\left({\widehat{b}}_1\right)}=\sum _{m=1}^k{w}_m^2\left(\frac{1}{y_{1m}}-\frac{1}{n_{1m}^v}+\frac{1}{{\alpha }_{1m}}+\frac{1}{y_{0m}}-\frac{1}{n_{0m}^v}+\frac{1}{{\alpha }_{0m}}\right),$$

(27)

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.

A simulation study

We conduct a simulation study to examine the finite sample performance of the IPW estimator ${\widehat{\beta }}_{I1}$ and the AIPW estimator ${\widehat{\beta }}_{A2}$. We consider the Poisson regression model (22). The covariates Z₁ and Z₂ 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₁,Z₂) and T, the outcome variable Y follows a Poisson distribution with mean T exp(b₀+b₁Z₁+θ Z₂) where b₀=−0.5, b₁=−0.8 and θ=−0.6, and W follows a Poisson distribution with mean T exp(a₀+a₁Z₁+γ Z₂) where a₀=−1.3, a₁=−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₁ through the logistic regression model logit{π _i (X,A)}=A−Z₁−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 ${\widehat{\beta }}_{I1}$ and ${\widehat{\beta }}_{A2}$ for β=(b₀,b₁,θ). We also compare the performance of the estimators ${\widehat{\beta }}_{I1}$ and ${\widehat{\beta }}_{A2}$ with the complete-case (CC) estimator ${\widehat{\beta }}_C$ 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 Simulation comparison of the IPW estimator ${\widehat{\beta }}_{I1}$ , the AIPW estimator ${\widehat{\beta }}_{A2}$ and the complete-case (CC) estimator ${\widehat{\beta }}_C$ under various sample sizes and selection probabilities

Table 2 Simulation comparison of the IPW estimator ${\widehat{\beta }}_{I1}$ , the AIPW estimator ${\widehat{\beta }}_{A2}$ and the complete-case (CC) estimator ${\widehat{\beta }}_C$ under various sample sizes and selection probabilities

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₀, b₁ 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.

Table 2 gives the results under Model 3 and Model 4. The bias remains small for ${\widehat{\beta }}_{I1}$ and ${\widehat{\beta }}_{A2}$. The empirical standard errors are also close to the corresponding estimated standard errors. The coverage probabilities remain close to the nominal level 95% for all IPW and AIPW estimators. However, the complete-case estimator yields larger bias and incorrect coverage probability because of the association between the validation probability and the auxiliary variable A and/or the covariate Z₁, in which case the missing is not missing completely at random. The AIPW performs better than IPW with smaller standard errors.

An Application

A community-based, nonrandomized, open-label influenza vaccine (CAIV-T) study was conducted in Temple-Belton, Texas during the 2000-2001 influenza season. The total 11,606 healthy children aged 18 months - 18 years were involved in this study and about 20% of them received a single dose of CAIV-T in 2000. The primary clinical outcome was based on an nonspecific case definition called medically attended acute respiratory infection (MAARI), which included all International Classification of Diseases, Ninth Revision, Clinical Modification diagnoses codes (ICD-9 codes 381-383, 460-487) for upper and lower respiratory tract infections, otitis media and sinusitis. MAARI outcomes and demographic data were extracted from the Scott & White Health Plan administrative database. For each visit, one or two International Classification of Diseases, Ninth Revision, Clinical Modification diagnosis codes were listed. Visits for which asthma diagnosis codes alone were noted, without another MAARI code, were excluded. More details about this study can be found in 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.

Table 3 Study data for influenza epidemic season 2000-01, by age and vaccine group (from Halloran et al.[2003])

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₁ be the treatment indicator (1=vaccine and 0=placebo). Let Z₂=(η₁,η₂) be the dummy variables indicating three age groups, where η₁=1 if the age is in the range 1.5–4, η₁=0, otherwise, and η₂=1 if the age is in the range 5–9, η₂=0, otherwise. The reference group is the age 10–18. The exposure time for all children is taken as T=1 year.

Consider a Poisson regression model with mean T exp(b₀+b₁Z₁+θ₁η₁+θ₂η₂). Using the IPW estimator ${\widehat{\beta }}_{I1}$, the estimates (standard errors) are ${\widehat{b}}_0=-0.7659$ (${\widehat{\sigma }}_{b_0}=0.1046$), ${\widehat{b}}_1=-1.5830$ (${\widehat{\sigma }}_{b_1}=0.5017$), ${\widehat{\theta }}_1=-0.5572$ (${\widehat{\sigma }}_{{\theta }_1}=0.2111$) and ${\widehat{\theta }}_2=-0.0199$ (${\widehat{\sigma }}_{{\theta }_2}=0.1472$). The age-adjusted relative rate (RR) in the vaccinated group compared with the unvaccinated group equals $exp\left({\widehat{b}}_1\right)=exp(-1.5830)=0.2054$, which means that the rate of developing MAARI for the vaccinated group is 20% of that for the unvaccinated group. In terms of the vaccine efficacy VE=1−RR=0.7946, this represents about 80% reduction in the risk of developing MAARI for the vaccinated group compared to the unvaccinated group. The 95% confidence interval of RR obtained by using the delta method is (0.0768,0.5490), showing clear evidence that the vaccinated children have less risk of influenza than the unvaccinated children. The 95% confidence interval for VE is (0.4510,0.9232).

Using the AIPW estimator ${\widehat{\beta }}_{A2}$, the estimates (standard errors) are ${\widehat{b}}_0=-2.0703$ (${\widehat{\sigma }}_{b_0}=0.0851$), ${\widehat{b}}_1=-1.8072$ (${\widehat{\sigma }}_{b_1}=0.3786$), ${\widehat{\theta }}_1=0.6452$ (${\widehat{\sigma }}_{{\theta }_1}=0.1966$) and ${\widehat{\theta }}_2=0.6235$ (${\widehat{\sigma }}_{{\theta }_2}=0.1265$). The age-adjusted relative rate (RR) is $exp\left({\widehat{b}}_1\right)=exp(-1.8072)=0.1641$. The estimated VE is 0.8359 and the 95% confidence interval is (0.6553,0.9219). The estimator ${\widehat{\beta }}_{A2}$ yields smaller standard errors and confidence intervals with more precision than using ${\widehat{\beta }}_{I1}$.

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-size-weighted 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 ${\widehat{\beta }}_{A2}$ 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.

Appendix A

Proof of Proposition 1.

Since

$$\begin{array}{l}\sum _{i=1}^n(1-{\xi }_i)Ê\left\{S_{\beta }\left(Y\right|X_i\left)\right|Z_i,A_i\right\}=\sum _{i\in \bar{V}}\sum _{j\in V(Z_i,A_i)}{S}_{\beta }\left(Y_j\right|X_j)/n^V(Z_i,A_i)\\ =\sum _{i\in V}\left\{n^{\bar{V}}(Z_i,A_i)/n^V(Z_i,A_i)\right\}{S}_{\beta }\left(Y_i\right|X_i),\end{array}$$

we have

$$\begin{array}{l}\sum _{i=1}^n{W}_i^{E1}=\sum _{i\in V}\left(1+\frac{n^{\bar{V}}(Z_i,A_i)}{n^V(Z_i,A_i)}\right)S_{\beta }\left(Y_i\right|X_i)=\sum _{i=1}^n\frac{{\xi }_i}{{\widehat{\pi }}_i^z}{S}_{\beta }(Y_i\left|X_i\right)=\sum _{i=1}^n{W}_i^{I1}.\end{array}$$

(A.1)

This shows that the mean score estimator ${\widehat{\beta }}_{E1}$ is the same as the IPW estimator ${\widehat{\beta }}_{I1}$. Further, since

$$\begin{array}{l}\sum _{i=1}^n\left(1-\frac{{\xi }_i}{{\widehat{\pi }}_i^z}\right)Ê\left\{S_{\beta }\right(Y\left|X_i\right)|Z_i,A_i\}\\ =\sum _{i\in \bar{V}}\sum _{j\in V(Z_i,A_i)}\frac{S_{\beta }\left(Y_j\right|X_j)}{n^V(Z_i,A_i)}-\sum _{i\in V}\frac{n^{\bar{V}}(Z_i,A_i)}{n^V(Z_i,A_i)}\sum _{j\in V(Z_i,A_i)}\frac{S_{\beta }\left(Y_j\right|X_j)}{n^V(Z_i,A_i)}\\ =\sum _{i\in V}\frac{n^{\bar{V}}(Z_i,A_i)}{n^V(Z_i,A_i)}{S}_{\beta }\left(Y_i\right|X_i)-\sum _{i\in V}\frac{n^{\bar{V}}(Z_i,A_i)}{n^V(Z_i,A_i)}{S}_{\beta }(Y_i\left|X_i\right)=0,\end{array}$$

we have $\sum _{i=1}^n{W}_i^{A1}=\sum _{i=1}^n{W}_i^{I1}$. Thus the AIPW estimator ${\widehat{\beta }}_{A1}$, the IPW estimator ${\widehat{\beta }}_{I1}$ and the mean score estimator ${\widehat{\beta }}_{E1}$ are equivalent to each other.

Note that

$$\begin{array}{l}\sum _{i=1}^n\left(1-\frac{{\xi }_i}{{\widehat{\pi }}_i^z}\right)Ê\left\{S_{\beta }\right(Y\left|X_i\right)|X_i,A_i\}\\ =\sum _{i\in \bar{V}}\sum _{j\in V(X_i,A_i)}\frac{S_{\beta }\left(Y_j\right|X_j)}{n^V(X_i,A_i)}-\sum _{i\in V}\frac{n^{\bar{V}}(Z_i,A_i)}{n^V(Z_i,A_i)}\sum _{j\in V(X_i,A_i)}\frac{S_{\beta }\left(Y_j\right|X_j)}{n^V(X_i,A_i)}\\ =\sum _{i\in V}\frac{n^{\bar{V}}(X_i,A_i)}{n^V(X_i,A_i)}{S}_{\beta }\left(Y_i\right|X_i)-\sum _{i\in V}\frac{n^V(X_i,A_i)}{n^V(X_i,A_i)}\frac{n^{\bar{V}}(Z_i,A_i)}{n^V(Z_i,A_i)}{S}_{\beta }(Y_i\left|X_i\right)\\ =\sum _{i\in V}\left\{\frac{n^{\bar{V}}(X_i,A_i)}{n^V(X_i,A_i)}-\frac{n^{\bar{V}}(Z_i,A_i)}{n^V(Z_i,A_i)}\right\}{S}_{\beta }\left(Y_i\right|X_i),\end{array}$$

(A.2)

which is not zero unless $Z_i^c$ is linearly related to Z _i and in this case β is not identifiable. Hence the AIPW estimator ${\widehat{\beta }}_{A2}$ is different from the AIPW estimator ${\widehat{\beta }}_{A1}$.

By (A.1) and (A.2), we have

$$\begin{array}{lcr}\sum _{i=1}^n{W}_i^{A2}& =& \sum _{i\in V}\left\{\frac{n(Z_i,A_i)}{n^V(Z_i,A_i)}+\frac{n^{\bar{V}}(X_i,A_i)}{n^V(X_i,A_i)}-\frac{n^{\bar{V}}(Z_i,A_i)}{n^V(Z_i,A_i)}\right\}{S}_{\beta }\left(Y_i\right|X_i)\\ =& \sum _{i\in V}\left\{1+\frac{n^{\bar{V}}(X_i,A_i)}{n^V(X_i,A_i)}\right\}{S}_{\beta }\left(Y_i\right|X_i)=\sum _{i=1}^n{W}_i^{I2}.\end{array}$$

Following the same arguments leading to (A.1), we also have $\sum _{i=1}^n{W}_i^{E2}=\sum _{i=1}^n{W}_i^{I2}$. Hence, the estimators ${\widehat{\beta }}_{I2}$, ${\widehat{\beta }}_{E2}$ and ${\widehat{\beta }}_{A2}$ are equivalent. By following the steps in (A.2), we also have $\sum _{i=1}^n\left(1-\frac{{\xi }_i}{{\widehat{\pi }}_i}\right)Ê\left\{S_{\beta }\left(Y\right|X_i\left)\right|X_i,A_i\right\}=0$. Hence, ${\widehat{\beta }}_{A3}$ is the same as ${\widehat{\beta }}_{I2}$. Therefore, these are essentially two different estimators. □

Proof of Theorem 1.

Applying the first order Taylor expansion, ${\tilde{\pi }}_i-{\pi }_i={\left(\mathrm{\partial \pi }\right(X_i,A_i,{\psi }_0)/\mathrm{\partial \psi })}^{\prime }\left(\widehat{\psi }-{\psi }_0\right)+o_p\left(n^{-1/2}\right)$. From (13), we have

$$n^{-1/2}{U}_I=n^{-1/2}\sum _{i=1}^n\frac{{\xi }_i}{{\pi }_i}{S}_{\beta }\left(Y_i\right|X_i)-n^{-1/2}\sum _{i=1}^n\frac{{\tilde{\pi }}_i-{\pi }_i}{{\tilde{\pi }}_i{\pi }_i}{\xi }_i{S}_{\beta }(Y_i\left|X_i\right)$$

(A.3)

The second term of (A.3) is

$$\begin{array}{l}{n}^{-1/2}\sum _{i=1}^n\frac{{\tilde{\pi }}_i-{\pi }_i}{{\tilde{\pi }}_i{\pi }_i}{\xi }_i{S}_{\beta }\left(Y_i\right|X_i)\\ =n^{-1/2}\sum _{i=1}^n{{\pi }_i}^{-2}{\xi }_i{S}_{\beta }\left(Y_i\right|X_i){\left(\frac{\mathrm{\partial \pi }(X_i,A_i,{\psi }_0)}{\mathrm{\partial \psi }}\right)}^{\prime }\left(\widehat{\psi }-{\psi }_0\right)\\ =E\left\{{{\pi }_i}^{-2}{\xi }_i{S}_{\beta }\left(Y_i\right|X_i){\left(\frac{\mathrm{\partial \pi }(X_i,A_i,{\psi }_0)}{\mathrm{\partial \psi }}\right)}^{\prime }\right\}{n}^{1/2}\left(\widehat{\psi }-{\psi }_0\right)+o_p\left(1\right)\end{array}$$

(A.4)

By (11), (A.3) and (A.4), we have

$$\begin{array}{l}{n}^{-1/2}{U}_I=n^{-1/2}\sum _{i=1}^n\left(\frac{{\xi }_i}{{\pi }_i}{S}_{\beta }\left(Y_i\right|X_i)-O_i\right)+o_p\left(1\right)=n^{-1/2}\sum _{i=1}^n{Q}_i^I+o_p\left(1\right).\end{array}$$

Now consider the AIPW estimator ${\widehat{\beta }}_A$ based on solving the estimating equation (14). For simplicity, we denote E _a {S _β (Y|X _i )|X _i ,A _i } by E _i and $\tilde{E}\left\{S_{\beta }\left(Y\right|X_i\left)\right|X_i,A_i\right\}$ by ${\tilde{E}}_i$. We note that

$$\begin{array}{lcr}{n}^{-1/2}{U}_A& =& n^{-1/2}\sum _{i=1}^n\left[\frac{{\xi }_i}{{\pi }_i}{S}_{\beta }\left(Y_i\right|X_i)+\left(1-\frac{{\xi }_i}{{\pi }_i}\right)E_i\right]\\ +n^{-1/2}\sum _{i=1}^n\left[\left(\frac{{\xi }_i}{{\tilde{\pi }}_i}-\frac{{\xi }_i}{{\pi }_i}\right)\left\{S_{\beta }\left(Y_i\right|X_i)-E_i\right\}+\left(1-\frac{{\xi }_i}{{\tilde{\pi }}_i}\right)\left({\tilde{E}}_i-E_i\right)\right].\end{array}$$

Suppose that ${\tilde{\pi }}_i$ and ${\tilde{E}}_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

$$\begin{array}{l}{n}^{-1/2}{U}_A=n^{-1/2}\sum _{i=1}^n\left[\frac{{\xi }_i}{{\pi }_i}{S}_{\beta }\left(Y_i\right|X_i)+\left(1-\frac{{\xi }_i}{{\pi }_i}\right)E_i\right]+o_p\left(1\right).\end{array}$$

It can be shown that under MAR I, $n^{-1}\partial U_I/\mathrm{\partial \beta }\stackrel{P}{\to }I\left(\beta \right)$ and $n^{-1}\partial U_A/\mathrm{\partial \beta }\stackrel{P}{\to }I\left(\beta \right)$. By routine derivations, we have

$$\begin{array}{lcr}{n}^{1/2}\left({\widehat{\beta }}_I-\beta \right)& =& I^{-1}\left(\beta \right)n^{-1/2}{U}_I+o_p\left(1\right)=I^{-1}\left(\beta \right)n^{-1/2}\sum _{i=1}^n{Q}_i^I+o_p\left(1\right),\\ n^{1/2}\left({\widehat{\beta }}_A-\beta \right)& =& I^{-1}\left(\beta \right)n^{-1/2}{U}_A+o_p\left(1\right)=I^{-1}\left(\beta \right)n^{-1/2}\sum _{i=1}^n{Q}_i^A+o_p\left(1\right).\end{array}$$

By the central limit theorem, both $n^{1/2}\left({\widehat{\beta }}_I-\beta \right)$ and $n^{1/2}\left({\widehat{\beta }}_A-\beta \right)$ have asymptotically normal distributions with mean zero and covariances equal to $I^{-1}\left(\beta \right)\text{Var}\left(Q_i^I\right)I^{-1}\left(\beta \right)$ and $I^{-1}\left(\beta \right)\text{Var}\left(Q_i^A\right)I^{-1}\left(\beta \right)$, respectively.

Next, we examine the covariance matrices $\text{Var}\left(Q_i^I\right)$ and $\text{Var}\left(Q_i^A\right)$ to understand the efficiency gain of ${\widehat{\beta }}_A$ over ${\widehat{\beta }}_I$. Note that $Q_i^I={\xi }_i/{\pi }_i{S}_{\beta }\left(Y_i\right|X_i)-O_i$ and $Q_i^A={\xi }_i/{\pi }_i{S}_{\beta }\left(Y_i\right|X_i)+(1-{\xi }_i/{\pi }_i)E_i$. Denote A _i =ξ _i /π _i S _β (Y _i |X _i ) and B _i =(1−ξ _i /π _i )E _i . Then $Q_i^I=Q_i^A-B_i-O_i$. Under MAR I, $\text{Cov}\left(Q_i^A,O_i\right)=\mathrm{E}\left(Q_i^A{O}_i\right)=E\left\{E\left(Q_i^A\right|{\xi }_i,X_i,A_i)O_i\right\}=E\left\{E_i{O}_i\right\}=E\left\{E_{i}E\right(O_i|X_i,A_i)\}=0,$ and

$$\begin{array}{lcr}\text{Cov}(Q_i^A,B_i)& =& E\left[\left\{\frac{{\xi }_i}{{\pi }_i}{S}_{\beta }\left(Y_i\right|X_i)+\left(1-\frac{{\xi }_i}{{\pi }_i}\right)E_i\right\}\left(1-\frac{{\xi }_i}{{\pi }_i}\right)E_i\right]\\ =& E\left[{\left(1-\frac{{\xi }_i}{{\pi }_i}\right)}^2{E}_i^2\right]-E\left[{\left(1-\frac{{\xi }_i}{{\pi }_i}\right)}^2{S}_{\beta }\left(Y_i\right|X_i)E_i\right]\\ +E\left[\left(1-\frac{{\xi }_i}{{\pi }_i}\right)S_{\beta }\left(Y_i\right|X_i)E_i\right]=0.\end{array}$$

Hence, $\text{Cov}\left(Q_i^A,B_i+O_i\right)=0$. It follows that $\text{Var}\left(Q_i^I\right)=\text{Var}\left(Q_i^A\right)+\text{Var}(B_i+O_i)$. Since $Q_i^A=S_{\beta }\left(Y_i\right|X_i)-\left(1-\frac{{\xi }_i}{{\pi }_i}\right)\left\{S_{\beta }\left(Y_i\right|X_i)-E_i\right\}$ and the two terms are uncorrelated under MAR I, we have $\text{Var}\left(Q_i^A\right)=\text{Var}\left(S_{\beta }\left(Y_i\right|X_i)\right)+\text{Var}\left((1-{\xi }_i/{\pi }_i)\left\{S_{\beta }\left(Y_i\right|X_i)-E_i\right\}\right)$, where the first term equals I(β). This completes the proof of Theorem 1. □

Proof of Corollary 1.

Let $Q_i^{A1}={\xi }_i/{\pi }_i^z{S}_{\beta }\left(Y_i\right|X_i)+(1-{\xi }_i/{\pi }_i^z)E\left\{S_{\beta }\left(Y\right|X_i\left)\right|Z_i,A_i\right\}$, and $Q_i^{A2}={\xi }_i/{\pi }_i^z{S}_{\beta }\left(Y_i\right|X_i)+\left(1-{\xi }_i/{\pi }_i^z\right)E\left\{S_{\beta }\left(Y\right|X_i\left)\right|X_i,A_i\right\}$. By (16),

$$\begin{array}{lcr}\text{Var}\left(Q_i^{A1}\right)& =& \text{Var}\left(S_{\beta }\left(Y_i\right|X_i)\right)+\text{Var}\left(\left(1-\frac{{\xi }_i}{{\pi }_i^z}\right)\left[S_{\beta }\left(Y_i\right|X_i)-E\left\{S_{\beta }\left(Y\right|X_i\left)\right|Z_i,A_i\right\}\right]\right),\\ \text{Var}\left(Q_i^{A2}\right)& =& \text{Var}\left(S_{\beta }\left(Y_i\right|X_i)\right)+\text{Var}\left(\left(1-\frac{{\xi }_i}{{\pi }_i^z}\right)\left[S_{\beta }\left(Y_i\right|X_i)-E\left\{S_{\beta }\left(Y\right|X_i\left)\right|X_i,A_i\right\}\right]\right).\end{array}$$

The second term of $\text{Var}\left(Q_i^{A1}\right)$ equals Σ_{A 1}(β) and the second term of $\text{Var}\left(Q_i^{A2}\right)$ equals Σ_{A 2}(β). 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 $Q_i^{A1}$ and $Q_i^{A2}$ contributes to the difference in the asymptotic variances of ${\widehat{\beta }}_{A1}$ and ${\widehat{\beta }}_{A2}$. Since $E\left\{{\left(1-{\xi }_i/{\pi }_i^z\right)}^2|Y_i,X_i,A_i\right\}=E\left\{{\left(1-{\xi }_i/{\pi }_i^z\right)}^2|X_i,A_i\right\}=(1-{\pi }_i^z)/{\pi }_i^z$ under MAR I,

$$\begin{array}{lcr}{\Sigma }_{A2}\left(\beta \right)& =& E\left(\frac{1-{\pi }_i^z}{{\pi }_i^z}\left\{E\left\{S_{\beta }^2\left(Y_i\right|X_i\left)\right|X_i,A_i\right\}-{\left[E\left\{S_{\beta }\left(Y\right|X_i\left)\right|X_i,A_i\right\}\right]}^2\right\}\right)\\ =& E\left(\frac{1-{\pi }_i^z}{{\pi }_i^z}\left\{E\left\{S_{\beta }^2\left(Y_i\right|X_i\left)\right|Z_i,A_i\right\}-{\left[E\left\{S_{\beta }\left(Y\right|X_i\left)\right|Z_i,A_i\right\}\right]}^2\right\}\right)\\ -E\left(\frac{1-{\pi }_i^z}{{\pi }_i^z}\left\{{\left[E\left\{S_{\beta }\left(Y\right|X_i\left)\right|X_i,A_i\right\}\right]}^2-{\left[E\left\{S_{\beta }\left(Y\right|X_i\left)\right|Z_i,A_i\right\}\right]}^2\right\}\right)\\ =& {\Sigma }_{A1}\left(\beta \right)-E\left(\frac{1-{\pi }_i^z}{{\pi }_i^z}\text{Var}\left(\left[E\left\{S_{\beta }\left(Y\right|X_i\left)\right|X_i,A_i\right\}\right]|Z_i,A_i\right)\right),\end{array}$$

which is less than Σ_{A 1}(β) 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). We note that

$$\begin{array}{l}E\left\{{{\pi }_i}^{-1}{S}_{\beta }\left(Y_i\right|X_i){\left(\frac{\mathrm{\partial \pi }(Z_i,A_i,{\psi }_0)}{\mathrm{\partial \psi }}\right)}^{\prime }\right\}\\ =E\left[{{\pi }_i}^{-1}E\left\{S_{\beta }\left(Y_i\right|X_i\left)\right|Z_i,A_i\right\}{\left(\frac{\mathrm{\partial \pi }(Z_i,A_i,{\psi }_0)}{\mathrm{\partial \psi }}\right)}^{\prime }\right]\\ =\sum _{z,a}\rho (z,a){\pi }^{-1}(z,a,{\psi }_0)E\left\{S_{\beta }\left(Y\right|X\left)\right|Z=z,A=a\right\}{\left(\frac{\mathrm{\partial \pi }(z,a,{\psi }_0)}{\mathrm{\partial \psi }}\right)}^{\prime },\end{array}$$

and

$$\begin{array}{lcr}{n}^{-1/2}\sum _{i=1}^n{S}_i^{\psi }& =& n^{-1/2}\sum _{i=1}^n\frac{({\xi }_i-\pi (Z_i,A_i,{\psi }_0\left)\right)}{\pi (Z_i,A_i,{\psi }_0)(1-\pi (Z_i,A_i,{\psi }_0\left)\right)}\frac{\mathrm{\partial \pi }(Z_i,A_i,{\psi }_0)}{\mathrm{\partial \psi }}\\ =& \sum _{z,a}{n}^{-1/2}\sum _{j=1}^n\frac{\left({\xi }_j-\pi (z,a,{\psi }_0)\right)I\left(Z_j=z,A_j=a\right)}{\pi (z,a,{\psi }_0)(1-\pi (z,a,{\psi }_0\left)\right)}\frac{\mathrm{\partial \pi }(z,a,{\psi }_0)}{\mathrm{\partial \psi }}.\end{array}$$

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, $\frac{\mathrm{\partial \pi }(z,a,{\psi }_0)}{\mathrm{\partial \psi }}$ 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.

We have

$$\begin{array}{lcr}{n}^{-1/2}\sum _{i=1}^n{O}_i& =& E\left\{{{\pi }_i}^{-1}{S}_{\beta }\left(Y_i\right|X_i){\left(\frac{\mathrm{\partial \pi }(Z_i,A_i,{\psi }_0)}{\mathrm{\partial \psi }}\right)}^{\prime }\right\}{\left(I^{\psi }\right)}^{-1}{n}^{-1/2}\sum _{i=1}^n{S}_i^{\psi }\\ =& \sum _{z,a}E\left\{S_{\beta }\left(Y\right|X\left)\right|Z=z,A=a\right\}{n}^{-1/2}\sum _{j=1}^n\frac{{\xi }_j-\pi (z,a,{\psi }_0)}{\pi (z,a,{\psi }_0)}I\left(Z_j=z,A_j=a\right),\end{array}$$

and

$$\begin{array}{l}{n}^{-1/2}\sum _{i=1}^n{B}_i=-n^{-1/2}\sum _{i=1}^n\frac{({\xi }_i-\pi (Z_i,A_i\left)\right)}{\pi (Z_i,A_i)}{E}_i\\ =-n^{-1/2}\sum _{i=1}^n\frac{({\xi }_i-\pi (Z_i,A_i\left)\right)}{\pi (Z_i,A_i)}E\left\{S_{\beta }\left(Y_i\right|X_i\left)\right|X_i,A_i\right\}\\ =-\sum _{z,a}{n}^{-1/2}\sum _{j=1}^n\frac{{\xi }_j-\pi (z,a,{\psi }_0)}{\pi (z,a,{\psi }_0)}I\left(Z_j=z,A_j=a\right)E\left\{S_{\beta }\left(Y_j\right|X_j\left)\right|Z_j=z,Z_j^c,A_j=a\right\}.\end{array}$$

Then (21) holds. It follows that ${\widehat{\beta }}_A$ is more efficient than ${\widehat{\beta }}_I$ unless V a r{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. □

Appendix B

Proof of the simplified variance formula ( 27 )

The information matrix I(β) is a (k+1)×(k+1) symmetric matrix given by

$$\begin{array}{l}I(b_0,b_1,\theta )=\left[\begin{array}{lllll}\sum _{l,m}{q}_{\mathit{\text{lm}}}& \sum _{m=1}^k{q}_{1m}& q_{11}+q_{01}& \cdots & q_{1r}+q_{0r}\\ \sum _{m=1}^k{q}_{1m}& \sum _{m=1}^k{q}_{1m}& q_{11}& \cdots & q_{1r}\\ q_{11}+q_{01}& q_{11}& q_{11}+q_{01}& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ q_{1r}+q_{0r}& q_{1r}& 0& \cdots & q_{1r}+q_{0r}\end{array}\right],\end{array}$$

where r=k−1 and $q_{\mathit{\text{lm}}}=E\left(T_i{e}^{b_{\mathit{\text{lm}}}}I\right\{\text{individual}i\text{in category}(l,m)\left\}\right)$. For ease of presentation in the following, we drop the augments (b₀,b₁,θ) and use I for I(b₀,b₁,θ). Let I _ij be the cofactor of the (i,j)th element of I.

First we need to find the elements on the second row of the information matrix I(b₀,b₁,θ). Note that for a matrix A, ((a _ij )_n×n)⁻¹=(A _ji )_n×n/|A|, where A _ij is the cofactor of a _ij in the matrix A and |A| is the determinant of A. Also note that for a block matrix

$$B=\left(\begin{array}{ll}{B}_{11}& B_{12}\\ B_{21}& B_{22}\end{array}\right),$$

the determinant $\left|B\right|=\left|B_{22}\right||B_{11}-B_{12}{B}_{22}^{-1}{B}_{21}|$. We have

$$\begin{array}{lcr}\left|I\right|& =& \prod _{m=1}^{k-1}(q_{1m}+q_{0m})\\ \left|\left(\begin{array}{ll}\sum _{\mathit{\text{lm}}}{q}_{\mathit{\text{lm}}}& \sum _{m=1}^k{q}_{1m}\\ \sum _{m=1}^k{q}_{1m}& \sum _{m=1}^k{q}_{1m}\end{array}\right)-\left(\begin{array}{ll}\sum _{l=0}^1\sum _{m=1}^{k-1}{q}_{\mathit{\text{lm}}}& \sum _{m=1}^{k-1}{q}_{1m}\\ \sum _{m=1}^{k-1}{q}_{1m}& \sum _{m=1}^{k-1}\frac{{\left(q_{1m}\right)}^2}{q_{1m}+q_{0m}}\end{array}\right)\right|\\ =& \prod _{m=1}^{k-1}(q_{1m}+q_{0m})\left|\begin{array}{ll}{q}_{1k}+q_{0k}& q_{1k}\\ q_{1k}& q_{1k}+\sum _{m=1}^{k-1}\frac{q_{0m}{q}_{1m}}{q_{1m}+q_{0m}}\end{array}\right|\\ =& \prod _{m=1}^k(q_{1m}+q_{0m})\sum _{m=1}^k\frac{q_{0m}{q}_{1m}}{q_{1m}+q_{0m}},\end{array}$$

and

$$\begin{array}{lcr}{I}_{21}& =& -\left|\begin{array}{llll}\sum _{m=1}^k{q}_{1m}& q_{11}+q_{01}& \cdots & q_{1r}+q_{0r}\\ q_{11}& q_{11}+q_{01}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ q_{1r}& 0& \cdots & q_{1r}+q_{0r}\end{array}\right|\\ =& -\prod _{m=1}^{k-1}(q_{1m}+q_{0m})\left[\sum _{m=1}^k{q}_{1m}-\sum _{m=1}^{k-1}{q}_{1m}\right]\\ =& -q_{1k}\prod _{m=1}^{k-1}(q_{1m}+q_{0m}).\end{array}$$

Hence, the (2,1)th element of I⁻¹ is

$${\left(I^{-1}\right)}_{21}=\frac{I_{21}}{\left|I\right|}=-\frac{q_{1k}}{W(q_{1k}+q_{0k})},$$

where $W=\sum _{m=1}^k{q}_{0m}{q}_{1m}/(q_{1m}+q_{0m})$.

To calculate (2,2)th element of I⁻¹, we have

$$\begin{array}{lcr}{I}_{22}& =& \left|\begin{array}{llll}\sum _{l,m}^k{q}_{\mathit{\text{lm}}}& q_{11}+q_{01}& \cdots & q_{1r}+q_{0r}\\ q_{11}+q_{01}& q_{11}+q_{01}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ q_{1r}+q_{0r}& 0& \cdots & q_{1r}+q_{0r}\end{array}\right|\\ =& \prod _{m=1}^{k-1}(q_{1m}+q_{0m})\left[\sum _{l,m}{q}_{\mathit{\text{lm}}}-\sum _{l=0}^1\sum _{m=1}^{k-1}{q}_{\mathit{\text{lm}}}\right]=\prod _{m=1}^k(q_{1m}+q_{0m}).\end{array}$$

Hence, the (2,2)th element of I⁻¹ is (I⁻¹)₂₂=I₂₂/|I|=1/W.

To calculate (2,3)th element of I⁻¹, we have

$$\begin{array}{lc}{I}_{23}& =-\left|\begin{array}{lllll}\sum _{l,m}{q}_{\mathit{\text{lm}}}& \sum _{m=1}^k{q}_{1m}& q_{12}+q02& \cdots & q_{1r}+q_{0r}\\ q_{11}+q_{01}& q_{11}& 0& \cdots & 0\\ q_{12}+q_{02}& q_{12}& q_{12}+q_{02}& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ q_{1r}+q_{0r}& q_{1r}& 0& \cdots & q_{1r}+q_{0r}\end{array}\right|\\ =-\prod _{m=2}^{k-1}(q_{1m}+q_{0m})\left|\left(\begin{array}{ll}\sum _{\mathit{\text{lm}}}{q}_{\mathit{\text{lm}}}& \sum _{m=1}^k{q}_{1m}\\ q_{11}+q_{01}& q_{11}\end{array}\right)\right.\\ \left.-\left(\begin{array}{ll}\sum _{m=2}^{k-1}(q_{1m}+q_{0m})& \sum _{m=2}^{k-1}{q}_{1m}\\ 0& 0\end{array}\right)\right|\\ =-\prod _{m=2}^{k-1}(q_{1m}+q_{0m})(q_{0k}{q}_{11}-q_{1k}{q}_{01}).\end{array}$$

Hence

$$\begin{array}{lcr}{\left(I^{-1}\right)}_{23}& =& \frac{I_{23}}{\left|I\right|}=-\frac{q_{11}}{W(q_{11}+q_{01})}+\frac{q_{1k}}{W(q_{1k}+q_{0k})}.\end{array}$$

To calculate (2,4)th element of I⁻¹, we have

$$\begin{array}{l}{I}_{24}=\left|\begin{array}{llllll}\sum _{l,m}{q}_{\mathit{\text{lm}}}& \sum _{m=1}^k{q}_{1m}& q_{11}+q_{01}& q_{13}+q_{03}& \cdots & q_{1r}+q_{0r}\\ q_{11}+q_{01}& q_{11}& q_{11}+q_{01}& 0& \cdots & 0\\ q_{12}+q_{02}& q_{12}& 0& 0& \cdots & 0\\ q_{13}+q_{03}& q_{13}& 0& q_{13}+q_{03}& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ q_{1r}+q_{0r}& q_{1r}& 0& 0& \cdots & q_{1r}+q_{0r}\end{array}\right|.\end{array}$$

By switching the 2nd and the 3rd row, we have

$$\begin{array}{l}{I}_{24}=-n^k\left|\begin{array}{llllll}\sum _{l,m}{q}_{\mathit{\text{lm}}}& \sum _{m=1}^k{q}_{1m}& q_{11}+q_{01}& q_{13}+q_{03}& \cdots & q_{1r}+q_{0r}\\ q_{12}+q_{02}& q_{12}& 0& 0& \cdots & 0\\ q_{11}+q_{01}& q_{11}& q_{11}+q_{01}& 0& \cdots & 0\\ q_{13}+q_{03}& q_{13}& 0& q_{13}+q_{03}& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ q_{1r}+q_{0r}& q_{1r}& 0& 0& \cdots & q_{1r}+q_{0r}\end{array}\right|.\end{array}$$

Hence

$$\begin{array}{l}{\left(I^{-1}\right)}_{24}=\frac{I_{24}}{\left|I\right|}=-\frac{q_{12}}{W(q_{12}+q_{02})}+\frac{q_{1k}}{W(q_{1k}+q_{0k})}.\end{array}$$

In general, to calculate I_2(m+2) for m=1,⋯,k−1, we can obtain a matrix with a (k−2)×(k−2) diagonal right lower block by switching rows of I_2(m+2) even number of times when m is odd and by switching rows odd number of times when m is even. Then similar to calculating (I⁻¹)₂₄, we have

$${\left(I^{-1}\right)}_{2(m+2)}=-\frac{q_{1m}}{W\left(q_{1m}+q_{0m}\right)}+\frac{q_{1k}}{W(q_{1k}+q_{0k})}.$$

For l=1 and 1≤m≤k, the (i,j)th element of G _lm is g _ij =1 for (i,j)=(1,1),(1,2),(1,m+2),(2,1),(2,2),(2,m+2),(m+2,1),(m+2,2),(m+2,m+2), and g _ij =0 elsewhere. We have the (2,2)th element of I⁻¹G _lm I⁻¹ as

$$\begin{array}{l}{\left(I^{-1}\right)}_{21}^2+{\left(I^{-1}\right)}_{22}^2+{\left(I^{-1}\right)}_{2(m+2)}^2+2{\left(I^{-1}\right)}_{21}{\left(I^{-1}\right)}_{2(m+2)}\\ +2{\left(I^{-1}\right)}_{22}{\left(I^{-1}\right)}_{2(m+2)}+2{\left(I^{-1}\right)}_{21}{\left(I^{-1}\right)}_{22}\\ =\frac{1}{W^2}{\left(\frac{q_{0m}}{q_{1m}+q_{0m}}\right)}^2.\end{array}$$

Since for l=0 and 1≤m≤k, g _ij =1 for (i,j)=(1,1),(1,m+2),(m+2,1),(m+2,m+2), and g _ij =0 elsewhere, in this case the (2,2)th element of I⁻¹G _lm I⁻¹ is

$$\begin{array}{l}{\left(I^{-1}\right)}_{21}^2+{\left(I^{-1}\right)}_{2(m+2)}^2+2{\left(I^{-1}\right)}_{21}{\left(I^{-1}\right)}_{2(m+2)}\\ ={({\left(I^{-1}\right)}_{21}+{\left(I^{-1}\right)}_{2(m+2)})}^2\\ =\frac{1}{W^2}{\left(\frac{q_{1m}}{q_{1m}+q_{0m}}\right)}^2.\end{array}$$

Hence the (2,2)th element of the asymptotic covariance matrix of ${\widehat{\beta }}_{A2}$ is given by ${\sigma }_{b_1}^2=1/W+U/W^2$, where

$$U=\sum _{\alpha ,l,m}\rho (\alpha ,l,m)\frac{1-{\rho }^V(\alpha ,l,m)}{{\rho }^V(\alpha ,l,m)}\alpha p_{\mathit{\text{lm}}}(1-p_{\mathit{\text{lm}}}){\left(\frac{q_{(1-l)m}}{q_{0m}+q_{1m}}\right)}^2.$$

Note that P(Y _i =1|A _i =1,i∈ category (l,m)) can be estimated by $y_{\mathit{\text{lm}}}/n_{\mathit{\text{lm}}}^v$and ρ(l,m) by α _lm /n. Thus q _lm can be estimated by $({\alpha }_{\mathit{\text{lm}}}/n)(y_{\mathit{\text{lm}}}/n_{\mathit{\text{lm}}}^v)$. Then we can estimate W by$Ŵ=n^{-1}\sum _{m=1}^k{\alpha }_{0m}{\alpha }_{1m}{y}_{0m}{y}_{1m}/({\alpha }_{0m}{y}_{0m}{n}_{1m}^v+{\alpha }_{1m}{y}_{1m}{n}_{0m}^v)$. By replacing ${\rho }_{l,m}^v$ with $n_{\mathit{\text{lm}}}^v/{\alpha }_{\mathit{\text{lm}}}$ and p _lm with $y_{\mathit{\text{lm}}}/n_{\mathit{\text{lm}}}^v$, we can estimate U by

$$\begin{array}{lc}Û& =n^{-1}\sum _{l,m}{\alpha }_{\mathit{\text{lm}}}\frac{{\alpha }_{\mathit{\text{lm}}}-n_{\mathit{\text{lm}}}^v}{n_{\mathit{\text{lm}}}^v}\frac{y_{\mathit{\text{lm}}}}{n_{\mathit{\text{lm}}}^v}\frac{n_{\mathit{\text{lm}}}^v-y_{\mathit{\text{lm}}}}{n_{\mathit{\text{lm}}}^v}\left({\alpha }_{(1-l)m}\frac{y_{(1-l)m}}{n_{(1-l)m}^v}\right){\left({\alpha }_{1m}\frac{y_{1m}}{n_{1m}^v}+{\alpha }_{0m}\frac{y_{0m}}{n_{0m}^v}\right)}^{-2}\\ =n^{-1}\sum _{m=1}^k\left[\frac{{\alpha }_{1m}-n_{1m}^v}{n_{1m}^v}\frac{n_{1m}^v-y_{1m}}{n_{1m}^v}{\alpha }_{0m}\frac{y_{0m}}{n_{0m}^v}+{\alpha }_{0m}-n_{0m}^v{n}_{0m}^v\frac{n_{0m}^v-y_{0m}}{n_{0m}^v}{\alpha }_{1m}\frac{y_{1m}}{n_{1m}^v}\right]\\ \times \left({\alpha }_{1m}\frac{y_{1m}}{n_{1m}^v}{\alpha }_{0m}\frac{y_{0m}}{n_{0m}^v}\right){\left({\alpha }_{1m}\frac{y_{1m}}{n_{1m}^v}+{\alpha }_{0m}\frac{y_{0m}}{n_{0m}^v}\right)}^{-2}\\ =n^{-1}\sum _{m=1}^k\left[\left(\frac{1}{y_{1m}}-\frac{1}{n_{1m}^v}+\frac{1}{{\alpha }_{1m}}-\frac{n_{1m}^v}{{\alpha }_{1m}{y}_{1m}}\right)+\left(\frac{1}{y_{0m}}-\frac{1}{n_{0m}^v}+\frac{1}{{\alpha }_{0m}}-\frac{n_{0m}^v}{{\alpha }_{0m}{y}_{0m}}\right)\right]\\ \times {\left(\frac{{\alpha }_{1m}{y}_{1m}{\alpha }_{0m}{y}_{0m}}{{\alpha }_{1m}{y}_{1m}{n}_{0m}^v+{\alpha }_{0m}{y}_{0m}{n}_{1m}^v}\right)}^2.\end{array}$$

From $Ŵ$ and $Û$, we obtain an estimate of ${\sigma }_{b_1}^2$ as follows

$$\begin{array}{l}{\widehat{\sigma }}_{b_1}^2\\ =\frac{1}{{Ŵ}^2}\sum _{m=1}^k\left[\left(\frac{1}{y_{1m}}-\frac{1}{n_{1m}^v}+\frac{1}{{\alpha }_{1m}}\right)+\left(\frac{1}{y_{0m}}-\frac{1}{n_{0m}^v}+\frac{1}{{\alpha }_{0m}}\right)\right]{\left(\frac{{\alpha }_{1m}{y}_{1m}{\alpha }_{0m}{y}_{0m}}{{\alpha }_{1m}{y}_{1m}{n}_{0m}^v+{\alpha }_{0m}{y}_{0m}{n}_{1m}^v}\right)}^2\\ +\frac{1}{Ŵ}+\frac{1}{{Ŵ}^2}\sum _{m=1}^k\left[-\frac{n_{1m}^v}{{\alpha }_{1m}{y}_{1m}}-\frac{n_{0m}^v}{{\alpha }_{0m}{y}_{0m}}\right]{\left(\frac{{\alpha }_{1m}{y}_{1m}{\alpha }_{0m}{y}_{0m}}{{\alpha }_{1m}{y}_{1m}{n}_{0m}^v+{\alpha }_{0m}{y}_{0m}{n}_{1m}^v}\right)}^2\\ =\frac{1}{{Ŵ}^2}\sum _{m=1}^k\left[\left(\frac{1}{y_{1m}}-\frac{1}{n_{1m}^v}+\frac{1}{{\alpha }_{1m}}+\frac{1}{y_{0m}}-\frac{1}{n_{0m}^v}+\frac{1}{{\alpha }_{0m}}\right)\right]{\left(\frac{{\alpha }_{1m}{y}_{1m}{\alpha }_{0m}{y}_{0m}}{{\alpha }_{1m}{y}_{1m}{n}_{0m}^v+{\alpha }_{0m}{y}_{0m}{n}_{1m}^v}\right)}^2.\end{array}$$

The variance of ${\widehat{b}}_1$ is estimated by $\widehat{\text{Var}}\left({\widehat{b}}_1\right)=n^{-1}{\widehat{\sigma }}_{b_1}^2$. □