# 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))