Estimation of a zero-inflated Poisson regression model with missing covariates via nonparametric multiple imputation methods

Abstract

Zero-inflated Poisson (ZIP) regression is widely applied to model effects of covariates on an outcome count with excess zeros. In some applications, covariates in a ZIP regression model are partially observed. Based on the imputed data generated by applying the multiple imputation (MI) schemes developed by Wang and Chen (Ann Stat 37:490–517, 2009), two methods are proposed to estimate the parameters of a ZIP regression model with covariates missing at random. One, proposed by Rubin (in: Proceedings of the survey research methods section of the American Statistical Association, 1978), consists of obtaining a unified estimate as the average of estimates from all imputed datasets. The other, proposed by Fay (J Am Stat Assoc 91:490–498, 1996), consists of averaging the estimating scores from all imputed data sets to solve the imputed estimating equation. Moreover, it is shown that the two proposed estimation methods are asymptotically equivalent to the semiparametric inverse probability weighting method. A modified formula is proposed to estimate the variances of the MI estimators. An extensive simulation study is conducted to investigate the performance of the estimation methods. The practicality of the methodology is illustrated with a dataset of motorcycle survey of traffic regulations.

Appendix

1.1 Proof of Theorem 1

It can be obtained from the empirical CDF $\hat{F}(x|Y_i,\varvec{V}_i)$ in (10) that for $i=1,\dots ,n$ and $v=1,\dots ,M$,

$$\begin{aligned} E_{\hat{F}}(\tilde{S}_{iv}({\varvec{\theta }})|Y_i,\varvec{V}_i) =\sum _{k=1}^{n}\dfrac{\delta _kI(Y_k=Y_i,\varvec{V}_k=\varvec{V}_i)S_k({\varvec{\theta }})}{\sum _{r=1}^{n}\delta _rI(Y_r=Y_i,\varvec{V}_r=\varvec{V}_i)}. \end{aligned}$$

(18)

By using the expression of $U_v({\varvec{\theta }})$ in (11), the expression of $E_{\hat{F}}(\tilde{S}_{iv}({\varvec{\theta }})|Y_i,\varvec{V}_i)$ in (18), and the fact that

$$\begin{aligned} \sum _{i=1}^{n}(1-\delta _i)E_{\hat{F}}(\tilde{S}_{iv}({\varvec{\theta }})|Y_i,\varvec{V}_i) =\sum _{k=1}^n\delta _kS_k({\varvec{\theta }})\left[ \dfrac{1}{\hat{\pi }(Y_k,\varvec{V}_k)}-1\right] , \end{aligned}$$

(19)

$v=1,\ldots ,M$, we can have $E_{\hat{F}}(U_v({\varvec{\theta }})|\mathcal {O})=n^{-1/2}\sum _{i=1}^n[\delta _i/\hat{\pi }(Y_i,\varvec{V}_i)]S_i({\varvec{\theta }}) =U_w({\varvec{\theta }},\hat{\varvec{\pi }})$, $v=1,\ldots ,M$. Similarly, it can be shown that $E_{\hat{F}}(\partial {U}_v({\varvec{\theta }})/{\varvec{\theta }}|\mathcal {O}) =\partial {U}_w({\varvec{\theta }},\hat{\varvec{\pi }})/\partial {\varvec{\theta }}$ and, hence, $E(\partial {U}_v({\varvec{\theta }})/\partial {\varvec{\theta }}) =E(\partial {U}_w({\varvec{\theta }},\hat{\varvec{\pi }})/\partial {\varvec{\theta }})$, $v=1,\ldots ,M$.

Recall that $S_i^*({\varvec{\theta }})=E(S_i({\varvec{\theta }})|Y_i,\varvec{V}_i)$, $i=1,\dots ,n$. As given in (16), $U_{m2}({\varvec{\theta }})$ can be expressed as follows:

