library(survival)
setwd("D:\\TEMP")
dat <- read.table("addict.csv", header = TRUE, sep = ",")
Set contrasts for clinic and prison. Set the reference category for clinic, making Clinic 1 (base = 1) the reference category:
dat$clinic <- factor(dat$clinic, levels = c(1, 2), labels = c("1", "2"))
contrasts(dat$clinic) <- contr.treatment(2, base = 1, contrasts = TRUE)
levels(dat$clinic)
Same for prison, making absence of a prison record the reference category:
dat$prison <- factor(dat$prison, levels = c(0, 1), labels = c("0", "1"))
contrasts(dat$prison) <- contr.treatment(2, base = 1, contrasts = TRUE)
levels(dat$prison)
Cox proportional hazards model (for comparison):
addict.cph01 <- coxph(Surv(stop, status, type = "right") ~ clinic + prison + dose, method = "breslow", data = dat)
summary(addict.cph01)
Exponential model:
addict.exp01 <- survreg(Surv(stop, status, type = "right") ~ clinic + prison + dose, dist = "exponential", data = dat)
summary(addict.exp01)
shape.exp = 1 / addict.exp01$scale
shape.exp
R outputs the parameter estimates for the accelerated failure time form of the exponential model. We multiply the estimated regression coefficients by -1 to get estimates consistent with the Cox proportional hazards parameterisation. For example, the estimated accleration factor comparing a patient with a prison record to one without a prison record is exp(-0.2526) = 0.78. Having a prison record accelerates time to event by a factor of 0.78. The estimated hazard ratio comparing a patient with a prison record to one without a prison record is exp(-1 × 0.2526) = 1.29. Having a prison record increases the daily hazard of relapse by a factor of 1.29. We now run a Weibull model:
addict.wei01 = survreg(Surv(stop, status, type = "right") ~ clinic + prison + dose, dist = "weibull", data = dat)
summary(addict.wei01)
The Weibull shape parameter is the reciprocal of the variable called scale in the model output:
shape.wei = 1 / addict.wei01$scale
shape.wei
The acceleration factor comparing clinic 2 with clinic 1 patients is exp(0.8806) = 2.41. The estimated median survival time (time off heroin) is double for patients from clinic 2 compared with those from clinic 1. We now plot the estimated survival for patients attending clinic 1 and clinic 2. Both have a prison record and received a maximum dose of 50 mg methadone per day:
addict.km01 <- survfit(Surv(stop, status) ~ clinic, type = "kaplan-meier", data = dat)
dat.clin01 <- data.frame(clinic = as.factor(c(1)), prison = as.factor(c(0)), dose = 50)
dat.clin02 <- data.frame(clinic = as.factor(c(2)), prison = as.factor(c(0)), dose = 50)
days.clin01 <- predict(addict.wei01, newdata = dat.clin01, type = "quantile", p = 0:100/100)
days.clin02 <- predict(addict.wei01, newdata = dat.clin02, type = "quantile", p = 0:100/100)
surv.wei01 <- 1 - 0:100/100
plot(addict.km01, xlab = "Days to relapse", ylab = "Cumulative proportion to experience event", lty = c(1,2), mark.time = FALSE, main = "")
lines(days.clin01, surv.wei01, type = "s", col = "red", lty = 1)
lines(days.clin02, surv.wei01, type = "s", col = "red", lty = 2)
legend(x = "topright", legend = c("Clinic 1: Kaplan Meier","Clinic 2: Kaplan Meier", "Clinic 1: Weibull model", "Clinic 2: Weibull model"), lty = c(1,2,1,2), col = c("black", "black", "red", "red"), bty = "n", cex = 0.80)
Now compare the three models using AIC:
extractAIC(addict.cph01)
extractAIC(addict.exp01)
extractAIC(addict.wei01)
The AIC for the Cox model is the smallest, indicating that this model provides the best fit with the data (this is consistent with the diagnostics we ran earlier to assess how consistent the data was with the exponential and Weibull distributions).
Here we use the psm function in the rms library to develop an AFT model. The psm function is a modification of survreg and is used for fitting the accelerated failure time family of parametric survival models.
library(rms)
addict.aft01 <- psm(Surv(stop, status) ~ clinic + prison + dose, dist = "weibull", data = dat)
addict.aft01
Compare addict.aft01 with addict.wei01:
addict.aft01
summary(addict.wei01)
What is the effect of Clinic 2 on retention time (after adjusting for the effect of presence of a prison record and maximum daily methadone dose)?
exp(addict.aft01$coefficients[2])
Treatment at Clinic 2 doubles patient retention time. What does this mean in terms of calendar time?
log.t <- as.numeric(addict.aft01$coefficients[1]) + (as.numeric(addict.aft01$coefficients[2]) * 1)
exp(log.t)
Treatment at Clinic 2 results in patients remaining on the program for an additional 250 days (compared with those treated at Clinic 1).