In borealexander/tslrt: Function for treatment switching logrank test

knitr::opts_chunk$set(echo = TRUE)

Installation

The 'tslrt' package is available on GitHub and can be downloaded and installed by using the 'devtools' package

#install.packages("devtools")
library(devtools)
#install_github("borealexander/tslrt")
library(tslrt)

Load packages

To use the 'tslrt' package there is also other packages that are used

library(data.table)
library(survival)
library(ggplot2)
library(dplyr)

Example to simulate data and use tslrt

For now it exists 2 ways to simulate data with the 'tslrt' package. The first is to simulate the control and experimental arm as exponential distribution, but that at a specific delay time a proportion of patients in control switches to the experimental treatment. The other way the control and experimental arms are still exponentially distributed, but the control arm now has a probability that they switch arm after progression, where the distribution of progressions are exponential.

Example with delay

# Set parameters

n_c <- 100
n_e <- 100
median_c_1 <- 10
median_c_2 <- 20
median_e <- 20
p_change <- 0.9
delay <- 3
end_event <- 150
rec_period <- 12
rec_power <- 1

# Model as input
model_delay_switch <- list(n_c = n_c,
                           n_e = n_e,
                           median_c_1 = median_c_1,
                           median_c_2 = median_c_2,
                           delay = delay,
                           median_e = median_e,
                           p_change = p_change,
                           end_event = end_event,
                           rec_period = rec_period,
                           rec_power = rec_power)

# Simulate data
ex_sim_delay <- sim_crossover_delay(model = model_delay_switch)

# print out the top of the simulated data
head(ex_sim_delay)

We also plot the data to show how it looks like

# Plot data
fit_delay <- survfit( Surv(time, event) ~ group, data = ex_sim_delay)
plot(fit_delay, col = c("black", "blue"), mark.time = TRUE, pch = c(1, 2))
legend('topright', legend = c("Control", "Experimental"),
       col = c("black", "blue"), lty = c(1, 1), pch = c(1, 2))

There is also a function to plot create the analytical survival and hazard ratio function for this model

t <- seq(1,50)

# Hazard ratio
HR_switch <- HR_delay_switch(lambda_c_1 = log(2)/median_c_1,
                             lambda_c_2 = log(2)/median_c_2,
                             lambda_e = log(2)/median_e,
                             delay = delay,
                             p = p_change)

plot(t, HR_switch$HR(t), type = "l",
     xlab = "Time (months)", ylab = "HR")

We also plot the analytical survival curve in comparison to the data. Of course with a small sample size the data and the model will look a bit different, but the figure shows that when looking at the control model comparing switching and no switching there is a big change in how the survival curve looks like.

plot(fit_delay, col = c("black", "blue"), mark.time = TRUE, pch = c(1, 2))
lines(t, exp(-t*log(2)/median_c_1), col = "red", lty = 2)
lines(t, HR_switch$S(t), col = "black", lty = 2)
lines(t, exp(-t*log(2)/median_e), col = "blue", lty = 2)
legend('topright', legend = c("Control", "Experimental",
                              "Control model - no switching",
                              "Control model - with switching",
                              "Experimental model"),
       col = c("black", "blue", "red", "black", "blue"),
       lty = c(1, 1, 2, 2, 2), 
       pch = c(1, 2, NA, NA, NA))

We now want to use the simulated data to perform the (weighted) logrank test. The first step is to create a risk table

# Create risk table
delay_risk_table <- get_risk_table(ex_sim_delay)

head(delay_risk_table)

With the risk table we can then calculate the weights and plot the result. There is 3 (4?) types of ways to calculate the weights at the moment. The "logrank" is a standard logrank test with weight 1. The "theta" is the new weighted test and needs to have a hazard ratio function depending on time as input. The "fh", which is the Flemming-Harrington test which needs rho and gamma as input parameters.

# create risk table
delay_risk_table <- get_risk_table(ex_sim_delay)

