---
title: "GLM Examples"
author: "Nancy"
date: "01/12/2021"
output:
  pdf_document: default
  html_document: default
---

I'm sure there are nicer ways to set things up, but I'm using a default version, except I added `scipen`, to force the output not to use exponential notation, which I hate. 

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
options(scipen=10, digits=4)
```

```{r libraries, include=FALSE}
library(faraway)
library(ggplot2)
library(dplyr)
library(kableExtra)
library(SMPracticals)
```

## Example 1: binomial data (Exercise 2.1 in ELM)

Faraway asks first to fit a model with all interactions, so here it is:

```{r Example 1}
data(esoph)
head(esoph)
help(esoph)
xtabs(ncases ~ agegp+tobgp, data = esoph)
```
There could be more elegant ways to tabulate the data, but this is at least feasible to look at. In the help file for `xtabs` this dataset is used as an example. And I learned about "flat tables". 
```{r}
ftable(xtabs(cbind(ncases,ncontrols) ~ agegp, data = esoph))
ftable(xtabs(cbind(ncases,ncontrols) ~ agegp + tobgp, data = esoph))
```

And now I'll try the full model as requested in (a): 
```{r}
farmod1 <- glm(cbind(ncases, ncontrols) ~ agegp * tobgp * alcgp, data = esoph, family = binomial)

summary(farmod1)
```
So we probably didn't really want this. We've fit a parameter to every observation, not helpful. What is with all the .L, .C etc? 

```{r}
is.ordered(esoph$agegp)
```
Ordered factors use polynomial contrasts to fit the various levels. This is a kind of compromise between treating the values as continuous and treating them as completely unordered. In the `help` file for this dataset you'll see models with covariates `unclass(tobgp)` and `unclass(alcgp)`. 

Now Faraway suggests using backward elimination to simplify the model "as far as is reasonable". 

```{r}
farmod2 <- step(farmod1, direction = "backward")
summary(farmod2)
```
This eliminated all the interactions, but still used orthogonal polynomials for the variables. And up to 5th degree for age! Doesn't quite make sense. We could either treat the variables as continuous, or as unordered factors. 
```{r}
farmod3 <- glm(cbind(ncases, ncontrols) ~ unclass(agegp) + unclass(alcgp) + unclass(tobgp), family = binomial, data = esoph)
summary(farmod3)
```
This seems a bit too simple, treating each covariate as continuous and just fitting a linear model. Myabe we could try linear in age, but with tobacco and alcohol as factors. Below are several versions, they can be compared using `anova`, or their residual deviances can be compared. 
```{r}
esoph2 <- mutate(esoph,age = factor(agegp,ordered=FALSE), tob = factor(tobgp, ordered = F), alc = factor(alcgp, ordered = F) )
# new dataset has unordered factors 
farmod4 <- glm(cbind(ncases, ncontrols) ~ unclass(agegp) + tob*alc, data = esoph2, family = binomial)
farmod5 <- glm(cbind(ncases, ncontrols) ~ unclass(agegp) + tob + alc, data = esoph2, family = binomial)
farmod6 <- glm(cbind(ncases, ncontrols) ~ agegp + tob + alc, data = esoph2, family = binomial)
```

## Example 2 Poisson 
This example in Chapter 3 of ELM is used to introduce Poisson regression. Faraway first considers ordinary linear regression, and notes that the residuals don't seem to have constant variance. Back in the day we'd take logs or square root of the response and then fit linear regression again, but a more modern approach is to fit a Poisson glm.

The residual deviance in a Poisson model is a test of fit of the model, as it is with a Binomial.  The second fit below allows an inflation factor on the variance, while still using the simplicity of the Poisson for fitting. 
```{r poisson gala}
data(gala)
head(gala)
help(gala)

gal.glm <- glm(Species ~ . - Endemics,
             family = poisson,
             data = gala)

summary(gal.glm)

gal.quasi <- glm(Species ~ . - Endemics,
                 family = quasipoisson,
                 data = gala)
summary(gal.quasi)
plot(gal.quasi$linear.predictors, gal.quasi$residuals)
```
The canonical link for the Poisson is the log link. Identity links are used in rate models, see ELM \S 3.2.

In SM (p. 511) the suggestion is made that one route to an overdispersed Poisson would be to model $Y$ as Poisson with mean $\mu\epsilon$, where $\epsilon$ is a latent random variable with expected value 1. If a gamma distribution with shape $\nu$ and mean 1 is postulated, then the marginal distribution of $Y$ is negative binomial. The negative binomial model can be fit by maximum likelihood using `MASS::glm.nb`; see ELM 3.3.

## Example 3 Gamma 
I will use the `chimp` data set analysed in SM, Example 10.16.

![](../SMtab105.png)

```{r gamma}
data(chimps)
head(chimps)
```

The response is time to learn each of 10 words, in minutes, and there are four chimpanzees. See p.485 of SM. 

Davison's first analysis is a simple `lm(y ~ chimp + word, data = chimps)`, and then he applies Tukey's 1 df test for non-additivity (see Example 8.24) to confirm that the additive model is not a good fit. 
```{r gamma fits}
model1 <- glm(y ~ chimp + word, family=Gamma(link = "log"), data = chimps)
options(digits=6)
anova(model1)
```
This reproduces the first half of Table 10.6. The other side is the same thing with the factors introduced in the opposite order.
```{r gamma fits more}
model2 <- glm(y ~ word + chimp, family = Gamma(link = "log"), data = chimps)
summary(model2)
anova(model2)
```
It makes a bit of difference, but not substantial. 
To test for a `chimp` effect we need an estimate of the dispersion parameter, which is obtained from `summary` as 0.434 (Davison reports .432).  

The `chimp` test statistic is $(6.22/3)/0.434 = 4.77$ which gives a $p$-value of $0.008$ when referred to an $F_{3, 27}$ distribution. The `word` test statistic is $(39.19/9)/0.434$ about $10^{-6}$. 
```{r tests}
test1 <- (6.22/3)/0.434
pf(test1, 3, 27, lower.tail = F)
test2 <- (38.46/9)/0.434
pf(test2, 9, 27, lower.tail = F)
```
Davison suggests two other approaches to this data: an ordinary normal theory model fit to $\log(y)$, and a gamma model with the canonical link $1/\mu$. I'll let you explore those on your own.