Evaluation of a multi-arm multi-stage Bayesian design with stan

[Rani]

Hello all,
In this article, article (Evaluation of a multi-arm multi-stage Bayesian design for phase II drug selection trials – an example in hemato-oncology), they study the Evaluation of a multi-arm multi-stage Bayesian design, for example, let us we want to compare two new treatments {trt=1,2} to a common reference {trt=0}

> event = sample(0:1, size = 1000, replace =  T)
> trt = as.factor( sample(c(0,1,2), size = 1000, replace = T) )
> mydata = data.frame(event, trt)
> str(mydata)
'data.frame':	1000 obs. of  2 variables:
 $ event: int  1 1 0 0 0 0 0 1 1 1 ...
 $ trt  : Factor w/ 3 levels "0","1","2": 3 3 1 1 2 3 2 2 2 1 ...

It’s not clear for me how they did the computations, but with rstanarm package, I run the model as follows:

> library(rstanarm)
> post <- stan_glm(event~trt-1, 
+                      data = mydata,
+                      family =  binomial("logit"), 
+                      prior = normal(location = 0, scale = 10, autoscale = FALSE),
+                      prior_intercept =  NULL, 
+                      prior_aux = NULL,
+                      QR= FALSE,
+                      mean_PPD = FALSE,
+                      seed=123456,
+                      cores=2,
+                      chains = 3,
+                      refresh=0)
> simdata = as.data.frame(post)
> head(simdata)
        trt0        trt1        trt2
1 -0.2271590 -0.15856774 -0.06239190
2 -0.2255432  0.13231237  0.03647478
3  0.2200402  0.09135786 -0.04482848
4  0.1804914  0.11293285 -0.08651130
5  0.0746649  0.07895283 -0.11162836
6 -0.1859310 -0.13899772 -0.11444129

Then they said, the inefficacy of each drug was assessed by three approaches, as shown in this picture,

which really I do not know how to do this with rsatnarm, could anyone help and thanks in advance.

[martinmodrak]

Sorry, short on time, but those quantities should be derivable directly from the posterior distribution of the coefficients and/or posterior predictions. Maybe @imadmali is not busy and can provide a bit more details?

[jatotterdell]

As stated by Martin these probabilities can be calculated directly from your draws.
In simdata you have the posterior draws representing the log-odds of response under each treatment which can be transformed to draws for the probability of response \pi_k.
Then you just need to calculate the quoted probabilities by taking the mean over these draws for the comparisons of interest, e.g. using your code to generate a simdat, (1) can be calculated as

p <- plogis(as.matrix(simdata))
p_0 <- 0.5
P_1 <- colMeans(p < p_0)
g_1 <- 0.95
P_1
     trt0      trt1      trt2 
0.7290000 0.0450000 0.4776667 

P_1 > g_1
 trt0  trt1  trt2 
FALSE FALSE FALSE 

So that P_1 is the value of P(\pi_k<p_0|y_{ki},n_{ki}) for k=0,1,2, and P_1 > g_1 the comparison with the decision threshold \gamma_1.

Similarly for (2) and (3),

Delta <- 0.15
delta <- 0.10
g_2 <- 0.05
g_3 <- 0.95
p_d <- p[, -1] - p[, 1]
P_2 <- colMeans(p_d > Delta)
P_3 <- colMeans(p_d > delta)
P_2
       trt1        trt2 
0.009666667 0.001333333

P_3
      trt1       trt2 
0.15733333 0.01666667 

P_2 < g_2
trt1 trt2 
TRUE TRUE 

P_3 > g_3
 trt1  trt2 
FALSE FALSE 

[Rani]

many thanks, looks great