Zero-and-one-inflated Poisson regression model

Abstract

In this paper, a zero-and-one-inflated Poisson (ZOIP) regression model is proposed. The maximum likelihood estimation (MLE) and Bayesian estimation for this model are investigated. Three estimation methods of the ZOIP regression model are obtained based on data augmentation method which is expectation-maximization (EM) algorithm, generalized expectation-maximization (GEM) algorithm and Gibbs sampling respectively. A simulation study is conducted to assess the performance of the proposed estimation for various sample sizes. Finally, an accidental deaths data set is analyzed to illustrate the practicability of the proposed method.

Appendix

1.1 Appendix A.1

Proof of Theorem 3.1

The proof depends on the smoothness of $Q(\varvec{\eta },\varvec{\eta }^{(k)}) = E[\ell _{c}(\varvec{\eta }|{\varvec{Y}}, {\varvec{B}}_{1},{\varvec{B}}_{2})|{\varvec{Y}}, \varvec{\eta }^{(k)}]$. We rewrite

$$\begin{aligned} Q(\varvec{\eta },\varvec{\eta }^{(k)}) = Q_{1}(\varvec{\eta },\varvec{\eta }^{(k)}) + Q_{2}(\varvec{\eta },\varvec{\eta }^{(k)}) + Q_{3}(\varvec{\eta },\varvec{\eta }^{(k)}), \end{aligned}$$

where

$$\begin{aligned}&Q_{1}(\varvec{\eta },\varvec{\eta }^{(k)})=\sum _{i=1}^{n} E\bigg [(1-B_{1i})\big (Y_{i}{\varvec{z}}_{i}^{\mathrm {T}} {\varvec{\beta }}-\mathrm {e}^{{\varvec{z}}_{i}^{\mathrm {T}} {\varvec{\beta }}}\big )\big |{\varvec{Y}},\varvec{\eta }^{(k)}\bigg ],\\&Q_{2}(\varvec{\eta },\varvec{\eta }^{(k)})=\sum _{i=1}^{n} E\bigg [B_{1i}\varvec{\omega }_{1i}^{\mathrm {T}} {\varvec{\gamma }}_{1}-\ln \big (1+ \mathrm {e}^{\varvec{\omega }_{1i}^{\mathrm {T}} {\varvec{\gamma }}_{1}}\big )-(1-B_{1i})\ln (Y_{i}!)| {\varvec{Y}},\varvec{\eta }^{(k)}\bigg ], \end{aligned}$$

and

$$\begin{aligned} Q_{3}(\varvec{\eta },\varvec{\eta }^{(k)})=\sum _{i=1}^{n} E\bigg [B_{2i}\varvec{\omega }_{2i}^{\mathrm {T}} {\varvec{\gamma }}_{2}-\ln \big (1+\mathrm {e}^ {\varvec{\omega }_{2i}^{\mathrm {T}}{\varvec{\gamma }}_{2}}\big )\bigg | {\varvec{Y}},\varvec{\eta }^{(k)}\bigg ]. \end{aligned}$$

$\square$

Simple calculation yields

$$\begin{aligned} \begin{aligned}&Q_{1}(\varvec{\eta },\varvec{\eta }^{(k)})\\&\quad = \sum _{i=1}^{n}I\{Y_{i}=0\}\left[ \frac{(1-p_{1i}^{(k)}) P(V=0)}{p_{1i}^{(k)}p_{2i}^{(k)}+(1-p_{1i}^{(k)}) P(V_{i}^{(k)}=0)}\right] \bigg (Y_{i}{\varvec{z}}_{i}^{\mathrm {T}} {\varvec{\beta }}-\mathrm {e}^{{\varvec{z}}_{i}^{\mathrm {T}} {\varvec{\beta }}}\bigg )\\&\qquad +\sum _{i=1}^{n}I\{Y_{i}=1\}\left[ \frac{(1-p_{1i}^{(k)}) P(V=1)}{p_{1i}^{(k)}(1-p_{2i}^{(k)})+(1-p_{1i}^{(k)}) P(V_{i}^{(k)}=1)}\right] \bigg (Y_{i}{\varvec{z}}_{i}^{\mathrm {T}} {\varvec{\beta }}-\mathrm {e}^{{\varvec{z}}_{i}^{\mathrm {T}} {\varvec{\beta }}}\bigg )\\&\qquad +\sum _{i=1}^{n}I\{Y_{i}\ge 2\} \big (Y_{i}{\varvec{z}}_{i}^{\mathrm {T}} {\varvec{\beta }}-\mathrm {e}^{{\varvec{z}}_{i}^{\mathrm {T}} {\varvec{\beta }}}\big ), \end{aligned} \end{aligned}$$

which is continuous in $\varvec{\eta }$ and $\varvec{\eta }^{(k)}$. A similar conclusion can be applied to $Q_{2}(\varvec{\eta },\varvec{\eta }^{(k)})$ and $Q_{3}(\varvec{\eta },\varvec{\eta }^{(k)})$. It is easy to find that the Eq. (10) in Wu (1983) for $Q(\varvec{\eta },\varvec{\eta }^{(k)})$ is correct. According to the GEM algorithm, it is easy to that