# calculate weights
logrank_weights <- calculate_weights(delay_risk_table, method = "logrank")
delay_weights <- calculate_weights(delay_risk_table, method = "theta", hr_fun = HR_switch$HR)
fh_weights <- calculate_weights(delay_risk_table, method = "fh", rho = 1, gamma = 0)


# data table for plotting the weights
weight_dt <- data.table(t  = c(logrank_weights$t, delay_weights$t, fh_weights$t),
                        w = c(logrank_weights$w, delay_weights$w, fh_weights$w),
                        Test = rep(c('logrank', 'new weights', 'fh(1,0)'),
                                   c(length(logrank_weights$w),length(delay_weights$w),length(fh_weights$w))))


ggplot(weight_dt, aes(x = t, y = w, color = Test)) +
  geom_line(size = 1.2) +
  theme_bw() +
  labs(x = "Time (months)", y = "w", title = "Comparison of weights for different tests") +
  theme(text = element_text(size = 18),
        panel.grid.minor = element_blank())

When the weights have been calculated, the last step is to caulculate the Z-scores for the tests

# Calculate Z-scores

# Logrank
calculate_zs(logrank_weights)

# New weights
calculate_zs(delay_weights)

# Flemming-Harrington
calculate_zs(fh_weights)

Example with exponential progression switching

For the model with treatment switching after progression, where the progression is exponentially distributed, we can use the 'tslrt' pacakge in the same ways as above, but with different functions for simulating data and create the HR function.

# Set parameters

n_c <- 100
n_e <- 100
median_c <- 10
median_e <- 20
median_prog <- 5
p_change <- 0.9
end_event <- 150
rec_period <- 12
rec_power <- 1

# Model as input
model_prog_switch <- list(n_c = n_c,
                           n_e = n_e,
                           median_c = median_c,
                           median_progression = median_prog,
                           median_e = median_e,
                           p_change = p_change,
                           end_event = end_event,
                           rec_period = rec_period,
                           rec_power = rec_power)

# Simulate data
ex_sim_prog <- sim_exp_crossover(model = model_prog_switch)

head(ex_sim_prog)

We plot the data

# Plot data
fit_prog <- survfit( Surv(time, event) ~ group, data = ex_sim_prog)
plot(fit_prog, col = c("black", "blue"), mark.time = TRUE, pch = c(1, 2))
legend('topright', legend = c("Control", "Experimental"),
       col = c("black", "blue"), lty = c(1, 1), pch = c(1, 2))

For this model we have the survival and hazard ratio function

t <- seq(1,50)

# Hazard ratio
HR_switch <- exp_prog_switch_fun(lambda_c = log(2)/median_c,
                                 lambda_e = log(2)/median_e,
                                 lambda_p = log(2)/median_prog,
                                 p = p_change)

plot(t, HR_switch$HR(t), type = "l",
     xlab = "Time (months)", ylab = "HR")

We can also plot the model in comparison to the data. Using a small sample size the data and model will differ a bit, but note the difference in control when we add switching.

plot(fit_prog, col = c("black", "blue"), mark.time = TRUE, pch = c(1, 2))
lines(t, exp(-t*log(2)/median_c), col = "red", lty = 2)
lines(t, HR_switch$S(t), col = "black", lty = 2)
lines(t, exp(-t*log(2)/median_e), col = "blue", lty = 2)
legend('topright', legend = c("Control", "Experimental",
                              "Control model - no switching",
                              "Control model - with switching",
                              "Experimental model"),
       col = c("black", "blue", "red", "black", "blue"),
       lty = c(1, 1, 2, 2, 2), 
       pch = c(1, 2, NA, NA, NA))

To use the weighted logrank test we start by creating a risk table

# Create risk table
prog_risk_table <- get_risk_table(ex_sim_prog)

head(prog_risk_table)

With the risk table we can then calculate the different versions of the logrank tests in the sam was as above. We also plot how the weights will differ between the different methods.

# calculate weights
logrank_weights <- calculate_weights(prog_risk_table, method = "logrank")
prog_weights <- calculate_weights(prog_risk_table, method = "theta", hr_fun = HR_switch$HR)
fh_weights <- calculate_weights(prog_risk_table, method = "fh", rho = 1, gamma = 0)

