#+ knitr_setup, include = FALSE
# Whether to cache the intensive code sections. Set to FALSE to recalculate
# everything afresh.
cacheoption <- TRUE
# 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)
#' # Modelling with age
#'
#' It seems that, no matter what I do, the C-index of a model, random forest or
#' otherwise, is about 0.78. I decided to try to some simpler models, intially
#' based purely on age which is clearly the biggest factor in this dataset.
#' (And, I assume, most datasets where a reasonably broad range of ages is
#' present.)
#'
#' The sanity check works: giving the model more data does indeed result in a
#' better fit. However, on top of that, I was surprised by just how good the
#' performance can be when age alone is considered!
#+ user_variables, message=FALSE
data.filename <- '../../data/cohort-sanitised.csv'
n.trees <- 500
continuous.vars <-
c(
'age', 'total_chol_6mo', 'hdl_6mo', 'pulse_6mo', 'crea_6mo',
'total_wbc_6mo', 'haemoglobin_6mo'
)
untransformed.vars <- c('anonpatid', 'time_death', 'imd_score', 'exclude')
source('../lib/shared.R')
require(ggrepel)
#' ## Load and prepare data
#+ load_and_prepare_data
# Load the data and convert to data frame to make column-selecting code in
# prepData simpler
COHORT.full <- data.frame(fread(data.filename))
# Define process settings; nothing for those to not transform, and missingToBig
# for the continuous ones...
process.settings <-
list(
var = c(untransformed.vars, continuous.vars),
method =
c(
rep(NA, length(untransformed.vars)),
rep('missingToBig', length(continuous.vars))
),
settings = rep(NA, length(untransformed.vars) + length(continuous.vars))
)
COHORT.prep <-
prepData(
# Data for cross-validation excludes test set
COHORT.full,
cols.keep,
process.settings,
surv.time, surv.event,
surv.event.yes,
extra.fun = caliberExtraPrep
)
n.data <- nrow(COHORT.prep)
# Define indices of test set
test.set <- sample(1:n.data, (1/3)*n.data)
#' ## Models
#'
#' ### The normal model
#'
#' All the variables, as in the vector `surv.predict`.
#+ normal_model, cache=cacheoption
# Fit random forest
surv.model.fit <-
survivalFit(
surv.predict,
COHORT.prep[-test.set,],
model.type = 'rfsrc',
n.trees = n.trees,
split.rule = 'logrank',
n.threads = 7,
nsplit = 20
)
print(surv.model.fit)
# Get C-index
c.index.test <-
cIndex(surv.model.fit, COHORT.prep[test.set, ], model.type = 'ranger')
#' The C-index on the held-out test set is **`r round(c.index.test, 3)`**.
#' ### Just age
#'
#' What if all we had to go on was age?
#+ just_age_model, cache=cacheoption
# Fit random forest
surv.model.fit <-
survivalFit(
c('age'),
COHORT.prep[-test.set,],
model.type = 'ranger',
n.trees = n.trees,
split.rule = 'logrank',
n.threads = 8,
respect.unordered.factors = 'partition'
)
print(surv.model.fit)
# Get C-index
c.index.test <-
cIndex(surv.model.fit, COHORT.prep[test.set, ], model.type = 'ranger')
#' The C-index on the held-out test set is **`r round(c.index.test, 3)`**.
#' ### No model, literally just age
#'
#' What if we constructed the C-index based purely on patients' ages?
#+ just_age_cindex
c.index.age <-
as.numeric(
survConcordance(
Surv(time_death, surv_event) ~ age,
COHORT.prep
)$concordance
)
#' The C-index on the whole dataset based purely on age is
#' **`r round(c.index.test, 3)`**. That's most of our predictive accuracy right
#' there! Reassuringly, it's also equal to the value predicted by the random
#' forest model based purely on age...
#' ### Age and gender
#'
#' OK, age and gender.
#+ age_gender_model, cache=cacheoption
# Fit random forest
surv.model.fit <-
survivalFit(
c('age', 'gender'),
COHORT.prep[-test.set,],
model.type = 'ranger',
n.trees = n.trees,
split.rule = 'logrank',
n.threads = 8,
respect.unordered.factors = 'partition'
)
print(surv.model.fit)
# Get C-index
c.index.test <-
cIndex(surv.model.fit, COHORT.prep[test.set, ], model.type = 'ranger')
#' The C-index on the held-out test set is **`r round(c.index.test, 3)`**.
#' ### Age, gender and history of liver disease
#'
#' Let's add a third variable. In the replication of the Cox model with missing
#' data included, liver disease was the most predictive factor after age, so
#' it's a reasonable next variable to add.
#+ age_gender_liver_model, cache=cacheoption
# Fit random forest
surv.model.fit <-
survivalFit(
c('age', 'gender', 'hx_liver'),
COHORT.prep[-test.set,],
model.type = 'ranger',
n.trees = n.trees,
split.rule = 'logrank',
n.threads = 8,
respect.unordered.factors = 'partition'
)
print(surv.model.fit)
# Get C-index
c.index.test <-
cIndex(surv.model.fit, COHORT.prep[test.set, ], model.type = 'ranger')
#' The C-index on the held-out test set is **`r round(c.index.test, 3)`**.
#' ### Age, gender and heart failure
#'
#' A different third variable: heart failure, the second most important variable
#' (after age) from random forest modelling.
#+ age_gender_hf_model, cache=cacheoption
# Fit random forest
surv.model.fit <-
survivalFit(
c('age', 'gender', 'heart_failure'),
COHORT.prep[-test.set,],
model.type = 'ranger',
n.trees = n.trees,
split.rule = 'logrank',
n.threads = 8,
respect.unordered.factors = 'partition'
)
print(surv.model.fit)
# Get C-index
c.index.test <-
cIndex(surv.model.fit, COHORT.prep[test.set, ], model.type = 'ranger')
#' The C-index on the held-out test set is **`r round(c.index.test, 3)`**.
#' ### Just gender
#'
#' Just gender, as a sanity check.
#+ just_gender_model, cache=cacheoption
# Fit random forest
surv.model.fit <-
survivalFit(
c('gender'),
COHORT.prep[-test.set,],
model.type = 'ranger',
n.trees = n.trees,
split.rule = 'logrank',
n.threads = 8,
respect.unordered.factors = 'partition'
)
print(surv.model.fit)
# Get C-index
c.index.test <-
cIndex(surv.model.fit, COHORT.prep[test.set, ], model.type = 'ranger')
#' The C-index on the held-out test set is **`r round(c.index.test, 3)`**.
#' ### Everything except age
#'
#' How do we do if we use all the variables _except_ age?
#+ no_age_model, cache=cacheoption
# Fit random forest
surv.model.fit <-
survivalFit(
surv.predict[surv.predict != 'age'],
COHORT.prep[-test.set,],
model.type = 'ranger',
n.trees = n.trees,
split.rule = 'logrank',
n.threads = 8,
respect.unordered.factors = 'partition'
)
print(surv.model.fit)
# Get C-index
c.index.test.all.not.age <-
cIndex(surv.model.fit, COHORT.prep[test.set, ], model.type = 'ranger')
#' The C-index on the held-out test set is
#' **`r round(c.index.test.all.not.age, 3)`**.
#' ### Predicting age
#'
#' So why does the model which doesn't include age do so well? Clearly the other
#' variables allow you to predict age with reasonable accuracy... So let's try
#' just that as a final test.
#+ predict_age, cache=cacheoption
options(rf.cores = n.threads)
age.model <-
rfsrc(
formula(
paste0(
# Predicting just the age
'age ~ ',
# Predictor variables then make up the other side
paste(surv.predict[!(surv.predict %in% c('age', 'most_deprived'))], collapse = '+')
)
),
COHORT.use[-test.set, ],
ntree = n.trees,
splitrule = 'mse',
na.action = 'na.impute',
nimpute = 3
)
age.predictions <- predict(age.model, COHORT.prep[test.set, ])
age.cor <- cor(age.predictions$predictions, COHORT.prep[test.set, 'age'])
to.plot <-
data.frame(
age = COHORT.prep[test.set, 'age'],
predicted = age.predictions$predictions
)
ggplot(sample.df(to.plot, 10000), aes(x = age, y = predicted)) +
geom_point(alpha = 0.2)
#' It doesn't look that great, but there is some correlation...
#' r^2 = `r age.cor^2` which is unremarkable, but OK. A more relevant measure
#' would be the pure-age C-index, ie if I gave you a pair of patients, how
#' often could you predict who was older?
c.index.age.on.age <-
1 - as.numeric(
survConcordance(
age ~ predicted,
to.plot
)$concordance
)
#' This comes out as **`r round(c.index.age.on.age, 3)`**, as compared to
#' **`r round(c.index.test.all.not.age, 3)`**, which was the C-index for the
#' survival model based on all other variables.
#'
#' This makes sense. The maths of C-indices is a bit tricky (how much they add
#' to one-another depends in large part of how correlated the variables are, as
#' well as probably being somewhat non-linear anyway),
#' but clearly some significant fraction of the all-except-age model's
#' predictive power comes from its ability to infer age from the remaining
#' variables, and then use _that_ (implicitly) to predict time to death.
#'
#' The mechanism for this could be that people are likely to
#' get more diseases as they get older, so if you have a, b _and_ c you're
#' likely to be older. Second-order predictivity may occur if particular
#' combinations of disorders or test results are common in certain age groups.
#'
#' ## Conclusion
#'
#' So, in conclusion, the sanity check has worked: giving the
#' random forest model more data to work with improves its performance. Age
#' alone is less predictive, adding gender makes it slightly more predictive,
#' and so on.
#'
#' Further, though, a huge amount of the model's predictivity arises from just
#' a patien's age. Not only is that alone a good predictor, but the
#' reasonable performance on the model of all factors except age is explained in
#' part by those factors' ability to act as a proxy for age.
