Deriving Bayes Factors from Marginal Likelihoods

Why Integrate Likelihoods?

Two PK/PD models can each fit the observed data well at some parameter value, yet differ in how much of their parameter space gives a good fit. A likelihood evaluated at the best-fitting parameter describes that single point. Bayesian model comparison instead averages the likelihood over the parameter values allowed by a model's prior. This average is its marginal likelihood, and the ratio of two marginal likelihoods is a Bayes factor.

The integration matters because a model that fits only in a narrow region of parameter space receives less average support than one that predicts the data well across a substantial part of its prior distribution. Thus a Bayes factor reflects both fit and parameter uncertainty. It is conditional on the candidate models, their likelihoods, and their parameter priors; it is not a prior-free measure of fit.

Definition: For model $M_k$, let $\theta_k$ denote its parameters, $p(y\mid\theta_k,M_k)$ the likelihood for observed data $y$, and $\pi(\theta_k\mid M_k)$ a proper prior density. The marginal likelihood is $$ m_k(y)=p(y\mid M_k)=\int p(y\mid\theta_k,M_k)\,\pi(\theta_k\mid M_k)\,d\theta_k. $$ For two models, the Bayes factor favoring $M_1$ over $M_0$ is $$ BF_{10}=\frac{m_1(y)}{m_0(y)}. $$ A value above 1 means the observed data have greater marginal likelihood under $M_1$ than under $M_0$, given the specified priors.

The prior density integrates to 1, so the marginal likelihood is a prior-weighted average of the likelihood, not an unweighted integral. The two models may have different parameters and different parameter spaces; each likelihood is integrated with respect to its own model-specific prior. The same observed data and data representation must be used for both models. For continuous outcomes, marginal likelihoods are densities and depend on the measurement scale, but a common change of units multiplies both by the same data-only factor, which cancels in their ratio.

Deriving the Ratio

The marginal likelihood follows by integrating the joint density of the data and parameters over the parameters. The joint density factors into a likelihood and a prior:

p(y,\theta_k\mid M_k) =p(y\mid\theta_k,M_k)\,\pi(\theta_k\mid M_k).

Since the possible parameter values are not observed, sum or integrate over them to obtain the probability density of the data under the model. Dividing the resulting data densities gives the Bayes factor. This operation incorporates the uncertainty about parameters under each model rather than selecting a single estimate.

p(y\mid M_k)=\int p(y,\theta_k\mid M_k)\,d\theta_k =\int p(y\mid\theta_k,M_k)\pi(\theta_k\mid M_k)\,d\theta_k, \qquad BF_{10}=\frac{\int p(y\mid\theta_1,M_1)\pi(\theta_1\mid M_1)\,d\theta_1} {\int p(y\mid\theta_0,M_0)\pi(\theta_0\mid M_0)\,d\theta_0}.

This expression also clarifies why substituting maximum-likelihood estimates into the likelihood does not, by itself, produce a Bayes factor: it replaces an average over the prior with a value at one point. The integration can be analytic for some models, or approximated numerically when the likelihood and prior do not yield a closed form.

Worked Example: A Normal Response With an Uncertain Drug Effect

Suppose a hypothetical PK/PD analysis summarizes a response at a specified exposure as a deviation $y=4$ response units from baseline. For illustration, treat the residual standard deviation as known at 1 response unit. Model $M_0$ has no exposure effect and predicts $y\sim N(0,1)$. Model $M_1$ predicts $y\mid\beta\sim N(\beta,1)$, where the uncertain exposure effect has prior $\beta\sim N(0,4)$. Calculate $BF_{10}$.

Under $M_0$, the marginal likelihood is just the normal density at $y=4$:

m_0(4)=\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{4^2}{2}\right) =\frac{1}{\sqrt{2\pi}}\exp(-8) \approx 0.0001338.

Under $M_1$, integrating over $\beta$ is equivalent to adding the variance of the prior effect to the residual variance. Thus $y\sim N(0,1+4)=N(0,5)$, and

m_1(4)=\frac{1}{\sqrt{2\pi(5)}}\exp\left(-\frac{4^2}{2(5)}\right) =\frac{1}{\sqrt{10\pi}}\exp(-1.6) \approx 0.03602.

Therefore,

BF_{10} =\frac{m_1(4)}{m_0(4)} =\sqrt{\frac{1}{5}}\exp(8-1.6) =\sqrt{0.2}\,\exp(6.4) \approx 269.2.

For these hypothetical data and priors, the marginal likelihood under the uncertain-effect model is about 269 times that under the no-effect model. The alternative does not receive this support simply because it can fit $y=4$: its prior assigns probability across many possible effects, including values that predict the observed response less well. The integrated result accounts for that uncertainty.

A PK/PD Example With a Binary Response

For a binary pharmacodynamic endpoint, a simple model can represent the probability $q$ of response at a specified exposure. The likelihood for one observed responder is $q$. The marginal likelihood under an uncertain $q$ is therefore the prior average of $q$. This setup is deliberately simple: it isolates the integration that also appears within more elaborate exposure-response models.

Worked Example: Averaging a Prior for an Exposure-Response Probability

In a hypothetical study, one subject is observed at a specified exposure and has a binary response ($y=1$). Under $M_0$, the response probability is fixed at $q_0=0.01111$. Under $M_1$, the prior assigns probability $0.5$ to $q=0.2000$ and probability $0.5$ to $q=0.6472$. These two prior masses sum to 1, so the prior is proper. Calculate the Bayes factor favoring $M_1$.

The likelihood for the observed response is $p(y=1\mid q)=q$. Under $M_0$, there is no parameter to integrate:

