chapter_2.Rmd

---
title: "Chapter 2: Mixed Models"
author: "Brian Callander"
date: "February 27, 2018"
output: html_document
---

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

library(tidyverse)
library(broom)
library(lme4)
library(nlme)
```

## Exercise 1

a. The mixed model can be specified as $y_i = \alpha + b_{j[i]} + \epsilon_i$, where $b_{j[i]} \sim \mathcal N (0, \sigma_b^2)$ and $\epsilon_i \sim \mathcal N (0, \sigma^2)$.

b. The `anova` call tests the hypothesis that all $b_{j[i]} = 0$. In the mixed-model setup, this is equivalent to $\sigma_b^2 = 0$.

c. There is strong evidence for between-sow variation, i.e. that $\sigma_b^2 > 0$.

d. An unbiased estimator for between-sow variation is $\hat\sigma_b^2 = \hat\sigma_0^2 - \hat\sigma^2 / J = 0.6351067 - 0.2131667 / 6 = `r 0.6351067 - 0.2131667 / 6`$.

## Exercise 2


We have $a = \alpha + J^{-1} \sum_1^J c_j$ and $e_i = b_i + J^{-1} \sum_1^J \epsilon_{ij}$. The errors $e_i$ are mutually independent since the constituent random variables are. Since $\mathbb V (e_i) = \sigma_b^2 + J^{-1} \sigma^2$, we have 

$$
\newcommand{\rss}{\mathop{RSS}}
\hat \sigma_b^2 = \frac{\rss_0}{I - 1} - \frac{\hat\sigma^2}{J}
$$

where $\rss_0$ is the residual sum of squares for $\bar y_{i.} = a + e_i$.

The above argument is symmetric for $\hat\sigma_c^2$, which can be estimated by

$$
\hat \sigma_c^2 = \frac{\rss_0}{J - 1} - \frac{\hat\sigma^2}{I}
$$

where $\rss_0$ is now the residual sum of squares for $\bar y_{.j} = a' + e_j'$.

## Exercise 3

a. The final row of `anova(lm(y ~ 1 + x + a))` would test the hypothesis.

b. We could fit the following two models

```r
m1 <- lm(y ~ 1 + x + a, data = df)
dfa <- df %>% 
  group_by(a) %>%
  summarise(
    y = mean(y),
    x = mean(x)
  ) 
