chapter_1.Rmd

---
title: "Chapter 1"
date: "February 19, 2018"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, comment = NA, warning = FALSE, message = FALSE)

library(tidyverse)
library(broom)
```

## Exercise 1

We know that distance should increase linearly with time, where zero time corresponds to zero distance. 

```{r}
df1 <- data.frame(
  length = c(1, 3, 4, 5), 
  hours  = c(0.1, 0.4, 0.5, 0.6)
)

summary(lm(length ~ 0 + hours, data = df1))$coefficients
```

Although `length / hours = speed`, the following model yields different results.

```{r}
summary(lm(length / hours ~ 1, data = df1))$coefficients
```

This latter model corresponds to the mean speed of each journey.

```{r}

df1 %>% 
  summarise(
    mean(length / hours)
  )
```

## Exercise 2

The linear model $y_i = \beta + \epsilon_i$ does not depend on the $x_i$, so $\beta = \bar y$.

## Exercise 3

Using linearity of $\mathbb E$, we can show that $\hat \beta$ is unbiased:

$$
\begin{align}
  \mathbb E (\hat\beta)
  &=
  \mathbb E ( (X^T X)^{-1} X^T \hat y)
  \\
  &=
  (X^T X)^{-1} X^T \mathbb E (\hat y)
  \\
  &=
  (X^T X)^{-1} X^T \mathbb E (X \beta + \epsilon)
  \\
  &=
  (X^T X)^{-1} X^T (X \beta + \mathbb E (\epsilon))
  \\
  &=
  (X^T X)^{-1} X^T X \beta
  \\
  &=
  \beta.
\end{align}
$$

This derivation does not depend on

1. independence of $Y_i$;
2. homoscedasticity of $Y_i$; or
3. normality of $Y_i$.

## Exercise 4

a. The model $y_{ij} = \alpha + \beta_i + \epsilon_{ij}$ can be written as $y_{ij} = X_{ji} \beta_i + \epsilon_{ij}$, where

$$
X
=
\begin{pmatrix}
  1 & 0 & 0 \\
  1 & 0 & 0 \\
  0 & 1 & 0 \\
  0 & 1 & 0 \\
  0 & 0 & 1 \\
  0 & 0 & 1
\end{pmatrix}
,
\qquad
\beta
=
\begin{pmatrix}
  \beta_1 \\
  \beta_2 \\
  \beta_3
\end{pmatrix}
$$

and $i = 1, 2, 3$ and $j = 1, 2$. The columns of $X$ are clearly linearly independent and the model is therefore identifiable.

b. The model $y_{ij} = \alpha + \beta_i + \gamma_j + \epsilon_{ij}$ can be written as $y_{ij} = X_{ji} b_i + \epsilon_{ij}$, where $\alpha = 0$ and

$$
X 
=
\begin{pmatrix}
  1 & 0 & 0 & 0 & 0 & 0 \\
  1 & 0 & 0 & 1 & 0 & 0 \\
  1 & 0 & 0 & 0 & 1 & 0 \\
  1 & 0 & 0 & 0 & 0 & 1 \\
  1 & 1 & 0 & 0 & 0 & 0 \\
  1 & 1 & 0 & 1 & 0 & 0 \\
  1 & 1 & 0 & 0 & 1 & 0 \\
  1 & 1 & 0 & 0 & 0 & 1 \\
  1 & 0 & 1 & 0 & 0 & 0 \\
  1 & 0 & 1 & 1 & 0 & 0 \\
  1 & 0 & 1 & 0 & 1 & 0 \\
  1 & 0 & 1 & 0 & 0 & 1 
\end{pmatrix}
,
\qquad
b
=
\begin{pmatrix}
  \alpha \\
  \beta_2 \\
  \beta_3 \\
  \gamma_2 \\
  \gamma_3 \\
  \gamma_4 \\
\end{pmatrix}
$$

In this setup, we are contrasting $i = j = 1$ against all other groups.

## Exercise 5

This can be modelled as `lm(deformation ~ 0 + load + I(load^2):alloy)`.

## Exercise 6

We have

$$
\begin{align}
  X^T \hat\mu 
  &=
  X^T X \hat\beta
  \\
  &=
  X^T X (X^T X)^{-1} X^T y
  \\
  &=
  X^T y
\end{align}
$$

from which we can conclude that the sum of residuals, $\sum_1^n (y_i - \hat\mu_i) = 0$, is zero.

## Exercise 7

Assuming homoscadisticity and independence of observations, we can show

$$ 
\newcommand{\tr}{\mathop{tr}}
\begin{align}
  \mathbb E (\vert r \vert^2)
  &=
  \tr \mathbb E \left(r r^T \right)
  \\
  &=
  \tr \mathbb E \left((I - A) y y^T (I - A)^T \right)
  \\
  &=
  \tr (I - A) \mathbb E  (y y^T) (I - A)
  \\
  &=
  \tr (I - A) \left( \mathbb V (y) + \mathbb E (y) \mathbb E (y)^T \right) (I - A)
  \\
  &=
  \tr (I - A) \sigma^2 + \tr \mathbb E (r) \mathbb E (r)^T
  \\
  &=
  (n - p) \sigma^2 + 0
\end{align}
$$

since $\tr A = \tr X (X^T X)^{-1} X^T = \tr X^T X (X^T X)^{-1} = \tr I$. 

## Exercise 8

Let's start by loadding the `MASS` library and fitting the full model.

```{r}
library(MASS)

