<style type="text/css">

body{ /* Normal  */
      font-size: 12px;
  }
td {  /* Table  */
  font-size: 12px;
}
h1.title {
  font-size: 20px;
  color: DarkBlue;
}
h1 { /* Header 1 */
  font-size: 18px;
}
h2 { /* Header 2 */
    font-size: 16px;
}
h3 { /* Header 3 */
  font-size: 14px;
}
code.r{ /* Code block */
    font-size: 12px;
}
pre { /* Code block - determines code spacing between lines */
    font-size: 12px;
}
</style>

---

# Aim
In this document we explore the confidence validity of bootstrap estimates. The goal of this investigation is to study the possibility of using the bootstrap to simulate a monte carlo experiment on a single (small) data set. 

---

# Set up
This section covers the global parameters and the necessary packages for executing the simulation study. 

---

## Fixing the random seed
Later on, we will be using random numbers. We fix the seed of the random number generator, so that we can always reproduce the results obtained in this document. 
```{r}
set.seed(123)
```
The random seed is a number used to initialize the pseudo-random number generator. The initial value (the seed) itself does not need to be random. The resulting process is random because the seed value is not used to generate randomness - it merely forms the starting point of the algorithm for which the results are random.

---

## Loading the required packages
```{r message=FALSE, warning=FALSE}
library(dplyr)    # Data manipulation
library(magrittr) # Pipes
library(purrr)    # Functional programming
library(mice)     # Data imputation
library(psych)    # Descriptives
library(knitr)    # 
library(kableExtra) #Cool tables
```

---

## Number of simulations
In this document we perform a Monte Carlo experiment. Here we define the number of replications that makes up this experiment. 
```{r}
nsim = 10000
```

---

# Establish the origin
We define a completed version of the `mice::boys` data set as the origin - conventionally we would consider this to be an infinitely large population - against which we compare our estimates. 
```{r}
origin <- mice(boys, m=1, print = FALSE) %>% # impute only once
  mice::complete() # extract completed data
```

---

## Model of interest
We define the following model to be evaluated in simulation
$$\text{wgt} = \beta_0 + \beta_1\text{hgt} + \beta2\text{age} + \epsilon$$
For the `origin` set, we extract the linear model as follows:
```{r}
truemodel <- origin %$%
  lm(wgt ~ hgt + age)
```

---

# Draw the bootstrap sets
We use the `replicate()` function to draw `nsim =`$10^4$ bootstrap sets from the 748 rows in the `origin`.
```{r cache=TRUE}
simdata <- replicate(nsim, 
                     origin[sample(1:748, 748, replace = TRUE), ], 
                     simplify = FALSE)
```

---

# Bias
First we evaluate the model of interest on each of the drawn bootstrap samples. We do this with the `purrr::map()` function, which *maps* the evaluation on each of the listed elements in the `simdata` object.
```{r}
model <- simdata %>% 
  map(~ lm(wgt ~ hgt + age, data = .x))
```
Then, we extract the estimates from the `model` evaluations.
```{r cache=TRUE}
estimates <- model %>% 
  map(coef) %>% 
  do.call(rbind, args = .) # bind rows into matrix
```

The obtained average estimates are 
```{r cache=TRUE}
estimates %>% 
  describe(quant = c(.025, .975)) %>%
  .[, c(2:4, 8:9, 11:12, 14:15)] %>%
  kable("html") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = F)
```

We find that the estimates display on average a low bias. 
```{r}
estimates %>% 
  colMeans
```


---

# Confidence Validity

---

## Approach 1: CI's with base `R`

We use the `confint()` function to extract the confidence intervals for each of the linear models in the `model` object. 
```{r cache=TRUE}
ci <- model %>% 
  map(confint)

```
We then calculate the proportion of confidence intervals that cover the `truemodel` parameters
```{r cache=TRUE}
cov <- ci %>% 
  map(function(x) x[, 1] <= coef(truemodel) & coef(truemodel) <= x[, 2]) %>%
  do.call(rbind, args = .) # bind rows into matrix
```

---

