---
title: "Lecture 12: log-linear models/Poisson regression"
author: "Bob O'Hara"
date: "bob.ohara@ntnu.no"
output: beamer_presentation
---

```{r setup, include=FALSE}
# April 11th
knitr::opts_chunk$set(echo = FALSE)
```

# Before we start...


- This week's exercises wioll be available soon (sorry)
- I will try to get a syllabus up this week.
- Last week's exercises don't have to be handed in until Friday

# Data Generating Models


Modern statistics deals much more with mechanisms

One major part: how the data were collected

# Poisson Processes

Assume events happen at a constant rate, $\lambda$. If we observe for a time $t$ then the expected number of events is $\mu = \lambda t$. The actual number varies around this, and follows a Poisson distribution:

$$
Pr(N = r) = \frac{e^{-\mu} \mu^r}{r!}
$$

# A Poisson Distribution

```{r PlotPoisson, comment=NA, results='hide'}
par(mfrow=c(2,3), mar=c(4.1,3,2,1))
sapply(c(0.1, 0.5, 1, 2, 10, 100), function(mu, N) { 
  dat <- table(rpois(N, mu))
  plot(dat, type="h", lwd=6, lend=3, col="grey30", main = paste0("Mean: ", mu), cex.main=2)
#  hist(rpois(N, mu), col="grey", main = paste0("Mean: ", mu))
}, N=1e3)

```


# Propoerties of the Poisson

if $r_1 \sim Poisson(mu_1)$ and $r_2 \sim Poisson(mu_2)$

- $E(r_1) = \mu_1$
- $Var(r_1) = \mu_1$
- $r_1 + r_2 \sim Poisson(\mu_1 + \mu_2)$
- $r_1 | r_1 + r_2 \sim Binomial(r_1 + r_2, \mu_1/(\mu_1 + \mu_2))$
- If $s \sim Binomial(N, p)$ with large $N$ and small $p$ then $s \approx Poisson(np)$

Poisson and binomial distributions are closely linked

# Inference

Suppose we observe $n$ counts from a Poisson with unknown mean $\mu$, what is the maximum lilelihood estimate?

$$
Pr(N=r | \mu) = \frac{e^{-\mu} \mu^r}{r!}
$$
so 
$$
l(\mu | r) = -\mu + r log(\mu) - log(r!)
$$

DIffernetiate & set to 0:

$$
0 = -1 + r \frac{1}{\mu}
$$

so $\hat{\mu} = r$

# Is this a GLM?

Remember that a GLM has a likelihood with the form 

$$
l(\theta | y) = \frac{y \theta -b(\theta)}{a(\phi)} + c(y, \phi)
$$

For the Poisson we have $l(\mu | r) = -\mu + r log(\mu) - log(r!)$, so 

- $\theta = log(\mu)$, 
- $a(\phi) = 1$,
- and $b(\theta) = -e^\theta$, $c(y, \phi) = - log(r!)$

# What This Means I

- $\theta = log(\mu)$

We have, naturally, a log link.

- this is the _canonical link_ 

Makes sense: if we are counting, the process is multiplicative (double the effort, double the counts)

This is additive on the log scale.

# What This Means II

- $a(\phi) = 1$,

The dispersion is fixed, same as saying $Var(r) = \mu$

The amount of variation is determined by the mean

(we will see how to relax this later)

# Interpretation

The log link means that the model is multiplicative

$$
\begin{aligned}
\log(\mu) = \alpha + \beta x \\
\mu = e^{\alpha+\beta x} = e^\alpha e^{\beta x}
\end{aligned}
$$

So the effect is multiplicative. For example, let $x$ be 0 or 1, and $\beta = 0.01$. The means are

$$
\begin{aligned}
\mu_0 &= e^{\alpha+0.01 0} = e^\alpha \\
\mu_1 &= e^{\alpha+0.01 1} = e^{\alpha + 0.01}
\end{aligned}
$$

So the ratio $\mu_1/\mu_0$ is $e^{\alpha + 0.01}/e^{\alpha} = e^{0.01} \approx 1.01$ 

If a coefficient is small, it is (approximately) the percent increase

# Symmetry

The coefficients are symmetrical

e.g. if $\beta = -0.01$ then

