#setwd("C:/Users/kejima/Dropbox/Andrew_murine/JAX/type-I-error") alpha_level = 0.05 # significance level tiff("Figure3.tif", res = 300, width = 8, height = 4, units = 'in') layout(matrix(c(1,2), 2, 1, byrow = T)) for(case in c(1,2)){ # case 1 or 2 if(alpha_level==0.05&case==1){Y=c(0.0,0.20);Y1=c(0.0,0.10)}; if(alpha_level==0.05&case==2){Y=c(0.0,0.20);Y1=c(0.0,0.10)} if(alpha_level==0.01&case==1){Y=c(0.0,0.20);Y1=c(0.0,0.04)}; if(alpha_level==0.01&case==2){Y=c(0.0,0.20);Y1=c(0.0,0.04)} if(alpha_level==0.005&case==1){Y=c(0.0,0.08);Y1=c(0.0,0.02)}; if(alpha_level==0.005&case==2){Y=c(0.0,0.08);Y1=c(0.0,0.02)} if(alpha_level==0.001&case==1){Y=c(0.0,0.035);Y1=c(0.0,0.01)}; if(alpha_level==0.001&case==2){Y=c(0.0,0.035);Y1=c(0.0,0.01)} data = NULL for(n in c(3,4,5,8,12,16,20)){ # sample size for(s in c('f','m')){ # sex for(stid in c(140,142,141,144,143,162)){ # strain if(s=='f'&stid==144) next data = rbind(data, read.csv(paste(stid,"-",n,"-",s,"-case-",case,"-type-I-error-",alpha_level,".csv",sep=""),header=F)) } } } mu <- data[,seq(1,13,3)]; lc <- data[,seq(2,14,3)]; uc <- data[,seq(3,15,3)] par(mgp=c(1.5,1.0,0), mai=c(0.5,0.5,0.2,0.3), ps=9) for(i in 1:4){ # for each test plot(c(1:77)+(i-1)*77,mu[,i],xlim=c(0,77*5),pch=c(rep(1,5),rep(16,6)),col=c(rep(2,5),rep(1,6)),axes=F, ylim=Y, xlab="", ylab="",cex=1.0) arrows(c(1:77)+(i-1)*77, lc[,i], c(1:77)+(i-1)*77,col=c(rep(2,5),rep(1,6)), uc[,i], angle = 90, length = 0.02) arrows(c(1:77)+(i-1)*77, uc[,i], c(1:77)+(i-1)*77,col=c(rep(2,5),rep(1,6)), lc[,i], angle = 90, length = 0.02) par(new=T) } axis(2, pos = -1, at = seq(min(Y),max(Y),max(Y)/5), adj = 0, las = 2, cex.axis=1.0) # y axis segments(-1,alpha_level,77*4,alpha_level, lty=2); i=5 # bootstrap par(new=T) plot(c(1:77)+(i-1)*77,mu[,i],xlim=c(0,77*5),pch=c(rep(1,5),rep(16,6)),col=c(rep(2,5),rep(1,6)),axes=F, ylim=Y1, xlab="", ylab="",cex=1.0) arrows(c(1:77)+(i-1)*77, lc[,i], c(1:77)+(i-1)*77,col=c(rep(2,5),rep(1,6)), uc[,i], angle = 90, length = 0.02) arrows(c(1:77)+(i-1)*77, uc[,i], c(1:77)+(i-1)*77,col=c(rep(2,5),rep(1,6)), lc[,i], angle = 90, length = 0.02) abline(v=c(1:3)*77+0.5, lty=1, lwd=1);abline(v=4*77+0.5, lty=1, lwd=2); segments(77*4,alpha_level,77*5,alpha_level, lty=2); mtext(paste("Significance level = ",alpha_level,", Case ",case,sep="")) mtext(side=1, c("Student's t test","Welch's t test","Wilcoxon test","Permutation test","Bootstrap test"), line=1, at = c(0:4)*77+40) mtext(side=1, rep(c(3,4,5,8,12,16,20),5), line = 0.0, at = rep(0:4,each=7)*77+seq(5.5,77,11)) # x axis mtext(side=2, "Type I error rate", line=1.0) axis(4, at = seq(min(Y1),max(Y1),max(Y1)/5), adj = 0, las = 2, pos=5*77+2, cex.axis=1.0) # y axis } dev.off() tiff("Figure4.tif", res = 300, width = 8, height = 4, units = 'in') layout(matrix(c(1,2), 2, 1, byrow = T)) for(case in c(1,2)){ # case 1 or 2 if(alpha_level==0.05&case==1){Y=c(0.0,0.20);Y1=c(0.0,0.10)}; if(alpha_level==0.05&case==2){Y=c(0.0,0.20);Y1=c(0.0,0.10)} if(alpha_level==0.01&case==1){Y=c(0.0,0.20);Y1=c(0.0,0.04)}; if(alpha_level==0.01&case==2){Y=c(0.0,0.20);Y1=c(0.0,0.04)} if(alpha_level==0.005&case==1){Y=c(0.0,0.08);Y1=c(0.0,0.02)}; if(alpha_level==0.005&case==2){Y=c(0.0,0.08);Y1=c(0.0,0.02)} if(alpha_level==0.001&case==1){Y=c(0.0,0.035);Y1=c(0.0,0.01)}; if(alpha_level==0.001&case==2){Y=c(0.0,0.035);Y1=c(0.0,0.01)} data = NULL for(n in c(3,4,5,8,12,16,20)){ # sample size for(s in c('f','m')){ # sex for(stid in c(140,142,141,144,143,162)){ # strain if(s=='f'&stid==144) next data = rbind(data, read.csv(paste(stid,"-",n,"-",s,"-case-",case,"-type-I-error-",alpha_level,".csv",sep=""),header=F)) } } } mu <- data[,seq(1,13,3)]; lc <- data[,seq(2,14,3)]; uc <- data[,seq(3,15,3)] temp1 = (lc>alpha_level) # count if the lower bound of 95% CI is above significance level temp2 = matrix(0,7,5) temp3 = (ucalpha_level) # count if the lower bound of 95% CI is above significance level temp2 = matrix(0,7,5) temp3 = (uc