Implementation and analysis design of an adaptive-outcome trial in R

• Platform trial

• ProBio
  • Rationale
    
  • Study design and implementation
• Main Statistical aspects
  • Randomization
    
  • Evaluating therapies
    
  • Operating characteristics via simulations
    
  • Bayesian survival analysis

• Tackle different problems with R

Slides available at http://rpubs.com/alecri/useR2019

• Gold standard for evaluating one therapy for one disease (one pre-established statistical hypothesis)
  
• Traditional clinical trials are ineffective for addressing multiple questions
  
• Conventional 2-group clinical trials are simple but inefficient (separate control population, different protocols, too large or not sufficiently powered)
  
• High costs, slow progress, and a high failure rate

• Platform trial may overcome those problems

• Evaluate several targeted therapies for one disease perpetually, and further accept additions or exclusions of new therapies. Focus is on the disease rather than any particular experimental therapy
  
• One master protocol and consent process, common IT infrastructure, logistics, DSMB

An outcome-adaptive randomized multi-arm biomarker driven trial in patients with metastatic castrate resistant prostate cancer (mCRPC)

Background

• Multiplicity of available treatments for mCRPC and new therapies are expected to soon enrich the landscape of available therapies

• Very heterogeneous response rates and increasing costs

• Retrospective analyses have identified prognostic biomarkers. Could they be predictive as well?

• How to optimize the treatment selection and identify the optimal sequencing of available therapies?

• Patients are stratified based on their biomarker subgroup combination, and then randomized to either the control (standard of care) or one of the active arms

• Outcome-adaptive randomization is implemented to assign more patients to more promising (effective) treatments within the biomarker subgroup combinations

• Treatments are constantly (monthly) evaluated within the biomarker signatures

• Highly effective treatments will graduate from the platform trial and enter into a validation trial (fixed randomization 1:1)

• Patients who progress, will be re-genotyped and re-randomized (max 2 randomizations)

• Fixed before enrolling 50 patients in the active arms

• After, randomization probabilities will be updated monthly based on the accumulated data (PFS)

• Proportional to $\pi_{ij}$, the (Bayesian) probability of superiority for a treatment $i$ in the biomarker subgroup combination $j$:

• Treatments are compared within the biomarker signatures of interest using the control group as comparator

• the main outcome is a survival time. We will use Bayesian parametric model to contrast the distributions of PFS

• Monthly, we will decide if continuing enrollment of new patients to a treatment signature combination, or to early stop (graduation, futility, max patients)

• graduation:
  • $n_{is} \ge 20$
    
  • $\pi_{is} \ge .85$
    
  • $\pi_{ij} \ge .65$, with $j \in s$
• stop futility:
  • $n_{is} \ge 20$
    
  • $\pi_{is} \le .15$
    
  • $\pi_{ij} \le .3$, with $j \in s$
• stop max patients:
  • $n_{is} \ge 150$

• Complicated design – complicated assessment of operating characteristics
  
• Statistical simulations
  
• Different types of errors
  • Graduation when there is no effect
    
  • Fail to graduate when there is an effect
    
  • Graduate drug-signature combination where the biomarker subgroup is too large
    
  • Fail to graduate drug-signature combination where the biomarker subgroup is too small
    
• Discovery trial, relatively liberal graduation criteria
  
• Graduated combinations will be validated in a side trial nested within the ProBio platform

• Simulation Scenarios:
  1. no treatment works better in any signature
     
  2. one treatment (Carboplatin) works better (HR $\approx$ 2)
     
  3. one treatment (Carboplatin) works better (HR $\approx$ 2.4) only in the biomarker subgroup combination belonging to the signature DRD+
     
  4. hormonal therapies work better in the signature TP53- & AR-, Carboplatin in the signature DRD+, while chemotherapies work better in the signature TEfus+ with HR ranging from 2 to 3.
• Results:
  1. Assumed PFS times (median)
  2. Error rate and power
     
  3. Average number of enrolled participants
     
  4. Probabilities of superiority
  5. Graduation times

• Parametric survival models: Weibull
  
• $\log(T_i) = \mu + \gamma^\top z_i + \sigma W$, $W$ has the extreme value distribution
  

• $T_i \sim \textrm{Weibull}(\exp(\mu + \gamma^\top z_i), \sigma)$

• with $\alpha = \sigma^{-1}$, the priors are:
  • $\alpha \sim \textrm{exponential}(1)$
    
  • $\gamma_i \sim \textrm{N}(0, 100)$
    

