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

library("survival")

## Recall levels of the status variable dc: 2 is censored, 1 and 3 are deaths.
coxmodel <- coxph(Surv(melanoma$years, melanoma$dc != 2) ~ sex + ulc + thick + age,
                  data = melanoma)

sf <- survfit(coxmodel, newdata =
              data.frame(sex = as.factor(0),
                         age = 50,
                         thick = 3,
                         ulc = as.factor(0)),
              conf.type = "log-log"
              )
S <- sf$surv
H <- -log(S)
plot(log(H)~log(sf$time),
     ylim = c(-9.5, -1),
     type = "s",
     xlab = "log(years)",
     ylab = "log(cumulative hazard)")
lines(log(-log(sf$upper))~log(sf$time),
      type = "s",
      lty = "21")
lines(log(-log(sf$lower))~log(sf$time),
      lty = "21",
      type = "s")

## inner, extra x-axis in year-scale
ticks <- c(0.05,0.1,0.5,1,5)
axis(1, at = log(ticks), labels = ticks, tcl = 0.3, cex.axis=0.8, mgp = c(-3,-1.2,0))