m <- lm(
  loss ~ 1 + 
         tens + hard + 
         I(tens^2) + tens:hard + I(hard^2) + 
         I(tens^3) + I(tens^2):hard + tens:I(hard^2) + I(hard^3), 
  data = Rubber
)

mfull <- m
summary(m)
```

This has an AIC of `r AIC(m)`. There are no features that are significant ($\alpha = 0.05$) and we try removing the 'least significant' feature: `hard:I(tens^2)`.

```{r}
m <- lm(
  loss ~ 1 + 
         tens + hard + 
         I(tens^2) + tens:hard + I(hard^2) + 
         I(tens^3) + tens:I(hard^2) + I(hard^3), 
  data = Rubber
)
summary(m)
```

This has an AIC of `r AIC(m)`. The $R^2$ has increased and some features are now closer to significance. Let's again drop the least significant feature: `hard`.

```{r}
m <- lm(
  loss ~ 1 + 
         tens + 
         I(tens^2) + tens:hard + I(hard^2) + 
         I(tens^3) + tens:I(hard^2) + I(hard^3), 
  data = Rubber
)
summary(m)
```

The AIC is now `r AIC(m)`. The $R^2$ has not dropped and two features have become significant. Let's try dropping another feature: `I(hard^2)`.

```{r}
m <- lm(
  loss ~ 1 + 
         tens + 
         I(tens^2) + tens:hard +
         I(tens^3) + tens:I(hard^2) + I(hard^3), 
  data = Rubber
)
msig <- m
summary(m)
```

The AIC is `r AIC(m)`, about 1 higher than the previous model, which is consistent with the intepretation that the dropped feature is noise. There was still no drop in $R^2$ and now every feature is significant. Out of curiosity, we can drop one final feature: `tens`.

```{r}
m <- lm(
  loss ~ 1 + 
         I(tens^2) + tens:hard +
         I(tens^3) + tens:I(hard^2) + I(hard^3), 
  data = Rubber
)
summary(m)
```

The AIC `r AIC(m)` has increased several units, indicating that we have dropped a decent signal. The $R^2$ has also decreased, so we we keep `tens` in the model.

```{r}
step(mfull)
```

```{r}
maic <- lm(
  loss ~
    1 + 
    tens + hard + 
    I(tens^2) + tens:hard + I(hard^2) + 
    I(tens^3) + tens:I(hard^2) + I(hard^3), 
  data = Rubber
)

summary(maic)
```

The model chosen by `step` is more complex than that chosen by significance selection.

```{r}
resolution <- 100
hard <- seq(min(Rubber$hard), max(Rubber$hard), length.out = resolution)
tens <- seq(min(Rubber$tens), max(Rubber$tens), length.out = resolution)

