Computation for Commensurate Precision Models

Why computation needs careful parameterization

A commensurate model connects a current-trial parameter to a historical parameter through a precision: large precision concentrates the current parameter near the historical one, while small precision allows greater separation. When that precision is uncertain, computation must represent both the data contributions and the uncertainty in borrowing correctly. In particular, the precision is positive, and changing its computational scale changes the density through a Jacobian.

The model below uses continuous outcomes and known sampling variances to make each contribution visible. The same computational principles apply when the endpoint has a different likelihood: use that endpoint’s correct likelihood, retain the commensurate link, and transform the positive precision consistently. The model assumes the historical and current data are conditionally independent given their respective parameters, and that the specified sampling variances are known.

Definition: Let \(\bar{x}_H\) and \(\bar{x}_C\) be historical and current sample means, based on sample sizes \(n_H\) and \(n_C\). For known within-study variances \(\sigma_H^2\) and \(\sigma_C^2\), specify \(\bar{x}_H\mid\theta_H\sim N(\theta_H,\sigma_H^2/n_H)\), \(\bar{x}_C\mid\theta_C\sim N(\theta_C,\sigma_C^2/n_C)\), \(\theta_H\sim N(m_0,v_0)\), \(\theta_C\mid\theta_H,\tau\sim N(\theta_H,1/\tau)\), and \(\tau\sim\operatorname{Gamma}(a_\tau,b_\tau)\). The Gamma distribution uses shape \(a_\tau\) and rate \(b_\tau\), so its density is proportional to \(\tau^{a_\tau-1}\exp(-b_\tau\tau)\), for \(\tau>0\). Here \(\tau\) is the commensurate precision, and \(1/\tau\) is the conditional variance of \(\theta_C-\theta_H\).

This specification treats the historical estimate as data about \(\theta_H\), not as a fixed prior mean for \(\theta_C\). The current estimate informs \(\theta_C\), and the commensurate distribution links the two parameters. The variance parameters \(v_0,\sigma_H^2,\sigma_C^2\) are variances, not standard deviations. If the original observations are available, their full likelihood can be used; with normal observations and known variance, the sample mean is sufficient for inference about the corresponding mean once any parameter-independent factors are removed.

Assembling the posterior computation

Write \(v_H=\sigma_H^2/n_H\) and \(v_C=\sigma_C^2/n_C\). Conditional independence gives a joint likelihood equal to the product of the historical and current sampling densities. Multiplying by the prior for \(\theta_H\), the commensurate density, and the Gamma prior gives the joint posterior up to its normalizing constant. A convenient computational parameter is \(\eta=\log\tau\), which ranges over the whole real line and guarantees \(\tau=\exp(\eta)>0\).

\begin{aligned} \log \pi(\theta_H,\theta_C,\eta\mid\bar{x}_H,\bar{x}_C) ={}&-\frac{(\bar{x}_H-\theta_H)^2}{2v_H} -\frac{(\bar{x}_C-\theta_C)^2}{2v_C} -\frac{(\theta_H-m_0)^2}{2v_0}\\ &+\left(a_\tau+\frac12\right)\eta -\left(b_\tau+\frac{(\theta_C-\theta_H)^2}{2}\right)\exp(\eta) +\text{constant}. \end{aligned}

The constant does not depend on \(\theta_H,\theta_C,\eta\), so it can be omitted for posterior sampling or for ratios of posterior densities. The term \(\tfrac12\eta-\tfrac12\tau(\theta_C-\theta_H)^2\) comes from the normal commensurate density, whose precision-dependent factor is \(\tau^{1/2}\). The Gamma prior contributes \((a_\tau-1)\eta-b_\tau\exp(\eta)\). Finally, changing variables from \(\tau\) to \(\eta\) contributes a further \(+\eta\), producing the coefficient \(a_\tau+\tfrac12\).

Key computational rule: If the posterior density is evaluated on the \(\eta=\log\tau\) scale, evaluate the density with respect to \(\eta\), including the Jacobian \(\left|d\tau/d\eta\right|=\exp(\eta)=\tau\). Do not add this Jacobian when evaluating a density with respect to \(\tau\) itself.

Worked Example: Building a log-posterior contribution

Suppose a hypothetical historical study has \(n_H=10\), \(\bar{x}_H=10.8\), and known \(\sigma_H=2.5\), while a hypothetical current trial has \(n_C=8\), \(\bar{x}_C=12.4\), and known \(\sigma_C=2\). Let \(m_0=10\), \(v_0=25\), and \(\tau\sim\operatorname{Gamma}(2,0.5)\), with the second Gamma parameter a rate. Evaluate the log posterior kernel at \(\theta_H=11\), \(\theta_C=12\), and \(\eta=0\).