# data table for plotting the weights
weight_dt <- data.table(t = c(logrank_weights$t, prog_weights$t, fh_weights$t),
                        w = c(logrank_weights$w, prog_weights$w, fh_weights$w),
                        Test = rep(c('logrank', 'new weights', 'fh(1,0)'),
                                   c(length(logrank_weights$w),length(prog_weights$w),length(fh_weights$w))))

ggplot(weight_dt, aes(x = t, y = w, color = Test)) +
  geom_line(size = 1.2) +
  theme_bw() +
  labs(x = "Time (months)", y = "w", title = "Comparison of weights for different tests") +
  theme(text = element_text(size = 18),
        panel.grid.minor = element_blank())

The last step is to calculate the Z-scores for the different tests.

# Calculate Z-scores

# Logrank
calculate_zs(logrank_weights)

# New weights
calculate_zs(prog_weights)

# Flemming-Harrington
calculate_zs(fh_weights)

Use tslrt on real data

To use the 'tslrt' package on real data it is needed a data set that contains time for event, an indicator for event of censoring and group (either "control" or "experimental"), and a hazard ratio function depending on time for calculating the weights.

example_data <- data.table(OS_time = rexp(100, rate = 0.1),
                           event_indicator = sample(c(FALSE, TRUE), 100, replace = TRUE),
                           arm = sample(c("placebo", "active"), 100, replace = TRUE))

As an example we have data set that looks like

head(example_data)

To be able to use this with the 'tslrt' package we have to use the "correct" variable names which can be done by

arms <- ifelse(example_data$arm == "placebo", 1,2)

ex_dt <- data.table(time = example_data$OS_time,
                    event = example_data$event_indicator,
                    group = ifelse(example_data$arm == "placebo", "control", "experimental"))

This is only a simple example and with a real data set there might of course be some additional data preparation to get the data in the correct format.

We also plot the data to show how it looks like

# Plot data
fit <- survfit( Surv(time, event) ~ group, data = ex_dt)
plot(fit, col = c("black", "blue"), mark.time = TRUE, pch = c(1, 2))
legend('topright', legend = c("Control", "Experimental"),
       col = c("black", "blue"), lty = c(1, 1), pch = c(1, 2))

To use the new weighted logrank test we also need a hazard ratio function (this function needs to be a one variable function on time to work as input) for calculating the weights

HR_data <- function(t){HR <- pmin(0.5 + 0.01*t, .99) }

Again this is just a simple example on how to create a hazard ratio function. When actually using the weighted logrank test it is best to create the hazard ratio function based on assumptions on the trial or calculating the hazard ratio from the data.

Now to perform the different types logrank tests, we do as with the simulated examples.

# risk table
risk_table <- get_risk_table(ex_dt)

# calculate weights
logrank_weights <- calculate_weights(risk_table, method = "logrank")
delay_weights <- calculate_weights(risk_table, method = "theta", hr_fun = HR_data)
fh_weights <- calculate_weights(risk_table, method = "fh", rho = 1, gamma = 0)

# data table for plotting the weights
weight_dt <- data.table(t  = c(logrank_weights$t, delay_weights$t, fh_weights$t),
                        w = c(logrank_weights$w, delay_weights$w, fh_weights$w),
                        Test = rep(c('logrank', 'new weights', 'fh(1,0)'),
                                   c(length(logrank_weights$w),length(delay_weights$w),length(fh_weights$w))))


ggplot(weight_dt, aes(x = t, y = w, color = Test)) +
  geom_line(size = 1.2) +
  theme_bw() +
  labs(x = "Time (months)", y = "w", title = "Comparison of weights for different tests") +
  theme(text = element_text(size = 18),
        panel.grid.minor = element_blank())

# Calculate Z-scores

# Logrank
calculate_zs(logrank_weights)

# New weights
calculate_zs(delay_weights)

# Flemming-Harrington
calculate_zs(fh_weights)