I am working on the Continual Reassessment Method for phase I clinical trials which models the relationship of drug dose and toxicitiy probability, using a simple power function. I must say, I am not a mathematician, but still want to be able to implement a basic CRM algorithm in R.

Assume we have some inital dose-toxicity model $F(d_k,β)=d_k^{exp⁡(β)}$, where $F$ represents the probability of toxicity, $d$ represents the dose level and $\beta$ is the model parameter which is going to be updated based on the likelihood of the data and a prior distribution. After observing $i-1$ patients, the binomial likelihood of the data is $$L_{i-1}(\beta)=\prod ^{i-1} _{j=1}[{F(x_j,\beta)}]^{Y_j}[1-F(x_j,\beta)]^{1-Y_j}$$ with the outcome $Y$ either be 1 or 0, and the prior distribution $G(β)$, usually a normal distribution. Given these, we want to update $\beta$ using the posterior mean $$\hat{\beta}_{i-1}=\frac{\int{\beta L_{i-1}(\beta )dG(\beta)}}{{\int{L_{i-1}(\beta )dG(\beta)}}}$$ which then can be inserted in $F(d_k,β)=d_k^{exp⁡(β)}$ to define the current shape of the model (which we apply for the drug dose for the current $i$th patient).

My problem is that I do not understand how I am supposed to calculate the posterior mean $\hat{\beta}_{i-1}$. In the first place, I am confused by the fact that there is a function $G(β)$ in the differential which is not part of the integral (please correct me if I am wrong).

All the notation is based on the book on CRM by Cheung (2011). Any help is appreciated! Thanks a lot!

  • $\begingroup$ I believe you are confused here in part by a simple matter of notation. The integral with respect to $\mathrm{d}G$ is meant to be read as "with respect to the probability measure generated by the distribution $G$," and you might just as well replace it with $g(\beta)\,\mathrm{d}\beta$, where $g$ is the density corresponding to $G$. $\endgroup$ – David C. Norris Nov 4 '17 at 1:55
  • $\begingroup$ I should think, however, that a sum running from $j=1$ to $i$ would describe the likelihood corresponding to $i$ observed patients (not $i-1$, as you indicate). Are all your $i$'s and $i-1$'s verbatim from Cheung? $\endgroup$ – David C. Norris Nov 4 '17 at 2:00
  • $\begingroup$ Thank you David! That indeed clarifies the formula to me. I still have to figure out how to implement this in R or JAGS, but you definitely helped me with this step. Regarding your second comment, you are absolutely right: I edited the product in the likelihood, now it should match Cheung's notation and display the correct things. $\endgroup$ – LuckyPal Nov 4 '17 at 10:45
  • $\begingroup$ Glad this helped. Given that you now have $i-1$ uniformly everywhere, you might like to effect the simplification of replacing it throughout by $i$. Always a good idea to simplify the 'math' as much as possible, before starting to write code! $\endgroup$ – David C. Norris Nov 4 '17 at 13:40
  • 1
    $\begingroup$ If you do not want to implement it by yourself: Did you take a look at the R-packages CRM, bcrm? $\endgroup$ – Márcio Augusto Diniz Nov 4 '17 at 15:39

In the end, I used OpenBUGS and specified the model the following way:

model {
for (j in 1:J){
tox[j] ~ dbin(theta[j],n[j])
theta[j] <- pow(skeleton[j],beta)
beta ~ dgamma(rate,scale)

Here, as prior distribution is used a Gamma distribution with rate and scale parameter required to be specified by the user. The variable skeleton represents the vector $d_k$ in the question, the variable theta is the probability of toxicity. With OpenBUGS, samples from the posterior distribution are drawn, hence all values of interest can be calculated.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.