Setting up the Scenarios

Setting up the simulations to be run.createDesigncreates atibble with every combination of the parameters.Each line corresponds to one simulation to be run.

The functiongenerate_delayed_effect needs the columns:n_trt: number of patients in the treatment arm,n_ctrl: number of patients in the control arm,delay: delay until onset of treatment effect,hazard_ctrl: hazard under control and before onset oftreatment effect,hazart_trt: hazard under treatment,t_max: maximum time of observation.

An example designtibble with all parameters filled outcan be created withdesing_skeleton_delayed_effect. Use thefunction to output an example function call that you can copy and modifyas needed, or assign the result to a variable to obtain a designtibble with some default parameters.

By default this will create a simulation design skeleton forsimulations of 50 patients in each arm, a constant hazard of 0.2 undercontrol and a hazard of 0.02 under treatment after the effect onsetvarying from 0 to 10.

N_sim<-100Assumptions<-assumptions_delayed_effect()#> expand.grid(#>   delay=m2d(seq(0, 10, by=2)), # delay of 0, 1, ..., 10 months#>   hazard_ctrl=m2r(24),         # median survival control of 24 months#>   hazard_trt=m2r(36),          # median survival treatment of 36 months#>   random_withdrawal=m2r(120)   # median time to random withdrawal 10 years#> )Options<-design_fixed_followup()#> expand.grid(#>   n_trt=150,         # 150 patients in the treatment arm#>   n_ctrl=150,        # 150 patients in the control arm#>   followup=m2d(24),  # followup 2 years#>   recruitment=m2d(6) # recruitment time 6 months#>#> )Design<-merge(Assumptions, Options,by=NULL)knitr::kable(Design)

Defining the ‘Generate’ funcion

Define the data generating ‘generate’ function, that besidesimulating the time to event also applies the different censoringprocesses.

my_generator<-function(condition,fixed_objects=NULL){generate_delayed_effect(condition, fixed_objects)|>recruitment_uniform(condition$recruitment)|>random_censoring_exp(condition$random_withdrawal)|>admin_censoring_time(condition$followup)}

Defining the ‘Summarise’ function

Next, we need to specify a summary function that computes the desiredoperating characteristics for each simulation scenario, and eachanalysis method. In the example below, we use a summary function thatcomputes the power (as we consider scenarios under the alternative inthe example) of the log-rank test and the max-combo test acrosssimulated scenarios. For each scenario the function just averages thenumber of times the computed p-value is below the significance levelobtain the power.

The results object contains the results of all replications of thecorresponding method for each row of theDesign object. Inthis exampleresults$p contains allN_simp-values returned by theanalyse_maxcombo oranalyse_logrank functions respectively. The Summary willcontain columns with the rejection rate and some other summarystatistics for both methods.

alpha<-0.05Summarise<-create_summarise_function(maxcombo =summarise_test(alpha),logrank  =summarise_test(alpha))

Putting it all together

Now we put it all together: in Design we give the scenarios definedbefore. We want to run 100 replications for each scenario. We want togenerate data using thegenerate_delayed_effect functionusing the parameters fromDesign and analyse eachreplication of each scenario with the two functionsanalyse_logrank andanalyse_maxcombo. Theoutput should be summarised with theSummarise functiondefined before and the simulations should be run in parallel.

res<-runSimulation(  Design,replications = N_sim,generate = my_generator,analyse =list(logrank  =analyse_logrank(),maxcombo =analyse_maxcombo()  ),summarise = Summarise,cl = cl,save=FALSE)#>#> Number of parallel clusters in use: 8#>#> Simulation complete. Total execution time: 10.61s#> Warning in system2("quarto", "-V", stdout = TRUE, env = paste0("TMPDIR=", : Ausführung von Kommando '"quarto" TMPDIR=C:/Users/Tobias11/AppData/Local/Temp/RtmpGKzPDq/file530456a56ad2#> -V' ergab Status 1

Finally we select the interesting columns from the output. Since allother parameters are the same for each scenario we just select delay.And we are interested in the rejection rate of the tests.

res|>subset(select=c("delay","maxcombo.rejection_0.05","logrank.rejection_0.05"))|>  knitr::kable()

Options<-design_group_sequential()#> expand.grid(#>   n_trt=200,          # 200 patients in the treatment arm#>   n_ctrl=200,         # 200 patients in the control arm#>   followup=m2d(48),   # maximum followup time 4 years#>   recruitment=m2d(6), # recruitment time 6 months#>   interim_events=150, # interim analysis after 150 events#>   final_events=300    # final analysis after 300 events#> )Design<-merge(Assumptions, Options,by=NULL)knitr::kable(Design)

‘Analyse’ functions with an interim analysis

Theanalyse_group_sequential function allows to combinetwo or more analyse functions to create an analysis functioncorresponding to a group sequential design. The arguments are the timesor events after which the analyses are done, the nominal alpha at eachstage and the analyse functions to be used at each stage.

