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


```{r setup, warning=F, message=F}
library(tidyverse)
dat <- read.table("data/ratliver.txt", header=T)
knitr::opts_chunk$set(echo = TRUE, tidy = FALSE, fig.align = 'center', fig.pos = 'H')
```

## Question 1

### (a)

```{r}
fit <- lm(y ~ x1 + x2 + x3, data=dat)
coef(fit)
deviance(fit)
sigma(fit)^2
```

The fitted model is 
$$Y = 0.2659 -0.0212 x_1 + 0.0143 x_2 +4.178 x_3.$$

$\text{RSS} = 0.0896$ and $\hat{\sigma}^2 = 0.07729^2 = 0.0060.$

### (b)

```{r}
n <- nrow(dat)
p <- 4
a <- qchisq(0.95, n - p)
b <- qchisq(0.05, n - p)
RSS <- deviance(fit)
c(sqrt(RSS / a), sqrt(RSS / b))
```

A 90% CI for $\sigma$ is $(0.060, 0.111)$.

### (c)

```{r}
predict(fit, data.frame(x1=165, x2=8, x3=0.85), level=0.95, interval="confidence")
```

A 95% CI for the expected value of $Y$ (i.e.\ for $E(X)$) given the stated explanatory
variables is given through $$(0.343, 0.509).$$

### (d)

```{r}
summary(fit)
```

From the output $\beta_2$ has a large associated $p$-value and so $x_2$ can be dropped from the model.

### (e)

```{r}
which(lm.influence(fit)$hat > 2 * p /n)
lm.influence(fit)$hat[3]
```

The third point [3] is a high leverage point as $h_{33} > 2\times 4/19 = 0.421.$ Note that the $x_3$ value is high but this is not the only point with $x_3 =1.00$. It is
the combination of coordinates that make this point a high leverage point.

### (f) 

```{r}
which(cooks.distance(fit) > qf(0.5, p, n - p))
cooks.distance(fit)[3]
c(p, n - p)
qf(0.5, p, n - p)
```

The $F_{4,15}$ median is 0.878 so [3] is an outlier because $D_3 = 0.9296>F^{-1}_{4,15}(0.5)$.

### (g)

```{r}
fit1 <- lm(y ~ x1 + x2 + x3, data=dat[-3,])
summary(fit1)
```

Dropping [3] the model becomes
$$Y=0.3114 - 0.0078x_1 + 0.0090 x_2 + 1.4849 x_3$$
The variables only explain 2.1% of the $Y$ variability in this case.

### (h)

```{r}
confint(fit1, level=0.95)
```

All intervals contain 0 so any of the parameters individually can be set to 0.

### (i) and (j)

```{r}
sd(dat$y[-3])
```

The $y$-sample standard deviation ignoring [3] is 0.071768 whereas the
estimate from the multiple regression model is 0.07825. Thus the simple model
$Y_i = \beta_0 + \epsilon_i$ has smaller s.d than the multiple regression
model. Ignoring the outlier, it looks like the three variables are not useful 
in explaining $Y$.

### (k)

```{r}
fit0 <- lm(y ~ 1, data=dat[-3,])
(fobs <- (deviance(fit0) - deviance(fit1)) / 3 / (deviance(fit1) / (18 - 4)))
1 - pf(fobs, 3, 18 - 4)
```

The $p$-value from the F-test is 0.958 which suggests that all three covariates are not important to explaining the response.

### Aside

If you had a look a the pairwise scatter plot, you would have noticed that `x1` and `x3` are highly correlated.

```{r}
GGally::ggpairs(dat, progress=FALSE)
```