During each bootstrap iteration, new incidence and survival models are fitted to the bootstrapped registry data, generating new parameters. Based on these parameters, each model is then used for prediction, either predicting an incident population or predicting the survival of this population.
This section provides details on how to provide custom models for both of these 2 roles.
Fitting models
Both incidence and survival objects must contain a call object that holds the initial function call used to build the model; this is obtained through match.call(). This call must contain an argument (name not important) which is passed the value data, as it is this argument which is changed to provide the bootstrapped data during simulation.
For example:
build_my_survival_object <- function(formula, input_data) {
# Build a survival model using the specified formula and input data
model <- ...
object <- list(model=model,
call=match.call()) # the function call must be included as an item 'call'
class(object) <- "myobj"
object
}
It is crucial that the parameter passing in the data to fit the model to is labelled data, as below.
data <- data.frame(...)
myobj <- build_my_survival_object(Surv(time, status) ~ sex, data)
prevalence(...
surv_model=myobj, # This will work
...)
The example below will not work.
some_data <- data.frame(...)
myobj <- build_my_survival_object(Surv(time, status) ~ sex, some_data)
prevalence(...
surv_model=myobj, # This WON'T work, since the data parameter was called 'some_data' instead
...)
Once the models have been built, generic methods are used to generate the estimated incident cases and survival probabilities, which is achieved by specifying an appropriate S3 class method. The following sections describes these methods and their parameterisation. See Hadley Wickham’s guide to S3 objects for further support on object-oriented programming in R.
An additional source of support is the source code for the existing objects that have been provided with the package which is freely available on CRAN and the development code is hosted on GitHub. For example, homogeneous_poisson.R contains the necessary methods for the default incidence model, and survregmin.R and flexsurv.R provide survival objects for the default survival implementation and flexsurv objects respectively.
Incidence
In this example a homogeneous Poisson process will still be used to demonstrate how to provide custom incidence objects. This process is parameterised by a single value: the rate \(\lambda\), which will need to be saved in the model object along with the function call.
The function below builds an object that contains both the process parameter (\(\lambda\)) and the function call.
build_poisson <- function(input_data) {
rate <- nrow(input_data) / as.numeric(difftime(max(input_data$entrydate), min(input_data$eventdate)))
# Build a survival model using the specified formula and input data
object <- list(rate=rate,
call=match.call()) # the function call must be included as an item 'call'
class(object) <- "pois"
object
}
When building the object, remember that the input data frame needs to be called data.
data <- prevsim
pois_mod <- build_poisson(input_data=data)
Printing the object shows that the requirements are met:
- Any required parameters are saved (rate)
- The call is saved and the input data was passed in as
data
- The object has a class (pois) in this case
## $rate
## [1] 0.2721829
##
## $call
## build_poisson(input_data = data)
##
## attr(,"class")
## [1] "pois"
Incidence objects need to implement a method called draw_incident_population that will generate the incident population specified by their incident times and any covariates that are required for the survival modelling. The incident times are relative to the origin, which in the simulation is the index minus N years, where N is max(num_years_to_estimate). The required parameterisation of draw_incident_population is shown below.
# object: The incidence model that will have been created on the bootstrapped data
# data: The data available to fit the model on, will be the registry data set provided to prevalence as this acts as the attribute prior distribution.
# timeframe: A single number specifying how long to generate incident cases for.
# covars: A character vector specifying the names of individual covariates that must be included in the returned data.table (or data frame)
# Returns a data.table (or data frame but data.table is preferred) where each row represents an incident case with:
# - The first column being the time since the origin, i.e. index date - N year prevalence
# - Any subsequent columns holding covariates that must be provided as specified in the 'covars' argument
draw_incident_population <- function(object, data, timeframe, covars, ...)
For this example using a homogeneous Poisson process, inter-arrival times are exponentially distributed, so simulating an incident population simply requires sampling inter-arrival times, turning these into arrival times, and then sampling individual attributes from the prior information (the registry data that is passed into the function).
draw_incident_population.pois <- function(object, data, timeframe, covars, ...) {
# Firstly draw inter-arrival times in the period [0, timeframe].
# The expected number is simply timeframe * rate so we'll take this amount + a margin for error.
expected <- 1.5 * (timeframe * object$rate)
# Now draw inter-arrival times
inter_arrival <- rexp(expected, object$rate)
# Determine absolute incident times
incident_times <- cumsum(inter_arrival)
# Truncate to those within the timeframe
incident_times <- incident_times[incident_times < timeframe]
num_incident <- length(incident_times)
# Sample individual attributes into a matrix. The required attributes are given by 'covars' argument
attrs <- do.call('cbind', lapply(covars, function(x) sample(data[[x]], num_incident, replace=T)))
# Now add the incident times as the first column
attrs <- cbind(incident_times, attrs)
# Convert to data frame and add column names
attrs <- data.frame(attrs)
colnames(attrs) <- c('incident_time', covars)
# Return this data frame
attrs
}
To validate that an incidence model has been correctly specified, the validate_incidence_model function has been provided. It accepts the incidence model itself and the registry data that it has been designed to work with. It will verify that the required attributes and methods are available, and that draw_incident_population successfully simulates individuals.
If any issues are found then the function stops execution and displays an error message detailing the fault, otherwise it returns a list of simulated arrival times, allowing for further diagnostics to be performed.
inc_times <- validate_incidence_model(pois_mod, prevsim)
## Incidence model passed all tests.
## incident_time
## 1 4.738095
## 2 5.492447
## 3 12.004696
## 4 14.582234
## 5 20.676506
## 6 26.975532
Once the object has been validated, it can be used in prevalence through the inc_model argument. Note that an additional argument is required to provide the name of the column in the data set that provides the incident dates, since this is no longer provided by the unused inc_formula option.
prevalence(index='2013-01-30',
num_years_to_estimate=c(3, 5, 10, 20),
data=prevsim,
inc_model = pois_mod,
surv_formula = Surv(time, status) ~ age + sex,
dist='weibull',
population_size = 1e6,
incident_column = 'entrydate',
death_column = 'eventdate')
## Warning in `[.data.table`(results, , `:=`(c("time_to_index", "time_to_entry"), :
## Column 'time_to_entry' does not exist to remove
## Estimated prevalence at 2013-01-30:
## 3 years: 208 (20.8 per 1e+05)
## 5 years: 304 (30.4 per 1e+05)
## 10 years: 514 (51.4 per 1e+05)
## 20 years: 839.14 (83.91 per 1e+05)
Survival
For this example a Weibull model with age as the sole covariate will be used. The survival model will be fitted using the flexsurv package, this functionality is already implemented as discussed earlier, but it is an appropriate demonstration.
A Weibull survival model is parameterised by its coefficients for each covariate and the distribution specific parameters. The function below builds an object of class mysurv that contains these coefficients, as well as the function call (which is also saved by the flexsurv object).
library(flexsurv)
build_wei <- function(data) {
mod <- flexsurvreg(Surv(time, status) ~ age, data=data, dist='weibull')
obj <- list(coefs=coef(mod),
call=match.call())
class(obj) <- 'mysurv'
obj
}
With just these two attributes, a fully specified survival model has been generated. It has the required saved information:
- Any required parameters are saved
- The call is saved and the input data was passed in as
data
- The object has a class (mysurv)
survobj <- build_wei(data=data)
survobj
## $coefs
## shape scale age
## -0.64951710 10.86692915 -0.04472098
##
## $call
## build_wei(data = data)
##
## attr(,"class")
## [1] "mysurv"
Survival objects have two methods that need to be implemented:
extract_covarspredict_survival_probability
extract_covars simply returns a character string detailing which of the covariates passed into prevalence through the data argument are used in the survival model. This allows the simulation to know how to generate the incident population as described above in draw_incident_population. In fact, the output of extract_covars is passed directly into draw_incident_population through the covars parameter.
# object: The survival model
# Returns a character vector detailing the covariates required to fit this model. All of
# these values should be columns in the data that is passed in the main 'prevalence' function.
extract_covars <- function(object)
For this survival model age is the only patient-level covariate being used.
extract_covars.mysurv <- function(object) {
"age"
}
predict_survival_probability estimates survival probability at the index date. It is specified as follows:
# object: The survival object
# newdata: A data frame (or data.table) with the incident population stored with their
# required covariates for the model.
# times: A vector of times to estimate survival probability at for individuals in
# corresponding rows of 'newdata'. This should be the same length as there are
# rows in 'newdata' since each individual has their survival probability estimated once.
# Returns:
# A vector of survival probabilities of length equal to the length of 'times' and the
# number of rows in 'newdata'.
predict_survival_probability <- function(object, newdata, times)
For this simple Weibull model these estimates are simply provided by \(1-F(x)\), with the CDF already implemented in base R as pweibull.
predict_survival_probability.mysurv <- function(object, newdata, times) {
# Calculate linear predictor, this will form the shape parameter
shape <- exp(object$coefs[1] + newdata$age*object$coefs[3])
scale <- exp(object$coefs[2])
1 - pweibull(times, shape, scale)
}
While more in-depth testing would be required to validate the predictions output by predict_survival_probability, from a programming perspective at least it is outputting numbers.
predict_survival_probability(survobj, newdata=data.frame(age=c(50, 70)), times=c(100, 100))
## [1] 0.4941064 0.4202817
There is also a function validate_survival_model that checks the survival model contains the required attributes and methods, and that predict_survival_probability outputs sensible probabilities that are monotonically decreasing with time. If the test passes, it returns survival probabilities taken from random individuals in the supplied registry data at random time-points.
probs <- validate_survival_model(survobj, prevsim)
## Survival model passed all tests.
## [1] 0.4211790 0.4148378 0.3953903 0.3866155 0.4464323 0.3957396
Plugging this model into the prevalence function now works.
prevalence(index='2013-01-30',
num_years_to_estimate=c(3, 5, 10, 20),
data=prevsim,
inc_formula = entrydate ~ 1,
surv_model = survobj,
population_size = 1e6,
death_column = 'eventdate',
N_boot = 100)
## Estimated prevalence at 2013-01-30:
## 3 years: 208 (20.8 per 1e+05)
## 5 years: 304 (30.4 per 1e+05)
## 10 years: 514 (51.4 per 1e+05)
## 20 years: 898.82 (89.88 per 1e+05)
