#' Codes used in the manuscript
#'
#' reReg version 1.4.4 requires the latest reda, which can be installed via
#' remotes::install_github('wenjie2wang/reda', upgrade = 'never')

## Section 2
library("reReg")
library("Hmisc")
subset(simDat, id %in% c(5, 7, 12))

## Section 5.1
args(reda::Recur)
(reObj <- with(simDat, Recur(t.start %2% t.stop, id, event, status)))
summary(reObj)

## Section 5.2
args(plotEvents)
args(reReg::plotEvents.control)
plot(reObj)
plotEvents(Recur(t.start %to% t.stop, id, event, status) ~ 1, data = simDat)
plotEvents(Recur(t.start %to% t.stop, id, event, status) ~ x1, data = simDat)
simDat$event2 <- with(simDat, ifelse(t.stop > 10 & event > 0, 2, event))
plotEvents(Recur(t.start %to% t.stop, id, event2, status) ~ x1, data = simDat)
simDat2 <- simDat
simDat2$t.start <- as.Date(simDat2$t.start + simDat2$x2 * 5, origin = "20-01-01")
simDat2$t.stop <- as.Date(simDat2$t.stop + simDat2$x2 * 5, origin = "20-01-01")
plotEvents(Recur(t.start %to% t.stop, id, event, status) ~ x1, data = simDat2, calendarTime = TRUE)

## Section 5.3
mcf1 <- reReg(Recur(t.start %to% t.stop, id, event, status) ~ 1, data = simDat, B = 200)
plot(mcf1)
args(basebind)
mcf2 <- reReg(Recur(t.start %to% t.stop, id, event, status) ~ 1, subset = x1 == 0,
  data = simDat, B = 200)
mcf3 <- update(mcf2, subset = x1 == 1)
g1 <- plot(mcf2)
g2 <- plot(mcf3)
basebind(g1, g2, legend.title = "X1", legend.labels = 0:1)
plot(reObj, mcf = TRUE, mcf.conf.int = TRUE)
mcf0 <- mcf(Recur(t.start %to% t.stop, id, event, status) ~ x1, data = simDat)
plot(mcf0, conf.int = TRUE)

## Section 5.4
args(reReg)
args(reReg.control)
fm <- Recur(t.stop, id, event, status) ~ x1 + x2
set.seed(0)
system.time(fit1 <- reReg(fm, data = simDat, model = "cox", B = 200))
summary(fit1)
set.seed(0)
system.time(fit2 <- reReg(fm, data = simDat, model = "cox|cox", B = 200))
summary(fit2)
argsAnywhere(plot.reReg)
plot(fit2)
newdata <- expand.grid(x1 = 0:1, x2 = mean(simDat$x2))
plot(fit2, newdata = newdata, showName = TRUE)
set.seed(0)
system.time(fit3 <- reReg(fm, data = simDat, model = "gsc", se = "sand", B = 200))
summary(fit3, test = TRUE)

## Section 5.5
summary(reReg(fm, data = simDat, model = "cox.LWYY"))
library("survival")
summary(coxph(Surv(t.start, t.stop, event) ~ x1 + x2 + cluster(id), data = simDat))$coef

## Section 5.6
args(simGSC)
data(simDat, package = "reReg")
set.seed(0)
dat <- simGSC(200, summary = TRUE)
identical(dat, simDat)

fm <- Recur(t.stop, id, event, status) ~ x1 + x2
sizes <- c(100, 200, 400, 600, 800, 1000)
set.seed(0)
times <- sapply(sizes, function(n) {
  rowMeans(replicate(50, {
    dat <- simGSC(n)
    c(system.time(reReg(fm, data = dat, model = "cox|cox", B = 0))[3], system.time(reReg(fm,
      data = dat, model = "am|am", B = 0))[3], system.time(reReg(fm, data = dat,
      model = "gsc", B = 0))[3])
  }))
})
rownames(times) <- c("\\texttt{model = cox|cox}", "\\texttt{model = am|am}", "\\texttt{model = gsc}")
knitr::kable(round(times, 3), col.names = paste0("n = ", sizes), caption = "Runtime (in seconds) under different sample sizes.",
  row.names = TRUE, booktabs = TRUE, escape = FALSE)

## Section 6
data(colorectal, package = "frailtypack")
str(colorectal)
(colreObj <- with(colorectal, Recur(time0 %to% time1, id, new.lesions, state)))
summary(colreObj)
fm <- Recur(time0 %to% time1, id, new.lesions, state) ~ treatment
plotEvents(fm, data = colorectal, base_size = 8, cex = 0.8, recurrent.name = "New lesions",
  terminal.name = "Death")
plot(mcf(fm, data = colorectal), conf.int = T)
set.seed(0)
fm <- update(fm, ~. + age + who.PS + prev.resection)
fitCox <- reReg(fm, data = colorectal, model = "cox|cox", B = 200)
summary(fitCox)
set.seed(0)
fitGSC <- reReg(fm, data = colorectal, model = "gsc", se = "sand", B = 200)
summary(fitGSC, test = TRUE)