---
title: "Bootstrap sampling and Bootstrap Confidence Intervals (Take 2)"
author: ''
output:
  html_document:
    toc: yes
  pdf_document: null
  pdf_presentation: null
  slidy_presentation:
    css: slidy.css
    widescreen: yes
  word_document: null
---



# review of the bootstrap receipe

This is slighlty updated version of our previous lecture on the bootstrap.

The goal is to review the concepts and practice coding.


## Receipe

<ul>
<li>Formulate your problem
<li> Get (or simulate) the 'real' dataset $D$ you want to analyze
<li>Initialize an empty list (this will store your bootstrap distribution)
<li>START THE BOOTSTRAP REPLICATES LOOP (repeat $B$ times)
<li>generate a new bootstrap sample (i.e. a new dataset, resampled from $D$)
<li>Calculate the statistics of interest: $S$ for this sample
<li>Add this estimated stats of interest to your list
<li>END OF BOOTSTRAP REPLICATES LOOP
<li>Analyze the bootstrap distribution (your list): plot, ci, p-values etc.
</ul>



To generate a synthetic  dataset is easy in R:

What is the distribution I sample from here? What are the parameters of this distribution? How many sample do I request? Why do I use the function set.seed?


```{r}
set.seed(87612345)
data=rnorm(20,mean=4,sd=.75)
```

Same questions:


```{r}
set.seed(117426)
data=rchisq(20,1,nc=3)
```

And if you are not familiar with the chi-square distribution, read:

<a href="https://astrostatistics.psu.edu/su07/R/html/stats/html/Chisquare.html">Chi-square in R</a>

We can enerate a new bootstrap sample ('fake new data') by resampling the data with replacement. Note how I choose to generate a sample (called bdata here) of the same length as the original sample.

```{r}
bdata=sample(data,20,replace=T)
```

Now I can calculate statistics of interest, which I note with the generic letter $S$, for this sample. The choice of $S$ depends on what I want to estimate/ my problem.


## Review some common statistics



<b>Calculate statistics of interest: $S$ for this sample</b>

Can you define a statistic?

A statistic is a quantity that I can calculate just based on the data, i.e. a given
function of the sample value.

We can note this $S=g(x)$.

Some examples of choices for $S$:

<ul>
<li>Sample mean (what is the function $g$?)
<li>Sample variance
<li>Sample median
<li>Sample quantiles
<li>Sample standard deviation
<li>Sample min
<li>Sample max
<li>Sample order statistics
<li>Sample gap between top and next top value
<li>Sample number of dogs with more than $k$ fleas on them (salmon -lice)
<li>Topology of sampled tree
<li>Degree distribution of sampled graphs
<li>Whatever you want to calculate from the data
</ul>


SAMPLE STATISTICS ARE RANDOM VARIABLES AND HAVE THEIR OWN DISTRIBUTION

We can say some things about sample distributions EVEN IF we dont know the distribution the data comes from

EXAMPLES
<ul>
<li>mean and variance of S=sample mean

What is the variance of the sample mean statistic?
$$Var(\bar{x}) = Var(\frac{1}{n} \sum_{i=1}^n x_i) = Var(x_i)/n = V/n$$
   

IF we know which distribution data is drawn for, we can theoretically predict the distribution of sample statistics

Example: if $x_i \sim N(\mu,\sigma^2)$, then what is the distribution of $bar(x}$?

<li>mean and variance of S=sample variance

</ul>

Mean and medians ?

difference between mean and median?

```{r}
v=c(1,2,3,4)
mean(v)
median(v)
v=c(1,2,3,4000)
mean(v)
median(v)
#v = c(1,2,90,10000)
#median(v)
# 46
#mean(v)
v = c(1,2,90,100000000)
median(v)
```

Try to memorize the mean and variance of common distributions