$$
\mu_1 = e^{\alpha - 0.01 1} = e^{\alpha} e^{-0.01} = e^{\alpha}/e^{0.01}
$$

- reduces the mean by $e^{0.01}$ rather than increasing it by $e^{0.01}$

# Hypothesis Testing and Deviance

We can use AIC/BIC just like before. But ANOVA is a bit different


- AIC = Deviance + 2*Number of parameters
- lowest is best

# Deviance

From before Easter:

-2 times difference in likelihood between 2 (nested) models follows a $\chi^2$ distribution

We call $-2l(\theta | Y)$ the _deviance_

So we can test whether a term shoud be in the model 

# Fitting a GLM in R

This is easy:

```{r FitGLM, comment=NA, echo=TRUE}
X <- 1:100
SimR <- rpois(100, lambda = exp(1.5 + 0.0001*X))
mod <- glm(SimR ~ 1, family=poisson)
mod1 <- glm(SimR ~ X, family=poisson)

# More formally
mod <- glm(SimR ~ 1, family=poisson("log"))

```

- we use `glm` not `lm` 
- `family=poisson` says to use the Poisson distribution
- `family=poisson("log")` says to use the Poisson distribution with a log link
- if we do not specify a link function, R will use the canonical link
    - i.e. the log link for the Poisson

# Analysis of Deviance

```{r AnovaGLM, comment=NA, echo=TRUE}
anova(mod, mod1, test="LRT")
```

# Looking at a GLM in R

```{r SummGLM, comment=NA, echo=TRUE}
summary( mod1)
```


# An Example: The Hastings Rarities

British birders were worried that a lot of observations from around Hastings betwen 1890 and 1930 were frauds

John Nelder (co-inventor of GLMs) took a look at the data, and compared Hastings with two areas next to Hastings

https://en.wikipedia.org/wiki/Hastings_Rarities

# Hastings Rarities Data

- Year (1895 to 1954)
- Area (Hastings, Sussex, Kent)
     - Hasting is a town in Sussex: Kent is next door
- Era (A, B: A is before 1925)
    - A is when the frauds were thought to occur
- Class (National rarity of species)
- Count: number of records

We will only look at the rarest species (Class I)

The problem: were there more rarities recorded around Hastings before 1925 (in Era A)?

# Hastings Data

