---
title: "Lecture 5: How Good is Our Regression"
author: "(Grafen & hails Chapter 9)"
date: "bob.ohara@ntnu.no"
output: beamer_presentation
---

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


## Regression: How good is the model?

```{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.)

```

Last week we learned how to draw straight lines


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


```{r PlotLines1}
mod <- lm(logEggMass ~ logFemaleMass, data=Data)
ss <- var(Data[,c("logFemaleMass", "logEggMass")]); delta <- 1
slope <- (ss[2,2] - delta*ss[1,1] + sqrt((ss[2,2]-delta*ss[1,1])^2 +4*delta*ss[1,2]^2))/(2*ss[1,2])
intercept <- mean(Data$logEggMass) - slope*mean(Data$logFemaleMass)

par(mar=c(4.1,4.1,1,1), lwd=2)
plot(Data$logFemaleMass, Data$logEggMass, xlim=c(4,8), ylim=c(0,6.2),
     xlab="Log body mass", ylab="Log egg mass", 
     main="", lwd=3, col="lightgreen")
abline(a=intercept, b=slope)

```


## How well does it explain the data?

How could we summarise these?

```{R Variation, fig.height=6}

x <- 1:100
yGood <- rnorm(length(x), x, 5)
yPoor <- rnorm(length(x), x, 20)
yBad <- rnorm(length(x), x, 100)
yCrap <- rnorm(length(x), 0.0005*x^3, 50)

par(mfrow=c(2,2), mar=c(2,2,1,1), oma=c(2,2,0,0))
plot(x, yGood)
plot(x, yPoor)
plot(x, yBad)
plot(x, yCrap)
mtext("x", 1, outer=TRUE)
mtext("y", 2, outer=TRUE)
```

## How well does it explain the data? R^2

One way is to ask what proportion of the total variation is explained by the model

Informally

$$
\frac{Var(\alpha + \beta X)}{Var(Y)} = 1 - \frac{Var(\varepsilon)}{Var(Y)}
$$

## More formally 

$$
Var(y) = \frac{1}{N-1} \sum(y_i - \bar{y})^2
$$

$$
Var(E(y)) = \frac{1}{N-1} \sum(\alpha + \beta x_i - \bar{y})^2 = \frac{1}{N-1} \beta^2 \sum(x_i - \bar{x})^2
$$

$$
Var(\varepsilon) = \frac{1}{N-1} \sum(y_i - (\alpha + \beta x_i))^2
$$

In practice, we calculate it as

$$
R^2 = 1 - \frac{(\sum y_i - \mu_i)^2}{(\sum y_i - \bar{y})^2}
$$

## What is a good $R^2$?

It depends!

```{R DifferentR2, fig.height=6}

x <- scale(1:100)
eps <- scale(rnorm(length(x), 0, 1))
prop <- c(0.1, 0.5, 0.7, 0.9)
y10 <- sqrt(prop[1])*x + sqrt(1-prop[1])*eps
y50 <- sqrt(prop[2])*x + sqrt(1-prop[2])*eps
y70 <- sqrt(prop[3])*x + sqrt(1-prop[3])*eps
y90 <- sqrt(prop[4])*x + sqrt(1-prop[4])*eps

par(mfrow=c(2,2), mar=c(2,2,3,1), oma=c(2,2,0,0), cex.main=2)
plot(x, y10, main=expression(R^2==10~'%'), xlab="", ylab="")
plot(x, y50, main=expression(R^2==50~'%'), xlab="", ylab="")
plot(x, y70, main=expression(R^2==70~'%'), xlab="", ylab="")
plot(x, y90, main=expression(R^2==90~'%'), xlab="", ylab="")
mtext("x", 1, outer=TRUE)
mtext("y", 2, outer=TRUE)
```

## Model Checking

How could the model go wrong?

There are lots of ways the model go wrong 

Discuss and list at least 3

- draw a plot of how it might look

The assumptions:

- linear in covariate
- error is normally distributed, with equal variance

## Graphical Checks

The assumption is that the error term $\varepsilon_i$ is unstructured error. So if we see structure, we might have a problem.

Easier to see what is going on with graphs: not just if there is a problem, but where it is


## Residuals

An important statistic we use is the _residual_: the difference between the observed and expected values:

$$
r_i = y_i - E(y_i) = y_i - (\hat{\alpha} + \hat{\beta} X_i)
$$

If $\hat{\alpha} = \alpha$ and $\hat{\beta} = \beta$, then $r_i = \varepsilon_i$

We also sometimes standardise them:

$$
t_i = \frac{y_i - E(y_i)}{\sqrt{var(r_i)}}
$$