$$\begin{aligned} U_{m2}({\varvec{\theta }})= & {} \dfrac{1}{\sqrt{n}}\sum _{i=1}^{n}\left[ \delta _iS_i({\varvec{\theta }})+(1-\delta _i)S_i^*({\varvec{\theta }})\right] \nonumber \\&+\dfrac{1}{\sqrt{n}}\sum _{i=1}^{n}(1-\delta _i) \left[ \tilde{S}_{i}({\varvec{\theta }})-E_{\hat{F}}(\tilde{S}_{i1}({\varvec{\theta }})|\mathcal {O})\right] \nonumber \\&+\dfrac{1}{\sqrt{n}}\sum _{i=1}^{n}(1-\delta _i) \left[ E_{\hat{F}}(\tilde{S}_{i1}({\varvec{\theta }})\big |\mathcal {O})-S_i^*({\varvec{\theta }})\right] . \end{aligned}$$

(20)

Note that the second term of the expression of $U_{m2}({\varvec{\theta }})$ in (20) can be reformulated as

$$\begin{aligned} \dfrac{1}{\sqrt{n}}\sum _{i=1}^{n}(1-\delta _i) [\tilde{S}_{i}({\varvec{\theta }})-E_{\hat{F}}(\tilde{S}_{i1}({\varvec{\theta }})|\mathcal {O})] =O_p(M^{-1/2}), \end{aligned}$$

(21)

where $\tilde{S}_i({\varvec{\theta }})=\sum _{v=1}^{M}\tilde{S}_{iv}({\varvec{\theta }})/M$. The third term of the expression of $U_{m2}({\varvec{\theta }})$ in (20) can be rewritten as follows:

$$\begin{aligned}&\dfrac{1}{\sqrt{n}}\sum _{i=1}^{n}(1-\delta _i)\left[ E_{\hat{F}}(\tilde{S}_{i1}({\varvec{\theta }})|\mathcal {O})-S_i^*({\varvec{\theta }})\right] \nonumber \\&= \dfrac{1}{\sqrt{n}}\sum _{i=1}^{n}(1-\delta _i)\left[ \sum _{k=1}^{n}\dfrac{\delta _kS_k({\varvec{\theta }})I(Y_k=Y_i,\varvec{V}_k=\varvec{V}_i)}{\sum _{r=1}^{n}\delta _rI(Y_r=Y_i,\varvec{V}_r=\varvec{V}_i)} \nonumber \right. \\&\quad \left. -S_i^*({\varvec{\theta }})\sum _{k=1}^{n}\dfrac{\delta _kI(Y_k=Y_i,\varvec{V}_k=\varvec{V}_i)}{\sum _{r=1}^{n}\delta _rI(Y_r=Y_i,\varvec{V}_r=\varvec{V}_i)}\right] \nonumber \\&= \dfrac{1}{\sqrt{n}}\sum _{i=1}^{n}(1-\delta _i)\left[ \sum _{k=1}^{n}\dfrac{\delta _k[S_k({\varvec{\theta }})-S_i^*({\varvec{\theta }})] I(Y_k=Y_i,\varvec{V}_k=\varvec{V}_i)}{\sum _{r=1}^{n}\delta _rI(Y_r=Y_i,\varvec{V}_r=\varvec{V}_i)}\right] \nonumber \\&= \dfrac{1}{\sqrt{n}}\sum _{i=1}^{n}(1-\delta _i) \left[ \sum _{k=1}^n\dfrac{\delta _k[S_k({\varvec{\theta }})-S_i^*({\varvec{\theta }})] I(Y_k=Y_i,\varvec{V}_k=\varvec{V}_i)}{\sum _{s=1}^{n}I(Y_s=Y_i,\varvec{V}_s=\varvec{V}_i)}\right] \nonumber \\&\quad \times \,\left[ \dfrac{\sum _{s=1}^{n}I(Y_s=Y_i,\varvec{V}_s=\varvec{V}_i)}{\sum _{r=1}^{n}\delta _r I(Y_r=Y_i,\varvec{V}_r=\varvec{V}_i)}\right] \nonumber \\&= \dfrac{1}{\sqrt{n}}\sum _{i=1}^{n}\dfrac{(1-\delta _i)}{\hat{\pi }(Y_i,\varvec{V}_i)} \left[ \dfrac{\sum _{k=1}^{n}\delta _k[S_k({\varvec{\theta }})-S_k^*({\varvec{\theta }})]I(Y_k=Y_i,\varvec{V}_k=\varvec{V}_i)}{\sum _{s=1}^{n}I(Y_s=Y_i,\varvec{V}_s=\varvec{V}_i)}\right] \nonumber \\&=\dfrac{1}{\sqrt{n}}\sum _{k=1}^n\delta _k[S_k({\varvec{\theta }})-S_k^*({\varvec{\theta }})] \left\{ \sum _{i=1}^{n}\dfrac{(1-\delta _i)}{\hat{\pi }(Y_i,\varvec{V}_i)} \left[ \dfrac{I(Y_k=Y_i,\varvec{V}_k=\varvec{V}_i)}{\sum _{s=1}^{n}I(Y_s=Y_i,\varvec{V}_s=\varvec{V}_i)}\right] \right\} \nonumber \\&=\dfrac{1}{\sqrt{n}}\sum _{k=1}^{n}\delta _k[S_k({\varvec{\theta }})-S_k^*({\varvec{\theta }})] \left[ \dfrac{1-\hat{\pi }(Y_k,\varvec{V}_k)}{\hat{\pi }(Y_k,\varvec{V}_k)}\right] \nonumber \\&= \dfrac{1}{\sqrt{n}}\sum _{k=1}^{n}\delta _k\left[ \dfrac{1-\pi (Y_k,\varvec{V}_k)}{\pi (Y_k,\varvec{V}_k)}\right] \left[ S_k({\varvec{\theta }})-S_k^*({\varvec{\theta }})\right] +o_p(1). \end{aligned}$$