```{r HastingsData, echo=FALSE}
HastingsFile <- "../Data/HastingsData.csv"
if(file.exists(HastingsFile)) { 
  HastingsYears <- read.csv(HastingsFile, stringsAsFactors = FALSE) 
} else {
  YearA <- 1895:1924
  XA.I <- c(2,4,2,0,1,4,7,7,10,7,19,9,13,8,15,6,14,14,8,27,22,15,6,4,3,5,0,3,6,2)
  XA.II <- c(3,2,0,2,1,1,7,3,4,8,7,4,3,3,3,2,3,2,7,9,7,5,1,5,2,3,2,4,2,3)
  XA.III <- c(3,2,5,4,3,6,3,9,11,9,16,9,3,6,4,7,10,5,5,5,5,4,4,4,4,4,1,5,2,7)
  
  YSA.I <- c(0,0,0,1,0,1,0,1,1,0,1,0,0,1,3,1,1,0,2,0,0,1,0,0,0,0,0,1,0,0)
  YSA.II <- c(0,0,1,1,0,0,1,0,0,3,3,1,2,0,0,0,0,0,0,0,0,0,0,1,0,0,1,2,0,0)
  YSA.III <- c(0,1,1,0,4,1,3,1,3,1,0,3,3,4,5,2,1,3,0,1,0,0,1,1,1,1,0,2,0,2)

  YKA.I <- c(0,1,1,0,1,0,0,2,0,1,1,1,2,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0)
  YKA.II <- c(0,0,0,1,2,0,0,0,0,1,0,0,3,1,0,0,2,0,0,0,0,0,0,0,0,0,1,0,0,0)
  YKA.III <- c(1,3,1,1,0,0,0,0,2,0,0,4,1,1,1,0,0,01,1,0,1,1,0,0,1,2,0,0,0,0)

  YearB <- 1925:1954
  XB.I <- c(0,1,0,0,0,0,2,1,2,3,3,6,3,1,0,0,0,1,0,1,0,3,3,3,2,4,5,4,3,3)
  XB.II <- c(2,1,1,1,1,0,2,2,0,1,1,0,2,2,0,1,2,0,2,0,0,0,2,3,2,5,7,4,2,5)
  XB.III <- c(5,1,0,3,1,0,1,0,2,3,3,6,5,4,2,1,0,2,0,1,1,2,6,1,6,11,7,20,7,12)

  YSB.I <- c(0,1,0,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,2,0,4,2,1,2,2)
  YSB.II <- c(0,0,0,0,0,0,1,0,0,2,0,0,0,1,0,0,0,1,0,0,0,1,0,3,1,2,0,1,0,0)
  YSB.III <- c(0,1,1,2,3,0,2,2,2,0,0,1,4,0,0,0,0,0,0,1,0,0,1,2,3,2,1,2,1,1)

  YKB.I <- c(0,0,0,0,0,0,0,0,1,1,0,1,0,1,1,0,0,2,2,0,1,0,0,0,1,4,2,1,4,4)
  YKB.II <- c(0,1,1,0,0,0,0,1,0,2,1,1,1,0,0,1,0,0,1,0,0,2,0,3,0,2,1,1,4,5)
  YKB.III <- c(0,0,0,1,1,0,0,2,0,0,2,2,1,0,1,0,2,2,2,3,1,1,1,0,1,5,6,3,10,6)

  HastingsYears <- data.frame(
    Year = c(rep(YearA, 9), rep(YearB, 9)),
    Cat = rep(c("XA.I", "XA.II", "XA.III", "YSA.I", "YSA.II", "YSA.III", 
                "YKA.I", "YKA.II", "YKA.III",
                "XB.I", "XB.II", "XB.III", "YSB.I", "YSB.II", "YSB.III", 
                "YKB.I", "YKB.II", "YKB.III"), each=length(XA.I))
    )
  HastingsYears$AreaCode <- gsub("[AB]\\..*", "", HastingsYears$Cat)
  HastingsYears$Era <- c("A", "B")[1+grepl("B", HastingsYears$Cat)]
  HastingsYears$Class <- gsub(".*\\.", "", HastingsYears$Cat)

  HastingsYears$Area <- sapply(HastingsYears$AreaCode, function(x) switch(x,
                                                                          X = "Hastings",
                                                                          YS = "Sussex",
                                                                          YK = "Kent"))
  HastingsYears$AreaCode <- NULL
  HastingsYears$Cat <- NULL
  HastingsYears$Count <- c(XA.I, XA.II, XA.III, YSA.I, YSA.II, YSA.III, YKA.I, YKA.II, YKA.III, 
                           XB.I, XB.II, XB.III, YSB.I, YSB.II, YSB.III, YKB.I, YKB.II, YKB.III)
  write.csv(HastingsYears, file=HastingsFile)
}

```

Concerned about Hastings before abnout 1925....

```{r PlotHastingsData, echo=FALSE}
HastingsYearsI <- HastingsYears[HastingsYears$Class=="I",]

HastingsYearsI$Era <- factor(HastingsYearsI$Era)
HastingsYearsI$Era <- relevel(HastingsYearsI$Era, ref = "B")

HastingsYearsI$Area <- factor(HastingsYearsI$Area)
HastingsYearsI$Area <- relevel(HastingsYearsI$Area, ref = "Sussex")


# Make sure baselines are correct
HastingsYearsI$Colour <- c("red2", "blue", "red4")[as.numeric(HastingsYearsI$Area)]

plot(HastingsYearsI$Year + 0.2*(as.numeric(HastingsYearsI$Area)-2), 
     HastingsYearsI$Count, type="n", xlab="Year", ylab="Frequency")
rect(1850, -5, 1925, 59, col="grey90", border=NA)
points(HastingsYearsI$Year + 0.2*(as.numeric(HastingsYearsI$Area)-2), 
     HastingsYearsI$Count, col = HastingsYearsI$Colour, pch=18, cex=2)
legend(1935, 20, levels(HastingsYearsI$Area), col=c("red2", "blue", "red4"), pch=18, cex=2)
box(bty="o")

```

# Fitting & Testing the Model

