매우 중요
Suppose \(Y\sim \text{poisson}(\lambda)\). Then, it is possible that this model includes only some of explanatory variables, i.e., there exists some latent variables we cannot figure out. Because of this, for \(Y\sim \text{poisson}(\lambda)\), \(\lambda\) can vary as well.
Let \(E(\lambda)=\mu\), \(\text{Var}(\lambda)=\sigma^2\). Then, by the Tower property, \[\begin{eqnarray*} E(Y)&=&E_\lambda\{E(Y|\lambda) \}= E\{\lambda \}=\lambda,\\ \text{Var}(Y)&=&E_\lambda\{ \text{Var}(Y|\lambda) \}+\text{Var}_\lambda\{ E(Y|\lambda) \}\\ &=&E(\lambda)+\text{Var}(\lambda)=\mu+\sigma^2>\mu. \end{eqnarray*}\] Thus, we can say, there is overdispersion relative to the poisson model.
We have three approaches dealing with the overdispersion:
two-parameter family that permits \(\text{Var}(Y)>E(Y)\), for example, negative binomial.
quasi-likelihood approach assuming \(\text{Var}(Y)=cE(Y)\), without assuming a particular parametric distribution function for \(Y\).
mixture modeling approach. For example, let \(Y|\lambda\sim \text{poisson}(\lambda)\), \(\lambda\sim \Gamma(\alpha, \beta)\).
Suppose that given \(\lambda\), \(Y|\lambda\sim \text{poisson}(\lambda)\) and \(\lambda\sim\Gamma(k,k/\mu).\) Then, \[ f(\lambda;k,\mu)= \frac{(k/\mu)^k}{\Gamma(k)}\exp\left(-\frac{k\lambda }{\mu}\right), \mbox{ }\lambda>0. \]
We have \(E(\lambda)=\mu\) and \(\text{Var}(\lambda)=\mu^2/k\).
Then, \(Y\) has probability mass function \[\begin{eqnarray*} P(Y=y)&=& E[P(Y=y|\lambda)]\\ &=&\int_0^\infty \frac{e^{-\lambda}\lambda^y}{y!}\frac{(k/\mu)^k}{\Gamma(k)}\exp\left(-\frac{k\lambda }{\mu}\right)d\lambda\\ &=&\cdots= \frac{\Gamma(y+k)}{\Gamma(k)\Gamma(y+1)}\left(\frac{\mu}{\mu+k}\right)^y \left(1-\frac{\mu}{\mu+k}\right)^k, \end{eqnarray*}\] for \(y=0,1,2,\ldots.\). This is p.m.f of negative binomial distribution.
Note that
\[ E(Y)= \mu, \mbox{ } \text{Var}(Y)=\mu+\mu^2/k. \] As \(n\rightarrow \infty\), \(\text{Var}(\lambda)\rightarrow 0\) and \(\text{Var}(Y)\rightarrow \mu\) (Negative bionmial \(\rightarrow\) poisson).
Note that \[\begin{eqnarray*} P(Y=y; k, \mu)&=& \frac{\Gamma(y+k)}{\Gamma(k)\Gamma(y+1)}\left(\frac{\mu}{\mu+k}\right)^y \left(1-\frac{\mu}{\mu+k}\right)^k,\\ &=&\exp\left\{y \log \frac{\mu}{\mu+k} + k \log \frac{k}{\mu+k} +\log\frac{\Gamma(y+k)}{\Gamma(k)\Gamma(y+1)} \right\}\\ &=&\exp\left\{y \theta + k \log (1-e^\theta) +\log\frac{\Gamma(y+k)}{\Gamma(k)\Gamma(y+1)} \right\}. \end{eqnarray*}\] Thus, for known \(k\), the distribution of \(Y\) has the exponential dispersion form \(\exp\left\{\frac{y \theta- b(\theta)}{a(\phi)} +c(y;\phi) \right\}\) with
\(\theta=\log\frac{\mu}{\mu+k}\), \(b(\theta)=- k\log (1-e^\theta)\), \(a(\phi)=1\).
When \(Y\) has this distribution as p.m.f, we will write \(Y\sim \text{NegBin}(\mu,k)\).
Suppose now that \(Y_1,\ldots,Y_N\) are independent, where \(Y_i\sim \text{NegBin}(\mu_i,k)\) with \[ \log\mu_i=x_i'\beta. \] For fixed \(k\), we can estimate \(\beta\) by iterative reweighted least squares (IRLS). For fixed \(\beta\), the log-likelihood for \(k\) is \[\begin{eqnarray*} l(k)&=&\sum_i \left\{y_i \log \frac{\mu_i}{\mu_i+k} + k \log \frac{k}{\mu_i+k} +\log\Gamma(y_i+k)-\log{\Gamma(k)}- \log{\Gamma(y_i+1)} \right\}\\ &=& \sum_i\left\{\log \Gamma(y_i+k)-(y_i+k)\log(\mu_i+k)\right\}+n\{k\log k-\log \Gamma(k) \} +c(y;\mu). \end{eqnarray*}\] This expression can be maximized over \(k\)(differentiate, set to zero, and solve numerically for \(k\)). We then iterate these steps, alternatively maximizing over \(\beta\) for the current \(k\), and then \(k\) for the current \(\beta\). Because \[ E\left[ \frac{(Y_i-\mu_i)^2}{\mu_i}\right]= 1+\frac{\mu_i}{k}, \]
a reasonable starting value for \(k\) can be obtained by solving the equation \[ \frac{1}{n}X^2=1+\frac{\bar\mu}{k} \] where \(X^2\) is the Peason statistic obtained by fitting a poisson model, and \(\bar\mu\) is the average of the fitted mean for that model.