---
title: "Lecture 6: More Regression"
author: "Bob O'Hara"
date: "bob.ohara@ntnu.no"
output: beamer_presentation
---

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

## Regression & R

Read in the data (and sub-sample)
```{R ReadData, echo=TRUE, fig.height=4}
BirdEggs <- read.csv(file="../Data/BirdEggs.csv", 
                     stringsAsFactors = FALSE)
# sub-sample 30 observations, just to make things clearer later
BirdEggs <- BirdEggs[sample.int(nrow(BirdEggs), size=20),]
```


## Plot the data

```{R PlotData, echo=TRUE, fig.height=4}
par(mar=c(4.1,4.1,1,1))
plot(BirdEggs$logFemaleMass, BirdEggs$logEggMass, 
     xlab="Log body mass", ylab="Log egg mass", 
     main="", lwd=2, col="black", pch=4)

```

## The Model

The model is a straight line: we are regressing y against x 

 _____

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

where $x_i$ is the log body mass, and $y_i$ is the log egg mass, and 

$$
\varepsilon_i \sim N(0, \sigma^2)
$$

## Fitting the Model in R

The code to fit the model is simple:

```{r FitModel, echo=TRUE}
mod <- lm(logEggMass ~ logFemaleMass, data=BirdEggs)
```

- logEggMass ~ logFemaleMass is the formula that describes the model: Y~X means we regress Y against X 
- data=BirdEggs just gives the data frame where the data are

## The Parameter Estimates

```{r GetCoefs, echo=TRUE}
round(coef(mod), 2)
round(sigma(mod), 2)

```

So the model is 

$$
\begin{aligned}
y_i &= `r round(coef(mod)[1], 2)` + `r round(coef(mod)[2], 2)` x_i + \varepsilon_i \\
\varepsilon_i &\sim N(0, `r round(sigma(mod), 2)`^2)
\end{aligned}
$$

## Plotting the Model

```{R PlotModel, echo=TRUE, fig.height=4}
par(mar=c(4.1,4.1,1,1))
plot(BirdEggs$logFemaleMass, BirdEggs$logEggMass, 
     xlab="Log body mass", ylab="Log egg mass", 
     main="", lwd=2, col="black", pch=4)
abline(mod, col="blue", lwd=3)
```

## Uncertainty in the parameters

```{r GetConfints, echo=TRUE}
round(confint(mod), 2)
```

## Looking at uncertainty in the intercept

If we fix the slope, then we can look at variation in the intercept. These are the upper & lower 95% confidence limits

```{r PlotIntercepts, echo=FALSE, results='hide'}
ConfIntervals <- confint(mod)

ConfPoly.x <- c(-100,100,100,-100)
ConfPoly.y <- rep(ConfIntervals[1,], each=2) + coef(mod)[2]*ConfPoly.x
par(mar=c(4.1,4.1,1,1))
plot(BirdEggs$logFemaleMass, BirdEggs$logEggMass, 
     xlab="Log body mass", ylab="Log egg mass", 
     main="", lwd=2, col="black", pch=4)
polygon(ConfPoly.x, ConfPoly.y, col="wheat", border=NA)
abline(mod, col="darkblue", lwd=2)

# replicate(20, abline(a=rnorm(1, coef(mod)[1], summary(mod)$coefficients[1,2]), b=coef(mod)[2], col=4))
points(BirdEggs$logFemaleMass, BirdEggs$logEggMass, col="black", pch=4, lwd=2)

```

## Looking at uncertainty in the intercept

If we fix the slope, then we can look at variation in the intercept. These are the upper & lower 95% confidence limits

```{r PlotIntercepts2, echo=FALSE, results='hide'}
ConfIntervals <- confint(mod)

par(mar=c(4.1,4.1,1,1))
plot(BirdEggs$logFemaleMass, BirdEggs$logEggMass, 
     xlab="Log body mass", ylab="Log egg mass", 
     main="", lwd=2, col="black", pch=4)
polygon(ConfPoly.x, ConfPoly.y, col="wheat", border=NA)
abline(mod, col="darkblue", lwd=2)

replicate(20, abline(a=rnorm(1, coef(mod)[1], summary(mod)$coefficients[1,2]), b=coef(mod)[2], col=4))
points(BirdEggs$logFemaleMass, BirdEggs$logEggMass, col="black", pch=4, lwd=2)

```