(22)

Hence, from (21) and (22), $U_{m2}({\varvec{\theta }})$ can be re-expressed as

$$\begin{aligned} U_{m2}({\varvec{\theta }})= & {} \dfrac{1}{\sqrt{n}}\sum _{i=1}^{n}\left\{ \dfrac{\delta _i}{\pi (Y_i,\varvec{V}_i)}S_i({\varvec{\theta }})+ \left[ 1-\dfrac{\delta _i}{\pi (Y_i,\varvec{V}_i)}\right] S_i^*({\varvec{\theta }})\right\} \nonumber \\&+ O_p(M^{-1/2})+o_p(1). \end{aligned}$$

(23)

Because the first term is the sum of independent and identically distributed random vectors, it can be shown by the multivariate central limit theorem that $U_{m2}({\varvec{\theta }}){\mathop {\rightarrow }\limits ^{d}}\mathcal {N}(\varvec{0},M({\varvec{\theta }},\varvec{\pi }))$ as $n,M\rightarrow \infty$, where $M({\varvec{\theta }},\varvec{\pi })$ is given in (15). In addition, $U_{m2}({\varvec{\theta }})$ in (23) can be expressed as

$$\begin{aligned} U_{m2}({\varvec{\theta }})= & {} U_w({\varvec{\theta }},\varvec{\pi }) +\dfrac{1}{\sqrt{n}}\sum _{i=1}^{n}\left[ 1-\dfrac{\delta _i}{\pi (Y_i,\varvec{V}_i)}\right] S_i^*({\varvec{\theta }})\nonumber \\&+ O_p(M^{-1/2}) +o_p(1). \end{aligned}$$

(24)

Because $\hat{\varvec{\theta }}_{m2}^{(M)}$ is the solution of $U_{m2}({\varvec{\theta }})=\varvec{0}$, it follows by a Taylor’s expansion of $U_{m2}(\hat{\varvec{\theta }}_{m2}^{(M)})$ at ${\varvec{\theta }}$ and the expression of $U_{m2}({\varvec{\theta }})$ in (23) that

$$\begin{aligned} \varvec{0}=U_{m2}(\hat{\varvec{\theta }}_{m2}^{(M)})= & {} \dfrac{1}{\sqrt{n}}\sum _{i=1}^{n}\left\{ \dfrac{\delta _i}{\pi (Y_i,\varvec{V}_i)}S_i({\varvec{\theta }}) +\left[ 1-\dfrac{\delta _i}{\pi (Y_i,\varvec{V}_i)}\right] S_i^*({\varvec{\theta }})\right\} \\&- G({\varvec{\theta }},\varvec{\pi })\sqrt{n}(\hat{\varvec{\theta }}_{m2}^{(M)}-{\varvec{\theta }})+O_p(M^{-1/2})+o_p(1), \end{aligned}$$

