---
title: "Lecture 11: Review/Intro to GLMs"
author: "Bob O'Hara"
date: "bob.ohara@ntnu.no"
output: beamer_presentation
---

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

# Before we start...

- By my reckoning we'll have 4 lectures & exercises (plus revision)

- Number of Exercises you'll need to hand in 5 (out if 9 in total)

- This week's exercises don't have to be handed in until next Friday

# Today

A re-cap, to try to put everything in context

Start on Generalised Linear Models

- like the linear models we have done, but with different types of data

# Example: Trying to Explain Bird Brain Size

The bird brain data we have been using was collected to look at the influence of life span on brain size in birds.

Easy, right? But...

```{r getBirdBrains, echo=FALSE}
BirdBrains <- read.csv("../Data/BirdBrains.csv")
BirdBrains$Mode.of.development <- factor(BirdBrains$Mode.of.development)
BirdBrains$M.of.Dev <- c("precocial", "semi-precocial", "semi-altricial", "altricial")[BirdBrains$Mode.of.development]
BirdBrains$M.of.Dev. <- factor(BirdBrains$M.of.Dev, levels = 
                                 c("altricial", "semi-altricial",
                                   "semi-precocial", "precocial"))

Covars <- c("Maximum.lifespan", "Age.at.first.reprodction",
            "Incubation.length", "Clutch.size", "Mean.latitude",
            "logBodyMass", "M.of.Dev.")
# Scale the continuous covariates (messily)
ToScale <- Covars[sapply(Covars, function(wh, df)
  is.numeric(df[,wh]), df=BirdBrains)]
BirdBrains[,ToScale] <- scale(BirdBrains[,ToScale])
```

# Example: Trying to Explain Bird Brain Size

Body mass is correlated with everything

```{r PlotBirdBrainsLiveSpan, echo=FALSE, fig.height=3.5, fig.width=4}
# plot(BirdBrains$Maximum.lifespan, BirdBrains$logBrainMass)
pairs(BirdBrains[,c("logBrainMass", "logBodyMass", "Maximum.lifespan")])
```

May be other correlations too


# Before the modelling...

- Check the data 
    - plot it
    - no obvious problems, except some horrible skew: 
    
```{r HistBirdBrainsMasses, echo=FALSE, fig.height=2, fig.width=5}
par(mfrow=c(1,2), mar=c(4.1,2,1,1))
hist(BirdBrains$Body.mass, xlab="", main="Body Mass (g)")
hist(BirdBrains$Brain.mass, xlab="", main="Brain Mass (g)")
# plot(BirdBrains$logBodyMass, BirdBrains$Maximum.lifespan)
# pairs(BirdBrains[,c("logBrainMass", "logBodyMass", "Maximum.lifespan")])
```
 
# Before the modelling...
    
- Think about what what we want to do

- Ask if brain size is explained by lifespan (maximum lifespan)
- Correct for confounders, e.g. body size

Here: Age at first reproduction, Incubation length, Clutch size, Mean latitude, Mode of Development

Mode of Development

- how developed a chick is when it hatches
- precocial is most developed

# The Strategy

Fit a model with lifespan & other confounders explaining brain size

We want to use the confounders that make a difference, but discard those that don't

Use model selection to chose the confounders

Check the model


# Bird Brains: how I would approach it

For the bird brains, the data were collected to ask if life span was related to brain size. Other variables were included because they might have an effect (e.g. body size)

So, we have a specific hypothesis to test

But we want to find the other covariates that might have an effect

So, we use AIC to find these other covariates


# Model

```{r GetGoodModel, echo=TRUE, message=FALSE, comment=NA}
# Don't use lifespan here, select best model with the rest
library(bestglm)
UseData <- cbind(BirdBrains[,Covars[-1]], 
                 y=BirdBrains$logBrainMass)
AllSubsets <- bestglm(Xy=UseData, IC="AIC")
```

Now add max. lifespan

```{r AddLifespan, echo=TRUE, message=FALSE, comment=NA}

TestLife <- update(AllSubsets$BestModel, 
                   .~.+Maximum.lifespan, data=BirdBrains)
```