$$\begin{aligned}&\ell _{c}\bigg ({\varvec{\gamma }}_{1}^{(k+1)}\bigg |{\varvec{Y}}, {\varvec{B}}_{1}^{(k+1)},{\varvec{B}}_{2}^{(k+1)}\bigg )-\ell _{c} \bigg ({\varvec{\gamma }}_{1}^{(k)}\bigg |{\varvec{Y}},{\varvec{B}}_{1}^{(k+1)}, {\varvec{B}}_{2}^{(k+1)}\bigg )\nonumber \\&\quad \ge \bigg |\bigg |{\varvec{\gamma }}_{1}^{(k+1)} -{\varvec{\gamma }}_{1}^{(k)}\bigg |\bigg |^{2}/2, \end{aligned}$$

(13)

$$\begin{aligned} \nonumber \\&\ell _{c}\bigg ({\varvec{\gamma }}_{2}^{(k+1)}\bigg |{\varvec{Y}}, {\varvec{B}}_{1}^{(k+1)},{\varvec{B}}_{2}^{(k+1)}\bigg )-\ell _{c} \bigg ({\varvec{\gamma }}_{2}^{(k)}\bigg |{\varvec{Y}},{\varvec{B}}_{1}^{(k+1)}, {\varvec{B}}_{2}^{(k+1)}\bigg )\nonumber \\&\quad \ge \bigg |\bigg |{\varvec{\gamma }}_{2}^{(k+1)} -{\varvec{\gamma }}_{2}^{(k)}\bigg |\bigg |^{2}/2 \end{aligned}$$

(14)

and

$$\begin{aligned} \ell _{c}\bigg ({\varvec{\beta }}^{(k+1)}\bigg |{\varvec{Y}}, {\varvec{B}}_{1}^{(k+1)},{\varvec{B}}_{2}^{(k+1)}\bigg )&-\ell _{c}\bigg ({\varvec{\beta }}^{(k)}\bigg |{\varvec{Y}}, {\varvec{B}}_{1}^{(k+1)},{\varvec{B}}_{2}^{(k+1)}\bigg )\nonumber \\&\ge \bigg |\bigg |\bigg ({\varvec{\beta }}^{(k+1)}-{\varvec{\beta }}^{(k)}\bigg )\bigg |\bigg |^{2}/2 \end{aligned}$$

(15)

respectively. And according to Eqs. (13), (14), and (15), it is easy to find that $||\varvec{\eta }^{(k+1)}-\varvec{\eta }^{(k)}||\rightarrow 0$ as $k \rightarrow \infty$. So according to Theorem 2 and Theorem 5 in Wu (1983), the conclusion is obvious.

1.2 Appendix A.2

The full conditional distribution of ${\varvec{\beta }}, {\varvec{\gamma }}_1$ and ${\varvec{\gamma }}_2$ with Normal prior. After the data augmentation step, we have the joint posterior of ${\varvec{\beta }}$, ${\varvec{\gamma }}_{1}$ and ${\varvec{\gamma }}_{2}$ as follows:

$$\begin{aligned} \begin{aligned}&{\pi [{\varvec{\beta }},{\varvec{\gamma }}_{1}, {\varvec{\gamma }}_{2}|\text {rest},{\varvec{Y}}]} \propto \exp \left\{ \left( \sum _{i=1}^{n}(1-B_{1i})V_{i} {\varvec{z}}_{i}\right) ^{\mathrm {T}}{\varvec{\beta }}\right\} \exp \\&\quad \left\{ \left( \sum _{i=1}^{n}B_{1i}\varvec{\omega _{1i}}\right) ^{\mathrm {T}}\varvec{\gamma _{1}}\right\} \exp \left\{ \left( \sum _{i=1}^{n}B_{2i}\varvec{\omega }_{2i}\right) ^{\mathrm {T}}{\varvec{\gamma }}_{2}\right\} \\&\quad \prod _{i=1}^{n}\frac{\exp \{-(1-B_{1i})\exp {({\varvec{z}}_{i}^ {\mathrm {T}}{\varvec{\beta }})}\}}{ (1+\exp {(\varvec{\omega }_{1i}^{\mathrm {T}} {\varvec{\gamma }}_{1})})(1+\exp {(\varvec{\omega }_{2i} ^{\mathrm {T}}{\varvec{\gamma }}_{2})})} \times \Pi ({\varvec{\beta }},{\varvec{\gamma }}_{1}, {\varvec{\gamma }}_{2}), \end{aligned} \end{aligned}$$

