---
title: "Lecture 7: Multiple Regression"
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...

Exercises to be handed in by 17:00 

- if we get a folder set up


## 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.

## 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, so it had to be included

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

```

## Another Example: Health care

The health care data isn't linear. We can try transforming, or we can fit a polynomial (i.e. include $x^2$, $x^3$ etc. terms)

```{R GetHealthData, fig.height=4.5, fig.height=4.5}
rawdata <- read.csv("../Data/LifeExpectancy.csv")
par(mar=c(4.1,4.1,1,1))
plot(rawdata$Health.Spending.per.capita, rawdata$Life.exectancy, type="n",
     xlab="Health Spending per capita", ylab="Life expectancy (years)")
text(rawdata$Health.Spending.per.capita, rawdata$Life.exectancy, rawdata$Country)
```

## 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_0 \\
   \beta_1 \\
   \beta_2 \\
\vdots \\
   \beta_p \\
  \end{array}  \right)
 +   \left( \begin{array}{c}
   \varepsilon_1 \\
   \varepsilon_2 \\
   \varepsilon_3 \\
\vdots \\
   \varepsilon_n \\
  \end{array}  \right)
$$

- $\beta_0$ is the intercept

## 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:

- you won't have to calculate this: the computer does it, and 
- the computer actually doesn't use this



## Fitting the Model

We will try to explain (log) brain mass with:

- Maximum lifespan
- logBodyMass

We can write the model as

logBrainMass ~ Maximum lifespan + logBodyMass

## Fitting the Model

Parameter Estimates

```{R FitMVModel, echo=TRUE}
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)[3], 2)` x_{i2} + 
\varepsilon_i
$$

Where the X's are maximum lifespan and logBodyMass

```{R SummMVModel, echo=FALSE}
round(summary(Mod)$coefficients, 3)
```

## 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(x=BirdBrains$logBodyMass, 
                     y=BirdBrains$Maximum.lifespan, 
       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

```{R StdData, echo=TRUE}
BirdBrains$Max.lifespan.std <-
  scale(BirdBrains$Maximum.lifespan)
BirdBrains$logBodyMass.std <- 
  scale(BirdBrains$logBodyMass)
```

## Fit the Standardised Model

```{R FitStdModel, echo=TRUE}
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 covariates 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

```{r BirdEggsData, echo = FALSE}
BirdEggs <- read.csv(file="../Data/BirdEggs.csv",
                     stringsAsFactors = FALSE)

```

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

```{r BirdEggs, echo = TRUE}
BirdEggs$lgFM.std <- scale(BirdEggs$logFemaleMass) # standardise
Mod.quad <- lm(logEggMass ~ lgFM.std + I(lgFM.std^2),
               data=BirdEggs)
```

The model formula is 

logEggMass ~ lgFM.std + I(lgFM.std^2)

We need the I() to tell R to use the quadratic

## Bird Eggs Summary

```{r SummBirdEggs, echo = TRUE}
round(summary(Mod.quad)$coefficients, 2)
```

We still see a positive linear term, but the quadratic is negative

- so, what does the curve look like?

## 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

Categorical Variables


## 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),]
mod <- lm(logEggMass ~ logFemaleMass, data=BirdEggs)

```

## 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)


## Add a Title

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