# The Result

```{r AnovaLifespan, echo=TRUE, message=FALSE, comment=NA}
anova(AllSubsets$BestModel, TestLife) # could use anova(TestLife)

```

Yep, it still has an effect

# The Model

... but the effect is small compared to the other effects (note: these are standardised)

```{r SummLifespan, echo=TRUE, message=FALSE, comment=NA}
round(summary(TestLife)$coefficients,2)
```

# Model Fit 1

```{r ResidSummLifespan, echo=TRUE, message=FALSE, comment=NA}
plot(fitted(TestLife), resid(TestLife))
```
Looks OK

# Model Fit 2

```{r ResidSummLifespan2, echo=TRUE, message=FALSE, comment=NA, fig.height=4.5}
par(mfrow=c(3,1), mar=c(2,4,1,1))
plot(jitter(BirdBrains$Age.at.first.reprodction), resid(TestLife))
plot(BirdBrains$logBodyMass, resid(TestLife))
plot(BirdBrains$Maximum.lifespan, resid(TestLife))

```


# Model Fit 3

```{r QQNormSummLifespan, echo=TRUE, message=FALSE, comment=NA}
qqnorm(resid(TestLife))
qqline(resid(TestLife))
```
Looks OK, tails a but thick?
 - leptokurtic?

# Model Fit 3

```{r CooksDNormSummLifespan, echo=TRUE, message=FALSE, comment=NA}
plot(fitted(TestLife), cooks.distance(TestLife))



```

Looks OK

# Conclusion

Yes, life span does help to explain brain size

But effect is small: the change from smallest to largest lifespan is `r round(diff(range(BirdBrains$Maximum.lifespan))*coef(TestLife)["Maximum.lifespan"], 2)` which is a bit less than hte change from one extreme to the other of mode of development, and is less than half the effect of changing body mass by 1 standard deviation.

- in part this is because body mass and lifespan are correlated

# Break

![Eric](Eric.jpeg)

# Generalised Linear Models

Not everything is normal

Question: what other data types are there?
(e.g. counts of things)


# A slight diversion

We have 2 colour morphs, with an initial frequency of the white morph of $p$. 

For every individual the red morph produces, the white morph produces $1+s$

How does the frequency of the white morph change?

![](J._B._S._Haldane.jpg)

# Asexual population genetics

The change in frequency:

$$
p_{1} = \frac{(1+s) p_0}{(1+s) p_0 + 1(1-p_0)} = \frac{(1+s) p_0}{ 1 + sp_0 }
$$


```{R Selection}
p <- c(0.01, rep(NA,99))
s <- 0.1
Time <- 1:length(p)
for(t in Time[-1]) { p[t] <- ((1+s)*p[t-1])/(1 + s*p[t-1]) }
plot(Time, p, type="l", xlab="Generations")
```


# Asexual population genetics

Another way of making the same calculation...

$$
\frac{p_{1}}{1 - p_{1}} = \frac{(1 + s)p_{0}}{1 - p_{0}}
$$

So...

$$
\frac{p_{2}}{1 - p_{2}} = (1 + s)\frac{p_{1}}{1 - p_{1}} = (1 + s)(1 + s)\frac{p_{0}}{1 - p_{0}}
$$

or

$$
\log{\frac{p_{t}}{1 - p_{t+r}}} = \mathbf{t} \log(1 + s) + \log{\frac{p_{0}}{1 - p_{0}}}
$$

# We can get a straight Line!

This is linear in $t$

$$
\log{\frac{p_{t}}{1 - p_{t+r}}} = \mathbf{t} \log(1 + s) + \log{\frac{p_{0}}{1 - p_{0}}}
$$

```{R StraightLine}
lgt.p <- log(p/(1-p))
plot(Time, lgt.p, type="l", xlab="Generations")

```

# The (modelling) point

We can transform the model so that it is linear

$$
log{\frac{p_i}{1-p_i}} = \alpha + \beta x_i
$$

# Inference

Suppose we know $p_0$, and then sample $N$ individuals to estimate $p_1$: of these $r$ are white.

