Example code for time-to-event analysis in R, as in whether repeated ad viewings lead to a sale

timetoevent.R

##
## Example code for time-to-event analysis in R
## [email protected]
## Dec 28, 2012
##
## joineR package: analyzing longitudinal data where the response
## from each person is a time-sequence of repeated measurements 
## and we are interested in a possibly censored time-to-event outcome
##
## example: repeated ad viewings leading to a sale

library(joineR)   # main functions
library(survival) # for plotting


data(mental) # from joineR

mental[1:5,] # examine

## convert balanced (wide) to unbalanced (tall) format

mental.unbalanced <- to.unbalanced(mental, id.col = 1, times = c(0, 1, 2, 4, 6, 8),
Y.col = 2:7, other.col = 8:11)
dim(mental.unbalanced)
names(mental.unbalanced)
names(mental.unbalanced)[3] <- "Y"

mental.balanced <- to.balanced(mental.unbalanced, id.col = 1, time.col = 2,
Y.col = 3, other.col = 4:7)
dim(mental.balanced)
names(mental.balanced)

## get summary for a bal object

summarybal(mental, Y.col = 2:7, times = c(0, 1, 2, 4, 6, 8), na.rm = TRUE)

## explore covariance structure for balanced data

y <- as.matrix(mental[, 2:7]) 

# convert mental from list format to numeric matrix format
means <- matrix(0,3,6)
for (trt in 1:3) {
   ysub <- y[mental$treat == trt,]
   means[trt,] <- apply(ysub, 2, mean, na.rm = TRUE)
   }
residuals <- matrix(0, 150, 6)
for (i in 1:150) {
   residuals[i,] <- y[i,] - means[mental$treat[i],]
   }
V <- cov(residuals, use = "pairwise")
R <- cor(residuals, use = "pairwise")
round(cbind(diag(V), R), 3)


## explore covariance structure for UNbalanced data

vgm <- variogram(indv = mental.unbalanced[, 1], time = mental.unbalanced[, 2], 
 Y = mental.unbalanced[, 3])
vgm$sigma2
par(mfrow = c(1, 3))
plot(vgm$svar[, 1], vgm$svar[, 2], pch = 19, cex = 0.5, xlab = "u", ylab = "V(u)")
plot(vgm, points = FALSE, ylim = c(0, 200)) 
plot(vgm)


###################################################
# Formal modelling  


# make joindata object

mental.long <- mental.unbalanced[, 1:3]
mental.surv <- UniqueVariables(mental.unbalanced, 6:7, id.col = 1)
mental.baseline <- UniqueVariables(mental.unbalanced, 4, id.col = 1)
mental.jd <- jointdata(mental.long, mental.surv, mental.baseline, id.col = "id", 
 time.col = "time")


# make model with joint()

# joint model with random latent association, based on Wulfsohn and Tsiatis (1997) 
# allowing for longitudinal and survival covariates
# three choices for latent process: 
#    model = c("intslope",  # joint rand effects model with join intercept and slope 
#                           #      as per Wulfsohn and Tsiatis (default) 
#              "int",       # intercept only
#              "quad")      # random quadratic for repeated meas submodel

model.jointrandom <- joint(data=mental.jd, 
                          Y ~ 1 + time + treat,  # longitudinal submodel
                          Surv(surv.time, cens.ind) ~ treat, # hazard submodel
                          model = "int")

summary(model.jointrandom)


# get standard error and CIs for mean response in random efffects 
# this is SLOW (2 min?)

model.jointrandom.se <- jointSE(model.jointrandom, n.boot = 100) # use 100+
model.jointrandom.se

######################
# un-connected samples of plotting

#fit a Kaplan-Meier and plot it 

fit <- survfit(Surv(time, status) ~ x, data = survival::aml) 
plot(fit, lty = 2:3) 
legend(70, .8, c("Maintained", "Nonmaintained"), lty = 2:3) 

#fit a Cox proportional hazards model and plot the  
#predicted survival for a 60 year old 

fit <- coxph(Surv(futime, fustat) ~ age, data = ovarian) 
plot(survfit(fit, newdata=data.frame(age=60)),
     xscale=365.25, xlab = "Years", ylab="Survival")

still need to make plots to visualize results