$$ m_0(1)=0.01111. $$

Under $M_1$, average the likelihood across the two prior values:

$$ m_1(1) =(0.5)(0.2000)+(0.5)(0.6472) =0.1000+0.3236 =0.4236. $$

The resulting Bayes factor is

BF_{10}=\frac{0.4236}{0.01111}\approx 38.13.

The observed response is about 38 times as probable under the prior-averaged alternative as under the specified null model. This calculation uses the likelihood at each prior value, not the likelihood at whichever value best fits the response. A discrete prior is used here to make the averaging transparent; in continuous PK/PD models the same definition uses an integral.

Continuous Integration and Parameter Uncertainty

A continuous prior makes the averaging explicit as an integral. For a binary response probability, a Beta prior is convenient because its density is proportional to a power of $q$ times a power of $1-q$. The integral of that form is a beta function, which gives a direct calculation without choosing a single fitted value of $q$.

Worked Example: Ten Hypothetical Subjects at One Exposure

Suppose 4 of 10 hypothetical subjects respond at a specified exposure. Compare $M_0$, which fixes $q=0.2000$, with $M_1$, which assigns $q\sim\operatorname{Beta}(2,2)$. For the observed pattern of four responses and six nonresponses, the likelihood is $q^4(1-q)^6$. The binomial coefficient for treating the data as a count would be the same under both models and cancels from the Bayes factor.

Under $M_0$, the marginal likelihood for that response pattern is

m_0(y)=(0.2)^4(0.8)^6 =0.0016\times0.262144 =0.0004194304.

The $\operatorname{Beta}(2,2)$ prior density is $q(1-q)/B(2,2)$, where $B(a,b)=\int_0^1 q^{a-1}(1-q)^{b-1}\,dq$. Integrating the likelihood gives

m_1(y) =\frac{1}{B(2,2)}\int_0^1 q^{4+1}(1-q)^{6+1}\,dq =\frac{B(6,8)}{B(2,2)} =\frac{1/10296}{1/6} =\frac{1}{1716} \approx 0.0005828.

Therefore,

BF_{10} =\frac{0.0005828}{0.0004194304} \approx 1.389.

The integrated evidence modestly favors the model with uncertain response probability, rather than overwhelming support. Although the maximum-likelihood estimate under $M_1$ is $4/10=0.4$, the Bayes factor is not a comparison of the likelihood at $q=0.4$ with the null likelihood. It averages over the entire $\operatorname{Beta}(2,2)$ prior, including values that do not fit the observed proportion as closely.

Computation, Diagnostics, and Alternatives

For a nonlinear PK/PD model, the marginal likelihood may require numerical integration over several parameters, such as an exposure-response maximum, an exposure scale, and residual variability. A useful computational outline is:

Specify the candidate model and prior.
Write the likelihood for the observed data and a normalized prior over every parameter being integrated. Confirm that the prior is proper.

Evaluate the log likelihood and log prior.
Use the same observations and response definition across models. Log-scale calculations help prevent numerical underflow when likelihoods are small.

Approximate each integral and check it.
Use an appropriate method, such as quadrature or a validated marginal-likelihood estimator. Check numerical stability, integration coverage, and agreement across adequate independent calculations.

Form the ratio and state its direction.
Calculate $m_1(y)/m_0(y)$, identify which model is in the numerator, and interpret the result only for the specified models, likelihoods, and priors.

Warning signs include a marginal likelihood that changes materially when numerical settings are tightened, failure to cover important tails or modes, and disagreement between independent calculations. Parameter transformations also require care: if integration is performed in transformed coordinates, the appropriate Jacobian must be included so the prior density represents the intended distribution. Improper priors generally do not define normalized prior predictive densities, so their arbitrary normalizing constants can make a Bayes factor undefined.

Bayes factors answer a different question from common alternatives. A likelihood-ratio test compares maximized likelihoods and, for nested models, uses a reference distribution under its assumptions. AICc and BIC summarize fit with a complexity adjustment; they are not the same as integrating the likelihood over a specified prior. Subject-level cross-validation assesses predictive performance on held-out subjects. These approaches can inform different scientific goals, but their numerical criteria should not be treated as interchangeable.

Common pitfall: “The more flexible model has the larger maximum likelihood, so its Bayes factor must exceed 1.” This is incorrect: a maximum is not a prior average. The correction is to integrate the likelihood over the model's proper parameter prior, then compare the resulting marginal likelihoods. A model can have a high peak but a low average if good fit occurs only in a small prior region.

A Bayes factor cannot repair a scientifically inadequate candidate set or an implausible likelihood. Its interpretation is conditional on the assumed response model, residual structure, and parameter priors. It is evidence for one specified model relative to another, not proof that either model is true, and it does not by itself establish a causal drug effect.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

In the definition of a marginal likelihood, what is integrated and what weights the likelihood values?
If $m_1(y)=0.12$ and $m_0(y)=0.03$, what is $BF_{10}$, and what does it describe?
Why does evaluating the alternative likelihood only at its maximum not generally give a Bayes factor?
Why can an improper parameter prior cause a problem for Bayes-factor calculation?

Solutions. (1) Integrate the likelihood over the model's parameter space, weighted by its normalized prior density. (2) $BF_{10}=0.12/0.03=4$; the observed data have four times the marginal likelihood under $M_1$ as under $M_0$, given the stated models and priors. (3) A maximum is one likelihood value, whereas a marginal likelihood averages over parameter values under the prior. (4) An improper prior has no finite normalizing constant, so the prior predictive density may be undefined and arbitrary constants may fail to cancel across models.