## O'Brien-Fleming Bounds for GSD with interim analysis at information time 1/2nominal_alpha<- ldbounds::ldBounds(c(0.5,1))$nom.alphaclusterExport(cl,"nominal_alpha")Analyse<-list(logrank_seq  =analyse_group_sequential(followup =c(condition$interim_events, condition$final_events),followup_type =c("event","event"),alpha = nominal_alpha,analyse_functions =analyse_logrank()  ),maxcombo_seq =analyse_group_sequential(followup =c(condition$interim_events, condition$final_events),followup_type =c("event","event"),alpha = nominal_alpha,analyse_functions =analyse_maxcombo()  ))

A ‘Summarise’ function for the more complex scenario

The output of the function created withanalyse_group_sequential contains additional columns.rejected_at_stage includes the stage at which the null wasfirst rejected orInf if the null was not rejected,N_pat andN_evt contain the number of patientsrecruited and the number of events observed before the null was rejectedandfollowup contains the time after study start at whichthe last analysis was done.

The results object also includes the results returned by each stageinresults_stages, but here we only use the overalltest-decision.

Summarise<-create_summarise_function(maxcombo_seq =summarise_group_sequential(),logrank_seq =summarise_group_sequential())

Putting it all together

The call torunSimulation looks almost the same as abovebut now the additional columns we defined in ourSummarisefunctions are included in the result.

res<-runSimulation(  Design,replications = N_sim,generate = my_generator,analyse = Analyse,summarise = Summarise,cl = cl,save=FALSE)#>#> Number of parallel clusters in use: 8#>#> Simulation complete. Total execution time: 21.83s#> Warning in system2("quarto", "-V", stdout = TRUE, env = paste0("TMPDIR=", : Ausführung von Kommando '"quarto" TMPDIR=C:/Users/Tobias11/AppData/Local/Temp/RtmpGKzPDq/file5304116837a0#> -V' ergab Status 1

In the case of a group sequential design we are also interested inthe average running time of the study in terms of patients recruited,number of events and running time of the study.

res|>subset(select=c("delay","maxcombo_seq.rejection","logrank_seq.rejection","maxcombo_seq.n_pat","logrank_seq.n_pat","maxcombo_seq.n_evt","logrank_seq.n_evt","maxcombo_seq.followup","logrank_seq.followup"    ))|>  knitr::kable()

Calculating the true values of the summary statistics

To evaluate the performance of an estimator, we first compute thevalues of some true summary statistics to which the estimates will becompared. The most relevant true summary statistics can be computed by aconvenience function for each scenario. Just pipe the Design data.frameto the function and the values of the statistics are added ascolumns.

Options<-design_fixed_followup()#> expand.grid(#>   n_trt=150,         # 150 patients in the treatment arm#>   n_ctrl=150,        # 150 patients in the control arm#>   followup=m2d(24),  # followup 2 years#>   recruitment=m2d(6) # recruitment time 6 months#>#> )Design<-merge(Assumptions, Options,by=NULL)Design<- Design|>true_summary_statistics_delayed_effect(cutoff_stats =20)knitr::kable(Design)

Defining theSummarise function

In the Summarise function the true value against which the estimatorshould be compared has to be specified. If coverage and average width ofthe confidence intervals should be estimated, the CI bounds should alsobe specified.

The arguments to the functions are left un-evaluated and are laterevaluated in theresults andconditiondatasets respectively. So any expressions using variables from resultscan be used for the estimated value and the CI bounds and expressionsusing variables from condition can be used in the argument for the realvalue.

In this case we want to compare the hazard ratio estimated by the Coxmodel to the geometric average hazard ratio as well as to the hazardratio after onset of treatment, calculated from the two respectivecolumns of the Design data.frame.

Note that one name can be used twice to summarise the output of oneanalysis method two times, like in this case, comparing it to twodifferent summary statistics.

Summarise<-create_summarise_function(coxph=summarise_estimator(hr, gAHR_20, hr_lower, hr_upper,name="gAHR"),coxph=summarise_estimator(hr, hazard_trt/hazard_ctrl, hr_lower, hr_upper,name="HR")  )

Putting it all together

Analyse<-list(coxph =analyse_coxph())res<-runSimulation(  Design,replications = N_sim,generate = my_generator,analyse = Analyse,summarise = Summarise,cl = cl,save=FALSE)#>#> Number of parallel clusters in use: 8#>#> Simulation complete. Total execution time: 0.77s#> Warning in system2("quarto", "-V", stdout = TRUE, env = paste0("TMPDIR=", : Ausführung von Kommando '"quarto" TMPDIR=C:/Users/Tobias11/AppData/Local/Temp/RtmpGKzPDq/file530433727d3d#> -V' ergab Status 1

res|>subset(select=c("delay","coxph.HR.bias","coxph.gAHR.bias","coxph.HR.mse","coxph.gAHR.mse","coxph.HR.coverage","coxph.gAHR.coverage"  ))|>  knitr::kable()