
plot.reliability.data.diagram <- function(x,sort=TRUE,...) {
  if(is.numeric(x)==TRUE) {
    n <- length(x)
    if(sort==TRUE) {
      x.time <- x[sort.list(x,decreasing=TRUE)]
      x.status <- rep(1,n)
      id <- c(1:n)
    } else {
      x.time <- x
      x.status <- rep(1,n)
      id <- c(1:n)
    }
  } else if (is.data.frame(x)==TRUE) {
    n <- dim(x)[1]
    if(sort==TRUE) {
      x.time <- x[sort.list(x[,1],decreasing=TRUE),1]
      x.status <- x[sort.list(x[,1],decreasing=TRUE),2]
      id <- row.names(x)[sort.list(x[,1],decreasing=TRUE)]
    } else {
      x.time <- x[,1]
      x.status <- x[,2]
      id <- row.names(x)
    }
  } else if (is.matrix(x)==TRUE) {
    n <- dim(x)[1]
    if(sort==TRUE) {
      x.time <- x[sort.list(x[,1],decreasing=TRUE),1]
      x.status <- x[sort.list(x[,1],decreasing=TRUE),2]
      id <- rawnames(x)[sort.list(x[,1],decreasing=TRUE)]
    } else {
      x.time <- x[,1]
      x.status <- x[,2]
      id <- rawnames(x)
    }
  }
  plot(x.time,c(1:n),
       xlim=c(0,max(x.time)),
       type="n",
       ylab="#",
       xlab="Lifetime data",
       yaxt="n")
  axis(2,1:n,labels=c(n:1))
  points(x.time,c(1:n),pch=c(1,4)[x.status+1],...)
  for( i in c(1:n) ) {
    lines(c(0,x.time[i]),c(i,i),...)
  }
}

compare.histgram.and.densities <- function(X,hist=TRUE,plot=TRUE) {
  require(MASS)
  est.weibull <- fitdistr(X$time,densfun="weibull")
  est.lnorm <- fitdistr(X$time,densfun="lognormal")
  est.gamma <- fitdistr(X$time,densfun="gamma")
  if(plot==TRUE) {
    if(hist==TRUE) {
      hist(X$time,freq=FALSE,ylim=c(0,0.3),
           xlab="Lifetime")
      plot(function(x){return(dlnorm(x,meanlog=est.lnorm$estimate[1],
                                     sdlog=est.lnorm$estimate[2]))},
           xlim=c(0,12),
           col="blue",
           add=TRUE)
    } else {
      plot(function(x){return(dlnorm(x,meanlog=est.lnorm$estimate[1],
                                     sdlog=est.lnorm$estimate[2]))},
           xlim=c(0,12),
           col="blue",xlab="Lifetime",ylab="Density")
    }
    plot(function(x){return(dweibull(x,shape=est.weibull$estimate[1],
                                     scale=est.weibull$estimate[2]))},
         xlim=c(0,12),
         col="red",
         add=TRUE)
    plot(function(x){return(dgamma(x,shape=est.gamma$estimate[1],
                                   rate=est.gamma$estimate[2]))},
         xlim=c(0,12),
         col="green",
         add=TRUE)
    legend("topright",
           col=c("blue","red","green"),
           lty=c(1,1,1),
           legend=c("lognormal","weibull","gamma"))
  }
  return(list(weibull=est.weibull,lnorm=est.lnorm,gamma=est.gamma))
}

compare.KM.and.distributions <- function(fits,data) {
  require(survival)
  km <- survfit(Surv(time,status)~1,data=data)
  plot(km,
       xlab="Time",
       ylab="Survival Function")
  plot(function(x) pweibull(x,shape=fits$weibull$estimate[1],scale=fits$weibull$estimate[2],
                            lower.tail=FALSE),
       xlim=c(0,max(km$time)),add=TRUE,col="green")
  plot(function(x) plnorm(x,meanlog=fits$lnorm$estimate[1],sdlog=fits$lnorm$estimate[2],
                          lower.tail=FALSE),
       xlim=c(0,max(km$time)),add=TRUE,col="blue")
  plot(function(x) pgamma(x,shape=fits$gamma$estimate[1],rate=fits$gamma$estimate[2],
                          lower.tail=FALSE),
       xlim=c(0,max(km$time)),add=TRUE,col="orange")
  legend("topright",lty=c(1,2,1,1,1),col=c("black","black","green","blue","orange"),
         legend=c("Kaplan-Meier Estimate","95% Confidence Interval",
                  paste("Weibull",round(fits$weibull$loglik,5)),
                  paste("Lognormal",round(fits$lnorm$loglik,5)),
                  paste("Gamma",round(fits$gamma$loglik,5))),
         cex=0.5)
}

compare.KM.and.distributions.2 <- function(data) {
  require(survival)
  km <- survfit(Surv(time,status)~1,data=data)
  fits <- compare.histgram.and.densities(data,hist=FALSE,plot=FALSE)
  survreg.weibull <- survreg(Surv(time, status)~1, data=data, dist="weibull")
  fits$weibull$estimate[1] <- 1/survreg.weibull$scale
  fits$weibull$estimate[2] <- exp(survreg.weibull$coefficients[1])
  survreg.lognormal <- survreg(Surv(time, status)~1, data=data, dist="lognormal")
  fits$lnorm$estimate[1] <- survreg.lognormal$coefficients[1]
  fits$lnorm$estimate[2] <- survreg.lognormal$scale
  plot(km,
       xlab="Time",
       ylab="Survival Function")
  plot(function(x) pweibull(x,shape=fits$weibull$estimate[1],scale=fits$weibull$estimate[2],
                            lower.tail=FALSE),
       xlim=c(0,max(km$time)),add=TRUE,col="green")
  plot(function(x) plnorm(x,meanlog=fits$lnorm$estimate[1],sdlog=fits$lnorm$estimate[2],
                          lower.tail=FALSE),
       xlim=c(0,max(km$time)),add=TRUE,col="blue")
  legend("topright",lty=c(1,2,1,1,1),col=c("black","black","green","blue","orange"),
         legend=c("Kaplan-Meier Estimate","95% Confidence Interval",
                  paste("Weibull",round(fits$weibull$loglik,5)),
                  paste("Lognormal",round(fits$lnorm$loglik,5))),
         cex=0.5)
}