---
title: "Lab Sheet 3 Solution"
author: "Jennifer Chan"
date: "06/03/2020"
output: html_document
---


```{r setup, warning=F, message=F}
library(tidyverse)
library(broom)
library(GGally)
knitr::opts_chunk$set(echo = TRUE, tidy = FALSE, fig.align = 'center', fig.pos = 'H')
```

## Question 1

### (a)
```{r}
x1 <- c(0.4, 0.4, 3.1, 0.6, 4.7, 1.7, 9.4, 10.1, 11.6, 12.6, 13.8, 10.9, 23.1, 29.9)
x2 <- c(53, 23, 19, 34, 24, 65, 44, 31, 29, 58, 55, 37, 46, 51)
x3 <- c(158, 163, 37, 157, 59, 123, 46, 117, 173, 112, 117, 111, 114, 124)
y <- c(64, 60, 71, 61, 54, 77, 81, 93, 93, 51, 60, 76, 96, 99)
dat <- data.frame(y, x1, x2, x3)
var(dat)
cor(dat)
ggpairs(dat)
```

The scatterplot does not show any outliers. Variables $x_1$ and $y$ as well as $x_1$ and $x_2$ show some positive correlation (linear dependence) -- this is also seen in the empirical correlation matrix. The scatterplot between $x_1$ and $y$ shows some fanning in and between $x_1$ and $x_2$ some fanning out.

### (b)

```{r}
fit1 <- lm(y ~ ., data=dat)
summary(fit1)
```
The fitted value and residual of the first observation are `r fit1$fitted[1]` and `r fit1$res[1]`, respectively.

### (c)

The 95\% confidence interval for $\beta_3$ is 
```{r}
confint(fit1, 4, level=0.95)
```

### (d)

The squared value of the multiple correlation coefficient is `r glance(fit1)[1,1]`.

### (e) 

```{r}
n <- nrow(dat)
p <- length(coef(fit1))
h <- lm.influence(fit1)$hat
which(rank(h)==n - 4 + 1)
any(h > 2 * p / n)
```
Thus, observation number `r which(rank(h)==n-4+1)` has fourth largest leverage value. There are no high leverage points in this data set.

### (f) 

$x_3$ can be dropped because it has largest p-value, which is larger than 5\%. We cannot exclude $x_2$ directly by looking at the summary output in (b) but can answer this through the following general F-test
```{r}
fit2 = lm(y ~ x1, data=dat)
anova(fit2, fit1)
``` 
Because the p-value of this general F-test, which tests $H_0: \beta_2=\beta_3=0$ against $H_1: \beta_2 \neq 0 \text{ and/or } \beta_3 \neq 0$ is greater than 0.05 we retain $H_0$, implying that both $x_2$ and $x_3$ can be simultaneously dropped.

### (g)

Model `fit2` in (h) is a reasonable model with $x_1$ alone. Its $R^2$ value is `r glance(fit2)[1,1]`, only a little smaller than the full model in (b). The estimated error standard deviation of the final model is `r sigma(fit2)`, which is smaller than the `r sigma(fit1)` in model `fit1`.

### (h)

```{r}
predict(fit2, data.frame(x1=40), interval="prediction", level=0.94)
```
Thus the 94% prediction interval for $Y|x_1=40$ calculates to $[69.75,148.55]$ (2dp).

### (i)

```{r}
max(cooks.distance(fit1))
max(cooks.distance(fit2))
```

The values differ because Cook's distance is a function of the response variable and all explanatory variables, i.e. the full model uses for the calculation of Cook's distance $y,x_1,x_2,x_3$, whereas Cook's distance for the final model is only a function of $y$ and $x_1$.