# misor examples -- c.f. Govindarajulu et al. (submitted) # Peter Kraft # 13 March 2006 source("misor.txt") ######################################################### # Plots comparing twolocus bias to three locus bias when third locus is redundant m <- matrix(c(0,0,0,0,1,1,1,0,0,1,1,1),4,3,byrow=T) outtab <- array(0,c(99,4,6)) for (i1 in 1:99) { q1 <- i1/100 q <- c(q1,(1-q1)/3,(1-q1)/3,(1-q1)/3) outtab[i1,,] <- misor(m[,1:2],q,e=.01)[,c(9,16:20)] } plot((1:99)/100,outtab[,1,2],ylim=c(1,2),type="n",ylab="OR estimate",xlab="haplotype A-B frequency",main="Two-locus Haplotype - OR estimates") lines((1:99)/100,outtab[,1,2]) outtab <- array(0,c(99,4,6)) for (i1 in 1:99) { q1 <- i1/100 q <- c(q1,(1-q1)/3,(1-q1)/3,(1-q1)/3) outtab[i1,,] <- misor(m,q,e=.01)[,c(9,16:20)] } lines((1:99)/100,outtab[,1,2],lty=2) ######################################################### # Plots for the paper m1 <- matrix(c(0,0,0,0,1,1,1,0,0,1,1,1),4,3,byrow=T) outtab1 <- array(0,c(99,4,6)) for (i1 in 1:99) { q1 <- i1/100 q <- c(q1,(1-q1)/3,(1-q1)/3,(1-q1)/3) outtab1[i1,,] <- misor(m1[,1:2],q,e=.01)[,c(9,16:20)] } m2 <- matrix(c(0,0,0,0,1,0,1,0,0,1,1,0, 0,0,1,0,1,1,1,0,1,1,1,1),8,3,byrow=T) outtab2 <- array(0,c(99,8,6)) for (i1 in 1:99) { q1 <- i1/100 q <- c(q1,rep((1-q1)/7,7)) outtab2[i1,,] <- misor(m2[,1:2],q,e=.01)[,c(9,16:20)] } par(mfrow=c(3,1)) plot((1:99)/100,outtab1[,1,2],ylim=c(1,2.1),type="n",ylab="OR estimate",xlab="haplotype frequency",main="OR estimates") lines((1:99)/100,outtab1[,1,2]) lines((1:99)/100,outtab2[,1,2],lty=2) abline(h=2,lty=3) plot((1:99)/100,outtab1[,1,5],ylim=c(0,1),type="n",ylab="PVP",xlab="haplotype frequency",main="Predictive value positive") lines((1:99)/100,outtab1[,1,5]) lines((1:99)/100,outtab2[,1,5],lty=2) plot((1:99)/100,outtab1[,1,6],ylim=c(0,1),type="n",ylab="PVN",xlab="haplotype frequency",main="Predictive value negative") lines((1:99)/100,outtab1[,1,6]) lines((1:99)/100,outtab2[,1,6],lty=2) legend(locator(1),lty=c(1,2),legend=c("Two-locus haplotypes","Three-locus haplotypes")) ############################################################### # Sensitivity specificity et cetera par(mfrow=c(1,1)) plot((1:99)/100,outtab[,1,3],type="n",ylim=c(min(outtab[,,3:4]),max(outtab[,,3:4])),ylab=" ",xlab="haplotype A-B frequency",main="Sensitivity and specificity") lines((1:99)/100,outtab[,1,3]) lines((1:99)/100,outtab[,2,3],col="blue") lines((1:99)/100,outtab[,4,3],col="red") lines((1:99)/100,outtab[,1,4],lty=2) lines((1:99)/100,outtab[,2,4],col="blue",lty=2) lines((1:99)/100,outtab[,4,4],col="red",lty=2) legend(locator(1),legend=c("AB sensitivity","Ab and aB sensitivity","ab sensitivity","AB specificity","Ab and aB specificity","ab specificity"),col=c("black","blue","red"),lty=c(1,1,1,2,2,2)) plot((1:99)/100,outtab[,1,3],type="n",ylim=c(min(outtab[,,5:6]),max(outtab[,,5:6])),ylab=" ",xlab="haplotype A-B frequency",main="Predictive value positive and negative") lines((1:99)/100,outtab[,1,5]) lines((1:99)/100,outtab[,2,5],col="blue") lines((1:99)/100,outtab[,4,5],col="red") lines((1:99)/100,outtab[,1,6],lty=2) lines((1:99)/100,outtab[,2,6],col="blue",lty=2) lines((1:99)/100,outtab[,4,6],col="red",lty=2) legend(locator(1),legend=c("AB PVP","Ab and aB PVP","ab PVP","AB PVN","Ab and aB PVN","ab PVN"),col=c("black","blue","red"),lty=c(1,1,1,2,2,2))