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
newbili <- seq(min(pbc3$bili), max(pbc3$bili),length.out = 500)
newlogbili <- seq(min(pbc3$logbili), max(pbc3$logbili),length.out = 500)

## function for calulating pseudo residuals for a model using pbc3 data
pseudoresiduals <- function(t, fit){
  ## pseudobservations
  p <- as.vector(pseudosurv(pbc3$followup, pbc3$event, tmax = t)[,3])
  ## predictions
  s <- as.vector(summary(fit, times = t)$surv)
  ## pseudo-residuals
  (p-s)/sqrt(s*(1-s))
}

fit <- survfit(coxph(Surv(followup, event) ~ tment, data=pbc3),
               newdata = data.frame(tment = pbc3$tment),
               type = "kaplan-meier")

offset <- c(5/4, 7/4)
## common range on y-axis
range <- range(mapply(function(t){pseudoresiduals(t, fit)}, times))

par(mfrow=c(2,2))
for (t in times){
  residuals <- pseudoresiduals(t, fit)
  stripchart(residuals ~ pbc3$tment,
             vertical = TRUE,
             pch = 1,
             ylim = range,
             ylab = paste("Pseudo-residual, time=", round(t,2)),
             group.names = c("Placebo", "CyA"),
             at = offset,
             col = "grey"
             )
  abline(h = 0,lty = "21")
  points(offset[1], mean(residuals[pbc3$tment==0]), pch=4, cex=2)
  points(offset[2], mean(residuals[pbc3$tment==1]), pch=4, cex=2)
}