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

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)
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.

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?

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

2 Likes

many thanks, looks great