## Simple Linear Regression Quiz

Topics:

- Visualizing data.
- Fitting a linear model.
- Assessing whether the assumptions are reasonable.
- Inference (CIs, PIs)

The quiz will be posted this evening and due Friday evening.

You may utlize resources such as the book, the internet, notes, slides, etc.
You must do the work yourself however; no asking questions of other students,
forums, etc.  Please email me if you any questions come up.


## Consider fitted models to four different data sets


##

```{r}
summary(m1 <- lm(y1 ~ x1, data = anscombe))
```


##

```{r}
summary(lm(y2 ~ x2, data = anscombe))
```

##

```{r}
summary(lm(y3 ~ x3, data = anscombe))
```


##

```{r}
summary(lm(y4 ~ x4, data = anscombe))
```


## Anscombe's Quartet

```{r, echo=FALSE}
ff <- y ~ x
mods <- setNames(as.list(1:4), paste0("lm", 1:4))
for(i in 1:4) {
  ff[2:3] <- lapply(paste0(c("y","x"), i), as.name)
  mods[[i]] <- lmi <- lm(ff, data = anscombe)
}
op <- par(mfrow = c(2, 2), mar = 0.1+c(4,4,1,1), oma =  c(0, 0, 2, 0))
for(i in 1:4) {
  ff[2:3] <- lapply(paste0(c("y","x"), i), as.name)
  plot(ff, data = anscombe, col = "red", pch = 21, bg = "orange", cex = 1.2,
       xlim = c(3, 19), ylim = c(3, 13))
  abline(mods[[i]], col = "blue")
}
par(op)
```


## For a valid model we need

1. The conditional mean of Y|X is a linear function of X.
2. The variance of Y|X is the same for any X.
3. The errors (and thus the Y|X are independent of one another).
4. The errors are normally distributed with mean zero.

5*. No "outliers".

These can all be assessed with residual plots.


## Basic residual plots

How to construct:

1. Calculate $\hat{e}_i = y_i - \hat{y}_i$ for each point in your data set (also
available as `m1$res`).
2. Create a scatter plot with the residuals on the y axis.  On the x-axis you
can plot either the x-coordinates or the fitted values (\hat{y}_i).

```{r, eval=FALSE}
plot(anscombe$x1, m1$res) # x versus residuals
plot(m1$fit, m1$res) # x versus fitted
plot(m1, 1) # built-in function
```


## 

```{r, echo=FALSE}
ff <- y ~ x
mods <- setNames(as.list(1:4), paste0("lm", 1:4))
for(i in 1:4) {
  ff[2:3] <- lapply(paste0(c("y","x"), i), as.name)
  mods[[i]] <- lmi <- lm(ff, data = anscombe)
}
op <- par(mfrow = c(2, 2), mar = 0.1+c(4,4,1,1), oma =  c(0, 0, 2, 0))
for(i in 1:4) {
  plot(mods[[i]], 1)
  #ff[2:3] <- lapply(paste0(c("y","x"), i), as.name)
  #plot(ff, data = anscombe, col = "red", pch = 21, bg = "orange", cex = 1.2,
  #     xlim = c(3, 19), ylim = c(3, 13))
  #abline(mods[[i]], col = "blue")
}
par(op)
```


## Anscombe I {.build}

1. The conditional mean of Y|X is a linear function of X.
*OK!*

2. The variance of Y|X is the same for any X.
*OK!*

3. The errors (and thus the Y|X are independent of one another).
*OK!*

4. The errors are normally distributed with mean zero.
*Probably ok?*

5. No "outliers".
*OK!*


## 

```{r, echo=FALSE}
ff <- y ~ x
mods <- setNames(as.list(1:4), paste0("lm", 1:4))
for(i in 1:4) {
  ff[2:3] <- lapply(paste0(c("y","x"), i), as.name)
  mods[[i]] <- lmi <- lm(ff, data = anscombe)
}
op <- par(mfrow = c(2, 2), mar = 0.1+c(4,4,1,1), oma =  c(0, 0, 2, 0))
for(i in 1:4) {
  plot(mods[[i]], 1)
  #ff[2:3] <- lapply(paste0(c("y","x"), i), as.name)
  #plot(ff, data = anscombe, col = "red", pch = 21, bg = "orange", cex = 1.2,
  #     xlim = c(3, 19), ylim = c(3, 13))
  #abline(mods[[i]], col = "blue")
}
par(op)
```


## Anscombe II {.build}

1. The conditional mean of Y|X is a linear function of X.
*No way!  Looks quadratic*

2. The variance of Y|X is the same for any X.
*OK!*

3. The errors (and thus the Y|X are independent of one another).
*Could be an issue, but probably not.*

4. The errors are normally distributed with mean zero.
*Doesn't look like it.*

5. No "outliers".
*Mmm, 8 and 6 maybe?*


## 

