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

# Background

This week we will see how to fit linear models in **R**, and estimate parameter uncertainty and confidence intervals. The data come from a paper by Johnson & Raven (1973, *Science* 179:893-895) and consist of the diversity of plants species across 30 islands in the Galapagos Archipelago. The data set also contains measurements of 5 possible explanatory variables:

* `Area`: the area of the island (km<sup>2</sup>)  
* `Elevation`: the highest elevation of the island (m)  
* `Nearest`: the distance from the nearest island (km)  
* `Scruz`: the distance from Santa Cruz Island (km)  
* `Adjacent`: the area of the adjacent island (km<sup>2</sup>)

The `gala` dataset is included in the **faraway** package. Let's load the package and inspect the data.

```{r load_gala_dataset}
library(faraway)
head(gala)
```

The second column `Endemics` is the number of plant species that are endemic to each island, but we will just use the total number of plant species recorded in `Species`.


# A simple model for plant diversity

We seek a linear model that would help us explain the diversity (number) of plant species as a function of some predictors, which are sometimes referred to as *covariates*. From the theory of island biogeography, it would seem reasonable that plant diversity is related to the area of an island, such that

<center>species diversity = $f$(area)</center><br>

More formally, we can write our explicit linear regression model as

$$
\text{species}_i = \beta_0 + \beta_1 \text{area}_i + e_i
$$

where $i$ indicates an island.

Let's plot the repsonse and predictor to see what the relationship looks like.

```{r plot_gala, echo=FALSE, fig.height=4.5, fig.width=5, fig.align="center"}
par(mai = c(0.9, 0.9, 0.3, 0.1), omi = rep(0, 4))
plot(gala$Area, gala$Species, pch = 16, col = "dark green",
     ylab = "Number of species", xlab = expression(paste("Area of island (", km^2, ")")))
```

## Estimating $\hat{\boldsymbol{\beta}}$

We saw in lecture that we can estimate the regression parameters as