## Approach 2: manually calculated CI's
For unknown $\sigma^2_\beta$, the $1-\alpha$ confidence interval for $\beta$ is defined as  
\[\mu = \hat{\beta} \pm t_{n-1, a/2}\frac{S}{\sqrt{N}}.\]. The following code extracts the estimates and the standard error
```{r cache=TRUE}
manual <- model %>% 
  map(function(x) cbind(vcov(x) %>% diag %>% sqrt(.), coef(x))) %>% # est and resp. vars
  map(function(x) cbind(x[, 2] - qt(.975, 747) * x[, 1],  x[, 2] + qt(.975, 747) * x[, 1])) %>%
  map(function(x) x[, 1] <= coef(truemodel) & coef(truemodel) <= x[, 2]) %>%
  do.call(rbind, args = .) # bind rows into matrix
```

---

## Approach 3: bootstrap CI's cf. bootstrap SE
Instead of drawing the standard errors from every model, we can also make use of the bootstrap estimate of the standard error - i.e. the square root of the variance of the vector of estimates. This gives us the empirical standard deviation of the realized bootstrap sampling distribution of the estimates. 
```{r}
manual2 <- model %>% 
  map(function(x) cbind(sqrt(diag(var(estimates))), coef(x))) %>% # sd of 
  map(function(x) cbind(x[, 2] - qt(.975, 747) * x[, 1],  x[, 2] + qt(.975, 747) * x[, 1])) %>%
  map(function(x) x[, 1] <= coef(truemodel) & coef(truemodel) <= x[, 2]) %>%
  do.call(rbind, args = .)
```

---

## Approach 4: bootstrap CI's cf average CI
Another approach is to calculate the bootstrap CI based on the quantiles of the vectors of estimates. Simply dividing the confidence interval width would then yield an empirical equivalent to $\pm t_{n-1, a/2}\frac{S}{\sqrt{N}}$. 
```{r cache = TRUE}
# confidence intervals and widths
intercept <- data.frame(est = mean(estimates[, 1]), 
                        ciw = diff(quantile(estimates[, 1], probs = c(.025, .975)))) %>%
  mutate(low = est - ciw/2,
         up = est + ciw/2)
hgt <- data.frame(est = mean(estimates[, 2]), 
                  ciw = diff(quantile(estimates[, 2], probs = c(.025, .975)))) %>%
  mutate(low = est - ciw/2,
         up = est + ciw/2)
age <- data.frame(est = mean(estimates[, 3]), 
                  ciw = diff(quantile(estimates[, 3], probs = c(.025, .975)))) %>%
  mutate(low = est - ciw/2,
         up = est + ciw/2)
# coverages
covint <- data.frame(est = estimates[, 1], 
                     low = estimates[, 1] - intercept$ciw/2,
                     up = estimates[, 1] + intercept$ciw/2) %>%
  mutate(cov = low <= coef(truemodel)[1] & coef(truemodel)[1] <= up)
covhgt <- data.frame(est = estimates[, 2], 
                     low = estimates[, 2] - hgt$ciw/2,
                     up = estimates[, 2] + hgt$ciw/2) %>%
  mutate(cov = low <= coef(truemodel)[2] & coef(truemodel)[2] <= up)
covage <- data.frame(est = estimates[, 3], 
                     low = estimates[, 3] - age$ciw/2,
                     up = estimates[, 3] + age$ciw/2) %>%
  mutate(cov = low <= coef(truemodel)[3] & coef(truemodel)[3] <= up)
```

---

## Results
Below are the coverages for each of the investigated scenarios. Results are depicted for the coverage of the confidence intervals over the regression estimates for $10^4$ simulations.
```{r}
output <- rbind(colMeans(cov), 
                colMeans(manual), 
                colMeans(manual2), 
                colMeans(cbind(covint$cov, covhgt$cov, covage$cov)))
rownames(output) <- c("confint", 
                      "manual", 
                      "bootstrap SE", 
                      "bootstrap CI")
kable(output, "html") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = F)
```

---

# Replication
The replication of these results [can be found here](BootstrapCI_Validity_replication.html)

---

END OF DOCUMENT

---

```{r}
sessionInfo()
```


