Chapter 28 A Generalized Linear Model for Poison Response Data

For all \(i = 1, \ldots, n\), \(y_i \sim \text{Possion}(\lambda_i)\), \(\log(\lambda_i) = x_i'\beta\). The Poisson log likelihood is \[ \ell(\beta\mid y) = \sum_{i=1}^n \left[y_ix_i'\beta - exp(x_i'\beta) - \log(y_i!) \right] \] The \(\ell(\beta\mid y)\) can be maximized using Fisher’s scoring method to obtain the MLE.

Let \(\lambda = \exp(x'\beta)\) and \(\tilde{\lambda} = exp(\tilde{x}'\beta)\) where \(\tilde{x} = [x_1, \ldots, x_{j-1}, x_{j} + 1, x_{j+1}, \ldots, x_p ]'\), we have \(\tilde{\lambda}/\lambda = \exp(\beta_j)\). This means that all other explanatory variables held constant, the mean response at \(x_j + 1\) is \(\exp(\beta_j)\) times the mean response at \(x_j\).

o = glm(y ~ x, family = poisson(link = "log"))
summary(o) 

# likelihood ratio test 
anova(o, test = "Chisq")

Lack of Fit: Under saturated model, \(\lambda_i = y_i\). Then the likelihood ratio statistic for testing the Poisson GLM as the reduced model vs. the saturated model as the full model is \[ 2 \sum_{i = 1}^n \left[y_i\log\left(\frac{y_i}{\hat\lambda_i}\right) - (y_i - \hat\lambda_i) \right] \] which is the Deviance Statistic for the Poisson case.

The deviance residuals are given by \[ d_i \equiv \text{sign}(y_i - \hat\lambda_i)\sqrt{2\left[y_i\log\left(\frac{y_i}{\hat\lambda_i} \right) - (y_i - \hat\lambda_i)\right]} \] The Pearson’s Chi-square statistic is \(X^2 = \sum_{i=1}^{n}\left(\frac{y_{i}-\widehat{E}\left(y_{i}\right)}{\sqrt{\widehat{\operatorname{Var}}\left(y_{i}\right)}}\right)^{2} = \sum_{i=1}^n \left(\frac{y_i - \hat\lambda_i}{\sqrt{\hat\lambda_i}}\right)^2\). The Pearson residure \(r_i = (y_i - \hat\lambda_i)/\sqrt{\hat\lambda_i}\).

d = resid(o, type = "deviance")
r = resid(o, type = "pearson")

For the Poisson case, \(Var(y) = E(y) = \lambda\) is a function of \(E(y)\). If either the Deviance Statistic or the Pearson Chi-Square Statistic suggests a lack of fit that cannot be explained by other reasons (e.g., poor model for the mean or a few extreme outliers), overdispersion may be the problem.

Quasi-likelihood: Suppose \(Var(y_i) = \phi \lambda_i\) for some unknown dispersion parameter \(\phi > 1\). \(\phi\) can be estimated by \(\hat\phi = \sum_{i = 1}^n d_i^2/(n-p)\) or \(\hat\phi = \sum_{i=1}^n r_i^2/(n-p)\).

• The estimated variance of \(\hat\beta\) is multiplied by \(\hat\phi\).
• For Wald type inferences, the standard normal null distribution is replaced by \(t\) with \(n-p\) degrees of freedom.
• Any test statistic \(T\) that was assumed \(\chi_q^2\) under \(H_0\) is replaced with \(T/(q\hat\phi)\) and compared to an \(F\) distribution with \(q\) and \(n-p\) degrees of freedom.

# estimates of the dispersion parameter 
deviance(o)/df.residual(o)

sum(r^2)/df.residual(o)