########################################
## conditional probability of an event
## occuring before patient time 'min_f' and
## before calendar time 'r_period + f_period',
## given that a patient is recruited at 'r_time'.
## Note: 'r_time' is a vector here
cond_p_event = function(r_time,
min_f,
r_period,
f_period,
change_points,
lambdas,
mu,
dropouts){
f_time = pmin(min_f, r_period + f_period - r_time)
event_prob = numeric(length(f_time))
for (i in seq_along(event_prob)){
event_prob[i] = 1 - surv_pieces(f_time[i], change_points, lambdas, mu, dropouts)
}
event_prob
}
########################################
## probability density of recruitment time
p_r = function(r_time, r_period, k) k * (r_time / r_period) ^ (k - 1) / r_period
########################################
## cond_p_event times p_r
p_event_r = function(r_time,
min_f,
r_period,
f_period,
change_points,
lambdas,
mu,
dropouts,
k){
cond_p_event(r_time,
min_f,
r_period,
f_period,
change_points,
lambdas,
mu,
dropouts) * p_r(r_time, r_period, k)
}
########################################
## probability of an event occuring
## prior to patient-time c(change_points, Inf) AND
## before calendar time 'r_period + f_period'
p_event = function(r_period,
f_period,
change_points,
lambdas,
mu,
dropouts,
k){
min_fs = c(change_points, Inf)
integrals = numeric(length(min_fs))
for (i in seq_along(integrals)){
integrals[i] = integrate(p_event_r,
lower = 0,
upper = r_period,
min_f = min_fs[i],
r_period = r_period,
f_period = f_period,
change_points = change_points,
lambdas = lambdas,
mu = mu,
dropouts = dropouts,
k = k)$value
}
integrals
}
#' Expected number of events
#'
#' \code{expected_events_single_arm} returns the expected number of events at a given data-cut-off time.
#' @param{dco} Time of data cut-off.
#' @param{recruitment} List of recruitment information.
#' Containing \enumerate{
#' \item Sample size, \code{n}
#' \item Recruitment period, \code{r_period}
#' \item Recruitment parameter for power model, \code{k}
#' }
#' @param{model} The piecewise hazard model.
#' A list containing the \code{change_points} and \code{lambdas}.
#' @param{mu} The dropout rate. Assumed independent to event hazard.
#' @param{dropouts} Is dropout the event of interest? Default is FALSE.
#' @param{total_only} By default, the function returns only the expected number of events.
#' To see more detailed output set \code{total_only = FALSE}.
#' @return The expected number of events at \code{dco}.
#' @export
expected_events_single_arm = function(dco = 28,
recruitment,
model,
mu = 0,
dropouts = FALSE,
total_only = TRUE){
n = recruitment$n
r_period = recruitment$r_period
k = recruitment$k
change_points = model$change_points
lambdas = model$lambdas
if (length(dco) > 1) stop("dco must be length 1")
if (length(n) > 1) stop("n must be length 1")
# assume dco is greater than r_period.
# what if dco < r_period
if (dco < r_period){
n = n * (dco / r_period) ^ k
r_period = dco
}
f_period = dco - r_period
expected_events = n * p_event(r_period,
f_period,
change_points,
lambdas,
mu,
dropouts,
k)
total_events = expected_events[length(expected_events)]
if (total_only) return(total_events)
events_per_period = diff(c(0, expected_events))
prop_events = events_per_period / total_events
list(total_events = total_events,
events_per_period = cbind(c(change_points, Inf), events_per_period),
prop_events = prop_events,
n_recruited = n)
}
#' Expected number of events
#'
#' \code{expected_events_two_arm} returns the expected number of events at a given data-cut-off time.
#' @param{dco} Time of data cut-off.
#' @param{recruitment} List of recruitment information.
#' Containing \enumerate{
#' \item Sample size on control, \code{n_0}
#' \item Sample size on experimental, \code{n_1}
#' \item Recruitment period, \code{r_period}
#' \item Recruitment parameter for power model, \code{k}
#' }
#' @param{model} The piecewise hazard model.
#' A list containing the \code{change_points} and \code{lambdas}.
#' @param{mu} The dropout rate. Assumed independent to event hazard.
#' @param{dropouts} Is dropout the event of interest? Default is FALSE.
#' @param{total_only} By default, the function returns only the expected number of events.
#' To see more detailed output set \code{total_only = FALSE}.
#' @return The expected number of events at \code{dco}.
#' @export
expected_events_two_arm = function(dco = 28,
recruitment,
model,
mu = 0,
dropouts = FALSE,
total_only = TRUE){
n_0 = recruitment$n_0
n_1 = recruitment$n_1
r_period = recruitment$r_period
k = recruitment$k
change_points = model$change_points
lambdas_0 = model$lambdas_0
lambdas_1 = model$lambdas_1
if (length(dco) > 1) stop("dco must be length 1")
if (length(n_0) > 1) stop("n_0 must be length 1")
if (length(n_1) > 1) stop("n_1 must be length 1")
# assume dco is greater than r_period.
# what if dco < r_period
if (dco < r_period){
n_0 = n_0 * (dco / r_period) ^ k
n_1 = n_1 * (dco / r_period) ^ k
r_period = dco
}
f_period = dco - r_period
expected_events_0 = n_0 * p_event(r_period,
f_period,
change_points,
lambdas_0,
mu,
dropouts,
k)
expected_events_1 = n_1 * p_event(r_period,
f_period,
change_points,
lambdas_1,
mu,
dropouts,
k)
total_events_0 = expected_events_0[length(expected_events_0)]
total_events_1 = expected_events_1[length(expected_events_1)]
total_events = total_events_0 + total_events_1
if (total_only) return(c(total_events_0 = total_events_0,
total_events_1 = total_events_1,
total_events = total_events))
events_per_period = cbind(events_0 = diff(c(0, expected_events_0)),
events_1 = diff(c(0, expected_events_1)))
prop_events = rowSums(events_per_period) / total_events
hrs = lambdas_1 / lambdas_0
av_hr = exp(sum(log(hrs) * prop_events))
list(total_events_0 = total_events_0,
total_events_1 = total_events_1,
total_events = total_events,
events_per_period = cbind(c(change_points, Inf), events_per_period),
prop_events = prop_events,
hrs = hrs,
av_hr = av_hr,
n_recruited_0 = n_0,
n_recruited_1 = n_1,
n_recruited = n_0 + n_1)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.