#######################################################################
# scanone_cph.R
# For "A guide to QTL mapping with R/qtl", Karl W Broman and Saunak Sen
#
# scanone via Cox proportional hazards model, using H-K regression
#######################################################################
scanone.cph <-
function(cross, pheno.col=1)
{
  require(survival)
  
  pheno <- pull.pheno(cross, pheno.col)
  cross <- subset(cross, ind = !is.na(pheno))
  pheno <- pheno[!is.na(pheno)]

  if(class(pheno) != "Surv")
    stop("Need the phenotype to be of class \"Surv\".")
  
  chrtype <- sapply(cross$geno, class)
  if(any(chrtype=="X")) {
    warning("Dropping X chromosome.")
    cross <- subset(cross, chr=(chrtype != "X"))
  }
  chr <- names(cross$geno)

  result <- NULL
  for(i in 1:nchr(cross)) {
    if(!("prob" %in% names(cross$geno[[i]]))) {
      warning("First running calc.genoprob.")
      cross <- calc.genoprob(cross)
    }
    p <- cross$geno[[i]]$prob

    # pull out map; drop last column of probabilities
    map <- attr(p, "map")
    p <- p[,,-dim(p)[3],drop=FALSE]
      
    lod <- apply(p, 2, function(a,b)
                 diff(coxph(b ~ a)$loglik)/log(10), pheno)

    z <- data.frame(chr=chr[i], pos=map, lod=lod)

    # special names for rows
    w <- names(map)
    o <- grep("^loc-*[0-9]+", w)
    if(length(o) > 0) # inter-marker locations cited as "c*.loc*"
      w[o] <- paste("c",names(cross$geno)[i],".",w[o],sep="")
    rownames(z) <- w

    result <- rbind(result, z)
  }

  class(result) <- c("scanone", "data.frame")
  result
}

# end of scanone_cph.R