Bayesian Model Averaging Across CRM Skeletons

Why Average Across CRM Skeletons?

A CRM skeleton encodes a proposed pattern of toxicity probabilities across the dose levels. In practice, more than one pattern may be plausible, and the data may not identify one skeleton decisively. Bayesian model averaging (BMA) represents this uncertainty by fitting each candidate skeleton and combining its posterior toxicity probabilities, rather than committing to one curve alone.

This tutorial uses a power-model CRM with a finite set of candidate skeletons. The same averaging principle applies to other well-defined CRM models, but the likelihood and parameter updating must match the model actually used. Dose selection below means choosing the available dose whose model-averaged posterior mean toxicity probability is closest to the target, subject to any protocol-specified eligibility rules.

Definition: Let \(M_m\) denote candidate skeleton model \(m\), with prior model probability \(\pi_m=P(M_m)\). For dose level \(i\), let \(a_{mi}\) be its skeleton probability and let \(q_{mi}(\gamma)\) be its toxicity probability at CRM parameter value \(\gamma\). The model-averaged posterior toxicity probability is \( \bar q_i= \sum_m P(M_m\mid D)\,E\{q_{mi}(\gamma)\mid D,M_m\}\), where \(D\) is the observed trial data. Thus, both the parameter and the skeleton-model uncertainty are represented in the average.

Model, Likelihood, and Assumptions

In the power model, the skeleton values are between zero and one and are ordered by dose. The transformation \(q_{mi}(\gamma)=a_{mi}^{\gamma}\), with \(\gamma>0\), preserves that ordering. A value of \(\gamma\) below one shifts probabilities upward, while a value above one shifts them downward. The skeletons and the prior distribution for \(\gamma\) are specified before updating with the trial outcomes.

q_{mi}(\gamma)=a_{mi}^{\gamma},\qquad 0<a_{m1}<a_{m2}<\cdots<a_{mK}<1

Suppose \(n_i\) participants have evaluable toxicity outcomes at dose \(i\), of whom \(y_i\) experience the toxicity. Under a binomial likelihood, conditional on the assigned doses and \(\gamma,M_m\), the observed outcomes are independent and each participant at dose \(i\) has toxicity probability \(q_{mi}(\gamma)\). The likelihood, omitting binomial coefficients that do not depend on \(m\) or \(\gamma\), is

L_m(\gamma;D)\ \propto\ \prod_{i=1}^{K}q_{mi}(\gamma)^{y_i}\{1-q_{mi}(\gamma)\}^{n_i-y_i}.

This likelihood presumes outcomes are observed and correctly classified, and that the model represents the relevant dose-toxicity relationship. If outcomes are delayed or censored, the likelihood must account for that rather than treating unevaluated participants as non-toxicities. Given a prior density or mass function \(g_m(\gamma)\), the within-model posterior is proportional to \(L_m(\gamma;D)g_m(\gamma)\). The marginal likelihood \(m_m(D)\) averages the likelihood over that prior. Posterior model probabilities then follow from Bayes' rule:

m_m(D)=\int L_m(\gamma;D)g_m(\gamma)\,d\gamma,\qquad P(M_m\mid D)=\frac{\pi_m m_m(D)}{\sum_{\ell}\pi_\ell m_\ell(D)}.

For a discrete prior, replace the integral with a sum over the supported values of \(\gamma\). The marginal likelihood is the link between the data and the model weights: models that assigned greater prior-predictive probability to the observed data receive greater posterior weight. The weights are nonnegative and sum to one.

Why the Averaging Formula Holds

Let \(T_i\) denote the event that a new participant treated at dose \(i\) experiences toxicity. Conditional on model \(M_m\), the posterior predictive probability of \(T_i\) is the posterior mean of \(q_{mi}(\gamma)\). Since the candidate models form a partition of the model uncertainty, the law of total probability gives

P(T_i\mid D)=\sum_m P(T_i\mid D,M_m)P(M_m\mid D) =\sum_m E\{q_{mi}(\gamma)\mid D,M_m\}P(M_m\mid D).

