Nothing
fit_exponential_incidence <- function(inc_form, data) {
entry_col <- all.vars(update(inc_form, .~0))
if (is.null(entry_col) || length(entry_col) != 1) {
stop("Error: Please provide only a single LHS term in the inc_formula, representing the column holding incident date.")
}
if (! entry_col %in% colnames(data)) {
stop("Error: Column '", entry_col, "' not found in '", quote(data), "'.")
}
strata <- all.vars(update(inc_form, 0~.))
# Obtain rate for combination of each of these factors and throw error if any instances < required
if (length(strata) == 1) {
if (! strata %in% colnames(data)) {
stop("Error: Column '", strata, "' not found in '", quote(data), "'.")
}
if (!all(sapply(strata, function(x) is.factor(data[[x]])))) {
stop("Error: Please format any incidence strata as factors.")
}
strata_counts <- table(data[, strata])
if (any(strata_counts < MIN_INCIDENCE)) {
stop("Error: Less than ", MIN_INCIDENCE, " incident cases. Please ensure data has sufficient incident cases.")
}
# Calculate time between first and last patient for each strata
durations <- vapply(names(strata_counts), function(s) {
sub <- data[[strata]] == s
as.numeric(difftime(max(data[sub, entry_col]), min(data[sub, entry_col]), units='days'))
}, numeric(1))
rate <- data.frame(strata_counts / durations)
names(rate)[1:(ncol(rate)-1)] <- strata
strata_names <- strata
} else if (length(strata) > 1) {
stop("Default incidence model can currently only handle one stratification.")
} else {
rate <- nrow(data) / as.numeric(difftime(max(data[[entry_col]]), min(data[[entry_col]]), units='days'))
strata_names <- NULL
}
# Formulate call with the formula evaluated in this environment with the data argument as a promise
# that will be evaluated later on the bootstrapped data
func_call <- match.call()
func_call$inc_form <- eval(inc_form)
obj <- list(rate=rate, call=func_call, strata=strata_names)
class(obj) <- "expinc"
obj
}
#' @export
draw_incident_population.expinc <- function(object, data, timeframe, covars) {
if (timeframe < 0) {
stop("Error: timeframe is negative with a value of ", timeframe, ".")
}
for (c in covars) {
if (!c %in% colnames(data)) {
stop("Error: Column '", c, "' not found in '", quote(data), "'.")
}
}
if (is.null(object$strata)) {
initial_num_inds <- INCIDENCE_MARGIN * object$rate * timeframe
time_to_entry <- cumsum(rexp(initial_num_inds, object$rate))
time_to_entry <- time_to_entry[time_to_entry < timeframe]
new_data <- data.table::data.table(time_to_entry)
} else {
foo <- apply(object$rate, 1, function(row) {
rate <- as.numeric(row['Freq'])
initial_num_inds <- INCIDENCE_MARGIN * rate * timeframe
time_to_entry <- cumsum(rexp(initial_num_inds, rate))
time_to_entry <- time_to_entry[time_to_entry < timeframe]
# Form data frame with these factor levels too
new_data <- data.table::data.table(time_to_entry)
new_data[, names(row[-length(row)]) := row[-length(row)]]
new_data
})
new_data <- data.table::rbindlist(foo)
}
if (length(covars) > 0) {
# draw covariate values from model that are in survival formula and not in data frame
missing_covars <- setdiff(covars, colnames(new_data))
higher_order <- missing_covars[grep("^I\\(.+\\)$", missing_covars)]
covars_to_sample <- setdiff(missing_covars, higher_order)
# Possible to have length covars > 0 and no missing variables, i.e. if the variable is used as
# stratification in incidence model
if (length(covars_to_sample) > 0) {
# Draw new covariates from prior with appropriate stratification
if (is.null(object$strata)) {
new_covars <- setNames(lapply(covars_to_sample, function(x) sample(data[[x]], nrow(new_data), replace=T)), covars_to_sample)
data.table::setDT(new_covars)
} else {
strata_vals <- object$rate[[object$strata]]
freqs <- table(new_data[[object$strata]])
covars_list <- setNames(lapply(strata_vals, function(sx) {
priors <- data[data[[object$strata]] == sx, ]
num <- freqs[as.character(sx)]
setNames(lapply(covars_to_sample, function(x) sample(priors[[x]], num, replace=T)), covars_to_sample)
}), strata_vals)
new_covars <- data.table::rbindlist(covars_list)
}
# Add in higher order terms
# NB This method of evaluating higher order won't work if the covariate isn't in this new_covars data, however,
# this is acceptable as the covariates missing from this data set will be those used in the incidence modelling and
# so will be categorical and thus not not have higher order terms
for (term in higher_order) {
new_covars[, (term) := eval(parse(text=term))]
}
# Combine together to form a population that has an entry time and the required covariates for modelling survival
new_data <- cbind(new_data, new_covars)
}
}
new_data
}
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.