BLANKET=function(ranks1,ranks2,include.a=F){ n=length(ranks1) temp.matrix=outer(ranks1,ranks2,FUN="==") num.list1=outer(ranks1,rep(1,n),FUN="*") num.list2=outer(rep(1,n),ranks1,FUN="*") temp1=apply(temp.matrix,1,cumsum) a=apply(temp1,1,cumsum) b=num.list1-a c=num.list2-a d=n-a-b-c exp.a=(num.list1/n)*num.list2 exp.b=num.list1-exp.a exp.c=num.list2-exp.a exp.d=n-exp.a-exp.b-exp.c z=sqrt((a-exp.a)^2/(exp.a)+(b-exp.b)^2/(exp.b)+ (c-exp.c)^2/(exp.c)+(d-exp.d)^2/(exp.d))*sign(a-exp.a) pstat=1-pnorm(z) pfisher=1-phyper(a-1, num.list2, n-num.list2, num.list1) if (include.a==T){ output=list(pstat,pfisher,a) names(output)=c("pstat","pfisher","a") } else { output=list(pstat,pfisher) names(output)=c("pstat","pfisher") } return(output) } ####### run an example n=100 ranks1=1:n ranks2=ranks1 lowerset=1:(n*.1) ranks2[lowerset]=sample(ranks2[lowerset]) ranks2[-lowerset]=sample(ranks2[-lowerset]) date() result=BLANKET(ranks1,ranks2) date() image(log10(result\$pfisher)) ### get empirical distribution of minima num.perms=1000 num.common=c(400,500) pmin.matrix=matrix(1,num.perms,length(num.common)) pmin20.matrix=matrix(1,num.perms,length(num.common)) pmin.a10.matrix=matrix(1,num.perms,length(num.common)) for (k in (1:length(num.common))){ pmin=rep(0,num.perms) pmin20=rep(0,num.perms) pmin.a10=rep(0,num.perms) for (i in (1:num.perms)){ ranks1=c(1:num.common[k]) ranks2=sample(ranks1) result=BLANKET(ranks1,ranks2,include.a=T) pmin[i]=min(result\$pfisher) pmin20[i]=min(result\$pfisher[1:20,1:20]) pmin.a10[i]=min(result\$pfisher[result\$a<=10]) print(c(i,k)) } pmin.matrix[,k]=pmin pmin20.matrix[,k]=pmin20 pmin.a10.matrix[,k]=pmin.a10 } y1=apply(pmin.matrix,2,quantile,.05) y2=apply(pmin20.matrix,2,quantile,.05) y3=apply(pmin.a10.matrix,2,quantile,.05) plot(num.common,y1,type="b",ylim=c(0,max(y1,y3,y2))) lines(num.common,y2) lines(num.common,y3,lty=3) plot(num.common,apply(pmin20.matrix,2,quantile,.05)) plot(num.common,apply(pmin.a10.matrix,2,quantile,.05)) cbind(num.common,apply(pmin.matrix,2,quantile,.05),apply(pmin.a10.matrix,2,quantile,.05)) temp=read.table("C:\\Documents and Settings\\Administrator\\Desktop\\for_FWright_test_case.txt",sep="\t") p1=temp[,2] p2=temp[,4] ranks1=rank(p1) ranks2=rank(p2) plot(ranks1,ranks2) result=BLANKET(ranks1,ranks2) ########### temp1=read.table("pmin.matrix") temp2=read.table("pmin20.matrix") temp3=read.table("pmin.a10.matrix") temp1=temp1[,1:3] temp2=temp2[,1:3] temp3=temp3[,1:3] threshold.all=c(apply(temp1,2,quantile, .05),apply(pmin.matrix,2,quantile,.05)) threshold.20=c(apply(temp2,2,quantile, .05),apply(pmin20.matrix,2,quantile,.05)) threshold.a10=c(apply(temp3,2,quantile, .05),apply(pmin.a10.matrix,2,quantile,.05)) plot(c(100,200,300,400,500),threshold.all,type="b",ylim=c(0,max(threshold.20))) lines(c(100,200,300,400,500),threshold.20) lines(c(100,200,300,400,500),threshold.a10,lty=2) cbind(c(100,200,300,400,500),threshold.all,threshold.20,threshold.a10) plot(c(100,200,300,400,500),log10(threshold.all),type="b") lines(c(100,200,300,400,500),log10(threshold.20)) lines(c(100,200,300,400,500),log10(threshold.a10),lty=2) **** #examples BB=read.table("C:\\Users\\Kim\\Desktop\\BB_lists.txt",sep="\t") BS=read.table("C:\\Users\\Kim\\Desktop\\BS_lists.txt",sep="\t") Hammond=read.table("C:\\Users\\Kim\\Desktop\\Hammond_lists.txt",sep="\t") par(mfrow=c(2,2)) BB.n=dim(BB)[1] BB.ranks1=c(1:BB.n) BB.ranks2=c(1:BB.n)[order(BB[,2])] BB.ranks2=BB.ranks2[BB[,1]] #cbind(BB.ranks1,BB.ranks2)[1:20,] plot(BB.ranks1,BB.ranks2) BS.n=dim(BS)[1] BS.ranks1=c(1:BS.n) BS.ranks2=c(1:BS.n)[order(BS[,2])] BS.ranks2=BS.ranks2[BS[,1]] plot(BS.ranks1,BS.ranks2) Hammond.n=dim(Hammond)[1] Hammond.ranks1=c(1:Hammond.n) Hammond.ranks2=c(1:Hammond.n)[order(Hammond[,2])] Hammond.ranks2=Hammond.ranks2[Hammond[,1]] plot(Hammond.ranks1,Hammond.ranks2) ###### par(mfrow=c(2,2)) BB.result=BLANKET(BB.ranks1,BB.ranks2,include.a=T) BB.p=c(min(BB.result\$pfisher), min(BB.result\$pfisher[BB.result\$a<=10]), min(BB.result\$pfisher[1:20,1:20])) image(log10(BB.result\$pfisher)) BS.result=BLANKET(BS.ranks1,BS.ranks2,include.a=T) min(BS.result\$pfisher[BS.result\$a<=10]) BS.p=c(min(BS.result\$pfisher), min(BS.result\$pfisher[BS.result\$a<=10]), min(BS.result\$pfisher[1:20,1:20])) image(log10(BS.result\$pfisher)) Hammond.result=BLANKET(Hammond.ranks1,Hammond.ranks2,include.a=T) min(Hammond.result\$pfisher[Hammond.result\$a<=10]) Hammond.p=c(min(Hammond.result\$pfisher), min(Hammond.result\$pfisher[Hammond.result\$a<=10]), min(Hammond.result\$pfisher[1:20,1:20])) image(log10(Hammond.result\$pfisher)) pvalue.matrix=rbind(BB.p,BS.p,Hammond.p) dimnames(pvalue.matrix)[[2]]=c("min.all","min.a<=10","min.20X20")