James E. Pustejovsky
July 12, 2026
For a new project, I’ve been thinking about model evaluation and model selection in meta-analysis. You can’t get very far into the literature on this without running into Kullback-Leibler divergence, which is a metric for quantifying discrepancies between distributions. Among other things, KL divergence is the basis of information criteria, including the Akaike (Akaike, 1973), Bayesian (Schwarz, 1978), Deviance (Spiegelhalter et al., 2002), Watanabe (Watanabe, 2010), and LOO (Vehtari & Lampinen, 2002) information criteria. For two distributions \(P\) and \(Q\) with probability density functions \(p(x)\) and \(q(x)\), respectively, the KL divergence of Q from P is \[ KL(P||Q) = \mathbb{E}_P\left[\log p(X) - \log q(X)\right], \] where \(\mathbb{E}_P\) means that the expectation is taken with respect to \(P\). Note that KL-divergence is not, strictly speaking, a distance metric because it is not symmetric: \(KL(P||Q) \neq KL(Q || P)\) because the expectation is taken under different distributions in each case.
One of the vexing things about working with KL divergence is interpreting its magnitude. If two models differ by 4 in KL divergence, is that big and meaningful? Or trivial? How should an analyst make judgements on the basis of KL divergence, AICs, or LOOICs? There are surely many answers to this question. In this post, I’ll explore a technique that’s been called calibration.
Suppose we are examining two different models fitted to the same data, trying to decide which one fits better. The (estimated) KL divergences of the two models differ by \(\Delta\). Calibration involves constructing an analogy to these two models by finding simpler-to-conceptualize models with the same KL divergence, \(KL = \Delta\). McCulloch (1989) suggested calibrating based on the KL divergence of a Bernoulli\((\pi)\) distribution from a Bernoulli\((\frac{1}{2})\) distribution. Domingue et al. (2025) recently described a relative model fit index called the InterModel Vigorish, which is derived along essentially the same lines as McCulloch (1989)’s calibration model. Goutis & Robert (1998) gave examples of calibrations based on other single-parameter distributions including the divergence of Poisson\((\lambda)\) from Poisson\((1)\) and divergence of \(N(\mu,1)\) from \(N(0,1)\). The diversity of potential calibration models is effectively endless, constrained only by one’s imagination and computing power.
Here is a slight generalization of the calibration model proposed in McCulloch (1989), where I use Binomial instead of Bernoulli distributions. Suppose that the true model \(P\) is Binomial\(\left(N,\frac{1}{2}\right)\) but our prediction model \(Q\) is \(\left(N,\pi\right)\). The KL divergence of the latter from the former is \[ \begin{aligned} KL\left(Bin\left(N,\frac{1}{2}\right), Bin\left(N,\pi\right)\right) &= N \times \left[\log\left(\frac{1}{2}\right) -\frac{1}{2} \left[\log \pi + \log(1 - \pi)\right]\right] \\ &= -\frac{K}{2} \log\left(4 \pi(1 - \pi)\right). \end{aligned} \] Given a KL-divergence of \(\Delta\), the prediction model with corresponding KL-divergence would have \(\pi = \frac{1}{2}\left[1 + \sqrt{1 - \exp\left(-\frac{2\Delta}{N}\right)}\right]\).
In the interactive Figure 1, you can specify a degree of divergence \(\Delta\) (along with the number of trials \(N\)) and see the binomial distribution with probability \(\pi\) that has KL-divergence of \(\Delta\) from a binomial distribution with probability \(\frac{1}{2}\). You can think of this as representing the divergence (or lack of fit) that corresponds to assuming you have a biased \(\pi\) coin across \(N\) flips, when the coin is actually fair.
math = require("mathjs")
binomial = require( 'https://cdn.jsdelivr.net/gh/stdlib-js/stats-base-dists-binomial@umd/browser.js' )
pi = (1 + math.sqrt(1 - math.exp(-2 * Delta / N))) / 2
pts = N + 1
dat1 = Array(pts).fill().map((element, index) => {
let x = index;
let p = binomial.pmf(index, N, 0.5);
let q = binomial.pmf(index, N, pi);
return ({
x: x,
p: p,
q: q
})
})McCulloch’s calibration model feels backwards to me—why would one predict that a coin is biased when it is actually perfectly fair? The other way to set up this problem would be to imagine that you are flipping \(N\) coins and that you expect the coin to be fair, but it is actually biased. This corresponds to finding the probability for \(Q\) such that the KL divergence of \(P\) from \(Q(\phi)\) is equal to \(\Delta\): \[ KL\left(Bin\left(N,\phi\right), Bin\left(N,\frac{1}{2}\right)\right) = N \times \left[\log 2 + \phi \log \phi + (1 - \phi) \log(1 - \phi)\right]. \] This is a bit harder to solve for \(\phi\), but no problem to work out numerically. Figure 2 displays the calibrated models corresponding to a given difference \(\Delta\) in KL divergence, for whatever value of \(N\) that you like. For small or moderate \(N\), this calibration model is distinct from McCulloch’s because of the asymmetry of KL divergence.
function newtonRaphson(f, df, c, x0, xmin, xmax, tolerance, maxIterations) {
let x = x0;
let iter = 0;
let error = Number.POSITIVE_INFINITY;
let fmin = f(xmin) - c;
let fmax = f(xmax) - c;
if (math.sign(fmin) == math.sign(fmax)) {
console.warn("The function does not have a solution within the specified range.");
if (math.abs(fmin) < math.abs(fmax)) {
return(xmin);
} else {
return(xmax);
}
} else {
while (error > tolerance && iter < maxIterations) {
let fx = f(x) - c;
let dfx = df(x);
let step = -fx / dfx;
let p = 0;
let xNext = x + step;
while (xNext <= xmin || xNext >= xmax) {
xNext = x + step / 2**p;
p++;
}
error = Math.abs(fx);
x = xNext;
iter++;
}
if (iter === maxIterations) {
console.warn("Newton-Raphson method did not converge within the maximum number of iterations.");
}
return x;
}
}
function KL_binomial(phi) {
return N * (math.log(2) + phi * math.log(phi) + (1 - phi) * math.log(1 - phi));
}
function dKL_binomial(phi) {
return N * (math.log(phi) - math.log(1 - phi))
}
phi = newtonRaphson(KL_binomial, dKL_binomial, Delta, 0.75, 0.50, 0.9999, 1e-5, 100)
dat2 = Array(pts).fill().map((element, index) => {
let x = index;
let p = binomial.pmf(index, N, phi);
let q = binomial.pmf(index, N, 0.5);
return ({
x: x,
p: p,
q: q
})
})The appeal of coin flips (or binomial distributions) for this sort of calibration exercise is that it they are familiar and just about as simple as you can get—one bit of information per flip—as far as probability models go. However, this isn’t necessarily the most intuitive calibration model for meta-analysis. In a conventional random effects meta-analysis, the effect size parameters are typically assumed to be normally distributed about some average effect and with some standard deviation that measures the degree of heterogeneity in the results across studies. Given that meta-analysts are accustomed to thinking about Gaussian distributions, perhaps a more natural calibration model would be one that involves the standard normal distribution. Here I’ll treat the standard normal \(N(0,1)\) distribution as the inferior prediction model and consider a variety of superior models for the unknown true distribution. For greater flexibility (as well as consistency with the coin flip models), I’ll always write KL divergences for samples of \(N\) observations, which is just \(N\) times the KL divergence of a single observation.1
A simple starting point would be to compare two normal distributions to each other. The KL divergence of a \(N(\mu,\sigma^2)\) distribution from the standard normal distribution is quite simple: \[ KL_N\left(N\left(\mu,\sigma^2\right), N\left(0,1\right)\right) = \frac{N}{2}\left(\mu^2 + \sigma^2 - 1 - \log \sigma^2\right). \] Even simpler, Goutis & Robert (1998) considered a calibration model involving normal distributions with unit variance, so the KL divergence is just \(N \times \mu^2 / 2\). Thus, taking \(\mu = \sqrt{\frac{2\Delta}{N}}\) leads to two unit-Gaussian distributions for which an \(N\)-sample has a divergence of \(\Delta\).
But why focus only on location shifts? Instead, what if the true model \(P\) had mean zero but differing variance? Then \[ KL_N\left(N\left(0,\sigma^2\right), N\left(0,1\right)\right) = \frac{N}{2}\left(\sigma^2 - 1 - \log \sigma^2\right). \] To obtain a KL divergence of \(\Delta\), set \(\sigma^2 = - W\left(-\exp\left[-\left(\frac{2\Delta}{N} + 1\right)\right]\right)\), where \(W()\) is the Lambert W function.
Neither of these calibration models seems especially salient in the meta-analysis context. There, it seems rather implausible that a model would be badly mis-specified in terms of the center or spread of the distribution of effects. I would find it more plausible that a model might be mis-specified in terms of its shape somehow—perhaps missing some skewness or kurtosis in the distribution of effects that a (Gaussian) random effects model would not capture.
What if the superior model was a shifted gamma distribution with shape parameter \(\alpha\), scale parameter \(\beta\), and shift \(c\)? If \(X + c \sim \Gamma(\alpha, \beta)\), then setting \(\beta = \sqrt{\alpha}\) and \(c = -\sqrt{\alpha}\) leads to a distribution with mean zero, unit variance, skewness \(2 / \sqrt{\alpha}\), and differential entropy \[ \mathcal{H}(P) = \alpha - \frac{1}{2}\log (\alpha) + \log \Gamma(\alpha) + (1 - \alpha) \psi(\alpha), \] where \(\Gamma()\) is the gamma function and \(\psi()\) is the digamma function. The KL divergence of an \(N\)-sample from this shifted gamma distribution to the unit normal distribution is therefore \[ \begin{aligned} KL_N\left(\Gamma\left(\alpha,\sqrt\alpha\right) + \sqrt\alpha, N\left(0,1\right)\right) &= \frac{1}{2}\log(2\pi) + \frac{1}{2} - \mathcal{H}(P) \\ &= \frac{1}{2}\log(2\pi) + \frac{1}{2} - \alpha + \frac{1}{2}\log (\alpha) - \log \Gamma(\alpha) - (1 - \alpha) \psi(\alpha). \end{aligned} \] Although there’s not a closed-form inverse, it’s easy enough to solve numerically for \(\alpha\).
If mis-specification could manifest as skewness, then why not kurtosis too? Suppose the true model is a Student \(t\) distribution with degrees of freedom \(\nu\). For simplicity, I won’t worry about forcing the \(t\) to have unit variance. For an \(N\)-sample, the KL-divergence of the standard normal from a \(t(\nu)\) distribution is \[ \begin{aligned} KL_N\left(t(\nu), N\left(0,1\right)\right) &= \frac{N}{2}\left[\log(2\pi) + \frac{\nu}{\nu - 2} \right. \\ & \qquad \qquad \left. - (\nu + 1)\left(\psi\left(\frac{\nu + 1}{2}\right) - \psi\left(\frac{\nu}{2}\right)\right) \right. \\ & \qquad \qquad \qquad \qquad \qquad\left. - \log(\nu) - 2 \log B\left(\frac{\nu}{2}, \frac{1}{2}\right) \right], \end{aligned} \] where \(\psi()\) is the digamma function and \(B()\) is the beta function. For a given \(\Delta\), the expression can be solved numerically for \(\nu\).
The interactive Figure 3 displays all four of these calibration models, where the superior distributions differ by location, scale, skew, or kurtosis, respectively. One thing that jumps out to me about this visualization is that the divergence from the \(t(\nu)\) distribution seems much subtler than the others. With this distribution, all of the action is in the tails. Put differently, the KL-divergence is driven by differences in the probability assigned to rare extreme outcomes.
normal = require( 'https://cdn.jsdelivr.net/gh/stdlib-js/stats-base-dists-normal@umd/browser.js' )
student = require( 'https://cdn.jsdelivr.net/gh/stdlib-js/stats-base-dists-t@umd/browser.js' )
gamma = require( 'https://cdn.jsdelivr.net/gh/stdlib-js/stats-base-dists-gamma@umd/browser.js' )
trigamma = require( 'https://cdn.jsdelivr.net/gh/stdlib-js/math-base-special-trigamma@umd/browser.js' )
// Location-shift
mu = math.sqrt(2 * Delta / N)
// Scale-shift
function KL_sigma(s) {
return N * (s**2 - 1 - 2 * math.log(s)) / 2;
}
function dKL_sigma(s) {
return (s - 1 / s) * N;
}
sigma = newtonRaphson(KL_sigma, dKL_sigma, Delta, 1.1, 1, 100, 1e-5, 100)
// Skewness
function logKL_gamma(la) {
let a = math.exp(la);
let c1 = (math.log(2 * math.pi) + 1) / 2;
let c2 = gamma.entropy(a, math.sqrt(a));
return math.log(N) + math.log(c1 - c2);
}
function dlogKL_gamma(la) {
let a = math.exp(la);
const c1 = (math.log(2 * math.pi) + 1) / 2;
let c2 = gamma.entropy(a, math.sqrt(a));
let c3 = 1 / (2 * a) - 1;
let c4 = (1 - a) * trigamma(a);
return a * (c3 - c4) / (c1 - c2);
}
log_Delta = math.log(Delta)
log_alpha = newtonRaphson(logKL_gamma, dlogKL_gamma, log_Delta, -log_Delta, -10, 10, 1e-5, 100)
alpha = math.exp(log_alpha)
// Kurtosis
function KL_t(v) {
let c1 = math.log(2 * math.pi) / 2;
let c2 = v / (2 * (v - 2));
let c3 = student.entropy(v);
return N * (c1 + c2 - c3);
}
function dKL_t(v) {
let c1 = (v - 2)**(-2);
let c2 = 1 / (2 * v);
let c3 = (trigamma((v + 1) / 2) - trigamma(v / 2)) * (v + 1) / 4;
return -N * (c1 + c2 + c3);
}
nu = newtonRaphson(KL_t, dKL_t, Delta, 3, 2, 100, 1e-5, 100)
// Compute densities
pts2 = 201
dat3 = Array(pts2).fill().map((element, index) => {
let x = - scale + index * 2 * scale / (pts2 - 1);
let p = normal.pdf(x,0,1);
let q = normal.pdf(x, mu, 1);
let r = normal.pdf(x,0, sigma);
let s = gamma.pdf(x + math.sqrt(alpha), alpha, math.sqrt(alpha));
let t = student.pdf(x, nu);
return ({
x: x,
p: p,
q: q,
r: r,
s: s,
t: t
})
})
colormap = [{v:"q",c:"blue"},{v:"r",c:"green"},{v:"s",c:"orange"},{v:"t",c:"red"}]
function density_plot(v, c) {
return Plot.plot({
height: 200,
width: 400,
y: {
grid: false,
axis: false,
domain: [0, 0.45]
},
x: {
label: null
},
marks: [
Plot.lineY(dat3, {x: "x", y: "p", stroke: "black"}),
Plot.areaY(dat3, {x: "x", y: v, fill: c, fillOpacity: 0.5}),
Plot.lineY(dat3, {x: "x", y: v, stroke: c})
]
});
}
plot_list = Object.values(colormap).map(x => density_plot(x.v, x.c))How useful are these calibration models? I’m not entirely sure yet.2 I find them appealing because they turn the very abstract concept of KL-divergence into something that can be visualized. However, having to make arbitrary choices of distributions to construct the calibration model does seem a bit weird to me. I will just leave this up and let it marinate a bit.
All of the math that follows amounts to applications of the same general derivation. For a distribution \(P\), \[ \begin{aligned} KL_N\left(P, N\left(0,1\right)\right) &= \mathbb{E}_P\left[\log p(X) - \log \phi(X)\right] \\ &= \mathbb{E}_P\left[\log p(X) + \frac{1}{2}\log(2\pi) + \frac{1}{2} X^2\right] \\ &= \frac{1}{2}\log(2\pi) + \frac{1}{2} \left(\left[\mathbb{E}_P (X) \right]^2 + \mathbb{V}_P(X)\right) - \mathcal{H}(P), \end{aligned} \] where \(\mathcal{H}(P)\) is the differential entropy of the distribution \(P\). Thus, the KL divergence depends only on the entropy, mean, and variance of the distribution. If \(P\) is constructed to have the same mean and variance as the standard normal, then things simplify even further to \(KL_N\left(P, N\left(0,1\right)\right) = \frac{1}{2}\log(2\pi) + \frac{1}{2} - \mathcal{H}(P)\).↩︎
If you have thoughts on these models, please feel free to comment!↩︎
---
title: Calibration models for Kullback-Leibler divergence
date: '2026-07-12'
categories:
- distribution-theory
- Kullback-Leibler-divergence
code-fold: true
code-tools: true
toc: true
bibliography: "../information-criteria-references.bib"
csl: "../apa.csl"
css: styles.css
link-citations: true
---
For a [new project](/presentations/SRSM-2026-effect-metric-choice.html), I've been thinking about model evaluation and model selection in meta-analysis.
You can't get very far into the literature on this without running into Kullback-Leibler divergence, which is a metric for quantifying discrepancies between distributions.
Among other things, KL divergence is the basis of information criteria, including the Akaike [@Akaike1973information], Bayesian [@Schwarz1978estimating], Deviance [@Spiegelhalter2002bayesian], Watanabe [@Watanabe2010asymptotic], and LOO [@Vehtari2002bayesian] information criteria.
For two distributions $P$ and $Q$ with probability density functions $p(x)$ and $q(x)$, respectively, the KL divergence of Q from P is
$$
KL(P||Q) = \mathbb{E}_P\left[\log p(X) - \log q(X)\right],
$$
where $\mathbb{E}_P$ means that the expectation is taken with respect to $P$.
Note that KL-divergence is not, strictly speaking, a distance metric because it is not symmetric: $KL(P||Q) \neq KL(Q || P)$ because the expectation is taken under different distributions in each case.
One of the vexing things about working with KL divergence is interpreting its magnitude. If two models differ by 4 in KL divergence, is that big and meaningful? Or trivial?
How should an analyst make judgements on the basis of KL divergence, AICs, or LOOICs?
There are surely many answers to this question.
In this post, I'll explore a technique that's been called _calibration_.
Suppose we are examining two different models fitted to the same data, trying to decide which one fits better.
The (estimated) KL divergences of the two models differ by $\Delta$.
Calibration involves constructing an analogy to these two models by finding simpler-to-conceptualize models with the same KL divergence, $KL = \Delta$.
@McCulloch1989local suggested calibrating based on the KL divergence of a Bernoulli$(\pi)$ distribution from a Bernoulli$(\frac{1}{2})$ distribution.
@Domingue2025intermodel recently described a relative model fit index called the InterModel Vigorish, which is derived along essentially the same lines as @McCulloch1989local's calibration model.
@Goutis1998model gave examples of calibrations based on other single-parameter distributions including the divergence of Poisson$(\lambda)$ from Poisson$(1)$ and divergence of $N(\mu,1)$ from $N(0,1)$.
The diversity of potential calibration models is effectively endless, constrained only by one's imagination and computing power.
## Coin flip calibrations
Here is a slight generalization of the calibration model proposed in @McCulloch1989local, where I use Binomial instead of Bernoulli distributions.
Suppose that the true model $P$ is Binomial$\left(N,\frac{1}{2}\right)$ but our prediction model $Q$ is $\left(N,\pi\right)$.
The KL divergence of the latter from the former is
$$
\begin{aligned}
KL\left(Bin\left(N,\frac{1}{2}\right), Bin\left(N,\pi\right)\right) &= N \times \left[\log\left(\frac{1}{2}\right) -\frac{1}{2} \left[\log \pi + \log(1 - \pi)\right]\right] \\
&= -\frac{K}{2} \log\left(4 \pi(1 - \pi)\right).
\end{aligned}
$$
Given a KL-divergence of $\Delta$, the prediction model with corresponding KL-divergence would have $\pi = \frac{1}{2}\left[1 + \sqrt{1 - \exp\left(-\frac{2\Delta}{N}\right)}\right]$.
In the interactive @fig-mcculloch, you can specify a degree of divergence $\Delta$ (along with the number of trials $N$) and see the binomial distribution with probability $\pi$ that has KL-divergence of $\Delta$ from a binomial distribution with probability $\frac{1}{2}$.
You can think of this as representing the divergence (or lack of fit) that corresponds to assuming you have a biased $\pi$ coin across $N$ flips, when the coin is actually fair.
::::: {.grid .column-page}
:::: {.g-col-9 .center}
```{ojs}
//| echo: false
math = require("mathjs")
binomial = require( 'https://cdn.jsdelivr.net/gh/stdlib-js/stats-base-dists-binomial@umd/browser.js' )
pi = (1 + math.sqrt(1 - math.exp(-2 * Delta / N))) / 2
pts = N + 1
dat1 = Array(pts).fill().map((element, index) => {
let x = index;
let p = binomial.pmf(index, N, 0.5);
let q = binomial.pmf(index, N, pi);
return ({
x: x,
p: p,
q: q
})
})
```
```{ojs}
//| echo: false
//| label: fig-mcculloch
//| fig-cap: "Binomial$(K,\\pi)$ distribution (black lollipops) that diverges by $\\Delta$ from a Binomial$(K,\\frac{1}{2})$ distribution (blue bars)"
Plot.plot({
height: 350,
width: 700,
y: {
grid: false,
label: "Probability"
},
x: {
label: "X"
},
marks: [
Plot.barY(dat1, {x: "x", y: "p", fill: "blue", fillOpacity: 0.2}),
Plot.ruleX(dat1, {x: "x", y: "q", strokeWidth: 3}),
Plot.dot(dat1, {x: "x", y: "q", fill: "black", r: 5})
]
})
```
::::
:::: {.g-col-3}
```{ojs}
//| echo: false
//| output: false
Nmax = 50
Delta_max = 20
Delta_step = 0.1
viewof N = Inputs.input(10)
viewof Delta = Inputs.input(4)
```
```{ojs}
//| panel: input
//| echo: false
Inputs.bind(
Inputs.range(
[1, Nmax],
{step: 1, label: tex`N`}
),
viewof N
)
Inputs.bind(
Inputs.range(
[0, Delta_max],
{step: Delta_step, label: tex`\Delta`}
),
viewof Delta
)
pi_toprint = pi.toFixed(3)
tex`\pi = ${pi_toprint}`
```
::::
:::::
McCulloch's calibration model feels backwards to me---why would one predict that a coin is biased when it is actually perfectly fair?
The other way to set up this problem would be to imagine that you are flipping $N$ coins and that you expect the coin to be fair, but it is actually biased.
This corresponds to finding the probability for $Q$ such that the KL divergence of $P$ from $Q(\phi)$ is equal to $\Delta$:
$$
KL\left(Bin\left(N,\phi\right), Bin\left(N,\frac{1}{2}\right)\right) = N \times \left[\log 2 + \phi \log \phi + (1 - \phi) \log(1 - \phi)\right].
$$
This is a bit harder to solve for $\phi$, but no problem to work out numerically. @fig-binomial displays the calibrated models corresponding to a given difference $\Delta$ in KL divergence, for whatever value of $N$ that you like.
For small or moderate $N$, this calibration model is distinct from McCulloch's because of the asymmetry of KL divergence.
::::: {.grid .column-page}
:::: {.g-col-9 .center}
```{ojs}
//| echo: false
function newtonRaphson(f, df, c, x0, xmin, xmax, tolerance, maxIterations) {
let x = x0;
let iter = 0;
let error = Number.POSITIVE_INFINITY;
let fmin = f(xmin) - c;
let fmax = f(xmax) - c;
if (math.sign(fmin) == math.sign(fmax)) {
console.warn("The function does not have a solution within the specified range.");
if (math.abs(fmin) < math.abs(fmax)) {
return(xmin);
} else {
return(xmax);
}
} else {
while (error > tolerance && iter < maxIterations) {
let fx = f(x) - c;
let dfx = df(x);
let step = -fx / dfx;
let p = 0;
let xNext = x + step;
while (xNext <= xmin || xNext >= xmax) {
xNext = x + step / 2**p;
p++;
}
error = Math.abs(fx);
x = xNext;
iter++;
}
if (iter === maxIterations) {
console.warn("Newton-Raphson method did not converge within the maximum number of iterations.");
}
return x;
}
}
function KL_binomial(phi) {
return N * (math.log(2) + phi * math.log(phi) + (1 - phi) * math.log(1 - phi));
}
function dKL_binomial(phi) {
return N * (math.log(phi) - math.log(1 - phi))
}
phi = newtonRaphson(KL_binomial, dKL_binomial, Delta, 0.75, 0.50, 0.9999, 1e-5, 100)
dat2 = Array(pts).fill().map((element, index) => {
let x = index;
let p = binomial.pmf(index, N, phi);
let q = binomial.pmf(index, N, 0.5);
return ({
x: x,
p: p,
q: q
})
})
```
```{ojs}
//| echo: false
//| label: fig-binomial
//| fig-cap: "Binomial$(K,\\frac{1}{2})$ distribution (black lollipops) that diverges by $\\Delta$ from a Binomial$(K,\\pi)$ distribution (green bars)"
Plot.plot({
height: 350,
width: 700,
y: {
grid: false,
label: "Probability"
},
x: {
label: "X"
},
marks: [
Plot.barY(dat2, {x: "x", y: "p", fill: "green", fillOpacity: 0.2}),
Plot.ruleX(dat2, {x: "x", y: "q", strokeWidth: 3}),
Plot.dot(dat2, {x: "x", y: "q", fill: "black", r: 5})
]
})
```
::::
:::: {.g-col-3}
```{ojs}
//| panel: input
//| echo: false
Inputs.bind(
Inputs.range(
[1, Nmax],
{step: 1, label: tex`N`}
),
viewof N
)
Inputs.bind(
Inputs.range(
[0, Delta_max],
{step: Delta_step, label: tex`\Delta`}
),
viewof Delta
)
phi_toprint = phi.toFixed(3)
tex`\phi = ${phi_toprint}`
```
::::
:::::
## Standard Normal Calibrations
The appeal of coin flips (or binomial distributions) for this sort of calibration exercise is that it they are familiar and just about as simple as you can get---one bit of information per flip---as far as probability models go.
However, this isn't necessarily the most intuitive calibration model for meta-analysis.
In a conventional random effects meta-analysis, the effect size parameters are typically assumed to be normally distributed about some average effect and with some standard deviation that measures the degree of heterogeneity in the results across studies.
Given that meta-analysts are accustomed to thinking about Gaussian distributions, perhaps a more natural calibration model would be one that involves the standard normal distribution.
Here I'll treat the standard normal $N(0,1)$ distribution as the inferior prediction model and consider a variety of superior models for the unknown true distribution.
For greater flexibility (as well as consistency with the coin flip models), I'll always write KL divergences for samples of $N$ observations, which is just $N$ times the KL divergence of a single observation.[^cross-entropy]
[^cross-entropy]: All of the math that follows amounts to applications of the same general derivation. For a distribution $P$,
$$
\begin{aligned}
KL_N\left(P, N\left(0,1\right)\right) &= \mathbb{E}_P\left[\log p(X) - \log \phi(X)\right] \\
&= \mathbb{E}_P\left[\log p(X) + \frac{1}{2}\log(2\pi) + \frac{1}{2} X^2\right] \\
&= \frac{1}{2}\log(2\pi) + \frac{1}{2} \left(\left[\mathbb{E}_P (X) \right]^2 + \mathbb{V}_P(X)\right) - \mathcal{H}(P),
\end{aligned}
$$
where $\mathcal{H}(P)$ is the differential entropy of the distribution $P$. Thus, the KL divergence depends only on the entropy, mean, and variance of the distribution. If $P$ is constructed to have the same mean and variance as the standard normal, then things simplify even further to $KL_N\left(P, N\left(0,1\right)\right) = \frac{1}{2}\log(2\pi) + \frac{1}{2} - \mathcal{H}(P)$.
A simple starting point would be to compare two normal distributions to each other. The KL divergence of a $N(\mu,\sigma^2)$ distribution from the standard normal distribution is quite simple:
$$
KL_N\left(N\left(\mu,\sigma^2\right), N\left(0,1\right)\right) = \frac{N}{2}\left(\mu^2 + \sigma^2 - 1 - \log \sigma^2\right).
$$
Even simpler, @Goutis1998model considered a calibration model involving normal distributions with unit variance, so the KL divergence is just $N \times \mu^2 / 2$. Thus, taking $\mu = \sqrt{\frac{2\Delta}{N}}$ leads to two unit-Gaussian distributions for which an $N$-sample has a divergence of $\Delta$.
But why focus only on location shifts? Instead, what if the true model $P$ had mean zero but differing variance? Then
$$
KL_N\left(N\left(0,\sigma^2\right), N\left(0,1\right)\right) = \frac{N}{2}\left(\sigma^2 - 1 - \log \sigma^2\right).
$$
To obtain a KL divergence of $\Delta$, set $\sigma^2 = - W\left(-\exp\left[-\left(\frac{2\Delta}{N} + 1\right)\right]\right)$, where $W()$ is the [Lambert W function](https://en.wikipedia.org/wiki/Lambert_W_function).
Neither of these calibration models seems especially salient in the meta-analysis context.
There, it seems rather implausible that a model would be badly mis-specified in terms of the center or spread of the distribution of effects.
I would find it more plausible that a model might be mis-specified in terms of its shape somehow---perhaps missing some skewness or kurtosis in the distribution of effects that a (Gaussian) random effects model would not capture.
What if the superior model was a shifted gamma distribution with shape parameter $\alpha$, scale parameter $\beta$, and shift $c$? If $X + c \sim \Gamma(\alpha, \beta)$, then setting $\beta = \sqrt{\alpha}$ and $c = -\sqrt{\alpha}$ leads to a distribution with mean zero, unit variance, skewness $2 / \sqrt{\alpha}$, and differential entropy
$$
\mathcal{H}(P) = \alpha - \frac{1}{2}\log (\alpha) + \log \Gamma(\alpha) + (1 - \alpha) \psi(\alpha),
$$
where $\Gamma()$ is the [gamma function](https://en.wikipedia.org/wiki/Gamma_function) and $\psi()$ is the [digamma function](https://en.wikipedia.org/wiki/Digamma_function).
The KL divergence of an $N$-sample from this shifted gamma distribution to the unit normal distribution is therefore
$$
\begin{aligned}
KL_N\left(\Gamma\left(\alpha,\sqrt\alpha\right) + \sqrt\alpha, N\left(0,1\right)\right) &= \frac{1}{2}\log(2\pi) + \frac{1}{2} - \mathcal{H}(P) \\
&= \frac{1}{2}\log(2\pi) + \frac{1}{2} - \alpha + \frac{1}{2}\log (\alpha) - \log \Gamma(\alpha) - (1 - \alpha) \psi(\alpha).
\end{aligned}
$$
Although there's not a closed-form inverse, it's easy enough to solve numerically for $\alpha$.
If mis-specification could manifest as skewness, then why not kurtosis too?
Suppose the true model is a Student $t$ distribution with degrees of freedom $\nu$.
For simplicity, I won't worry about forcing the $t$ to have unit variance.
For an $N$-sample, the KL-divergence of the standard normal from a $t(\nu)$ distribution is
$$
\begin{aligned}
KL_N\left(t(\nu), N\left(0,1\right)\right) &= \frac{N}{2}\left[\log(2\pi) + \frac{\nu}{\nu - 2} \right. \\
& \qquad \qquad \left. - (\nu + 1)\left(\psi\left(\frac{\nu + 1}{2}\right) - \psi\left(\frac{\nu}{2}\right)\right) \right. \\
& \qquad \qquad \qquad \qquad \qquad\left. - \log(\nu) - 2 \log B\left(\frac{\nu}{2}, \frac{1}{2}\right) \right],
\end{aligned}
$$
where $\psi()$ is the [digamma function](https://en.wikipedia.org/wiki/Digamma_function) and $B()$ is the [beta function](https://en.wikipedia.org/wiki/Beta_function).
For a given $\Delta$, the expression can be solved numerically for $\nu$.
The interactive @fig-normal displays all four of these calibration models, where the superior distributions differ by location, scale, skew, or kurtosis, respectively.
One thing that jumps out to me about this visualization is that the divergence from the $t(\nu)$ distribution seems much subtler than the others.
With this distribution, all of the action is in the tails.
Put differently, the KL-divergence is driven by differences in the probability assigned to rare extreme outcomes.
::::: {.grid .column-page}
:::: {.g-col-9 .center}
```{ojs}
//| echo: false
normal = require( 'https://cdn.jsdelivr.net/gh/stdlib-js/stats-base-dists-normal@umd/browser.js' )
student = require( 'https://cdn.jsdelivr.net/gh/stdlib-js/stats-base-dists-t@umd/browser.js' )
gamma = require( 'https://cdn.jsdelivr.net/gh/stdlib-js/stats-base-dists-gamma@umd/browser.js' )
trigamma = require( 'https://cdn.jsdelivr.net/gh/stdlib-js/math-base-special-trigamma@umd/browser.js' )
// Location-shift
mu = math.sqrt(2 * Delta / N)
// Scale-shift
function KL_sigma(s) {
return N * (s**2 - 1 - 2 * math.log(s)) / 2;
}
function dKL_sigma(s) {
return (s - 1 / s) * N;
}
sigma = newtonRaphson(KL_sigma, dKL_sigma, Delta, 1.1, 1, 100, 1e-5, 100)
// Skewness
function logKL_gamma(la) {
let a = math.exp(la);
let c1 = (math.log(2 * math.pi) + 1) / 2;
let c2 = gamma.entropy(a, math.sqrt(a));
return math.log(N) + math.log(c1 - c2);
}
function dlogKL_gamma(la) {
let a = math.exp(la);
const c1 = (math.log(2 * math.pi) + 1) / 2;
let c2 = gamma.entropy(a, math.sqrt(a));
let c3 = 1 / (2 * a) - 1;
let c4 = (1 - a) * trigamma(a);
return a * (c3 - c4) / (c1 - c2);
}
log_Delta = math.log(Delta)
log_alpha = newtonRaphson(logKL_gamma, dlogKL_gamma, log_Delta, -log_Delta, -10, 10, 1e-5, 100)
alpha = math.exp(log_alpha)
// Kurtosis
function KL_t(v) {
let c1 = math.log(2 * math.pi) / 2;
let c2 = v / (2 * (v - 2));
let c3 = student.entropy(v);
return N * (c1 + c2 - c3);
}
function dKL_t(v) {
let c1 = (v - 2)**(-2);
let c2 = 1 / (2 * v);
let c3 = (trigamma((v + 1) / 2) - trigamma(v / 2)) * (v + 1) / 4;
return -N * (c1 + c2 + c3);
}
nu = newtonRaphson(KL_t, dKL_t, Delta, 3, 2, 100, 1e-5, 100)
// Compute densities
pts2 = 201
dat3 = Array(pts2).fill().map((element, index) => {
let x = - scale + index * 2 * scale / (pts2 - 1);
let p = normal.pdf(x,0,1);
let q = normal.pdf(x, mu, 1);
let r = normal.pdf(x,0, sigma);
let s = gamma.pdf(x + math.sqrt(alpha), alpha, math.sqrt(alpha));
let t = student.pdf(x, nu);
return ({
x: x,
p: p,
q: q,
r: r,
s: s,
t: t
})
})
colormap = [{v:"q",c:"blue"},{v:"r",c:"green"},{v:"s",c:"orange"},{v:"t",c:"red"}]
function density_plot(v, c) {
return Plot.plot({
height: 200,
width: 400,
y: {
grid: false,
axis: false,
domain: [0, 0.45]
},
x: {
label: null
},
marks: [
Plot.lineY(dat3, {x: "x", y: "p", stroke: "black"}),
Plot.areaY(dat3, {x: "x", y: v, fill: c, fillOpacity: 0.5}),
Plot.lineY(dat3, {x: "x", y: v, stroke: c})
]
});
}
plot_list = Object.values(colormap).map(x => density_plot(x.v, x.c))
```
```{ojs}
//| echo: false
//| label: fig-normal
//| fig-cap: "Standard normal distribution (black line) that diverges by $\\Delta$ from selected non-normal distributions, including $N(\\mu, 1)$ (blue density), $N(0,\\sigma^2)$ (green density), shifted $\\Gamma(\\alpha, \\sqrt{\\alpha})$ (orange density), or $t(\\nu)$ (red density)."
html`<div style="display: grid; grid-template-columns: 1fr 1fr; column-gap: 5px; row-gap: 5px;">
${plot_list}
</div>`
```
::::
:::: {.g-col-3}
```{ojs}
mu_toprint = mu.toFixed(3)
D_mu = N * mu**2 / 2
D_mu_toprint = D_mu.toFixed(3)
sigma_toprint = sigma.toFixed(3)
D_sigma = KL_sigma(sigma).toFixed(3)
alpha_toprint = alpha.toFixed(3)
D_alpha = math.exp(logKL_gamma(log_alpha)).toFixed(3)
nu_toprint = nu.toFixed(3)
D_nu = KL_t(nu).toFixed(3)
```
```{ojs}
//| panel: input
//| echo: false
Inputs.bind(
Inputs.range(
[1, Nmax],
{step: 1, label: tex`N`}
),
viewof N
)
Inputs.bind(
Inputs.range(
[0, Delta_max],
{step: Delta_step, label: tex`\Delta`}
),
viewof Delta
)
viewof scale = Inputs.range(
[1., 10],
{value: 5, step: 0.5, label: "Scale"}
)
tex`\mu = ${mu_toprint}`
tex`\sigma = ${sigma_toprint}`
tex`\alpha = ${alpha_toprint}`
tex`\nu = ${nu_toprint}`
```
::::
:::::
How useful are these calibration models? I'm not entirely sure yet.[^comment]
I find them appealing because they turn the very abstract concept of KL-divergence into something that can be visualized.
However, having to make arbitrary choices of distributions to construct the calibration model does seem a bit weird to me.
I will just leave this up and let it marinate a bit.
[^comment]: If you have thoughts on these models, please feel free to comment!