First, \(v_H=2.5^2/10=0.625\) and \(v_C=2^2/8=0.5\). At \(\eta=0\), \(\tau=\exp(0)=1\). The historical likelihood’s quadratic contribution is \(-\tfrac12(10.8-11)^2/0.625=-0.0320\), and the current likelihood’s is \(-\tfrac12(12.4-12)^2/0.5=-0.160\).

The historical prior contributes \(-\tfrac12(11-10)^2/25=-0.0200\). The commensurate quadratic contribution is \(-\tfrac12(12-11)^2(1)=-0.500\). The transformed precision-prior terms are \((2+\tfrac12)(0)-[0.5+\tfrac12(12-11)^2]\exp(0)=-1.00\). Thus the log posterior kernel, omitting constants, is \(-0.0320-0.160-0.0200-0.500-1.00=-1.712\).

This is an unnormalized log density at one parameter state, not a posterior probability. The calculation also shows why each component should be checked separately: the two likelihood terms use the sampling variances of the means, whereas the commensurate term uses the unknown variance \(1/\tau\).

Positive precision: transformed and conditional calculations

The log transformation avoids proposing negative precision values in unconstrained numerical algorithms. Its Jacobian is not optional: probability densities are defined relative to a variable, and the density changes when the variable changes. The next calculation isolates the difference between the density on the precision scale and the density on the log-precision scale.

Worked Example: Applying the log-precision Jacobian

For the hypothetical prior \(\tau\sim\operatorname{Gamma}(2,0.5)\), the normalized density on the \(\tau\) scale is \(f_\tau(\tau)=0.25\tau\exp(-0.5\tau)\), for \(\tau>0\). Because \(\tau=\exp(\eta)\), the density on the \(\eta\) scale is \(f_\eta(\eta)=f_\tau(\exp(\eta))\exp(\eta)\).

At \(\tau=0.25\), the density on the original scale is \(f_\tau(0.25)=0.25(0.25)\exp(-0.125)\approx0.05516\). Since \(\eta=\log(0.25)\), the transformed density is \(f_\eta(\eta)=0.05516(0.25)\approx0.01379\). At \(\tau=4\), \(f_\tau(4)=0.25(4)\exp(-2)\approx0.1353\), whereas \(f_\eta(\log 4)=0.1353(4)\approx0.5413\).

These are densities with respect to different scales, so their numerical values need not match and a density can exceed one. Omitting the factor \(\exp(\eta)\) would give the wrong distribution for the transformed variable, even though the transformed precision remains positive.

This model also has useful conditional distributions, which can support Gibbs sampling or checks on a more general implementation. Conditional on \(\theta_H,\theta_C\), the commensurate normal density contributes a factor proportional to \(\tau^{1/2}\exp[-\tau(\theta_C-\theta_H)^2/2]\). Multiplication by the Gamma prior gives another Gamma kernel:

\tau\mid\theta_H,\theta_C,\text{data} \sim\operatorname{Gamma}\left( a_\tau+\frac12,\; b_\tau+\frac{(\theta_C-\theta_H)^2}{2} \right).

The rate parameter increases with squared disagreement. Holding the other quantities fixed, a larger difference between the current and historical parameters lowers the conditional mean precision, \( (a_\tau+\tfrac12)/(b_\tau+(\theta_C-\theta_H)^2/2)\). This is the model’s direct computational expression of reduced commensurability when the parameter values are far apart.

Worked Example: Deriving and using conditional updates

Use the hypothetical settings above, and suppose a computation is currently at \(\theta_H=11\), \(\theta_C=12.2\). Their difference is \(1.2\). The conditional precision distribution is \(\operatorname{Gamma}(2+\tfrac12,\;0.5+1.2^2/2)=\operatorname{Gamma}(2.5,1.22)\). Its conditional mean is \(2.5/1.22\approx2.049\), and its conditional variance is \(2.5/1.22^2\approx1.679\), using the shape-rate Gamma convention.

The current parameter also has a normal full conditional distribution. Its precision is the sum of the current data precision and the commensurate precision: \(1/v_C+\tau=1/0.5+1=3\). At \(\bar{x}_C=12.4\), \(\theta_H=11\), and \(\tau=1\), its conditional mean and variance are \[ \frac{\bar{x}_C/v_C+\tau\theta_H}{1/v_C+\tau} =\frac{12.4/0.5+11}{3} =\frac{35.8}{3}\approx11.93, \qquad \frac{1}{3}\approx0.3333. \] Thus \(\theta_C\mid\theta_H,\tau,\text{data}\sim N(11.93,0.3333)\), with the second normal parameter a variance.