```{r hastingsModel, comment=NA, echo=TRUE}
Hast.mod <- glm(Count ~ Area*Era, family="poisson", 
                data=HastingsYearsI)
anova.hast <- anova(Hast.mod, test = "Chisq")
signif(data.frame(anova.hast), 3)

```

- The test statistic for the interaction is `r round(anova(Hast.mod)["Area:Era","Deviance"], 2)` with `r anova(Hast.mod)["Area:Era","Df"]` DF
- If we test this against a $\chi^2$ distribution, we get p = `r signif(data.frame(anova.hast)["Area:Era", "Pr..Chi."], 2)`.

# Parameter Estimates

```{r hastingsPars, comment=NA, echo=TRUE}
round(summary(Hast.mod)$coefficients, 2)
```

We can see that there were about exp(`r round(coef(Hast.mod)["AreaHastings"], 2)`) =  `r round(exp(coef(Hast.mod)["AreaHastings"]),2)` times as many rare species around Hastings in Era B

But exp(`r round(coef(Hast.mod)["AreaHastings:EraA"], 2)`) =  `r round(exp(coef(Hast.mod)["AreaHastings:EraA"]),2)` times more than _that_ before 1925


# 95% Confidence intervals

All parameters (on log scale):
```{r Conint, echo=TRUE, message=FALSE, comment=NA}
round(CI <- confint(Hast.mod), 2)
```

For Hasting:AreaA on couint scale

```{r Conint2, echo=TRUE, message=FALSE, comment=NA}
round(exp(CI["AreaHastings:EraA",]),2)
```


# The Full Summary 1

```{r HastingsSumm1, comment=NA, echo=TRUE}
summ.Hast <- paste(capture.output(print(summary(Hast.mod), digits=2)), "\n")
cat(summ.Hast[2:7])

```

The model (useful reminder) and summary of residuals (ignore!)

# The Full Summary 2

```{r HastingsSumm2, comment=NA, echo=FALSE}
cat(summ.Hast[9:18])

```

Parameter estimates: very useful.

- also use coef() and confint()

# The Full Summary 3

```{r HastingsSumm3, comment=NA, echo=FALSE}
cat(summ.Hast[19:26])
```

- Dispersion was mentioned last week; more on this and use of deviance later...
- AIC not useful without more models
- Ignore Fisher scorings: relate to efficiency in model fitting


# Residuals

Our raw residuals are $y_i - E(y_i)$, i.e. Observed - Expected

We can plot these

```{r PlotResids, echo=TRUE, fig.height=4}
x <- 1:100; y <- rpois(length(x), 0.5*x)
mod <- glm(y~log(x), family="poisson")
Resid <- y - fitted(mod)
par(mfrow=c(1,2))
plot(x, y); plot(fitted(mod), Resid)
```

Remember, $Var(x) = E(x)$



# Better Residuals

There are several solutions to this problem (none perfect)

- Pearson residuals: $(x - \mu_x)/\sigma_x$
- Deviance Residuals: $sgn(y_i - E(y_i))\sqrt{D_i}$
    - deviance for one datum is $D_i = -2*l(y_i|\theta_i)$
    - $sgn(x)$ is 1 if $x>0$ and -1 if $x<0$

Deviance residuals are the default: they control for some of the variation in shape, but aren't perfect


# Better Residuals

```{r PlotBetterResids, echo=TRUE, fig.height=4}
par(mfrow=c(1,3))
plot(x, y); plot(fitted(mod), resid(mod, type="deviance")); 
plot(fitted(mod), resid(mod, type="pearson")); 
```


# Residuals for Discrete Data

We get lines, from when $y=0,1,2....$

```{r DiscRes}
y <- rpois(length(x), 0.05*x)
mod.d <- glm(y ~ x, family="poisson")
par(mfrow=c(1,2))
plot(x, y)
plot(fitted(mod.d), resid(mod.d))
```


This is normal (if annoying)

# Model Checking

Deviance Residuals can still be informative

```{r PlotSimResids, echo=TRUE, fig.height=4}
x <- 1:100; 
y <- rpois(length(x), exp(2-0.015*x + 0.0004*x^2))
mod <- glm(y~x, family="poisson")
par(mfrow=c(1,2))
plot(x, y); plot(x, resid(mod))
```

Look, curvature!

# Hastings Residuals