where the $\Pi ({\varvec{\beta }},{\varvec{\gamma }}_{1},{\varvec{\gamma }}_{2})$ is the prior of the parameters ${\varvec{\beta }}, {\varvec{\gamma }}_{1}$ and ${\varvec{\gamma }}_{2}$. Let $N({\varvec{\beta }}_{0},\sigma _{\beta }{\varvec{I}}_{q})$, $N({\varvec{\gamma }}_{01},\sigma _{\gamma _{1}}{\varvec{I}}_{r1})$ and $N({\varvec{\gamma }}_{02},\sigma _{\gamma _{2}}{\varvec{I}}_{r2})$ be the priors for parameters ${\varvec{\beta }}, {\varvec{\gamma }}_{1}$ and ${\varvec{\gamma }}_{2}$ respectively and assume that they are mutually independent. The full conditional distributions of ${\varvec{\beta }}$, ${\varvec{\gamma }}_{1}$ and ${\varvec{\gamma }}_{2}$ are not standard distributions. Their densities are as follows:

$$\begin{aligned}&{\pi [{\varvec{\beta }}|\text {rest},{\varvec{Y}}]}\propto (\sigma _{\beta })^{-q}\exp \left\{ -\left( \frac{1}{2}\sigma _{\beta }^{-2}\right) ({\varvec{\beta }}-{\varvec{\beta }}_{0})^{\mathrm {T}} ({\varvec{\beta }}-{\varvec{\beta }}_{0})\right\} \nonumber \\&\quad \times \exp \left\{ \left( \sum _{i=1}^{n}(1-B_{1i})V_{i}{\varvec{z}}_{i} \right) ^{\mathrm {T}} {\varvec{\beta }}\right\} \prod _{i=1}^{n}\exp {\left\{ -(1-B_{1i})\exp {\bigg ({\varvec{z}}_{i}^ {\mathrm {T}}{\varvec{\beta }}\big )}\right\} }, \end{aligned}$$

(16)

$$\begin{aligned}&{\pi [{\varvec{\gamma }}_{1}|\text {rest},{\varvec{Y}}]}\propto (\sigma _{\gamma _{1}})^{-r_{1}}\exp \left\{ -\left( \frac{1}{2} \sigma _{\gamma _{1}} ^{-2}\right) ({\varvec{\gamma }}_{1}-{\varvec{\gamma }}_{01})^{\mathrm {T}} ({\varvec{\gamma }}_{1}-{\varvec{\gamma }}_{01})\right\} \nonumber \\&\quad \times \exp \left\{ \left( \sum _{i=1}^{n}B_{1i}\varvec{\omega }_{1i} \right) ^{\mathrm {T}} {\varvec{\gamma }}_{1}\right\} \prod _{i=1}^{n}\bigg (1+\exp {\bigg (\varvec{\omega }_{1i}^{\mathrm {T}} {\varvec{\gamma }}_{1}\bigg )}\bigg )^{-1} \end{aligned}$$

(17)

and

$$\begin{aligned} \begin{aligned}&{\pi [{\varvec{\gamma }}_{2}|\text {rest},{\varvec{Y}}]}\propto (\sigma _{\gamma _{2}})^{-r_{2}}\exp \left\{ -\left( \frac{1}{2}\sigma _{\gamma _{2}}^ {-2}\right) ({\varvec{\gamma }}_{2}-{\varvec{\gamma }}_{02})^{\mathrm {T}} ({\varvec{\gamma }}_{2}-{\varvec{\gamma }}_{02})\right\} \\&\quad \times \exp \left\{ \left( \sum _{i=1}^{n}B_{2i}\varvec{\omega }_{2i}\right) ^{\mathrm {T}} {\varvec{\gamma }}_{2}\right\} \prod _{i=1}^{n}\bigg (1+\exp {\bigg (\varvec{\omega }_{2i}^{\mathrm {T}} {\varvec{\gamma }}_{2}\bigg )}\bigg )^{-1}. \end{aligned} \end{aligned}$$

(18)

It is also easy to get the second order partial derivative of ${\pi [{\varvec{\beta }}|\text {rest},{\varvec{Y}}]}$,

$$\begin{aligned} \frac{\partial ^{2} {\pi [{\varvec{\beta }}|\text {rest},{\varvec{Y}}]}}{\partial {\varvec{\beta }}\partial {\varvec{\beta }}^{T}} =\frac{1}{2}\sigma _{\beta }^{-2}{\varvec{I}}_{q}+\sum _{i=1}^{n} (1-B_{1i})\exp ({\varvec{z}}_{i}^{T}{\varvec{\beta }}) {\varvec{z}}_{i}{\varvec{z}}_{i}^{T}>0. \end{aligned}$$

So ${\pi [{\varvec{\beta }}|\text {rest},{\varvec{Y}}]}$ are log-concave. The log-concave property of conditional densities ${\pi [{\varvec{\gamma }}_{1}|\text {rest},{\varvec{Y}}]}$ and ${\pi [{\varvec{\gamma }}_{2}|\text {rest},{\varvec{Y}}]}$ can be proved similarly. So ARS can be used to sample ${\varvec{\beta }}$, ${\varvec{\gamma }}_{1}$ and ${\varvec{\gamma }}_{2}$ from their respective full conditional distributions.