z <- matrix(
  predict(msig, expand.grid(hard = hard, tens = tens)),
  nrow = length(hard)
)

contour(x = hard, y = tens, z = z)
```

```{r}
z <- matrix(
  predict(maic, expand.grid(hard = hard, tens = tens)),
  nrow = length(hard)
)

contour(x = hard, y = tens, z = z)
```

## Exercise 9

```{r}
m9 <- lm(breaks ~ tension + wool + tension:wool, data = warpbreaks)
summary(m9)
```

```{r}
anova(m9)
```

```{r}
interaction.plot(warpbreaks$tension, warpbreaks$wool, warpbreaks$breaks)
```

## Exercise 10

The scatterplot offers some evidence of a non-linear relationship.

```{r, echo=FALSE}
cars %>%
  ggplot(aes(x = speed, y = dist)) +
  geom_point(alpha = 0.4) +
  geom_smooth(se = FALSE) +
  geom_smooth(method = 'lm', colour = 'red', se = FALSE) +
  labs(
    x = 'Speed (mph)',
    y = 'Distance (feet)',
    title = 'Stopping distance as a function of speed',
    subtitle = 'Loess vs. linear smooths'
  )
```


```{r}
m10 <- lm(dist ~ 1 + speed + I(speed^2), data = cars)
summary(m10)
```

```{r}
step(m10)
```

```{r}
m10_aic <- lm(dist ~ 1 + I(speed^2), data = cars)
summary(m10_aic)
```

```{r}
m10_pval <- lm(dist ~ 1 + speed, data = cars)
summary(m10_pval)
```

Both methods of selection give models with similar performance with respect to $R^2$, AIC, and significance.
The diagnostic plots are also similar, with the AIC-selected model perhaps having the advantage. One should be careful with the interpretation of either model since distance must be positive but the linear model allows distance to be negative. This is particularly clear with `m10_pval`, which would predict a negative stopping distance for a stationary car. The positive intercept for `m10_aic` is not so problematic since we would expect the intercept to be positive for any positive speed due to reaction time. We interpret the constant of the `m10_aic` model as the distance travelled between receiving the stop signal and actually applying the breaks.

We could also deal with the zero-speed problem by removing the constant.

```{r}
m10_0 <- lm(dist ~ 0 + speed + I(speed^2), data = cars)
summary(m10_0)
```

Now both linear and quadrtic terms are significant, as expected from theoretical considerations. Moreover, `m10_0` has the best AIC score of the three.

```{r}
AIC(m10_pval, m10_aic, m10_0)
```

To calculate the reaction time in seconds, we convert distance from feet to miles and time from hours to seconds.

```{r}
c <- summary(m10_aic)$coefficients[1, 1]
cars %>%
  mutate(
    reaction_seconds = ((c / 5280) / speed) * (60 * 60)
  ) %>%
  summarise(mean(reaction_seconds))
```

```{r}
c <- summary(m10_0)$coefficients[1, 1]
cars %>%
  mutate(
    reaction_seconds = (c / 5280) * (60 * 60)
  ) %>%
  summarise(mean(reaction_seconds))
```

The above considerations indicate that we should exclude `m10` and `m10_pval` as appropriate models. The remaining models exclude either the contant term (`m10_0`) or the linear term (`m10_aic`) whilst keeping the quadratic term.

## Exercise 11

```{r}
cars11 <- model.matrix(dist ~ 1 + speed + I(speed^2), cars)

ols <- function(y, X) {
  d <- qr(X) # the decomposition
  r <- d$rank
  Rinv <- backsolve(qr.R(d), diag(1, nrow = r))
  return(Rinv %*% t(qr.Q(d)) %*% y)
}

b <- ols(cars$dist, cars11)
b
```

```{r}
m10_summary <- summary(m10)
m10_summary$coefficients[,1]
```

```{r}
residual_variance <- function(y, X) {
  b <- ols(y, X)
  r <- y - X %*% b
  p <- ncol(X)
  n <- length(y)
  return(sum(r^2) / (n - p))
}

v <- residual_variance(cars$dist, cars11)
se <- sqrt(v)
c(v, se)
```

```{r}
m10_summary$sigma
```