# Misor - calculates odds ratios for phase-known haplotypes measured with error
# See Govindarajulu et al. (submitted) for details
# Peter Kraft (Fall 2005)

misor <- function(haps,q,e=0,mode="DOM",a=log(.001/.999),b=log(2),cut=.001)
{

  f1 <- exp(a+b)/(1+exp(a+b))
  f0 <- exp(a)/(1+exp(a))

  miss <- as.matrix(dist(haps,method="manhattan"))
  miss <- e^miss * (1-e)^(NCOL(haps)-miss)

  qstar <- as.vector(q %*% miss)
  J <- sum(q>cut)
  output <- matrix(0,J,20)

  for (g in (1:NROW(haps))[q>cut])
  {
    qr <- q[g]                                 # Define risk haplotype frequency
    if (mode=="DOM")
    {
      k <- f1*(qr^2+2*qr*(1-qr)) + f0*(1-qr)^2   # Pr(case) = population prevalence -- DOMINANT
      PrDZ <- f1*(qr^2+2*qr*(1-qr))
      PrdZ <- (1-f1)*(qr^2+2*qr*(1-qr))
      p11 <- PrDZ/k                              # True CaCo cell counts
      p10 <- 1-p11
      p01 <- PrdZ/(1-k)
      p00 <- 1-p01
      or  <- p11*p00/(p10*p01)                   # True OR
    }
    if (mode=="REC")
    {
      k <- f1*qr^2 + f0*(1-qr^2)                 # Pr(case) = population prevalence -- RECESSIVE
      PrDZ <- f1*qr^2
      PrdZ <- (1-f1)*qr^2
      p11 <- PrDZ/k                              # True CaCo cell counts
      p10 <- 1-p11
      p01 <- PrdZ/(1-k)
      p00 <- 1-p01
      or  <- p11*p00/(p10*p01)                   # True OR
    }
    qrs <- qstar[g]                              # This is observed frequency of risk haplotype
    PrDZs <- 0
    PrdZs <- 0
    PrZZ <- 0
    Przz <- 0
    PrZ <- 0
    PrZs <- 0
    Przz <- 0
    mm1 <- miss[g,g]
    mm0 <- sum((1-miss[-g,g]) * q[-g]) / sum(q[-g])
    m <- c(mm0,mm1)
    for (i1 in 0:1)
    for (i2 in 0:1)
    for (m1 in 0:1)
    for (m2 in 0:1)
    {
      PrHM <- qr^(i1+i2) * (1-qr)^(2-i1-i2) * m[i1+1]^m1 * (1-m[i1+1])^(1-m1) * m[i2+1]^m2 * (1-m[i2+1])^(1-m2)
      if (mode == "DOM") {isZs <- (i1*m1 + i2*m2 + (1-i1)*(1-m1) + (1-i2)*(1-m2)) > 0}
      if (mode == "DOM") {isZ  <- (i1 + i2) > 0}
      if (mode == "REC") {isZs <- (i1*m1 + i2*m2 + (1-i1)*(1-m1) + (1-i2)*(1-m2)) == 2}
      if (mode == "REC") {isZ  <- (i1 + i2) == 2}
      PrDZs <- PrDZs + f1*PrHM*isZs*isZ + f0*PrHM*isZs*(!isZ)
      PrdZs <- PrdZs + (1-f1)*PrHM*isZs*isZ + (1-f0)*PrHM*isZs*(!isZ)
      PrZZ <- PrZZ + PrHM*isZs*isZ
      Przz <- Przz + PrHM*(!isZs)*(!isZ)
      PrZ <- PrZ +  PrHM*isZ
      PrZs <- PrZs +  PrHM*isZs
    }
    p11s <- PrDZs/k
    p10s <- 1-p11s
    p01s <- PrdZs/(1-k)
    p00s <- 1-p01s
    ors <- p11s*p00s / (p10s*p01s)
    sensZ <- PrZZ/PrZ
    specZ <- Przz/(1-PrZ)
    pvpZ <- PrZZ/PrZs
    pvnZ <- Przz/(1-PrZs)
    output[g,] <-c(qr,PrDZ,PrdZ,p11,p10,p01,p00,or,qrs,PrDZs,PrdZs,p11s,p10s,p01s,p00s,ors,sensZ,specZ,pvpZ,pvnZ)
    colnames(output) <- c("q","PrDZ","PrdZ","p11","p10","p01","p00","or","qs","PrDZs","PrdZs",
                          "p11s","p10s","p01s","p00s","ors","sensZ","specZ","pvpZ","pvnZ")
  }
  output
}