## Uncertainty in the slope

If we fix the intercept, then we can look at variation in the slope. These are the uppr & lower 95% confidence limits

```{r PlotSlopes0, echo=FALSE, results='hide'}
Centre <- mean(BirdEggs$logEggMass) -
  ConfIntervals[2,]*mean(BirdEggs$logFemaleMass)

ConfPoly.y2 <- rep(Centre, each=2) + rep(ConfIntervals[2,], each=2)*ConfPoly.x
par(mar=c(4.1,4.1,1,1))
plot(BirdEggs$logFemaleMass, BirdEggs$logEggMass, 
     xlab="Log body mass", ylab="Log egg mass", 
     main="", lwd=2, col="black", pch=4)
polygon(ConfPoly.x, ConfPoly.y2, col="wheat", border=NA)
abline(mod, col="darkblue", lwd=2)

AddLines <- function(mod) {
  slope <- rnorm(1, coef(mod)[2], summary(mod)$coefficients[2,2])
  int <- mean(mod$model$logEggMass) - slope*mean(mod$model$logFemaleMass)
  abline(a=int, b=slope, col=4)
}
# replicate(10, AddLines(mod))
points(BirdEggs$logFemaleMass, BirdEggs$logEggMass, col="black", pch=4, lwd=2)


```


## Uncertainty in the slope

If we fix the intercept, then we can look at variation in the slope. These are the uppr & lower 95% confidence limits

```{r PlotSlopes, echo=FALSE, results='hide'}
Centre <- mean(BirdEggs$logEggMass) -
  ConfIntervals[2,]*mean(BirdEggs$logFemaleMass)

ConfPoly.y2 <- rep(Centre, each=2) + rep(ConfIntervals[2,], each=2)*ConfPoly.x
par(mar=c(4.1,4.1,1,1))
plot(BirdEggs$logFemaleMass, BirdEggs$logEggMass, 
     xlab="Log body mass", ylab="Log egg mass", 
     main="", lwd=2, col="black", pch=4)
polygon(ConfPoly.x, ConfPoly.y2, col="wheat", border=NA)
abline(mod, col="darkblue", lwd=2)

AddLines <- function(mod) {
  slope <- rnorm(1, coef(mod)[2], summary(mod)$coefficients[2,2])
  int <- mean(mod$model$logEggMass) - slope*mean(mod$model$logFemaleMass)
  abline(a=int, b=slope, col=4)
}
replicate(10, AddLines(mod))
points(BirdEggs$logFemaleMass, BirdEggs$logEggMass, col="black", pch=4, lwd=2)


```

## Uncertainty 

In reality both slope & intercept are uncertain

```{r PlotAll, echo=FALSE, results='hide'}

PredData <- data.frame(logFemaleMass=seq(-20, 20, length=50))
PredData$logEggMass <- predict(mod, newdata = PredData, interval="confidence")
par(mar=c(4.1,4.1,1,1))
plot(BirdEggs$logFemaleMass, BirdEggs$logEggMass, 
     xlab="Log body mass", ylab="Log egg mass", 
     main="", lwd=2, col="black", pch=4)
polygon(c(PredData$logFemaleMass, rev(PredData$logFemaleMass)), 
        c(PredData$logEggMass[,"lwr"], rev(PredData$logEggMass[,"upr"])), col="wheat", border=NA)
abline(mod, col="darkblue", lwd=2)

AddLinesBS <- function(mod) {
  dat <- mod$model[sample.int(nrow(mod$model), nrow(mod$model), replace = TRUE), ]
  abline(lm(logEggMass~logFemaleMass, data=dat), col=4)
}
replicate(10, AddLinesBS(mod))
points(BirdEggs$logFemaleMass, BirdEggs$logEggMass, col="black", pch=4, lwd=2)

```

## Summaries

```{r Summary, echo=TRUE, comment=NA}
summary(mod)
```

## Summaries: Coefficients

