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 <-".../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 reconstruct 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.9939753