• $P(\mu_{i} > \mu_{0})$ is computed through MCMC

Simulated example dataset

Kaplan-Meier

km <- survfit(Surv(time_m, delta) ~ trt, data = dat_month)
km
Call: survfit(formula = Surv(time_m, delta) ~ trt, data = dat_month)

                  n events median 0.95LCL 0.95UCL
trt=Control      44     26   5.92    4.61      NA
trt=Enzalutamide 12      7   9.71    5.76      NA
trt=Abiraterone   8      5   6.11    2.46      NA
trt=Carboplatin  31     15  17.86   10.47      NA
trt=Cabazitaxel   8      7   3.19    1.28      NA
trt=Docetaxel    13      8   6.97    3.41      NA

ggsurvplot(km)

Frequentist models

Cox model

fit_cox <- coxph(Surv(time_m, delta) ~ trt, data = dat_month)
fit_cox
Call:
coxph(formula = Surv(time_m, delta) ~ trt, data = dat_month)

                coef exp(coef) se(coef)    z     p
trtEnzalutamide -0.3       0.7      0.4 -0.7 0.466
trtAbiraterone   0.2       1.2      0.5  0.4 0.681
trtCarboplatin  -1.1       0.3      0.3 -3.3 0.001
trtCabazitaxel   0.5       1.7      0.4  1.2 0.228
trtDocetaxel    -0.1       0.9      0.4 -0.3 0.768

Likelihood ratio test=17  on 5 df, p=0.004
n= 116, number of events= 68 

Parametric Weibull model

fit_w <- flexsurvreg(Surv(time_m, delta) ~ trt, data = dat_month, dist = "weibull")
fit_w
Call:
flexsurvreg(formula = Surv(time_m, delta) ~ trt, data = dat_month, 
    dist = "weibull")

Estimates: 
                 data mean  est      L95%     U95%     se       exp(est)
shape                 NA     0.9405   0.7670   1.1532   0.0978       NA 
scale                 NA     8.7122   5.7824  13.1266   1.8221       NA 
trtEnzalutamide   0.1034     0.2755  -0.6129   1.1639   0.4533   1.3172 
trtAbiraterone    0.0690    -0.2560  -1.2739   0.7620   0.5194   0.7742 
trtCarboplatin    0.2672     1.1231   0.4398   1.8064   0.3486   3.0743 
trtCabazitaxel    0.0690    -0.6130  -1.5024   0.2765   0.4538   0.5417 
trtDocetaxel      0.1121     0.0913  -0.7514   0.9341   0.4300   1.0956 
                 L95%     U95%   
shape                 NA       NA
scale                 NA       NA
trtEnzalutamide   0.5418   3.2024
trtAbiraterone    0.2797   2.1425
trtCarboplatin    1.5524   6.0882
trtCabazitaxel    0.2226   1.3185
trtDocetaxel      0.4717   2.5450

N = 116,  Events: 68,  Censored: 48
Total time at risk: 831.8
Log-likelihood = -227.9, df = 7
AIC = 469.8

data {
      // number of gamma parameters
      int<lower=0> P; 
      // data for censored subjects
      int<lower=0> N_m;
      matrix[N_m,P] X_m;
      vector[N_m] y_m;
      // data for observed subjects
      int<lower=0> N_o;
      matrix[N_o,P] X_o;
      vector[N_o] y_o;
    }
    
    parameters {
      vector[P] gamma; 
      real<lower=0> alpha;
    }
    
    transformed parameters{
      vector[N_m] gamma_m;
      vector[N_o] gamma_o;
      gamma_m = exp(X_m*gamma);
      gamma_o = exp(X_o*gamma);
    }
    
    model {
      gamma ~ normal(0, 100);
      alpha ~ exponential(1);
      
      // evaluate likelihood for censored and uncensored subjects
      target += weibull_lpdf(y_o | alpha, gamma_o);
      target += weibull_lccdf(y_m | alpha, gamma_m);
    }
    
    // generate posterior quantities of interest
    generated quantities{
      real sigma;
      real lambda;
      vector[P-1] beta;
      vector[P-1] hr;
      vector[P-1] af;
      vector[P] mu;

      sigma = 1/alpha;
      mu[1] = exp(gamma[1])*tgamma(1 + alpha);
      lambda = exp(-gamma[1]*alpha);
      for (n in 1:(P-1)){
        beta[n] = -gamma[n+1]*alpha;
        hr[n] = exp(-gamma[n+1]*alpha);
        af[n] = exp(gamma[n+1]);
        mu[n+1] = exp(gamma[1] + gamma[n+1])*tgamma(1 + alpha);
      }
    }

