#+ knitr_setup, include = FALSE

# Whether to cache the intensive code sections. Set to FALSE to recalculate
# everything afresh.
cacheoption <- FALSE
# Disable lazy caching globally, because it fails for large objects, and all the
# objects we wish to cache are large...
opts_chunk$set(cache.lazy = FALSE)

#' # Performance of an age-only model
#' 
#+ setup, message=FALSE

bootstraps <- 200

data.filename <- '../../data/cohort-sanitised.csv'
n.threads <- 16

source('../lib/shared.R')
require(rms)

COHORT.full <- data.frame(fread(data.filename))

COHORT.use <- subset(COHORT.full, !exclude)

n.data <- nrow(COHORT.use)

# Define indices of test set
test.set <- testSetIndices(COHORT.use, random.seed = 78361)

COHORT.use <- prepSurvCol(COHORT.use, surv.time, surv.event,surv.event.yes)

#' ## Concordance index
#' 
#' The c-index can be calculated taking age itself as a risk score. Since it's
#' purely rank-based, it doesn't matter that age is nonlinearly related to true
#' risk.
#' 
#' We bootstrap this based purely on the test set to make the bootstrap
#' variability commensurate with other models tested on the test set. You can
#' imagine that this model was 'trained' on the training set, even though it's
#' so simple that we actually did no such thing...
#' 
#+ c_index

# Define a trivial function used for bootstrapping which simply returns the
# c-index of a 'model' based purely on age.

ageOnly <- function(df, indices, df.test) {
  # Create a Kaplan-Meier curve from the bootstrap sample
  km.by.age <-
    survfit(
      Surv(surv_time, surv_event) ~ age,
      data = df[indices, ],
      conf.type = "log-log"
    )

  # Return the calibration score
  c(
    calibration.score =
      calibrationScore(
        calibrationTable(km.by.age, df.test)
      )$area
  )
}

# C-index is by definition non-variable because we do not bootstrap or otherwise
# permute the test set, and there's no 'training' phase
age.c.index <- 
  as.numeric(
    survConcordance(
      Surv(surv_time, surv_event) ~ age,
      COHORT.use[test.set,]
    )$concordance
  )

# Bootstrap to establish variability of calibration
age.only.boot <-
  boot(
    data = COHORT.use[-test.set,],
    statistic = ageOnly,
    R = bootstraps,
    parallel = 'multicore',
    ncpus = n.threads,
    df.test =  COHORT.use[test.set,]
  )

age.only.boot.stats <- bootStats(age.only.boot, uncertainty = '95ci')

#' C-index is
#' **`r round(age.c.index, 3)`** 
#' on the held-out test set (not that it really matters, the model isn't
#' 'trained' as such for the discrimination test...it's just oldest patient dies
#' first).
#' 
#' 
#' ## Calibration
#' 
#' 
#+ calibration

km.by.age <-
  survfit(
    Surv(surv_time, surv_event) ~ age,
    data = COHORT.use[-test.set,],
    conf.type = "log-log"
  )

calibration.table <- calibrationTable(km.by.age, COHORT.use[test.set, ])

print(calibrationScore(calibration.table))

calibrationPlot(calibration.table)

#' Calibration score is
#' **`r round(age.only.boot.stats['calibration.score', 'val'], 3)`
#' (`r round(age.only.boot.stats['calibration.score', 'lower'], 3)` -
#' `r round(age.only.boot.stats['calibration.score', 'upper'], 3)`)** 
#' on the held-out test set.

varsToTable(
  data.frame(
    model = 'age',
    imputation = FALSE,
    discretised = FALSE,
    c.index = age.c.index,
    c.index.lower = NA,
    c.index.upper = NA,
    calibration.score = age.only.boot.stats['calibration.score', 'val'],
    calibration.score.lower = age.only.boot.stats['calibration.score', 'lower'],
    calibration.score.upper = age.only.boot.stats['calibration.score', 'upper']
  ),
  performance.file,
  index.cols = c('model', 'imputation', 'discretised')
)