I have been trying to create code to calculate the power for an adaptive design
with a survival endpoint according to the method of Schafer and Muller
('Modification of the sample size and the schedule of interim analyses in
survival trials based on interim inspections,' Stats in Med, 2001). This
design allows for the sample size to be increased (if necessary) based on an
interim look at the data, while still preserving the overall type I error of
the trial. I have included my code below. The computed power seems to be
lower than I would expect. It seems to me that the code is failing to add
sufficient additional sample size at the interim look, although I can't
understand why. Perhaps someone will see something I don't . . .
Best regards,
-Cody Hamilton
#############################################################################################
# Program: Adaptive - no futility 2
# Author: Cody Hamilton
# Date: 04.04.07
#
# This program computes the average probability of two outcomes for an adaptive
# trial with a survival outcome - fail to show a significant improvement of
# treatment over control, or succeed to show a significant improvement of
treatment
# over control. No quitting for futility is allowed in this program. If the
trial appears
# underpowered at the interim look, the sample size allowed is increased up to
the maximum
# allowable size specified by the user. The program calculates the average
sample size at the
# trial end. The adaptive design is based off of the paper by Schafer and
Muller in Stats in
# Med (2001).
#
# For now, the program assumes only one interim look that must occur BEFORE the
end
# of the planned accrual period. In other words, the program will only allow
you to
# extend the accrual, not to restart it. The program also assumes exponentially
# distributed survival times. The program requires you to specify the maximum
time
# to which you are willing to extend the accrual period and the maximum time
you wish to
# let the trial run if necessary. The program calculates the expected number of
deaths at
# this time and lets the trial run until this kth death is observed. Thus, the
actual trial
# length my be less or more than the specified max time.
#
# PARAMETERS:
# p.0 = 12 month survival for control group
# p.a = 12 month survival for treatment group
# n.0 = number of patients in the control group
# n.a = number of patients in the treatment group
# t.first = time of first (and only) interim look (must be < t.A)
# t.A = length of planned accrual period
# t.end = planned end of trial
# alpha = overall type I error for the trial
# cond.pow = desired conditional power for the trial (power given the observed
# information at the the first interim look)
# t.maxA = maximum allowable accrual time
# t.maxEnd = desired maximum allowable total trial time
# num.sim = number of simulations
#
#############################################################################################
#set parameters for simulation
p.0 <- .4
p.a <- .6
n.0 <- 225
n.a <- 225
t.first <- 18
t.A <- 12
t.end <- 30
alpha <- 0.05
cond.pow <- .8
t.maxA <- t.A + 12
t.maxEnd <- t.maxA + 12
num.sim <- 1000
## create vectors to hold results
# result holds the final result for the trial (success or fail)
# final.size holds the final sample size for the trial
# final.t holds the final stop time for the trial
# final.stat holds the value for T(k) at the end of the trial
# z holds the value for the final test statistic
# theta is theta from Schafer Muller 2001
result<-rep(0,num.sim)
final.size<-rep(0,num.sim)
final.t<-rep(0,num.sim)
final.stat<-rep(0,num.sim)
z<-rep(0,num.sim)
# k1 is the number of expected events at the end of the trial under current
design
# muhat is the median survival in the combined sample under alternative
hypothesis for PARTNER
muhat<- (-12*log(2)/log(1-p.a)-12*log(2)/log(1-p.0))/2
n<- n.0 + n.a
k1<- (n/t.A)*( exp(-0.69*t.end/muhat) + 0.69*t.A/muhat -
exp(-0.69*(t.end-t.A)/muhat) )/(0.69/muhat)
# start the simulation
for (i in 1:num.sim) {
# if t.first >= t.A then print an error
if (t.first >= t.A) {print ('ERROR - t.first must be < t.A!') }
# calculate the true medians based on 12 month survival rates assuming
exponential hazards
med.0<- -12*log(2)/log(1-p.0)
med.a<- -12*log(2)/log(1-p.a)
# generate start times for initial trial (no adaption yet)
# y.0 holds the times to event from trial start time for the control group
# y.a holds the times for the treatment group
start.0<-runif(n.0,0,t.A)
start.a<-runif(n.a,0,t.A)
# Generate survival data for initial trial design (no adaption yet)
y.0<-start.0 + rexp(n.0,(log(2)/med.0))
y.a<-start.a + rexp(n.a,(log(2)/med.a))
### Compute number of events in generated sample at interim look ###
# Compute number of events in control and treatment groups at interim look
eventfirst.0<-sum(y.0<t.first)
eventfirst.a<-sum(y.a<t.first)
### Calculate conditional power at interim look under the original design
# estimate median survival for both groups at first look = number of events
# divided by sum of patient time at interim look (linearized rate)
patTime.0<-
sum(pmin((y.0-start.0)[start.0<t.first],t.first-start.0[start.0<t.first]))
patTime.a<-
sum(pmin((y.a-start.a)[start.a<t.first],t.first-start.a[start.a<t.first]))
#patTime.0<- sum(t.first-start.0[start.0<t.first])
#patTime.a<- sum(t.first-start.a[start.a<t.first])
med0.hat <-log(2)/(eventfirst.0/patTime.0)
meda.hat <-log(2)/(eventfirst.a/patTime.a)
# compute total number of events at interim look (k0)
k0<-eventfirst.0 + eventfirst.a
# calculate T(k0) from page 3744 of Schafer and Muller.
# time.0 is the time on study for all patients experiencing an event before
t.first
# c.j is the number of patients randomized to the control group before
t.first
# that have been on study as long or longer than the to,e fpr the jth event.
e.j is
# the equivalent number in the experimental group.
T.k<-0
t<-max(y.0[y.0<t.first],y.a[y.a<t.first])
time.0<-(y.0-start.0)[y.0<t]
time.a<-(y.a-start.a)[y.a<t]
for (j in 1:length(time.0)) {
c.j<-sum(
pmin((t.first-start.0[start.0<t]),(y.0-start.0)[start.0<t])>=time.0[j] )
e.j<-sum(
pmin((t.first-start.a[start.a<t]),(y.a-start.a)[start.a<t])>=time.0[j] )
T.k <- T.k + 1 - c.j/(c.j+e.j)
}
for (j in 1:length(time.a)) {
c.j<-sum(
pmin((t.first-start.0[start.0<t]),(y.0-start.0)[start.0<t])>=time.a[j] )
e.j<-sum(
pmin((t.first-start.a[start.a<t]),(y.a-start.a)[start.a<t])>=time.a[j] )
T.k <- T.k - c.j/(c.j+e.j)
}
# compute epsilon from Muller and Shafer
tau <- T.k
eps <- pnorm( (qnorm(alpha/2)*sqrt(k1) + tau)/sqrt(k1 - k0) )
# compute expected number of events at the end of the trial (exp.k) based on
info available at the
# interim look using formula (and notation) from p. 3746 of Schafer and
Muller (2001)
# A is the time from the interim look to the end of the scheduled accrual
(t.A)
# t is the time from the interim look to the scheduled end of trial (t.end)
# muhat is median survival in the combined sample
# n is the number of patients yet to be enrolled (under the initial design)
# B is the time from the first enrollment to the interim look
# r is the set of start times - first randomization for patients that are
randomized and still alive at the interim look
# delta.k is the number of additional deaths expected by scheduled end of
trial (based on Schoenfeld and Richter)
# exp.k is the expected number of deaths at t.end under the original trial
design (no extended accrual)
A<-t.A - t.first
t<-t.end - t.first
muhat<-log(2)/((eventfirst.0 + eventfirst.a)/(sum(patTime.0) +
sum(patTime.a)))
n<- n.0 + n.a - (sum(start.0<t.first) + sum(start.a<t.first))
B<-t.first-min(c(start.0,start.a))
r<-c(start.0[(start.0<t.first)&(y.0>t.first)],start.a[(start.a<t.first)&(y.a>t.first)])-min(c(start.0,start.a))
delta.k<-(n/A)*( exp(-0.69*t/muhat) + 0.69*A/muhat - exp(-0.69*(t-A)/muhat)
)/(0.69/muhat) + sum( 1 - exp(-0.69*(t + B - r)/muhat) )
exp.k<-delta.k + k0
# Calculate number of necessary events to get desired conditional power
(formula on p. 3745)
zeps<-qnorm(1-eps)
zpow<-qnorm(cond.pow)
theta <- min(max(0,log(meda.hat/med0.hat)),log(med.a/med.0))
needK<-(2*(zpow + zeps)/theta)^2 + k0
### If conditional power at interim look is acceptable, let trial finish and
compute final test statistic
### after the exp.k death
if(exp.k > needK){
# final death is exp.k (stop trial after the exp.k event)
k <- exp.k
# final sample size equals the original planned size
final.size[i] <- n.0 + n.a
# Determine when exp.k is reached
final.t[i]<-sort(c(y.0,y.a))[exp.k]
# calculate T(k)
# time.0, c.j, e.j are as defined above
T.k<-0
time.0<-(y.0-start.0)[y.0<final.t[i]]
time.a<-(y.a-start.a)[y.a<final.t[i]]
for (j in 1:length(time.0)) {
c.j<-sum(
pmin((final.t[i]-start.0[start.0<final.t[i]]),(y.0-start.0)[start.0<final.t[i]])>=time.0[j]
)
e.j<-sum(
pmin((final.t[i]-start.a[start.a<final.t[i]]),(y.a-start.a)[start.a<final.t[i]])>=time.0[j]
)
T.k <- T.k + 1 - c.j/(c.j+e.j)
}
for (j in 1:length(time.a)) {
c.j<-sum(
pmin((final.t[i]-start.0[start.0<final.t[i]]),(y.0-start.0)[start.0<final.t[i]])>=time.a[j]
)
e.j<-sum(
pmin((final.t[i]-start.a[start.a<final.t[i]]),(y.a-start.a)[start.a<final.t[i]])>=time.a[j]
)
T.k <- T.k - c.j/(c.j+e.j)
}
#T.k<-0
#time.0<-(y.0-start.0)[y.0<final.t[i]]
#time.a<-(y.a-start.a)[y.a<final.t[i]]
#for (j in 1:length(time.0)) {
#c.j<-sum( (y.0-start.0)[start.0<final.t[i]]>=time.0[j] )
#e.j<-sum( (y.a-start.a)[start.a<final.t[i]]>=time.0[j] )
#T.k <- T.k + 1 - c.j/(c.j+e.j)
#}
#for (j in 1:length(time.a)) {
#c.j<-sum( (y.0-start.0)[start.0<final.t[i]]>=time.a[j] )
#e.j<-sum( (y.a-start.a)[start.a<final.t[i]]>=time.a[j] )
#T.k <- T.k - c.j/(c.j+e.j)
#}
# calculate test statistic - note that tau = T(k0)
z[i] <- (T.k - tau)/sqrt(k - k0)
pval <- 1 - pnorm(z[i])
# if p-value is less than epsilon, then trial is a success (otherwise
its a failure)
result[i] <- ifelse(pval<eps,'Success','Fail')
}
### If conditional power at interim look is less than necessary, then see if
additional patients could
### raise the conditional power to an acceptable level. If not, enroll the
max.
if(exp.k < needK){
# Calculate number of expected number of events if max number of
patients is enrolled.
# n.max is the total number of patients under the maximum accrual time
(assuming accrual remains constant)
# n is the number of patients to be enrolled between t.first and end of
trial
# delta.k expected number of additional events
# other variables are defined above
A <- t.maxA - t.first
t <- t.maxEnd - t.first
n.max <- n.0 + n.a + round( (t.maxA - t.A)*(n.0 + n.a)/t.A )
n <- n.max - (sum(start.0<t.first) + sum(start.a<t.first))
delta.k <- (n/A)*( exp(-0.69*t/muhat) + 0.69*A/muhat -
exp(-0.69*(t-A)/muhat) )/(0.69/muhat) + sum( 1 - exp(-0.69*(t + B - r)/muhat) )
# if number of expected events under maximum sample size is less than
the number
# necessary to gain the desired conditional power, then set new accrual
period to max
# accrual period
# if number of expected events under maximum sample size is > the number
# necessary to gain the desired conditional power, then find the number
# of months of extra accrual (j) necessary to get the required number
of events
# and set new.tA = t.A + j
if((delta.k + k0) < needK) {new.tA <- t.maxA}
else {
new.tA<-0
j<-0
while (new.tA==0) {
j<-j+1
A<-t.A + j - t.first
t<-t.end + j - t.first
# n is the extra number of patients to be enrolled if accrual
period is lengthened
# to time t.A + j
# delta.k is the expected number of additional events by time
t.A + j
n<-round(n.0 + n.a + j*(n.0 + n.a)/t.A - (sum(start.0<t.first)
+ sum(start.a<t.first)))
delta.k<-(n/A)*( exp(-0.69*t/muhat) + 0.69*A/muhat -
exp(-0.69*(t-A)/muhat) )/(0.69/muhat) + sum( 1 - exp(-0.69*(t + B - r)/muhat) )
# if adding j to the accrual time provides enough events then
stop and set new.tA = t.A + j (otherwise try a larger j)
if ( (delta.k + k0)>needK ) {new.tA<-t.A + j}
# some times delta.k + k0 doesnt quite reach needK at tmaxA due
to round off error - in this case
# set new.tA to t.maxA
if ( (j>(t.maxA-t.A))&new.tA==0 ) {new.tA<-t.maxA}
}
}
#k is the expected number of deaths at new.tEnd
k<- k0 + delta.k
# create additional set of patients accrued until time new.tA (new end
of accrual)
new.tEnd <- t.end + new.tA - t.A
# number of extras is the extra accrual time times the average rate of
accrual
nextra.0 <- round( (0.5*(n.0+n.a)/t.A)*(new.tA - t.A) )
nextra.a <- round( (0.5*(n.0+n.a)/t.A)*(new.tA - t.A) )
# final sample size = original sample plus the extra sample
final.size[i] <- nextra.0 + nextra.a + n.0 + n.a
#generate start times for additional patients
startextra.0<-runif(nextra.0,t.A,new.tA)
startextra.a<-runif(nextra.a,t.A,new.tA)
#generate additional sample
yextra.0<- startextra.0 + rexp(nextra.0,(log(2)/med.0))
yextra.a<- startextra.a + rexp(nextra.a,(log(2)/med.a))
# Determine when k is reached (its the kth element of the sorted list
of all event times)
final.t[i]<-sort(c(y.0,y.a,yextra.0,yextra.a))[k]
#combine original and extra data to compute T(k)
yall.0<-c(y.0,yextra.0)
yall.a<-c(y.a,yextra.a)
startall.0<-c(start.0,startextra.0)
startall.a<-c(start.a,startextra.a)
# calculate T(k)
# time.0 is the time on study for all control events before final.t[j]
(time.a is
# the equivalent number in the treatment group)
# c.j is the number of patients randomized to the control group before
final.t[i]
# that have been on study longer than the time on study for the jth
death. e.j
# is the equivalent number in the experimental group
T.k<-0
time.0<-(yall.0-startall.0)[yall.0<final.t[i]]
time.a<-(yall.a-startall.a)[yall.a<final.t[i]]
for (j in 1:length(time.0)) {
c.j<-sum(
pmin((final.t[i]-startall.0[startall.0<final.t[i]]),(yall.0-startall.0)[startall.0<final.t[i]])>=time.0[j]
)
e.j<-sum(
pmin((final.t[i]-startall.a[startall.a<final.t[i]]),(yall.a-startall.a)[startall.a<final.t[i]])>=time.0[j]
)
T.k <- T.k + 1 - c.j/(c.j+e.j)
}
for (j in 1:length(time.a)) {
c.j<-sum(
pmin((final.t[i]-startall.0[startall.0<final.t[i]]),(yall.0-startall.0)[startall.0<final.t[i]])>=time.a[j]
)
e.j<-sum(
pmin((final.t[i]-startall.a[startall.a<final.t[i]]),(yall.a-startall.a)[startall.a<final.t[i]])>=time.a[j]
)
T.k <- T.k - c.j/(c.j+e.j)
}
# calculate test statistic - note that tau = T(k0)
z[i] <- (T.k - tau)/sqrt(k - k0)
pval <- 1 - pnorm(z[i])
# if p-value is less than epsilon, then trial is a success (otherwise
its a failure)
result[i] <- ifelse(pval<eps,'Success','Fail')
}
final.stat[i] <- T.k
}
print('% for each outcome')
table(result)/num.sim
print('Average total sample size')
mean(final.size)
This message contains information which may be confident...{{dropped:8}}
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