create_data_list <- function(event, time, formula, data){
  return(
    list(
      N_m = sum(data[event] == 0),
      X_m = model.matrix(formula, data = data[data[, event] == 0, ]),
      y_m = data[data[, event] == 0, time],
      N_o = sum(data[event] == 1),
      X_o = model.matrix(formula, data = data[data[, event] == 1, ]),
      y_o = data[data[, event] == 1, time],
      P = ncol(model.matrix(formula, data = data[data[, event] == 1, ]))
    )
  )
}

d_list_x <- create_data_list(event = "delta", time = "time_m", 
                             formula = ~ trt, data = dat_month)
str(d_list_x)
List of 7
 $ N_m: int 48
 $ X_m: num [1:48, 1:6] 1 1 1 1 1 1 1 1 1 1 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:48] "3" "4" "6" "9" ...
  .. ..$ : chr [1:6] "(Intercept)" "trtEnzalutamide" "trtAbiraterone" "trtCarboplatin" ...
  ..- attr(*, "assign")= int [1:6] 0 1 1 1 1 1
  ..- attr(*, "contrasts")=List of 1
  .. ..$ trt: chr "contr.treatment"
 $ y_m: num [1:48] 20 20 19 19 19 19 18 17 17 17 ...
 $ N_o: int 68
 $ X_o: num [1:68, 1:6] 1 1 1 1 1 1 1 1 1 1 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:68] "1" "2" "5" "7" ...
  .. ..$ : chr [1:6] "(Intercept)" "trtEnzalutamide" "trtAbiraterone" "trtCarboplatin" ...
  ..- attr(*, "assign")= int [1:6] 0 1 1 1 1 1
  ..- attr(*, "contrasts")=List of 1
  .. ..$ trt: chr "contr.treatment"
 $ y_o: num [1:68] 6.73 13.39 4.79 7.25 17.86 ...
 $ P  : int 6

fit_v <- sampling(weibull_mod_x, data = d_list_x, chains = 4, iter = 6000, 
                  warmup = 1000, pars= c("gamma", "sigma", "lambda", "alpha", 
                                         "beta", "hr", "af", "mu"),
                  seed = 6494684)
fit_v
Inference for Stan model: 7ba39673c8164060800454620cc11679.
4 chains, each with iter=6000; warmup=1000; thin=1; 
post-warmup draws per chain=5000, total post-warmup draws=20000.

            mean se_mean    sd    2.5%     25%     50%     75%   97.5%
gamma[1]    2.20    0.00  0.23    1.78    2.04    2.19    2.35    2.68
gamma[2]    0.36    0.00  0.51   -0.59    0.00    0.33    0.68    1.40
gamma[3]   -0.15    0.00  0.59   -1.21   -0.55   -0.18    0.22    1.11
gamma[4]    1.18    0.00  0.38    0.46    0.93    1.17    1.43    1.98
gamma[5]   -0.57    0.00  0.51   -1.53   -0.92   -0.59   -0.24    0.47
gamma[6]    0.15    0.00  0.48   -0.75   -0.17    0.14    0.46    1.14
sigma       1.15    0.00  0.13    0.93    1.06    1.14    1.23    1.43
lambda      0.15    0.00  0.04    0.08    0.12    0.15    0.17    0.24
alpha       0.88    0.00  0.10    0.70    0.81    0.88    0.94    1.07
beta[1]    -0.31    0.00  0.44   -1.21   -0.59   -0.29    0.00    0.51
beta[2]     0.13    0.00  0.51   -0.95   -0.19    0.16    0.48    1.06
beta[3]    -1.03    0.00  0.33   -1.70   -1.25   -1.03   -0.81   -0.40
beta[4]     0.50    0.00  0.44   -0.40    0.21    0.52    0.80    1.32
beta[5]    -0.13    0.00  0.41   -0.98   -0.40   -0.12    0.15    0.63
hr[1]       0.81    0.00  0.35    0.30    0.55    0.75    1.00    1.67
hr[2]       1.29    0.00  0.65    0.39    0.82    1.18    1.62    2.89
hr[3]       0.38    0.00  0.13    0.18    0.29    0.36    0.44    0.67
hr[4]       1.81    0.01  0.79    0.67    1.24    1.68    2.23    3.73
hr[5]       0.95    0.00  0.39    0.37    0.67    0.89    1.16    1.89
af[1]       1.64    0.01  1.13    0.55    1.00    1.40    1.98    4.07
af[2]       1.04    0.01  0.87    0.30    0.58    0.83    1.25    3.02
af[3]       3.52    0.01  1.48    1.59    2.52    3.23    4.17    7.26
af[4]       0.65    0.00  0.39    0.22    0.40    0.55    0.78    1.60
af[5]       1.31    0.01  0.73    0.47    0.84    1.14    1.58    3.11
mu[1]       8.89    0.02  2.13    5.65    7.40    8.60   10.03   13.90
mu[2]      13.85    0.07  8.52    5.52    8.96   12.01   16.48   32.47
mu[3]       8.82    0.06  6.66    2.86    5.12    7.16   10.49   24.12
mu[4]      29.73    0.07 10.10   16.28   22.86   27.72   34.37   55.05
mu[5]       5.44    0.02  2.92    2.14    3.56    4.76    6.48   12.92
mu[6]      11.08    0.04  5.50    4.74    7.53    9.81   13.15   24.95
lp__     -232.71    0.02  2.03 -237.54 -233.80 -232.35 -231.21 -229.84
         n_eff Rhat
