# Calculating the posterior mean of the model parameter in the Continual Reassessement Method

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!

• 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$. – David C. Norris Nov 4 '17 at 1:55
• 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? – David C. Norris Nov 4 '17 at 2:00
• 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. – LuckyPal Nov 4 '17 at 10:45
• 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! – David C. Norris Nov 4 '17 at 13:40
• If you do not want to implement it by yourself: Did you take a look at the R-packages CRM, bcrm? – Márcio Augusto Diniz Nov 4 '17 at 15:39

## 1 Answer

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.