```{r PlotHastingsResids}
par(mfrow=c(1,2))
plot(fitted(Hast.mod), resid(Hast.mod))
plot(HastingsYearsI$Year, resid(Hast.mod), col=HastingsYearsI$Colour)
legend(1925, 4.6, levels(HastingsYearsI$Area), col=c("red2", "red4", "blue"), pch=1, cex=1.2)

```

Rise in reports after 1950?


# Overdispersion

We assume that the mean controls the variance

But this is not always true: there might be extra variation

- e.g. for the Hastings data there might be more variation between years

We can check this!

# Overdispersion

If the mean controls the variance, it then controls the amount of residual deviance

It turns out that the residual deviance should (asymptotically) follow a $\chi^2$ distribution

- If there are lots of DF then residual deviance $\approx$ Deviance
- so we can test this!

# Testing Overdispersion

Remember the summary?

```{r CheckHastingsOD, echo=TRUE, comment=NA}
cat(summ.Hast[20:24])
```

We can use this to test for overdispersion

# Testing Overdispersion

The residual deviance is `r round(deviance(Hast.mod), 2)`, with `r df.residual(Hast.mod)` degrees of freedom. The p-value is 

```{r testDisp, echo=TRUE}
pchisq(deviance(Hast.mod), 
       df=df.residual(Hast.mod), 
       lower.tail = FALSE)
```

So it is unlikely that the data come from a Poisson distribution

- could be that there is another variable that should be in the model
- or there is just a lot more variation

# Estimating Overdispersion

The ratio of deviance degrees of freedom is `r round(deviance(Hast.mod), 2)`/`r df.residual(Hast.mod)` = `r round(deviance(Hast.mod)/df.residual(Hast.mod), 2)`

This is more useful: it acts like a residual variance.

If there is no overdispersion, this should be 1.

# Dealing With Overdispersion

There are a few ways to deal with overdispersion

- Correct in the likelihood
- Use a mixed model (later?)
- Use a different distribution

# Correct in the likelihood

The likelihood is 

$$
l(\theta | y) = \frac{y \theta -b(\theta)}{a(\phi)} + c(y, \phi)
$$

So we can estimate $\phi$, the dispersion. We can use the deviance ratio.

Deviance/Degrees of Freedom

```{r EstDisp, echo=TRUE}
Dispersion <- deviance(Hast.mod)/df.residual(Hast.mod)
```

We can plug that into the summary:

```{r SummDisp, echo=TRUE}
summary(Hast.mod, dispersion = Dispersion)
```

# Effect of Overdispersion

Effect is to increase standard errors by sqrt(Dispersion):

```{r SummDispStdErr, echo=TRUE}
round(summary(Hast.mod)$coefficients[1:3,"Std. Error"],2)
round(summary(Hast.mod, dispersion = 
                Dispersion)$coefficients[1:3,"Std. Error"],2)
```

# Effect of Overdispersion

Similar effect on anova():

```{r AnovaDispStdErr, echo=TRUE}

cat(paste(capture.output(anova(Hast.mod, test="LRT")), "\n")[10:13])
cat(paste(capture.output(anova(Hast.mod, dispersion = 100, test = "LRT")), "\n")[10:13])

```

# Use a different distribution

The Negative Binomial distribution assumes that there is over-dispersion

```{r GLMNB, echo=TRUE}
Hast.NB <- MASS::glm.nb(Count ~ Area*Era, data=HastingsYearsI)

round(summary(Hast.NB)$coefficients, 2)
```

# Use a different distribution: long version

Our model is $\log(\mu_i) = \sum_j X_{ij} \beta_j$. But we could add a random term, so it becomes $\log(\mu_i) = \sum_j X_{ij} \beta_j + \varepsilon_i$

If we use $\varepsilon_i \sim N(0, \sigma^2)$ this is like a regression

- need a Generalised Linear Mixed Model to estimate it

We could also use $e^{\varepsilon_i} \sim \chi^2_\nu$. This is the same as assuming a negative binomial distribution.

# Summary

- GLMs are like LMs
- anova() is (almost) the same
- summary() is the same
    - but parameter interpretation is important
- Because the Poisson assumes the variance, we might have to deal with that
- We should check for overdispersion & correct if necessary