# Bret Larget # October 28, 2010, November 6, 2011 # Function for bootstrap # Function carries out the nonparametric bootstrap # First argument is a single sample # Second argument is function to apply to each sample; default is mean # Third argument is number of bootstrap replicates; default is 10,000 # Fourth argument is alternative; default is "two.sided", # other choices are "less" and "greater". Partial match is okay. boot.test = function(x,fun=mean,B=10000,alternative="two.sided") { n = length(x) boot.sample = apply(matrix(sample(x,n*B,replace=T),nrow=B,ncol=n),1,fun) alternative = match.arg(alternative, c("two.sided","less","greater")) test.stat = fun(x) if (alternative=="two.sided") p.value = sum(abs(boot.sample) >= abs(test.stat))/B else if (alternative=="less") p.value = sum(boot.sample <= test.stat)/B else if (alternative=="greater") p.value = sum(boot.sample >= test.stat)/B return(list(p.value=p.value,test.stat=test.stat,boot.sample=boot.sample)) } # Function for a bootstrap percentile confidence interval # Note that BCa method is more complicated, but better; see library(boot) # Example: # > sockeye = read.table("sockeye.txt", header=T) # > str(sockeye) # > sockeye.boot = boot.ci(sockeye$mass) # > sockeye.boot$ci # > library(lattice) # > densityplot(sockeye.boot$boot.sample,plot.points=F) boot.ci = function(x,fun=mean,B=10000,conf=0.95) { n = length(x) boot.sample = apply(matrix(sample(x,n*B,replace=T),nrow=B,ncol=n),1,fun) if ( conf <= 0 || conf >= 1 ) stop("Must have confidence level between 0 and 1") ci = quantile(boot.sample,c((1-conf)/2,1-(1-conf)/2)) return(list(ci=ci,boot.sample=boot.sample)) } # Function for a two-independent sample permutation test (difference in means). perm.test = function(x,y,B=10000,alternative="two.sided") { n1 = length(x) n2 = length(y) f = function(z,n1,n2) { return(mean(z[1:n1])-mean(z[(n1+1):(n1+n2)])) } out = rep(0,B) z = c(x,y) for( i in 1:B ) out[i] = f(sample(z),n1,n2) test.stat = f(z,n1,n2) alternative = match.arg(alternative, c("two.sided","less","greater")) if (alternative=="two.sided") p.value = sum(abs(out) >= abs(test.stat))/B else if (alternative=="less") p.value = sum(out <= test.stat)/B else if (alternative=="greater") p.value = sum(out >= test.stat)/B return(list(p.value=p.value,test.stat=test.stat,test.sample=out)) }