where $G({\varvec{\theta }},\varvec{\pi })=E[-\partial {U}_w({\varvec{\theta }},\varvec{\pi })/(\sqrt{n}\partial {\varvec{\theta }})]$. Therefore, it can be obtained that

$$\begin{aligned} \sqrt{n}(\hat{\varvec{\theta }}_{m2}^{(M)}-{\varvec{\theta }})= & {} \dfrac{1}{\sqrt{n}}G^{-1}({\varvec{\theta }},\varvec{\pi })\left\{ \sum _{i=1}^n \left[ \dfrac{\delta _i}{\pi (Y_i,\varvec{V}_i)}S_i({\varvec{\theta }})+\left( 1-\dfrac{\delta _i}{\pi (Y_i,\varvec{V}_i)}\right) S_i^*({\varvec{\theta }})\right] \right\} \\&+\, O_p(M^{-1/2})+o_p(1). \end{aligned}$$

This implies that $\sqrt{n}(\hat{\varvec{\theta }}_{m2}^{(M)}-{\varvec{\theta }}){\mathop {\rightarrow }\limits ^{d}}\mathcal {N}(\varvec{0},\Delta _m({\varvec{\theta }}))$ as $n,M\rightarrow \infty$, where $\Delta _m({\varvec{\theta }})=G^{-1}({\varvec{\theta }},\varvec{\pi })M({\varvec{\theta }},\varvec{\pi })[G^{-1}({\varvec{\theta }},\varvec{\pi })]^T$.

Let $\hat{\varvec{\theta }}_v$ be the solution to the estimating equations $U_v({\varvec{\theta }})=\varvec{0}$. We have by a Taylor’s expansion of $U_v(\hat{\varvec{\theta }}_v)$ at ${\varvec{\theta }}$ that

$$\begin{aligned} \varvec{0}=U_v(\hat{\varvec{\theta }}_v) =U_v({\varvec{\theta }})-G({\varvec{\theta }},\varvec{\pi })\sqrt{n}(\hat{\varvec{\theta }}_v-{\varvec{\theta }})+o_p(1). \end{aligned}$$

Hence, it follows that $\sqrt{n}(\hat{\varvec{\theta }}_v-{\varvec{\theta }})=G^{-1}({\varvec{\theta }},\varvec{\pi })U_v({\varvec{\theta }})+o_p(1)$. Because $\hat{\varvec{\theta }}_{m1}^{(M)}=\sum _{v=1}^{M}\hat{\varvec{\theta }}_v/M$, using the above result and the expressions for $U_{m2}({\varvec{\theta }})$ in (13) and (23), we can have

$$\begin{aligned} \sqrt{n}(\hat{\varvec{\theta }}_{m1}^{(M)}-{\varvec{\theta }})= & {} G^{-1}({\varvec{\theta }},\varvec{\pi })\left( \dfrac{1}{M}\sum _{v=1}^MU_v({\varvec{\theta }})\right) +o_p(1)\nonumber \\= & {} \dfrac{1}{\sqrt{n}}G^{-1}({\varvec{\theta }},\varvec{\pi })\sum _{i=1}^n\left\{ \dfrac{\delta _i}{\pi (Y_i,\varvec{V}_i)} S_i({\varvec{\theta }})+\left[ 1-\dfrac{\delta _i}{\pi (Y_i,\varvec{V}_i)}\right] S_i^*({\varvec{\theta }})\right\} \nonumber \\&+\, O_p(M^{-1/2})+o_p(1), \end{aligned}$$

(25)

and it is shown easily that $\sqrt{n}(\hat{\varvec{\theta }}_{m1}^{(M)}-{\varvec{\theta }}){\mathop {\rightarrow }\limits ^{d}}{N}(\varvec{0},\Delta _{m}({\varvec{\theta }}))$ as $n,M\rightarrow \infty$.