## Why they are useful

The assumption is that the error term $\varepsilon_i$ is unstructured error. So if we see structure, we might have a problem.

So we should check this

## Plot against the predicted value

Plot $r_i$ against the predicted value, $\hat{\mu_i} = \hat{\alpha} + \hat{\beta} x_i$

```{r ResidPlot1, fig.height=4}
par(mar=c(4.1,4.1,1,1), lwd=2)
plot(predict(mod), resid(mod), 
     xlab="predicted log egg mass", ylab="residual", 
     main="", lwd=3, col="lightblue")

```


## Plot against the covariate

Plot $r_i$ against $x_i$

```{r ResidPlot2, fig.height=4}
par(mar=c(4.1,4.1,1,1), lwd=2)
plot(Data$logFemaleMass[!is.na(Data$logFemaleMass)], resid(mod), 
     xlab="Log body mass", ylab="residual", 
     main="", lwd=3, col="lightblue")

```

(same as last plot, but this will become useful with multiple regression)


## Anscombe's Quartet: when things gow rong

```{r AnscombeData, fig.height=5}
par(mfrow=c(2,2), mar=c(2,2,2,2), oma=c(0,0,0,0), pch=16)
plot(anscombe$x1, anscombe$y1, xlab="x", ylab="y")
plot(anscombe$x2, anscombe$y2, xlab="x", ylab="y")
plot(anscombe$x3, anscombe$y3, xlab="x", ylab="y")
plot(anscombe$x4, anscombe$y4, xlab="x", ylab="y")

```

## Anscombe's Residuals: 1

```{r AnscombeResid1, fig.height=4}

anscombe1 <- lm(y1~x1, data=anscombe)
anscombe2 <- lm(y2~x2, data=anscombe)
anscombe3 <- lm(y3~x3, data=anscombe)
anscombe4 <- lm(y4~x4, data=anscombe)

par(mfrow=c(1,2), mar=c(4,4,3,1), oma=c(0,0,0,0), pch=16)
plot(anscombe$x1, anscombe$y1, xlab="x", ylab="y", main="Data")
abline(lm(y1~x1, data=anscombe), col=2)
plot(fitted(anscombe1), resid(anscombe1), xlab="Fitted", ylab="Residuals", main="Residuals")

```

These look OK: no structure

## Anscombe's Residuals: 2

```{r AnscombeResid2, fig.height=4}
par(mfrow=c(1,2), mar=c(4,4,3,1), oma=c(0,0,0,0), pch=16)
plot(anscombe$x2, anscombe$y2, xlab="x", ylab="y", main="Data")
abline(lm(y2~x2, data=anscombe), col=2)
plot(fitted(anscombe2), resid(anscombe2), xlab="Fitted", ylab="Residuals", main="Residuals")

```