```{r SummaryCoef, echo=FALSE, comment=NA}
summ2 <- paste(capture.output(print(summary(mod), digits=2)), "\n")
cat(summ2[9:14])

# round(summary(mod)$coefficients, 2)
```


## Summaries: Last stuff

```{r SummaryLastStuff, echo=FALSE, comment=NA}
cat(summ2[16:18])
```
## Summaries: Last stuff

```{r SummaryLastStuff2, echo=FALSE, comment=NA}
knitr::asis_output("\\footnotesize")
cat(summ2[16])
```
- Estimate of $\sigma$
- Degrees of freedom: how many "spare" data points we have

```{r SummaryLastStuff3, echo=FALSE, comment=NA}
knitr::asis_output("\\footnotesize")
cat(summ2[17])
```
- Multiple R-squared: $R^2$ (see last lecture). What proportion of the variation are we explaining?
- Adjusted R-squared: ignore (at least for now)

```{r SummaryLastStuff4, echo=FALSE, comment=NA}
knitr::asis_output("\\footnotesize")
cat(summ2[18])
```
- Test of if the data explains anything. Usually very silly to test this.

## Pause


## Model Checking

Lots more plots!

## Building a Plot

```{r FittedResidPlot, echo=TRUE, fig.height=4, fig.width=6}
plot(fitted(mod), resid(mod))
```

(for your own plots, do what you feel comfortable with)

## Add axis labels

```{r FittedResidPlotAddLabels, echo=TRUE, fig.height=4, fig.width=6}
plot(fitted(mod), resid(mod), xlab="Fitted log egg mass", 
     ylab="Residual")
```

## Change the plot colour

```{r FittedResidPlotAddColour, echo=TRUE, fig.height=4, fig.width=6}
plot(fitted(mod), resid(mod), xlab="Fitted log egg mass", 
     ylab="Residual", col="red")
```

(col= is documented in ?par, as are a lot of other options)

## Quantile Plots

Special function!

```{r QQPlot, fig.height=3.5, fig.width=6}
par(mar=c(4.1,4.1,1,1), lwd=2)
qqnorm(resid(mod)) 
qqline(resid(mod)) 

```

(mar= sets margin sizes)

## Quantile Plots

Special function!

```{r QQPlotTitle, fig.height=3.5, fig.width=6}
par(mar=c(4.1,4.1,1,1), lwd=2)
qqnorm(resid(mod), main="A Plot") 
qqline(resid(mod)) 

```

(main= gives title)

## Now Into More Dimensions


## Multiple Regression

```{r GetData}
Data.all <- read.csv("../Data/EggBodySize.csv")
Data <- Data.all[!apply(Data.all, 1, function(x) any(is.na(x))),c("species", "development.category", "Egg.mass..g.", "female.mass..g.", "mm..n..mean.femur")]
Data$logEggMass <- log(Data$Egg.mass..g.)
Data$logFemaleMass <- log(Data$female.mass..g.)

```

We have learned how to draw straight lines

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

But this is limited as a model: we cannot draw more complicated curves, and we cannot explain or predict the effects of more than one covariate.

## The Model

This is out basic model

$$
\begin{aligned}
y_i &= \alpha + \sum_{j=1}^p \beta_j x_{ij} + \varepsilon_i \\
y_i &= \alpha + \beta_1 x_{i1} + \beta_2 x_{i2} + \beta_3 x_{i3} + \dots + \beta_p x_{ip} + \varepsilon_i
\end{aligned}
$$

So we replace $\beta x_j$ with $\sum_{j=1}^p \beta_j x_{ij}$. 

- we have $p$ covariates, labelled from $j=1$ to $p$
- we have $p$ covariate effects
- the j$^{th}$ covariate values for the i$^{th}$ individual is $x_{ij}$



## Design Matrices

We can write this more compactly. First, we turn the intercept into a covariate by using a covariate with a value of 1 for every data point. Then we write all of the covariates in a matrix, $X$:

$$
X = \left( \begin{array}{ccc}
   1 & 2.3 & 3.0 \\
   1 & 4.9 & -5.3 \\
   1 & 1.6 & -0.7 \\
\vdots & \vdots & \vdots \\
   1 & 8.4 & 1.2 \\
  \end{array}  \right)
