library(msage) # figure 6 s1a<-runScoreTestComparison(tags=1000,n=2,phi=1,samples=100,b.mult=8,lib.size=c(10000,100000,10000,100000)) v<-calcRankDiff(s1a) png("exact_diffranks.png",width=800,height=500) boxplot(v~col(v),names=c("ML","CML","PL","QL","CR","qCML"),ylab="mean rank difference (qCML-others)") abline(h=0) dev.off() # figure 5 auca<-calcAUCs(s1a,NULL) d<-(1:6-1)*4+1 ind<-c(d,d+1,d+2,d+3) png("aucs.png",width=900,height=400) boxplot(auca[,ind]~col(auca),ylab="AUC",names=rep(c("ML","CML","PL","QL","CR","qCML"),4),cex.axis=.7) mtext("Wald",1,adj=.14,padj=3.5) mtext("Score",1,adj=.38,padj=3.5) mtext("LR",1,adj=.62,padj=3.5) mtext("Exact",1,adj=.865,padj=3.5) abline(v=c(6.5,12.5,18.5)) dev.off() # figure 3 s1000.1.1.30<-runScoreTestComparison(tags=1000,n=2,phi=1,samples=100,b.mult=1,lib.size.range=c(30000,90000)) s1000.1.1.30.5<-runScoreTestComparison(tags=1000,n=5,phi=1,samples=100,b.mult=1,lib.size.range=c(30000,90000)) s2<-s1000.1.1.30 s2$est[1:100,1:6]<-.5 s2$est<-matrix(s2$est[,1],ncol=1) n2<-calcNominal(s2) s5<-s1000.1.1.30.5 s5$est[1:100,1:6]<-.5 s5$est<-matrix(s5$est[,1],ncol=1) n5<-calcNominal(s5) png("tests_truephi.png",height=500,width=900) par(mfrow=c(1,2)) boxplot(n2[,1:4]~col(n2[,1:4]),names=c("Wald","Score","LR","Exact"),main="n1=n2=2",ylab="FP rate"); abline(h=.05); boxplot(n5[,1:4]~col(n5[,1:4]),names=c("Wald","Score","LR","Exact"),main="n1=n2=5",ylab="FP rate"); abline(h=.05); dev.off() # figure 4 s1<-runScoreTestComparison(tags=1000,n=2,phi=1,samples=30,b.mult=1,lib.size.range=c(30000,90000)) aa<-fisherExactScore(s1) res<-fisherNominal(aa) png("fisher_diffests.png",height=600,width=800) boxplot(res~col(res),names=c("ML", "CML", "PL", "QL", "CR", "qCML")) abline(h=.05) dev.off() # figure 1, 2 a.0001.1<-runPhiEstimateComparison1o(samples=1000,n=3,lib.p=.0001,phi=1) a.0005.1<-runPhiEstimateComparison1o(samples=1000,n=3,lib.p=.0005,phi=1) a.0001.25<-runPhiEstimateComparison1o(samples=1000,n=3,lib.p=.0001,phi=0.25) a.0005.25<-runPhiEstimateComparison1o(samples=1000,n=3,lib.p=.0005,phi=0.25) c.0001.1<-runPhiEstimateComparison(samples=1000,tags=100,n=3,lib.p=.0001,phi=1) c.0005.1<-runPhiEstimateComparison(samples=1000,tags=100,n=3,lib.p=.0005,phi=1) c.0001.25<-runPhiEstimateComparison(samples=1000,tags=100,n=3,lib.p=.0001,phi=0.25) c.0005.25<-runPhiEstimateComparison(samples=1000,tags=100,n=3,lib.p=.0005,phi=0.25) png("1sample_w_offsets.png",750,600) par(mfrow=c(2,2)) plotPhiEstimateComparison(a.0001.1) plotPhiEstimateComparison(a.0005.1) plotPhiEstimateComparison(a.0001.25) plotPhiEstimateComparison(a.0005.25) dev.off() png("100tags_w_offsets.png",750,600) par(mfrow=c(2,2)) plotPhiEstimateComparison(c.0001.1) plotPhiEstimateComparison(c.0005.1) plotPhiEstimateComparison(c.0001.25) plotPhiEstimateComparison(c.0005.25) dev.off()