Stable Deviance Contributions Near Saturation

Why Deviance Needs Special Care Near Saturation

A fitted model can predict a mean extremely close to the observed count. In that situation, the deviance contribution should be small: the fitted model and the saturated model make nearly the same prediction for that observation. But the usual formulas subtract logarithmic terms whose leading parts nearly cancel. Floating-point arithmetic can therefore lose much of the small difference that remains.

The remedy is to preserve the small relative change rather than form and subtract nearly equal quantities. The function log1p(x) computes $\log(1+x)$ accurately when $x$ is close to zero. A further rearrangement expresses the deviance as a sum of nonnegative terms; for very small arguments, those terms can be evaluated with a convergent series.

Deviance Contributions and Stable Functions

Definition: For a binomial observation with $m$ trials, $y$ successes, and fitted mean $\mu=mp$, assume $0<y<m$ and $0<\mu<m$. Its deviance contribution is $$ d_{\mathrm{B}}(y,\mu)=2\left[y\log\left(\frac{y}{\mu}\right)+(m-y)\log\left(\frac{m-y}{m-\mu}\right)\right]. $$ For a Poisson observation with $y>0$ and fitted mean $\mu>0$, the contribution is $$ d_{\mathrm{P}}(y,\mu)=2\left[y\log\left(\frac{y}{\mu}\right)-(y-\mu)\right]. $$ In both cases, the saturated mean is $y$, the mean that reproduces the observed count exactly.

For the binomial case, set $t=\mu-y$, the difference between the fitted and saturated means. The two logarithms in the half-deviance become $-\log1p(t/y)$ and $-\log1p(-t/(m-y))$, respectively. Their first-order terms cancel. The Poisson expression has the same issue: writing $x=t/y$, its half-deviance is $y\{x-\log1p(x)\}$.

Definition: For $x>-1$, define $$ h(x)=x-\log1p(x). $$ For $|x|<1$, its convergent power series is $$ h(x)=\sum_{k=2}^{\infty}\frac{(-1)^k x^k}{k} =\frac{x^2}{2}-\frac{x^3}{3}+\frac{x^4}{4}-\cdots. $$ For small $|x|$, evaluate this series rather than subtracting $\log1p(x)$ from $x$. Otherwise, the direct definition using $\log1p$ is generally suitable.

Worked Example: A Binomial Mean Just Above Saturation

Suppose $m=100$, $y=40$, and the fitted mean is $\mu=40.1$. Then $t=0.1$, $y=40$, and $m-y=60$. The half-deviance can be evaluated as $$ \frac{d_{\mathrm{B}}}{2} =40h(0.1/40)+60h(-0.1/60). $$

Using the series for $h$, the combined term of order $k$ is $$ \frac{t^k}{k}\left(\frac{(-1)^k}{40^{k-1}}+\frac{1}{60^{k-1}}\right). $$ The terms for $k=2$ through $k=5$ are, respectively, $$ 0.0002083333333333,\quad -0.0000001157407407,\quad 0.0000000005063657,\quad -0.0000000000006269. $$ The $k=6$ term is approximately $1.84\times10^{-15}$, and subsequent terms are smaller. Summing the full series gives a half-deviance of approximately $0.0002082180983332$, so $$ d_{\mathrm{B}}\approx 0.0004164361967. $$ The fourth-order term changes the displayed precision; stopping after the cubic term would not justify reporting these digits.

Both series arguments are small, and their stable contributions are added rather than obtained by subtracting two nearly cancelling first-order quantities. The positive fitted-mean difference also gives $h(0.0025)$ and $h(-1/600)$, both nonnegative.

Worked Example: A Poisson Mean Just Above Saturation

Let $y=20$ and $\mu=20.02$, so $t=0.02$ and $x=t/y=0.001$. The stable half-deviance is $$ \frac{d_{\mathrm{P}}}{2}=20h(0.001). $$

The series gives $$ h(0.001)=\frac{0.001^2}{2}-\frac{0.001^3}{3} +\frac{0.001^4}{4}-\frac{0.001^5}{5}+\cdots $$ $$ =5.0000000000\times10^{-7} -3.3333333333\times10^{-10} +2.5000000000\times10^{-13} -2.0000000000\times10^{-16}+\cdots. $$ Thus $h(0.001)\approx4.9966691647\times10^{-7}$, giving $$ \frac{d_{\mathrm{P}}}{2}\approx9.9933383293\times10^{-6}, \qquad d_{\mathrm{P}}\approx1.9986676659\times10^{-5}. $$ The leading quadratic term alone would give $d_{\mathrm{P}}\approx2.0\times10^{-5}$. It indicates the scale, but the additional series terms are needed for the more precise value shown.

Deriving the Stable Forms

The algebraic identities below justify the stable expressions and explain why they remain informative near saturation. They do not change the definition of deviance; they rearrange it so that small contributions are not computed as differences of large, nearly equal terms.

Result: With $t=\mu-y$, $q=m-y$, and the interior conditions in the definition, $$ \frac{d_{\mathrm{B}}(y,\mu)}{2} =yh\left(\frac{t}{y}\right)+qh\left(-\frac{t}{q}\right), \qquad \frac{d_{\mathrm{P}}(y,\mu)}{2} =yh\left(\frac{t}{y}\right). $$ Each expression is nonnegative and equals zero at saturation, $\mu=y$.