This establishes the model-averaged quantity as a posterior predictive toxicity probability, not just a convenient average of curves. It also explains why the weights are posterior model probabilities, not automatically the prior model probabilities. Within each model, uncertainty in \(\gamma\) is integrated first; uncertainty about which model applies is then integrated using the posterior model weights.

Worked Example: Updating the Parameter Within Skeleton A

Consider a hypothetical trial with one participant treated at dose 1 and no observed toxicity. Candidate skeleton A is \((0.10,0.30,0.60)\), and the prior puts probability \(0.5\) on each of \(\gamma=0.5\) and \(\gamma=1.5\). At dose 1, the likelihood of this outcome is \(1-q_{A1}(\gamma)\).

At \(\gamma=0.5\), \(q_{A1}=0.10^{0.5}=0.3162\), so the likelihood is \(1-0.3162=0.6838\). At \(\gamma=1.5\), \(q_{A1}=0.10^{1.5}=0.03162\), so the likelihood is \(1-0.03162=0.9684\). The unnormalized posterior masses are \(0.5(0.6838)\) and \(0.5(0.9684)\); the common factor \(0.5\) cancels when normalizing. Therefore,

P(\gamma=0.5\mid D,M_A) =\frac{0.6838}{0.6838+0.9684}\approx0.4139,\qquad P(\gamma=1.5\mid D,M_A)\approx0.5861.

The likelihood and posterior calculation are based on the stated single non-toxicity at dose 1. The posterior shifts toward \(\gamma=1.5\), because that value assigned a higher probability to observing no toxicity at the starting dose.

Compute Each Model's Posterior Toxicity Curve

For each dose and each candidate model, calculate the toxicity probability at every supported \(\gamma\), multiply by its posterior probability, and sum. The result is a posterior mean toxicity probability under that model. It is not generally equal to the toxicity probability evaluated at the posterior mean of \(\gamma\), because \(a_{mi}^{\gamma}\) is nonlinear in \(\gamma\).

Worked Example: The Posterior Curve Under Skeleton A

Use skeleton A and the posterior parameter probabilities just calculated. At dose 1, the two probabilities are \(0.10^{0.5}=0.3162\) and \(0.10^{1.5}=0.03162\). Weighting them by \(0.4139\) and \(0.5861\) gives

\bar q_{A1}=0.4139(0.3162)+0.5861(0.03162)\approx0.1494.

At dose 2, the corresponding probabilities are \(0.30^{0.5}=0.5477\) and \(0.30^{1.5}=0.1643\), giving \(0.4139(0.5477)+0.5861(0.1643)\approx0.3230\). At dose 3, they are \(0.60^{0.5}=0.7746\) and \(0.60^{1.5}=0.4648\), giving \(0.4139(0.7746)+0.5861(0.4648)\approx0.5930\). Thus the within-model posterior mean curve is approximately \((0.1494,0.3230,0.5930)\). The values are probabilities of toxicity, not proportions observed in the trial.

A different skeleton can yield a different parameter posterior, because its likelihood for the same outcomes differs. For the single non-toxicity in this example, consider skeleton B \((0.15,0.40,0.70)\), again with equal prior masses on \(\gamma=0.5\) and \(1.5\). The likelihoods are \(1-\sqrt{0.15}=0.6127\) and \(1-0.15^{1.5}=0.9419\). Normalizing gives posterior parameter weights approximately \(0.3941\) and \(0.6059\). Applying those weights to the transformed skeleton values gives the model-specific posterior curve below.

\begin{aligned} \bar q_{B1}&=0.3941(0.3873)+0.6059(0.05809)\approx0.1878,\\ \bar q_{B2}&=0.3941(0.6325)+0.6059(0.2530)\approx0.4025,\\ \bar q_{B3}&=0.3941(0.8367)+0.6059(0.5857)\approx0.6846. \end{aligned}

Worked Example: Average the Two Candidate Curves

For these hypothetical data, the marginal likelihood under A is the prior-weighted likelihood average \(\{0.6838+0.9684\}/2=0.8261\). Under B it is \(\{0.6127+0.9419\}/2\approx0.7773\). Suppose the prior model probabilities are equal. The posterior probability of A is therefore \(0.8261/(0.8261+0.7773)\approx0.5152\), and that of B is approximately \(0.4848\).