:-(

## Anscombe's Residuals: 3

```{r AnscombeResid3, fig.height=5}
par(mfrow=c(1,2), mar=c(4,4,3,1), oma=c(0,0,0,0), pch=16)
plot(anscombe$x3, anscombe$y3, xlab="x", ylab="y", main="Data")
abline(lm(y3~x3, data=anscombe), col=2)
plot(fitted(anscombe3), resid(anscombe3), xlab="Fitted", ylab="Residuals", main="Residuals")

```

The outlier stands out, with a large residual


## Anscombe's Residuals: 4

```{r AnscombeResid4, fig.height=5}
par(mfrow=c(1,2), mar=c(4,4,3,1), oma=c(0,0,0,0), pch=16)
plot(anscombe$x4, anscombe$y4, xlab="x", ylab="y", main="Data")
abline(lm(y4~x4, data=anscombe), col=2)
plot(fitted(anscombe4), resid(anscombe4), xlab="Fitted", ylab="Residuals", main="Residuals")

thing <- resid(anscombe4)[anscombe$x4>10]
```

The outlier actually has a small residual (`r as.character(thing)`).

## Influential Points

The fitted line goes through the point on the right

```{r AnscombeData4a, fig.height=5}
par(mfrow=c(1,1), mar=c(2,2,2,2), oma=c(0,0,0,0), pch=16)
plot(anscombe$x4, anscombe$y4, xlab="x", ylab="y")
abline(lm(y4~x4, data=anscombe))
```


## Influential Points

Here the outlier pulls the whole line up towards it

```{r AnscombeData4b, fig.height=5}
par(mfrow=c(1,1), mar=c(2,2,2,2), oma=c(0,0,0,0), pch=16)
plot(anscombe$x3, anscombe$y3, xlab="x", ylab="y")
abline(lm(y3~x3, data=anscombe))
```

## What is an influnetial point?

If we remove it, the model will change a lot


```{r AnscombeData3Influence, fig.height=5, results='hide'}

RemovePoint <- function(which, x, y) {
  if(length(x)!=length(y)) stop("x & y should be same length")
#  if(!is.integer(which)) stop("which should be an integer")
  if(which>length(x)) stop("which should not be higher than length(x)")
  if(which>length(y)) stop("which should not be higher than length(y)")
  mod <- lm(y[-which]~ x[-which])
  c(intercept=unname(coef(mod)["(Intercept)"]),
    slope=unname(coef(mod)["x[-which]"]), confint(mod)["x[-which]",])
}
# RemovePoint(which = 1, x=anscombe$x3, y=anscombe$y3)

Fits <- sapply(seq_along(anscombe$x3),  RemovePoint, x=anscombe$x3, y=anscombe$y3)

par(mfrow=c(2,1), mar=c(2,4,2,2), oma=c(2,0,0,0), pch=16)
plot(anscombe$x3, anscombe$y3, xlab="x", ylab="y")
apply(Fits, 2, function(x) 
  abline(a=x["intercept"], b=x["slope"], col="red"))
abline(lm(y3~x3, data=anscombe), lwd=2)
plot(anscombe$x3, Fits["slope",], xlab="x", ylab="slope")
```


## Influence and Leverage

We can generalise this idea by asking how much the fitted values for the other points change if we remove a data point

$$
D_i = \frac{\sum_{j=1}^{n} (\hat{y_j} - \hat{y}_{j(-i)})^2}{s^2}
$$
where $s^2 = 1/(n-1) \sum r_i^2$ is the _mean squared error_

## What is influential?

Large values of $D_i$ mean a large influence

- $D_i$ > 1, or 4/n?

```{r AnscombeData3CooksD, fig.height=5, results='hide'}
par(mfrow=c(2,1), mar=c(2,4,2,2), oma=c(2,0,0,0), pch=16)
plot(anscombe$x3, anscombe$y3, xlab="x", ylab="y")
plot(anscombe$x3, cooks.distance(anscombe3), xlab="x", ylab="Cook's D")
```


## Probability Plots

We are assuming that the residuals are normally distributed

We can check this with normal probability plots

Plot the residual against its expected value (i.e. if the model is correct)

## Probability Plots

If we sort the data (smallest to largest), we can plot them against their expected values, i.e. plot $r_i$ against the normal quantile

```{r QQResidPlot, fig.height=4}
par(mar=c(4.1,4.1,1,1), lwd=2)
qqnorm(resid(mod), main="", lwd=3, col="lightblue")
qqline(resid(mod))

```

## Constructing Probability Plots

```{r ContructQQPlot, fig.height=4}
Resids <- resid(mod)
NormQuants <- qnorm(1:length(Resids)/(1+length(Resids)))

par(mfrow=c(1,3), mar=c(4.1,2.1,3,1), oma=c(0,2,0,0), lwd=2)
plot(Resids, lwd=3, col="lightblue", main="Residuals", ylab="")
plot(sort(Resids), lwd=3, col="lightblue", 
     main="Sorted Residuals", ylab="", xlab="Rank")
plot(NormQuants, sort(Resids), lwd=3, col="lightblue",
     main="Residuals vs Normal quantiles", ylab="", xlab="Normal Quantiles")
mtext("Resisuals", 2, outer=TRUE)
```

## What you can see

```{r UseQQPlot, fig.height=4}
Res.norm <- rnorm(100)
Res.outl <- c(7, Res.norm[-1])
Res.skew <- Res.norm^1.2
Res.t <- Res.norm*rgamma(length(Res.norm), 3,3)

par(mfrow=c(1,4), mar=c(4.1,2.1,3,1), oma=c(0,2,0,0), lwd=2)
qqnorm(Res.norm, main="Normal", ylab="", col=4)
qqline(Res.norm)

qqnorm(Res.outl, main="Outliers", ylab="", col=4)
qqline(Res.outl)

qqnorm(Res.skew, main="Skewness", ylab="", col=4)
qqline(Res.skew)

qqnorm(Res.t, main="Thick tails", ylab="", col=4)
qqline(Res.t)
```



## What to Do when your model is wrong

Break

## Individual Data Points

Is your data point wrong?

- typos?
- real but unique

If it is wrong, correct, if it is right, might want to remove it & see if that makes a big difference

- if it does, be careful!

## Assumptions

![](LifeExpectancy.gif)

## Assumptions

The bigger picture

- curve

- skew

- heteroscedasticity


## Problems

```{R AssumptionProblems}
Data.res <- data.frame(
  x = scale(Data$logFemaleMass),
  y = scale(Data$logEggMass)
)
mod0 <- lm(y~x, data=Data.res)
Data.res$Resid <- resid(mod0)
Data.res$Resid[Data.res$Resid<0.65] <- Data.res$Resid[Data.res$Resid<0.65] + 0.5
Data.res$E1 <- coef(mod0)["(Intercept)"] +
  Data.res$x*coef(mod0)["x"]
Data.res$y1 <- Data.res$E1 + Data.res$Resid
mod1 <- lm(y1~x, data=Data.res)

# Quadratic
Data.res$E2 <- Data.res$E1 + 0.3*Data.res$x^2
Data.res$y2 <- Data.res$E2 + Data.res$Resid
mod2 <- lm(y2~x, data=Data.res)
Data.res$Resid2 <- resid(mod2)

# Heteroscedastic
Data.res$E3 <- fitted(mod1)
Data.res$y3 <- Data.res$E3 + 
  (Data.res$x-min(Data.res$x)+0.1)*Data.res$Resid
mod3 <- lm(y3~x, data=Data.res)
Data.res$Resid3 <- resid(mod3)

Data.res$E4 <- fitted(mod1)
Data.res$y4 <- Data.res$E4 + 10*(Data.res$Resid^2)
mod4 <- lm(y4~x, data=Data.res)
Data.res$Resid4 <- resid(mod4)

par(mfcol=c(2,4), mar=c(2,2,4,1), ann=FALSE)
plot(Data.res$x, Data.res$y1)
mtext("Normal Data", 3, line=1, cex=1.5)
plot(Data.res$x, Data.res$Resid)

plot(Data.res$x, Data.res$y2)
mtext("Curved Response", 3)
plot(Data.res$x, Data.res$Resid2)

plot(Data.res$x, Data.res$y3)
mtext("Heteroscedastic", 3)
plot(Data.res$x, Data.res$Resid3)

plot(Data.res$x, Data.res$y4)
mtext("Skewed Residuals", 3)
plot(Data.res$x, Data.res$Resid4)

```

## Possible Solutions

Transform the covariate

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

e.g. $\sqrt(x_i)$, $x_i^2$, $\log(x_i)$, 

Add more terms

- quadratic

$$
y_i = \alpha + \beta x_i + \gamma x_i^2  + \varepsilon_i
$$

More about this later


## Transformations

Transform the response

e.g. $\sqrt(x_i)$, $x_i^2$, $\log(x_i)$

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

## Box-Cox transformations

General Class of transformations

$$
y_i \rightarrow y_i^p 
$$

if $p = 0$, use $log(y_i)$

## Using Box-Cox transformations

```{r UseBoxCox, fig.height=4}
Data.res$y.bc <- Data.res$y + 3
Data.res$y.bc2 <- Data.res$y.bc^2
Data.res$y.bc0.5 <- Data.res$y.bc^0.5
Data.res$y.bclog <- log(Data.res$y.bc)

par(mfrow=c(1,4), mar=c(4.1,2.1,3,1), oma=c(0,2,0,0), lwd=2)
plot(Data.res$x, Data.res$y.bc, main="Untransformed")
 abline(lm(y.bc~x, data = Data.res), col=2)
plot(Data.res$x, Data.res$y.bc2, main="Squared Transformation")
 abline(lm(y.bc2~x, data = Data.res), col=2)
plot(Data.res$x, Data.res$y.bc0.5, main="Square Root Transformation")
 abline(lm(y.bc0.5~x, data = Data.res), col=2)
plot(Data.res$x, Data.res$y.bclog, main="log Transformation")
 abline(lm(y.bclog~x, data = Data.res), col=2)

```

## Heteroscedasticity

Variance changes with the mean

- Box-Cox can also solve this

```{r Heteroscedasticity, fig.height=4}

het.x <- rep(1:100,2)
a=2
b = 1
sigma2 = het.x^1.3
eps = rnorm(length(het.x),mean=0,sd=sqrt(sigma2))
het.y=a+b*het.x + eps
mod.het <- lm(het.y ~ het.x)
mod.het2 <- lm(sqrt(het.y) ~ het.x)

par(mfrow=c(1,3), mar=c(4.1,2.1,3,1), oma=c(0,2,0,0), lwd=2)
plot(het.y, het.x, main="Data")
abline(mod.het, col=2)

plot(fitted(mod.het), resid(mod.het), main="Residuals")
plot(fitted(mod.het2), resid(mod.het2), main="Residuals, sqrt transformation")

```

## Next Week

Multiple Regression

- more than 1 explanatory variable

- Grafen & Hails Chapter 4




