--- title: "Time-to-Event Data MCP-Mod" output: rmarkdown::html_vignette bibliography: '`r system.file("REFERENCES.bib", package = "DoseFinding")`' csl: '`r system.file("american-statistical-association.csl", package = "DoseFinding")`' link-citations: yes vignette: > %\VignetteIndexEntry{Time-to-Event Data MCP-Mod} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, child = "children/settings.txt"} ``` In this vignette we illustrate how to use the `DoseFinding` package with time-to-event (survival) endpoints. We show how to simulate a dose-finding trial, how to apply the generalized MCP-Mod methodology to the log hazard ratios estimated from a Cox proportional hazards model, and how to perform power and sample size calculations following @deng2019. For continuously distributed data see the [analysis of normal data vignette][v2], and for binary data the [binary data vignette][v3]. [v2]: analysis_normal.html [v3]: binary_data.html ```{r, message = FALSE} library(DoseFinding) library(ggplot2) library(survival) ## set seed and ensure reproducibility across R versions set.seed(2010, kind = "Mersenne-Twister", sample.kind = "Rejection", normal.kind = "Inversion") ``` ## Example trial We will first generate an appropriate dataset that we can use to illustrate the analysis of time-to-event data. Assume a dose-finding study is planned for a hypothetical investigational treatment. The endpoint of interest is the time to a certain event. The treatment is tested with doses 0 (placebo), 0.05, 0.2, 0.5 and 1, with an allocation ratio of 2:1:1:1:2. The larger allocation to placebo and highest dose is chosen to increase the power. It is assumed that the hazard rate for the placebo group is $\lambda_0 = 0.6$ and that the hazard ratio for the maximum effect over the dose range is $\text{HR} = 0.5$. This results in a hazard rate of $\lambda_1 = 0.3$. We define the parameters for the example just described: ```{r} doses <- c(0, 0.05, 0.2, 0.5, 1) # placebo and 4 treatment groups lambda0 <- 0.6 # yearly hazard rate for the placebo group HR1 <- 0.5 # hazard ratio for the treatment effect at the highest dose lambda1 <- lambda0 * HR1 # hazard rate for the highest dose alRatio <- c(2, 1, 1, 1, 2) # allocation ratio across the dose groups ``` We use the following candidate set of dose-response models for the mean response: an $\text{E}_{\text{max}}$, a sigmoidal $\text{E}_{\text{max}}$, a quadratic and an exponential model. In the survival context a *decreasing* response (a lower hazard) typically corresponds to a beneficial treatment effect, so we specify `direction = "decreasing"`. The dose-response models are placed on the *log hazard* scale, because dose-response models in the `DoseFinding` package assume additive effects. Consequently: - `placEff` is the log hazard rate of the placebo arm, i.e. $\log(\lambda_0)$. - `maxEff` is the maximum treatment effect within the dose range, i.e. the difference between the response at the highest dose and `placEff`. On the log hazard scale this equals the log hazard ratio, since \[ \log(\text{HR}) = \log(\lambda_1) - \log(\lambda_0). \] ```{r, fig.width = 8, out.width = '100%'} mods <- Mods( emax = c(0.1, 0.5), sigEmax = c(0.5, 4), quadratic = -0.9, exponential = 0.2, doses = doses, direction = "decreasing", # decreasing (lower hazard) is beneficial placEff = log(lambda0), # log hazard rate of the placebo arm maxEff = log(lambda1) - log(lambda0) # log hazard ratio at the highest dose ) plotMods(mods, superpose = TRUE, ylab = "log(hazard rate)") ``` To illustrate the methodology we simulate a survival data set based on the trial assumptions and one of our assumed candidate models. The function below draws survival times from an exponential distribution under one of the candidate models. It supports - *event-driven censoring*: the study is run until a total of `etotal` events have been observed, and - optional *staggered study entry*: an entry (calendar) time is added to each patient's event time, the data are censored on the calendar-time scale at the time of the `etotal`-th event, and the observed follow-up time is then expressed on the study-time scale. ```{r} #' @param mods `Mods` object defining the candidate models to simulate from #' @param model name or index of the model whose hazards are used to simulate #' (refers to column names or column indices of getResp(mods)) #' @param alRatio allocation ratio across the dose groups #' @param n total sample size #' @param etotal number of events at which event-driven censoring is applied #' @param staggered logical, whether to use staggered study entry #' @return data.frame with columns group, time, status simulate_surv_data <- function(mods, model, alRatio, n, etotal, staggered = FALSE) { doses <- attr(mods, "doses") y <- getResp(mods) lambda0 <- exp(attr(mods, "placEff")) ## assign each patient to a dose group (log hazard rate) according to alRatio log_rate_resp <- rep(y[, model], rndDesign(alRatio/sum(alRatio), n)) ## entry (calendar) times entry <- if (staggered) { runif(n, 0, qexp(etotal / n, lambda0)) } else { rep(0, n) } ## event times on the calendar-time scale event <- entry + rexp(n, exp(log_rate_resp)) ## event-driven censoring at the time of the etotal-th event t_cut <- sort(event)[etotal] status <- as.numeric(event <= t_cut) obs_time <- pmin(event, t_cut) - entry # follow-up on the study-time scale data.frame( group = factor(names(log_rate_resp), levels = as.character(doses)), time = obs_time, status = status ) } ``` We simulate a trial of 300 patients using the hazards from the $\text{E}_{\text{max}}$ model (`model = 1`) and assuming we will require 130 events. ```{r} n <- 300 dat <- simulate_surv_data( mods, model = 1, alRatio = alRatio, n = n, etotal = 130, staggered = TRUE ) head(dat) ``` We can use a Kaplan-Meier plot to visualize the survival times by dose. ```{r, fig.width = 6, out.width = '100%'} km <- survfit(Surv(time, status) ~ group, data = dat) plot(km, col = seq_along(doses), lwd = 2, xlab = "Years", ylab = "Survival probability", main = "Kaplan-Meier estimator by dose") legend("topright", legend = paste("dose", doses), col = seq_along(doses), lwd = 2, bty = "n") ``` ## Analysis with MCP-Mod We will now walk through the typical steps of a MCP-Mod analysis using our simulated trial data. ### MCP step The generalized MCP-Mod approach requires an estimate $\hat\mu$ of the (in this case placebo-adjusted) dose-group effects together with its covariance matrix $\hat S$. We obtain both from a Cox proportional hazards model, which expresses the hazard as \[ \lambda(t \mid \text{dose}) = \lambda_0(t)\exp(\beta_1 x_1 + \dots + \beta_p x_p), \] where the coefficients $\beta$ are the log hazard ratios of each dose group relative to placebo. ```{r} coxfit <- coxph(Surv(time, status) ~ group, data = dat) summary(coxfit) coef_cox <- coef(coxfit) # placebo-adjusted log hazard ratios S_cox <- vcov(coxfit) # covariance matrix of the estimates ``` We now test for a dose-response signal with `MCTtest()`. Since the coxfit model returns placebo-adjusted estimates (there is no placebo coefficient), we use `placAdj = TRUE` and drop the placebo dose from `dose`, we test one-sided. ```{r} mct <- MCTtest( dose = doses[-1], # placebo-adjusted, so drop placebo resp = coef_cox, S = S_cox, models = mods, type = "general", alternative = "one.sided", placAdj = TRUE ) mct ``` The output reports, per candidate model, the optimal-contrast $t$-statistics together with multiplicity-adjusted $p$-values. Small adjusted $p$-values give evidence for a dose-response relationship. If we would have specified a significance level of 0.025 one-sided, then we can consider the null hypothesis of no dose-response effect versus placebo is rejected since several of the adjusted p-values are smaller than 0.025. ### Dose-response estimation Dose-response modeling can then proceed with a combination of bootstrapping and model averaging. For detailed explanations refer to the [vignette for analysis of continuous data][v2]. Fitting is done on the log-hazard ratio scale so we are using placebo-adjusted models. Since `maFitMod` does not support placebo-adjusted dose-response estimation we need to define a couple of helper functions to perform the model averaging. ```{r} one_bootstrap_prediction <- function(mu_hat, S_hat, doses, bounds, dose_seq) { sim <- drop(mvtnorm::rmvnorm(1, mu_hat, S_hat)) fit <- lapply(c("emax", "sigEmax", "exponential", "quadratic"), function(mod) fitMod(doses, sim, model = mod, S = S_hat, type = "general", bnds = bounds[[mod]], placAdj = TRUE)) index <- which.min(sapply(fit, gAIC)) pred <- predict(fit[[index]], doseSeq = dose_seq, predType = "effect-curve") return(pred) } ## bs_predictions is a doses x replications matrix, ## probs is a 4-element vector of increasing probabilities for the quantiles summarize_predictions <- function(bs_predictions, probs) { stopifnot(length(probs) == 4) med <- apply(bs_predictions, 1, median) quants <- apply(bs_predictions, 1, quantile, probs = probs) bs_df <- as.data.frame(cbind(med, t(quants))) names(bs_df) <- c("median", "low_out", "low_in", "high_in", "high_out") return(bs_df) } predict_and_plot <- function(mu_hat, S_hat, doses, dose_seq, n_rep) { bs_rep <- replicate( n_rep, one_bootstrap_prediction(mu_hat, S_hat, doses, defBnds(max(doses)), dose_seq)) bs_summary <- summarize_predictions(bs_rep, probs = c(0.025, 0.25, 0.75, 0.975)) bs_summary <- as.data.frame(exp(bs_summary)) # back to hazard ratio scale ci_half_width <- qnorm(0.975) * sqrt(diag(S_hat)) cox_summary <- data.frame(dose = doses, mu_hat = exp(mu_hat), low = exp(mu_hat - ci_half_width), high = exp(mu_hat + ci_half_width)) gg <- ggplot(cbind(bs_summary, dose_seq = dose_seq)) + geom_line(aes(dose_seq, median)) + geom_ribbon(aes(x = dose_seq, ymin = low_in, ymax = high_in), alpha = 0.2) + geom_ribbon(aes(x = dose_seq, ymin = low_out, ymax = high_out), alpha = 0.2) + geom_point(aes(dose, mu_hat), cox_summary) + scale_y_continuous(breaks = seq(0.5, 1.5, 0.25)) + geom_errorbar(aes(dose, ymin = low, ymax = high), cox_summary, width = 0, alpha = 0.5) + xlab("Dose") + ylab("Hazard ratio") + labs(title = "Bootstrap estimates for population hazard ratio", subtitle = "confidence levels 50% and 95%") return(gg) } ``` We then use these functions to plot the estimated dose-response curve together with the Cox model estimates on the original hazard ratio scale. ```{r} predict_and_plot(coef_cox, S_cox, doses[-1], dose_seq = seq(0, 1, 0.01), 1000) ``` ## Power and sample size calculations at design stage For power and sample size calculations in generalized MCP-Mod we need to make an assumption about the covariance matrix $S$ of the placebo-adjusted estimates under the alternative. @deng2019 derive an approximation for the covariance matrix of the estimated log hazard ratios for a Cox model which we will use in the following. Let $n_k$ be the number of patients in dose group $k$ ($k = 0$ for placebo), $D$ the total number of events, and $\beta_k$ the log hazard ratio of dose group $k$ vs. placebo (with $\beta_0 = 0$). Then, with \[ p_k = \frac{(n_k/n_0)\exp(\beta_k)}{1 + \sum_{i=1}^{K}(n_i/n_0)\exp(\beta_i)}, \] the covariance matrix $S$ of $\hat\beta = (\hat\beta_1, \dots, \hat\beta_K)$ is approximately \[ S \approx \frac{1}{D} \begin{pmatrix} p_1^{-1} + p_0^{-1} & p_0^{-1} & \cdots & p_0^{-1}\\ p_0^{-1} & p_2^{-1} + p_0^{-1} & \cdots & p_0^{-1}\\ \vdots & \vdots & \ddots & \vdots\\ p_0^{-1} & p_0^{-1} & \cdots & p_K^{-1} + p_0^{-1} \end{pmatrix}. \] Since this matrix depends on the total number of events rather than the total number of patients, this facilitates the use of event-driven dose-finding designs. Under the null hypothesis of no dose-response effect versus placebo ($\beta_k = 0$ for all $k$) the covariance matrix simplifies to a matrix $S_0$ that depends only on the allocation ratios to each arm and the total number of events, and is used to obtain the critical value for the multiple contrast test. The following function calculates both $S$ (one per candidate model) and $S_0$. Following @deng2019 we evaluate $S$ at half the assumed effect, `0.5 * beta`, an intermediate value between the null and the alternative. Then using $S$ and $S_0$ we compute the power under each candidate model. The function below builds the optimal contrast for each model (using its own $S$), derives the critical value from $S_0$, and evaluates the power for a given number of events with `powMCT()`. We also define a function that allows us to calculate the required number of events to achieve a target average power using binary search. ```{r} #' Power of the multiple contrast test for a time-to-event endpoint #' #' @param mods Set of dose-response models; an object of class `Mods`. #' @param alpha One-sided significance level. #' @param alRatio Positive allocation weights for all treatment groups, #' including placebo. #' @param etotal Positive number of events. #' #' @return A list containing: #' - `power_table`: power under each candidate model and the average power; #' - `power`: unrounded power under each candidate model; #' - `power_av`: unrounded average power; #' - `contMat`: matrix of optimal contrasts; #' - `critV`: critical value for the multiple contrast test. powMCT_TTE <- function(mods, alpha, alRatio, etotal) { y <- getResp(mods) doses <- attr(mods, "doses") n_group <- nrow(y) n_model <- ncol(y) n_active <- n_group - 1 model_names <- colnames(y) if (is.null(model_names)) { model_names <- paste0("model", seq_len(n_model)) } # Placebo-adjusted log hazard ratios. beta <- sweep(y, 2, y[1, ], FUN = "-") # Group probabilities evaluated at half the assumed treatment effect. weighted_hazards <- exp(0.5 * beta) * alRatio p <- sweep( weighted_hazards, 2, colSums(weighted_hazards), FUN = "/" ) # Covariance matrices under each candidate alternative. S <- array( 0, dim = c(n_active, n_active, n_model), dimnames = list( as.character(doses[-1]), as.character(doses[-1]), model_names ) ) for (i in seq_len(n_model)) { p0 <- p[1, i] S_i <- matrix( 1 / p0, nrow = n_active, ncol = n_active ) diag(S_i) <- 1 / p[-1, i] + 1 / p0 S[, , i] <- S_i / etotal } # Covariance matrix under the null hypothesis. p_null <- alRatio / sum(alRatio) p0_null <- p_null[1] S0 <- matrix( 1 / p0_null, nrow = n_active, ncol = n_active ) diag(S0) <- 1 / p_null[-1] + 1 / p0_null S0 <- S0 / etotal # Obtain the optimal contrast for each candidate model using the # covariance matrix associated with that model. contMat <- vapply( seq_len(n_model), FUN = function(i) { contrast_i <- optContr( mods, doses = doses[-1], S = S[, , i], placAdj = TRUE )[[1]] contrast_i[, i] }, FUN.VALUE = numeric(n_active) ) dimnames(contMat) <- list( as.character(doses[-1]), model_names ) # The null correlation matrix and critical value do not depend on which # candidate model is treated as the true model. covMat0 <- crossprod(contMat, S0 %*% contMat) corMat0 <- cov2cor(covMat0) integration_control <- mvtnorm.control() critV <- critVal( corMat0, alpha, df = 0, alternative = "one.sided", control = integration_control ) model_power <- numeric(n_model) for (i in seq_len(n_model)) { model_power[i] <- powMCT( contMat, alpha = alpha, altModels = mods, S = S[, , i], placAdj = TRUE, alternative = "one.sided", df = Inf, critV = critV, control = integration_control )[i] } average_power <- mean(model_power) power_table <- data.frame( model = c(model_names, "average"), power = round(c(model_power, average_power), 3), row.names = NULL ) list( power_table = power_table, power = stats::setNames(model_power, model_names), power_av = average_power, contMat = contMat, critV = critV ) } #' Required number of events for a target power #' #' @param mods Set of dose-response models; an object of class `Mods`. #' @param alpha One-sided significance level. #' @param alRatio Positive allocation weights for all treatment groups, #' including placebo. #' @param targPow Target average power. #' @param lowerE Lower bound for the number of events. #' @param upperE Upper bound for the number of events. #' #' @return A list containing: #' - `events`: smallest integer event count attaining the target power; #' - `power`: average power at that event count; #' - `power_table`: model-specific powers at that event count; #' - `target_power`: requested average power. eventsMCT_TTE <- function(mods, alpha, alRatio, targPow, lowerE = 1, upperE = 1000) { # alRatio = 1 makes sampSize() iterate the total number of events directly; # targFunc still uses the alRatio argument of eventsMCT_TTE (lexical scoping) events <- sampSize( upperN = upperE, lowerN = lowerE, targFunc = function(e) { powMCT_TTE(mods = mods, alpha = alpha, alRatio = alRatio, etotal = e)$power_av }, target = targPow, alRatio = 1, Ntype = "total" )$samp.size final_result <- powMCT_TTE( mods = mods, alpha = alpha, alRatio = alRatio, etotal = events ) list( events = events, power = final_result$power_av, power_table = final_result$power_table, target_power = targPow ) } ``` For example, the power for 120 events at a one-sided significance level of 0.025 is: ```{r} powMCT_TTE(mods, alpha = 0.025, alRatio = alRatio, etotal = 120) ``` And we can calculate the number of required events to achieve an average power of 80%: ```{r} eventsMCT_TTE(mods, alpha = 0.025, alRatio = alRatio, targPow = 0.8, lowerE = 100, upperE = 300) ``` ### Calculating required sample size for a trial For trial planning purposes we may want to determine the required sample size for the trial rather than the number of events. To connect the number of events to a sample size, @deng2019 use the relationship between the sample size and the expected number of events under exponentially distributed survival times. More specifically, assume that patients are recruited uniformly over $[0, T_r]$ (end of recruitment) with study end $T_s$, and that survival times in group $k$ are exponential with hazard $\lambda_k$. Then the probability that a patient in group $k$ experiences an event by the end of the study is $$q_k = 1 - \frac{\exp\{-\lambda_k(T_s-T_r)\}-\exp(-\lambda_kT_s)}{\lambda_kT_r}.$$ If $w_k=r_k/\sum_i r_i$ is the allocation proportion, the expected total number of events for a sample size $N$ is $$E(D\mid N) = N\sum_{k=0}^{K}w_kq_k.$$ We can use these formulas to calculate the expected number of events for each candidate model for a given sample size and determine the required sample size to achieve a specified target average power. Note that this function could also be easily modified to use the minimum or maximum power instead of the mean across candidate models. ```{r} #' Required sample size for target average power with a time-to-event endpoint #' #' For each candidate model, the expected number of events is calculated from #' that model's hazards. Model-specific powers are then evaluated at the #' corresponding expected event counts and averaged across candidate models. #' #' @param mods Set of dose-response models; an object of class `Mods`. #' @param alpha One-sided significance level. #' @param alRatio Positive allocation ratios for all treatment groups, #' including placebo. #' @param targPow Target average power across candidate models. #' @param end_rec End of the recruitment period, in years. #' @param end_stu End of the study, in years. #' @param lowerN Lower bound for the total sample size. #' @param upperN Upper bound for the total sample size. #' #' @return A list containing: #' - `sampSize`: smallest integer sample size attaining the target average power; #' - `power`: average power at the returned sample size; #' - `model_power`: model-specific powers at the returned sample size; #' - `expected_events`: expected events under each candidate model; #' - `event_probability`: overall event probability under each candidate model; #' - `event_probability_by_group`: group-specific event probabilities; #' - `allocation`: allocation proportions; #' - `power_table`: model-specific expected events and powers; #' - `target_power`: requested average power. sampSizeMCT_TTE <- function(mods, alpha, alRatio, targPow, end_rec, end_stu, lowerN = 1, upperN = 1000) { y <- as.matrix(getResp(mods)) n_group <- nrow(y) n_model <- ncol(y) model_names <- colnames(y) if (is.null(model_names)) { model_names <- paste0("model", seq_len(n_model)) } allocation <- alRatio / sum(alRatio) # Each column contains the group-specific hazards under one candidate model. lambda <- exp(y) # Group-specific probability of observing an event by the end of the study. # # This is equivalent to # # 1 - ( # exp(-lambda * (end_stu - end_rec)) - # exp(-lambda * end_stu) # ) / (lambda * end_rec) # # but expm1() is more stable when lambda * end_rec is small. event_probability_by_group <- 1 - exp(-lambda * (end_stu - end_rec)) * (-expm1(-lambda * end_rec)) / (lambda * end_rec) # Overall event probability under each candidate model. event_probability <- colSums( sweep( event_probability_by_group, MARGIN = 1, STATS = allocation, FUN = "*" ) ) names(event_probability) <- model_names colnames(event_probability_by_group) <- model_names # Evaluate average power for a specified integer sample size. evaluate_sample_size <- function(n) { expected_events <- n * event_probability model_power <- vapply( seq_len(n_model), FUN = function(j) { power_result <- powMCT_TTE( mods = mods, alpha = alpha, alRatio = alRatio, etotal = expected_events[j] ) unname(power_result$power[j]) }, FUN.VALUE = numeric(1) ) names(expected_events) <- model_names names(model_power) <- model_names list( expected_events = expected_events, model_power = model_power, power_av = mean(model_power) ) } # alRatio = 1 makes sampSize() iterate the total sample size directly; # evaluate_sample_size() still uses the alRatio argument of sampSizeMCT_TTE required_n <- sampSize( upperN = upperN, lowerN = lowerN, targFunc = function(n) evaluate_sample_size(n)$power_av, target = targPow, alRatio = 1, Ntype = "total" )$samp.size final_result <- evaluate_sample_size(required_n) power_table <- data.frame( model = c(model_names, "average"), expected_events = c( unname(final_result$expected_events), mean(final_result$expected_events) ), power = c( unname(final_result$model_power), final_result$power_av ), row.names = NULL ) list( sampSize = required_n, power = final_result$power_av, model_power = final_result$model_power, expected_events = final_result$expected_events, event_probability = event_probability, event_probability_by_group = event_probability_by_group, allocation = allocation, power_table = power_table, target_power = targPow ) } ``` Now we can calculate how many patients we would need to recruit to achieve an average power of 80%, assuming we would recruit for 2 years and we would end the study after 4 years. ```{r} sampSizeMCT_TTE(mods, alpha = 0.025, alRatio = alRatio, end_rec = 2, end_stu = 4, targPow = 0.8, lowerN = 100, upperN = 300) ``` ## References