source("http://192.38.117.59/~linearpredictors/datafiles/readPbc3.R")

library(survival)

plot.integr.intensities <- function(fit, pch=1, ref=1, xlab = "", ylab = "", add=TRUE,...){
  extended.integr.intensities <- function(fit, i, ...){
    times <- sort(fit$time)
    -log(
         approx(fit[i]$time, fit[i]$surv, 
                times,
                rule=2, # use the closest data extreme outside range 
                method="constant", ...)$y
         )
  }

  draw.base <- function(fit,i,...){
    ## calculating default ranges
    xrange <- c(0, max(-log(fit[ref]$surv)))
    yrange <- c(0, max(-log(fit$surv)))
    
    ## basis plot
    plot(0,0, type="n",
         xlim = xrange,       
         ylim = yrange,
         xlab = "Cumulative hazard, 1. quintile group",
         ylab = paste("Cumulative hazard,", i,". quintile group")
    )
  }
  ngroups <- length(fit$strata)


  i <- 1
  j <- 1 # counter for plotting symbol
  while (i <= ngroups)
    {

      if (i != ref) {
              draw.base(fit,i,...)
          points(extended.integr.intensities(fit,ref),
                 extended.integr.intensities(fit,i),
                 pch = pch[j], ...)
        j <- j + 1
      }
      i <- i + 1
    }
}


## Now, the integrated intensities for Pbc3 data.

par(mfrow = c(2,2))
dataset <- subset(pbc3, followup>0)
pvec <- seq(0,1,0.2)
dataset$biligrp <- cut(dataset$bili,quantile(dataset$bili,pvec))

fit <- survfit(Surv(dataset$followup, dataset$status != 0) ~ dataset$biligrp,
               type = "fleming-harrington")

plot.integr.intensities(fit, type="s", add=FALSE,
                        xlab = "Cumulated hazard, 1. quintile group",
                        ylab = paste("Cumulative hazard", )
                        )