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


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

# Question 1

#### (a)

```{r}
dat <- read.table(file="data/baseball.txt", header=TRUE)
plot(dat$EBP, dat$OBA)
poly.6 <- lm(OBA ~ poly(EBP, 6), data=dat)
summary(poly.6)
poly.2 <- lm(OBA ~ poly(EBP, 2), data=dat)
summary(poly.2)

RSS_H1 <- deviance(poly.6)
RSS_H0 <- deviance(poly.2)
p <- length(coef(poly.6))
q <- length(coef(poly.2))
n <- nrow(dat)
fobs <- (RSS_H0 - RSS_H1) / (p - q) / (RSS_H1 / (n - p))
1 - pf(fobs, p - q, n - p)

poly.2 <- lm(OBA ~ EBP + I(EBP^2), data=dat)
summary(poly.2)
```
The question says to look at polynomial model of order up to 6. So perhaps a good starting point is to look at a polynomial model of order 6 and look at the $p$-value from the summary output. The highest order which is significant is 2. A formal comparison of the polynomial model of order 2 vs order 6 using a F-test yields a $p$-value of 0.6265. Thus, the second-order polynomial model 
$$ Y_i = \beta_0 + \beta_1x_i + \beta_2 x_i^2 + \epsilon_i$$
where $x_i$ and $Y_i$ is the EBP and OBP for the $i$-th observation may be a good fit to the data. The $t$-test for $H_0: \beta_2$ yields a $p$-value of 0.00499 suggesting that we should not simplify the model further (using `base` plot in R). 

The following is just a sanity check.
```{r}
plot(dat$EBP, dat$OBA)
beta2 <- coef(poly.2)
poly.6 <- lm(OBA ~ poly(EBP, 6, raw=TRUE), data=dat)
beta6 <- coef(poly.6)
poly2 <- function(x) beta2[1] + beta2[2]*x + beta2[3]*x^2
poly6 <- function(x) beta6[1] + beta6[2]*x + beta6[3]*x^2 + beta6[4]*x^3 + beta6[5]*x^4 + beta6[6]*x^5 + beta6[7]*x^6
curve(poly2, from=min(dat$EBP), to=max(dat$EBP), add=T)
curve(poly6, from=min(dat$EBP), to=max(dat$EBP), add=T, col="red")
```
In the above plot, the red and black lines correspond to the fitted lines for `poly.6` and `poly.2`, respectively. It's not clear whether `poly.2` is the best model however `poly.6` is clearly an overfitted model. 

#### (b)
The fitted second order polynomial model is 
$$ \hat{Y}_i = 625.729571 - 2.679007 x_i + 0.008120x_i^2. $$

#### (c)

```{r}
predict(poly.2, newdata=data.frame(EBP=250, EBP2=62500), interval="prediction", level=0.9)
```

The predicted OBA at EBP of 250 is 463.4753. The associated 90% prediction interval is (417.9762, 508.9744). 

## Question 2

#### (a)

```{r}
dat <- read.csv("data/diabetes.csv")
skimr::skim(dat)
```
There data contains all numerica variables (or to be more specific 5 variables including the response are integers).

#### (b)

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

The explanatory variable that is most highly correlated with the response is the body mass index (correlation coefficient of 0.586).

#### (c)

```{r}
M0 <- lm(y ~ 1, data=dat)
Mf <- lm(y ~ ., data=dat)
```

#### (d)

The first model can be found by using:
```{r}
MforwardAIC <- step(M0, direction="forward", trace=0, k=2, scope= list(lower=M0, upper=Mf))
summary(MforwardAIC)
```

For the second model using F-test you must do it by "hand":
```{r}
drop1(Mf, test="F")
MstepFtest <- update(Mf, . ~ . - age)
drop1(MstepFtest, test="F")
MstepFtest <- update(MstepFtest, . ~ . - hdl)
add1(MstepFtest, test="F", scope=Mf)
drop1(MstepFtest, test="F")
MstepFtest <- update(MstepFtest, . ~ . - glu)
add1(MstepFtest, test="F", scope=Mf)
drop1(MstepFtest, test="F")
MstepFtest <- update(MstepFtest, . ~ . - tch)
add1(MstepFtest, test="F", scope=Mf)
drop1(MstepFtest, test="F")
summary(MstepFtest)
```

For the third model:

```{r}
MforwardAIC <- step(Mf, direction="backward", trace=0, k=2, scope= list(lower=M0, upper=Mf))
summary(MforwardAIC)
```

And for the third model:

```{r}
MstepBIC <- step(Mf, direction="both", trace=0, k=log(nrow(dat)), scope= list(lower=M0, upper=Mf))
summary(MstepBIC)
```

In this case all four models choose the same (remember this is not always the case). The chosen model includes `bmi`, `ltg`, `map`, `tc`, `sex` and `ldl`. 