m0 <- lm(y ~ x, data = dfa)
```

then estimate $\beta$ as the coefficient of $x$ in `m0`, and $\sigma_a^2$ as described in the chapter.

## Exercise 4

a.

$$
\begin{pmatrix}
  y_{11} \\
  y_{12} \\
  y_{21} \\
  y_{22} \\
  y_{31} \\
  y_{32} \\
  y_{41} \\
  y_{42}
\end{pmatrix}
=
\begin{pmatrix}
  1 & x_{11} \\
  1 & x_{12} \\
  1 & x_{21} \\
  1 & x_{22} \\
  1 & x_{31} \\
  1 & x_{32} \\
  1 & x_{41} \\
  1 & x_{42}
\end{pmatrix}
\begin{pmatrix}
  \alpha \\
  \beta
\end{pmatrix}
+
\begin{pmatrix}
  1 & 0 & 0 & 0 \\
  1 & 0 & 0 & 0 \\
  0 & 1 & 0 & 0 \\
  0 & 1 & 0 & 0 \\
  0 & 0 & 1 & 0 \\
  0 & 0 & 1 & 0 \\
  0 & 0 & 0 & 1 \\
  0 & 0 & 0 & 1 
\end{pmatrix}
\begin{pmatrix}
  a_1 \\
  a_2 \\
  a_3 \\
  a_4
\end{pmatrix}
+
\begin{pmatrix}
  \epsilon_{11} \\
  \epsilon_{12} \\
  \epsilon_{21} \\
  \epsilon_{22} \\
  \epsilon_{31} \\
  \epsilon_{32} \\
  \epsilon_{41} \\
  \epsilon_{42}
\end{pmatrix}
$$


b. There are 2 treatment levels $\alpha$, 4 trees $b$ with 2 per treatment level, and 4 measurements per tree.

$$
\begin{pmatrix}
  y_{11} \\
  y_{12} \\
  y_{13} \\
  y_{14} \\
  y_{21} \\
  y_{22} \\
  y_{23} \\
  y_{24} \\
  y_{31} \\
  y_{32} \\
  y_{33} \\
  y_{34} \\
  y_{41} \\
  y_{42} \\
  y_{43} \\
  y_{44}
\end{pmatrix}
=
\begin{pmatrix}
  1 & 0 \\
  1 & 0 \\
  1 & 0 \\
  1 & 0 \\
  1 & 0 \\
  1 & 0 \\
  1 & 0 \\
  1 & 0 \\
  0 & 1 \\
  0 & 1 \\
  0 & 1 \\
  0 & 1 \\
  0 & 1 \\
  0 & 1 \\
  0 & 1 \\
  0 & 1
\end{pmatrix}
\begin{pmatrix}
  \alpha_0 \\
  \alpha_1
\end{pmatrix}
+
\begin{pmatrix}
  1 & 0 & 0 & 0 \\
  1 & 0 & 0 & 0 \\
  1 & 0 & 0 & 0 \\
  1 & 0 & 0 & 0 \\
  0 & 1 & 0 & 0 \\
  0 & 1 & 0 & 0 \\
  0 & 1 & 0 & 0 \\
  0 & 1 & 0 & 0 \\
  0 & 0 & 1 & 0 \\
  0 & 0 & 1 & 0 \\
  0 & 0 & 1 & 0 \\
  0 & 0 & 1 & 0 \\
  0 & 0 & 0 & 1 \\
  0 & 0 & 0 & 1 \\
  0 & 0 & 0 & 1 \\
  0 & 0 & 0 & 1 
\end{pmatrix}
\begin{pmatrix}
  b_0 \\
  b_1 \\
  b_2 \\
  b_3
\end{pmatrix}
+
\begin{pmatrix}
  \epsilon_{11} \\
  \epsilon_{12} \\
  \epsilon_{13} \\
  \epsilon_{14} \\
  \epsilon_{21} \\
  \epsilon_{22} \\
  \epsilon_{23} \\
  \epsilon_{24} \\
  \epsilon_{31} \\
  \epsilon_{32} \\
  \epsilon_{33} \\
  \epsilon_{34} \\
  \epsilon_{41} \\
  \epsilon_{42} \\
  \epsilon_{43} \\
  \epsilon_{44}
\end{pmatrix}
$$



c.

$$
\begin{pmatrix}
  y_{111} \\
  y_{112} \\
  y_{113} \\
  y_{121} \\
  y_{122} \\
  y_{123} \\
  y_{131} \\
  y_{132} \\
  y_{133} \\
  y_{212} \\
  y_{211} \\
  y_{212} \\
  y_{221} \\
  y_{222} \\
  y_{221} \\
  y_{232} \\
  y_{231} \\
  y_{232}
\end{pmatrix}
=
\begin{pmatrix}
  1 & 0 \\
  1 & 0 \\
  1 & 0 \\
  1 & 0 \\
  1 & 0 \\
  1 & 0 \\
  1 & 0 \\
  1 & 0 \\
  1 & 0 \\
  1 & 1 \\
  1 & 1 \\
  1 & 1 \\
  1 & 1 \\
  1 & 1 \\
  1 & 1 \\
  1 & 1 \\
  1 & 1 \\
  1 & 1
\end{pmatrix}
\begin{pmatrix}
  \mu \\
  \alpha_2
\end{pmatrix}
+
\begin{pmatrix}
  1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\
  1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\
  1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\
  0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\
  0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\
  0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\
  0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 \\
  0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 \\
  0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 \\
  1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\
  1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\
  1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\
  0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\
  0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\
  0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\
  0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 \\
  0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 \\
  0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 
\end{pmatrix}
\begin{pmatrix}
  b_1 \\
  b_2 \\
  b_3 \\
  c_{11} \\
  c_{12} \\
  c_{23} \\
  c_{21} \\
  c_{32} \\
  c_{33}
\end{pmatrix}
+
\begin{pmatrix}
  \epsilon_{111} \\
  \epsilon_{112} \\
  \epsilon_{113} \\
  \epsilon_{121} \\
  \epsilon_{122} \\
  \epsilon_{123} \\
  \epsilon_{131} \\
  \epsilon_{132} \\
  \epsilon_{133} \\
  \epsilon_{212} \\
  \epsilon_{211} \\
  \epsilon_{212} \\
  \epsilon_{221} \\
  \epsilon_{222} \\
  \epsilon_{221} \\
  \epsilon_{232} \\
  \epsilon_{231} \\
  \epsilon_{232}
\end{pmatrix}
$$

## Exercise 5

a. With $\mu := \mathbb E (X + Z) = \mu_X + \mu_Z$ and 

$$
\newcommand{\cov}{\mathop{Cov}}
\begin{align}
  \cov (X + Z)
  &=
  \cov X + \cov Z + \cov (X, Z) \cov (X, Z)^T
\end{align}
$$
  
it is sufficient to show that $\cov (X, Z) = 0$. Since $X, Z$ are independent, it follows that $\mathbb E (XZ^T) = \mathbb E X \mathbb E Z^T = \mu_X \mu_Z^T$. Therefore,

$$
\begin{align}
  \cov (X, Z)
  &=
  \mathbb E \left( (X - \mu_X) (Z - \mu_Z)^T \right)
  \\
  &=
  \mathbb E \left( XZ^T + \mu_X \mu_Z^T - X\mu_Z^T - \mu_X Z^T \right)
  \\
  &=
  \mu_X \mu_Z^T + \mu_X \mu_Z^T - \mu_X \mu_Z^T - \mu_X \mu_Z^T 
  \\
  &=
  0
\end{align}
$$

b. The matrix form of the model is

$$
\begin{pmatrix}
  y_{11} \\
  y_{12} \\
  y_{13} \\
  y_{21} \\
  y_{22} \\
  y_{23} \\
  y_{31} \\
  y_{32} \\
  y_{33}
\end{pmatrix}
=
\begin{pmatrix}
  1 & x_{11} \\
  1 & x_{12} \\
  1 & x_{13} \\
  1 & x_{21} \\
  1 & x_{22} \\
  1 & x_{23} \\
  1 & x_{31} \\
  1 & x_{32} \\
  1 & x_{33}
\end{pmatrix}
\begin{pmatrix}
  \alpha \\
  \beta
\end{pmatrix}
+
\begin{pmatrix}
  1 & 0 & 0 \\
  1 & 0 & 0 \\
  1 & 0 & 0 \\
  0 & 1 & 0 \\
  0 & 1 & 0 \\
  0 & 1 & 0 \\
  0 & 0 & 1 \\
  0 & 0 & 1 \\
  0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
  b_1 \\
  b_2 \\
  b_3
\end{pmatrix}
+
\begin{pmatrix}
  \epsilon_{11} \\
  \epsilon_{12} \\
  \epsilon_{13} \\
  \epsilon_{21} \\
  \epsilon_{22} \\
  \epsilon_{23} \\
  \epsilon_{31} \\
  \epsilon_{32} \\
  \epsilon_{33}
\end{pmatrix}
$$

$$
\begin{align}
  \cov y
  &=
  \cov (Zb + \epsilon)
  \\
  &=
  \cov (Zb) + \cov \epsilon
  \\
  &=
  Z\cov (b) Z^T + \sigma^2
  \\
  &=
  \sigma_b^2 ZZ^T + \sigma^2
\end{align}
$$


$$
ZZ^T
=
\begin{pmatrix}
  1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
  1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
  1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
  0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 \\
  0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 \\
  0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 \\
  0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 \\
  0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 \\
  0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 
\end{pmatrix}
$$

$$
\cov y
=
\begin{pmatrix}
  \sigma_b^2 + \sigma & \sigma_b^2 & \sigma_b^2   & 0 & 0 & 0 & 0 & 0 & 0 \\
  \sigma_b^2 & \sigma_b^2 + \sigma^2 & \sigma_b^2 & 0 & 0 & 0 & 0 & 0 & 0 \\
  \sigma_b^2 & \sigma_b^2 & \sigma_b^2 + \sigma^2 & 0 & 0 & 0 & 0 & 0 & 0 \\
  0 & 0 & 0 & \sigma_b^2 + \sigma & \sigma_b^2 & \sigma_b^2   & 0 & 0 & 0 \\
  0 & 0 & 0 & \sigma_b^2 & \sigma_b^2 + \sigma^2 & \sigma_b^2 & 0 & 0 & 0 \\
  0 & 0 & 0 & \sigma_b^2 & \sigma_b^2 & \sigma_b^2 + \sigma^2 & 0 & 0 & 0 \\
  0 & 0 & 0 & 0 & 0 & 0 & \sigma_b^2 + \sigma & \sigma_b^2 & \sigma_b^2   \\
  0 & 0 & 0 & 0 & 0 & 0 & \sigma_b^2 & \sigma_b^2 + \sigma^2 & \sigma_b^2 \\
  0 & 0 & 0 & 0 & 0 & 0 & \sigma_b^2 & \sigma_b^2 & \sigma_b^2 + \sigma^2 
\end{pmatrix}
$$


## Exercise 6

The `Source` should be modelled as a fixed effect since this is the main factor of interest. Since `Site` is suspected to have a systemic effect, this can also be moddeled as a fixed effect. We model `Lot` and `Wafer` as nested random effects to account for additional variation.

### EDA

```{r}
Oxide %>%
  ggplot(aes(Thickness, fill = Source)) +
  geom_density(alpha = 0.3)