$$

So, the first column is the intercept, the second is the first covariate, and the third is the second covariate.

This is called the *Design Matrix*: it is helpful for writing down the model

## Writing the Model 

Using marix algebra, the regression model becomes

$$
\mathbf{Y} = X \mathbf{\beta} + \mathbf{\varepsilon}
$$

where $\mathbf{Y}$, $\mathbf{\beta}$ and $\mathbf{\varepsilon}$ are now all vectors of length $n$, where there are $n$ data points. $X$ is am $n \times p$ matrix.

We will not look at the mathematics in any detail: the point here is that the model for the effect of covariates can be written in the design matrix.

## Writing the Model 

$$
\mathbf{Y} = X \mathbf{\beta} + \mathbf{\varepsilon}
$$

is 

$$
\left( \begin{array}{c}
   y_1 \\
   y_2 \\
   y_3 \\
\vdots \\
   y_n \\
  \end{array}  \right) = \left( \begin{array}{ccc}
   1 & 2.3 & 3.0 \\
   1 & 4.9 & -5.3 \\
   1 & 1.6 & -0.7 \\
\vdots & \vdots & \vdots \\
   1 & 8.4 & 1.2 \\
  \end{array}  \right)
  \left( \begin{array}{c}
   \beta_1 \\
   \beta_2 \\
   \beta_3 \\
\vdots \\
   \beta_n \\
  \end{array}  \right)
 +   \left( \begin{array}{c}
   \varepsilon_1 \\
   \varepsilon_2 \\
   \varepsilon_3 \\
\vdots \\
   \varepsilon_n \\
  \end{array}  \right)
$$

## The Solution (just so you can see it)

After a bit of matric algebra, one can find the ML solution:

$$
\mathbf{b} = (X^T X)^{-1}X^T \mathbf{Y}
$$
where $\mathbf{b}$ is the MLE for $\mathbf{\beta}$.

In practice, (a) you won't have to calculate this: the computer does it, and (b) the computer actually doesn't use this


## An Example: Bird Brains

The data wehave been using was collected to look at the effects of longevity on brain size. 

But size is a counfounder, i.e. it also has an effect

```{R GetBirdBrainData, fig.height=4.5, fig.height=4.5}
BirdBrains <- read.csv(file="../Data/BirdBrains.csv")
pairs(BirdBrains[,c("Maximum.lifespan",
                    "logBodyMass", "logBrainMass")])

```

## Fitting the Model

We wil ltry to explain (log) brain mass with:

- Maximum lifespan
- Age at first reprodction
- logBodyMass

We can write the model as

logBrainMass ~ Maximum lifespan +  Age at first reprodction + logBodyMass

## Fitting the Model

```{R FitMVModel}
Mod <- lm(logBrainMass ~ Maximum.lifespan +  logBodyMass,
          data=BirdBrains)

round(summary(Mod)$coefficients, 3)
```

## The Model

The model is

$$
y_i = `r round(coef(Mod)[1], 2)` + `r round(coef(Mod)[2], 2)` x_{i1} + 
`r round(coef(Mod)[2], 2)` x_{i2} + 
\varepsilon_i
$$

Where the X's are maximum lifespan and logBodyMass

## The Model

With one covariate we have a straight line, with two we have a plane.
```{R PlotBirdBrainResults, fig.width=4.5, fig.height=4.5}
  library(scatterplot3d)
Mod2 <- lm(logBrainMass ~ logBodyMass + Maximum.lifespan,
          data=BirdBrains)

s3d <- scatterplot3d(y=BirdBrains$Maximum.lifespan, 
       x=BirdBrains$logBodyMass, 
       z=BirdBrains$logBrainMass, bg=2,
       ylab="Maximum lifespan", xlab="log Body Mass",
       zlab="log Brain Mass", angle=55, pch=20)
  s3d$plane3d(Mod2, lty.box = "solid")
  
```

## Model Interpretation

The parameters say what happens if we increase a parameter by 1 if we hold the other covariates constant

e.g. if we increase maximum lifespan by 1 year, log brain mass increases by `r round(coef(Mod)['Maximum.lifespan'], 3)`.

