Time to Event Dose Finding Trial
Berry Consultants, LLC
1 Introduction
We consider here a dose-finding trial with a time to event primary endpoint. The primary analysis of the primary endpoint is a Bayesian piecewise exponential model (Ibrahim, Chen, and Sinha 2001; Christensen et al. 2010) of the time to event, with an underlying dose-response model to borrow information across the treatment arms. For illustration purposes, we consider a trial with a control arm and four active doses.
2 Piecewise Constant Proportional Hazard Model
Let \(T_{i}\) be the event time (for modeling purposes, measured in weeks) for the primary endpoint for the \(i^{th}\) patient. Without loss of generality, we can take the initial time after which the event can occur to be \(t = 0\) as it is usually done in survival analysis. We also make the more restrictive assumption that the event will always occur so that event times are always finite. In other words, we often assume that event times are distributed across the open time interval \(t \in (0, +\infty)\). When follow-up for patients is finite, as in this case study, we will assume that \(t \in (0, T_{max})\)
Let us now consider a set of \(J\) pre-determined cutpoints that partition the range of the data, i.e. \(0 < s_{1} < \dots < s_{J} = T_{max}\). These divide the time domain in \(J\) intervals \(\underbrace{[0, s_{1}]}_{I_{1}}, \underbrace{(s_{1}, s_{2}]}_{I_{2}}, \dots, \underbrace{(s_{J-1}, s_{J}]}_{I_{J}}\).
In our motivating example, we consider \(J = 4\) cutpoints \(\mathbf{s} = (5, 10, 15, 20)\).
2.1 Model for Control Hazard Rates
In many cases, in survival analysis the occurrence of an event is promoted by the accumulation of some “stimulus.” We refer to this stimulus as a hazard. A hazard function, \(\lambda(t)\), quantifies the differential amount of hazard introduced to the system at time \(t\).
A relatively simple model considers as baseline hazard (hazard function for the control arm) a constant function on each interval, i.e. \[ \lambda(t) = \lambda_{j} I_{(s_{j-1}, s_{j}]}(t), \] where \(I_{(s_{j-1}, s_{j}]}(t)\) is the indicator function that is equal to \(1\) if \(t \in (s_{j-1}, s_{j}]\), and \(0\) otherwise.
The cumulative hazard function \[ H(t) = \int_{0}^{t} \lambda(s) \ \mathrm{d} s = \sum_{j = 1}^{J} \lambda_{j} \int_{0}^{t} I_{j}(s) \ \mathrm{d} s \] quantifies the total hazard accumulated up to time \(t\).
In the survival modeling context, \(S(t)\) quantifies the probability that the event occurs anytime after time \(t\).
This can also be interpreted as the probability that the event does not occur before time \(t\), or in other words the probability that it survives until time \(t\).
Because of this \(S(t)\) is referred to as the survival function.
The more hazard that accumulates by a given time the less likely the event should have survived to that time.
In other words, the survival function should decay with increasing hazard.
Survival models make the particular assumption of an exponential decay,
\[
S(t) = \exp \left\{ - \sum_{j = 1}^{J} \lambda_{j} \int_{0}^{t} I_{j}(s) \ \mathrm{d} s \right\}.
\]
This expression can be made more explicit, for example, considering an event occurring at \(t \in I_{k}\), i.e., \(s_{k-1} < t < s_{k}\).
Then we have
\[\begin{equation}
S(t) = \exp \left\{ - \sum_{j = 1}^{k-1} \lambda_{j} (s_{j} - s_{j-1}) - \lambda_{k} (t - s_{k-1}) \right\}.
\tag{2.1}
\end{equation}\]
In other words, assuming a piecewise constant hazard function is equivalent to assuming a piecewise exponential model for the event time.
We will discuss the prior distribution on the parameters \(\lambda_{1}, \dots, \lambda_{J}\) in Section 4.1.
Continuing our practical example, consider the weekly hazard rates \(\mathbf{\lambda} = (0.06, 0.06, 0.03, 0.03)\). With this value of \(\mathbf{\lambda}\), we can visualize the hazard function, the cumulative hazard function, and the survival function.
2.2 Model for Treatment Hazard Rates
In the presence of multiple arms, we still model the survival distribution for the event times as piecewise exponential: \[T_{i} \sim PE(\boldsymbol \eta_{i}),\] where \(\boldsymbol \eta_{i} = (\eta_{i1}, \dots, \eta_{iJ})\) represents the set of hazard rates (events per week) within each segment for the \(i^{th}\) patient.
The hazard rates are modeled as follows: \[\eta_{ij} = \lambda_{j} \exp\{ X_{i,d} \theta_{d} \}, \quad j = 1, \dots, J, \ i = 1, \dots, n, \ d = 1, \dots, D,\] where \(\lambda_{j}\) is the baseline hazard rate for time segment \(j\), \(\theta_{d}\) is the log hazard ratio (HR) for dose \(d\), and \(X_{i,d}\) is the treatment indicator variable, where \(X_{i,d} = 1\) if patient \(i\) is randomized to dose \(d\), and zero otherwise. We will discuss the prior distribution on the parameters \(\theta_{1}, \dots, \theta_{D}\) in Section 4.1.
In our example, let us consider \(D = 4\) doses with hazard ratios \(e^{\mathbf{\theta}} = (1.5, 2, 2.5, 1.75)\).
2.3 Dose Response Model
A dose-response model is used to borrow information for the treatment effects across each of the four active doses, leading to increased estimation efficiency. The log hazard ratio (HR) for each dose, \(\theta_{d}\), is modeled using a normal dynamic linear model (NDLM), i.e. \[ \theta_{1} \sim \text{N}(0, 1), \quad \theta_{d} \sim \text{N}(\theta_{d-1}, \tau^{2}), \ d = 2, \dots, D, \] where dose 1 has an independent standard normal prior distribution, and each subsequent dose has a prior distribution centered on the previous dose, with a variance term \(\tau^{2}\).
The parameter \(\tau^{2}\) controls the degree of smoothing between doses. To visualize its effect, we can repeatedly sample from the NDLM using some fixed values of \(\tau^{2}\).
Small values of \(\tau^{2}\) indicate that the responses of successive doses are likely to be very similar, and that there is more borrowing (smoothing) across doses. Large values of \(\tau^{2}\) indicate that the responses of successive doses are highly variable, and that there is less smoothing across doses.
Note: we sampled dose-response curves for fixed values of \(\tau^{2}\). In the model, we will put a prior on \(\tau^{2}\), therefore the implied prior on the dose-response curve will have an even larger variability (it incorporates uncertainty on \(\tau^{2}\) as well). We will discuss the prior distribution on the parameter \(\tau^{2}\) in Section 4.1.
3 Data Generation
To demonstrate how we can make inference on this model, we first need to simulate some data.
Let us load some useful functions and create the quantities corresponding to the underlying true model parameter.
Click to explore the functions in utils_tte.R.
installer <- function(package){
installed_packages <- package %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
install.packages(packages[!installed_packages])
}
}
# List of dependent packages
packages = c(
'tibble',
'dplyr',
'tidyr',
'readr',
'stringr',
'knitr',
'ggplot2',
'latex2exp',
'RColorBrewer',
'rstan',
'bayesplot',
'coda',
'survival',
'survminer',
'reshape2',
'parallel',
'kableExtra'
)
# Install + source the packages
installer(packages)
invisible(lapply(packages, library, character.only = TRUE))
# Initialize global variables
if ('rstan' %in% packages) {
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores()) # use multiple cores if available, speeds up simulations
}
theme_set(theme_bw(base_size = 14))
options(dplyr.summarise.inform = FALSE)
# Auxiliary ---------------------------------------------------------------
#' Extract the legend from a ggplot
#' @param plt: plot whose legend is to be extracted
#' @return legend object
ext_legend = function(plt){
tmp = ggplot_gtable(ggplot_build(plt))
leg = which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
legend = tmp$grobs[[leg]]
return(legend)
}
#' Empty theme
custom_theme = function() {
theme_survminer() %+replace%
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
panel.grid.major = element_line(colour = "grey90"),
panel.border = element_blank(),
panel.background = element_blank()
)
}
#' Rename the columns of a dataframe
#' @param df: data.frame whose column are to be renamed
#' @param names: names of the columns
#' @return data.frame with new column names
colnames_inplace = function(df, names){
colnames(df) = names
df
}
#' Function to calculate the integrated accrual profile
#' @param t: time
#' @param peak_rate: peak accrual rate (weeks)
#' @param ramp_up: ramp up complete at week
#' @param t0: date of last enrolled patient (in weeks, since start of trial)
#' @return data.frame with new column names
Lambda = function(t, peak_rate, ramp_up, t0 = 0) {
if (ramp_up == 0){
peak_rate * (t + t0)
}
else {
0.5 * peak_rate/ramp_up * (t + t0)^2 * as.numeric((t + t0) < ramp_up) +
(peak_rate * ((t + t0) - ramp_up) + 0.5 * peak_rate * ramp_up) * as.numeric((t + t0) >= ramp_up)
}
}
#' Calculate the hazard function for piecewise constant hazards
#' @param t: time grid where to evaluate the cumulative hazard function
#' @param s_j: cut-points of the J intervals where the hazards are constant (length J + 1)
#' @param lambda: hazard rates in each of the J intervals (if a matrix, the hazard rates are on each row)
#' @return hazard function for each possible row of lambda
ht = function(t, s_j, lambda){
J = length(s_j) - 1
if (any(t < 0)) {
stop("Cannot have negative event times!\n")
}
s_j[J+1] = +Inf
if (is.vector(lambda)){
J = length(lambda)
out = array(NA, dim = length(t))
for (j in 1:J){
idx = which( (t >= s_j[j]) & (t < s_j[j+1]))
out[idx] = lambda[j]
}
}
else if (is.matrix(lambda)){
J = ncol(lambda)
out = array(NA, dim = c(nrow(lambda), length(t)))
for (j in 1:J){
idx = which( (t >= s_j[j]) & (t < s_j[j+1]))
out[,idx] = lambda[,j]
}
}
else {
stop("The hazard rates lambda should be a vector or a matrix!\n")
}
return (out)
}
#' Calculate the cumulative hazard function for piecewise constant hazards
#' @param t: time grid where to evaluate the cumulative hazard function
#' @param s_j: cut-points of the J intervals where the hazards are constant (length J + 1)
#' @param lambda: hazard rates in each of the J intervals (if a matrix, the hazard rates are on each row)
#' @return cumulative hazard function for each possible row of lambda
Ht = function(t, s_j, lambda){
J = length(s_j) - 1
s_j[J+1] = +Inf
if (is.vector(lambda)){
J = length(lambda)
out = array(0, dim = length(t))
}
else if (is.matrix(lambda)){
J = ncol(lambda)
out = array(0, dim = c(length(t), nrow(lambda)))
}
else {
stop("The hazard rates lambda should be a vector or a matrix!\n")
}
EXP = t(outer(t, s_j[2:(J+1)], FUN = function(x, y){pmin(x, y)})) - s_j[1:J]
EXP[which(EXP < 0)] = 0
if (is.vector(lambda)){
out = as.numeric(colSums(EXP * lambda))
}
else if (is.matrix(lambda)){
out = t(lambda %*% EXP)
}
return (out)
}
#' Calculate the survival function for piecewise constant hazards
#' @param t: time grid where to evaluate the survival function
#' @param s_j: cut-points of the J intervals where the hazards are constant (length J + 1)
#' @param lambda: hazard rates in each of the J intervals (if a matrix, the hazard rates are on each row)
#' @return survival function for each possible row of lambda
St = function(t, s_j, lambda){
out = exp(-Ht(t, s_j, lambda))
return (out)
}
# Data generation ---------------------------------------------------------
#' Normalized integrated accrual profile
#' @param x: time
#' @param T_max: maximum follow up time
#' @param peak_rate: peak accrual rate (weeks)
#' @param ramp_up: ramp up complete at week
#' @param t0: date of last enrolled patient (in weeks, since start of trial)
#' @return data.frame with new column names
Ft = function(x, T_max, peak_rate, ramp_up, t0 = 0) {
(Lambda(x, peak_rate, ramp_up, t0) - Lambda(0, peak_rate, ramp_up, t0)) / (Lambda(T_max, peak_rate, ramp_up, t0) - Lambda(0, peak_rate, ramp_up, t0))
}
#' Inverse of the normalized integrated accrual profile
#' @param Ft: Normalized integrated accrual profile
#' @param u: uniform draw (to use the inverse cdf method)
#' @param T_max: maximum follow up time
#' @param peak_rate: peak accrual rate (weeks)
#' @param ramp_up: ramp up complete at week
#' @param t0: date of last enrolled patient (in weeks, since start of trial)
#' @return data.frame with new column names
Ftinv = function(Ft, u, T_max, peak_rate, ramp_up, t0 = 0) {
N = length(u)
a = rep(0, N)
b = rep(T_max, N)
binf = rep(0, N)
bsup = rep(0, N)
for (j in 1:20) {
idx_low = which(Ft((a + b) / 2, T_max, peak_rate, ramp_up, t0) <= u)
binf[idx_low] = (a[idx_low] + b[idx_low])/2
bsup[idx_low] = b[idx_low]
idx_upp = which(Ft((a + b) / 2, T_max, peak_rate, ramp_up, t0) >= u)
bsup[idx_upp] = (a[idx_upp] + b[idx_upp]) / 2
binf[idx_upp] = a[idx_upp]
a = binf
b = bsup
}
return ( (a + b) / 2 )
}
#' Sample observations from a piecewise constant hazard model
#' @param n: sample size
#' @param X: design matrix (without the column of leading ones)
#' @param lambda: hazard rates in each of the J intervals
#' @param theta: log HR for the different active doses (length D)
#' @param s_j: cut-points of the J intervals where the hazards are constant (length J + 1)
#' @param T_max: maximum follow up time after the last patient randomized
#' @param accr_profile: accrual profile
#' @param t0: date of last enrolled patient (in weeks, since start of trial)
#' @return (time, event, dates): event times (or censoring times, whichever is sooner), event indicators, and accrual dates
sample_pwc_data = function(n, X, lambda, theta, s_j, T_max, accr_profile, t0){
# Sample times from start of accrual according to a Poisson Process
n_max = rpois(1, Lambda(2*n/accr_profile$peak_rate, accr_profile$peak_rate, accr_profile$ramp_up, t0))
if (n_max < n){
stop('Cannot simulate n individuals. Consider increasing accrual!\n')
}
# Use the inverse cdf method to sample n_max accrual dates
dates_temp = Ftinv(Ft, runif(n_max), 2*n/accr_profile$peak_rate, accr_profile$peak_rate, accr_profile$ramp_up, t0)
# Take only n accrual dates, sort them and shift them by the initial date
dates = t0 + sort(dates_temp)[1:n]
# Now we use the inverse cdf method to sample from the likelihood model
u = runif(n)
XtBeta = as.numeric(X %*% theta)
time = event = dropout = array(NA, dim = n)
J = length(lambda)
s_j[J+1] = max(s_j[J+1], dates[length(dates)] - dates[1] + T_max + 1)
for (i in 1:n){
# Define the max censoring time as the gap time between ith enrollment and
# last enrollment plus the maximum follow up after last randomized (T_max)
obs_window_i = dates[length(dates)] - dates[i] + T_max
# Calculate the hazard rate for individual i
lambda_i = lambda * exp(XtBeta[i])
# Check which time interval the observation should fall into
H_temp = c(0, cumsum(lambda_i * diff(s_j)))
interval_idx = head(which(log(u[i]) >= - H_temp), 1) - 1
if (length(interval_idx) == 0){ # censored at the observation window endpoint
time[i] = obs_window_i
event[i] = 0
}
else { # sampled from the model in that interval
event[i] = 1
time[i] = s_j[interval_idx] - 1/lambda_i[interval_idx] * (log(u[i]) + H_temp[interval_idx])
if (time[i] > obs_window_i){ # it still could be censored in that interval
time[i] = obs_window_i
event[i] = 0
}
}
}
return (list('time' = time, 'event' = event, 'dates' = dates))
}
# Predictive probabilities ------------------------------------------------
#' Compute the summary statistics for running the piecewise constant hazards model
#' The methodology is described here: https://myweb.uiowa.edu/pbreheny/7210/f19/notes/10-24.pdf
#' @param y: event times or censoring times, whichever is sooner (vector of length N)
#' @param nu: event indicators (vector of length N)
#' @param s_j: cut-points of the J intervals where the hazards are constant (length J + 1)
#' @return (delta_i, Hij): interval number in which subject i failed, cumulative hazards
get_summary_stats = function(y, nu, s_j){
N = length(y)
J = length(s_j) - 1
s_j[J+1] = max(y) + 1
# Calculate the sufficient statistics for the STAN model to run
delta_i = array(0, dim = N) # interval number in which subject i failed
Hij = array(NA, dim = c(N, J))
for (i in 1:N){
if (nu[i] == 1){
delta_i[i] = head(which(y[i] < s_j), 1) - 1
}
Hij[i,] = pmin(y[i], s_j[2:(J+1)]) - s_j[1:J]
Hij[i,which(Hij[i,] < 0)] = 0
}
# Check if we have enough events in each interval
if (!(all(table(factor(delta_i, levels = 0:J)) > 4))) {
warning ("Some intervals have less than 5 observations. Consider collapsing into larger intervals.\n")
}
stopifnot(rowSums(Hij) == y) # Just an additional check
return (list('delta_i' = delta_i, 'Hij' = Hij))
}
#' Followup the current patients for T_max weeks using a piecewise constant
#' hazard model
#' @param subj_id: subject IDs for the currently enrolled patients
#' @param time: vector indicating exposure time for each patient
#' @param event: vector indicating if each patient has already had the event or not
#' @param dropout: vector indicating if each patient has already dropped out or not
#' @param lambda: hazard rates in each of the J intervals
#' @param s_j: cut-points of the J intervals where the hazards are constant (length J + 1)
#' @param X: design matrix (without the column of leading ones)
#' @param beta: regression coefficients (length p)
#' @param T_max: maximum follow up time (censoring window)
#' @param cens_rate: weekly rate corresponding to the desired censoring rate
#' @return imputed dataset
sample_followup_data_curr = function(subj_id, time, event, dropout, lambda, s_j,
X, beta, T_max, cens_rate){
# We use the inverse cdf method to impute using the piecewise constant hazard
# likelihood
idx_not_to_impute = which(event == 1 | dropout == 1) # if already observed, nothing to impute
idx_to_impute = which(event == 0 & dropout == 0) # if not observed, events could occur
# Now we use the inverse cdf method to sample from the likelihood model
n = length(time)
J = length(lambda)
s_j[J+1] = +Inf
XtBeta = as.numeric(X %*% beta)
imp_time = imp_event = array(NA, dim = n)
imp_dropout = array(0, dim = n)
imp_time[idx_not_to_impute] = time[idx_not_to_impute]
imp_event[idx_not_to_impute] = event[idx_not_to_impute]
imp_dropout[idx_not_to_impute] = dropout[idx_not_to_impute]
u = runif(n)
if (length(idx_to_impute) > 0){ # continue followup for T_max weeks
for (i in idx_to_impute){
lambda_i = lambda * exp(XtBeta[i])
H_temp = c(0, cumsum(lambda_i * diff(pmax(0, s_j - time[i]))))
interval_idx = head(which(log(u[i]) >= - H_temp), 1) - 1
if (length(interval_idx) == 0){ # event is censored
imp_time[i] = time[i] + T_max
imp_event[i] = 0
}
else { # event is sampled in one of the intervals defined by s_j
imp_time[i] = time[i] - 1/lambda_i[interval_idx] * (log(u[i]) + H_temp[interval_idx])
imp_event[i] = 1
if (imp_time[i] > time[i] + T_max){ # it still could be censored in that interval
imp_time[i] = time[i] + T_max
imp_event[i] = 0
}
}
# Additionally, a censoring time is sampled from an exponential distribution
# with a rate that results in the assumed dropout rate
c_time = time[i] + rexp(1, cens_rate) # using the memoryless property of Exp
# If this censoring time is less than the event time for the subject, then
# the patient will be censored at this time and will not count as an event.
if (c_time < imp_time[i]) {
imp_dropout[i] = 1
imp_event[i] = 0
imp_time[i] = c_time
}
}
}
imp_dat = tibble(ID = subj_id,
dropout = imp_dropout,
event = imp_event,
time = imp_time)
return (imp_dat)
}
#' Sample new observations from a piecewise constant hazard model
#' @param n: sample size to be imputed
#' @param lambda: hazard rates in each of the J intervals
#' @param s_j: cut-points of the J intervals where the hazards are constant (length J + 1)
#' @param T_max: maximum follow up time (censoring window)
#' @param cens_rate: weekly rate corresponding to the desired censoring rate
#' @return imputed dataset
sample_pwc_data_new = function(n, lambda, s_j, T_max, cens_rate){
# Now we use the inverse cdf method to sample from the likelihood model
u = runif(n)
time = event = array(NA, dim = n)
dropout = array(0, dim = n)
J = length(lambda)
for (i in 1:n){
# Check which time interval the observation should fall into
H_temp = c(0, cumsum(lambda * diff(s_j)))
interval_idx = head(which(log(u[i]) >= - H_temp), 1) - 1
if (length(interval_idx) == 0){ # censored at the observation window endpoint
time[i] = T_max
event[i] = 0
}
else { # sampled from the model in that interval
event[i] = 1
time[i] = s_j[interval_idx] - 1/lambda[interval_idx] * (log(u[i]) + H_temp[interval_idx])
if (time[i] > T_max){ # it still could be censored in that interval
time[i] = T_max
event[i] = 0
}
}
# Additionally, a censoring time is sampled from an exponential distribution
# with a rate that results in the assumed dropout rate at 1 year.
c_time = rexp(1, cens_rate) # censoring time SINCE TIME 0 (not post-blanking)
# If this censoring time is less than the event time for the subject, then
# the patient will be censored at this time and will not count as an event.
if (c_time < time[i]) {
dropout[i] = 1
event[i] = 0
time[i] = c_time
}
}
imp_dat = tibble(ID = paste0('imp_', 1:n),
DPOUTFL = dropout,
peefl = event,
EXPTIME = time)
return (imp_dat)
}source('utils_tte.R')
key_arms = c('Control', 'Dose 1', 'Dose 2', 'Dose 3', 'Dose 4')
key_palette = c('black', brewer.pal(9, 'Set1')[1:(length(key_arms) - 1)])
data_seed = 12345
s_j_true = c(0, 5, 10, 15, 20) # cutpoints
lambda_true = c(0.06, 0.06, 0.03, 0.03) # baseline hazard rates
theta_true = log(c(1.5, 2, 2.5, 1.75)) # hazard ratiosWe first have to generate the allocation to each arm. In our example, we use an allocation \(2:1:1:1:1\) in blocks of \(12\).
N = 250 # sample size
D = length(theta_true) # number of doses
rnd_ratio = c(4, 2, 2, 2, 2) # randomization ratio
blk_size = sum(rnd_ratio) # block size
# Generate the random allocations
set.seed(data_seed)
dose_label = tibble(BLOCK = c(rep(1:floor(N/blk_size), each = blk_size),
rep(floor(N/blk_size) + 1, N - length(rep(1:floor(N/blk_size),
each = blk_size))))) %>%
group_by(BLOCK) %>%
mutate(TRTNUM = sample(rep(0:D, times = rnd_ratio), length(BLOCK))) %>%
ungroup() %>%
pull(TRTNUM)We then create the allocation matrix \(X\) where \(X_{i,d} = 1\) if patient \(i\) is randomized to the active dose \(d\), and zero otherwise. In the example below, we show the first rows of \(X\), where the allocations correspond to (Control, Dose 3, Dose 2, Dose 4, Control, Dose 1).
# Put the allocations into the dose matrix format
X = matrix(0, N, D)
for (s in 1:D){
X[which(dose_label == s),s] = 1
}
# Show X in output
X %>%
colnames_inplace(sprintf('Dose %s', 1:D)) %>%
head(n = 6) %>%
knitr::kable(booktabs = TRUE,
caption = 'First rows of the dose allocation matrix X.') %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"))| Dose 1 | Dose 2 | Dose 3 | Dose 4 |
|---|---|---|---|
| 0 | 0 | 0 | 0 |
| 0 | 0 | 1 | 0 |
| 0 | 1 | 0 | 0 |
| 0 | 0 | 0 | 1 |
| 0 | 0 | 0 | 0 |
| 1 | 0 | 0 | 0 |
We can also setup the accrual parameters. For example, consider a peak rate accrual of \(3\) patients per week with no ramp-up period.
accr_profile = NULL
accr_profile$peak_rate = 3 # in weeks
accr_profile$ramp_up = 0 # in weeksGiven a survival function, exact simulation of event times can be achieved using a variant of the inverse cumulative distribution function method. In fact, solving \(S(t)\) for \(t\), we get \[\begin{equation} t = s_{k-1} - \dfrac{\log S(t)}{\lambda_{k}} - \dfrac{\sum_{j=1}^{k-1} \lambda_{j} (s_{j} - s_{j-1})}{\lambda_{k}} \tag{3.1} \end{equation}\]
Hence we can use the following method:
- Sample \(u \sim U(0, 1)\)
- Find the correct interval using the conditions \(s_{k-1} < t \leq s_{k}\), i.e. \[ \begin{cases} s_{k-1} < t & \Leftrightarrow && \log u < - \sum_{j = 1}^{k - 1} \lambda_{j} (s_{j} - s_{j-1}) \\ s_{k} \geq t & \Leftrightarrow && \log u \geq - \sum_{j = 1}^{k} \lambda_{j} (s_{j} - s_{j-1}) \end{cases} \]
- Calculate \(t\) using (3.1)
The figures below illustrate this algorithm, implemented in the function sample_pwc_data().
Let us generate some data, and then we can censor them the maximum follow-up time of \(20\) weeks.
set.seed(data_seed)
data = sample_pwc_data(n = N, # sample size
X = X, # dose allocation matrix
lambda = lambda_true, # baseline hazard rates
theta = theta_true, # log hazard ratios
s_j = s_j_true, # cutpoints
T_max = 0, # follow-up after last randomized patient
accr_profile = accr_profile, # accrual profile
t0 = 0 # current week (0 for beginning of time)
)
# Censor the observations at the maximum follow-up
data$event[which(data$time > 20)] = 0
data$time[which(data$time > 20)] = 20
y = as.numeric(data$time) # either censoring time (if event == 0) or death time (if event == 1)
nu = as.numeric(data$event) # event indicatorWe can visualize the data we just generated with a Kaplan-Meier plot, and compare it with the (theoretical) survival curves.
# True survival curves
tibble(X = rep(t_grid, length(key_arms)),
y = c(St(t_grid, s_j_true, lambda_true),
St(t_grid, s_j_true, exp(theta_true[1]) * lambda_true),
St(t_grid, s_j_true, exp(theta_true[2]) * lambda_true),
St(t_grid, s_j_true, exp(theta_true[3]) * lambda_true),
St(t_grid, s_j_true, exp(theta_true[4]) * lambda_true)),
TRTPN = factor(rep(key_arms, each = length(t_grid)), levels = key_arms)) %>%
ggplot() +
geom_line(aes(x = X, y = y, col = TRTPN), size = 1) +
geom_vline(xintercept = s_j_true, lty = 3, size = 1) +
scale_x_continuous(name = 'Follow-up time (weeks)', expand = c(0, 0)) +
scale_y_continuous(name = TeX("$S(t)$"), limits = c(0, 1), expand = c(0, 0)) +
scale_color_manual(name = 'Arm', values = key_palette) +
theme(legend.position = 'top', text = element_text(size = 16))
# KM Plot
dat_KM = tibble(Duration = y,
Outcome = nu,
Dose = factor(dose_label))
fit_KM <- survfit(Surv(Duration, Outcome) ~ Dose, data = dat_KM)
p_km = ggsurvplot(fit_KM,
data = dat_KM,
size = 1,
legend.title = "Arm",
legend.labs = key_arms,
palette = key_palette,
axes.offset = F,
ggtheme = theme_bw(base_size = 16)
)
p_km$plot +
geom_vline(xintercept = s_j_true, size = 1, lty = 3) +
labs(x = 'Follow-up time (weeks)', y = TeX("$S(t)$"))
We can also inspect the accrual of the trial that we just simulated.
tibble(TIMESINCERAND = data$dates) %>%
mutate(Subject = row_number()) %>%
mutate(EXPACCR = Lambda(t = TIMESINCERAND, peak_rate = accr_profile$peak_rate,
ramp_up = accr_profile$ramp_up)) %>%
ggplot() +
geom_step(aes(x = TIMESINCERAND, y = Subject, col = 'observed'), size = 1) +
geom_line(aes(x = TIMESINCERAND, y = EXPACCR,
col = paste0(accr_profile$peak_rate, ' per week')),
size = 1) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
labs(x = 'Weeks since start of accrual',
y = 'Cumulative number of subjects',
title = 'Observed vs simulated accrual') +
scale_color_brewer(name = '', palette = 'Dark2') +
theme(legend.position = "top")
4 Model Fitting
When fitting the data, the first choice we are confronted with is the specification of the time intervals \(s_{1}, \dots, s_{J}\) in which the hazard rates are constant. At this stage we do not know the underlying mechanism that generated the data, and our guess is to use two intervals with a single cutpoint at \(t = 10\). Note that, even though we used four intervals in the data generation, the hazard rates were constant across the first two intervals and across the last two intervals. Therefore, choosing only two intervals falls within the same family that generated the data (model is not misspecified) and we should be able to recover the underlying true parameters.
model_seed = 11111
# Define interval breaks (in weeks)
eps = 1E-6
s_j = c(0, eps + c(10, 20))
J = length(s_j) - 1Note that we have added a small perturbation to s_j to ensure that no observation falls at the cutpoints (this can give rise to computational issues).
4.1 Prior Elicitation
We first outline a potential strategy to elicit the priors in this case study.
The baseline hazard rates \(\lambda_{j}\) are assumed to be constant within each segment \(j\), and modeled with independent weakly-informative prior distributions: \[\lambda_{j} \sim \text{Gamma}(a_{j}, b_{j}),\] where the \(\text{Gamma}(a, b)\) distribution has the following density function: \[f(\lambda) = \frac{b^{a}}{\Gamma(a)} \lambda^{a-1}\text{exp}(-b \lambda), \quad \lambda > 0\] with mean \(a / b\) and weight \(a\).
We can choose the prior on the baseline hazard rates using previous studies. For example, say that it is expected that \(40\%\) of patients in the control group will experience a primary endpoint event within the first \(20\) weeks of exposure. Assuming an exponential distribution, the weekly hazard rate corresponding to \(40\%\) cumulative event rate by \(20\) weeks is equal to \[ \begin{align*} &S(20) = 1 - 0.4 = \exp\{-H(20)\} = \exp\{-20 \lambda\} \\ \Rightarrow \quad &\lambda = - \dfrac{\log (1 - 0.4)}{20} = 0.0255. \end{align*} \] Therefore the prior for the event rate should be centered at \(0.0255\) for the first \(20\) weeks. If we are not very confident in this prior mean, the event rates could be given very little weight. For example, choosing \(a_{j} = 1\) and \(b_{j} = 39.2156863\) is equivalent to observing \(1\) event with \(39.1523\) weeks of exposure. The piecewise exponential model provides the flexibility to allow for the control hazard rate to change over the time period, but a priori the distributions are the same for each time segment.
The parameter \(\tau^{2}\) controls the degree of smoothing between doses, and is estimated using an inverse gamma prior: \[\tau^{2} \sim \text{inv-Gamma}\left(\dfrac{\alpha}{2} , \dfrac{\alpha \beta^{2}}{2} \right),\] where the \(\text{inv-Gamma}(\alpha / 2, \alpha \beta^{2} / 2)\) distribution has the following density function: \[f(\tau^{2}) = \frac{(\alpha / 2)^{\alpha / 2} \beta^{\alpha}}{\Gamma(\alpha / 2)} \left( \frac{1}{\tau^{2}} \right)^{\alpha/2+1}\text{exp}\left( -\frac{\alpha \beta^{2}}{2 \tau^{2}} \right), \quad \tau^{2} > 0\] with central value \(\beta\) and weight \(\alpha\).
We can choose a non-informative prior so that the observed data very quickly overcomes its effect. For example, \(\alpha = 1, \beta = 1\) assumes a weight of 1 and central value for the variance parameter of 1.
# Prior hazard hyperparameters
a_j = rep(1, J)
b_j = rep(1/round(- log (1 - 0.4) / 20, 4), J)
alpha = 1
beta = 14.2 Posterior Sampling
When dealing with censored data, we observe the minimum time between the event time and the censoring time, i.e. \(y_{i} = \min\{T_{i}, C_{i}\}\), and the censoring indicator \(\delta_{i} = 1(T_{i} \leq C_{i})\).
In general, the likelihood corresponding to the \(i^{th}\) observation in survival analysis is
\[ L_{i} = \lambda(y_{i} )^{\delta_{i}} S(y_{i}). \]
In our specific model, we have an analytic expression for the hazard function and for the survival function. Therefore, the likelihood can be rewritten as \[ \begin{align*} L_{i} &= \left\{ e^{\mathbf{x}_{i}^T \boldsymbol{\theta}} \lambda_{0} (y_{i}) \right\}^{\delta_{i}} S_{0}(y_{i})^{e^{\mathbf{x}_{i}^T \boldsymbol{\theta}}} \\ &= \left( e^{\mathbf{x}_{i}^T \boldsymbol{\theta}} \lambda_{i^{\star}} \right)^{\delta_{i}} \exp\left\{ - e^{\mathbf{x}_{i}^T \boldsymbol{\theta}} \sum_{j = 1}^{i^{\star}} H_{ij} \lambda_{j} \right\}, \end{align*} \] where \(i^{\star}\) denotes the interval in which the time \(y_{i}\) (either censoring time or event time) falls, and \(H_{ij}\) denotes the exposure of subject \(i\) in each segment \(j\), i.e. \[H_{ij} = \begin{cases} s_{j} - s_{j-1} & \text{if } y_{i} \notin I_{j} \\ y_{i} - s_{j-1} & \text{if } y_{i} \in I_{j}. \end{cases}\] In the last equality, we used the expression for the survival function that we had previously described in (2.1).
We can use the function get_summary_stats() to calculate the important summary statistics \((H_{ij}, i^{\star})\), i.e. exposure in each segment for each patient and interval event indicator for each patient.
These will be passed into STAN and will be used to calculate the likelihood.
# Calculate summary statistics
summ_stats = get_summary_stats(y, nu, s_j)
# Put data in list for STAN
stan_data <- list(N = N,
X = X,
D = D,
J = J,
delta = summ_stats$delta_i,
H = summ_stats$Hij,
a_j = a_j,
b_j = b_j,
alpha = alpha,
beta = beta
)We can run the following STAN model.
writeLines(readLines("./TTE_pwc.stan"))functions{
// Likelihood function
real pc_haz(vector lambda, vector beta, int[] delta, matrix H, matrix X){
vector[num_elements(delta)] loglik_i;
real loglik;
vector[num_elements(delta)] Lambda;
vector[num_elements(delta)] eta;
eta = X * beta;
Lambda = - (H * lambda) .* exp(eta);
loglik_i = Lambda;
for (i in 1:num_elements(delta)) {
if (delta[i] > 0){
loglik_i[i] = loglik_i[i] + log(lambda[delta[i]]) + eta[i];
}
}
loglik = sum(loglik_i);
return loglik;
}
}
data{
int<lower=1> N; // number of patients
int<lower=1> D; // number of doses
matrix[N,D] X; // design matrix (without column of ones)
int<lower=1> J; // number of intervals
int delta[N]; // vector indicating in which interval the event happens for each patient (0 if it does never happen)
matrix[N,J] H; // matrix with the sufficient statistics to calculate the cumulative hazard
vector<lower=0>[J] a_j; // Weight hyperparameters for each component of the gamma prior on lambda
vector<lower=0>[J] b_j; // Mean hyperparameters for each component of the gamma prior on lambda
real<lower=0> alpha; // Weight hyperparameter for the inv-gamma prior on tau^2
real<lower=0> beta; // Central value hyperparameter for the inv-gamma prior on tau^2
}
parameters{
vector<lower=0>[J] lambda; // piecewise constant baseline hazards
vector[D] theta; // covariate effects
real<lower=0> tau2;
}
transformed parameters {
real<lower=0> tau = sqrt(tau2);
}
model{
// Priors
lambda ~ gamma(a_j, b_j); // prior for the piecewise constant baseline hazards
tau2 ~ scaled_inv_chi_square(alpha, beta);
theta[1] ~ normal(0, 1);
for (d in 2:D){
theta[d] ~ normal(theta[d-1], tau);
}
// Likelihood
target += pc_haz(lambda, theta, delta, H, X);
}The main part of this model is the likelihood calculation via the function pc_haz().
This simply performs the likelihood update on the log scale, as described above.
We can now run the model and put the results (MCMC chains) into a matrix whose columns correspond to the different model parameters \(\lambda_{1}, \dots, \lambda_{J}, \theta_{1}, \dots, \theta_{D}, \tau^{2}\).
# Run STAN model
n_iters = 5000
burnin = 2500
n_chains = 2
set.seed(model_seed)
fit <- stan(file = './TTE_pwc.stan',
data = stan_data,
chains = n_chains,
warmup = burnin,
iter = n_iters,
seed = model_seed) %>%
extract()
# Extract model output
samples = fit$lambda %>%
cbind(fit$theta) %>%
cbind(fit$tau2) %>%
colnames_inplace(c(sprintf('lambda[%s]', 1:J), sprintf('theta[%s]', 1:D), 'tau^2'))4.3 MCMC Diagnostics
The traceplots for the model parameters do not show any concerning behavior.
# Traceplots
mcmc_trace(samples, facet_args = list(labeller = ggplot2::label_parsed,
ncol = 2))
The total sample size of the chains was 5000, and the effective sample sizes are very close to their optimal values.
# Effective sample size
print(round(coda::effectiveSize(samples)))lambda[1] lambda[2] theta[1] theta[2] theta[3] theta[4] tau^2
5240 5000 5205 5000 5559 5017 5000 For some parameters, the the effective sample sizes are larger than the actual sample sizes. As detailed in the STAN guide, this can happen for parameters which have close to Gaussian posterior and little dependency on other parameters.
5 Posterior Checks
Once we have the MCMC draws, we can calculate any functional of interest, e.g. posterior mean, median, quantiles, etc.
samples %>%
as_tibble %>%
pivot_longer(cols = everything(), names_to = 'Parameter') %>%
group_by(Parameter) %>%
summarize(Mean = mean(value),
Sd = sd(value),
Median = median(value),
Low = quantile(value, probs = 0.025),
Upp = quantile(value, probs = 0.975),
ESS = round(coda::effectiveSize(value))) %>%
mutate(Mean = sprintf("$%0.3f$", Mean),
Sd = sprintf("$%0.3f$", Sd),
MCI = sprintf("$%0.3f$ ($%0.3f$, $%0.3f$)", Median, Low, Upp),
ESS = sprintf("$%s$", ESS)) %>%
select(Parameter, Mean, `Std. Dev.` = Sd, `Median (95% CI)` = MCI,
`Effective Sample Size` = ESS) %>%
mutate(Parameter = str_replace_all(Parameter, '\\[', '\\_'),
Parameter = str_remove_all(Parameter, '\\]'),
Parameter = paste0("$\\", Parameter, "$")) %>%
knitr::kable(booktabs = TRUE,
caption = 'Posterior estimates of the model parameters.') %>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"))| Parameter | Mean | Std. Dev. | Median (95% CI) | Effective Sample Size |
|---|---|---|---|---|
| \(\lambda_1\) | \(0.066\) | \(0.010\) | \(0.066\) (\(0.048\), \(0.087\)) | \(5240\) |
| \(\lambda_2\) | \(0.029\) | \(0.006\) | \(0.028\) (\(0.018\), \(0.043\)) | \(5000\) |
| \(\tau^2\) | \(0.693\) | \(1.429\) | \(0.395\) (\(0.105\), \(2.921\)) | \(5000\) |
| \(\theta_1\) | \(0.340\) | \(0.224\) | \(0.342\) (\(-0.102\), \(0.773\)) | \(5205\) |
| \(\theta_2\) | \(0.737\) | \(0.215\) | \(0.738\) (\(0.312\), \(1.160\)) | \(5000\) |
| \(\theta_3\) | \(0.724\) | \(0.214\) | \(0.724\) (\(0.305\), \(1.129\)) | \(5559\) |
| \(\theta_4\) | \(0.596\) | \(0.221\) | \(0.597\) (\(0.159\), \(1.025\)) | \(5017\) |
We can calculate the survival curves associated to each MCMC draw. Averaging across these sampled survival curves will give us an estimate (posterior mean) which hopefully matches the data well.
# Survival curves for each treatment arm
post_surv_ctr = St(t = t_grid,
s_j = s_j,
lambda = samples[,1:J])
post_surv_trt1 = St(t = t_grid,
s_j = s_j,
lambda = exp(samples[,J+1]) * samples[,1:J])
post_surv_trt2 = St(t = t_grid,
s_j = s_j,
lambda = exp(samples[,J+2]) * samples[,1:J])
post_surv_trt3 = St(t = t_grid,
s_j = s_j,
lambda = exp(samples[,J+3]) * samples[,1:J])
post_surv_trt4 = St(t = t_grid,
s_j = s_j,
lambda = exp(samples[,J+4]) * samples[,1:J])
post_survival = post_surv_ctr %>%
reshape2::melt(varnames = c('t', 'iteration')) %>%
mutate(t = t_grid[t],
TRTPN = key_arms[1]) %>%
bind_rows(post_surv_trt1 %>%
reshape2::melt(varnames = c('t', 'iteration')) %>%
mutate(t = t_grid[t],
TRTPN = key_arms[2])) %>%
bind_rows(post_surv_trt2 %>%
reshape2::melt(varnames = c('t', 'iteration')) %>%
mutate(t = t_grid[t],
TRTPN = key_arms[3])) %>%
bind_rows(post_surv_trt3 %>%
reshape2::melt(varnames = c('t', 'iteration')) %>%
mutate(t = t_grid[t],
TRTPN = key_arms[4])) %>%
bind_rows(post_surv_trt4 %>%
reshape2::melt(varnames = c('t', 'iteration')) %>%
mutate(t = t_grid[t],
TRTPN = key_arms[5])) %>%
group_by(t, TRTPN) %>%
summarize(mean = mean(value),
low = quantile(value, probs = 0.025),
upp = quantile(value, probs = 0.975)) %>%
ungroup() %>%
mutate(TRTPN = factor(TRTPN, levels = key_arms))
# KM Estimate
dat_KM = tibble(ExpTime = y, PRIM = nu, Dose = factor(key_arms[dose_label + 1]))
fit_KM <- survfit(Surv(ExpTime, PRIM) ~ Dose, data = dat_KM)
p_km = ggsurvplot(fit_KM,
data = dat_KM,
palette = key_palette,
legend = "top",
legend.title = "Arm",
legend.labs = levels(dat_KM$Dose),
axes.offset = FALSE,
break.time.by = 10,
ggtheme = custom_theme()
)
# Posterior survival curves and KM plot
p_km = p_km$plot +
geom_line(data = post_survival, aes(x = t, y = mean, col = TRTPN), size = 1, lty = 4) +
# geom_ribbon(data = post_survival,
# aes(x = t, ymin = low, ymax = upp, fill = TRTPN),
# alpha = 0.25, inherit.aes = FALSE) +
geom_vline(xintercept = s_j, lty = 3) +
scale_y_continuous(expand = c(0.001, 0.001)) +
labs(x = 'Follow-up time (weeks)', y = 'Survival probability') +
theme(legend.position = 'top')
p_km
We see that the model fit seems to be in agreement with the Kaplan-Meier curves. We can also compare the distribution of the hazard ratios with the true underlying values that we used to simulate the data.
# Dose response curve (HR density for each dose)
samples %>%
as_tibble() %>%
select(contains('theta')) %>%
add_column('theta[0]' = 0, .before = 'theta[1]') %>%
pivot_longer(cols = everything(), names_to = 'TRTPN', values_to = 'HR',
names_prefix = 'theta') %>%
mutate(TRTPN = as.numeric(str_remove_all(TRTPN, "\\[|\\]")),
HR = exp(HR),
TRTPN = factor(key_arms[as.numeric(TRTPN) + 1], levels = key_arms)) %>%
group_by(TRTPN) %>%
summarize(low = quantile(HR, probs = 0.025),
upp = quantile(HR, probs = 0.975),
mean = mean(HR),
median = median(HR)) %>%
ggplot() +
geom_point(aes(x = TRTPN, y = mean), pch = 10, size = 4, stroke = 1) +
geom_line(aes(x = TRTPN, y = mean, group = 1), size = 1) +
geom_line(aes(x = TRTPN, y = low, group = 1), lty = 3) +
geom_line(aes(x = TRTPN, y = upp, group = 1), lty = 3) +
geom_point(aes(x = 1:5, y = c(1, exp(theta_true))), col = 'red') +
geom_ribbon(aes(x = TRTPN, ymin = low, ymax = upp, group = 1), alpha = 0.25) +
labs(x = '', y = 'Hazard ratios')
Lastly, we also show that the underlying baseline hazard rates can be recovered well by our model implementation.
# Baseline hazard function
reshape2::melt(ht(t_grid, s_j, samples[,1:J]), varnames = c('Iteration', 'X')) %>%
mutate(X = t_grid[X]) %>%
group_by(X) %>%
dplyr::summarize(ybar = mean(value),
UI = quantile(value, probs = 0.975),
LI = quantile(value, probs = 0.025)) %>%
ggplot() +
geom_line(aes(x = X, y = ybar), size = 1) +
geom_ribbon(aes(x = X, ymin = LI, ymax = UI), alpha = 0.3) +
geom_step(data = data.frame(s_j_true = s_j_true,
lambda_true = c(lambda_true, tail(lambda_true, 1))),
aes(x = s_j_true, y = lambda_true), size = 1, col = 'red') +
geom_vline(xintercept = s_j, lty = 3) +
scale_x_continuous(name = 'time (weeks)', expand = c(0, 0)) +
labs(y = TeX("$\\lambda_{j}$"))
6 Sensitivity Analyses
We describe here possible additional analyses that might be performed when implementing a trial.
6.1 Smoothing Prior Hyperparameters
We repeat the analysis described in Section 4 with a different prior for the smoothing parameter \(\tau^{2}\), favoring more borrowing of information across doses. We will be using \(\alpha = 10\), \(\beta = 0.1\).
# Define interval breaks (in weeks)
s_j = c(0, eps + c(10, 20))
J = length(s_j) - 1
# Prior hazard hyperparameters
a_j = rep(1, J)
b_j = rep(1/round(- log (1 - 0.4) / 20, 4), J)
alpha = 10
beta = 0.1
# Calculate summary statistics
summ_stats = get_summary_stats(y, nu, s_j)
# Put data in list for STAN
stan_data <- list(N = N,
X = X,
D = D,
J = J,
delta = summ_stats$delta_i,
H = summ_stats$Hij,
a_j = a_j,
b_j = b_j,
alpha = alpha,
beta = beta
)
# Run STAN model
set.seed(model_seed)
fit <- stan(file = './TTE_pwc.stan',
data = stan_data,
chains = n_chains,
warmup = burnin,
iter = n_iters,
seed = model_seed) %>%
extract()
# Extract model output
samples = fit$lambda %>%
cbind(fit$theta) %>%
cbind(fit$tau2) %>%
colnames_inplace(c(sprintf('lambda[%s]', 1:J), sprintf('theta[%s]', 1:D), 'tau^2'))
# Dose response curve (HR density for each dose)
samples %>%
as_tibble() %>%
select(contains('theta')) %>%
add_column('theta[0]' = 0, .before = 'theta[1]') %>%
pivot_longer(cols = everything(), names_to = 'TRTPN', values_to = 'HR',
names_prefix = 'theta') %>%
mutate(TRTPN = as.numeric(str_remove_all(TRTPN, "\\[|\\]")),
HR = exp(HR),
TRTPN = factor(key_arms[as.numeric(TRTPN) + 1], levels = key_arms)) %>%
group_by(TRTPN) %>%
summarize(low = quantile(HR, probs = 0.025),
upp = quantile(HR, probs = 0.975),
mean = mean(HR),
median = median(HR)) %>%
ggplot() +
geom_point(aes(x = TRTPN, y = mean), pch = 10, size = 4, stroke = 1) +
geom_line(aes(x = TRTPN, y = mean, group = 1), size = 1) +
geom_line(aes(x = TRTPN, y = low, group = 1), lty = 3) +
geom_line(aes(x = TRTPN, y = upp, group = 1), lty = 3) +
geom_point(aes(x = 1:5, y = c(1, exp(theta_true))), col = 'red') +
geom_ribbon(aes(x = TRTPN, ymin = low, ymax = upp, group = 1), alpha = 0.25) +
labs(x = '', y = 'Hazard ratios')
Clearly, the dose response curve looks smoother than the one obtained with a vague prior. In this case, we know the underlying true parameters and we realize that this extreme borrowing might be a bit excessive.
6.2 Number of Cutpoints
We would like to make as few assumptions as possible about the hazard function in order to avoid incorrect inference on the hazard ratios for the treatment effects. Towards a non-parametric approach, a possible way to do so is to restrict the intervals in which the hazard function is constant. Note that this additional flexibility does not come at no cost.
In particular:
- Numerical instability and unidentifiability issues are possible when no observations fall in some of the time intervals. This can be mitigated by the prior choice, but in that case the prior effect is not overwhelmed by the data
- The proportional hazard assumption, that in practice can be the most restrictive, still holds. This implies that survival curves for different doses cannot cross
We here repeat the experiment that we ran in Section 4 with \(8\) intervals between \(0\) and \(20\) weeks.
# Define interval breaks (in weeks)
s_j = c(0, eps + seq(2.5, 20, by = 2.5))
J = length(s_j) - 1
# Prior hazard hyperparameters
a_j = rep(1, J)
b_j = rep(1/round(- log (1 - 0.4) / 20, 4), J)
alpha = 1
beta = 1
# Calculate summary statistics
summ_stats = get_summary_stats(y, nu, s_j)
# Put data in list for STAN
stan_data <- list(N = N,
X = X,
D = D,
J = J,
delta = summ_stats$delta_i,
H = summ_stats$Hij,
a_j = a_j,
b_j = b_j,
alpha = alpha,
beta = beta
)
# Run STAN model
set.seed(model_seed)
fit <- stan(file = './TTE_pwc.stan',
data = stan_data,
chains = n_chains,
warmup = burnin,
iter = n_iters,
seed = model_seed) %>%
extract()
# Extract model output
samples = fit$lambda %>%
cbind(fit$theta) %>%
cbind(fit$tau2) %>%
colnames_inplace(c(sprintf('lambda[%s]', 1:J), sprintf('theta[%s]', 1:D), 'tau^2'))
# Survival curves for each treatment arm
post_surv_ctr = St(t = t_grid,
s_j = s_j,
lambda = samples[,1:J])
post_surv_trt1 = St(t = t_grid,
s_j = s_j,
lambda = exp(samples[,J+1]) * samples[,1:J])
post_surv_trt2 = St(t = t_grid,
s_j = s_j,
lambda = exp(samples[,J+2]) * samples[,1:J])
post_surv_trt3 = St(t = t_grid,
s_j = s_j,
lambda = exp(samples[,J+3]) * samples[,1:J])
post_surv_trt4 = St(t = t_grid,
s_j = s_j,
lambda = exp(samples[,J+4]) * samples[,1:J])
post_survival = post_surv_ctr %>%
reshape2::melt(varnames = c('t', 'iteration')) %>%
mutate(t = t_grid[t],
TRTPN = key_arms[1]) %>%
bind_rows(post_surv_trt1 %>%
reshape2::melt(varnames = c('t', 'iteration')) %>%
mutate(t = t_grid[t],
TRTPN = key_arms[2])) %>%
bind_rows(post_surv_trt2 %>%
reshape2::melt(varnames = c('t', 'iteration')) %>%
mutate(t = t_grid[t],
TRTPN = key_arms[3])) %>%
bind_rows(post_surv_trt3 %>%
reshape2::melt(varnames = c('t', 'iteration')) %>%
mutate(t = t_grid[t],
TRTPN = key_arms[4])) %>%
bind_rows(post_surv_trt4 %>%
reshape2::melt(varnames = c('t', 'iteration')) %>%
mutate(t = t_grid[t],
TRTPN = key_arms[5])) %>%
group_by(t, TRTPN) %>%
summarize(mean = mean(value),
low = quantile(value, probs = 0.025),
upp = quantile(value, probs = 0.975)) %>%
ungroup() %>%
mutate(TRTPN = factor(TRTPN, levels = key_arms))
# KM Estimate
dat_KM = tibble(ExpTime = y, PRIM = nu, Dose = factor(key_arms[dose_label + 1]))
fit_KM <- survfit(Surv(ExpTime, PRIM) ~ Dose, data = dat_KM)
p_km = ggsurvplot(fit_KM,
data = dat_KM,
palette = key_palette,
legend = "top",
legend.title = "Arm",
legend.labs = levels(dat_KM$Dose),
axes.offset = FALSE,
break.time.by = 2.5,
ggtheme = custom_theme()
)
# Posterior survival curves and KM plot
p_km = p_km$plot +
geom_line(data = post_survival, aes(x = t, y = mean, col = TRTPN), size = 1, lty = 4) +
# geom_ribbon(data = post_survival,
# aes(x = t, ymin = low, ymax = upp, fill = TRTPN),
# alpha = 0.25, inherit.aes = FALSE) +
geom_vline(xintercept = s_j, lty = 3) +
scale_y_continuous(expand = c(0.001, 0.001)) +
labs(x = 'Follow-up time (weeks)', y = 'Survival probability') +
theme(legend.position = 'top')
p_km
# Dose response curve (HR density for each dose)
samples %>%
as_tibble() %>%
select(contains('theta')) %>%
add_column('theta[0]' = 0, .before = 'theta[1]') %>%
pivot_longer(cols = everything(), names_to = 'TRTPN', values_to = 'HR',
names_prefix = 'theta') %>%
mutate(TRTPN = as.numeric(str_remove_all(TRTPN, "\\[|\\]")),
HR = exp(HR),
TRTPN = factor(key_arms[as.numeric(TRTPN) + 1], levels = key_arms)) %>%
group_by(TRTPN) %>%
summarize(low = quantile(HR, probs = 0.025),
upp = quantile(HR, probs = 0.975),
mean = mean(HR),
median = median(HR)) %>%
ggplot() +
geom_point(aes(x = TRTPN, y = mean), pch = 10, size = 4, stroke = 1) +
geom_line(aes(x = TRTPN, y = mean, group = 1), size = 1) +
geom_line(aes(x = TRTPN, y = low, group = 1), lty = 3) +
geom_line(aes(x = TRTPN, y = upp, group = 1), lty = 3) +
geom_point(aes(x = 1:5, y = c(1, exp(theta_true))), col = 'red') +
geom_ribbon(aes(x = TRTPN, ymin = low, ymax = upp, group = 1), alpha = 0.25) +
labs(x = '', y = 'Hazard ratios')
# Baseline hazard function
reshape2::melt(ht(t_grid, s_j, samples[,1:J]), varnames = c('Iteration', 'X')) %>%
mutate(X = t_grid[X]) %>%
group_by(X) %>%
dplyr::summarize(ybar = mean(value),
UI = quantile(value, probs = 0.975),
LI = quantile(value, probs = 0.025)) %>%
ggplot() +
geom_line(aes(x = X, y = ybar), size = 1) +
geom_ribbon(aes(x = X, ymin = LI, ymax = UI), alpha = 0.3) +
geom_step(data = data.frame(s_j_true = s_j_true,
lambda_true = c(lambda_true, tail(lambda_true, 1))),
aes(x = s_j_true, y = lambda_true), size = 1, col = 'red') +
geom_vline(xintercept = s_j, lty = 3) +
scale_x_continuous(name = 'time (weeks)', expand = c(0, 0)) +
labs(y = TeX("$\\lambda_{j}$"))
The estimates still look reasonable for this (slightly) overparametrized model.
7 Conclusions
Survival modeling is a powerful tool for modeling the occurrence of events, whether events are negative such as death or beneficial such as progression. The famous Cox proportional hazards model (Cox 1972) is an extension that we did not cover here which involves non-parametric modeling of the hazard function.
References
Original Computing Environment
writeLines(readLines(file.path(Sys.getenv("HOME"), ".R/Makevars")))
CXX14FLAGS += -O3 -mtune=native -arch x86_64 -ftemplate-depth-256sessionInfo()R version 4.1.2 (2021-11-01)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 10.16
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] parallel stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] kableExtra_1.3.4 reshape2_1.4.4 survminer_0.4.9
[4] ggpubr_0.4.0 survival_3.3-1 coda_0.19-4
[7] bayesplot_1.9.0 rstan_2.21.5 StanHeaders_2.21.0-7
[10] RColorBrewer_1.1-3 latex2exp_0.9.4 ggplot2_3.3.6
[13] knitr_1.39 stringr_1.4.0 readr_2.1.2
[16] tidyr_1.2.0 dplyr_1.0.9 tibble_3.1.7
loaded via a namespace (and not attached):
[1] matrixStats_0.62.0 webshot_0.5.3 httr_1.4.3
[4] tensorA_0.36.2 tools_4.1.2 backports_1.4.1
[7] bslib_0.3.1 utf8_1.2.2 R6_2.5.1
[10] DBI_1.1.2 colorspace_2.0-3 withr_2.5.0
[13] tidyselect_1.1.2 gridExtra_2.3 prettyunits_1.1.1
[16] processx_3.5.3 compiler_4.1.2 cli_3.3.0
[19] rvest_1.0.2 xml2_1.3.3 posterior_1.2.1
[22] labeling_0.4.2 bookdown_0.26 sass_0.4.1
[25] checkmate_2.1.0 scales_1.2.0 survMisc_0.5.6
[28] ggridges_0.5.3 callr_3.7.0 systemfonts_1.0.4
[31] digest_0.6.29 rmarkdown_2.14 svglite_2.1.0
[34] pkgconfig_2.0.3 htmltools_0.5.2 highr_0.9
[37] fastmap_1.1.0 rlang_1.0.2 rstudioapi_0.13
[40] farver_2.1.0 jquerylib_0.1.4 generics_0.1.2
[43] zoo_1.8-10 jsonlite_1.8.0 distributional_0.3.0
[46] car_3.0-13 inline_0.3.19 magrittr_2.0.3
[49] loo_2.5.1 Matrix_1.4-1 Rcpp_1.0.8.3
[52] munsell_0.5.0 fansi_1.0.3 abind_1.4-5
[55] lifecycle_1.0.1 stringi_1.7.6 yaml_2.3.5
[58] carData_3.0-5 pkgbuild_1.3.1 plyr_1.8.7
[61] grid_4.1.2 crayon_1.5.1 lattice_0.20-45
[64] splines_4.1.2 hms_1.1.1 ps_1.7.0
[67] pillar_1.7.0 ggsignif_0.6.3 codetools_0.2-18
[70] stats4_4.1.2 glue_1.6.2 evaluate_0.15
[73] data.table_1.14.2 RcppParallel_5.1.5 vctrs_0.4.1
[76] tzdb_0.3.0 gtable_0.3.0 purrr_0.3.4
[79] km.ci_0.5-6 assertthat_0.2.1 xfun_0.31
[82] xtable_1.8-4 broom_0.8.0 rstatix_0.7.0
[85] viridisLite_0.4.0 KMsurv_0.1-5 ellipsis_0.3.2