Nothing
## ----setup, message = FALSE, warning = FALSE----------------------------------
library(rprev)
library(survival)
data(prevsim)
## -----------------------------------------------------------------------------
summary(prevsim)
## ----basicsurvival, fig.height=4, fig.width=7, echo=F-------------------------
survf <- survfit(Surv(time, status) ~ sex, data=prevsim)
survf_df <- data.frame(t=survf$time, s=survf$surv, sex=rep(c('M', 'F'), survf$strata),
stringsAsFactors = TRUE)
ggplot2::ggplot(survf_df, ggplot2::aes(x=t, y=s, colour=sex)) +
ggplot2::geom_line() +
ggplot2::theme_bw() +
ggplot2::labs(x='Time (days)', y='Survival probability') +
ggplot2::ylim(0, 1)
## -----------------------------------------------------------------------------
prevalence_total <- prevalence(index='2013-01-30',
num_years_to_estimate=c(3, 5, 10, 20),
data=prevsim,
inc_formula = entrydate ~ sex,
surv_formula = Surv(time, status) ~ age + sex,
dist='weibull',
population_size = 1e6,
death_column = 'eventdate')
## -----------------------------------------------------------------------------
prevalence_total
## -----------------------------------------------------------------------------
summary(prevalence_total)
## -----------------------------------------------------------------------------
prevalence_total$estimates
## -----------------------------------------------------------------------------
head(prevalence_total$simulated)
## -----------------------------------------------------------------------------
llog <- flexsurv::flexsurvreg(Surv(time, status) ~ age + sex, data=prevsim, dist='llogis')
llog
## ---- cache=T-----------------------------------------------------------------
prev_llog <- prevalence(index='2013-01-30',
num_years_to_estimate=c(3, 5, 10, 20),
data=prevsim,
inc_formula = entrydate ~ sex,
surv_model=llog,
population_size = 1e6,
death_column = 'eventdate',
N_boot = 100)
## ---- cache=T-----------------------------------------------------------------
prev_llog
## -----------------------------------------------------------------------------
fix_cure_mod <- fixed_cure(Surv(time, status) ~ age + sex, data=prevsim, cure_time=5*365.25,
dist='weibull')
## -----------------------------------------------------------------------------
prevalence(index='2013-01-30',
num_years_to_estimate=20,
data=prevsim,
inc_formula = entrydate ~ sex,
surv_model=fix_cure_mod, # Pass in the cure model that was built above
population_size = 1e6,
death_column = 'eventdate',
N_boot = 30)
## -----------------------------------------------------------------------------
prevalence_total
## -----------------------------------------------------------------------------
mixture_cure <- flexsurvcure::flexsurvcure(Surv(time, status) ~ age + sex, data=prevsim, dist='weibull', link='logistic', mixture=TRUE)
## -----------------------------------------------------------------------------
prevalence(index='2013-01-30',
num_years_to_estimate=20,
data=prevsim,
inc_formula = entrydate ~ sex,
surv_model=mixture_cure,
population_size = 1e6,
death_column = 'eventdate',
N_boot = 30)
## ---- eval=F------------------------------------------------------------------
# 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
# }
## ---- eval=F------------------------------------------------------------------
# data <- data.frame(...)
# myobj <- build_my_survival_object(Surv(time, status) ~ sex, data)
# prevalence(...
# surv_model=myobj, # This will work
# ...)
## ---- eval=F------------------------------------------------------------------
# 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
# ...)
## -----------------------------------------------------------------------------
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
}
## -----------------------------------------------------------------------------
data <- prevsim
pois_mod <- build_poisson(input_data=data)
## -----------------------------------------------------------------------------
pois_mod
## ---- eval=F------------------------------------------------------------------
# # 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, ...)
## -----------------------------------------------------------------------------
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
}
## -----------------------------------------------------------------------------
inc_times <- validate_incidence_model(pois_mod, prevsim)
## -----------------------------------------------------------------------------
head(inc_times)
## -----------------------------------------------------------------------------
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')
## -----------------------------------------------------------------------------
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
}
## -----------------------------------------------------------------------------
survobj <- build_wei(data=data)
survobj
## ---- eval=F------------------------------------------------------------------
# # 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)
## -----------------------------------------------------------------------------
extract_covars.mysurv <- function(object) {
"age"
}
## ---- eval=F------------------------------------------------------------------
# # 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)
## -----------------------------------------------------------------------------
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)
}
## -----------------------------------------------------------------------------
predict_survival_probability(survobj, newdata=data.frame(age=c(50, 70)), times=c(100, 100))
## -----------------------------------------------------------------------------
probs <- validate_survival_model(survobj, prevsim)
head(probs)
## -----------------------------------------------------------------------------
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)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.