$$Var(\chi^2(k) = 2 k$$


Quantiles


```{r}
n=21
v = 1:n
v
quantile(v,probs=c(.3,.5,.7))
```


```{r}
w=v^2
w
quantile(w,probs=c(.3,.5,.7))
```


```{r}
w=abs(cos(v))
w
x <-1:21
quantile(w,probs=c(.3,.5,.7))
plot(x,w)
lines(x,w)
```


## building the bootstrap distribution

Add the stat of interest to your list


```{r}
# outside the bootstrap loop I will initialize an empty list
#
myListOfStats = NULL
# this line of code will append the new statistic value myStat to the existing list
# myListOfStats
# myListOfStats=c(myListOfStats,myStat)
```

Analyze the bootstrap distribution


Can you now see how the following code puts all these elements together?

```{r}
set.seed(117426)
data=rchisq(20,8,ncp=5)
# theoretical mean is k + lambda  df+ncp=13
mystatchoice = function(ichoice,mydata){
return(mean(mydata))
}

myListOfStats = NULL
for (i in 1:1000){
mybsample   =  sample(data,20,replace=T)
#mystat  =  mystatchoice(1, mybsample)
mystat = mean(mybsample)
myListOfStats = c(myListOfStats,mystat)
}

myListOfStats[1:20]
mean(myListOfStats)
```



Histograms are a way to see the distribution:

```{r}
hist(myListOfStats,nclass=20)
```



Quantiles are cutpoints that cut the values of the random variable (response) into zones of well-defined areas (i.e. probabilities). Here I ask which cut-point $c_1,c_2$ have the probabilities $P(X \leq c_1) = 0.025$ and $P(X \leq c_2) = 0.975$:

```{r}
# this will return the two cutpoints c_1 and c_2
# they are called the 2.5% percentile and 97.5% percentile of the empirical distribution
# provided in the vector myListOfStats (this vector stores random values of your random variable / response)
quantile(myListOfStats,probs=c(.025,.975))
```


Here I recall how to define a matrix and how to calculate the median of each row

```{r}
#apply
M=matrix(1:12,byrow=T,ncol=4)
M
apply(M,1,median)
```
Or the mean of each column

```{r}
#apply
M=matrix(1:12,byrow=T,ncol=4)
M
apply(M,2,mean)
```


Or the variance of each row

```{r}
#apply
M=matrix(1:12,byrow=T,ncol=4)
M
apply(M,2,var)
```



# Bootstrap sample

Suppose that a data set consists of the $n$ items $x_1, x_2, \ldots , x_n$.

A *bootstrap sample* is a sample of size $n$, drawn with $replacement$, from
the original sample.  Because the sampling is with replacement, the sample points
can be repeated in the bootstrap sample.  Here's an example.

```{r}
data=round(rnorm(6),1)
data
```

There are six items in the data set.  Now let's draw 10 bootstrap samples, each
of size 6, with replacement.  The R command to draw a sample of size $n$ from $x$,
with replacment, is $sample(x,n,replace=T)$.    We don't have to know $n$ ahead
of time, as we can just use $sample(x,length(x),replace=T)$.
Let's use $Nboot$ to denote the number of bootstrap samples.


```{r}
Nboot=10
for (i in 1:Nboot){
  bsample=sample(data,length(data),replace=T)
  print(bsample)
  }
```

#Bootstrap distribution

Generally we are interested in some statistic, which is nothing more than a function
of the sample values, say the mean, or a specified percentile. We can calculate
the statistic for each of our bootstrap samples, and thereby construct the  "bootstrap
distribution" of the  statistic.  The  following examples show the bootstrap distribution
of the mean and of the 90'th percentile, using 500 bootstrap samples.  Let's start with a sample of
size 15 from a $chi^2$ distribution with 2 degrees of freedom.


```{r}
Nboot=500
n=15; df=2
x=rchisq(n,df) #draw sample from a chi-squared distribution
means=NULL
p90=NULL
for (i in 1:Nboot){
  bsample=sample(x,n,replace=T)
  means=c(means,mean(bsample)) #calculate the mean 
  p90=c(p90,quantile(bsample,.9)) #calulate 90'th percentile
  }
hist(means,nclass=20,main="bootstrap distribution of the mean")
```

```{r}
hist(p90,nclass=20,main="bootstrap distribution of the 90'th percentile")
```

##Bootstrap confidence intervals

A common use of the bootstrap is to calculate a confidence interval
which doesn't depend on underlying distributional assumptions.
The $percentile method$ creates a $100(1-\alpha)$% confidence interval
using the upper and lower $\alpha/2$'th percentiles of the bootstrap
distribution.  Where **x** is some data, the upper and lower
$\alpha$'th percentiles of **x** are calculated in R using
using $quantile(x,c(\alpha/2,1-\alpha/2))$.

### Example:

draw a sample of size 10 from $N(2,1)$, and calculate
confidence intervals for the underlying population mean, using 
both the usual $t$-interval, and a bootstrap interval using 500 bootstrap
samples.


```{r}
nboot=500
n=10; mu=2; sd=1
alpha=.05 #for 95% interval
data=rnorm(n,mu,sd) #simulate a data set
# get the bootstrap distribution of the mean
 bmean=NULL
 for (i in 1:nboot){
  bsamp=sample(data,n,replace=T) #take a bootstrap sample from data
  bmean=c(bmean,mean(bsamp)) #accumlate the mean
  }
bint=quantile(bmean,c(alpha/2,1-alpha/2)) #calculate percentile interval
cat("bootstrap interval is ",bint)
```

As was indicated in the notes on the $t-interval$, 
the actual coverage of the $t-interval$ will be correct when the
the assumption that the underlying sample is from a normal 
distribution holds. Also, its empirical coverage is generally pretty
close to its nominal coverage if the underlying distribution is approximately
normal. (This is because the Central Limit Theorem says that sample means will typically
be close to normally distributed.)

The bootstrap percentile interval$ makes no assumption about the shape of the underlying distribution.  

## Assessing the empirical coverage of the bootstrap confidence interval

The following function bconflevel simulates $Nbatch$ samples of size $n$ from 
a distribution specifed by $FUN$.  Arguments to $FUN$ are passed in using the $...$ notation.
For each sample, a $100(1-\alpha)$% bootstrap confidence interval is constructed, and the function
counts the proportion of the samples where the confidence interval contains the true, but
assumed unknown, parameter. This proportion is the empirical coverage.

```{r}
bconflevel=function(Nbatch,Nboot=100,alpha=.05,n=20,truemean,FUN,...){
df=n-1 #degrees of freedom for t
ta2=qt(1-alpha/2,df) # 100(1-alpha/2)'th percentile of t
count=0
print(truemean)
for (i in 1:Nbatch){
data=FUN(n,...)
#generate the bootstrap distribution for the mean
  bmeans=NULL
  for (j in 1:Nboot){
    datab=sample(data,size=length(data),replace=T) #bootstrap sample
    bmeans=c(bmeans,mean(datab))}
  quants=quantile(bmeans,probs=c(alpha/2,1-alpha/2)) # upper and lower percentage points
#if the interval contains the true mean, increment the count
if(quants[1]<truemean && quants[2]>truemean)count=count+1
}
emplevel=count/Nbatch   #empirical level of the t-interval
paste("Nominal level = ",1-alpha," Empirical level = ",emplevel)
}
```
### Example:  sampling from a normal population.

$Nbatch=500$ samples of size 10 are drawn from a normal distribution
with mean 4 and standard deviation .3.  The empirical coverage
of the bootstrap interval is reported.
For the bootstrap interval, 1000 bootstrap samples are used for
each sample, so computing time for the bootstrap interval
is roughly 1000 times that for calculating a t-interval.

```{r}
bconflevel(Nbatch=500,Nboot=1000,n=10,truemean=4,FUN=rnorm,mean=4,sd=.3)
```

The empirical coverage here will generally be fairly 
close to the nominal coverage because we are sampling from a normal population.  Differences between the empirical
and the nominal level will depend on a number of things, including
the number of simulation batches, the sample size, the generated data, and the percentile method itself.
Keep in mind, though, that the percentile interval doesn't use any information about the underlying distribution.


### Example: sampling from a $\chi^2$ distribution

Here we are sampling from a distribution which is very non-normal.

```{r}
bconflevel(500,Nboot=1000,n=10,truemean=2,FUN=rchisq,df=2) #mean of chisq_2 is 2
```

## Boostrap interval for percentiles

There is no analytical form for a confidence interval 
for a population percentile. This differs from the situation of the population
mean, where the t-interval can be used.  However, it's easy to make a bootstrap interval for a percentile.

The following function carries out a simulation study assessing the actual coverage of the nominal 95\% bootstrap interval
for the 75'th percentile,  when the underlying distribution is $\chi^2_2$.  For this distribution, the 75'th 
percentile is

```{r}
qchisq(.75,2)
```

```{r}
bconflevel=function(Nbatch=500,Nboot=100,alpha=.05,n=20,FUN,...){
pct=.75; df=2
truep=qchisq(pct,2)
count=0
for (i in 1:Nbatch){
data=FUN(n,...) #simulate data
n=length(data)
#generate the bootstrap distribution for the 90th percentile
  bp90=NULL
  for (j in 1:Nboot){
    datab=sample(data,size=n,replace=T) #bootstrap sample
    bp90=c(bp90,quantile(datab,pct))} #accumulate the bootstrapped percentiles
    quants=quantile(bp90,probs=c(alpha/2,1-alpha/2)) # upper and lower percentage points
#if the interval contains the true percentile, increment the count
if(quants[1]<truep && quants[2]>truep)count=count+1
}
emplevel=count/Nbatch   #empirical level of the t-interval
paste("Nominal level = ",1-alpha," Empirical level = ",emplevel)
}
```

The following prints the empirical coverage of the percentile interval
with 5 replicate simulations. Each simulation uses 500 simulation batches,
but only 100 bootstrap samples for each generated data set, in order to
not be too computationally expensive.   This shows the variability in the bootstrap coverage.

```{r}
#run 5 replicate simulations, each with 500 batches, but only 100 bootstrap samples
for (j in 1:5){ 
print(bconflevel(Nbatch=500,Nboot=100,n=30,FUN=rchisq,df=2))}
```

The coverage will typically differ from the nominal level, here 95%,
but keep in mind that there is no analytical method to make a confidence
interval for a percentile, so the bootstrap has given a fairly simple way
to do something that we couldn't do otherwise.

In practice, one would want to use considerably more than 100 bootstrap
samples.  1000 might be reasonable.


centile interval
with 5 replicate simulations. Each simulation uses 500 simulation batches,
but only 100 bootstrap samples for each generated data set, in order to
not be too computationally expensive.   This shows the variability in the bootstrap coverage.

```{r}
#run 5 replicate simulations, each with 500 batches, but only 100 bootstrap samples
for (j in 1:5){ 
print(bconflevel(Nbatch=500,Nboot=100,n=30,FUN=rchisq,df=2))}
```

The coverage will typically differ from the nominal level, here 95%,
but keep in mind that there is no analytical method to make a confidence
interval for a percentile, so the bootstrap has given a fairly simple way
to do something that we couldn't do otherwise.

In practice, one would want to use considerably more than 100 bootstrap
samples.  1000 might be reasonable.


# try other distributions

## wrap things in a function?


```{r}

mymodel <- 1

Nboot=500
n=100; 

if(mymodel == 1){
mu <- 2
sd <- 0.25
x=rnorm(n,mu,sd) #draw sample from a normal distribution
}

if(mymodel == 2){
df=2
x=rchisq(n,df) #draw sample from a chi-squared distribution
}

means=NULL
p90=NULL
for (i in 1:Nboot){
  bsample=sample(x,n,replace = TRUE)
  means=c(means,mean(bsample)) #calculate the mean 
  p90=c(p90,quantile(bsample,.9)) #calulate 90'th percentile
  }
hist(means,nclass=20,main="bootstrap distribution of the mean")
```


```{r}
dumm <- function(nb=nb){
v <- array()


for (i in 1:nb){
v[i] <- mean(rchisq(i,df))
}
x <- 1:nb
y <- v
plot(x,y)
}
```


```{r}
dumm(nb=5)
```

```{r}
dumm(nb=50)
```



# Examples of bootstrapping


## Bootstrapping the standard error of the slope.


```{r}
x = c(1,1.5,2,3,4,4.5,5,5.5,6,6.5,7,8,9,10,11,12,13,14,15)
y = c(6.3,11.1,20,24,26.1,30,33.8,34.0,38.1,39.9,42,46.1,53.1,52,52.5,48,42.8,27.8,21.9)
n=length(y)
plot(x,y,xlab="hardwood %",ylab="tensile strength")
```

We would like to use $x$ to predict $y$.

   + Let's start with a **simple linear regression** model

$$y_i = \beta_0 + \beta_1 x_i + \epsilon_i$$


There is a builtin procedure "boot" which can be used to estimate
the standard error of any statistic.  In order to use it, you need
to create a function which returns the statistic you're interesedt in.
For the regression slope,

```{r}
# this is an extracting function that extracts the estimated slope  of a linear regression
b1=function(data,index){
  x=data$x[index]
  y=data$y[index]
  lmout=lm(y~x)
  b1=coef(lmout)[2]
 # this will return the estimated slope of the linear regression
  return(b1)}
```

"index" is a subset of 1:n which will be passed in by the boot command.


```{r}
library(boot)
data=data.frame(x=x,y=y)
bout <- boot(data,b1,R=1000)
# bout is an object that contains several attributes
# one of them called t contains the value (b1) returned by your extracting function
# bout$t therefore contains the values of the estimated slope for each bootstrap sample
# this is a vector of length R=1000 here -number of bootstrap sample-
# that contains the bootstrap distribution of the estimated slope
hist(bout)
```

The original estimate ($\hat \beta_1$), the standard error (estimated variance of $\hat \beta_1$)  and something called the bias of $\hat \beta_1$ are returned.


## Bootstrapping in population genetics

Once upon a time, there was ... linanthus parryae.

<img src="../images/blue2.png">

<img src="https://upload.wikimedia.org/wikipedia/commons/3/39/Yucca_brevifolia_bushy.jpg">

<a href="https://www.rwongphoto.com/gallery/mojave-desert-pictures/">More pictures of the Mojave desert</a>


<img src="../images/blue1.png">



$$GD=p(1-p)$$


## Bootstrapping relatedness coefficients



Microsatellites markers are used in forensic analysis


PCR  + electrophoresis on agarose gel

<img width=500 src="https://www.researchgate.net/profile/Prescott_Deininger2/publication/255797983/figure/fig1/AS:652244831707140@1532518694441/Agarose-gel-electrophoresis-and-codominant-segregation-of-polymorphic-Alu-PCR-products.png">




$$R_{ij}= \sum_l \sum_{i,j} z_{il}z_{jl}$$
## Bootstrapping abundances (species)

coral reef paper

so many others...


# assignment 5 help/notes



what does this code do?

```{r}
myfun=function(x){
if(length(x)==1){return(x)
} else
return(x[1]+myfun(x[-1]))
}

v = c(5,6,7)
myfun(v)

```


$$cox(x)+cos(2x)+cos(3x)+ \ldots$$

$$\sum_{i=1}^n g(i)$$

$$\sum_{i=1}^n (1/i)^p$$


# external links on the bootstrap

<ul>
<li><a href="https://arxiv.org/pdf/1411.5279.pdf">What teachers should know about the bootstrap</a>
<li><a href="https://wwwf.imperial.ac.uk/~agandy/teaching/ltcc/lecture4.pdf">Alex Gandy: Lecture on the bootstrap</a>
</ul>




