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

pbc3$logbili <- log(pbc3$bili)
pbc3$event <- ifelse(pbc3$status != 0, 1, 0)

## we select 4 time points corresponding to quintiles of observed event times
times <- quantile(pbc3$followup[pbc3$event == 1], c(.2,.4,.6,.8))  

## bili values in which to evaluate the smoother
newlogbili <- seq(min(pbc3$logbili), max(pbc3$logbili),length.out = 500)

## scatter plot: pseudo observations vs log of bilirubin
## Transformed by log(-log). Added piecewise constant linear predictor

## Quintile groups of bilirubin with group means as levels
pbc3$bgm <- pbc3$biligroup
biliQuintileMean <- tapply(pbc3$bili, pbc3$biligroup, mean)
levels(pbc3$bgm) <- biliQuintileMean

## Fitting Cox model
coxphModel <- coxph(Surv(followup, status != 0) ~ bgm, data=pbc3)
loghr <- c(0,coxphModel$coefficients) # added 0 as baseline level

## survival for first quintile group
sfit <- survfit(coxphModel,
                newdata = data.frame(bgm = as.factor(biliQuintileMean)[1])
                )

par(mfrow = c(2,2))      
for (t in times){
  
  ## survival for first quintile group calculated at time t
  ssfit <- summary(sfit, time = t)    
  cloglogbase <- log(-log(ssfit$surv))
  cloglog <- c(loghr, loghr[5]) + cloglogbase ## unfortunate notation
  
  plot(log(quintiles), cloglog, 
       type="s",
       ylim = c(-6,2),
       xlab = "log(bilirubin)",
       ylab = paste("time =",round(t,2)))
  
  ## smoother of pseudo observations
  pbc3$pseudo <- NULL
  pbc3$pseudo <- pseudosurv(pbc3$followup, pbc3$event, tmax = t)[,3]
  smoothed <- predict(loess(pseudo ~ logbili,
                            data = pbc3, degree = 1),
                      data.frame(logbili = newlogbili))
  
  ## smoothed observations added
  lines(log(-log(smoothed))~newlogbili) 
}