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

## Define log-variables, for convenience
pbc3$logbili <- log(pbc3$bili)
pbc3$logasptr <- log(pbc3$asptr)
pbc3$logalkph <- log(pbc3$alkph)

## Recoding the status variable for use with package "pseudo"
pbc3$event <- ifelse(pbc3$status != 0, 1, 0)

## Model of interest
cox.model <- coxph(Surv(followup, event) ~ tment + logbili + alb + logalkph + logasptr + age,
                   data = pbc3,
                   model = T) ## return the data, to which the model is actually fitted!

## The dataset actually used
used.data <- na.omit(get_all_vars(cox.model, pbc3))

fit <- survfit(cox.model,
               newdata = cox.model$model,
               type = "kaplan-meier")

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

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

  ## pseudo residuals vs log of bilirubin

par(mfrow = c(2,2))
range <- range(mapply(function(t){pseudoresiduals(t)}, times))

for (t in times){
  used.data$pseudores <- NULL
  used.data$pseudores <- pseudoresiduals(t)
  plot(used.data$pseudores ~ used.data$logalkph,
       ylim = range, 
       xlab = "log(alkaline phosphatase)",
       ylab = paste("Pseudo-residual, time=",round(t,2)),
       col = "grey")
  
  ## New alkph values in which to evaluate the smoother
  newlogalkph <- seq(min(used.data$logalkph),max(used.data$logalkph),length.out = 500)
  
  smoothed <- predict(loess(pseudores ~ logalkph,
                            data = used.data, degree = 1), 
                      data.frame(logalkph = newlogalkph))
  lines(smoothed ~ newlogalkph)
  abline(h = 0, lty = "21")
}