1.2 Proof of Theorem 2

Because $\hat{\varvec{\theta }}_{m2}^{(M)}$ is the solution of $U_{m2}({\varvec{\theta }})=M^{-1}\sum _{v=1}^{M}U_v({\varvec{\theta }})=\varvec{0}$, a Taylor’s expansion of $U_{m2}(\hat{\varvec{\theta }}_{m2}^{(M)})$ at ${\varvec{\theta }}$ can lead to $\varvec{0}=U_{m2}(\hat{\varvec{\theta }}_{m2}^{(M)}) =M^{-1}\sum _{v=1}^MU_v({\varvec{\theta }})-G({\varvec{\theta }},\varvec{\pi })\sqrt{n}(\hat{\varvec{\theta }}_{m2}^{(M)}-{\varvec{\theta }})+o_p(1)$, which implies that

$$\begin{aligned} \sqrt{n}(\hat{\varvec{\theta }}_{m2}^{(M)}-{\varvec{\theta }}) =G^{-1}({\varvec{\theta }},\varvec{\pi })\left( \dfrac{1}{M}\sum _{v=1}^{M}U_v({\varvec{\theta }})\right) +o_p(1). \end{aligned}$$

(26)

Thus, it follows from (25) and (26) that $\sqrt{n}(\hat{\varvec{\theta }}_{m2}^{(M)}-\hat{\varvec{\theta }}_{m1}^{(M)})=o_p(1)$. This shows that $\sqrt{n}(\hat{\varvec{\theta }}_{m2}^{(M)}-\hat{\varvec{\theta }}_{m1}^{(M)})$ converges in probability to $\varvec{0}$ as $n,M\rightarrow \infty$.

Next, we show that the semiparametric IPW estimator and the second MI-type estimator are asymptotically equivalent. $U_{m2}({\varvec{\theta }})$ can be expressed as

$$\begin{aligned} U_{m2}({\varvec{\theta }})= & {} \dfrac{1}{\sqrt{n}}\sum _{i=1}^{n} \left\{ \delta _iS_i({\varvec{\theta }})+(1-\delta _i)E_{\hat{F}}(\tilde{S}_{i1}({\varvec{\theta }})\big |Y_i,\varvec{V}_i) \right. \\&\left. +\,(1-\delta _i)\left[ \tilde{S}_i({\varvec{\theta }})-E_{\hat{F}}(\tilde{S}_{i1}({\varvec{\theta }})|Y_i,\varvec{V}_i)\right] \right\} . \end{aligned}$$

Using the fact given in (19), $n^{-1/2}\sum _{i=1}^n[\delta _iS_i({\varvec{\theta }})+(1-\delta _i)E_{\hat{F}}(\tilde{S}_{i1}({\varvec{\theta }})|Y_i,\varvec{V}_i)]$ can be expressed as

$$\begin{aligned}&\dfrac{1}{\sqrt{n}} \sum _{i=1}^{n}[\delta _iS_i({\varvec{\theta }})+(1-\delta _i)E_{\hat{F}}(\tilde{S}_{i1}({\varvec{\theta }})\big |Y_i,\varvec{V}_i)] \\&\quad = \dfrac{1}{\sqrt{n}}\sum _{i=1}^{n} \delta _iS_i({\varvec{\theta }})+\dfrac{1}{\sqrt{n}}\sum _{k=1}^{n}\delta _kS_k({\varvec{\theta }}) \left[ \dfrac{1}{\hat{\pi }(Y_k,\varvec{V}_k)}-1\right] =U_w({\varvec{\theta }},\hat{\varvec{\pi }}). \end{aligned}$$

Hence it can be obtained that

$$\begin{aligned} U_{m2}({\varvec{\theta }})=U_w({\varvec{\theta }},\hat{\varvec{\pi }})+\dfrac{1}{\sqrt{n}}\sum _{i=1}^n(1-\delta _i) \left[ \tilde{S}_i({\varvec{\theta }})-E_{\hat{F}}(\tilde{S}_{i1}({\varvec{\theta }})|Y_i,\varvec{V}_i)\right] . \end{aligned}$$

