# Bootstrap test example using mouse data from p. 11 of E&T.

set.seed(321)

survival <- c(94,197,16,38,99,141,23,52,104,146,10,51,30,40,27,46)

obs.diff.means <- mean( survival[1:7] ) - mean( survival[8:16] )

bss <- matrix( sample( survival, 16*40000, replace=T ), nrow=40000 )

diffmeans <- function( x ) { mean( x[1:7] ) - mean( x[8:16] ) }

diff.means <- apply( bss, 1, diffmeans )

# one-sided test p-value (bootstrap based on difference in sample means):
length( diff.means[ diff.means >= obs.diff.means ] )/40000 


# Now I'll do a bootstrap test based on Student's two-sample t
# statistic (as opposed to using the difference in sample means).

newt <- function( x ) {
    numerator <- mean( x[1:7] ) - mean( x[8:16] )
    var.pool <- ( var( x[1:7] )*6 + var( x[8:16] )*8 )/14
    denom <- sqrt( var.pool*( 1/7 + 1/9 ) )
    t <- numerator/denom
    t
  }

# observed Student's t statistic:
obs.new.t <- newt( survival )

new.t <- apply( bss, 1, newt )

# one-sided test p-value (bootstrap based on Student's t):
length( new.t[ new.t >= obs.new.t ] )/40000

 
# Now I'll do a bootstrap test based on Welch's statistic.

welch <- function( x ) {
    numerator <- mean( x[1:7] ) - mean( x[8:16] )
    denom <- sqrt( var( x[1:7] )/7 + var( x[8:16] )/9 )
    t <- numerator/denom
    t
  }

# observed Welch's statistic:
obs.welch <- welch( survival )

treat <- c(94,197,16,38,99,141,23)
control <- c(52,104,146,10,51,30,40,27,46)
 
treat.shift <- treat - mean( treat ) + mean( survival )
control.shift <- control - mean( control ) + mean( survival )

bss.treat <- matrix( sample( treat.shift, 7*40000, replace=T ), nrow=40000 )
bss.control <- matrix( sample( control.shift, 9*40000, replace=T ), nrow=40000 )
bss <- cbind( bss.treat, bss.control )

welch.reps <- apply( bss, 1, welch )

# one-sided test p-value (bootstrap based on Welch's statistic):
length( welch.reps[ welch.reps >= obs.welch ] )/40000

# For comparison purposes, I'll do Welch's test.

t.test( treat, control, alternative="greater", var.equal=F )