source("http://192.38.117.59/~linearpredictors/datafiles/readMelanoma.R")
library("timereg")

## levels of the status variable dc: 2 is censored, 1 and 3 are deaths.
additivemodelNonConst <- aalen(Surv(melanoma$years, melanoma$dc != 2) ~
                               sex + ulc + I(thick-3) + I(age-50),
                               data = melanoma)
confbands <- function(term){
  se <- sqrt(additivemodelNonConst$var.cum[,term])
  lines(additivemodelNonConst$cum[,term] - 1.96*se ~ eventtimes,
        type = "s",
        lty = "21")
  lines(additivemodelNonConst$cum[,term] + 1.96*se ~ eventtimes,
        type = "s",
        lty = "21")
}

par(mfrow = c(2,2))
eventtimes <- (additivemodelNonConst$cum)[,1]
plot(eventtimes, (additivemodelNonConst$cum)[,3],
     type = "s",
     xlab = "Years",
     ylab = "Cumulative regression function, Sex")
confbands(3)

eventtimes <- (additivemodelNonConst$cum)[,1]
plot(eventtimes, (additivemodelNonConst$cum)[,4],
     type = "s",
     xlab = "Years",
     ylab = "Cumulative regression function, Ulceration")
confbands(4)
  
eventtimes <- (additivemodelNonConst$cum)[,1]
plot(eventtimes, (additivemodelNonConst$cum)[,5],
     type = "s",
     xlab = "Years",
     ylab = "Cumulative regression function, Thickness")
confbands(5)

eventtimes <- (additivemodelNonConst$cum)[,1]
plot(eventtimes, (additivemodelNonConst$cum)[,6],
     type = "s",
     xlab = "Years",
     ylab = "Cumulative regression function, Age")
confbands(6)