## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", message = FALSE,
                      warning = FALSE, fig.width = 6.2, fig.height = 5,
                      fig.align = "center")
set.seed(1)

## ----install, eval = FALSE----------------------------------------------------
# install.packages("gcmrec")

## ----lib----------------------------------------------------------------------
library(gcmrec)

## ----quickstart, eval = FALSE-------------------------------------------------
# data(readmission)
# 
# fit <- gcmrec(Survr(id, time, event) ~ as.factor(dukes) + sex,
#               data = readmission, s = 3000)
# 
# summary(fit)     # hazard ratios with confidence intervals
# plot(fit)        # baseline survivor function

## ----data---------------------------------------------------------------------
data(readmission)
head(readmission, 4)

length(unique(readmission$id))              # patients
table(table(readmission$id) - 1)            # events per patient

## ----survr--------------------------------------------------------------------
head(Survr(readmission$id, readmission$time, readmission$event), 4)

## ----addcentime---------------------------------------------------------------
dat <- data.frame(id    = c(1, 1, 2, 2),
                  time  = c(5, 3, 7, 4),
                  event = c(1, 0, 1, 1))     # subject 2 ends on an event
addCenTime(dat)

## ----hydraulic----------------------------------------------------------------
data(hydraulic)
head(List.to.Dataframe(hydraulic), 3)

## ----explore, fig.height = 5, fig.cap = "Rehospitalisations of the first 40 patients. Each row is a patient, each dot a readmission, the cross the end of follow-up."----
graph.caltimes(readmission[readmission$id %in% 1:40, ])

## ----mcf, fig.height = 4.2, fig.cap = "Mean cumulative function by Dukes' stage: the expected number of readmissions per patient."----
m <- mcf(readmission, group = readmission$dukes)
m
plot(m)

## ----fit----------------------------------------------------------------------
fit <- gcmrec(Survr(id, time, event) ~ as.factor(dukes) + sex,
              data = readmission, s = 3000)
fit

## ----summary------------------------------------------------------------------
summary(fit)

## ----forest, fig.height = 3.2, fig.cap = "Hazard ratios with 95% confidence intervals. Intervals crossing the dashed line are compatible with no effect."----
plotForest(fit, labels = c("as.factor(dukes)2" = "Dukes C vs A-B",
                           "as.factor(dukes)3" = "Dukes D vs A-B",
                           "sex" = "Female vs male"))

## ----extractors---------------------------------------------------------------
coef(fit)
sqrt(diag(vcov(fit)))       # standard errors
logLik(fit)
AIC(fit)

## ----anova--------------------------------------------------------------------
anova(fit)

## ----anova-nested-------------------------------------------------------------
fit.small <- gcmrec(Survr(id, time, event) ~ as.factor(dukes),
                    data = readmission, s = 3000)
anova(fit.small, fit)

## ----plot-baseline, fig.height = 4.2, fig.cap = "Baseline survivor function on the effective age scale, with its 95% confidence band."----
plot(fit)

## ----plot-styled, fig.height = 4.2, fig.cap = "The same curve, restyled."-----
plot(fit, type.plot = "hazard", level = 0.99) +
  ggplot2::labs(title = "Baseline cumulative hazard",
                subtitle = "Colorectal cancer readmissions")

## ----predict, fig.height = 4.2, fig.cap = "Predicted survivor function of the next readmission, by Dukes' stage (men)."----
profiles <- data.frame(dukes = c(1, 2, 3), sex = 1)

predict(fit, profiles, type = "risk")      # relative risk of each profile

plotPredict(fit, profiles,
            labels = c("Dukes A-B", "Dukes C", "Dukes D"))

## ----effage-models------------------------------------------------------------
mod.per <- gcmrec(Survr(id, time, event) ~ as.factor(dukes) + sex,
                  data = readmission, s = 3000, typeEffage = "perfect")
mod.min <- gcmrec(Survr(id, time, event) ~ as.factor(dukes) + sex,
                  data = readmission, s = 3000, typeEffage = "minimal")

rbind(perfect = coef(mod.per), minimal = coef(mod.min))

## ----effage-compare, fig.height = 4.2, fig.cap = "Baseline survivor function under perfect and minimal repair."----
plotBaseline(list(perfect = mod.per, minimal = mod.min))

## ----lymphoma-----------------------------------------------------------------
data(lymphoma)
table(lymphoma$effage)

mod.can <- gcmrec(Survr(id, time, event) ~ as.factor(distrib),
                  data = lymphoma, s = 1000, cancer = lymphoma$effage)
coef(mod.can)

## ----effagedata---------------------------------------------------------------
data(GeneratedData)
dat.gen <- List.to.Dataframe(GeneratedData)

mod.eff <- gcmrec(Survr(id, time, event) ~ covar.1 + covar.2,
                  data = dat.gen, effageData = GeneratedData, s = 100)
coef(mod.eff)

## ----frailty------------------------------------------------------------------
mod.fra <- gcmrec(Survr(id, time, event) ~ as.factor(dukes) + sex,
                  data = readmission, s = 3000, Frailty = TRUE)
coef(mod.fra)
mod.fra$Xi                        # frailty parameter

## ----frailty-values, fig.height = 3.8, fig.cap = "Estimated frailties. Patients to the right accumulate events faster than their covariates predict."----
summary(mod.fra$frailties)

ggplot2::ggplot(data.frame(z = mod.fra$frailties),
                ggplot2::aes(x = z)) +
  ggplot2::geom_histogram(bins = 30, fill = "#0072B2", alpha = 0.85) +
  ggplot2::labs(x = "Estimated frailty", y = "Patients") +
  theme_gcmrec()

## ----jackknife----------------------------------------------------------------
sub <- readmission[readmission$id %in% unique(readmission$id)[1:60], ]

mod.info <- gcmrec(Survr(id, time, event) ~ as.factor(dukes),
                   data = sub, s = 3000)
mod.jack <- gcmrec(Survr(id, time, event) ~ as.factor(dukes),
                   data = sub, s = 3000, se = "Jacknife")

rbind(information = sqrt(diag(vcov(mod.info))),
      jackknife   = sqrt(diag(vcov(mod.jack))))

## ----sessioninfo, echo = FALSE------------------------------------------------
sessionInfo()