<\body> <\output> \; R : Copyright 2003, The R Foundation for Statistical Computing Version 1.8.1 \ (2003-11-21), ISBN 3-900051-00-3 \; R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. \; R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R in publications. \; Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for a HTML browser interface to help. Type 'q()' to quit R. \; --------------------------------------------------------- TeXmacs interface. Use v() to insert current graphic into TeXmacs buffer. Use q("yes") to quit and save. --------------------------------------------------------- +y> ||>|1>||>|||>|>||>>> <\input|>> > 3 <\output> 3 [1] 3 <\input|>> > 4 <\output> [1] 4 <\input| >> \; > The principle of the bootstrap method is that the data the we gathered can represent the external world. In order to get a distribution of possible values that we could have measured, we re-sample from the data. Let us say that we have a population of 2000 people, whose heights are the following: <\input|>> > population=c(rnorm(1500,mean=180,sd=10),rnorm(500,mean=140,sd=10)) <\input| >> hist(population,n=30);v() <\output> |ps>||||||> <\input| >> \; > Now let us assume that we measure the heights of 20 randomly chosen people exactly: <\input|>> > measurements=sample(population,20) <\input| >> measurements <\output> \ [1] 184.1049 186.0301 162.4606 176.1285 133.3787 135.8835 142.8405 175.5389 \ [9] 117.5551 187.8455 203.9851 141.7764 179.8468 176.3813 147.5747 119.2880 [17] 188.9668 173.2810 173.1390 184.8293 <\input| >> \; > We want to know the median of the distribution. To approximante that, we measure the median of the measurements: <\input|>> > mean(measurements) <\output> [1] 164.5417 <\input| >> @@@ <\output> [1] 164.5417 <\input| >> \; > We know now the mdeian of the mesurements. How wrong are we? If we had the whole population at our disposal, we could do the process of sampling many times, and see the distribution of the median: <\input|>> > pop.samples=t(sapply(1:10000,function(i) sample(population,20) )) <\input| >> dim(pop.samples) <\output> [1] 10000 \ \ \ 20 <\input| >> pop.samples.medians=apply(pop.samples,1,median) <\input| >> hist(pop.samples.medians,n=30);v() <\output> |ps>||||||> <\input| >> median(population) <\output> [1] 176.1048 <\input| >> 2 <\output> [1] 2 <\input| >> \; \; > We can see that we're likely to be about 10cm off the real median. But we don't have the real population. We only have our samples. To estimate the distribution of medians, we sample from the sample, instead of from the real population: <\input| >> sample.samples=t(sapply(1:10000,function(i) sample(measurements,20,rep=T) )) <\input| >> sample.samples.medians=apply(sample.samples,1,median) <\input| >> layout( matrix(1:2,2,1));par(cex=0.7);hist(pop.samples.medians,n=30,xlim=c(135,190) );hist(sample.samples.medians,n=30,xlim=c(135,190));v();layout(1) <\input|>> \; <\output> |ps>||||||> <\input| >> \; > As you see, we do not get the same distribution! The first distribution was created by sampling from the real population. The second distribution was sampled from our sample. Instead of using the sample as our model for the world, we can base a model on the sample, and then sample from it. For example, let us look at our sample. <\input| >> hist(measurements);rug(measurements);v() <\output> |ps>||||||> <\input| >> \; > It obviously does not look normal - it has too much of a tail on the left. We could still try to fit the best normal distrubution we can find, and then bootstrap. <\input| >> a=measurements-mean(measurements) <\input| >> a=a/sd(a) <\input| >> ks.test(a,pnorm) <\output> \; \ \ \ \ \ \ \ \ One-sample Kolmogorov-Smirnov test \; data: \ a\ D = 0.2335, p-value = 0.2254 alternative hypothesis: two.sided\ <\input| >> \; > You can see that because of the small sample size, the KS test is actually not significant. <\input| >> mm=mean(measurements);mm <\output> [1] 164.5417 <\input| >> msd=sd(measurements);msd <\output> [1] 25.19489 <\input| >> parametric=matrix( rnorm(10000*20,mean=mm,sd=msd), 10000, 20) <\input| >> parametric.medians=apply(parametric,1,median) <\input| >> layout( matrix(1:3,3,1));par(cex=0.7);hist(pop.samples.medians,n=30,xlim=c(135,190) );hist(sample.samples.medians,n=30,xlim=c(135,190) );hist(parametric.medians,n=30,xlim=c(135,190));v(width=5,height=5);layout(1) <\input|>> \; <\output> |ps>||||||> <\input| >> sd(sample.samples.medians) <\output> [1] 8.75055 <\input| >> sd(parametric.medians) <\output> [1] 6.861419 <\input| >> sd(pop.samples.medians) <\output> [1] 4.687592 <\input| >> \; > The Jackkinfe method was invented before the bootstrap method, and is very similar. Instead of taking random samples from the distribution, we drop each of our sample points. <\input| >> jack= t( \ sapply(1:length(measurements),function(i) measurements[-i] ) ) <\input| >> jack[1:4,1:4] <\output> \ \ \ \ \ \ \ \ \ [,1] \ \ \ \ [,2] \ \ \ \ [,3] \ \ \ \ [,4] [1,] 186.0301 162.4606 176.1285 133.3787 [2,] 184.1049 162.4606 176.1285 133.3787 [3,] 184.1049 186.0301 176.1285 133.3787 [4,] 184.1049 186.0301 162.4606 133.3787 <\input| >> jack.medians=apply( jack, 1, median ) <\input| >> layout( matrix(1:2,2,1));par(cex=0.7);hist(jack.medians,xlim=c(135,190) );hist(sample.samples.medians,n=30,xlim=c(135,190));v();layout(1) <\input|>> \; <\output> |ps>||||||> <\input| >> \; > We can see that the jackknife estimate of the error is too small. In fact it is too small by a factor of approxmimately >, or /n>>. <\input| >> sd(jack.medians) <\output> [1] 1.158305 <\input| >> sd( jack.medians)*sqrt(20) <\output> [1] 5.180097 \; > It is actually clear that the median will move between 1 2 or 3 points at most: <\input| >> median( 1:3) <\output> [1] 2 <\input| >> median(1:2) <\output> [1] 1.5 <\input| >> median(c(1,3)) <\output> [1] 2 <\input| >> median( 2:3)\ <\output> [1] 2.5 <\input| >> \; > And as we see in the above example, we get just 2 points.\ <\input| >> table( jack.medians ) <\output> jack.medians 173.280953077412 175.538905170459\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 10 \ \ \ \ \ \ \ \ \ \ \ \ \ \ 10\ <\input| >> \; > Because the median has this property that it ``jumps'', the jackknife is a bad choice in this case.\ Let us take another example: <\input| >> a=read.table("/home/dirk/data/regression.txt",head=T) <\input| >> plot(growth~tannin,data=a);v() <\output> |ps>||||||> <\input| >> reg=lm(growth~tannin,data=a) <\input| >> reg <\output> \; Call: lm(formula = growth ~ tannin, data = a) \; Coefficients: (Intercept) \ \ \ \ \ \ tannin \ \ \ \ \ \ 11.756 \ \ \ \ \ \ -1.217 \ <\input| >> names(reg) <\output> \ [1] "coefficients" \ "residuals" \ \ \ \ "effects" \ \ \ \ \ \ "rank" \ \ \ \ \ \ \ \ \ [5] "fitted.values" "assign" \ \ \ \ \ \ \ "qr" \ \ \ \ \ \ \ \ \ \ \ "df.residual" \ \ [9] "xlevels" \ \ \ \ \ \ "call" \ \ \ \ \ \ \ \ \ "terms" \ \ \ \ \ \ \ \ "model" \ \ \ \ \ \ \ <\input| >> reg$coefficients[2] <\output> \ \ \ tannin\ -1.216667\ <\input| >> \; > Let us say that we would like to know how exact this estimate of the slope is First, let us write a function that given growth and tannin tells us the slope: <\input| >> slope=function(gr,ta){reg=lm(gr~ta);reg$coeff[2]} <\input| >> slope(a$gro,a$tan) <\output> \ \ \ \ \ \ \ ta\ -1.216667\ <\input| >> \; > Now let us write a function that given idexes, calculates the slope for those points of data from a: <\input| >> index.slope=function(i) slope( a$gr[i], a$tan[i] ) <\input| >> a <\output> \ \ growth tannin 1 \ \ \ \ 12 \ \ \ \ \ 0 2 \ \ \ \ 10 \ \ \ \ \ 1 3 \ \ \ \ \ 8 \ \ \ \ \ 2 4 \ \ \ \ 11 \ \ \ \ \ 3 5 \ \ \ \ \ 6 \ \ \ \ \ 4 6 \ \ \ \ \ 7 \ \ \ \ \ 5 7 \ \ \ \ \ 2 \ \ \ \ \ 6 8 \ \ \ \ \ 3 \ \ \ \ \ 7 9 \ \ \ \ \ 3 \ \ \ \ \ 8 <\input| >> index.slope(1:9) <\output> \ \ \ \ \ \ \ ta\ -1.216667\ <\input| >> index.slope(1:3) <\output> ta\ -2\ <\input| >> index.slope(-1) <\output> \ \ \ \ \ \ \ ta\ -1.190476\ <\input| >> \; > Now the jackknife is very simple: <\input| >> jack.slopes=sapply(-(1:9), index.slope) <\input| >> jack.slopes <\output> \ \ \ \ \ \ \ ta \ \ \ \ \ \ \ ta \ \ \ \ \ \ \ ta \ \ \ \ \ \ \ ta \ \ \ \ \ \ \ ta \ \ \ \ \ \ \ ta \ \ \ \ \ \ \ ta \ \ \ \ \ \ \ ta\ -1.190476 -1.253133 -1.270270 -1.161359 -1.216667 -1.242038 -1.117117 -1.200501\ \ \ \ \ \ \ \ ta\ -1.321429\ <\input| >> \; > Bootstrap is a bit more complicated: <\input| >> boot=t( sapply(1:200,function(i) sample(1:9,rep=T) ) ) <\input| >> boot.slopes=apply(boot,1,index.slope) <\input| >> layout(matrix(1:2,2,1));par(cex=0.7);hist(boot.slopes,prob=T,xlim=c(-3,-0.5)) <\input| >> hist(jack.slopes,prob=T,xlim=c(-3,-0.5));v();layout(1) <\output> |ps>||||||> <\input| >> sd(boot.slopes) <\output> [1] 0.1943785 <\input| >> sd(jack.slopes) <\output> [1] 0.06090926 <\input| >> sd(jack.slopes)*sqrt(9) <\output> [1] 0.1827278 <\input| >> \; > A general rule of thumb says that to estimate the error, 200 bootstrap samples are enough. To estimate 95% confidence intervals, one needs around 1000. <\input| >> boot=t( sapply(1:1000,function(i) sample(1:9,rep=T) ) ) <\input| >> boot.slopes=apply(boot,1,index.slope) <\input| >> quantile(boot.slopes,0.025) <\output> \ \ \ \ \ 2.5%\ -1.609121\ <\input| >> quantile(boot.slopes,0.975) <\output> \ \ \ \ \ 97.5%\ -0.9084386\ <\input| >> \; > There are actually more acurate methods for calculating confidence intervals using quantiles. It can be done with the function bcanon in the library bootstrap: <\input| >> library(bootstrap) <\input| >> bcanon(1:9,1000,index.slope) <\output> $confpoints \ \ \ \ \ alpha \ bca point [1,] 0.025 -1.6774194 [2,] 0.050 -1.5909091 [3,] 0.100 -1.4782609 [4,] 0.160 -1.4000000 [5,] 0.840 -1.0920330 [6,] 0.900 -1.0520833 [7,] 0.950 -1.0061728 [8,] 0.975 -0.9530201 \; $z0 [1] -0.06521854 \; $acc [1] -0.001202818 \; $u [1] -1.190476 -1.253133 -1.270270 -1.161359 -1.216667 -1.242038 -1.117117 [8] -1.200501 -1.321429 \; $call bcanon(x = 1:9, nboot = 1000, theta = index.slope) \; Warning message:\ multi-argument returns are deprecated in: return(confpoints, z0, acc, u, call = call)\ <\input| >> \; > Permutation tests are very similar to bootstraps. They test if two distributions are the same: <\input| >> x=rnorm(10) <\input| >> y=rnorm(20)+1 <\input| >> mean(x)-mean(y) <\output> [1] -0.7245577 <\input| >> \; > In a permutation test, we permute the data, and ask how often it has such a difference between the samples: <\input| >> perm=t( sapply( 1:1000, function(i) sample( c(x,y) ) ) ) <\input| >> perm.diff = apply(perm, 1, function(xy){ mean(xy[1:10])-mean(xy[11:30]) } ) <\input| >> sum( abs(perm.diff) \ abs((mean(x)-mean(y))) ) <\output> [1] 36 <\input| >> 36/1000 <\output> [1] 0.036 <\input| >> \; > So this permutation test sees the data as significantly different <\input| >> t.test(x,y)$p.value <\output> [1] 0.07868077 <\input| >> \; > The permutation test tests the hypothesis that the two distributions are the same. With a bootstrap test, we could do a very similar test, except that we sample . \; <\initial> <\collection> <\references> <\collection> |?>> |?>> |?>> |?>> |?>> |?>> <\auxiliary> <\collection> <\associate|toc> Lecture 7 - Bootstrap and Jackknife |Example: Median>