---
title: "Lecture 10: Selecting Models"
author: "Bob O'Hara"
date: "bob.ohara@ntnu.no"
output: beamer_presentation
---

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

# Before we start...

Number of Exercises: ```ceiling(N/2)```

This week's exercises don't have to be handed in until after Easter (I will out them up tomorrow)

# Models

```{r GetCake, results='asis', echo=FALSE, message=FALSE}
library(lme4)
data("cake")
# Sugar <- cake[cake$recipe!="B" & cake$temp=="225",] # remove recipe B & use temperature of 175C
# Sugar$treatment <- c("Control", "Treatment")[1+grepl("A", Sugar$recipe)]
```

So far we have been fitting models, looking at how well they fit, and (hopefully) working out what they mean. But we also often want to compare models

- pick from a large number of models (QTLs)
- ask if one particular hypothesis explains the data

# Why select models?

We could fit a model with every covariate in it

- plus all interactions

But these get big & difficult to interpret. So we only want to important variables

Smaller models generally have smaller standard errors

- more precise



# Two types of problem: two solutions

Testing a specific hypothesis

- does the temperature at which a cake is baked affect how much you can break it?
- confirmatory

Finding a good model

- of these $10^6$ genetic markers, which ones explain a trait?
- exploratory

# Hypothesis Testing

e.g. does the temperature at which a cake is baked affect how much you can break it?

The basic approach: 

- fit a model with temperature
- fit a model without temperature

Compare them & see which fits better

# Hypothesis Testing

Hypothesis Testing is asymmetrical. We ask "does is model without the effect sufficient to explain the data?"


# How to Do Statistical Hypothesis Testing

- get a *null hypothesis* (i.e. without the effect)
- get an *alternative hypothesis* (i.e. with the effect)
- Chose a *test statistic* (e.g. the likelihood)
- calculate the distribution of the test statistic if the null hypothesis was true
- ask if the observed value of the statistic falls within the null distribution
- if it does not, declare the null hypothesis wrong

# Cake Testing

Null hypothesis: no factor explains the cake angle

Alternative hypothesis: angle can be explained by temperature

Test statistic: the likelihood


# The Test statistic

We could use anything, but the likelihood makes theoretical sense

- the probability of the data given the model and parameter estimates

Need to fit the model to get the parameter estimates:

```{r FitCake, echo=TRUE, comment=NA}
mod.null <- lm(angle~1, data=cake)
round(coef(mod.null), 2)
mod.alt <- lm(angle~temperature, data=cake)
round(coef(mod.alt), 2)

```


# The Test statistic


Question is whether the alternative model explains the data better than the null model, i.e. if the likelihood is higher

- us the ratio of likelihoods (or difference in log-likelihoods)

The likelihood is $P(Y | \hat{\theta})$ - the probability of the data given the maximum likelihood estimates of the parameters.

With maximum likelihood, the data are random. We ask

- if we had an inifinite number of samples for the data, how often would the alternative model be better than the null model 
- (given the MLEs)

# The Distribution of the Null Likelihood

simulate() simulates data from the model in mod.null

```{r PlotNullCakeCake, echo=TRUE, fig.height=3.5, fig.width=5}
CalcLhood <- function(y, x) {
  mod.null <- lm(y ~ 1)
  mod.alt <- lm(y ~ x)
  logLik(mod.alt) - logLik(mod.null)
}

# simulate the data from the (null) model 
SimNull <- simulate(mod.null, nsim=1e3)
Lhood.null <- apply(SimNull, 2, CalcLhood, 
                    x=cake$temperature)
```

# Compare to Likelihoods

```{r PlotNullCakeCake2, echo=TRUE, fig.height=3, fig.width=5}
hist(Lhood.null, col="light blue", border=NA, 
     main="", xlim=c(0, 17))
```

If temperature has an effect, the likelihood should be higher

# Compare to Likelihoods

