---
title: "Exercise 7.31"
author: "Per August Jarval Moen"
date: "26/10/2023"
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
\section*{Exercise 7.31}
We begin by loading the data:
```{r cache=TRUE}
url = "http://users.stat.ufl.edu/~aa/glm/data/Homicides.dat"
data = read.table(url, header=TRUE)
data = data[,c("race", "count")] #removing observation number
data$race=as.factor(data$race)
head(data)
```
Note that the data is un-grouped, while table 7.5 in the book shows the grouped representation. You can check that the dataset is the same as the table. For instance we have that:
```{r cache=TRUE}
sum( data[data$race==1,2]==0) #119, agrees with table
```
which agrees with table 7.5. 
\subsection*{a)}
We fit the Poisson GLM as specified in the exercise text:
```{r cache=TRUE}
poisson.model = glm(count~race, family=poisson(link="log"),data=data)
summary(poisson.model)
```
Interpretation: The rate of blacks knowing someone who has been murdered is estimated to be $\exp(\hat{\beta}) = 5.66$ times larger than the same rate for whites. 

\subsection*{b)}
There are arguably other characteristics that determine the rate at which one knows someone who is murdered, for instance socio-economic background, postal code, etc. Since these are not included in the model, the response variable $Y_i$ will show a larger variation/variance than the Poisson model suggests (remember that for the Poisson, the variance is equal to the mean/rate), making the Poisson model inadequate. 
\newline
\newline
To fit a Negative Binomial model, we use the MASS package:
```{r cache=TRUE}
library(MASS)
negbin.model = glm.nb(count~race,link="log",data=data)
summary(negbin.model)
```
Notice that the estimated parameters are almost the same as in the Poisson GLM. Also, note that the interpretation of $\hat{\beta}$ is the same as in the Poisson model, only that we swap the term "rate" with "mean". 
\newline
\newline
The estimated dispersion parameter is $\hat{\gamma} = \frac{1}{0.2023} \approx 4.94$, which is quite large. Note that R parametrizes the overdispersion parameter as $\theta = \frac{1}{\gamma}$.
\newline\newline
To formally compare the model fits, we use a likelihood-ratio test just as in exercise 5.30:
```{r cache=TRUE}
library(pscl)
odTest(negbin.model)
```
The p-value us practically zero, so we conclude there is strong evidence that the Negative Binomial model is a better fit (so the Poisson is overdispersed).
\subsection*{c)}
The Wald 95\% confidence intervals for $\beta$ are the usual confidence intervals based on asymptotic normality of the MLE: In both models, 
$$
\hat{\beta} \overset{d}{\approx}\text{N}(\beta, \hat{\text{se}}(\hat{\beta})),
$$
implying that a $95\%$ CI for $\beta$ using the Poisson GLM is
$$1.73314 \pm 1.96 \times 0.14657 = (1.44 \ , \ 2.02),$$
and a $95\%$ CI for $\beta$ using the Negative Binomial model is
$$\hat{\beta}\pm 1.96 \times 0.14657 = (1.27 \ , \ 2.2).$$
Exponentiating the above confidence intervals, we get that a 95\% confidence interval for the ratio of means using the Poisson GLM is
$$(4.25 \ , \ 7.55),$$
and a 95\% confidence interval for the ratio of means using the Negative Binomial model is
$$(3.54 \ , \ 9.03).$$
Which of these do we trust more? We have already seen strong evidence that the model fit of the Negative Binomial is better than that of the Poisson GLM, as the Poisson GLM is overdispersed. When GLMs are overdispersed, the model does not capture all of the variability of the data. The presence of overdispersion causes standard errors to be underestimated, and thus that confidence intervals are too narrow. A negative binomial model cannot be overdispersed (although it could of course have a poor model fit), so it is more robust.