#functions to generate bootstrap null distribution
#theta0 is the value of the test statistics under the complete null hypthesis
#tau0 is the scaling parameter (upper bound on variance of test statistics)
boot.null <- function(X,label,stat.closure,W=NULL,B=1000,test,nulldist,theta0=0,tau0=1,marg.null=NULL,marg.par=NULL,ncp=0,perm.mat,alternative="two.sided",seed=NULL,cluster=1,dispatch=0.05,keep.nulldist,keep.rawdist){
cat("running bootstrap...\n")
X<-as.matrix(X)
n<-ncol(X)
p<-nrow(X)
if(!(is.vector(W) | is.matrix(W) | is.null(W))) stop("W must be a vector or a matrix")
if(is.null(W))W<-matrix(1,nrow=p,ncol=n)
if(is.vector(W)){
if(length(W)==n) W<-matrix(W,nrow=p,ncol=n,byrow=TRUE)
if(length(W)==p) W<-matrix(W,nrow=p,ncol=n)
if(length(W)!=n & length(W)!=p) stop("Length of W does not match dim(X)")
}
if(is.matrix(W) & (dim(W)[1]!=p | dim(W)[2]!=n)) stop("W and X must have same dimension")
# Dispatch to cluster
if (is.numeric(cluster)) {
if(!is.null(seed)) set.seed(seed)
muboot <- boot.resample(X,label,p,n,stat.closure,W,B,test)
}
else {
autoload("clusterApply","snow")
autoload("clusterApplyLB","snow")
if(!is.null(seed)) clusterApply(cluster, seed, set.seed)
else clusterApply(cluster, runif(length(cluster), max=10000000), set.seed)
# Create vector of jobs to dispatch
if ((dispatch > 0) & (dispatch < 1)){
BtoNodes <- rep(B*dispatch, 1/dispatch)
} else {
BtoNodes <- rep(dispatch, B/dispatch)
}
FromCluster <- clusterApplyLB(cluster, BtoNodes, boot.resample,X=X,label=label,p=p,n=n,stat.closure=stat.closure,W=W, test=test)
muboot <- matrix(unlist(FromCluster), nrow=nrow(X))
}
Xnames<-dimnames(X)[[1]]
dimnames(muboot)<-list(Xnames,paste(1:B))
#fill in any nas by resampling some more
nas<-(is.na(muboot)|muboot=="Inf"|muboot=="-Inf")
count<-0
while(sum(nas)){
count<-count+1
if(count>1000) stop("Bootstrap null distribution computation terminating. Cannot obtain distribution without missing values after 1000 attempts. This problem may be resolved if you try again with a different seed.")
nascols<-unique(col(muboot)[nas])
for(b in nascols){
samp<-sample(n,n,replace=TRUE)
Xb<-X[,samp]
Wb<-W[,samp]
if(p==1){
Xb<-t(as.matrix(Xb))
Wb<-t(as.matrix(Wb))
}
Tb<-get.Tn(Xb,stat.closure,Wb)
muboot[,b]<-Tb[3,]*Tb[1,]/Tb[2,]
}
nas<-is.na(muboot)
}
rawboot <- matrix(nrow=0,ncol=0)
if(keep.rawdist) rawboot <- muboot
if(nulldist=="boot") muboot <- center.scale(muboot, theta0, tau0, alternative)
if(nulldist=="boot.cs") muboot <- center.scale(muboot, theta0, tau0, alternative)
if(nulldist=="boot.ctr") muboot <- center.only(muboot, theta0, alternative)
if(nulldist=="boot.qt") muboot <- quant.trans(muboot, marg.null, marg.par, ncp, alternative, perm.mat)
out <- list(muboot=muboot, rawboot=rawboot)
out
}
center.only <- function(muboot,theta0,alternative){
muboot<-(muboot-apply(muboot,1,mean))+theta0
if(alternative=="greater") muboot <- muboot
else if(alternative=="less") muboot <- -muboot
else muboot <- abs(muboot)
}
center.scale <- function(muboot, theta0, tau0, alternative){
muboot<-(muboot-apply(muboot,1,mean))*sqrt(pmin(1,tau0/apply(muboot,1,var)))+theta0
if(alternative=="greater") muboot <- muboot
else if(alternative=="less") muboot <- -muboot
else muboot <- abs(muboot)
}
quant.trans <- function(muboot, marg.null, marg.par, ncp, alternative, perm.mat){
### NB: Sanity checks occur outside this function at the beginning of MTP.
m <- dim(muboot)[1]
B <- dim(muboot)[2]
ranks <- t(apply(muboot,1,rank,ties.method="random"))
Z.quant <- switch(marg.null,
normal = marg.samp(marg.null="normal",marg.par,m,B,ncp),
t = marg.samp(marg.null="t",marg.par,m,B,ncp),
f = marg.samp(marg.null="f",marg.par,m,B,ncp),
perm = perm.mat)
Z.quant <- t(apply(Z.quant,1,sort))
### Left code like this for transparency. Could just as easily use quantile()
### for this first part, although it would be redundant.
if(marg.null!="perm"){
for(i in 1:m){
Z.quant[i,] <- Z.quant[i,][ranks[i,]]
}
}
else{
Z.quant <- t(apply(Z.quant,1,quantile,probs=seq(0,1,length.out=B),na.rm=TRUE))
for(i in 1:m){
Z.quant[i,] <- Z.quant[i,][ranks[i,]]
}
}
if(alternative=="greater") Z.quant <- Z.quant
else if(alternative=="less") Z.quant <- -Z.quant
else Z.quant <- abs(Z.quant)
Z.quant
}
boot.resample <- function (X, label, p, n, stat.closure, W, B, test){
muboot <- matrix(0, nrow = p, ncol = B)
samp <- sample(n, n * B, replace = TRUE)
if (any(test == c("t.twosamp.equalvar", "t.twosamp.unequalvar",
"f"))) {
label <- as.vector(label)
uniqlabs <- unique(label)
num.group <- length(uniqlabs)
groupIndex <- lapply(1:num.group, function(k) which(label ==
uniqlabs[k]))
if(sum(is.na(label))){
naindex<-c(1:num.group)[is.na(uniqlabs)]
groupIndex[[naindex]]<-which(is.na(label))
}
obs <- sapply(1:num.group, function(x) length(groupIndex[[x]]))
samp <- lapply(1:num.group, function(k) matrix(NA, nrow = B,
ncol = obs[k]))
for (j in 1:B) {
for (i in 1:num.group) {
uniq.obs <- 1
count <- 0
while (uniq.obs == 1) {
count <- count + 1
samp[[i]][j, ] <- sample(groupIndex[[i]], obs[i],
replace = TRUE)
uniq.obs <- length(unique(samp[[i]][j, ]))
if (count > 1000)
stop("Bootstrap null distribution computation terminating. Cannot obtain bootstrap sample with at least 2 unique observations after 1000 attempts. Sample size may be too small for bootstrap procedure but this problem may be resolved if you try again with a different seed.")
}
}
}
samp <- as.vector(t(matrix(unlist(samp), nrow = B, ncol = sum(obs))))
}
else if (test == c("f.twoway")) {
label <- as.vector(label)
utreat <- unique(label)
num.treat <- length(utreat)
num.block <- length(gregexpr("12", paste(label, collapse = ""))[[1]])
ublock <- 1:num.block
Breaks <- c(0, gregexpr(paste(c(num.treat, 1), collapse = ""),
paste(label, collapse = ""))[[1]], n)
BlockNum <- sapply(1:num.block, function(x) Breaks[x +
1] - Breaks[x])
block <- unlist(lapply(1:num.block, function(x) rep(x,
BlockNum[x])))
groupIndex <- lapply(1:num.block, function(j) sapply(1:num.treat,
function(i) which(label == utreat[i] & block == ublock[j])))
obs <- sapply(1:num.block, function(x) sapply(1:num.treat,
function(y) length(groupIndex[[x]][,y])))
samp <- lapply(1:(num.treat * num.block), function(k) matrix(NA,
nrow = B, ncol = obs[k]))
for (k in 1:B) {
for (i in 1:num.block) {
for (j in 1:num.treat) {
uniq.obs <- 1
count <- 0
while (uniq.obs == 1) {
count <- count + 1
samp[[(i - 1) * num.treat + j]][k, ] <- sample(groupIndex[[i]][,j],
obs[j, i], replace = TRUE)
uniq.obs <- length(unique(samp[[(i - 1) *
num.treat + j]][k, ]))
if (count > 1000)
stop("Bootstrap null distribution computation terminating. Cannot obtain bootstrap sample with at least 2 unique observations after 1000 attempts. Sample size may be too small for bootstrap procedure but this problem may be resolved if you try again with a different seed.")
}
}
}
}
samp <- as.vector(t(matrix(unlist(samp), nrow = B, ncol = sum(obs))))
}
cat("iteration = ")
muboot <- .Call("bootloop", stat.closure, as.numeric(X),
as.numeric(W), as.integer(p), as.integer(n), as.integer(B),
as.integer(samp), NAOK = TRUE)
cat("\n")
muboot <- matrix(muboot, nrow = p, ncol = B)
}