gamma[1] 10027    1
gamma[2] 15507    1
gamma[3] 15235    1
gamma[4] 13149    1
gamma[5] 14718    1
gamma[6] 14800    1
sigma    19633    1
lambda   12862    1
alpha    19453    1
beta[1]  16009    1
beta[2]  15745    1
beta[3]  14002    1
beta[4]  14956    1
beta[5]  15282    1
hr[1]    16149    1
hr[2]    17326    1
hr[3]    13578    1
hr[4]    15726    1
hr[5]    15524    1
af[1]    11651    1
af[2]    11508    1
af[3]    12749    1
af[4]    13398    1
af[5]    13363    1
mu[1]     9892    1
mu[2]    13171    1
mu[3]    13238    1
mu[4]    20244    1
mu[5]    16479    1
mu[6]    16883    1
lp__      7767    1

Samples were drawn using NUTS(diag_e) at Thu Jul 11 14:36:45 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Compared to a frequentist model

WeibullReg(Surv(time_m, delta) ~ trt, data = dat_month)
$formula
Surv(time_m, delta) ~ trt

$coef
                Estimate      SE
lambda           0.13056 0.03626
gamma            0.94049 0.09783
trtEnzalutamide -0.25910 0.42587
trtAbiraterone   0.24075 0.48970
trtCarboplatin  -1.05624 0.32996
trtCabazitaxel   0.57649 0.42693
trtDocetaxel    -0.08591 0.40430

$HR
                    HR     LB    UB
trtEnzalutamide 0.7717 0.3349 1.778
trtAbiraterone  1.2722 0.4872 3.322
trtCarboplatin  0.3478 0.1821 0.664
trtCabazitaxel  1.7798 0.7708 4.109
trtDocetaxel    0.9177 0.4155 2.027

$ETR
                   ETR     LB    UB
trtEnzalutamide 1.3172 0.5418 3.202
trtAbiraterone  0.7742 0.2797 2.142
trtCarboplatin  3.0743 1.5524 6.088
trtCabazitaxel  0.5417 0.2226 1.319
trtDocetaxel    1.0956 0.4717 2.545

$summary

Call:
survival::survreg(formula = formula, data = data, dist = "weibull")
                  Value Std. Error     z      p
(Intercept)      2.1647     0.2091 10.35 <2e-16
trtEnzalutamide  0.2755     0.4533  0.61 0.5433
trtAbiraterone  -0.2560     0.5194 -0.49 0.6221
trtCarboplatin   1.1231     0.3486  3.22 0.0013
trtCabazitaxel  -0.6130     0.4538 -1.35 0.1768
trtDocetaxel     0.0913     0.4300  0.21 0.8318
Log(scale)       0.0614     0.1040  0.59 0.5553

Scale= 1.06 

Weibull distribution
Loglik(model)= -227.9   Loglik(intercept only)= -236.7
    Chisq= 17.69 on 5 degrees of freedom, p= 0.0034 
Number of Newton-Raphson Iterations: 5 
n= 116 

traceplot(fit_v, c("gamma", "sigma", "alpha", "lambda",
                   "beta", "hr", "af", "mu"))