<\body> \; Till now we learned about vectors, matrices, and dataframes, and we also leared about classes. Let us look a bit deeper into R. <\input| >> a=c(x=1,y=5) <\input| >> a <\output> x y\ 1 5\ <\input| >> a=c(1:4,10:14) <\input| >> a <\output> [1] \ 1 \ 2 \ 3 \ 4 10 11 12 13 14 <\input| >> a=c(x=1:4,y=10:14) <\input| >> a <\output> x1 x2 x3 x4 y1 y2 y3 y4 y5\ \ 1 \ 2 \ 3 \ 4 10 11 12 13 14\ <\input| >> a=list(x=1:4,y=10:14) <\input| >> a <\output> $x [1] 1 2 3 4 \; $y [1] 10 11 12 13 14 <\input| >> a$x <\output> [1] 1 2 3 4 <\input| >> a$y <\output> [1] 10 11 12 13 14 <\input| >> a[["x"]] <\output> [1] 1 2 3 4 <\input| >> a[["y"]] <\output> [1] 10 11 12 13 14 <\input| >> a=list(x=1:4,y=c("a","b"),z=1:10) <\input| >> a <\output> $x [1] 1 2 3 4 \; $y [1] "a" "b" \; $z \ [1] \ 1 \ 2 \ 3 \ 4 \ 5 \ 6 \ 7 \ 8 \ 9 10 <\input| >> a[1] <\output> $x [1] 1 2 3 4 <\input| >> a[1:2] <\output> $x [1] 1 2 3 4 \; $y [1] "a" "b" <\input| >> a[[1]] <\output> [1] 1 2 3 4 <\input| >> \; > Using will give us a sublist of the list, which is still a list. Using will give us elements of the list. We already learned about , there is also a function called . <\input| >> lapply(1:10,function(x) x^2) <\output> [[1]] [1] 1 \; [[2]] [1] 4 \; [[3]] [1] 9 \; [[4]] [1] 16 \; [[5]] [1] 25 \; [[6]] [1] 36 \; [[7]] [1] 49 \; [[8]] [1] 64 \; [[9]] [1] 81 \; [[10]] [1] 100 <\input| >> a=lapply(1:10,function(x) x^2) <\input| >> a[1:3] <\output> [[1]] [1] 1 \; [[2]] [1] 4 \; [[3]] [1] 9 <\input| >> a[[3]] <\output> [1] 9 <\input| >> a[[1:3]] <\output> Error: recursive indexing failed at level 2 <\input| >> \; > sapply and lapply are very similar, except that sapply will turn its output into a matrix or vector in the end. <\input| >> a=lapply(31:40,c);a <\output> [[1]] [1] 31 \; [[2]] [1] 32 \; [[3]] [1] 33 \; [[4]] [1] 34 \; [[5]] [1] 35 \; [[6]] [1] 36 \; [[7]] [1] 37 \; [[8]] [1] 38 \; [[9]] [1] 39 \; [[10]] [1] 40 <\input| >> names(a)=c("a","b","c","d","e","f","g","h","i","j") <\input| >> a <\output> $a [1] 31 \; $b [1] 32 \; $c [1] 33 \; $d [1] 34 \; $e [1] 35 \; $f [1] 36 \; $g [1] 37 \; $h [1] 38 \; $i [1] 39 \; $j [1] 40 <\input| >> a$d <\output> [1] 34 <\input| >> a["e"] <\output> $e [1] 35 <\input| >> b=lapply(a,function(x) x-20) <\input| >> b <\output> $a [1] 11 \; $b [1] 12 \; $c [1] 13 \; $d [1] 14 \; $e [1] 15 \; $f [1] 16 \; $g [1] 17 \; $h [1] 18 \; $i [1] 19 \; $j [1] 20 <\input| >> sapply(a,function(x) x-20) <\output> \ a \ b \ c \ d \ e \ f \ g \ h \ i \ j\ 11 12 13 14 15 16 17 18 19 20\ <\input| >> \; > How to turn a list into a vector? <\input| >> b <\output> $a [1] 11 \; $b [1] 12 \; $c [1] 13 \; $d [1] 14 \; $e [1] 15 \; $f [1] 16 \; $g [1] 17 \; $h [1] 18 \; $i [1] 19 \; $j [1] 20 <\input| >> unlist(b) <\output> \ a \ b \ c \ d \ e \ f \ g \ h \ i \ j\ 11 12 13 14 15 16 17 18 19 20\ <\input| >> \; > Another extremely usefull function is <\input| >> setwd("~/R-course-2006/lecture8") <\input| >> list.files() <\output> [1] "factorial.txt" \ \ \ \ \ \ \ \ \ "oneway.txt" \ \ \ \ \ \ \ \ \ \ \ \ "Rcourse8_no_output.tm"\ [4] "Rcourse8_no_output.tm~" "Rcourse8.pdf" \ \ \ \ \ \ \ \ \ \ "Rcourse8.tm" \ \ \ \ \ \ \ \ \ \ [7] "Rcourse8.tm~" \ \ \ \ \ \ \ \ \ \ "regression.txt" \ \ \ \ \ \ \ \ "twoway.txt" \ \ \ \ \ \ \ \ \ \ \ <\input| >> fact=read.table("factorial.txt",sep="\\t",head=T) <\input| >> fact <\output> \ \ \ growth diet \ coat 1 \ \ \ \ 6.6 \ \ \ A light 2 \ \ \ \ 7.2 \ \ \ A light 3 \ \ \ \ 6.9 \ \ \ B light 4 \ \ \ \ 8.3 \ \ \ B light 5 \ \ \ \ 7.9 \ \ \ C light 6 \ \ \ \ 9.2 \ \ \ C light 7 \ \ \ \ 8.3 \ \ \ A \ dark 8 \ \ \ \ 8.7 \ \ \ A \ dark 9 \ \ \ \ 8.1 \ \ \ B \ dark 10 \ \ \ 8.5 \ \ \ B \ dark 11 \ \ \ 9.1 \ \ \ C \ dark 12 \ \ \ 9.0 \ \ \ C \ dark <\input| >> tapply(fact$growth,fact$diet,c) <\output> $A [1] 6.6 7.2 8.3 8.7 \; $B [1] 6.9 8.3 8.1 8.5 \; $C [1] 7.9 9.2 9.1 9.0 <\input| >> tapply(fact$growth,fact$coat,c) <\output> $dark [1] 8.3 8.7 8.1 8.5 9.1 9.0 \; $light [1] 6.6 7.2 6.9 8.3 7.9 9.2 <\input| >> tapply(fact$growth,fact$diet,median) <\output> \ \ \ A \ \ \ B \ \ \ C\ 7.75 8.20 9.05\ <\input| >> tapply(fact$growth,fact$coat,c) <\output> $dark [1] 8.3 8.7 8.1 8.5 9.1 9.0 \; $light [1] 6.6 7.2 6.9 8.3 7.9 9.2 <\input| >> tapply(fact$growth,fact$coat,median,simplify=F) <\output> $dark [1] 8.6 \; $light [1] 7.55 <\input| >> tapply(fact$growth,fact$coat,median,simplify=T) <\output> \ dark light\ \ 8.60 \ 7.55\ <\input| >> \; > dataframe are also lists, but special lists. <\input| >> fact <\output> \ \ \ growth diet \ coat 1 \ \ \ \ 6.6 \ \ \ A light 2 \ \ \ \ 7.2 \ \ \ A light 3 \ \ \ \ 6.9 \ \ \ B light 4 \ \ \ \ 8.3 \ \ \ B light 5 \ \ \ \ 7.9 \ \ \ C light 6 \ \ \ \ 9.2 \ \ \ C light 7 \ \ \ \ 8.3 \ \ \ A \ dark 8 \ \ \ \ 8.7 \ \ \ A \ dark 9 \ \ \ \ 8.1 \ \ \ B \ dark 10 \ \ \ 8.5 \ \ \ B \ dark 11 \ \ \ 9.1 \ \ \ C \ dark 12 \ \ \ 9.0 \ \ \ C \ dark <\input| >> fact[["growth"]] <\output> \ [1] 6.6 7.2 6.9 8.3 7.9 9.2 8.3 8.7 8.1 8.5 9.1 9.0 <\input| >> fact[["coat"]] <\output> \ [1] light light light light light light dark \ dark \ dark \ dark \ dark \ dark\ Levels: dark light <\input| >> fact[1:2] <\output> \ \ \ growth diet 1 \ \ \ \ 6.6 \ \ \ A 2 \ \ \ \ 7.2 \ \ \ A 3 \ \ \ \ 6.9 \ \ \ B 4 \ \ \ \ 8.3 \ \ \ B 5 \ \ \ \ 7.9 \ \ \ C 6 \ \ \ \ 9.2 \ \ \ C 7 \ \ \ \ 8.3 \ \ \ A 8 \ \ \ \ 8.7 \ \ \ A 9 \ \ \ \ 8.1 \ \ \ B 10 \ \ \ 8.5 \ \ \ B 11 \ \ \ 9.1 \ \ \ C 12 \ \ \ 9.0 \ \ \ C <\input| >> lapply(fact,c) <\output> $growth \ [1] 6.6 7.2 6.9 8.3 7.9 9.2 8.3 8.7 8.1 8.5 9.1 9.0 \; $diet \ [1] 1 1 2 2 3 3 1 1 2 2 3 3 \; $coat \ [1] 2 2 2 2 2 2 1 1 1 1 1 1 <\input| >> attributes(fact) <\output> $names [1] "growth" "diet" \ \ "coat" \ \; $class [1] "data.frame" \; $row.names \ [1] "1" \ "2" \ "3" \ "4" \ "5" \ "6" \ "7" \ "8" \ "9" \ "10" "11" "12" <\input| >> \; > The of an object can cause R to treat it differently. Two attributes that we already encountered are , and . fact is a list, but because the class is ``data.frame'', it is treated in a special way. <\input| >> a=matrix(1:12,3,4) <\input| >> a <\output> \ \ \ \ \ [,1] [,2] [,3] [,4] [1,] \ \ \ 1 \ \ \ 4 \ \ \ 7 \ \ 10 [2,] \ \ \ 2 \ \ \ 5 \ \ \ 8 \ \ 11 [3,] \ \ \ 3 \ \ \ 6 \ \ \ 9 \ \ 12 <\input| >> attributes(a) <\output> $dim [1] 3 3 <\input| >> \; > all that seperates a vector from a matrix is the attributes. <\input| >> attributes(a)$dim=c() <\input| >> a <\output> \ [1] \ 1 \ 2 \ 3 \ 4 \ 5 \ 6 \ 7 \ 8 \ 9 10 11 12 <\input| >> attributes(a)$dim=c(3,4) <\input| >> a <\output> \ \ \ \ \ [,1] [,2] [,3] [,4] [1,] \ \ \ 1 \ \ \ 4 \ \ \ 7 \ \ 10 [2,] \ \ \ 2 \ \ \ 5 \ \ \ 8 \ \ 11 [3,] \ \ \ 3 \ \ \ 6 \ \ \ 9 \ \ 12 <\input| >> attributes(a)$dim=c(4,3) <\input| >> a <\output> \ \ \ \ \ [,1] [,2] [,3] [1,] \ \ \ 1 \ \ \ 5 \ \ \ 9 [2,] \ \ \ 2 \ \ \ 6 \ \ 10 [3,] \ \ \ 3 \ \ \ 7 \ \ 11 [4,] \ \ \ 4 \ \ \ 8 \ \ 12 <\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] 182.6651 195.6872 182.5375 174.9788 183.1086 164.1614 127.4974 178.6637 \ [9] 184.6356 178.1290 179.3890 170.3763 138.7274 183.9492 182.9451 181.6166 [17] 191.0370 179.0268 124.8710 128.2074 <\input| >> \; > We want to know the median of the distribution. To approximante that, we measure the median of the measurements: <\input|>> > median(measurements) <\output> [1] 179.2079 <\input| >> \; > We know now the median 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| >> \; \; > 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) <\output> Error in library(bootstrap) : there is no package called 'bootstrap' <\input| >> install.packages("bootstrap") <\output> Warning in install.packages("bootstrap") : argument 'lib' is missing: using /home/michael/lib/R \; Warning message: unable to resolve 'cran.r-mirror.de'.\ Warning: unable to access index for repository http://cran.r-mirror.de/src/contrib \; Warning message: no package 'bootstrap' at the repositories in: download.packages(pkgs, destdir = tmpd, available = available, \ <\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| >> \; > \; <\initial> <\collection> <\references> <\collection> > > > > > > > > |?>> |?>> |?>> |?>> |?>> |?>> <\auxiliary> <\collection> <\associate|toc> |math-font-series||Lecture 7 - Bootstrap and Jackknife> |.>>>>|> |Example: Median |.>>>>|> > |Parametric bootstrap |.>>>>|> > |Jackknife |.>>>>|> > |Confidence intervals |.>>>>|> > |Permutation tests |.>>>>|> >