These updates follow by collecting the terms quadratic in the variable being updated. For \(\theta_C\), the quadratic terms are \(-\tfrac12(\bar{x}_C-\theta_C)^2/v_C-\tfrac12\tau(\theta_C-\theta_H)^2\); completing the square gives the precision-weighted mean above. The update calculations rely on the stated known-variance normal model and Gamma prior. A different likelihood or prior may not yield these same full conditionals.

Implementation choices and failure modes

A computation can use the conditional distributions when available, or evaluate the joint log posterior and use a numerical method that operates on unconstrained coordinates. A fixed value of \(\tau\) removes uncertainty about the degree of borrowing and can be useful as a reference fit, but it does not answer the uncertain-precision question. Numerical integration over \(\tau\) is another option when the remaining parameter dimension is small. In more complex models, sampling is often practical, but the target density must still match the intended model.

For other endpoint types, replace the normal sampling contributions rather than forcing a normal likelihood onto the data. For example, if \(x_C\) events are observed among \(n_C\) current participants with event probability \(p_C\), the binomial log-likelihood contribution, up to a constant, is \(x_C\log(p_C)+(n_C-x_C)\log(1-p_C)\). The historical contribution is written for its own parameter in the same way. The commensurate relationship must then be specified on the chosen parameter scale; a normal link on a probability scale is not interchangeable with one on a logit scale.

Several checks catch common implementation errors before results are interpreted:

Check the scale of every variance. The likelihood for a sample mean uses \(\sigma^2/n\), not \(\sigma^2\), while the commensurate variance is \(1/\tau\).
Check the Gamma convention. The formulas here use a rate. A software function parameterized by scale requires scale \(1/b_\tau\), not \(b_\tau\).
Check the transformed target. On the \(\eta\) scale, the log Jacobian is \(+\eta\). Without it, the implied prior and posterior for \(\tau\) are incorrect.
Check arithmetic on the log scale. Compute log densities and combine them by addition; multiplying very small densities directly can underflow.
Check that data enter once. A historical sample mean likelihood and a prior constructed from the same historical observations should not both be included unless the model explicitly justifies both contributions.
Check numerical extremes. Very large \(\eta\) can make \(\exp(\eta)\) overflow, and very negative \(\eta\) can effectively set the precision to zero. Avoid arbitrary clipping that changes the target; use numerically stable evaluations and investigate why the computation reaches problematic values.

Common mistake: “I sampled \(\eta=\log\tau\), so I only need to substitute \(\exp(\eta)\) for \(\tau\) in the original density.” This omits the change-of-variables factor and targets the wrong distribution. When the target is expressed as a density in \(\eta\), add \(\log|d\tau/d\eta|=\eta\); when it is expressed as a density in \(\tau\), do not add that term.

A useful implementation check is to compare a conditional update with the corresponding completed-square calculation, as in the example, and to compare log-posterior values computed from separate components with a single combined expression. When the model is small enough, numerical integration can provide an additional reference for posterior summaries. These checks assess whether the intended target has been coded; they do not establish that a simulation algorithm explores that target adequately.

Check Your Understanding

Try each question, then compare your reasoning with the solutions.

If \(\tau=\exp(\eta)\), what is the log-Jacobian contribution when evaluating a density on the \(\eta\) scale?
Suppose \(a_\tau=3\), \(b_\tau=2\), and \(\theta_C-\theta_H=2\). Give the conditional distribution of \(\tau\).
For \(v_C=0.25\), \(\tau=2\), \(\bar{x}_C=5\), and \(\theta_H=4\), calculate the conditional mean and variance of \(\theta_C\).
Why is a likelihood for a sample mean based on \(\sigma_C^2/n_C\), rather than \(\sigma_C^2\)?
A programmer substitutes \(\exp(\eta)\) everywhere but omits the Jacobian. What is wrong with this computation?

Solutions. (1) The contribution is \(+\eta\), since \(\log|d\tau/d\eta|=\log(\exp(\eta))=\eta\). (2) The shape is \(3+\tfrac12=3.5\), and the rate is \(2+2^2/2=4\), so \(\tau\mid\theta_H,\theta_C,\text{data}\sim\operatorname{Gamma}(3.5,4)\). (3) The conditional precision is \(1/0.25+2=6\). The mean is \((5/0.25+2(4))/6=28/6\approx4.667\), and the variance is \(1/6\approx0.1667\). (4) The sample mean of \(n_C\) independent observations with variance \(\sigma_C^2\) has variance \(\sigma_C^2/n_C\). (5) It omits the change-of-variables factor, so the density targeted on the log-precision scale is not the transformed density of the specified prior and posterior.

Computation for Commensurate Precision Models

What You'll Learn