R/sim_adtte.R
In iain-t-bennett-roche/flexsurvPlus: Provides Functions to Perform Survival Analysis for Economic models
Defines functions
sim_adtte
Documented in
sim_adtte
#' A function to simulate Survival data
#'
#' This function simulates survival data with correlated time to progression and
#' overall survival times. Optionally crossover from treatment arms can be simulated.
#'
#' The simulation times are derived from formulas in Austin 2012 adapted to enable correlations between endpoints.
#' Austin, P.C. (2012), Generating survival times to simulate Cox proportional hazards models with timeâvarying covariates. Statist. Med., 31: 3946-3958. https://doi.org/10.1002/sim.5452
#'
#' @param rc_siminfo Should simulation params be included in simulated dataframe (logical). Defaults to FALSE.
#' @param rc_origos Should OS without switching be included in simulated dataframe (logical). Defaults to FALSE.
#' @param id Identifer added to simulated dataframe. Defaults to 1.
#' @param seed Seed used for random number generator. Defaults to 1234.
#' @param rho correlation coefficient between TTP and OS. Defaults to 0.
#' @param pswitch proportion of patients who switch. Defaults to 0.
#' @param proppd proportion of patients with PFS before switch allowed. Defaults to 0.
#' @param beta_1a treatment effect (as log(Hazard Ratio)) for OS pre PFS. Defaults to log(0.7).
#' @param beta_1b treatment effect (as log(Hazard Ratio)) for OS post PFS. Defaults to log(0.7).
#' @param beta_2a treatment effect (as log(Hazard Ratio)) for OS pre PFS (switch). Defaults to log(0.7).
#' @param beta_2b treatment effect (as log(Hazard Ratio)) for OS post PFS (switch). Defaults to log(0.7).
#' @param beta_pd treatment effect on progression (as log(HR)). Defaults to log(0.4).
#' @param arm_n patients per arm. Defaults to 250.
#' @param enroll_start start of enrollment. Defaults to 0.
#' @param enroll_end end of enrollment. Defaults to 1.
#' @param admin_censor end of trial. Defaults to 2.
#' @param os_gamma weibull shape - for OS. Defaults to 1.2.
#' @param os_lambda weibull scale - for OS. Defaults to 0.3.
#' @param ttp_gamma weibull shape - for TTP. Defaults to 1.5.
#' @param ttp_lambda weibull scale - for TTP. Defaults to 2.
#' @keywords Survival, Simulation
#' @export
#' @examples
#' require(survival)
#' require(dplyr)
#'
#' ADTTE <- sim_adtte()
#' survfit(Surv(AVAL, event = CNSR == 0) ~ ARM, data = filter(ADTTE, PARAMCD == "PFS")) %>%
#' plot()
#'
#' survfit(Surv(AVAL, event = CNSR == 0) ~ ARM, data = filter(ADTTE, PARAMCD == "OS")) %>%
#' plot()
sim_adtte <- function(
rc_siminfo = FALSE, # should simulation info be included?
rc_origos = FALSE, # should original OS without crossover be included?
# params modified by sim
id = 1, # identifier
seed = 1234 , # seed
rho = 0, # correlation coefficient between switch time and os
pswitch = 0, # proportion switch
proppd = 0, # proportion of patients with PFS before switch allowed
beta_1a = log(0.7), # treatment effect (as log(Hazard Ratio))
beta_1b = log(0.7), # treatment effect (as log(Hazard Ratio))
beta_2a = log(0.7), # treatment effect (as log(Hazard Ratio))
beta_2b = log(0.7), # treatment effect (as log(Hazard Ratio))
# params not changed across all sims
beta_pd = log(0.4), # treatment effect on prog (as log(HR))
arm_n = 250, # patients per arm
enroll_start = 0, # start of enrollment
enroll_end = 1, # end of enrollment
admin_censor = 2, # end of trial
os_gamma = 1.2, # weibull shape - for OS
os_lambda = 0.3, # weibull scale - for OS
ttp_gamma = 1.5, # weibull shape - for TTP
ttp_lambda = 2 # weibull scale - for TTP
){
set.seed(seed)
# combine so can select based on x
n_switch <- arm_n * pswitch
# u1 is just used to setup then not used
# u2 and u3 are correlated uniform random variables used for OS and switch time
# u4 used for censoring time (entry time)
# u5 used for selecting switchers
u <- matrix(nrow = arm_n*2, ncol = 5,data = stats::runif(arm_n*2*5,0,1))
# setup correlations using normal variables then convert back to uniform
u[,3] <- stats::pnorm(rho* stats::qnorm(u[,2],0,1) + ((1-rho^2)^0.5)*stats::qnorm(u[,1],0,1),0,1)
# generate covariates
x.trt <- rep(c(0,1), each = arm_n)
# generate os without treatment
os.basis.t <- (-(log(u[,2]) / (os_lambda)))^(1/os_gamma)
# generate ttp
ttp.untreated.t <- (-(log(u[,3]) / (ttp_lambda)))^(1/ttp_gamma)
# apply treatment effect
ttp.treated.t <-(-(log(u[,3]) / (ttp_lambda*exp(beta_pd))))^(1/ttp_gamma)
# deive basis ttp
ttp.basis.t <- ifelse(x.trt==1, ttp.treated.t, ttp.untreated.t)
# generate pfs
pfs.basis.t <- pmin(os.basis.t, ttp.basis.t)
# generate entry time
t.entry <- enroll_start + stats::runif(u[,4])*(enroll_end-enroll_start)
# generate follow up time
t.censor <- admin_censor - t.entry
# define a time in trial after which switch happens
t.start.switch <- as.numeric(stats::quantile(t.entry + pfs.basis.t, proppd))
# who can switch - have ttp before os and censor
can.switch <- as.numeric(
ttp.basis.t < os.basis.t &
ttp.basis.t < t.censor &
t.entry+ttp.basis.t >= t.start.switch
)[x.trt==0]
rank <- rep(arm_n,arm_n)
# choose randomly from those who can switch to get the reqd number of switchers
rank[can.switch==1] <- order(u[x.trt==0 & can.switch==1,5])
x.switch <- c(as.numeric(rank<=n_switch), rep(0,arm_n))
# calculate how many actually switch
actual.p.switch <- sum(x.switch) / arm_n
# generate switch time
t.switch <- rep(NA,arm_n*2)
t.switch[x.switch==1] <- ttp.basis.t[x.switch==1]
# get control non switch os
os.orig.t <- rep(NA,arm_n*2)
os.orig.t[x.trt == 0] <- (-(log(u[x.trt==0,2])) / (os_lambda ))^(1/os_gamma)
# get experimental arm os
t0 <- ttp.basis.t
t0v <- t0 ^ os_gamma
c1 <- os_lambda * exp(beta_1a)
c2 <- os_lambda * exp(beta_1b)
v1 <- 1/os_gamma
os.orig.t[x.trt == 1] <- ifelse(-log(u[x.trt==1,2]) < c1 * t0v[x.trt==1],
(-log(u[x.trt==1,2]) / c1) ^ v1,
((-log(u[x.trt==1,2]) - c1*t0v[x.trt==1] + c2*t0v[x.trt==1]) / c2) ^ v1
)
# get switch treatment (start from orig)
os.cross.t <- os.orig.t
t1 <- ttp.basis.t
t2 <- ttp.basis.t + ttp.treated.t
t1v <- t1^os_gamma
t2v <- t2^os_gamma
c3 <- os_lambda*exp(beta_2a)
c4 <- os_lambda*exp(beta_2b)
d1 <- x.switch==1 & -log(u[,2]) < os_lambda*t1v
d3 <- x.switch==1 & -log(u[,2]) >= os_lambda*t1v + c3*(t2-t1)
d2 <- x.switch==1 & -log(u[,2]) >= os_lambda*t1v & !d3
os.cross.t[d1] <- (-log(u[d1,2]) / os_lambda )^v1
os.cross.t[d2] <- ((-log(u[d2,2])-os_lambda*t1v[d2] + c3*t1v[d2]) / c3 )^v1
os.cross.t[d3] <- ((-log(u[d3,2])-os_lambda*t1v[d3] - c3*(t2v[d3] - t1v[d3]) +c4*t2v[d3]) / c4 )^v1
# apply censoring
os.b.t <- pmin(os.basis.t, t.censor)
os.b.c <- as.numeric(os.b.t==t.censor)
os.o.t <- pmin(os.orig.t, t.censor)
os.o.c <- as.numeric(os.o.t==t.censor)
os.t <- pmin(os.cross.t, t.censor)
os.c <- as.numeric(os.t==t.censor)
# store ttp
ttp.t <- pmin(ttp.basis.t, t.censor)
ttp.c <- as.numeric(ttp.t==t.censor)
# derive PFS (min of os and ttp)
pfs.t <- pmin(os.t, ttp.t)
pfs.c <- as.numeric(pfs.t==t.censor)
# create a data frame
patid <- 1:(arm_n*2)
df_siminfo <- data.frame(
sim_seed = seed,
sim_rho = rho,
sim_proppd = proppd,
sim_pswitch = pswitch,
sim_beta_1a = beta_1a,
sim_beta_1b = beta_1b,
sim_beta_2a = beta_2a,
sim_beta_2b = beta_2b,
sim_id = id,
sim_actpswitch = actual.p.switch
)
df_pfs <- data.frame(
USUBJID = patid,
ARMCD = ifelse(x.trt == 1, "B", "A"),
ARM = ifelse(x.trt == 1, "Intervention Arm B", "Reference Arm A"),
PARAMCD = "PFS",
PARAM = "Progression Free Survival",
AVAL = round(pfs.t * 365),
AVALU = "DAYS",
CNSR = pfs.c
)
df_os <- data.frame(
USUBJID = patid,
ARMCD = ifelse(x.trt == 1, "B", "A"),
ARM = ifelse(x.trt == 1, "Intervention Arm B", "Reference Arm A"),
PARAMCD = "OS",
PARAM = "Overall Survival",
AVAL = round(os.t * 365),
AVALU = "DAYS",
CNSR = os.c
)
df_os_exsw <- data.frame(
USUBJID = patid,
ARMCD = ifelse(x.trt == 1, "B", "A"),
ARM = ifelse(x.trt == 1, "Intervention Arm B", "Reference Arm A"),
PARAMCD = "OSORIG",
PARAM = "Overall Survival Without Switching",
AVAL = round(os.o.t * 365),
AVALU = "DAYS",
CNSR = os.o.c
)
rc <- rbind(df_pfs, df_os)
if (rc_origos == TRUE){
rc <- rbind(rc, df_os_exsw)
}
if (rc_siminfo == TRUE){
rc <- cbind(rc, df_siminfo)
}
return(rc)
}
iain-t-bennett-roche/flexsurvPlus documentation built on Aug. 1, 2022, 10:10 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions sim_adtte
Documented in sim_adtte
#' A function to simulate Survival data
#'
#' This function simulates survival data with correlated time to progression and
#' overall survival times. Optionally crossover from treatment arms can be simulated.
#'
#' The simulation times are derived from formulas in Austin 2012 adapted to enable correlations between endpoints.
#' Austin, P.C. (2012), Generating survival times to simulate Cox proportional hazards models with timeâvarying covariates. Statist. Med., 31: 3946-3958. https://doi.org/10.1002/sim.5452
#'
#' @param rc_siminfo Should simulation params be included in simulated dataframe (logical). Defaults to FALSE.
#' @param rc_origos Should OS without switching be included in simulated dataframe (logical). Defaults to FALSE.
#' @param id Identifer added to simulated dataframe. Defaults to 1.
#' @param seed Seed used for random number generator. Defaults to 1234.
#' @param rho correlation coefficient between TTP and OS. Defaults to 0.
#' @param pswitch proportion of patients who switch. Defaults to 0.
#' @param proppd proportion of patients with PFS before switch allowed. Defaults to 0.
#' @param beta_1a treatment effect (as log(Hazard Ratio)) for OS pre PFS. Defaults to log(0.7).
#' @param beta_1b treatment effect (as log(Hazard Ratio)) for OS post PFS. Defaults to log(0.7).
#' @param beta_2a treatment effect (as log(Hazard Ratio)) for OS pre PFS (switch). Defaults to log(0.7).
#' @param beta_2b treatment effect (as log(Hazard Ratio)) for OS post PFS (switch). Defaults to log(0.7).
#' @param beta_pd treatment effect on progression (as log(HR)). Defaults to log(0.4).
#' @param arm_n patients per arm. Defaults to 250.
#' @param enroll_start start of enrollment. Defaults to 0.
#' @param enroll_end end of enrollment. Defaults to 1.
#' @param admin_censor end of trial. Defaults to 2.
#' @param os_gamma weibull shape - for OS. Defaults to 1.2.
#' @param os_lambda weibull scale - for OS. Defaults to 0.3.
#' @param ttp_gamma weibull shape - for TTP. Defaults to 1.5.
#' @param ttp_lambda weibull scale - for TTP. Defaults to 2.
#' @keywords Survival, Simulation
#' @export
#' @examples
#' require(survival)
#' require(dplyr)
#'
#' ADTTE <- sim_adtte()
#' survfit(Surv(AVAL, event = CNSR == 0) ~ ARM, data = filter(ADTTE, PARAMCD == "PFS")) %>%
#' plot()
#'
#' survfit(Surv(AVAL, event = CNSR == 0) ~ ARM, data = filter(ADTTE, PARAMCD == "OS")) %>%
#' plot()
sim_adtte <- function(
rc_siminfo = FALSE, # should simulation info be included?
rc_origos = FALSE, # should original OS without crossover be included?
# params modified by sim
id = 1, # identifier
seed = 1234 , # seed
rho = 0, # correlation coefficient between switch time and os
pswitch = 0, # proportion switch
proppd = 0, # proportion of patients with PFS before switch allowed
beta_1a = log(0.7), # treatment effect (as log(Hazard Ratio))
beta_1b = log(0.7), # treatment effect (as log(Hazard Ratio))
beta_2a = log(0.7), # treatment effect (as log(Hazard Ratio))
beta_2b = log(0.7), # treatment effect (as log(Hazard Ratio))
# params not changed across all sims
beta_pd = log(0.4), # treatment effect on prog (as log(HR))
arm_n = 250, # patients per arm
enroll_start = 0, # start of enrollment
enroll_end = 1, # end of enrollment
admin_censor = 2, # end of trial
os_gamma = 1.2, # weibull shape - for OS
os_lambda = 0.3, # weibull scale - for OS
ttp_gamma = 1.5, # weibull shape - for TTP
ttp_lambda = 2 # weibull scale - for TTP
){
set.seed(seed)
# combine so can select based on x
n_switch <- arm_n * pswitch
# u1 is just used to setup then not used
# u2 and u3 are correlated uniform random variables used for OS and switch time
# u4 used for censoring time (entry time)
# u5 used for selecting switchers
u <- matrix(nrow = arm_n*2, ncol = 5,data = stats::runif(arm_n*2*5,0,1))
# setup correlations using normal variables then convert back to uniform
u[,3] <- stats::pnorm(rho* stats::qnorm(u[,2],0,1) + ((1-rho^2)^0.5)*stats::qnorm(u[,1],0,1),0,1)
# generate covariates
x.trt <- rep(c(0,1), each = arm_n)
# generate os without treatment
os.basis.t <- (-(log(u[,2]) / (os_lambda)))^(1/os_gamma)
# generate ttp
ttp.untreated.t <- (-(log(u[,3]) / (ttp_lambda)))^(1/ttp_gamma)
# apply treatment effect
ttp.treated.t <-(-(log(u[,3]) / (ttp_lambda*exp(beta_pd))))^(1/ttp_gamma)
# deive basis ttp
ttp.basis.t <- ifelse(x.trt==1, ttp.treated.t, ttp.untreated.t)
# generate pfs
pfs.basis.t <- pmin(os.basis.t, ttp.basis.t)
# generate entry time
t.entry <- enroll_start + stats::runif(u[,4])*(enroll_end-enroll_start)
# generate follow up time
t.censor <- admin_censor - t.entry
# define a time in trial after which switch happens
t.start.switch <- as.numeric(stats::quantile(t.entry + pfs.basis.t, proppd))
# who can switch - have ttp before os and censor
can.switch <- as.numeric(
ttp.basis.t < os.basis.t &
ttp.basis.t < t.censor &
t.entry+ttp.basis.t >= t.start.switch
)[x.trt==0]
rank <- rep(arm_n,arm_n)
# choose randomly from those who can switch to get the reqd number of switchers
rank[can.switch==1] <- order(u[x.trt==0 & can.switch==1,5])
x.switch <- c(as.numeric(rank<=n_switch), rep(0,arm_n))
# calculate how many actually switch
actual.p.switch <- sum(x.switch) / arm_n
# generate switch time
t.switch <- rep(NA,arm_n*2)
t.switch[x.switch==1] <- ttp.basis.t[x.switch==1]
# get control non switch os
os.orig.t <- rep(NA,arm_n*2)
os.orig.t[x.trt == 0] <- (-(log(u[x.trt==0,2])) / (os_lambda ))^(1/os_gamma)
# get experimental arm os
t0 <- ttp.basis.t
t0v <- t0 ^ os_gamma
c1 <- os_lambda * exp(beta_1a)
c2 <- os_lambda * exp(beta_1b)
v1 <- 1/os_gamma
os.orig.t[x.trt == 1] <- ifelse(-log(u[x.trt==1,2]) < c1 * t0v[x.trt==1],
(-log(u[x.trt==1,2]) / c1) ^ v1,
((-log(u[x.trt==1,2]) - c1*t0v[x.trt==1] + c2*t0v[x.trt==1]) / c2) ^ v1
)
# get switch treatment (start from orig)
os.cross.t <- os.orig.t
t1 <- ttp.basis.t
t2 <- ttp.basis.t + ttp.treated.t
t1v <- t1^os_gamma
t2v <- t2^os_gamma
c3 <- os_lambda*exp(beta_2a)
c4 <- os_lambda*exp(beta_2b)
d1 <- x.switch==1 & -log(u[,2]) < os_lambda*t1v
d3 <- x.switch==1 & -log(u[,2]) >= os_lambda*t1v + c3*(t2-t1)
d2 <- x.switch==1 & -log(u[,2]) >= os_lambda*t1v & !d3
os.cross.t[d1] <- (-log(u[d1,2]) / os_lambda )^v1
os.cross.t[d2] <- ((-log(u[d2,2])-os_lambda*t1v[d2] + c3*t1v[d2]) / c3 )^v1
os.cross.t[d3] <- ((-log(u[d3,2])-os_lambda*t1v[d3] - c3*(t2v[d3] - t1v[d3]) +c4*t2v[d3]) / c4 )^v1
# apply censoring
os.b.t <- pmin(os.basis.t, t.censor)
os.b.c <- as.numeric(os.b.t==t.censor)
os.o.t <- pmin(os.orig.t, t.censor)
os.o.c <- as.numeric(os.o.t==t.censor)
os.t <- pmin(os.cross.t, t.censor)
os.c <- as.numeric(os.t==t.censor)
# store ttp
ttp.t <- pmin(ttp.basis.t, t.censor)
ttp.c <- as.numeric(ttp.t==t.censor)
# derive PFS (min of os and ttp)
pfs.t <- pmin(os.t, ttp.t)
pfs.c <- as.numeric(pfs.t==t.censor)
# create a data frame
patid <- 1:(arm_n*2)
df_siminfo <- data.frame(
sim_seed = seed,
sim_rho = rho,
sim_proppd = proppd,
sim_pswitch = pswitch,
sim_beta_1a = beta_1a,
sim_beta_1b = beta_1b,
sim_beta_2a = beta_2a,
sim_beta_2b = beta_2b,
sim_id = id,
sim_actpswitch = actual.p.switch
)
df_pfs <- data.frame(
USUBJID = patid,
ARMCD = ifelse(x.trt == 1, "B", "A"),
ARM = ifelse(x.trt == 1, "Intervention Arm B", "Reference Arm A"),
PARAMCD = "PFS",
PARAM = "Progression Free Survival",
AVAL = round(pfs.t * 365),
AVALU = "DAYS",
CNSR = pfs.c
)
df_os <- data.frame(
USUBJID = patid,
ARMCD = ifelse(x.trt == 1, "B", "A"),
ARM = ifelse(x.trt == 1, "Intervention Arm B", "Reference Arm A"),
PARAMCD = "OS",
PARAM = "Overall Survival",
AVAL = round(os.t * 365),
AVALU = "DAYS",
CNSR = os.c
)
df_os_exsw <- data.frame(
USUBJID = patid,
ARMCD = ifelse(x.trt == 1, "B", "A"),
ARM = ifelse(x.trt == 1, "Intervention Arm B", "Reference Arm A"),
PARAMCD = "OSORIG",
PARAM = "Overall Survival Without Switching",
AVAL = round(os.o.t * 365),
AVALU = "DAYS",
CNSR = os.o.c
)
rc <- rbind(df_pfs, df_os)
if (rc_origos == TRUE){
rc <- rbind(rc, df_os_exsw)
}
if (rc_siminfo == TRUE){
rc <- cbind(rc, df_siminfo)
}
return(rc)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.