We want to maximise the binomial likelihood

$$
l(p | r, N) = \log \left( \begin{array}{c} N \\ r  \end{array} \right) + 
r \log{p} + (N-r) log{(1-p)}
$$

# Tidying up

Re-arrange to collect terms in $r$:
$$
\begin{aligned}
l(p | r, N) &= \log \left( \begin{array}{c} N \\ r  \end{array} \right) + 
r \left( \log{p} - \log{(1-p)} \right)  + N \log{(1-p)} \\
            &= \log \left( \begin{array}{c} N \\ r  \end{array} \right) + 
            r \log \left(\frac{p}{1-p} \right)  + N \log{(1-p)}
\end{aligned}
$$

Back to $\log{p/(1-p)}$!

# Generalising Linear Models

It turn out that it's really convenient to model $\log{p/(1-p)}$:

$$
\log \frac{p_i}{1-p_i} = \sum_j X_{ij} \beta_j
$$

and we can work with the right hand side like we do with multiple regression.

# Components of a Generalized Linear Model

Random Part

 - the data (we assume a distribution, e.g. binomial)
 
 Systematic Part
 
 - the model for the mean for each data point: $\sum_j X_{ij} \beta_j$

The link Function

- transforms the model (which is linear) onto the scale of the data
- e.g. the _logit_ link function, $\log{p/(1-p)}$

# Binomial

We look at selection over time. We take samples of size $N_t$ at time $t$, and look at the colour, and record the number that are white, $r_t$. Then we have

Random part:
$$
r_t ~ Binom(N_t, p_t)
$$

Link function:
$$
\log \frac{p_t}{1-p_t} = \theta_t
$$
Systematic Part: 
$$
\theta_t = \alpha + \beta t
$$


(*) this is not the best model, but is OK for now

# Simulated data

```{r echo=TRUE, fig.height=3.5, fig.width=4}
N <- 5
r <- rbinom(length(p), N, p)
plot(Time, r)
lines(Time, p*N, col=2, lwd=1.6)
```

# Fitting the Model

```{r FitGLM, echo=TRUE, comment=NA}
glmod <- glm(cbind(r, N-r) ~ Time, 
             family=binomial("logit"))
round(summary(glmod)$coefficients, 3)

```

# The technical bit (i.e. not on the exam)

Nelder & Wedderburn showed that a lot of common models have likelihoods (for one data point) that look similar:

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

- $y$ is the data
- $\theta$ is the expected value (the mean)
- $\phi$ is the dispersion (the variance)
- $a()$, $b()$ & $c()$ are functions

Often $\phi$ is known, so this is also written as

$$
l(\theta | y) = a(y) t(\theta) + c(\theta) + d(y)
$$

# For example 

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

For the binomial we'll focus on the proportion, so $y = r/N$

$$
l(p | r, N) = \log \left( \begin{array} {c} N \\ r  \end{array} \right) + 
N \left( \frac{r}{N} \log \left(\frac{p}{1-p} \right)  + \log{(1-p)} \right)
$$

$$a(\phi) = 1/N$$
$$\theta = \log \left(\frac{p}{1-p} \right)$$

and $b(\theta) = \log{(1-p)}$ & $c = \log \left( \begin{array}{c} N \\ r  \end{array} \right)$

# Likelihood for everything

$$
l(\theta_i | y_i) = \sum {\frac{y_i \theta_i -b(\theta_i)}{a(\phi)} + c(y_i, \phi)}
$$

We can make $\theta_i$ a linear model (like we do for regression): 

$$\theta_i = \sum_j {X_{ij} \beta_j}$$

# The Workhorse of Modern Statistics

It turns out that this has some nice properties

- the likelihood is easy to maximise 
- the standard errors don't need any special tricks to calculate
- the model often behaves really nicely

Lots of models come under this framework, and many extensions come from it

# What's coming

There are several distributions we can use in GLMs. We will look at these:

- Poisson (for counts)
- Binomial (for counts of proportions)

But there are others (e.g. Gamma)

Distribution, link functions are new: the systematic part is the same

- but interpretation of the parameters needs care