Proof. For the binomial observation, $\mu=y+t$ and $m-\mu=q-t$. Substitution in the definition gives $$ \frac{d_{\mathrm{B}}}{2} =-y\log1p(t/y)-q\log1p(-t/q). $$ The identity $h(x)=x-\log1p(x)$ implies $-\log1p(x)=-x+h(x)$. Apply it to both terms. Their linear parts are $-y(t/y)=-t$ and $-q(-t/q)=t$, which cancel, leaving $yh(t/y)+qh(-t/q)$. For Poisson, substitution into its definition gives $$ \frac{d_{\mathrm{P}}}{2} =y[-\log1p(t/y)]+t =y\left[-\log1p(t/y)+t/y\right] =yh(t/y). $$ Finally, $h(0)=0$, and $h(x)\geq0$ for $x>-1$: since $\log(1+x)\leq x$ on that domain, $x-\log(1+x)\geq0$. The binomial arguments are greater than $-1$ because $0<\mu<m$; the Poisson argument is greater than $-1$ because $\mu>0$. This proves the identities and nonnegativity. $\square$

The series also gives a precise local approximation. For fixed interior $y,m$, with $t$ sufficiently close to zero, substituting the series for $h$ yields $$ \frac{d_{\mathrm{B}}}{2} =\frac{t^2}{2}\left(\frac{1}{y}+\frac{1}{m-y}\right)+O(t^3), \qquad \frac{d_{\mathrm{P}}}{2} =\frac{t^2}{2y}+O(t^3). $$ Indeed, the $k=2$ terms in the convergent series give the displayed quadratic terms; all remaining terms have powers at least three in $t$. This local expansion describes the scale of the deviance, but it is not a replacement for evaluating the full contribution when greater accuracy is needed.

Worked Example: Preserving a Tiny Poisson Contribution

Take $y=1{,}000{,}000$ and $\mu=1{,}000{,}000.01$, treating the displayed inputs as exact for the calculation. Then $t=0.01$ and $x=10^{-8}$. Directly computing $t-y\log(\mu/y)$ subtracts two numbers close to $0.01$. In the stable form, $$ \frac{d_{\mathrm{P}}}{2}=10^6 h(10^{-8}). $$

The first terms give $$ 10^6 h(10^{-8}) =10^6\left(\frac{10^{-16}}{2} -\frac{10^{-24}}{3} +\frac{10^{-32}}{4}-\cdots\right) \approx4.9999999667\times10^{-11}. $$ Therefore $d_{\mathrm{P}}\approx9.9999999333\times10^{-11}$. The direct expression asks floating-point arithmetic to recover a quantity around $10^{-10}$ from terms around $10^{-2}$. The series instead calculates directly at the scale of the answer.

In actual computation, the stored value of $\mu$ may differ slightly from the decimal value shown. The stable formula prevents avoidable loss after the inputs are represented; it cannot restore precision already lost while calculating or storing an input mean.

Practical Evaluation and Common Mistakes

A reliable implementation first forms $t=\mu-y$, then uses the relative changes in the stable identities. When $|x|$ is small, use a series for $h(x)$; when it is not small, compute $x-\log1p(x)$. The series has a quantitative stopping check: after retaining terms through $k=K$, for $|x|<1$, the absolute tail is bounded by $$ \sum_{k=K+1}^{\infty}\frac{|x|^k}{k} \leq\frac{|x|^{K+1}}{(K+1)(1-|x|)}. $$ Choose $K$ large enough that this bound is below the desired absolute error. For positive $x$, the series alternates; for negative $x$, its terms are positive. In either case, evaluate the small-argument function rather than relying on a difference of nearly equal numbers.

Using $\log(y/\mu)$ directly near equality. The ratio may round close to one before the logarithm is taken. Compute the relative change and use $\log1p$ instead.
Assuming that $\log1p$ alone solves every cancellation problem. In the Poisson formula, $x-\log1p(x)$ still subtracts nearly equal quantities. Use $h(x)$ with its series for small $|x|$.
Dropping higher-order terms while reporting excessive digits. The fourth-order term can affect the rounded result, as the binomial calculation shows. Bound or calculate the omitted tail before choosing the reported precision.
Interpreting the quadratic approximation as the exact deviance. It gives the leading local behavior, not the full contribution at an arbitrary fitted mean.

These expressions matter because deviance is used to monitor how a fitted model approaches the observed data. If near-saturated contributions are rounded to zero or computed inaccurately, objective comparisons and convergence diagnostics can be misleading even when the underlying fitted means are good. Stable evaluation makes a small deviance meaningful rather than merely a by-product of cancellation.

Check Your Understanding

Use the stable identities to explain both the algebra and the numerical choices.

For a binomial observation, what are the two relative-change arguments in the stable half-deviance when $t=\mu-y$?
Why can the expression $x-\log1p(x)$ still be inaccurate when $|x|$ is very small?
State the leading quadratic term of the Poisson deviance contribution near $\mu=y$.
For $y=30$, $m=80$, and a fitted mean $\mu=30.001$, write the stable binomial half-deviance in terms of $h$, without evaluating it.
Why is it unsafe to report many digits from a local series after omitting its fourth-order term?

Stable Deviance Contributions Near Saturation

What You'll Learn