## this test is a nonparametric test which extends friedman.test() to cases that allow replicates in blocks ## exact p value is calculated TwoSampleBlockRank.test=function(y,groups, blocks, alternative="two.sided") { ## calculate statsitics stat.rank=0 sizes=NULL ## remove blocks with missing cells index=NULL B=unique(blocks) G=unique(groups) for (i in B){ if(length(unique(groups[blocks==i]))==1) index=c(index, i) } if(!is.null(index)) { for (i in index) {blocks=ifelse(blocks==i, NA, blocks)} temp=!is.na(blocks) blocks=blocks[temp] groups=groups[temp] y=y[temp] } B=unique(blocks) G=unique(groups) lb=length(B) for (i in 1:lb){ stat.rank= sum(rank(y[blocks==B[i]]) [groups[blocks==B[i]]==G[1]])+stat.rank sizes=c(sizes, choose(sum(blocks==B[i]), sum(groups[blocks==B[i]]==G[1]))) } ## NULL Distribution ranksum.matrix=matrix(NA, ncol=max(sizes), nrow=lb) for (i in 1:length(B)) { temp=rowSums(choose.enumerate(rank(y[blocks==B[i]]), n=sum(groups[blocks==B[i]]==G[1]))) ranksum.matrix[i,1:length(temp) ]= temp } b.length=NULL for(i in 1:dim(ranksum.matrix)[[1]]) { b.length=c(b.length, sum(!is.na(ranksum.matrix[i, ]))) } temp=sort(b.length, index.return=T) min.length=min(b.length) b.length=b.length[temp$ix] ranksum.matrix=matrix(ranksum.matrix[temp$ix,], nrow=lb) ranksum=NULL continue1=T index=rep(1, length(B)) while(continue1){ ranksum=c(ranksum, colSums(matrix(ranksum.matrix[, 1:min.length], ncol=min.length,nrow=lb))) if(lb==1) continue1=F else if(index[lb]==b.length[lb]) { continue2=T temp=lb-1 while(continue2) { if(temp==1) {continue2=F; continue1=F} else if(index[temp]==b.length[temp]) {continue2=T; temp=temp-1} else{ continue2=F; index[temp]=index[temp]+1; a=ranksum.matrix[temp,1]; ranksum.matrix[temp, 1:(b.length[temp]-1)]=ranksum.matrix[temp, 2:(b.length[temp])] ranksum.matrix[temp, b.length[temp]]=a; index[(temp+1):length(index)]=1 } } } else{ index[lb]=index[lb]+1 a=ranksum.matrix[lb,1]; ranksum.matrix[lb, 1:(b.length[lb]-1)]=ranksum.matrix[lb, 2:(b.length[lb])] ranksum.matrix[lb, b.length[lb]]=a } } if(alternative=="two.sided") { p.value=2*mean(ranksum>=stat.rank) if(p.value>1) p.value=2*mean(ranksum<=stat.rank) if(p.value>1) p.value=1 } if(alternative=="greater") p.value=mean(ranksum>=stat.rank) if(alternative=="less") p.value=mean(stat.rank>=ranksum) output=list(p.value=p.value, stat=stat.rank) } choose.enumerate=function(x, n=NA){ if(is.na(n)) stop ("The number to choose is not specified!") lx=length(x) nrow=choose(lx, n) output=matrix(0, ncol=n, nrow=nrow) continue1= T count1=1 index=1:n while(continue1) { output[count1, ]=x[index] count1=count1+1; if(count1>nrow) continue1=F else { if(index[n]==length(x)) { continue2=T temp=n-1 max.limit=lx-1 while(continue2){ if(index[temp]==max.limit) {continue2=T; temp=temp-1; max.limit=max.limit-1} else{continue2=F;index[temp]=index[temp]+1; index[temp:n]=seq(index[temp],index[temp]+(n-temp) ) } } } else {index[n]=index[n]+1} } } output }