```{r PlotNullCakeCake3, echo=TRUE, fig.height=3, fig.width=5}
hist(Lhood.null, col="light blue", border=NA, main="", xlim=c(0, 17))
abline(v=CalcLhood(cake$angle, cake$temperature), col=2, lwd=3)

```

Here it is larger than all of the likelihoods

# Without simulation

We know from statistical theory that the distribution of the ratio of log-likelihoods should follow an F distribution, so we don't need to simulate it

The F distribution has 2 parameters, known a "degrees of freedom".

- numerator degrees of freedom: how many extra parameters are in the alternative model
- denominator degrees of freedom: how many parameters are used to estimate $\hat{\sigma}^2$
    - taken from the alternative model

# With R

We can get R to make the comparison:

```{r Ftest, comment=NA}
mod.null <- lm(angle ~ 1, dat=cake)
mod.alt <- lm(angle ~ temperature, dat=cake)

anova(mod.null, mod.alt)

```

The test statistic is the F-ratio (and the p-value)

# Degrees of Freedom

We have $N$ data points. Each is a "degree of freedom" that we can use in the estimation
Each df can be spent to estimate one parameter
The rest are used to estimate the residual variance

e.g. 

$$
y_i = \alpha + \beta x_i + \varepsilon_i
$$

2 parameters ($\alpha$ and $\beta$), so $N-2$ can be used to estimate $\sigma^2$

$N-2$ is the residual degrees of freedom

# Degrees of Freedom

If we compare 2 models, the difference in the residual degrees of freedom is the number of extra parameters in the alternative model

- this is the degrees of freedom.

(the same as used in a $\chi^2$ test)

# With R

We can get R to make the comparison:

```{r Ftest2, comment=NA}
mod.null <- lm(angle ~ 1, dat=cake)
mod.alt <- lm(angle ~ temperature, dat=cake)

(an <- anova(mod.null, mod.alt))

```

- df is the degrees of freedom for the model


# RSS

NOT the Royal Statistical Society

The Residual Sum of Squares

- the log-likelihood for a normal distributiuon is 

$$
l(\mathbf{x}|\mu, \sigma^2) = - \frac{n}{2} \log{\sigma^2} - \sum_{i=1}^n {\frac{(x_i - \mu)^2}{2\sigma^2} }
$$

And the main bit is $\sum_{i=1}^n {\frac{(x_i - \mu_i)^2}{2\sigma^2}}$

which is just a sum of squares, and a variance term

# Next

```{r Ftest3, comment=NA}
anova(mod.null, mod.alt)

```

So the RSS are (almost) likelihoods

- need to divide by an estimate of $\sigma^2$
- this is the RSS for the largest model divided by its residual degrees of freedom

# Why an F distribution?

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

when we divide by the residual stanadrd deviation we add more uncertainty

- but this also follows a $\chi^2$ distribution

The ratio of 2 independent $\chi^2$ distributions (divided by their degrees of freedom) is an F distibution

- later we will see $\chi^2$ distributions in similar tests


# The p-value

```{r PlotF}
curve(df(x, df1=an$Df[2], df2=an$Res.Df[2]), from=0, to=8,
      ylab="Density")
abline(v=an$F, col=2)
```

The p-value is $Pr(F>F_{obs})$. Here it is `r round(an[6,2],9)`

- so very unlikely

# ANOVA made easier

We have just used anova() to compare 2 models, but it has traditionally been used to compare several:

```{r ANOVA, comment=NA}
mod.all <- lm(angle~temperature+recipe+replicate, data=cake)
round(anova(mod.all), 2)
  
```


# ANOVA made easier

Each row is a test


```{r ANOVA.main, comment=NA, results="asis"}
mod.main <- lm(angle~temperature+recipe+replicate, data=cake)
print(xtable::xtable(anova(mod.main)), digits=2, comment=FALSE)
  
```

It compares a model with the terms above to one including that term e.g. the replicate line compares

```temperature + recipe```

to