For each dose, average the two within-model posterior means using these posterior model weights. Using unrounded values in the calculation and reporting the result to four significant figures gives

\begin{aligned} \bar q_1&\approx0.5152(0.1494)+0.4848(0.1878)=0.1680,\\ \bar q_2&\approx0.5152(0.3230)+0.4848(0.4025)=0.3616,\\ \bar q_3&\approx0.5152(0.5930)+0.4848(0.6846)=0.6374. \end{aligned}

The model-averaged curve is approximately \((0.1680,0.3616,0.6374)\). If the target toxicity probability is \(0.30\), the absolute differences are \(0.1320\), \(0.0616\), and \(0.3374\). Dose 2 is closest to target, so it is selected under the stated target-proximity rule, provided it is eligible under the trial's separate safety and protocol rules.

Interpretation, Checks, and Alternatives

The model-averaged curve is a weighted average of posterior predictive probabilities. It need not equal any one candidate curve. Nor is it generally obtained by averaging the skeleton entries first and then applying the CRM transformation: in general, transforming an average is not the same as averaging transformed values. In the example, the weights are applied to posterior toxicity probabilities at each dose, after the parameter has been updated within each skeleton.

Several checks help identify implementation or modeling problems:

Each candidate skeleton should be ordered and have entries strictly between zero and one; the parameter transformation should preserve dose ordering.
Within every model, posterior parameter weights should be nonnegative and sum to one. Posterior model weights should also sum to one, and each averaged probability should lie between zero and one.
With monotone candidate curves and the same nonnegative model weights at every dose, the model-averaged curve should be nondecreasing. A violation suggests a coding or indexing error, or a model specification that does not preserve monotonicity.
Compare observed toxicity counts with fitted probabilities, especially at doses with substantial data. A visibly poor fit across all candidate curves indicates that averaging them does not repair the misspecification.
Check sensitivity to the candidate skeletons and prior model probabilities. If the dose recommendation changes materially under reasonable alternatives, report that instability rather than presenting the averaged choice as definitive.

A single-skeleton CRM is simpler and may be appropriate when one skeleton is strongly justified, but it does not represent uncertainty across skeleton choices. Selecting one skeleton after inspection is also possible, but can make the final analysis depend heavily on a discrete choice. BMA retains multiple plausible curves and propagates their posterior uncertainty into the toxicity probabilities. Its limitations are equally important: it cannot include a plausible curve that was left out, posterior weights can concentrate on a model that is only the least inadequate candidate, and averaging does not guarantee good calibration or protect participants from overdose. Protocol eligibility rules and other safety procedures remain necessary.

Common pitfall: Averaging the skeleton entries, or using their prior weights, and then treating the result as the posterior toxicity curve. This ignores the data-dependent posterior probability of each model and the nonlinear CRM transformation. Update \(\gamma\) within each model, compute each model's posterior toxicity probability, and then average those probabilities using \(P(M_m\mid D)\).

Check Your Understanding

Answer each question, then compare with the solutions.

Why does the model-averaged toxicity probability use posterior model weights rather than prior model weights?
Suppose two models have posterior probabilities \(0.7\) and \(0.3\), with posterior mean toxicity probabilities \(0.20\) and \(0.50\) at a dose. What is the model-averaged probability?
With a monotone skeleton and \(\gamma>0\) in the power CRM, should the fitted toxicity probabilities decrease with dose? Explain.
Why is averaging skeleton values and then transforming generally not equivalent to averaging posterior toxicity probabilities?

Solutions. (1) The posterior weights condition on the observed data and represent the updated uncertainty about which model applies; they follow from Bayes' rule. (2) The average is \(0.7(0.20)+0.3(0.50)=0.14+0.15=0.29\), or a 29% posterior predictive toxicity probability. (3) No. Since \(a_{mi}\) increases with dose and \(\gamma>0\), \(a_{mi}^{\gamma}\) also increases, so the probabilities are nondecreasing. (4) The transformation \(a^\gamma\) is nonlinear in general, so transforming a weighted average does not commute with averaging transformed values; moreover, the proper calculation averages posterior probabilities using data-updated weights.