Examples
library(OneArm2stage)
library(survival)
library(IPDfromKM)Example 1 (Design a study with restricted follow-up)
Assuming we are designing a study with the 2-year event-free survival (EFS) as the primary endpoint to plan the sample size. The hypotheses are: H0: S0 = 62% versus H1: S1 = 80% at the landmark time point 2 years (x0=2). This is equivalent to a hazard ratio of 0.467 assuming the survival time follows the proportional hazard assumption (\(0.62^{0.467}=0.8\)). For this specific study, the historical accrual rate is assumed to be 5 patients per year and the restricted follow-up time for each patient is set to be 2 years (x=2). Given type I error rate and power set to be 0.05 and 0.8, and given the failure time follows an exponential distribution (“WB” distribution with shape = 1 is essentially an exponential distribution), we can use the function Optimal.rKJ to plan this study as shown below.
Design1 <- Optimal.rKJ(dist="WB", shape=1, S0=0.62, x0=2, hr=0.467, x=2, rate=5,
alpha=0.05, beta=0.2)Design1$Two_stage
# n1 c1 n c t1 MTSL ES PS
# 1 26 0.1349 46 1.6171 5.0946 11.2 34.6342 0.5537Based on the above output, we need to plan for an accrual time of 5.1 years (or equivalently enroll 26 patients assuming rate=5) in the 1st stage before we can conduct the interim analysis. If the decision is “Go”, we will enroll additional 20 patients (\(46 - 26 = 20\)) in the 2nd stage to conduct the final analysis. Critical values for the 1st and 2nd stages are 0.1349 and 1.6171, respectively. Overall, this study will take maximum 11.2 years to be completed and have an expected sample size of 34.6342 under H0. The probability of early stopping under H0 is 55.37%.
Example 2 (Conduct interim/final analysis)
Suppose we have followed the design in Example 1 and collected data from 26 patients for four variables: 1) Entry (when patients enter trial), 2) time (observed time of patients under the trial), 3) status (patient status indicator) and 4) Total (the sum of Entry and time) as shown below. For those who would like to try the example on their own, the data samples used can be found in the folder under /OneArm2stage/extdata.
#> Entry time status Total
#> 16 0.1751473 2.0000000 0 2.1751473
#> 15 0.2037084 1.9788860 1 2.1825944
#> 20 0.2057495 0.7585931 1 0.9643426
#> 2 0.4849195 0.3754064 1 0.8603259
#> 9 0.8156683 0.9108228 1 1.7264911
#> 12 0.9105992 2.0000000 0 2.9105992
#> 22 1.0186641 1.9943413 1 3.0130054
#> 26 1.0585012 2.0000000 0 3.0585012
#> 19 1.1806645 2.0000000 0 3.1806645
#> 8 1.4351793 1.4999062 1 2.9350855
#> 11 1.6707471 2.0000000 0 3.6707471
#> 1 1.7450642 1.8066005 1 3.5516647
#> 13 1.9559141 0.3306302 1 2.2865444
#> 14 2.0315492 2.0000000 0 4.0315492
#> 24 2.0733742 2.0000000 0 4.0733742
#> 25 2.5297900 2.0000000 0 4.5297900
#> 23 2.7490369 2.0000000 0 4.7490369
#> 17 3.0558739 2.0000000 0 5.0558739
#> 5 3.2362173 0.1812412 1 3.4174585
#> 10 3.3746414 1.0227582 1 4.3973996
#> 18 4.0819773 1.0126227 0 5.4623913
#> 7 4.3221594 0.7724406 0 6.3221594
#> 3 4.6492995 0.4453005 0 5.7776514
#> 21 4.7298983 0.3647017 0 6.7298983
#> 4 4.7536854 0.3409146 0 6.7536854
#> 6 4.9215797 0.1730203 0 6.9215797We can contruct a Kaplan-Meier survival curve for this hypothetical dataset (dat1, n=26). 
The survival rates estimated with the Kaplan-Meier method can also be otained with this hypothetical dataset (dat1, n=26).
#> Call: survfit(formula = Surv(time, status) ~ 1, data = dat1, conf.type = "log-log")
#>
#> time n.risk n.event survival std.err lower 95% CI upper 95% CI
#> 0.173 26 0 1.000 0.0000 NA NA
#> 0.181 25 1 0.960 0.0392 0.748 0.994
#> 0.331 24 1 0.920 0.0543 0.716 0.979
#> 0.341 23 0 0.920 0.0543 0.716 0.979
#> 0.365 22 0 0.920 0.0543 0.716 0.979
#> 0.375 21 1 0.876 0.0671 0.663 0.958
#> 0.445 20 0 0.876 0.0671 0.663 0.958
#> 0.759 19 1 0.830 0.0778 0.607 0.933
#> 0.772 18 0 0.830 0.0778 0.607 0.933
#> 0.911 17 1 0.781 0.0872 0.549 0.903
#> 1.013 16 0 0.781 0.0872 0.549 0.903
#> 1.023 15 1 0.729 0.0957 0.490 0.869
#> 1.500 14 1 0.677 0.1020 0.435 0.833
#> 1.807 13 1 0.625 0.1067 0.384 0.794
#> 1.979 12 1 0.573 0.1098 0.335 0.753
#> 1.994 11 1 0.521 0.1115 0.289 0.710
#> 2.000 10 0 0.521 0.1115 0.289 0.710We can use the function LRT (log-rank test) in package OneArm2stage to test H0: S0 = 0.62 against H1: S1 = 0.80 at the landmark time x0=2, using the data shown above (n = 26). The “Go”/“No-go” decision can be made based on the calculated critical value Z from the output below.
LRT(dist = "WB", shape = 1, S0=0.62, x0=2, data=dat1)
#> O E Z
#> 10.0000 8.1190 -0.6601Based on the above output, the test statistic Z = -0.6601, which is smaller than the critical value for the 1st stage c1 = 0.1349 (See Example 1). Therefore, the decision is “No-go”, and we need to stop the trial and declare futility.
Now let’s assume what we have collected is a different sample of 26 patients, and their data are shown below.
#> Entry time status Total
#> 16 0.1751473 2.0000000 0 2.175147
#> 15 0.2037084 2.0000000 0 2.203708
#> 20 0.2057495 1.6243964 1 1.830146
#> 2 0.4849195 0.8038681 1 1.288788
#> 9 0.8156683 1.9503701 1 2.766038
#> 12 0.9105992 2.0000000 0 2.910599
#> 22 1.0186641 2.0000000 0 3.018664
#> 26 1.0585012 2.0000000 0 3.058501
#> 19 1.1806645 2.0000000 0 3.180664
#> 8 1.4351793 2.0000000 0 3.435179
#> 11 1.6707471 2.0000000 0 3.670747
#> 1 1.7450642 2.0000000 0 3.745064
#> 13 1.9559141 0.7079877 1 2.663902
#> 14 2.0315492 2.0000000 0 4.031549
#> 24 2.0733742 2.0000000 0 4.073374
#> 25 2.5297900 2.0000000 0 4.529790
#> 23 2.7490369 2.0000000 0 4.749037
#> 17 3.0558739 2.0000000 0 5.055874
#> 5 3.2362173 0.3880967 1 3.624314
#> 10 3.3746414 1.7199586 0 5.374641
#> 18 4.0819773 1.0126227 0 6.081977
#> 7 4.3221594 0.7724406 0 6.322159
#> 3 4.6492995 0.4453005 0 6.649299
#> 21 4.7298983 0.3647017 0 6.729898
#> 4 4.7536854 0.3409146 0 6.753685
#> 6 4.9215797 0.1730203 0 6.921580We can contruct a Kaplan-Meier survival curve for this hypothetical dataset (dat2.1, n=26). 
The survival rates estimated with the Kaplan-Meier method can also be otained with this hypothetical dataset (dat2.1, n=26). We can see at time = 2 years, the survival rate = 0.758.
#> Call: survfit(formula = Surv(time, status) ~ 1, data = dat2.1, conf.type = "log-log")
#>
#> time n.risk n.event survival std.err lower 95% CI upper 95% CI
#> 0.173 26 0 1.000 0.0000 NA NA
#> 0.341 25 0 1.000 0.0000 NA NA
#> 0.365 24 0 1.000 0.0000 NA NA
#> 0.388 23 1 0.957 0.0425 0.729 0.994
#> 0.445 22 0 0.957 0.0425 0.729 0.994
#> 0.708 21 1 0.911 0.0601 0.688 0.977
#> 0.772 20 0 0.911 0.0601 0.688 0.977
#> 0.804 19 1 0.863 0.0736 0.632 0.954
#> 1.013 18 0 0.863 0.0736 0.632 0.954
#> 1.624 17 1 0.812 0.0850 0.572 0.925
#> 1.720 16 0 0.812 0.0850 0.572 0.925
#> 1.950 15 1 0.758 0.0951 0.510 0.892
#> 2.000 14 0 0.758 0.0951 0.510 0.892We can use LRT to test the hypotheses H0: s0 = 0.62 vs. H1: s1 = 0.80 again based on the alternative sample.
LRT(dist = "WB", shape = 1, S0=0.62, x0=2, data=dat2.1)
#> O E Z
#> 5.0000 9.1553 1.3733Since now the test statistic Z becomes 1.3733 and is larger than the 1st stage critical value c1 = 0.1349, the decision is “Go” and the trial will proceed to the 2nd stage. For that we will enroll additional 20 patients (46 - 26 = 20). We assume the final data after adding the additional 20 patients is as shown below.
#> Entry time status Total
#> 28 0.1751473 2.0000000 0 2.175147
#> 27 0.2037084 2.0000000 0 2.203708
#> 34 0.2057495 1.6243964 1 1.830146
#> 3 0.4849195 0.8038681 1 1.288788
#> 12 0.8156683 1.9503701 1 2.766038
#> 18 0.9105992 2.0000000 0 2.910599
#> 39 1.0186641 2.0000000 0 3.018664
#> 46 1.0585012 2.0000000 0 3.058501
#> 33 1.1806645 2.0000000 0 3.180664
#> 11 1.4351793 2.0000000 0 3.435179
#> 16 1.6707471 2.0000000 0 3.670747
#> 1 1.7450642 2.0000000 0 3.745064
#> 20 1.9559141 0.7079877 1 2.663902
#> 26 2.0315492 2.0000000 0 4.031549
#> 42 2.0733742 2.0000000 0 4.073374
#> 43 2.5297900 2.0000000 0 4.529790
#> 41 2.7490369 2.0000000 0 4.749037
#> 30 3.0558739 2.0000000 0 5.055874
#> 6 3.2362173 0.3880967 1 3.624314
#> 14 3.3746414 2.0000000 0 5.374641
#> 31 4.0819773 2.0000000 0 6.081977
#> 10 4.3221594 2.0000000 0 6.322159
#> 4 4.6492995 2.0000000 0 6.649299
#> 36 4.7298983 2.0000000 0 6.729898
#> 5 4.7536854 2.0000000 0 6.753685
#> 9 4.9215797 2.0000000 0 6.921580
#> 17 5.2925843 1.2395647 1 6.532149
#> 25 5.3239487 2.0000000 0 7.323949
#> 21 5.4715527 2.0000000 0 7.471553
#> 37 5.7529107 2.0000000 0 7.752911
#> 38 5.8474833 2.0000000 0 7.847483
#> 13 5.9075769 1.2108378 1 7.118415
#> 35 6.7189713 2.0000000 0 8.718971
#> 32 6.9993933 0.2237612 1 7.223154
#> 7 7.2432056 0.8315737 1 8.074779
#> 15 7.3439092 2.0000000 0 9.343909
#> 19 7.6240584 2.0000000 0 9.624058
#> 45 7.7256301 2.0000000 0 9.725630
#> 2 7.9269020 2.0000000 0 9.926902
#> 40 8.1280655 2.0000000 0 10.128066
#> 23 8.1931211 0.7256491 1 8.918770
#> 44 8.4734503 2.0000000 0 10.473450
#> 8 8.5440524 2.0000000 0 10.544052
#> 22 8.5786669 2.0000000 0 10.578667
#> 29 8.7362357 1.7539730 1 10.490209
#> 24 8.8650923 2.0000000 0 10.865092We can contruct a Kaplan-Meier survival curve for this combined dataset (n=46).
The survival rates estimated with the Kaplan-Meier method can also be otained with this hypothetical dataset (dat2.1, n=26). We can see at time = 2 years, the survival rate = 0.761.
#> Call: survfit(formula = Surv(time, status) ~ 1, data = dat2.2, conf.type = "log-log")
#>
#> time n.risk n.event survival std.err lower 95% CI upper 95% CI
#> 0.224 46 1 0.978 0.0215 0.856 0.997
#> 0.388 45 1 0.957 0.0301 0.837 0.989
#> 0.708 44 1 0.935 0.0364 0.811 0.978
#> 0.726 43 1 0.913 0.0415 0.785 0.966
#> 0.804 42 1 0.891 0.0459 0.758 0.953
#> 0.832 41 1 0.870 0.0497 0.732 0.939
#> 1.211 40 1 0.848 0.0530 0.707 0.924
#> 1.240 39 1 0.826 0.0559 0.682 0.909
#> 1.624 38 1 0.804 0.0585 0.658 0.893
#> 1.754 37 1 0.783 0.0608 0.634 0.877
#> 1.950 36 1 0.761 0.0629 0.610 0.860
#> 2.000 35 0 0.761 0.0629 0.610 0.860Let’s use LRT to conduct the final analysis with the data of total 46 patients.
LRT(dist="WB", shape=1, S0=0.62, x0=2, data=dat2.2)
#> O E Z
#> 11.0000 19.4704 1.9196Based on the above output, we can see the test statistic for final stage Z = 1.9196. It is larger than the critical value of final stage c = 1.6171 (See Example 1). Therefore, we can reject the null hypothesis and claim the success of the trial. Besides, the LRT output also provide information about the observed and expected number of events for the final stage (O = 11 and E = 19.4704).
Example 3 (Design a study with unrestricted follow-up)
Suppose we would like to design a study like the one in Example 1 but with unrestricted follow-up of 2 years (tf = 2). We can use the Optimal.KJ function to plan that.
Design2 <- Optimal.KJ(dist="WB", shape=1, S0=0.62, x0=2, hr=0.467, tf=2, rate=5,
alpha=0.05, beta=0.2)Design2$Two_stage
# n1 c1 n c t1 MTSL ES PS
# 1 16 -0.302 26 1.6135 3.0593 7.2 21.9187 0.3813Based on the above output, we need to have an accrual time of 3.1 years (or equivalently enroll 16 patients given accrual rate = 5) during the 1st stage. If the decision based on the result of the interim analysis is “Go”, then we will enroll additional 10 patients (\(26 - 16 = 10\)) during the 2nd stage before conducting the final analysis. Critical values for the 1st and 2nd stages are -0.302 and 1.6135, respectively. Overall, this study will take maximum 7.2 years to complete and with an expected sample size of 21.9187 under H0. The probability of early stopping under H0 is 38.13%.
Example 4 (Simulation)
The OneArm2stage package also includes the function Sim_KJ and Sim_rKJ that can be used to calculate empirical power and type-I error by simulation given the design parameters (e.g., t1, n1, n, c1, c) obtained from the optimal two-stage design with unrestricted or restricted follow-up, respectively.
As examples, we will conduct simulations and find the empirical power of the corresponding design based on results of Optimal.rkj using Sim_rKJ, assuming each of the four possible parametric distributions (‘WB’, ‘LN’, ‘GM’ and ‘LG’) of the baseline hazard function.
WB
Design3 <- Optimal.rKJ(dist="WB", shape=0.5, S0=0.3, x0=1, hr=0.65, x=1, rate=10,
alpha=0.05, beta=0.2)Design3$Two_stage
# n1 c1 n c t1 MTSL ES PS
# 1 38 0.1688 63 1.6306 3.7084 7.3 48.3058 0.567# calculate empirical power
Sim_rKJ(dist="WB", shape=0.5, S0=0.3, S1=0.3^(0.65), x0=1, x=1, rate=10, t1=3.7084,
c1=0.1688, c=1.6306, n1=38, n=63, N=10000, seed=5868)
#> [1] 0.796LN
Design6 <- Optimal.rKJ(dist="LN", shape=0.5, S0=0.3, x0=1, hr=0.65, x=1, rate=10,
alpha=0.05, beta=0.2)Design6$Two_stage
#> n1 c1 n c t1 MTSL ES PS
#> 1 39 0.1097 63 1.6281 3.8408 7.3 49.6291 0.5437# calculate empirical power
Sim_rKJ(dist="LN", shape=0.5, S0=0.3, S1=0.3^(0.65), x0=1, x=1, rate=10, t1=3.8408, c1=0.1097, c=1.6281, n1=39, n=63, N=10000, seed=20198)
#> [1] 0.799LG
Design7 <- Optimal.rKJ(dist="LG", shape=0.5, S0=0.3, x0=1, hr=0.65, x=1, rate=10,
alpha=0.05, beta=0.2)Design7$Two_stage
# n1 c1 n c t1 MTSL ES PS
# 1 38 0.1958 63 1.6306 3.7089 7.3 48.034 0.5776# calculate empirical power
Sim_rKJ(dist="LG", shape=0.5, S0=0.3, S1=0.3^(0.65), x0=1, x=1, rate=10, t1=3.7089, c1=0.1958, c=1.6306, n1=38, n=63, N=10000, seed=20198)
#> [1] 0.799GM
Design8 <- Optimal.rKJ(dist="GM", shape=0.5, S0=0.3, x0=1, hr=0.65, x=1, rate=10,
alpha=0.05, beta=0.2)Design8$Two_stage
# n1 c1 n c t1 MTSL ES PS
# 1 38 0.1511 63 1.6306 3.7079 7.3 48.4841 0.56# calculate empirical power
Sim_rKJ(dist="GM", shape=0.5, S0=0.3, S1=0.3^(0.65), x0=1, x=1, rate=10, t1=3.7079, c1=0.1511, c=1.6306, n1=38, n=63, N=10000, seed=20198)
#> [1] 0.8Example 5 (Fitting a historical data)
Suppose we would like to plan for a trial by assuming the failure time follows the Weibull distribution but we don’t know the parameters of the underlying distribution. In this situation we can use historical data from a previous literature to make an educated guess. The package IPDfromKM can extract individual patient data (IPD) from digitized Kaplan-Meier plots. The following example code shows how to do that using a study (Wang et al, 2018) in a group of manetel cell lymphoma patients. The plot we choose to extract the IPD is the KM curve of the overall survival (Fig.2D in Wang et al, 2018).
The below steps tell you how to extract IPD data yourself, which you are welcome to try.
## First, we take a screenshot of the plot and save it as a .png file. In this example it is named as "wang2018.png".
## Next, create a path directing to the .png file.
## e.g. path <-"./vignettes/wang2018.png"
## Last, type the following code in console, and follow the prompts to extract data points by clicking your mouse on the plot.
## df <- getpoints(path, 0, 24, 0, 100) # 0,24,0,100 are the ranges of x and y-axes in the imageBut for this vignette, we have saved a sample df dataset that can be used directly.
# a sample dataset that we already extracted from Wang et al, 2018.
df<- read.csv(system.file("extdata", "df.csv", package = "OneArm2stage"))
# risk time points
trisk <- c(0,2,4,6,8,10,12,14,16,18,20,22,24)
# number of patients at risk at each risk time point
nrisk.radio <- c(124,120,115,110,107,104,103,95,46,18,11,8,0)
# Preprocess the raw coordinates into an proper format for reco""nstruct IPD
pre_radio <- preprocess(dat=df, trisk=trisk, nrisk=nrisk.radio,totalpts=NULL,maxy=100)
#Reconstruct IPD
est_radio <- getIPD(prep=pre_radio,armID=0,tot.events=NULL)We can perform a survival analysis on the reconstructed IPD, and compare with the results from the original paper. The survival probability estimation at time = 12 is 0.86 with a 95% CI of (0.8, 0.92). This is very close to the result reported in Wang et al., 2018, which is 0.87 with a 95% CI of (0.79, 0.92).
plot_ex5 <- survreport(ipd1=est_radio$IPD,ipd2=NULL,arms=1,interval=2,s=c(0.75,0.5,0.25),showplots=F)
plot_ex5$arm1$survprob
#> Surv SE 0.95LCI 0.95UCI
#> time = 2 0.9677 0.0159 0.9371 0.9993
#> time = 4 0.9433 0.0208 0.9033 0.9850
#> time = 6 0.9184 0.0247 0.8712 0.9682
#> time = 8 0.8934 0.0280 0.8402 0.9499
#> time = 10 0.8681 0.0307 0.8099 0.9305
#> time = 12 0.8597 0.0316 0.8000 0.9239
#> time = 14 0.8254 0.0347 0.7602 0.8962
#> time = 16 0.7820 0.0379 0.7112 0.8599
#> time = 18 0.7650 0.0407 0.6893 0.8491
#> time = 20 0.6400 0.0914 0.4837 0.8467
#> time = 22 0.6400 0.0914 0.4837 0.8467Now, with the reconstructed IPD, we can fit a parametric regression model assuming Weibull distribution (which is the most widely used option). Before that, we need to first shift the IPD data into the proper format to be used with the FitDat function in our package.
ipd <- est_radio$IPD
dat3 <- as.data.frame(cbind(rep(0, nrow(ipd)),ipd$time, ipd$status))
colnames(dat3) <- c("Entry", "Time", "Cens")We use FitDat function to fit the historical data with the Weibull distribution. The function returns the AIC value for the fitted model and the estimated parameters. We can use these parameter estimates to design the trial as we have shown in Example 1 and 3.
modelWB<- FitDat(dat3)
#AIC
modelWB$AIC
#> Weibull
#> 301.7776
modelWB$parameter.estimates
#> $Weibull
#> shape scale
#> 0.1133671 3.9939753Example 6 (Design a study with pre-defined early stopping probability under H0)
The early stopping probability is the probability of ending a trial upon the completion of the interim stage, which happens when the test statistic \(Z_1\) at the interim analysis is less than \(c_1\). In situations where users wish to pre-define the early stopping probability under H0, the argument prStop can be used to achieve this goal. The following example shows how to set the stopping probability under H0 to be 0.65 while designing a two-stage trial with restricted follow-up. For unrestricted designs, it follows the same procedure using Optimal.KJ().
Optimal.rKJ(dist="WB", shape=1, S0=0.62, x0=2, hr=0.467, x=2, rate=5,
alpha=0.05, beta=0.2, prStop=0.65)
# $Two_stage
# n1 c1 n c t1 MTSL ES PS
# 1 30 0.386 46 1.6116 5.8247 11.2 35.0258 0.6503We can also choose the most desirable design by comparing multiple designs with varying stopping probabilities using the following codes:
# early stopping probability under H0
P <- t(seq(0.6, 0.85, l=6))
# a design without pre-defined stopping probability
o1 <- Optimal.rKJ(dist="WB", shape=1, S0=0.62, x0=2, hr=0.467, x=2, rate=5,
alpha=0.05, beta=0.2)
# designs with varying pre-defined stopping probabilities
res1 <- t(apply(P, 2, Optimal.rKJ,
dist="WB", shape=1, S0=0.62, x0=2, hr=0.467, x=2, rate=5, alpha=0.05, beta=0.2))
# compare results
r1 <- rbind(o1$Two_stage,
res1[[1]]$Two_stage, res1[[2]]$Two_stage,
res1[[3]]$Two_stage, res1[[4]]$Two_stage,
res1[[5]]$Two_stage, res1[[6]]$Two_stage)
# first row shows the design without pre-defined stopping probability (prStop=NULL)
# > r1
# n1 c1 n c t1 MTSL ES PS
# 1 26 0.1349 46 1.6171 5.0946 11.2 34.6350 0.5537
# 2 28 0.2544 46 1.6147 5.4413 11.2 34.7161 0.6004
# 3 30 0.3860 46 1.6116 5.8247 11.2 35.0258 0.6503
# 4 31 0.5247 47 1.6013 6.1331 11.4 35.5642 0.7001
# 5 33 0.6760 47 1.5952 6.5867 11.4 36.4433 0.7505
# 6 36 0.8436 46 1.5927 7.1199 11.2 37.6740 0.8005
# 7 39 1.0382 46 1.5768 7.6639 11.2 39.4686 0.8504Based on the above outputs, the sample size, critical value and length of the interim stage (\(n_1\), \(c_1\), and \(t_1\)) increase with the increase of the early stopping probability (\(PS\)). This trend also applies to the expected sample size under H0 (\(ES\)). However, the critical value at the final stage (\(c\)) gradually decreases as the \(PS\) increases. The total sample size (\(n\)) and the maximum total study length (\(MTSL\)) are relatively stable despite the increase of \(PS\). Based on the specific needs for any particular trial, a decision can be made from the above options based on the design parameters.