R/RandTree.R

In TBEST: Tree Branches Evaluated Statistically for Tightness

RandTree<-function(mydata, myinput, mystat, mymethod, mymetric, rand.fun = c("shuffle.column",
#function(mydata, mystat, mymethod, mymetric, rand.fun = c("shuffle.column", 
    "shuffle.block", "define.function"), by.block = NA, metric.args = list(), 
    rand.args = list())
{
    #myinput <- mydata$myinput
    ntest <- mydata$nperm
    indextable <- TreeStat(myinput, mystat = mystat, method = mymethod, 
        metric = mymetric, metric.args = metric.args)
    statnames <- mystat
    #nullstat <-vector("list",length(statnames))
    #names(nullstat) <- statnames
    if(any(statnames=="slb")){slbcount<-rep(0,ntest)}
    allcounts <- matrix(ncol = length(statnames), nrow = nrow(indextable), data = 0)
    colnames(allcounts) <- statnames
    for (i in 1:ntest) {
        if (rand.fun == "shuffle.column"){ 
            myrdata <- apply(myinput, 2, sample)
        }else if (rand.fun == "shuffle.block") {
	    if(is.na(by.block[1]))stop("by.block needs to be specified")
            myrdata <- t(myinput)
            myrlist <- by(myrdata, by.block, FUN = byfactor)
            for (j in 1:length(myrlist)) {
                if (j == 1) {
                  myrdata <- myrlist[[j]]
                }
                else {
                  myrdata <- rbind(myrdata, myrlist[[j]])
                }
            }
            myrdata <- t(myrdata)
        }else if (rand.fun == "define.function") {
            define.function <- match.fun(define.function)
            myrand.args <- vector("list", length(rand.args) + 
                1)
            myrand.args[[1]] <- myinput
            if (length(myrand.args) > 1) {
                myrand.args[2:length(myrand.args)] <- rand.args
            }
            myrdata <- do.call(define.function, myrand.args)
        }
        rindextable <- TreeStat(myrdata, mystat = mystat, method = mymethod, 
            metric = mymetric, metric.args = metric.args)
	if(any(statnames!="slb")){
        size <-rindextable[,"clustersize"]
        for (statname in statnames) {
                rstat <- rindextable[,statname]
                statmax <- max(rstat)
                randomX <- sort(size + 0.5 * rstat/statmax)
                rmatch <- bestmatch(rsize=sort(size),size = indextable[,
                        "clustersize"])
                rmatchl <- bestmatchl(rsize=sort(size),size = indextable[,
                        "clustersize"])
                data <- 2 * statmax * (randomX - sort(size))[rmatch]
                datal <- 2 * statmax * (randomX - sort(size))[rmatchl]
                mydata <- pmax(data, datal)
		allcounts[, statname] <- allcounts[, statname] + (mydata >=
            		indextable[, statname])
                #if(i==1){nullstat[[statname]]<-mydata}
                #else{nullstat[[statname]]<-c(nullstat[[statname]],mydata)}
        }
	}else if(any(statnames=="slb")){
		slbcount[i]<-sum(rindextable[nrow(rindextable),"slb"]>=indextable[nrow(rindextable),"slb"])
	}
    }
    if(any(statnames=="slb")){allcounts<-sum(slbcount)}
    return(list(allcounts, ntest))
}