• 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\):

\[P_{ij}(t) \propto P_{ij}^2(t-1) + \pi_{ij}^2\]

• 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 .25\), 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\) 0.5)
     
  3. one treatment (Carboplatin) works better (HR \(\approx\) 0.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 0.4 to 0.27.
• 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)
here are whatever error messages were returned
[[1]]
Stan model '7ba39673c8164060800454620cc11679' does not contain samples.

[[2]]
Stan model '7ba39673c8164060800454620cc11679' does not contain samples.

[[3]]
Stan model '7ba39673c8164060800454620cc11679' does not contain samples.

fit_v
Inference for Stan model: 7ba39673c8164060800454620cc11679.
1 chains, each with iter=6000; warmup=1000; thin=1; 
post-warmup draws per chain=5000, total post-warmup draws=5000.

            mean se_mean    sd    2.5%     25%     50%     75%   97.5%
gamma[1]    2.21    0.00  0.24    1.78    2.04    2.20    2.36    2.69
gamma[2]    0.35    0.01  0.51   -0.59    0.00    0.34    0.68    1.38
gamma[3]   -0.15    0.01  0.58   -1.18   -0.55   -0.19    0.22    1.05
gamma[4]    1.18    0.01  0.39    0.45    0.91    1.17    1.44    2.01
gamma[5]   -0.58    0.01  0.51   -1.55   -0.93   -0.61   -0.26    0.48
gamma[6]    0.14    0.01  0.49   -0.77   -0.18    0.12    0.46    1.14
sigma       1.15    0.00  0.13    0.92    1.07    1.14    1.23    1.44
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.87    0.94    1.08
beta[1]    -0.31    0.01  0.44   -1.18   -0.59   -0.30    0.00    0.51
beta[2]     0.13    0.01  0.50   -0.89   -0.19    0.16    0.48    1.03
beta[3]    -1.03    0.01  0.34   -1.72   -1.26   -1.03   -0.81   -0.39
beta[4]     0.51    0.01  0.44   -0.41    0.23    0.53    0.81    1.32
beta[5]    -0.13    0.01  0.42   -0.99   -0.40   -0.11    0.16    0.65
hr[1]       0.81    0.01  0.36    0.31    0.55    0.74    1.00    1.67
hr[2]       1.29    0.01  0.63    0.41    0.83    1.18    1.62    2.80
hr[3]       0.38    0.00  0.13    0.18    0.28    0.36    0.45    0.68
hr[4]       1.82    0.01  0.79    0.67    1.26    1.70    2.25    3.75
hr[5]       0.96    0.01  0.41    0.37    0.67    0.90    1.18    1.92
af[1]       1.63    0.02  1.00    0.56    1.00    1.40    1.98    3.99
af[2]       1.03    0.01  0.80    0.31    0.58    0.83    1.25    2.85
af[3]       3.53    0.03  1.56    1.56    2.50    3.22    4.20    7.43
af[4]       0.64    0.01  0.39    0.21    0.39    0.55    0.77    1.62
af[5]       1.31    0.01  0.73    0.46    0.83    1.13    1.58    3.13
mu[1]       8.93    0.04  2.16    5.64    7.37    8.64   10.12   14.07
mu[2]      13.84    0.12  7.64    5.31    9.01   12.01   16.53   32.76
mu[3]       8.78    0.11  6.48    2.91    5.14    7.21   10.54   23.47
mu[4]      29.83    0.16 10.56   16.10   22.90   27.91   34.36   55.31
mu[5]       5.42    0.05  2.93    2.14    3.53    4.72    6.39   13.09
mu[6]      11.05    0.09  5.60    4.75    7.43    9.77   13.16   24.74
lp__     -232.76    0.05  2.10 -237.84 -233.83 -232.38 -231.22 -229.85
         n_eff Rhat
gamma[1]  2633    1
gamma[2]  4113    1
gamma[3]  4388    1
gamma[4]  3378    1
gamma[5]  3780    1
gamma[6]  3674    1
sigma     4942    1
lambda    3392    1
alpha     4946    1
beta[1]   4249    1
beta[2]   4530    1
beta[3]   3547    1
beta[4]   3873    1
beta[5]   3769    1
hr[1]     4203    1
hr[2]     4984    1
hr[3]     3506    1
hr[4]     4161    1
hr[5]     3842    1
af[1]     3423    1
af[2]     3292    1
af[3]     3032    1
af[4]     3313    1
af[5]     3328    1
mu[1]     2544    1
mu[2]     4200    1
mu[3]     3768    1
mu[4]     4434    1
mu[5]     3733    1
mu[6]     3988    1
lp__      1867    1

Samples were drawn using NUTS(diag_e) at Tue Nov 19 15:28:54 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"))

mu_dat_v <- data.frame(rstan::extract(fit_v, "mu")$mu)
colnames(mu_dat_v) <- levels(dat_month$trt)

mu_dat_v %>% 
  gather(var, mu) %>%
  ggplot(aes(mu, col = var)) +
  geom_line(stat = "density") +
  xlim(c(0, 30)) +
  labs(col = "Treatment", x = "Mean")

Probability of superiority

diff_mat <- apply(mu_dat_v[, -1], 2, function(x) x - mu_dat_v[, 1])
colMeans(diff_mat > 0)
Enzalutamide  Abiraterone  Carboplatin  Cabazitaxel    Docetaxel 
      0.7490       0.3830       0.9994       0.1232       0.6088 

• Design of the study (simulations and operating characteristics)

• Implementation of adaptive randomization

• Interactive webapp for exemplification of ProBio

• Dynamin reports and presentations

• Interactive dashboard for continuous monitoring