Simulates data for studies involving treatment switching, incorporating time-dependent confounding. The generated data can be used to evaluate methods for handling treatment switching in survival analysis.
tssim(
tdxo = 0L,
coxo = 1L,
p_R = 0.5,
p_X_1 = NA_real_,
p_X_0 = NA_real_,
rate_T = NA_real_,
beta1 = NA_real_,
beta2 = NA_real_,
gamma0 = NA_real_,
gamma1 = NA_real_,
gamma2 = NA_real_,
gamma3 = NA_real_,
gamma4 = NA_real_,
zeta0 = NA_real_,
zeta1 = NA_real_,
zeta2 = NA_real_,
zeta3 = NA_real_,
alpha0 = NA_real_,
alpha1 = NA_real_,
alpha2 = NA_real_,
theta1_1 = NA_real_,
theta1_0 = NA_real_,
theta2 = NA_real_,
rate_C = NA_real_,
followup = NA_integer_,
days = NA_integer_,
n = NA_integer_,
NSim = 1000L,
seed = NA_integer_
)A list of data frames, each containing simulated longitudinal and event history data with the following variables:
id: Subject identifier.
trtrand: Randomized treatment assignment (0 = control,
1 = experimental)
bprog: Baseline prognosis (0 = good, 1 = poor).
tpoint: Treatment cycle index.
tstart: Start day of the treatment cycle.
tstop: End day of the treatment cycle.
L: Time-dependent covariate predicting survival and switching;
affected by treatment switching.
Llag: Lagged value of L.
Z: Disease progression status at tstop.
A: Treatment switching status at tstop.
Alag: Lagged value of A.
Y: Death indicator at tstop.
timeOS: Observed time to death or censoring.
died: Indicator of death by end of follow-up.
progressed: Indicator of disease progression by end of follow-up.
timePD: Observed time to progression or censoring.
xo: Indicator for whether treatment switching occurred.
xotime: Time of treatment switching (if applicable).
censor_time: Administrative censoring time.
Logical indicator for timing of treatment switching:
1: Treatment switching can occur at or after disease progression.
0: Treatment switching is restricted to the time of disease progression.
Logical indicator for arm-specific treatment switching:
1: Treatment switching occurs only in the control arm.
0: Treatment switching is allowed in both arms.
Probability of randomization to the experimental arm.
Probability of poor baseline prognosis in the experimental arm.
Probability of poor baseline prognosis in the control arm.
Baseline hazard rate for time to death.
Log hazard ratio for randomized treatment (R).
Log hazard ratio for baseline covariate (X).
Intercept for the time-dependent covariate model (L).
Coefficient for lagged treatment switching (Alag)
in the L model.
Coefficient for lagged L in the L model.
Coefficient for baseline covariate (X) in the
L model.
Coefficient for randomized treatment (R) in the
L model.
Intercept for the disease progression model (Z).
Coefficient for L in the Z model.
Coefficient for baseline covariate (X) in the
Z model.
Coefficient for randomized treatment (R) in the
Z model.
Intercept for the treatment switching model (A).
Coefficient for L in the A model.
Coefficient for baseline covariate (X) in the
A model.
Negative log time ratio for A (experimental arm).
Negative log time ratio for A (control arm).
Negative log time ratio for L.
Hazard rate for random (dropout) censoring.
Number of treatment cycles per subject.
Number of days in each treatment cycle.
Number of subjects per simulation.
Number of simulated datasets.
Random seed for reproducibility.
Kaifeng Lu, kaifenglu@gmail.com
The simulation algorithm is adapted from Xu et al. (2022), by
simulating disease progression status and by using the multiplicative
effects of baseline and time-dependent covariates on survival time.
The design options tdxo and coxo indicate
the timing of treatment switching and the study arm eligibility
for switching, respectively. Each subject undergoes followup
treatment cycles until administrative censoring.
At randomization, each subject is assigned treatment based on: $$R_i \sim \mbox{Bernoulli}(p_R)$$ and a baseline covariate is generated: $$X_i \sim \mbox{Bernoulli}(p_{X_1} R_i + p_{X_0} (1-R_i))$$
The initial survival time is drawn from an exponential distribution with hazard: $$rate_T \exp(\beta_1 R_i + \beta_2 X_i)$$ We define the event indicator at cycle \(j\) as $$Y_{i,j} = I(T_i \leq j\times days)$$
The initial states are set to \(L_{i,0} = 0\), \(Z_{i,0} = 0\), \(Z_{i,0} = 0\), \(Y_{i,0} = 0\). For each treatment cycle \(j=1,\ldots,J\), we set \(tstart = (j-1) \times days\).
Generate time-dependent covariates: $$\mbox{logit} P(L_{i,j}=1|\mbox{history}) = \gamma_0 + \gamma_1 A_{i,j-1} + \gamma_2 L_{i,j-1} + \gamma_3 X_i + \gamma_4 R_i$$
If \(T_i \leq j \times days\), set \(tstop = T_i\) and \(Y_{i,j} = 1\), which completes data generation for subject \(i\).
If \(T_i > j \times days\), set \(tstop = j\times days\), \(Y_{i,j} = 0\), and perform the following before proceeding to the next cycle for the subject.
Generate disease progression status: If \(Z_{i,j-1} = 0\), $$\mbox{logit} P(Z_{i,j}=1 | \mbox{history}) = \zeta_0 + \zeta_1 L_{i,j} + \zeta_2 X_i + \zeta_3 R_i$$ Otherwise, set \(Z_{i,j} = 1\).
Generate alternative therapy status: If \(A_{i,j-1} = 0\), determine switching eligibility based on design options. If switching is allowed: $$\mbox{logit} P(A_{i,j} = 1 | \mbox{history}) = \alpha_0 + \alpha_1 L_{i,j} + \alpha_2 X_i$$ If switching is now allowed, set \(A_{i,j} = 0\). If \(A_{i,j-1} = 1\), set \(A_{i,j} = 1\) (once switched to alternative therapy, remain on alternative therapy).
Update survival time based on changes in alternative therapy status and time-dependent covariates: $$T_i = j\times days + (T_i - j\times days) \exp\{ -(\theta_{1,1}R_i + \theta_{1,0}(1-R_i))(A_{i,j} - A_{i,j-1}) -\theta_2 (L_{i,j} - L_{i,j-1})\}$$
Additional random censoring times are generated from an exponential distribution with hazard rate \(rate_C\).
To estimate the true treatment effect in a Cox marginal structural model, one can set \(\alpha_0 = -\infty\), which effectively forces \(A_{i,j} = 0\) (disabling treatment switching). The coefficient for the randomized treatment can then be estimated using a Cox proportional hazards model.
Jessica G. Young, and Eric J. Tchetgen Tchetgen. Simulation from a known Cox MSM using standard parametric models for the g-formula. Statistics in Medicine. 2014;33(6):1001-1014.
NR Latimer, IR White, K Tilling, and U Siebert. Improved two-stage estimation to adjust for treatment switching in randomised trials: g-estimation to address time-dependent confounding. Statistical Methods in Medical Research. 2020;29(10):2900-2918.
Jing Xu, Guohui Liu, and Bingxia Wang. Bias and type I error control in correcting treatment effect for treatment switching using marginal structural models in Phse III oncology trials. Journal of Biopharmaceutical Statistics. 2022;32(6):897-914.
simulated.data <- tssim(
tdxo = 0, coxo = 0, p_R = 0.5, p_X_1 = 0.3, p_X_0 = 0.3,
rate_T = 0.002, beta1 = -0.5, beta2 = 0.3,
gamma0 = 0.3, gamma1 = -0.9, gamma2 = 0.7, gamma3 = 1.1, gamma4 = -0.8,
zeta0 = -3.5, zeta1 = 0.5, zeta2 = 0.2, zeta3 = -0.4,
alpha0 = 0.5, alpha1 = 0.5, alpha2 = 0.4,
theta1_1 = -0.4, theta1_0 = -0.4, theta2 = 0.2,
rate_C = 0.0000855, followup = 20, days = 30,
n = 500, NSim = 100, seed = 314159)
Run the code above in your browser using DataLab