[R] Adaptive design code

[Cody Hamilton]

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}}