---
subtitle: "TMA4268 Statistical Learning V2018"
title: "Module 5: RESAMPLING"
author: "Mette Langaas, Department of Mathematical Sciences, NTNU"
date: "week 7 2018 (Version 22.05.2018)"
output: 
  # html_document:
  #   toc: true
  #   toc_float: true
   # beamer_presentation:
   #   #theme: "AnnArbor"
   #   fig_width: 3
   #   fig_height: 3
   #   keep_tex: true
###    incremental: true # or >* for one at a time
  pdf_document:
#    latex_engine: xelatex
    # toc: true
   # toc_depth: 2
#editor_options: 
#  chunk_output_type: console
---

Latest changes: (22.05: added toc for html version, changed from k to j for better notation for CV-formulas due to double booking of k.)

```{r setup, include=FALSE}
showsol<-TRUE
library(knitr)
opts_chunk$set(results="hold",tidy=TRUE,warning=FALSE,error=FALSE,message=FALSE,echo=TRUE)
knitr::opts_chunk$set(fig.width=6, fig.height=6) 
```

# <a id="intro"> Introduction</a>

## Learning material for this module

* James et al (2013): An Introduction to Statistical Learning. Chapter 5.  
* [Classnotes 12.02.2018](https://www.math.ntnu.no/emner/TMA4268/2018v/notes/TMA4268ClassnotesM5.pdf) 

## Move to

* [Introduction](#intro)
* [Cross-validation](#cv) and [Recommended exercises on cross-validation](#recexcv)
* [Bootstrapping](#boot) and [Recommended exercises on bootstrapping](#recexboot)
* [Summing up](#summing)
* [Further reading](#further)
* [Packages to install before knitting this R Markdown file](#Rpackages)

---

## What will you learn?

* What is model assessment and model selection?
* Ideal solution in a data rich situation.
* Cross-validation: validation set - LOOCV and $k$-fold CV - what is the best?
* Bootstrapping - how and why.
* Summing up
* The plan for the interactive lesson on Wednesday/Friday.

---

# Generalization performance of learning method

* prediction capacity on independent test data
* inference and understanding 

This is important both for 

### Model selection

estimate the _performance_ of different models (often different order of complexity within one model class) to _choose the best model_.

### Model assessment 

having chosen a final model, estimating its performance (prediction error) on new data.

---

### Example 
from ComplsoryEx 1- Problem 1

We aim to do _model selection_ in KNN-regression, where true curve is
$f(x)=-x+x^2+x^3$ with $x \in [-3,3]$. $n=61$ for the training data.

![Regression problem: training data and true regression curve](Prob1f1.png)

---

### KNN-regression

$n=61$ both for training and for test data (using same x-grid).

![Regression problem: MSE for training and test set, M=1000 versions](Prob1f3.png)


$K$ small: high complexity (left) and $K$ large: low complexity (right). 
 
---

# The bias-variance trade-off

<!-- <figure><img src="Probl1fig4negK.png" height="250px" width="250px" align="right" /><figcaption>Regression problem: bias-variance traceoff</figure> -->
![Regression problem: bias-variance traceoff](Prob1f4.png)

Left (high complexity): low squared-bias (red) and high variance (green).
Right (low complexity): high squared-bias (red) and low variance (green). 

---

## Loss functions 

reminder - we will use

* Mean squared error (quadratic loss) for regression problems: 
$$Y_i=f({\bf x}_i)+\varepsilon_i \text{ for }i=1,\ldots,n \text{ and } \text{MSE}=\frac{1}{n}\sum_{i=1}^n (y_i-\hat{f}({\bf x}_i))^2$$

* 0/1 loss for classification problems:
$$P(Y=j\mid {\bf x}_0) \text{ for }j=1,\ldots,K $$
and classify to the class with the highest probability $\hat{y}_i$. Then $$\frac{1}{n}\sum_{i=1}^n \text{I}(y_i \neq \hat{y}_i)$$

## The challenge

our example was based on simulated data, so I had unlimited access to data. Now, let us move to real data!

---

# Data rich situation

If we had a large amount of data we could divide our data into three parts:

* training set: to fit the model
* validation set: to select the best model (aka model selection)
* test set: to assess how well the model fits on new independent data (aka model assessment)

**Q**: Why not enough with training and validation (why also test)?

**A**: We will be a bit optimistic if we report the error that we know we have chosen to be the smallest.


---

If this is the case - great - then you do not need Module 5. But, this is very seldom the case - so we will study other solutions based on efficient sample reuse with _resampling_ data.

An alternative strategy for model selection (using methods penalizing model complexity) is covered in Module 6.

First we look at _crossvalidation_, then at _bootstrapping_.

---

# Cross-validation (CV)

* the validation set approach
* leave one out cross validation (LOOCV)
* 5 and 10 fold crossvalidation (CV)
* selection bias - all elements of a model selection strategy need to be within the CV-loop
* recommended exercises

---

## The validation set approach

Consider the case when you have a data set consisting of $n$ observations. 

To fit a model and to test its predictive performance you randomly divide the data set into two parts ($n/2$ sample size each): 

* a _training set_ (to fit the model) and 
* a _validation set_ (to make predictions of the response variable for the observations in the validation set)

First: focus on model selection 

---

### Regression model selection example: validation set error
Auto data set (library `ISLR`): predict `mpg` (miles pr gallon) using polynomial function of `horsepower` (of engine), $n=392$.
<!-- ![One validation set](validset1.png) -->

```{r,echo=FALSE,eval=TRUE}
library(ISLR)
library(ggplot2)
library(ggpubr)
library(reshape2)
set.seed(123)
n=dim(Auto)[1]

testMSEmat=matrix(ncol=11,nrow=10)
for (newsample in 1:11)
{
  trainid=sample(n,n/2)
  for (polydeg in 1:10)
  {
    lm.fit=lm(mpg~poly(horsepower,polydeg),data=Auto,subset=trainid)
    lm.pred=predict(lm.fit,Auto)[-trainid]
    testMSEmat[polydeg,newsample]=mean((Auto$mpg[-trainid]-lm.pred)^2)
  }
}
yrange=c(15,28)

plotdf=data.frame("testMSE"=testMSEmat[,1],"degree"=1:10)
g0=ggplot(plotdf,aes(x=degree))+geom_line(y=testMSEmat[,1])+scale_y_continuous(limits = yrange)+scale_x_continuous(breaks=1:10)+labs(y="Validation set")
g0
```

---

### Regression example: validation set error for many random divisions

<!-- ![Original + 10 new validation sets](valid10.png) -->
```{r,eval=TRUE,echo=FALSE}
cols=rainbow(10)
g1=g0+geom_line(aes(x=1:10,y=testMSEmat[,2]),colour=cols[1])
g1=g1+geom_line(aes(x=1:10,y=testMSEmat[,3]),colour=cols[2])
g1=g1+geom_line(aes(x=1:10,y=testMSEmat[,4]),colour=cols[3])
g1=g1+geom_line(aes(x=1:10,y=testMSEmat[,5]),colour=cols[4])
g1=g1+geom_line(aes(x=1:10,y=testMSEmat[,6]),colour=cols[5])
g1=g1+geom_line(aes(x=1:10,y=testMSEmat[,7]),colour=cols[6])
g1=g1+geom_line(aes(x=1:10,y=testMSEmat[,8]),colour=cols[7])
g1=g1+geom_line(aes(x=1:10,y=testMSEmat[,9]),colour=cols[8])
g1=g1+geom_line(aes(x=1:10,y=testMSEmat[,10]),colour=cols[9])
g1=g1+geom_line(aes(x=1:10,y=testMSEmat[,11]),colour=cols[10])
g1
g1valid=g1
#ggarrange(g0,g1)
```

---

### Drawbacks with the validation set approach

* high variability of validation set error (which we think of as estimate for test set error) - since this is dependent on which observation are included in the training and validation set
* smaller sample size for model fit - since not all observations can be in the training set
* the validation set error may tend to overestimate the test set error for a model that is fit on the full data set (because - the more data the lower error, and here our training set is half of our data set).

---

## Leave-one-out cross-validation (LOOCV)

* If the data is very limited and the division of the data into two parts is unreasonable, leave-one-out cross-validation (LOOCV) can be used. 
* In LOOCV one observation at a time is left out and makes up the test set. 
* The remaining $n-1$ observations make up the training set. 
* The procedure of model fitting is repeated $n$ times, such that each of the $n$ observations is left out once. 
* The total prediction error is the mean across these $n$ models.

$$\text{MSE}_i=(y_i-\hat{y}_i)^2$$

$$CV_{n}=\frac{1}{n}\sum_{i=1}^n \text{MSE}_i$$

---

### Regression example: LOOCV

<!-- ![LOOCV](loocv.png) -->

```{r,echo=FALSE,eval=TRUE}
library(ISLR)
library(boot)
library(ggplot2)
set.seed(123)
n=dim(Auto)[1]
testMSEvec=NULL
start=Sys.time()
for (polydeg in 1:10)
  {
    glm.fit=glm(mpg~poly(horsepower,polydeg),data=Auto)
    glm.cv1=cv.glm(Auto, glm.fit,K=n)
    testMSEvec=c(testMSEvec,glm.cv1$delta[1])
}
stopp=Sys.time()
yrange=c(15,28)
plotdf=data.frame("testMSE"=testMSEvec,"degree"=1:10)
g0=ggplot(plotdf,aes(x=degree,y=testMSE))+geom_line()+geom_point()+scale_y_continuous(limits = yrange)+scale_x_continuous(breaks=1:10)+labs(y="LOOCV")
g0
```

---

\tiny
```{r,echo=TRUE,eval=FALSE}
library(ISLR)
library(boot)
library(ggplot2)
set.seed(123)
n=dim(Auto)[1]
testMSEvec=NULL
start=Sys.time()
for (polydeg in 1:10)
  {
    glm.fit=glm(mpg~poly(horsepower,polydeg),data=Auto)
    glm.cv1=cv.glm(Auto, glm.fit,K=n)
    testMSEvec=c(testMSEvec,glm.cv1$delta[1])
}
stopp=Sys.time()
yrange=c(15,28)
plotdf=data.frame("testMSE"=testMSEvec,"degree"=1:10)
g0=ggplot(plotdf,aes(x=degree,y=testMSE))+geom_line()+geom_point()+scale_y_continuous(limits = yrange)+scale_x_continuous(breaks=1:10)+labs(y="LOOCV")
g0
```
\normalsize

---

### Issues with leave-one-out cross-validation

* Good: 
    + no randomness in training/validation splits!
    + little bias since nearly the whole data set used for training (compared to half for validation set approach)
    
* Bad: 
    + expensive to implement - need to fit $n$ different models - however nice formula for linear model LOOCV - but not genereally so
    + high variance since: two training sets only differ by one observation - which makes estimates from each fold are highly correlated and this can lead to that their average can have high variance.

$$\text{Var}(\sum_{i=1}^na_iX_i+b)=\sum_{i=1}^n\sum_{j=1}^n a_ia_j\text{Cov}(X_i,X_j)$$
$$=\sum_{i=1}^na_i^2\text{Var}(X_i)+2\sum_{i=2}^n \sum_{j=1}^{i-1}
a_ia_j\text{Cov}(X_i,X_j).$$

---

### LOOCV for multiple linear regression

$$ CV_{n}=\frac{1}{n}\sum_{i=1}^n \left( \frac{y_i-\hat{y}_i}{1-h_{ii}} \right) ^2$$

where $h_i$ is the $i$th diagonal element (leverage) of the hat matrix ${\bf H}={\bf X}({\bf X}^T{\bf X})^{-1}{\bf X}^T$.

---

## $k$-fold cross-validation

See class notes for drawings!

---

### Formally

* Indices of observations - divided into $k$ folds: $C_1, C_2, \ldots, C_k$.
* $n_k$ elements in each fold, if $n$ is a multiple of $k$ then $n_k=n/k$.

$$\text{MSE}_k=\frac{1}{n_k}\sum_{i\in C_k}(y_i-\hat{y}_i)^2$$
where $\hat{y}_i$ is the fit for observation $i$ obtained from the data with part $k$ removed. 

$$CV_{k}=\frac{1}{n} \sum_{j=1}^k n_j \text{MSE}_j$$

Observe: setting $k=n$ gives LOOCV.

---

### Regression example: $5$ and $10$-fold cross-validation

```{r,echo=FALSE,eval=TRUE}
library(ISLR)
library(boot)
library(ggplot2)
set.seed(123)
n=dim(Auto)[1]
testMSEvec5=NULL
testMSEvec10=NULL
start=Sys.time()
for (polydeg in 1:10)
  {
    glm.fit=glm(mpg~poly(horsepower,polydeg),data=Auto)
    glm.cv5=cv.glm(Auto, glm.fit,K=5)
    glm.cv10=cv.glm(Auto, glm.fit,K=10)
    testMSEvec5=c(testMSEvec5,glm.cv5$delta[1])
    testMSEvec10=c(testMSEvec10,glm.cv10$delta[1])
}
stopp=Sys.time()
yrange=c(15,28)
plotdf=data.frame("testMSE5"=testMSEvec5,"degree"=1:10)
g0=ggplot(plotdf,aes(x=degree,y=testMSE5))+geom_line()+geom_point()+scale_y_continuous(limits = yrange)+scale_x_continuous(breaks=1:10)+labs(y="CV")+ggtitle("5 and 10 fold CV")
g0+geom_line(aes(y=testMSEvec10),colour="red")+geom_point(aes(y=testMSEvec10),colour="red")+ggtitle("5 fold (black), 10 fold (red)")
```

---

\tiny

```{r,echo=TRUE,eval=FALSE}
library(ISLR)
library(boot)
library(ggplot2)
set.seed(123)
n=dim(Auto)[1]
testMSEvec5=NULL
testMSEvec10=NULL
start=Sys.time()
for (polydeg in 1:10)
  {
    glm.fit=glm(mpg~poly(horsepower,polydeg),data=Auto)
    glm.cv5=cv.glm(Auto, glm.fit,K=5)
    glm.cv10=cv.glm(Auto, glm.fit,K=10)
    testMSEvec5=c(testMSEvec5,glm.cv5$delta[1])
    testMSEvec10=c(testMSEvec10,glm.cv10$delta[1])
}
stopp=Sys.time()
yrange=c(15,28)
plotdf=data.frame("testMSE5"=testMSEvec5,"degree"=1:10)
g0=ggplot(plotdf,aes(x=degree,y=testMSE5))+geom_line()+geom_point()+scale_y_continuous(limits = yrange)+scale_x_continuous(breaks=1:10)+labs(y="CV")+ggtitle("5 and 10 fold CV")
g0+geom_line(aes(y=testMSEvec10),colour="red")+geom_point(aes(y=testMSEvec10),colour="red")+ggtitle("5 fold (black), 10 fold (red)")
```
\normalsize

---

```{r,echo=FALSE,eval=TRUE}
library(ISLR)
library(boot)
library(ggplot2)
set.seed(123)
n=dim(Auto)[1]
start=Sys.time()
testMSEmat=matrix(ncol=11,nrow=10)

for (polydeg in 1:10)
  {
    glm.fit=glm(mpg~poly(horsepower,polydeg),data=Auto)
    for (newsample in 1:11)
    {
        glm.cv10=cv.glm(Auto, glm.fit,K=10)
      testMSEmat[polydeg,newsample]=glm.cv10$delta[1]
    }
}

stopp=Sys.time()
yrange=c(15,28)

plotdf=data.frame("testMSE"=testMSEmat[,1],"degree"=1:10)
g0=ggplot(plotdf,aes(x=degree))+geom_line(y=testMSEmat[,1])+scale_y_continuous(limits = yrange)+scale_x_continuous(breaks=1:10)+labs(y="CV10")
#g0

cols=rainbow(10)
g1=g0+geom_line(aes(x=1:10,y=testMSEmat[,2]),colour=cols[1])
g1=g1+geom_line(aes(x=1:10,y=testMSEmat[,3]),colour=cols[2])
g1=g1+geom_line(aes(x=1:10,y=testMSEmat[,4]),colour=cols[3])
g1=g1+geom_line(aes(x=1:10,y=testMSEmat[,5]),colour=cols[4])
g1=g1+geom_line(aes(x=1:10,y=testMSEmat[,6]),colour=cols[5])
g1=g1+geom_line(aes(x=1:10,y=testMSEmat[,7]),colour=cols[6])
g1=g1+geom_line(aes(x=1:10,y=testMSEmat[,8]),colour=cols[7])
g1=g1+geom_line(aes(x=1:10,y=testMSEmat[,9]),colour=cols[8])
g1=g1+geom_line(aes(x=1:10,y=testMSEmat[,10]),colour=cols[9])
g1=g1+geom_line(aes(x=1:10,y=testMSEmat[,11]),colour=cols[10])
g1
#ggarrange(g0,g1)
```

---

### Issues with $k$-fold cross-validation

1. As for the validation set, the result may vary according to how the folds are made, but the variation is in general lower than for the validation set approach.
2. Computational issues: less work with $k=5$ or 10 than LOOCV.
3. The training set is (k-1)/k of the original data set - this will estimated of prediction error that are biased upwards.
4. This bias is the smalles when $k=n$ (LOOCV), but we know that LOOCV has high variance.
5. This is way often $k=5$ or $k=10$ is used as a compromise between 3 and 4.

---

### $k$-fold cross-validation in classification
- what do we need to change from our regression set-up?

Other loss function - 0/1 loss, or maybe 1-AUC?

---

## The right and the wrong way to do cross-validation

[ISL book slides, page 17](https://lagunita.stanford.edu/c4x/HumanitiesScience/StatLearning/asset/cv_boot.pdf): model assessment.

* We have a two-class problem and would like to use a simple classification method, however, 
* we have many possible predictors $p=5000$ and not so big sample size $n=50$.

We use this strategy to produce a classifier:

1. We calculate the correlation between the class label and each of the $p$ predictors, and choose the $d=25$ predictors that have the highest (absolute value) correlation with the class label. (We need to have $d<n$ to fit the logistic regression uniquely.)
2. Then we fit our classifier (here: logistic regression) using only the $d=25$ predictors.

How can we use cross-validation to produce an estimate of the performance of this classifier? 

**Q**: Can we apply cross-validation only to step 2? Why or why not?

---

Can we apply cross-validation only to step 2?

**A:** No, step 1 is part of the training procedure (the class labels are used) and must be part of the CV to give an honest estimate of the performance of the classifier. 

* Wrong: Apply cross-validation in step 2.
* Right: Apply cross-validation to steps 1 and 2.

We will see in the Recommended Exercises that doing the wrong thing can give a misclassification error approximately 0 - even if the "true" rate is 50%.

---

### Selection bias in gene extraction on the basis of microarray gene-expression data

Article by [Christophe Ambroise and Geoffrey J. McLachlan, PNAS 2002](http://www.pnas.org/content/99/10/6562): Direct quotation from the abstract of the article follows.

* In the context of cancer diagnosis and treatment, we consider the
problem of constructing an accurate prediction rule on the basis of
a relatively small number of tumor tissue samples of known type
containing the expression data on very many (possibly thousands)
genes. 
* Recently, results have been presented in the literature
suggesting that it is possible to construct a prediction rule from
only a few genes such that it has a negligible prediction error rate.
* However, in these results the test error or the leave-one-out
cross-validated error is calculated without allowance for the selection
bias. 

---

* There is no allowance because the rule is either tested
on tissue samples that were used in the first instance to select the
genes being used in the rule or because the cross-validation of the
rule is not external to the selection process; that is, gene selection
is not performed in training the rule at each stage of the crossvalidation
process.
* We describe how in practice the selection bias
can be assessed and corrected for by either performing a crossvalidation
or applying the bootstrap external to the selection
process. 
* We recommend using 10-fold rather than leave-one-out
cross-validation, and concerning the bootstrap, we suggest using
the so-called .632 bootstrap error estimate designed to handle
overfitted prediction rules. 
* Using two published data sets, we
demonstrate that when correction is made for the selection bias,
the cross-validated error is no longer zero for a subset of only a few
genes.

---

# <a id="recexcv"> Recommended exercises on cross-validation</a>

**Problem 1**: Explain how $k$-fold cross-validation is implemented 

* draw a figure, 
* specify algorithmically what is done,
* and in particular how the "results" from each fold are aggregated,
* relate to one example from regression (maybe is complexity wrt polynomials of increasing degree in multiple linear regression or $K$ in KNN-regression?)
* relate to one example from classification (maybe is complexity wrt polynomials of increasing degree in logistic regression or $K$ in KNN-classification?)

Hint: the words "loss function", "fold", "training", "validation" are central.

---

**Problem 2**: What are the advantages and disadvantages of $k$-fold cross-validation relative to

* the validation set approach
* leave one out cross-validation (LOOCV)
* what are recommended values for $k$, and why?
    
Hint: the words "bias", "variance" and "computational complexity" should be included.

---

**Problem 3:**. Selection bias and the "wrong way to do CV".

The task here is to devise an algorithm to "prove" that the wrong way is wrong and that the right way is right. 

a. What are the steps of such an algorithm? Write down a suggestion. Hint: how do you generate data for predictors and class labels, how do you do the classification task, where is the CV in the correct way and wrong way inserted into your algorithm? Can you make a schematic drawing of the right and the wrong way? Hint: [ISL book slides, page 20+21](https://lagunita.stanford.edu/c4x/HumanitiesScience/StatLearning/asset/cv_boot.pdf) - but you can do better?

---

b. One possible version of this is presented in the R-code below. Go through the code and explain what is done in each step, then run the code and observe if the results are in agreement with what you expected. Make changes to the R-code if you want to test out different strategies.

```{r,eval=FALSE}
library(boot)
# GENERATE DATA
# reproducible
set.seed(4268)
n=50 #number of observations
p=5000 #number of predictors
d=25 #top correlated predictors chosen
kfold=10
#generating predictor data
xs=matrix(rnorm(n*p,0,4),ncol=p,nrow=n) #simple way to to uncorrelated predictors
dim(xs) # n times p
# generate class labels independent of predictors - so if all classifies as class 1 we expect 50% errors in general
ys=c(rep(0,n/2),rep(1,n/2)) #now really 50% of each
table(ys)

# WRONG CV - using cv.glm
# here select the most correlated predictors outside the CV
corrs=apply(xs,2,cor,y=ys)
hist(corrs)
selected=order(corrs^2,decreasing = TRUE)[1:d] #top d correlated selected
data=data.frame(ys,xs[,selected])
#apply(xs[,selected],2,cor,y=ys) yes, ave the most correlated
# then cv around the fitting of the classifier - use logistic regression and built in cv.glm function
logfit=glm(ys~.,family="binomial",data=data)
cost <- function(r, pi = 0) mean(abs(r-pi) > 0.5)
cvres=cv.glm(data=data,cost=cost,glmfit=logfit,K=kfold)
cvres$delta
# observe - near 0 misclassification rate

# WRONG without using cv.glm - should be similar (just added to see the similarity to the RIGHT version)
reorder=sample(1:n,replace=FALSE)
validclass=NULL
for (i in 1:kfold)
{
  neach=n/kfold
  trainids=setdiff(1:n,(((i-1)*neach+1):(i*neach)))
  traindata=data.frame(xs[reorder[trainids],],ys[reorder[trainids]])
  validdata=data.frame(xs[reorder[-trainids],],ys[reorder[-trainids]])
  colnames(traindata)=colnames(validdata)=c(paste("X",1:p),"y")
  data=traindata[,c(selected,p+1)]
  trainlogfit=glm(y~.,family="binomial",data=data)
  pred=plogis(predict.glm(trainlogfit,newdata=validdata[,selected]))
  print(pred)
  validclass=c(validclass,ifelse(pred > 0.5, 1, 0))
}
table(ys[reorder],validclass)
1-sum(diag(table(ys[reorder],validclass)))/n

# CORRECT CV
reorder=sample(1:n,replace=FALSE)
validclass=NULL
for (i in 1:kfold)
{
  neach=n/kfold
  trainids=setdiff(1:n,(((i-1)*neach+1):(i*neach)))
  traindata=data.frame(xs[reorder[trainids],],ys[reorder[trainids]])
  validdata=data.frame(xs[reorder[-trainids],],ys[reorder[-trainids]])
  colnames(traindata)=colnames(validdata)=c(paste("X",1:p),"y")
  foldcorrs= apply(traindata[,1:p],2,cor,y=traindata[,p+1]) 
  selected=order(foldcorrs^2,decreasing = TRUE)[1:d] #top d correlated selected
  data=traindata[,c(selected,p+1)]
  trainlogfit=glm(y~.,family="binomial",data=data)
  pred=plogis(predict.glm(trainlogfit,newdata=validdata[,selected]))
  validclass=c(validclass,ifelse(pred > 0.5, 1, 0))
}
table(ys[reorder],validclass)
1-sum(diag(table(ys[reorder],validclass)))/n
```


<!-- Here is one possibility to start from: -->

<!-- Generate data -->
<!-- a. Simulate high dimensional data (e.g. $p=500$? from independent or correlated normal variables) but with few samples (e.g. $n=50$).  -->
<!-- b. Randomly assign a class labels (how many classes? maybe only 2 to start with) - so the truth is that the misclassification rate can not get very small. What is the theoretical misclassification rate (for this random set)? -->

<!-- Classification task: -->
<!-- c. Decide on a way to choose a few of the predictors and  -->
<!-- d. perform a classification rule (logistic, LDA, naive Bayes?) on these predictors.  -->

<!-- Then, do CV ($k=5$?) on only d =wrong way, or c+d=right way. Report misclassification errors for both situations. -->

---

**Problem 4**: Trying out different versions of cross-validation with R. 

* Maybe you can devise a simulated data set, or use a data set from the `ILSR` library. 
* To see differences look at a large data set and a small data set (in the number of observations), and 
* focus on regression and on classification. 
* What if the optimal (lowest CV error) method has complexity $7$, but the CV error of complexity $2$ is not very different. What should we do then? Hint: rule of parsimony, Occams razor. Later (ch 6) we will look at the "1sd rule".

Use Section 5.3 in the ISL book as a starting point.

<!-- a. Compare candidate models -->

<!-- LOOCV with -->
<!-- $k$-fold CV with built in function: cv.glm -->
<!-- for loop -->

---

# The Bootstrap

* flexible and powerful statistical tool that
can be used to quantify _uncertainty_ associated with an estimator or statistical learning method
* we will look at getting an estimate for the standard error of a sample median and of a regression coefficient
* in Module 8 - bootstrapping is the core of the _ensemble method_ referred to at _bagging_=bootstrap aggregation,
* in TMA4300 Computation statistics - more on the bootstrap.

---

The inventor: Bradley Efron in 1979 - [see interview](https://www.youtube.com/watch?v=6l9V1sINzhE).

The name? _To pull oneself up by one's bootstraps_ from "The Surprising
Adventures of Baron Munchausen" by Rudolph Erich
Raspe:

_The Baron had fallen to the bottom of a deep lake. Just
when it looked like all was lost, he thought to pick himself
up by his own bootstraps._

**Idea: Use the data itself to get more information about a statistic (an estimator).**


<!-- ## Example: The mean $\bar{X}$ -->

<!-- Assume we have univariate random variables  $X_1, X_2, \ldots, X_n$, with some common mean $\mu$ and variance $\sigma^2$. -->

<!-- We can then form the average: $\bar{X}=\frac{1}{n} \sum_{i=1}^n X_i$, which we often use as an estimator for $\mu$. -->

<!-- If we want to estimate $\sigma^2$, a popular estimator is $S^2=\frac{1}{n-1}\sum_{i=1}^n (X_i-\bar{X})^2$, -->

<!-- If the observations are _independent_ (and identically distributed) we know that the variance of $\bar{X}$ is given by $\text{Var}(\bar{X})=\frac{\sigma^2}{n}$,  -->
<!-- which we may estimate as $\widehat{\text{Var}}(\bar{X})=\frac{S^2}{n}$.  -->

<!-- To form a confidence interval for $\mu$ we need to assume a distribution for the $X_i$s. If we assume that the $X_i$s are independent and normal $N(\mu,\sigma^2)$  and can get a 95% confidence interval from -->

<!-- $$ \bar{X}\pm t_{0.025,n-1}\frac{S}{\sqrt{n}}$$ -->

---

## Example: the standard deviation of the sample median?

Assume that we observe a random sample $X_1, X_2, \ldots, X_n$ from an unknown probability distribution $f$. We are interesting in saying something about the population median, and to do that we calculate the sample median $\tilde{X}$. But, how accurate is $\tilde{X}$ as an estimator?

The bootstrap was introduced as a computer-based method to estimate the standard deviation of an estimator, for example our estimator $\tilde{X}$.

But, before we look at the bootstrap method, first we assume that we know $F$ and can sample from $F$, and use simulations to answer our question.

---

```{r,eval=TRUE, echo=TRUE}
set.seed(123)
n=101
B=1000
estimator=rep(NA,B)
for (b in 1:B)
{
  xs=rnorm(n)
  estimator[b]=median(xs)
}
sd(estimator)
# approximation for large samples 
# (sd of median of standard normal)
1.253*1/sqrt(n)
```

---

```{r}
ggplot(data=data.frame(x=estimator),aes(x=x))+
  geom_density()
```

---

## Moving from simulation to bootstrapping

The bootstrap method is using the observed data to estimate the _empirical distribution_ $\hat{f}$, that is each observed value of $x$ is given probability $1/n$.

A _bootstrap sample_ $X^*_1,X^*_2,\ldots, X^*_n$ is a random sample drawn from $\hat{f}$.

A simple way to obtain the bootstrap sample is to _draw with replacement from $X_1, X_2, \ldots, X_n$_.

This means that our bootstrap sample consists of $n$ members of $X_1, X_2, \ldots, X_n$ - some appearing 0 times, some 1, some 2, etc.


---

```{r, eval=TRUE, echo=TRUE}
set.seed(123)
n=101
original=rnorm(n)
median(original)
boot1=sample(x=original,size=n,replace=TRUE)
table(table(boot1))
n-sum(table(table(boot1)))
median(boot1)
```

---

## The bootstrap algorithm for estimating standard errors

1. $B$ bootstrap samples: drawn with replacement from the original data 

2. Evaluate statistic: on each of the $B$ bootstrap samples to get $\tilde{X}^*_b$ for the $b$th bootstrap sample.

3. Estimate squared standard error by:
$$\frac{1}{B-1}\sum_{b=1}^B (\tilde{X}^*_b-\frac{1}{B}\sum_{b=1}^B \tilde{X}^*_b)^2$$


---

### with for-loop in R

```{r, eval=TRUE, echo=TRUE}
set.seed(123)
n=101
original=rnorm(n)
median(original)
B=1000
estimator=rep(NA,B)
for (b in 1:B)
  {
  thisboot=sample(x=original,size=n,replace=TRUE)
  estimator[b]=median(thisboot)
}
sd(estimator)
```

---

```{r}
ggplot(data=data.frame(x=estimator),aes(x=x))+
  geom_density()

```

---

### using built in `boot` function from library  `boot`

```{r, eval=TRUE, echo=TRUE}
library(boot)
set.seed(123)
n=101
original=rnorm(n)
median(original)
summary(original)
```

---

```{r}
boot.median=function(data,index) return(median(data[index]))
B=1000
boot(original,boot.median,R=B)
```

---

## With or without replacement?

In bootstrapping we sample _with replacement_ from our observations.

What if we instead sample _without replacement_?

(Then we would always get the same sample - given that the order of the sample points is not important for our estimator.)

(In permutation testing we sample without replacement to get samples under the null hypothesis - a separate field of research.)

---

### Example: multiple linear regression 

We assume, for observation $i$:
$$Y_i= \beta_0 + \beta_{1}  x_{i1} + \beta_2 x_{i2} + ... + \beta_p x_{ip} + \varepsilon_i,$$ where $i=1,2,...,n$. The model can be written in matrix form:
$${\bf Y}={\bf X} \boldsymbol{\beta}+{\boldsymbol{\varepsilon}}.$$

The least squares estimator: $\hat{\boldsymbol\beta}=({\bf X}^T{\bf X})^{-1} {\bf X}^T {\bf Y}$
has $\text{Cov}(\boldsymbol\beta)=\sigma^2({\bf X}^T{\bf X})^{-1}$.

We will in the recommended exercises look at how to use bootstrapping to estimate the covariance of the estimator. Why is that "needed" if we already know the mathematical formula for the standard deviation? Answer: not needed - but ok example.

We will not do this here - but our bootstrap samples can also be used to make confidence intervals for the regression coefficients or prediction intervals for new observations. This means that we do not have to rely on assuming that the error terms are normally distributed!

---

## Bagging 

Bagging is a special case of _ensemble methods_.

In Module 8 we will look at bagging, which is built on bootstrapping the the fact that it is possible to reduce the variance of a prediction by taking the average of many model fits.

* Assume that we have $B$ different predictors $X_1, X_2, \ldots, X_B$. We have built them on $B$ different bootstrap samples.
* All are predictors for some parameter $\mu$ and that all have some unknown variance $\sigma^2$. 
* We then decide that we want to use all the predictors together - equally weighted - and make
$\bar{X}=\frac{1}{n} \sum_{i=1}^n X_i$, which we often use to predict $\mu$.

---

We can therefore obtain a new model (our average of the individual models) that has a smaller variance than each of the individual model because

$$\text{Var}(\bar{X})=\frac{\sigma^2}{n}$$

Since this averaged predictor has smaller variance than each of the predictors we would assum that this is a more accurate prediction. However, the interpretation of this bagged prediction might be harder then for the separate predictors.

Models that have poor prediction ability (as we may see can happen with regression and classification trees) might benefit greatly from bagging. More in Module 8.


---

# <a id="recexboot"> Recommended exercises on bootstrapping</a>

**Problem 1**: We will calculate the probability that a given observation in our original sample is part of a bootstrap sample. This is useful for us to know in Module 8.

Our sample size is $n$.

a. We draw one observation from our sample. What is the probability of drawing observation $x_i$? And of not drawing observation $x_i$?
b. We make $n$ independent drawing (with replacement). What is the probability of not drawing observation $x_i$ in any of the $n$ drawings? What is then the probability that $x_i$ is in our bootstrap sample (that is, more than 0 times)?
c. When $n$ is large $(1-\frac{1}{n})^n=\frac{1}{\exp(1)}$. Use this to give a numerical value for the probability that a specific observation $x_i$ is in our bootstrap sample.
d. Write a short R code chunk to check your result. (Hint: An example on how to this is on page 198 in our ISLR book.) You may also study the result in c. - how fast this happens as a function of $n$.

```{r,echo=FALSE,eval=FALSE}
n=100
B=10000
j=1
res=rep(NA,B)
for (b in 1:B) res[b]=(sum(sample(1:n,replace=TRUE)==j)>0)
mean(res)
```

---

**Problem 2**: Explain with words and an algorithm how you would proceed to use bootstrapping to estimate the standard deviation of one of the regression parameters in multiple linear regression. Comment on which assumptions you make for your regression model.

**Problem 3**: Implement your algorithm from 2 both using for-loop and using the `boot` function. Hint: see page 195 of our ISLR book. Use our SLID data set and provide standard errors for the coefficient for age. Compare with the theoretical value $({\bf X}^T{\bf X})^{-1}\hat{\sigma}^2$ that you find in the output from the regression model.

```{r,eval=FALSE}
library(car)
library(boot)
SLID = na.omit(SLID)
n = dim(SLID)[1]
SLID.lm = glm(wages~., data = SLID)
summary(SLID.lm)$coeff[3,2]
```

---

# Summing up

* Use $k=5$ or $10$ fold cross-validation for model selection or assessment.
* Use bootstrapping to estimate the standard deviation of an estimator.

---

## Plan for interactive lecture 
(Wednesday 10.15-12.00 and Friday 12.15-14.00)

1.  Enter - participate in interactive lecture IL (with lecturer) or supervision of CompEx1 (TAs).
If supervision only - allocated table - or participate in IL and take breaks to ask for supervision of CompEx1.

2. If IL - answer: name+study programme+year+"I found CompEx1 ok or difficult"
Then groups are formed (by lecuturer) as homogeneous as possible (so, not random).

3. 10.15/12.15: Presentation round in groups: name+study programme+year+previous background in statistics+level of expertise in R+favorite hobby!
4. 10.20/12.20: Introduction to problems on cross-validation - work with problems 1-3 (4): [Recommended exercises on cross-validation](#recexcv)
5. 10.50/12.50: Summing up problems 1-2 (3) by lecturer. 
6. 11.00/13.00: Break - with fruits.

---

7. Maybe new groups, maybe not (we evaluate if some of the groups do not work well). If changes, need new presentation round.
8. 11.15/13.15: Introduction to problems on bootstrapping - work with problems 1-3: [Recommended exercises on bootstrapping](#recexboot)
9. 11.40/13.40: Summing up problems 1-2 (3) by lecturer. 
10. 11.45/13.45: Team Kahoot! (in ghostmode on Friday - to beat the Wednesday people).
11. 12.00/14.00: The end:-)

---

# <a id="further"> Further reading </a>

* [Videoes on YouTube by the authors of ISL, Chapter 5](https://www.youtube.com/playlist?list=PL5-da3qGB5IA6E6ZNXu7dp89_uv8yocmf), and corresponding [slides](https://lagunita.stanford.edu/c4x/HumanitiesScience/StatLearning/asset/cv_boot.pdf)
* [Solutions to exercises in the book, chapter 5](https://rstudio-pubs-static.s3.amazonaws.com/65561_43c0eaaa8565414eae333b47038f716c.html)

---

# <a id="Rpackages"> R packages</a>

\footnotesize
```{r,eval=FALSE}
# packages to install before knitting this R Markdown file
# to knit the Rmd
install.packages("knitr")
install.packages("rmarkdown")
# cool layout for the Rmd
install.packages("prettydoc") # alternative to github
#plotting
install.packages("ggplot2") # cool plotting
install.packages("ggpubr") # for many ggplots
#datasets
install.packages("ElemStatLearn") 
install.packages("ISLR")
# cross-validation and bootstrapping
install.packages("boot")
```


