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

We use the following packages:
```{r message=FALSE, warning=FALSE}
library(boot)
library(mice)
library(dplyr)
library(magrittr)
library(ggplot2)
```

I fix the random seed to `123`
```{r}
set.seed(123)
```

---

1. **Use the following function to calculate the bootstrapped estimate for the correlation between `age` and `weight` on the `boys` data set. Use `R=1000` bootstrap samples.**

```{r}
corfun <- function(data, i){
  data.sample <- data[i,]
  stat <- cor(data.sample$age, data.sample$wgt, use = "pairwise.complete.obs")
	return(stat)
}
```

To apply the `corfun()` function in the bootstrap function `boot()`, we run:
```{r cache = TRUE}
bootstr.cor <- 
  boys %>%
  boot(corfun, R = 1000)
```

Our object `bootstr.cor` contains all the bootstrap information about the correlation parameter

---

2. **Explore the contents of the `bootstr.cor` object. For example, use function `attributes()` to see the listed dimensions that are within the `bootstr.cor` object and the `class` of the object.**
```{r}
attributes(bootstr.cor)
```
we see that the `bootstr.cor` object has class `boot`. 
```{r}
head(bootstr.cor$data)
```

The most interesting listed dimensions are the following:

- The `$t0` contains the original estimate for the data 
```{r}
bootstr.cor$t0
```

- The `$t` contains the individual estimate for every bootstrapped sample
```{r}
head(bootstr.cor$t)
```

The average bootstrapped estimate would then be
```{r}
bootstr.cor$t %>% mean
```


- The `$data` has the original data set that has been served to function `boot()`:
```{r}
head(bootstr.cor$data)
```

---

3. **Plot the histogram of the individual bootstrapped estimates for every bootstrapped sample (i.e. `$t`).**
```{r}
plotdata <- data.frame(t = bootstr.cor$t)

plotdata %>%
  ggplot(aes(t)) + 
  geom_histogram(bins = 20) + 
  geom_density(col = "orange") + 
  geom_vline(xintercept = bootstr.cor$t0, col = "orange", linetype = "dashed", linewidth = 1.5)
```

---

4. **Add a column to boys that indicates whether boys are overweight by taking the boundary `bmi > 25`.** 

Remember that body mass index takes the ratio weight-over-height-squared into account and the validity of the interpretation does change for children vs. adults.

We can use function `mutate()` to assign values to a new column (variable) `ovw` for boys that have a `bmi <= 25` (smaller or equal) or `bmi > 25` (larger). The `if_else()` function can be directly applied to call `if_else(condition, true, false)`.
```{r}
boys2 <- boys %>%
    mutate(ovw = if_else(bmi <= 25, "Not Overweight", "Overweight")) 
```
There are now 20 `Overweight` boys:
```{r}
boys2 %$%
  table(ovw)
```

---

5. **Bootstrap a $X^2$-test that evaluates the distribution of overweight `ovw` boys over the regions `reg`. **
We can copy the previous function and replace the line for our test-statistic:
```{r}
X2fun <- function(data, i){
  data.sample <- data[i,]
  tab <- data.sample %$%
    table(ovw, reg) 
  X2 <- suppressWarnings(chisq.test(tab))
	return(X2$statistic)
}
```

I use `suppressWarnings()` to suppress the warnings the $X^2$ test is going to give us, because some of the expected cell-frequencies are below 5. I only return the chi-squared estimate by asking for `X2$statistic`. This is because function `boot` expects an estimate as return, not a whole list of estimates. 
 
```{r}

bootstr.X2 <- 
  boys2 %>%
  boot(X2fun, R = 1000)

bootstr.X2

```

---

6. **Plot the histogram of the individual bootstrapped estimates for every bootstrapped sample (i.e. `$t`).**
```{r}
plotdata <- data.frame(t = bootstr.X2$t)

plotdata %>%
  ggplot(aes(t)) + 
  geom_histogram(bins = 40) + 
  geom_vline(xintercept = bootstr.X2$t0, col = "orange", linetype = "dashed", linewidth = 1.5)
```

The sampling distribution of the estimate is quite skewed.

---
 
7. **Do a bootstrap of the regression estimates of the following model:** 

- `lm(wgt ~ age + hgt + I(hgt^2))`

The function that would give us the regression estimates for the parameters is e.g. :
```{r}
lmfun <- function(data, i){
  data.sample <- data[i,]
  coef <- data.sample %$%
    lm(wgt ~ age + hgt + I(hgt^2)) %>%
    coef()
	return(coef)
}
```

and if we plug that function into `boot()` for `R = 1000` bootstrap samples, we obtain
```{r}
bootstr.lm <- 
  boys %>%
  boot(lmfun, R = 1000)

bootstr.lm
```

Again, the observed data regression estimate is:
```{r}
bootstr.lm$t0
```

---

8. **Create a histogram of the individual bootstrapped estimates **
There are two ways to do this. 
1. Using the plot options of the class `boot`.
2. Construct the plots manually by extracting the regression coefficients from the boot class object.


Plot the distribution of the intercept by using the `index` argument (1 refers to the intercept and subsequent numbers to the remaining regression coefficients / slopes).
```{r}
plot(bootstr.lm, index = 1)
```
Plot the distribution of the slope of `wgt`:
```{r}
plot(bootstr.lm, index = 2)
```


Plot the distribution of the slope of `hgt`:
```{r}
plot(bootstr.lm, index = 3)
```
Plot the distribution of the slope of `hgt^2`:
```{r}
plot(bootstr.lm, index = 4)
```


**Option 2: manually**
Look at the code to see if you understand what it does. 
```{r}
est <- bootstr.lm$t
colnames(est) <- (c("Intercept", "wgt", "hgt", "hgt^2"))

reshape2::melt(est) %>%
  ggplot(aes(value)) + facet_wrap(~Var2, scales = 'free_x') +
  geom_histogram()
```


9. **Calculate the confidence intervals for the intercept and slopes**

We can use function `boot.ci()` to obtain the confidence intervals from our estimates. Remember that there we have 4 statistics: one intercept and three slopes. The bootstrap confidence interval of the intercept can be obtained with the argument `index = 1` and we choose the basic interval:
```{r}
boot.ci(bootstr.lm, type = "basic", index = 1)
```


```{r}
boot.ci(bootstr.lm, type = "basic", index = 2)
```
```{r}
boot.ci(bootstr.lm, type = "basic", index = 3)
```
```{r}
boot.ci(bootstr.lm, type = "basic", index = 4)
```

We can also use the `quantiles()` for the `$t` dimensions. For example, to obtain the cut-offs for `.025` and `.975` from the empirical distribution for the intercept, use:
```{r}
quantile(est[, "Intercept"], c(0.025,0.975))
```
or, for all estimates when using `apply()`:
```{r}
apply(est, 2, function(x) quantile(x, c(0.025,0.975)))
```

The `quantile()` function yields a slightly different result from the percentile bootstrap estimates, because the quantile function uses a quantile algorithm. Read all about these algorithms in the help. 

To calculate the percentile bootstrap CI by hand:
```{r}
sort(est[, "Intercept"])[c(25, 976)]
```

---

End of Practical