```{r, echo=FALSE}
ff <- y ~ x
mods <- setNames(as.list(1:4), paste0("lm", 1:4))
for(i in 1:4) {
  ff[2:3] <- lapply(paste0(c("y","x"), i), as.name)
  mods[[i]] <- lmi <- lm(ff, data = anscombe)
}
op <- par(mfrow = c(2, 2), mar = 0.1+c(4,4,1,1), oma =  c(0, 0, 2, 0))
for(i in 1:4) {
  plot(mods[[i]], 1)
  #ff[2:3] <- lapply(paste0(c("y","x"), i), as.name)
  #plot(ff, data = anscombe, col = "red", pch = 21, bg = "orange", cex = 1.2,
  #     xlim = c(3, 19), ylim = c(3, 13))
  #abline(mods[[i]], col = "blue")
}
par(op)
```


## Anscombe III {.build}

1. The conditional mean of Y|X is a linear function of X.
*Perfectly linear, but we've fit the wrong line!*

2. The variance of Y|X is the same for any X.
*Looks like a problem.*

3. The errors (and thus the Y|X are independent of one another).
*Hard to tell.*

4. The errors are normally distributed with mean zero.
*Hard to tell.*

5. No "outliers".
*Whoops!  Point 3 is a glaring problem.*


## 

```{r, echo=FALSE}
ff <- y ~ x
mods <- setNames(as.list(1:4), paste0("lm", 1:4))
for(i in 1:4) {
  ff[2:3] <- lapply(paste0(c("y","x"), i), as.name)
  mods[[i]] <- lmi <- lm(ff, data = anscombe)
}
op <- par(mfrow = c(2, 2), mar = 0.1+c(4,4,1,1), oma =  c(0, 0, 2, 0))
for(i in 1:4) {
  plot(mods[[i]], 1)
  #ff[2:3] <- lapply(paste0(c("y","x"), i), as.name)
  #plot(ff, data = anscombe, col = "red", pch = 21, bg = "orange", cex = 1.2,
  #     xlim = c(3, 19), ylim = c(3, 13))
  #abline(mods[[i]], col = "blue")
}
par(op)
```


## Anscombe IV {.build}

1. The conditional mean of Y|X is a linear function of X.
*Possibly, though the X doesn't appear to have much predictive power.*

2. The variance of Y|X is the same for any X.
*Hard to say.*

3. The errors (and thus the Y|X are independent of one another).
*Hard to tell.*

4. The errors are normally distributed with mean zero.
*Looks more uniform.*

5. No "outliers".
*The slope of the line is being completely determined by one point!*


## Assessing Normality {.build}

We can check the assumption that the errors are normal by looking at the distribution
of the residuals.  Difficult to do in a residual plot, so we use a QQ plot (for
quantile-quantile), aka normal probability plot.

*Quantile*: The $j^{th}$ quantile, $q_j$, is the value of a random variable $X$ 
that fulfills:

\[ P(X \le q_j) = j \]

For the standard normal distribution, $q_{.5} = 0$, $q_{.025} = -1.96$, 
$q_{.975} = 1.96$.


## Constructing a QQ plot {.build}

1. Standardize your residuals.
\[ \tilde{e}_i = \frac{\hat{e}_i - \bar{\hat{x}}}{s} \]

2. If you have $n$ standardized residuals, you can consider the lowest to be
the $1/n$ quantile, the second lowest, the $2/n$ quantile, the median to be the
$.5$ quantile, etc.

3. Look up these values for the standard normal distribution and find what the 
appropriate quantiles would be (this is what `qnorm()` does). These become your 
theoretical quantiles.

4. Plot the theoretical quantiles against the standardized residuals.


##

```{r, fig.height=5, fig.width = 5, fig.align='center'}
plot(m1, 2)
```


## Interpreting a QQ plot

- Perfectly normally distributed residuals would align along the identify line.
- Short tails will veer of the line horizontally.
- Long tails will veer off the line vertically.
- *Expect some variation, even from normal residuals!*


## Normal residuals

```{r fig.height=5, fig.width = 5, fig.align='center'}
x <- rnorm(40)
qqnorm(x)
qqline(x)
```


## Heavy tailed residuals

```{r fig.height=5, fig.width = 5, fig.align='center'}
x <- rt(40, 1)
qqnorm(x)
qqline(x)
```


## Returning to Anscome 1

```{r echo=FALSE}
par(mfrow = c(1, 2))
plot(m1, 1:2)
```


## Checking constant variance

Recall the quakes data:

```{r, echo=FALSE}
plot(stations ~ mag, data = quakes)
m1 <- lm(stations ~ mag, data = quakes)
abline(m1, col = "orange")
```


##

```{r, echo=FALSE}
par(mfrow = c(1, 2))
plot(m1, 1:2)
```


## 

```{r fig.height=5, fig.width = 5, fig.align='center'}
plot(m1, 3)
```


## Checking constant variance {.build}

- The scale location plot transforms the residuals to make non-constant
variance (heteroscedasticity) more apparent.
- The red line is a guide: flat = constant variance.
- The basic residual plot can also be used, but it's a bit more difficult to tell
(also, *that* red line refers to the linear trend).


## Activity 4

- Please work in pairs, trios.
- Submit whatever you have by the end of class via moodle and include all group
members in file name.

