Using R package to implement Bayesian phase I/II dose-finding design for three outcomes

Dear all,

I am trying to implement the work done by Suyu Liu, “A Bayesian Phase I/II Trial Design for Immunotherapy”, using R, since the code attached with that work takes a lot of time (more than 20 hours, and the code not complete). So that I tried to use trialr package since it used rstan, but this package allowed me to use two outcomes ( toxicity, efficacy ) and the work of Liu used three outcomes (immune response, toxicity, and efficacy).

I tried to use the R package trialr for two outcomes using the utility for sensitivity analysis written in R code below (table 1 below), I want to see if I used the correct utility and to see how to add a third outcome ( immune response ) to the model? (Question 1)

I attached to you the code and the output for 50 iterations and 60 patients ( I cannot do more, it took around 12 hours),

allocALL.txt (12.0 KB)
recommALL.txt (153 Bytes)
Rcode.txt (15.6 KB)
if you know how he produced the results, (Question 2)
I hope you can feed me back.

and thanks in advance.

### My code###
rm(list = ls())

Uti <- array(0,c(2,3,2)) # order: tox, eff, immuno
Uti[,,1] <- matrix(c(0,0,50,10,80,35),nrow=2)
Uti[,,2] <- matrix(c(5,0,70,20,100,45),nrow=2)
N.max= 60 # patients
outcomes <- '1NNN 2NNT 3NNT 4NNN 5NTN'
doses = c(.1,.3,.5,.7,.9)

fit <- stan_efftox(outcomes,
real_doses =doses,
efficacy_hurdle = 0.5, toxicity_hurdle = 0.3,
p_e = 0.1, p_t = 0.1,
eff0 = 0.5, tox1 = 0.65,
eff_star = 0.7, tox_star = 0.25,
alpha_mean = -7.9593, alpha_sd = 3.5487,
beta_mean = 1.5482, beta_sd = 3.5018,
gamma_mean = 0.7367, gamma_sd = 2.5423,
zeta_mean = 3.4181, zeta_sd = 2.4406,
eta_mean = 0, eta_sd = 0.2,
psi_mean = 0, psi_sd = 1,
seed = 123)

ndoses <- length(fit$prob_tox)
plot(1:ndoses, fit$prob_tox, type="b", pch=19, xlab="Dose level", ylab="Probability of toxicity", ylim=c(0,max(fit$prob_tox) + 0.15), col="green")
points(1:ndoses,fit$prob_eff, type="b", pch=18, col="blue")
abline(h=0.3, lwd=2, lty=4, col = "red")
legend(1, 0.4, legend=c("Toxicity", "Effecacy"),
col=c("green", "blue"), lty=1:2, cex=0.8)
1 Like

Tagging @brockk as the author of the trialr package if they have time to comment. If not, I’ll try to dig into the model to see what’s going on.

1 Like

I cannot access the article because paywall but other work by Liu and Yuan uses conjugate beta priors on binary outcomes for the obvious benefit that it greatly simplifies calculation. Is that the case here? If so, MCMC is not needed.


Hello all,
Here is the article,
nihms-1501127.pdf (549.9 KB)

the issue is not if he is using conjugate beta priors or not, the issue is how dose he produced the tables and the figures from the simulated data in the text file I attached (drive link). Yes, he used conjugate beta priors. But if there is another way to use MCMC by any R tool that designed for Bayesian dose-finding to reduce the time-consuming in his approach (20 hours to complete 100 simulated trials if you use a computer cluster, for me 5o iterations took 13 H on my simple laptop) it will be great. I would like if there is someone who can help to use rstan instead of that. I like rstan really.

Thanks for attaching the paper. From a quick glance though, my feeling is that nothing I have written in trialr would help fit this model, and that if you wanted to fit it, you’d be starting from scratch. I’m afraid I do not have the time to implement it myself.


Even for question one? you think trailr will not help? I mean go with two arms clinic trials directly with the best dose, by considering for instance the first stage by N= 20 and 25 and the second stage adding 20/30 patients.

What is CLT?

clinic trials, sorry, I corrected that above