```temperature + recipe + replicate```


# Why is ANOVA called ANOVA

ANOVA = Analysis of Variance

The Mean Sq is the mean square, i.e. Sum Sq/Df

- when the data cooperate, it is an estimate of the variance explained by that effect
- so we can sometimes use the Mean Square to eyeball how important a variable is

# Exploring for Good Models

Sometiumes we don't have strong hypotheses. Instead we might be exploring which variables might have an effect

- our aim is to get a good model overall
- e.g. for prediction

# What does a good model look like?

- Simple
- Fits the data well
- Understandable

We can measure simplicity and fit.

- Fit: likelihood
- Simplicity: number of parameters

# The problem

A model with always fit better if you add a parameter

Key issue: is adding the extra parameter worth it?

ANOVA answers this by asking if the improvement from the extra parameter can be explained as noise

# Penalisation

Another way of looking at the problem: we measure model adqequacy 

- is the model good enough for what we want?

We penalise complicated models

- measure complexity by number of parameters

Find the 'best' model as one with optimum between fit & complexity

# How to Penalise

There are several ways to penalise. Here I will mention two, which chose different criteria

- *AIC*: Akaike's Information Criterion
- *BIC*: Bayesian Information Criterion

AIC tries to find the model that best predicts thedata 

BIC trues to find the model most likely to be true

# AIC

Finds the model that would best predict replicate data

AIC = -2\*Likelihood + 2\*Number of Parameters

# BIC

Finds the model which is most likely to be "true"

BIC = -2\*Likelihood + log(N)\*Number of Parameters

- log(n) = log(sample size)
- penalises more than AIC

# Using AIC/BIC

Full Subset Selection

- calculate AIC/BIC for every model
- pick the best

Usually, if the values are within ~2 of each other, the models are pretty similar.

# Example: Trying to Explain Bird Brain Size

(we will get to the actual purpose of this study in a moment)

```{r getBirdBrains, echo=TRUE}
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])
```

Look at:

- logBodyMass
- Maximum lifespan
- Age at first reprodction
- Incubation length
- Clutch size
- Mode of development
- Mean latitude

# Fit the model

```{r FitFullModel, echo=TRUE, message=FALSE, comment=NA}
library(bestglm) # might need install.packages("bestglm")
UseData <- cbind(BirdBrains[,Covars], y=BirdBrains$logBrainMass)
AllSubsets <- bestglm(Xy=UseData, IC="AIC")
AllSubsets

```

# The Best model

Max. lifespan, body mass, mode of development (how developed a chick is when it hatches: precocial is most developed)

Model has $R^2$ of `r round(100*summary(AllSubsets$BestModel)$r.square)`

```{r BestModel, comment=NA}
round(summary(AllSubsets$BestModel)$coefficients, 2)
```

# Why we should't just use anova() Pt 1

```{r bestANOVA, echo=TRUE, message=FALSE, comment=NA}
anova(lm(logBrainMass~Mean.latitude+logBodyMass+Maximum.lifespan+M.of.Dev,
         data=BirdBrains))

```

Mean latitude has an effect!

# Why we should't just use anova() Pt 2

```{r bestANOVA2, echo=TRUE, message=FALSE, comment=NA}
anova(lm(logBrainMass~Maximum.lifespan+M.of.Dev+logBodyMass+Mean.latitude,
         data=BirdBrains))

```

Mean latitude has no effect!

# Why we should't just use anova()

Because the tests are sequential, the order usually matters

(unless the data are from a well-designed experiment)

# When to Use these Methods

ANOVA: testing specific hypotheses

AIC/Full subsets: finding the best 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 spcific hypothesis to test

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

So, we use AIC to find these other covariates

# Bird Brains: how I would approach it

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

# Now add max. lifespan
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)

```

# Next Week

I'll be jealous of you in Borneo

After Easter, I will wrap up where we are & how everything fits together

The exercises don't have to be handed in until after Easter.