But this is scale dependent: if we measure lifespan in decades, the coefficient changes to `r round(coef(Mod)['Maximum.lifespan']*10, 3)`

Makes it difficult to compare between different covariates: how do we compare 1 year to 1kg?

## Standardised Coefficients

We can look at the standardised coefficients

- standardise by the standard deviation
- also mean-centre

## Fit the Standardised Model

```{R FitStdModel, echo=TRUE}
BirdBrains$Max.lifespan.std <-
  scale(BirdBrains$Maximum.lifespan)
BirdBrains$logBodyMass.std <- 
  scale(BirdBrains$logBodyMass)
Mod.std <- lm(logBrainMass ~ logBodyMass.std + 
                Max.lifespan.std, data=BirdBrains)
round(summary(Mod.std)$coefficients, 3)
```

## Interpret the Standardised Model

```{R SummStdModel, echo=TRUE}
round(summary(Mod.std)$coefficients, 3)
```

- The intercept is now at the mean of Max. lifespan & log Body Mass
- the estimates are the effects of changing covriates by 1 standard deviation.
    - comparable as relative changes in the data
    - effect of log body mass about 10 times bigger than lifespan
    
## Plot the Standardised Model
    
```{R PlotStdBirdBrainResults, fig.width=4.5, fig.height=4.5}
s3d.std <- scatterplot3d(y=BirdBrains$Max.lifespan.std, 
       x=BirdBrains$logBodyMass.std, 
       z=BirdBrains$logBrainMass, bg=2,
       ylab="Maximum lifespan", xlab="log Body Mass",
       zlab="log Brain Mass", angle=55, pch=20)
  s3d.std$plane3d(Mod.std, lty.box = "solid")
  
```

## Another use of multiple regression

Approximating curves

We can approximate curves with a Taylor series:

$$
f(x) \approx \beta_0 + \beta_1 (x- \bar{x}) + \beta_2 (x- \bar{x})^2 + \beta_3 (x- \bar{x})^3 + \dots + \beta_p (x- \bar{x})^p
$$

So we can fit an approximate curve by regressing $Y$ against $X$, $X^2$, $x^3$ etc.

## Bird Eggs

For the bird eggs data, there might be a curve, so we can fit that

```{r BirdEggs, echo = TRUE}
BirdEggs <- read.csv(file="../Data/BirdEggs.csv",
                     stringsAsFactors = FALSE)
BirdEggs$lgFM.std <- scale(BirdEggs$logFemaleMass)
Mod.quad <- lm(logEggMass ~ lgFM.std + I(lgFM.std^2),
               data=BirdEggs)
round(summary(Mod.quad)$coefficients, 2)
```

## Plot Bird Eggs

```{r PlotBirdEggs, echo = FALSE}
Pred.eggs <- data.frame(
  logFemaleMass=seq(min(BirdEggs$logFemaleMass),
                    max(BirdEggs$logFemaleMass), length=50))
Pred.eggs$lgFM.std <- 
  scale(Pred.eggs$logFemaleMass, 
        center=attr(BirdEggs$lgFM.std, "scaled:center"),
        scale=attr(BirdEggs$lgFM.std, "scaled:scale"))
Pred.eggs <- cbind(Pred.eggs, 
                   predict(Mod.quad, newdata = Pred.eggs,
                        interval="prediction"))
plot(BirdEggs$logFemaleMass, BirdEggs$logEggMass, type="n",
     xlab="log Female Mass", ylab="log Egg Mass")
polygon(c(Pred.eggs$logFemaleMass,
          rev(Pred.eggs$logFemaleMass)),
        c(Pred.eggs$lwr, rev(Pred.eggs$upr)), col ="pink",
        border = NA)
points(BirdEggs$logFemaleMass, BirdEggs$logEggMass, pch=20)
lines(Pred.eggs$logFemaleMass, Pred.eggs$fit, col=2)

```

## Summary

We can fit models with more than 1 covariate

- comparison of the coefficients is a bit tricky
- with more than 2 covariates, plotting the model is a pain

We can use this to fit more complicated curves

## Next Week

How well does the model fit?

Do we need all of these parameters?