Permutation test vs. bootstrap for hypothesis testing

sigdiff.R

nA = 8
meanA = 48
sdA = 3
nB = 10
meanB = 53
sdB = 4

set.seed(42)

nBoot = 100000

datA = rnorm(nA,meanA,sdA)
datB = rnorm(nB,meanB,sdB)

dat = data.frame(obs=c(datA,datB),lab=c(rep("A",nA),rep("B",nB)))

lwd=3

bA = replicate(nBoot,boot(dat,"A"))
bB = replicate(nBoot,boot(dat,"B"))

curve(dnorm(x,meanA,sdA),lwd=lwd,xlim=c(40,65),ylim=c(0,0.5),ylab="density")
curve(dnorm(x,meanB,sdB),type="l",col="red",lwd=lwd,add=TRUE)
points(density(as.numeric(bA["mean",])),type="l",col="grey",lwd=lwd)
points(density(as.numeric(bB["mean",])),type="l",col="pink",lwd=lwd)
legend("topright",legend=c("Population distribution","Sampling distribution (mean)"),col=c("black","grey"),lwd=lwd)

# p-value using t.test function
ttest = t.test(dat$obs[dat$lab=="A"],dat$obs[dat$lab=="B"],var.equal=FALSE)
Tstat = ttest$statistic
ttp = ttest$p.value
ttvep = t.test(dat$obs[dat$lab=="A"],dat$obs[dat$lab=="B"],var.equal=TRUE)$p.value

# p-value using handmade t-test 
# https://github.com/wch/r-source/blob/trunk/src/library/stats/R/t.test.R
mean_A = mean(dat$obs[dat$lab=="B"])
mean_B = mean(dat$obs[dat$lab=="A"])
SD_A = sd(dat$obs[dat$lab=="A"])
SD_B = sd(dat$obs[dat$lab=="B"])
N_A = length(dat$obs[dat$lab=="A"])
N_B = length(dat$obs[dat$lab=="B"])

STDERR_A = SD_A/sqrt(N_A)
STDERR_B = SD_B/sqrt(N_B)
STDERR = sqrt(STDERR_A^2 + STDERR_B^2)

D = (mean_B-mean_A)/STDERR
dof = STDERR^4/(STDERR_A^4/(N_A-1) + STDERR_B^4/(N_B-1))
handp = 2*pt(-abs(D),dof)

# Test statistic: H0 is both datasets transformed to have identical means

# This bootstrap gives us uncertainty about summary (e.g. uncertainty about mean)
boot = function(dat,lab){
  alldat = dat$obs[dat$lab==lab]
  samp = sample(alldat,replace=TRUE)
  return(list(mean=mean(samp),sd=sd(samp),n=length(samp)))
}

Astar = dat$obs[dat$lab=="A"]-mean(dat$obs[dat$lab=="A"])
Bstar = dat$obs[dat$lab=="B"]-mean(dat$obs[dat$lab=="B"])
datstar = dat
datstar$obs = c(Astar,Bstar)
bAstar = replicate(nBoot,boot(datstar,"A"))
bBstar = replicate(nBoot,boot(datstar,"B"))
#Dstar = as.numeric(bBstar["mean",])-as.numeric(bAstar["mean",])/(sqrt((as.numeric(bAstar["sd",])^2)/as.numeric(bAstar["n",])+(as.numeric(bBstar["sd",])^2)/as.numeric(bBstar["n",])))
Dstar = as.numeric(bBstar["mean",])-as.numeric(bAstar["mean",])
D = mean_B - mean_A
bootp = 2*(1+sum(abs(Dstar)>=abs(D)))/(length(Dstar)+1)

plot(density(abs(Dstar)),lwd=lwd)
abline(v=abs(D),col="red",lty=2,lwd=lwd)

# Test statistic: H0 is when labels are shuffled

permute = function(dat){
  datstar = dat
  datstar$lab = sample(dat$lab,length(dat$lab),replace=FALSE)
  D = mean(datstar$obs[datstar$lab=="B"])-mean(datstar$obs[datstar$lab=="A"]) 
  return(D)
}

Dstar = replicate(nBoot,permute(dat))
Dstar = Dstar[!is.na(Dstar)]
permutep = 2*(1+sum(abs(Dstar)>=abs(D),na.rm=TRUE))/(length(Dstar)+1)

plot(density(Dstar),lwd=lwd)
abline(v=D,col="red",lty=2,lwd=lwd)

print("t.test (two-tailed, different var)")
print(ttp)
print("t.test (two-tailed, equal var)")
print(ttvep)
print("hand-coded t-test")
print(handp)
print("bootstrap t-test, data scaled for H0")
print(bootp)
print("bootstrap t-test, data shuffled for H0")
print(permutep)  # Why is this so much bigger?

permutep/bootp
permutep/ttp