Recall $\mathcal {B}({\varvec{\theta }};Y_i,\varvec{V}_i) =M^{-1/2}\sum _{v=1}^M\left[ \tilde{S}_{iv}({\varvec{\theta }})-E_{\hat{F}}(\tilde{S}_{i1}({\varvec{\theta }})|Y_i,\varvec{V}_i)\right]$, $i=1,\ldots ,n$. Because

$$\begin{aligned} \dfrac{1}{\sqrt{n}}\sum _{i=1}^{n}(1-\delta _i) \left[ \tilde{S}_i({\varvec{\theta }})-E_{\hat{F}}(\tilde{S}_{i1}({\varvec{\theta }})|Y_i,\varvec{V}_i)\right] M^{-1/2}\dfrac{1}{\sqrt{n}}\sum _{i=1}^{n}(1-\delta _i)\mathcal {B}({\varvec{\theta }};Y_i,\varvec{V}_i) \end{aligned}$$

and $(1-\delta _i)\mathcal {B}({\varvec{\theta }};Y_i,\varvec{V}_i)$ are independent and identically distributed random vectors with mean $\varvec{0}$ and covariance matrix $E\{(1-\delta _1)[\mathcal {B}({\varvec{\theta }},Y_1,\varvec{V}_1)]^{\otimes 2}\}$, it implies by the multivariate central limit theorem that $n^{-1/2}\sum _{i=1}^{n}(1-\delta _i)\mathcal {B}({\varvec{\theta }};Y_i,\varvec{V}_i)=O_p(1)$ and, hence,

$$\begin{aligned} U_{m2}({\varvec{\theta }})-U_w({\varvec{\theta }},\hat{\varvec{\pi }})= & {} \dfrac{1}{\sqrt{n}}\sum _{i=1}^{n}(1-\delta _i) \left[ \tilde{S}_i({\varvec{\theta }})-E_{\hat{F}}(\tilde{S}_{i1}({\varvec{\theta }})\big |Y_i,\varvec{V}_i)\right] \nonumber \\= & {} M^{-1/2}O_p(1)=O_p(M^{-1/2}). \end{aligned}$$

(27)

Let $\hat{\varvec{\theta }}_{m2}^{(M)}$ be the solution of $U_{m2}({\varvec{\theta }})=\varvec{0}$. Because by a Taylor’s expansion of $U_{m2}(\hat{\varvec{\theta }}_{m2}^{(M)})$ at ${\varvec{\theta }}$ and $U_w(\hat{\varvec{\theta }}_{ws},\hat{\varvec{\pi }})$ at ${\varvec{\theta }}$, respectively, we can have that

$$\begin{aligned} \varvec{0}=U_{m2}(\hat{\varvec{\theta }}_{m2}^{(M)}) =U_{m2}({\varvec{\theta }})-G({\varvec{\theta }},\varvec{\pi })\sqrt{n}(\hat{\varvec{\theta }}_{m2}^{(M)}-{\varvec{\theta }})+o_p(1) \end{aligned}$$

and

$$\begin{aligned} \varvec{0}=U_w(\hat{\varvec{\theta }}_{ws},\hat{\varvec{\pi }}) =U_{w}({\varvec{\theta }},\hat{\varvec{\pi }})-G({\varvec{\theta }},\varvec{\pi })\sqrt{n}(\hat{\varvec{\theta }}_{ws}-{\varvec{\theta }})+o_p(1), \end{aligned}$$

it can be shown from (27) that

$$\begin{aligned} \sqrt{n}(\hat{\varvec{\theta }}_{m2}^{(M)}-\hat{\varvec{\theta }}_{ws}) G^{-1}({\varvec{\theta }},\varvec{\pi })\left[ U_{m2}({\varvec{\theta }})-U_w({\varvec{\theta }},\hat{\varvec{\pi }})\right] +o_p(1) =o_p(1)+O_p(M^{-1/2}). \end{aligned}$$

Therefore, it follows that $\sqrt{n}(\hat{\varvec{\theta }}_{m2}^{(M)}-\hat{\varvec{\theta }}_{ws})$ converges in probability to $\varvec{0}$ as $n,M\rightarrow \infty$.