$$
\hat{\boldsymbol{\beta}} = (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{y}.
$$

Let's begin by creating the matrix of response values $\mathbf{y}$, which should have dimensions $n \times 1$.

```{r create_y_matrix}
nn <- nrow(gala)
yy <- matrix(gala$Species, nrow = nn, ncol = 1)
head(yy)
```

Next let's create the matrix of predictors $\mathbf{X}$. There is only one predictor in the model (`Area`), so our matrix should consist of a column of 1's and a column with the area of each island.

```{r create_design_matrix}
Intercept <- rep(1, nn)
Area <- gala$Area
XX <- cbind(Intercept, Area)
head(XX)
```

Matrix multiplication in **R** uses the `%*%` operator. You can transpose a matrix with `t()` and invert a matrix with `solve()`. So our estimate of the regression parameters in \hat{\boldsymbol{\beta}} is

```{r est_beta}
beta_hat <- solve(t(XX) %*% XX) %*% t(XX) %*% yy
beta_hat
```

We can interpret these coefficients to mean that an island with an area of 0 km<sup>2</sup> would be expected to have ~`r round(beta_hat[1], 0)` species of plants(!), and that you would add ~`r round(beta_hat[2]*100, 0)` species per 100 km<sup>2</sup> of island.

## Estimating $\hat{\mathbf{y}}$

Recall that the fitted values are given by

$$
\begin{align}
\hat{\mathbf{y}} &= \mathbf{X} \hat{\boldsymbol{\beta}} \\
  &= \mathbf{X} \left( (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{y} \right) \\
  &= \underbrace{\mathbf{X} (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^{\top}}_{\mathbf{H}} \mathbf{y} \\
  &= \mathbf{H} \mathbf{y}.
\end{align}
$$

So we can calculate the hat matrix $\mathbf{H}$ and solve for $\hat{\mathbf{y}}$

```{r get_hatmat}
HH <- XX %*% solve(t(XX) %*% XX) %*% t(XX)
y_hat <- HH %*% yy
```

Let's add the predicted values to our plot of the data.

```{r plot_gala_fit, fig.height=4.5, fig.width=5, fig.align="center"}
par(mai = c(0.9, 0.9, 0.3, 0.1), omi = rep(0, 4))
## plot the data
plot(gala$Area, gala$Species, pch = 16, col = "dark green",
     ylab = "Number of plant species", xlab = expression(paste("Area of island (", km^2, ")")))
## add the fitted values (line)
abline(a = beta_hat[1], b = beta_hat[2])
```

This is not a great fit, but it's a decent start. We'll expand upon our model later.

## A shortcut to model fitting

As we saw briefly in lecture, it's easy to fit a linear regression model using the function `lm()`, which takes arguments like

`lm(respone ~ predictor(s), data_to_use)`

To do so for our model above, we have

```{r fit_with_lm}
simple_model <- lm(Species ~ Area, gala)
simple_model
```

Note that these are the same estimates that we obtained earlier

We can obtain the estimates of the fitted values $\hat{\mathbf{y}}$ via one of two functions in **R**: `fitted()` and `predict()`. We'll see later that `predict()` has additional functionality that we can use.

```{r get_fitted_y}
## via fitted
y_hat_f <- fitted(simple_model)
## via predict
y_hat_p <- predict(simple_model, type = "response")
## compare these to each other
all.equal(y_hat_f, y_hat_p)
```

## Goodness-of-fit

We saw in lecture that we can use the *coefficient of determination* $(R^2)$ to measure how well a model fits the data. The formula is

$$
R^2 = 1 - \frac{SSE}{SSTO}
$$

We can easly compute these quantities from our model fit and data.

```{r compute_R2}
## SSE
resids <- yy - y_hat
SSE <- sum(resids^2) # = t(resids) %*% resids
## SSTO
SSTO <- sum((yy - mean(yy))^2)
## R^2
1 - SSE / SSTO
```

Let's compare this to the value in the ANOVA table that **R** produces

```{r check_R2}
summary(simple_model)
```

The estimates are indeed the same.

### Adjusted $R^2$

The ANOVA table above also reports the `Adjusted R-squared`, which we haven't discussed yet. It turns out that one can always increase the $SSR$ of a model by including more predictors, even if they are completed unrelated to the response. Recall our objective in least squares estimation is

$$
\min \sum SSE = \min \sum_i (y_i - X_i \beta).
$$

Increasing the dimension of $\mathbf{X}$ by adding more predictors (columns) necessitates an increased space within which to minimize this objective function, and hence you can get a better fit o the data. We can demonstrate this with a simple example. Suppose we added another predictor to our simple model that was nothing more than a vector of Gaussian white noise, which is unrelated to the response.

```{r demo_adj_R2}
set.seed(514)
## generate a vector of Gaussian white noise
WN <- rnorm(nn)
## add this to our Galapagos data frame
gala$WN <- WN 
## fit a model with Area & WN
summary(lm(Species ~ Area + WN, gala))
```

Notice that our $R^2$ value has indeed increased.

To account for this, we can use an *adjusted* $R^2$, which corrects for these additional degrees of freedom.

$$
\bar{R}^2 = 1 - \frac{SSE ~ / ~ df_{error}}{SSTO ~ / ~ df_{total}}
$$

Let's calculate $\bar{R}^2$ for our simple model and confirm that it matches the value in the ANOVA table for the model with only `Area` as a predictor.

```{r adjusted_R2}
1 - (SSE / (nn - 2)) / (SSTO / (nn - 1))
```


# A better model for plant diversity

Let's expand our model from above to include some more predictors. In particular, we'll add in the effects of an island's elevation and the distance to the nearest island, such that

<center>plant diversity = $f$(area, elevation, distance to nearest island)</center><br>  

More formally, we can write our explicit linear regression model as

$$
\text{species}_i = \beta_0 + \beta_1 \text{area}_i + \beta_2 \text{elevation}_i + \beta_3 \text{nearest}_i + e_i
$$

where $i$ indicates an island.

## Fit the model

We'll use `lm()` to estimate this model.

```{r fit_multiple_regr}
## larger model
full_mod <- lm(Species ~ Area + Elevation + Nearest, gala)
full_mod
```

Notice that our estimates for the intercept $\beta_0$ and the effect of island area $\beta_1$ changed when we accounted for additional predictors in the model.

## Hypothesis tests to compare models

We saw in lecture a number of different hypothesis tests we might wish to undertake. Let's run through examples of those here.

### Test of *all* predictors in a model

Suppose we wanted to test whether this collection of 3 predictors in our model were better than simply estimating the data by their mean.

$$
\Theta: \mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \mathbf{e} \\
\theta: \mathbf{y} = \boldsymbol{\mu} + \mathbf{e} \\
$$

We write the null hypothesis as

$$
H_0: \beta_1 = \beta_2 = \dots = \beta_k = 0
$$

and we would reject $H_0$ if $F > F^{(\alpha)}_{k_{\Theta} - k_{\theta}, n - k_{\Theta}}$

Recall that we wil base this test on an $F$-distribution, such that

$$
SSE_{\Theta} = \left( \mathbf{y} - \mathbf{X} \boldsymbol{\beta} \right)^{\top} \left( \mathbf{y} - \mathbf{X} \boldsymbol{\beta} \right) = \mathbf{e}^{\top} \mathbf{e} = SSE \\
SSE_{\theta} = \left( \mathbf{y} - \bar{y} \right)^{\top} \left( \mathbf{y} - \bar{y} \right) =  SSTO \\
\Downarrow \\
F = \frac{ \left( SSTO - SSE \right)  / (k - 1) } { SSE  / (n - k)}
$$

#### $F$-test by hand

Calculating the sums-of-squares in **R** is straightforward with matric operations.

```{r F_test_by_hand}
## get matrix of predictors
XX <- model.matrix(full_mod)
## estimate beta
beta_hat <- solve(t(XX) %*% XX) %*% t(XX) %*% yy
## total sum of squares
SSE <- t(yy - XX %*% beta_hat) %*% (yy - XX %*% beta_hat)
## error sum of squares
SSTO <- t(yy - mean(yy)) %*% (yy - mean(yy))
## F statistic
(F_stat <- ((SSTO - SSE) / (4 - 1)) / (SSE / (nn - 4)))
## F test
pf(F_stat, 4-1, nn-4, lower.tail = F)
```

This $F$-statistic is quite large and the $p$-value is very small, so we would reject the null hypothesis that we would be justified in dropping the 3 predictors from this model in favor of a mean-only model.

#### $F$-test with built-in functions

We can use the `anova()` function in **R** to conduct this $F$-test

```{r gala_full, echo = TRUE}
## null model; the '1' indicates an intercept-only model
null_mod <- lm(Species ~ 1, gala)
## use `anova('simple', 'complex')` to get the F-test results
anova(null_mod, full_mod)
```

We can see that we get the same $F$-statistic and $p$-value as earlier


### Testing one predictor

We might ask whether any one predictor could be dropped from our model. For example, can `nearest` be dropped from this model? One option is to fit this reduced model and compare to it to the full model via our $F$-test with $H_0: \beta_3 = 0$

$$
\begin{aligned}
\Theta: \text{species}_i &= \beta_0 + \beta_1 \text{area}_i + \beta_2 \text{elevation}_i + \beta_3 \text{nearest}_i + e_i \\
~ \\ 
\theta: \text{species}_i &= \beta_0 + \beta_1 \text{area}_i + \beta_2 \text{elevation}_i + e_i
\end{aligned}
$$

Here's our test in **R**.

```{r gala_one_predictor_F, echo = TRUE}
## full model as before
full_mod <- lm(Species ~ Area + Elevation + Nearest, gala)
## reduced model without `nearest`
reduced_mod <- lm(Species ~ Area + Elevation, gala)
## use `anova('reduced', 'full')` to get the F-test results
anova(reduced_mod, full_mod)
```

The $F$-statistic is small and the $p$-value is large so we cannot reject $H_0$, and hence we conclude that we would be justified in dropping `Nearest` from our model.

#### Using a $t$-test

Another option for testing the signficance of one predictor is to estimate a $t$-statistic as

$$
t_i = \frac{\beta_i}{\text{SE} \left( \beta_i \right)}
$$

and compare it to a $t$-distribution with $n - k$ degrees of freedom. We can get the results of the $t$-test with `summary()` from base **R** or `sumary()` from the **faraway** pkg, which is much less verbose.

```{r gala_one_predictor_t, echo = TRUE}
## via `base::summary()`
summary(full_mod)
## via `faraway::sumary()`
sumary(full_mod)
```

Note that we get $t^2 = F$ and the same $p$-value as before


### Testing 2+ predictors

Sometimes we might want to know whether we can drop 2+ predictors from a model. For example, could we drop both $\text{area}$ and $\text{nearest}$ from our full model?

$$
\begin{aligned}
\Theta: \text{species}_i &= \beta_0 + \beta_1 \text{area}_i + \beta_2 \text{elevation}_i + \beta_3 \text{nearest}_i + e_i \\
~ \\ 
\theta: \text{species}_i &= \beta_0 + \beta_2 \text{elevation}_i + e_i
\end{aligned}
$$

$H_0 : \beta_1 = \beta_3 = 0$

Again we can conduct this test by fitting both models and calculating $F$.

```{r gala_one_predictor_F2, echo = TRUE}
## reduced model without `area` and `nearest`
reduced_mod <- lm(Species ~ Elevation, gala)
## use `anova('reduced', 'full')` to get the F-test results
anova(reduced_mod, full_mod)
```

The $F$-statistic is rather small and the $p$-value is big so we cannot reject $H_0$, suggesting we would be justified in dropping area and elevation from the model.


#### Multiple $t$-tests versus one $F$-test

Could we instead evaluate our null hypothesis by examining our table of $t$-statistics and associated $p$-values?

```{r gala_one_predictor_t2, echo = TRUE}
## t-statistics from our full model
sumary(full_mod)
```

This creates the problem of having to make our decision based on 2 $p$-values rather than one

* Each $p$-value corresponds to a separate $t$-test with the other predictor in the model

* Therefore we must use the single $F$-test 


### Testing a subspace

Some tests cannot be expressed in terms of the inclusion or exclusion of predictors.

Consider a test of whether the areas of the current and adjacent island could be added together and used in place of the two separate predictors

$$
\text{species}_i = \beta_0 + \beta_1 \text{area}_i + \beta_2 \text{adjacent}_i + \dots + e_i \\
~ \\
\text{species}_i = \beta_0 + \beta_1 \text{(area + adjacent)}_i + \dots + e_i
$$

$H_0 : \beta_{\text{area}} = \beta_{\text{adjacent}}$

We can do this in **R** by using the `I()` function, which tells **R** to treat things inside the `()` as an operation.

```{r gala_subspace, echo = TRUE}
## full model as before
full_mod <- lm(Species ~ Area + Adjacent + Elevation + Nearest, gala)
## reduced model without `elevation + nearest`
comb_mod <- lm(Species ~ I(Area + Adjacent) + Elevation + Nearest, gala)
## use `anova('combined', 'full')` to get the F-test results
anova(comb_mod, full_mod)
```

The $F$-statistic is large and the $p$-value is reasonably small so we reject $H_0$.

<br>

What if we wanted to test whether a predictor had a specific (non-zero) value? For example, is there a 1:1 relationship between $\text{species}$ and $\text{elevation}$ after controlling for the other predictors?

$$
\text{species}_i = \beta_0 + \beta_1 \text{area}_i + \underline{1} \text{elevation}_i + \beta_3 \text{nearest}_i + e_i
$$

$H_0 : \beta_2 = 1$

We can do this in **R** by using the `offset()` function, which allows us to set a fixed beta.

```{r gala_subspace_elev, echo = TRUE}
## full model
full_mod <- lm(Species ~ Area + Elevation + Nearest, gala)
## model with effect of `elevation` = 1
fixed_mod <- lm(Species ~ Area + offset(1 * Elevation) + Nearest, gala)
## use `anova('comb', 'full')` to get the F-test results
anova(fixed_mod, full_mod)
```

The $F$-statistic is very large and the $p$-value is tiny so we reject $H_0$.


#### Using a $t$-test

We can also modify our $t$-test from before and use it for our comparison by including the hypothesized $\hat{\beta}_{H_0}$ as an offset

$$
t_i = \frac{(\hat{\beta_i} - \beta_{H_0})}{\text{SE} \left( \hat{\beta}_i \right)}
$$

We can get $\hat{\beta}_{\text{elevation}}$ and $\text{SE}(\hat{\beta}_{\text{elevation}})$ from our model summary

```{r gala_subspace_elev_2, echo = TRUE}
## full model
sumary(full_mod)
```

This leaves us with $\hat{\beta}_{\text{elevation}} \approx 0.17$ and $\text{SE}(\hat{\beta}_{\text{elevation}}) \approx 0.055$. 

We can calculate the $t$-test by hand and compare the results to the above $F$-test

```{r t-test_subspace, echo = TRUE, results = "hold"}
## t statistic
(t_value <- (0.171336 - 1) / 0.054519)
## p-value = t_alpha * Pr(t_value, df); `pt()` is the pdf for a t-dist
(p_value <- 1.96 * pt(t_value, 26))
## verify t^2 = F
all.equal(t_value^2, anova(fixed_mod, full_mod)$F[2], tolerance = 0.0001)
```

The $p$-value is the same and $t^2 = F$ as before


## Confidence intervals for $\beta$

We can also use confidence intervals (CI's) to express uncertainty in $\hat{\beta}_i$, which take the form

$$
100(1 - \alpha)\% ~ \text{CI}: \hat{\beta}_{i} \pm t_{n-p}^{(\alpha / 2)} \operatorname{SE}(\hat{\beta})
$$

where here $\alpha$ is our *predetermined* Type-I error rate. To do so, we need the information from the model summary.

```{r gala_CI, echo = TRUE}
## t-statistics from our full model
sumary(full_mod)
```

We can construct a 95% CI on $\hat{\beta}_{\text{area}}$ using the 2.5% and 97.5% percentiles of the $t$-distribution with $df = 30 - 4 = 26$.

```{r gala_CI_2, echo = TRUE}
## critical value for the t-dist
## `qt()` is the quantile function for the t-dist; `p` is the (1-alpha/2) value 
t_crit <- qt(p = 0.975, df = 30-4)
## 95% CI
CI95_beta <- 0.019085 + c(-1,1) * t_crit * 0.026764
round(CI95_beta, 3)
```


### Confidence intervals for all $\beta$

We can quickly get the 95% CI for all the $\beta_i$ using the function `confint()`.

```{r gala_CI_all, echo = TRUE}
## all of the 95% CI's
confint(full_mod)
```


## Bootstrap confidence intervals

The above $F$- and $t$-based CI's depend on the assumption of normality. Recall that the bootstrap method provides a way to construct CI's without this assumption

### Bootstrap procedure

Step 1: Fit your model to the data

```{r}
## we already fit `full_mod` above
```

Step 2: Calculate $\mathbf{e} = \mathbf{y} - \mathbf{X} \hat{\boldsymbol{\beta}}$ 

```{r bootstrap_CI_1-2, echo = TRUE}
## residuals from our full model
resids <- residuals(full_mod)
```

Step 3a: Generate $\mathbf{e}^*$ by sampling *with replacement* from $\mathbf{e}$  
Step 3b: Calculate $\mathbf{y}^* = \mathbf{X} \hat{\boldsymbol{\beta}} + \mathbf{e}^*$  
Step 3c: Estimate $\hat{\boldsymbol{\beta}}^*$ from $\mathbf{X}$ & $\mathbf{y}^*$  

```{r bootstrap_CI_3, echo = TRUE}
## number of boostrap samples
nb <- 1000
## empty matrix for beta estimates
beta_est <- matrix(NA, nb, 4)
## fitted values from our full model = X*beta
Xbeta <- fitted(full_mod)
## sample many times
for(i in 1:nb) {
  ## 3a: sample w/ replacement from e
  e_star <- sample(resids, rep = TRUE)
  ## 3b: calculate y_star
  y_star <- Xbeta + e_star
  ## 3c: re-estimate beta_star from X & y_star
  beta_star <- update(full_mod, y_star ~ .)
  ## save estimated betas
  beta_est[i,] <- coef(beta_star)
}
```

Step 4: Select the $\tfrac{\alpha}{2}$ and $(1 - \tfrac{\alpha}{2})$ percentiles from the saved $\hat{\boldsymbol{\beta}}^*$

```{r bootstrap_CI_4, echo = TRUE}
## extract 2.5% and 97.5% values
CI95 <- apply(beta_est, 2, quantile, c(0.025, 0.975))
colnames(CI95) <- c("Intercept", colnames(gala[,3:5]))
t(round(CI95, 3))
```


## Confidence interval for new predictions

Given a fitted model $\mathbf{y} = \mathbf{X} \hat{\boldsymbol{\beta}} + \mathbf{e}$, we might want to know the uncertainty around a new estimate $\mathbf{y}^*$ given some new predictor $\mathbf{X}^*$. Recall that this CI is given by 

$$
\hat{\mathbf{y}}^* \pm ~ t^{(\alpha / 2)}_{df} \sigma \sqrt{ {\mathbf{X}^*}^{\top} (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^* } 
$$


### CI for the mean response

Think back to our simple model for plant diversity as a function of island area, such that

$$
\text{species}_i = \beta_0 + \beta_1 \text{area}_i + e_i
$$

Suppose we wanted to estimate the uncertainty in the *average* diversity for an island with an area of 2000 km<sup>2</sup>. Recall that our CI on the mean response is given by


and here $\mathbf{X}^* = [1 ~~ 2000]$ (ie, a $1 \times 2$ row matrix).


#### CI by hand

We can calculate this by hand using matrix operations in **R**.

```{r CI_by_hand, warning=FALSE, message=FALSE}
## matrix of predictors
XX <- model.matrix(simple_model)
## new X; vector for now
X_star <- c(Intercept = 1, Area = 2000)
## inside sqrt
inner_X <- t(X_star) %*% solve(t(XX) %*% XX) %*% X_star
## critical t-value
t_crit <- qt(0.975, df = nn-2)
## estimated SD
sigma <- summary(simple_model)$sigma
## predicted y
y_star <- sum(X_star * coef(simple_model))
## 95% CI
y_star + c(-1,1) * t_crit * sigma * sqrt(inner_X)
```

#### CI with `predict()`

It's much easier to use the `predict()` function to give us a CI.

```{r CI_with_predict}
predict(simple_model, new = data.frame(t(X_star)),
        level = 0.95, interval = "confidence")
```

The estimated lower and upper bounds are the same as those we obtained by hand.

### Prediction interval for a new response

We saw in lecture that we can also estimate the uncertainty for a *specific* prediction rather than the average response. Let's now imagine we wanted to predict how many plant species we might find on a hypothetical island that might one day emerge from an underwater volcano and is ultimately 2000 km<sup>2</sup>. Recall that a PI is given by 

$$
\hat{\mathbf{y}}^* \pm ~ t^{(\alpha / 2)}_{df} \sigma \sqrt{1 + {\mathbf{X}^*}^{\top} (\mathbf{X}^{\top} \mathbf{X})^{-1} \mathbf{X}^* } 
$$

To calculate the PI by hand simply requires a different calculation for `inner_X` above.

```{r PI_by_hand, warning=FALSE, message=FALSE}
## new X_star
X_star <- c(Intercept = 1, Area = 2000)
## inside sqrt
inner_X <- 1 + t(X_star) %*% solve(t(XX) %*% XX) %*% X_star
## 95% CI
y_star + c(-1,1) * t_crit * sigma * sqrt(inner_X)
```

Notice how much wider this interval is than the CI we estimated above. We could instead use `predict()` as we did above for the CI by switching the `interval` argument to `"prediction"`.

```{r PI_with_predict}
predict(simple_model, new = data.frame(t(X_star)),
        level = 0.95, interval = "prediction")
```