#lm(Thickness ~ Source, data = Oxide)
```


```{r}
Oxide %>%
  ggplot(aes(Thickness, fill = Site)) +
  geom_density(alpha = 0.3)
```

### No random effects

For comparison with the mixed-effect models below, we check here what we can conclude without any random effects. The most basic linear model suggests that there is a significant `Source` effect, but there is insufficient evidence to suggest any systematic `Site` effect.

```{r}
m00 <- lm(Thickness ~ 0 + Source, data = Oxide)
m01 <- lm(Thickness ~ 0 + Source + Site, data = Oxide)
anova(m01)
```

```{r}
summary(m00)
```

There seems to be an increase in variance from site 1 to site 2, but there's otherwise no indication of problems with the model assumptions. 

```{r}
plot(m00)
```

### Mixed-effects

The full model with random effects is as follows.

```{r}
ctrl <- lmeControl(niterEM = 500, msMaxIter = 100)
options(contrasts = c("contr.treatment", "contr.treatment"))

m1 <- lme(Thickness ~ Source + Site, data = Oxide, ~ 1 | Lot/Wafer, control = ctrl)
summary(m1)
```

Neither MLE nor AIC indicate that `Site` has a systemic effect, so we can drop it from our model.

```{r}
m2 <- lme(Thickness ~ Source, data = Oxide, ~ 1 | Lot/Wafer, control = ctrl)
anova(m2, m1)
```

There is no strong evidence that any model assumptions have been violated.

```{r}
plot(m2)
plot(m2, Thickness ~ resid(.))
qqnorm(m2, ~ resid(.))
```

The random effects also look normal.

```{r}
qqnorm(m2, ~ranef(., level=1))
qqnorm(m2, ~ranef(., level=2))
```

There is insufficient evidence to reject the hypothesis that there is a systematic difference between the two sources. Any difference that might exist would by masked by the large variance between lots.

```{r}
summary(m2)
```


We can compare the model `m2` with the equivalent model without any random effects. There is strong evidence of both a `Lot` and `Wafer` random effect.

```{r}
m0 <- lme(Thickness ~ Source, data = Oxide)
m3 <- lme(Thickness ~ Source, data = Oxide, ~ 1 | Lot)
anova(m0, m3, m2)
```

I am unsure why `m0` would have 9 degrees of freedom. Since it uses only one binary factor and a constant, it should have 2 degrees of freedom.

## Exercise 9

Both `Method` and `Physique` are fixed effects since we could expect these to have systemic effects on `rounds`. The factor `Team` is a random effect because repeating the experiment with different teams could lead to different results.

A suitable starting model includes the interactions of the fixed effects. The anova shows that `Method` has a significant effect, but neither `Physique` nor the interaction term are significant.

```{r}
m <- lme(rounds ~ Method * Physique, data = Gun, ~1 | Team)
anova(m)
```

The `Physique` factor is still insignifiant after removing the interaction term.

```{r}
m1 <- lme(rounds ~ Method + Physique, data = Gun, ~1 | Team)
anova(m1, m)
```

Removing `Physique` altogether yields the best results.

```{r}
m2 <- lme(rounds ~ Method, data = Gun, ~1 | Team)
anova(m2, m1)
```

The diagnistics are mostly satisfactory. There is some evidence that `M2` has higher variance than `M1` and the random effects are not quite normal.

```{r}
diagnostic <- augment(m2) %>%
  mutate(
    .resid_student = residuals(m2, type = 'pearson', scaled = TRUE),
    .random = .fitted - .fixed
  )

diagnostic %>%
  ggplot(aes(.fitted, .resid, colour = Method)) +
  geom_point() 
  #geom_smooth(se = FALSE, colour = "black")

diagnostic %>%
  ggplot(aes(y = sqrt(abs(.resid_student)), x = .fitted, colour = Method)) +
  geom_point() +
  geom_smooth(se = FALSE, colour = "black", formula = y ~ x)

diagnostic %>%
  ggplot(aes(sample = .resid_student)) +
  geom_point(stat = 'qq') +
  geom_abline()

ranef(m2) %>%
  ggplot(aes(sample = `(Intercept)`)) +
  geom_point(stat = 'qq') +
  geom_abline()
```

The summary tells us that between-team variance is fairly small, smaller than within-team variance, which in turn is a smaller than the estimated effect size.

```{r}
summary(m2)
```

```{r}
intervals(m2)
```