Debugging Extreme-Predictor Fits

Why Trace an Extreme-Predictor Fit?

An initialized coefficient vector can be finite while its linear predictors are so large in magnitude that ordinary probability calculations fail. In a binomial model, a computed probability may round to exactly zero or one, even when its complementary probability is mathematically positive. In a Poisson model with a log link, exponentiating a large predictor may overflow. Those failures can then appear in the likelihood, score, or working weights and make an iteration look more conclusive than it is.

Debugging means following the same observation through each calculation, rather than treating a failed fit as one indivisible event. First check the initialized coefficients and the resulting predictor. Then check the mean and its probability-scale representation, the log-likelihood contribution, and the working weight. The key distinction is between a mathematically extreme quantity and a quantity that has become inaccurate or unrepresentable in finite-precision arithmetic.

Trace the Binomial Calculations

For observation \(i\), let \(Y_i\) be a binomial response with \(m_i\) trials and success probability \(p_i\). Under the logit link, the predictor is \(\eta_i=x_i^\mathsf{T}\beta\), and \(p_i=1/(1+\exp(-\eta_i))\). The binomial log-likelihood contribution, including its data-only combinatorial term, is

\ell_i(\beta)=\log {m_i \choose y_i}+y_i\log(p_i)+(m_i-y_i)\log(1-p_i).

When \(\eta_i\) is extreme, computing \(p_i\) and then subtracting it from one can discard important information. A stable calculation works with the two probabilities separately and computes the log-likelihood on the log scale.

Definition: Define the softplus function \(s(t)=\log(1+\exp(t))\). Its stable piecewise form is \(s(t)=t+\log(1+\exp(-t))\) for \(t>0\), and \(s(t)=\log(1+\exp(t))\) for \(t\leq 0\). For a logit predictor, compute \(\log(p_i)=-s(-\eta_i)\) and \(\log(1-p_i)=-s(\eta_i)\). These expressions avoid forming a rounded probability and then taking its logarithm.

In software, the \(\log(1+\exp(t))\) terms should be evaluated with a stable log-one-plus-exponential operation, often called log1p. For positive \(\eta_i\), the two probabilities can also be evaluated using different expressions: \(p_i=1/(1+\exp(-\eta_i))\) and \(1-p_i=\exp(-\eta_i)/(1+\exp(-\eta_i))\). The second expression preserves a small complement that may be lost by computing \(1-p_i\) from a rounded \(p_i\).

Result: For a binomial observation with \(m_i\) trials, the log-likelihood can be written as \(\ell_i=\log {m_i\choose y_i}+y_i\eta_i-m_i s(\eta_i)\). Its first and second derivatives with respect to \(\eta_i\) are \(y_i-m_i p_i\) and \(-m_i p_i(1-p_i)\), respectively.

To verify the result, the logit equation gives \(\log(p_i)=\eta_i-s(\eta_i)\) and \(\log(1-p_i)=-s(\eta_i)\). Substituting both identities into the binomial log-likelihood gives the stated expression. Also, \(s'(t)=\exp(t)/(1+\exp(t))\), which equals \(p_i\) at \(t=\eta_i\). Differentiating once therefore gives \(y_i-m_i p_i\); differentiating again, using \(dp_i/d\eta_i=p_i(1-p_i)\), gives \(-m_i p_i(1-p_i)\). \(\square\)

The binomial working weight for this predictor is \(w_i=m_i p_i(1-p_i)\), the negative second derivative just obtained. In an extreme tail this weight is genuinely small, but it should not be set to zero merely because a rounded probability has become exactly zero or one. For example, calculate \(p_i\) and \(1-p_i\) separately, then multiply them. An exactly zero weight produced by subtracting rounded values can falsely make the information contribution disappear.

Worked Example: A Positive Extreme Binomial Predictor

Consider one Bernoulli observation with \(m=1\), \(y=0\), predictor row \(x=(1)\), and initialized coefficient \(\beta=40\). The coefficient is finite, but the resulting predictor is \(\eta=x^\mathsf{T}\beta=40\).

The success probability is \(p=1/(1+\exp(-40))\), approximately \(1-4.2484\times10^{-18}\). In ordinary double-precision arithmetic, this can round to exactly 1. Subtracting that stored value from one then gives zero, so directly evaluating \(y\log(p)+(1-y)\log(1-p)\) may produce \(\log(0)=-\infty\).

The stable calculation uses \(s(40)=40+\log(1+\exp(-40))\), approximately 40, and the log-likelihood formula gives \(\ell=0\cdot40-s(40)\), approximately \(-40\). More precisely, the contribution is \(-\log(1+\exp(40))\), which is finite. The complementary probability is approximately \(4.2484\times10^{-18}\), so the likelihood is small but not zero.

The working weight is \(p(1-p)\), approximately \(4.2484\times10^{-18}\). Computing \(1-p\) after \(p\) has rounded to one instead gives a false weight of zero. The score is \(y-p\), approximately \(-1\); thus a tiny weight does not imply a tiny score. These are different quantities and should be checked separately.

Worked Example: A Negative Extreme Binomial Predictor

Now take \(m=1\), \(y=1\), \(x=(1)\), and initialized \(\beta=-40\), so \(\eta=-40\). Here \(p=\exp(-40)/(1+\exp(-40))\), approximately \(4.2484\times10^{-18}\), while \(1-p\) is close to one.

Since the observed outcome is a success, the log-likelihood contribution is \(\log(p)\). The stable expression gives \(\log(p)=-s(40)\), approximately \(-40\), rather than requiring the logarithm of a probability that may underflow for still more negative predictors. The likelihood itself is approximately \(4.2484\times10^{-18}\), but its log is a manageable finite number.

The working weight is again \(p(1-p)\), approximately \(4.2484\times10^{-18}\). It is positive mathematically, although it may underflow to zero at more extreme values. Record whether the zero arose in the probability calculation or in the multiplication: that detail helps distinguish a flawed computation from a value below the representable range. The score is \(y-p\), approximately 1.

Trace the Poisson Calculations

For a Poisson response \(Y_i\) with mean \(\mu_i\), the log link gives \(\mu_i=\exp(\eta_i)\). For a nonnegative integer observation \(y_i\), the probability mass and log-likelihood contribution are

P(Y_i=y_i)=\frac{\exp(-\mu_i)\mu_i^{y_i}}{y_i!}, \qquad \ell_i=y_i\eta_i-\exp(\eta_i)-\log(y_i!).

For computing a likelihood contribution, the log form is usually the useful one. Computing the probability mass first can underflow to zero even when the log-likelihood is finite. The working weight for a Poisson log-link observation is \(\mu_i\), so the same large predictor that threatens the likelihood calculation can also make the weight enormous.

Result: Under the Poisson log link, differentiating \(\ell_i=y_i\eta_i-\exp(\eta_i)-\log(y_i!)\) with respect to \(\eta_i\) gives the score contribution \(y_i-\mu_i\). The second derivative is \(-\mu_i\), so the working weight is \(\mu_i=\exp(\eta_i)\).

Indeed, \(\log(y_i!)\) is constant in \(\eta_i\), and the derivative of \(\exp(\eta_i)\) is \(\exp(\eta_i)=\mu_i\). A second differentiation gives \(-\mu_i\). Thus the log-likelihood, score, and weight all depend on the same exponential, but their numerical behavior differs: the log-likelihood may be finite when a probability mass underflows, while the mean and weight may eventually overflow. \(\square\)

Worked Example: A Large but Representable Poisson Mean

Suppose a Poisson observation has \(y=0\), \(x=(1)\), and initialized \(\beta=40\). Then \(\eta=40\) and \(\mu=\exp(40)\), approximately \(2.3539\times10^{17}\). This mean is extremely large, but it is representable in double-precision arithmetic.

The probability of observing zero is \(\exp(-\mu)\), too small to represent as a positive double-precision number; a direct probability calculation may return zero. Taking the logarithm of that computed zero gives \(-\infty\), which is not the correct finite log-likelihood. The log-likelihood formula instead gives \(\ell=0\cdot40-\exp(40)-\log(0!)=-2.3539\times10^{17}\), because \(0!=1\) and \(\log(1)=0\).

The score is \(y-\mu=-2.3539\times10^{17}\), and the working weight is \(\mu=2.3539\times10^{17}\). These are large but representable. This trace shows why a zero probability returned by exponentiation does not establish that the log-likelihood is negative infinity.

Worked Example: When the Poisson Mean Overflows

Keep \(y=0\) and \(x=(1)\), but let the initialized coefficient be \(\beta=710\), so \(\eta=710\). In double precision, \(\exp(710)\) exceeds the largest finite representable number and is returned as positive infinity.

The mathematical mean is finite, and the mathematical log-likelihood is \(-\exp(710)\), but neither value can be stored as a finite double. The computed weight \(\mu\) is also infinite. This is not the same failure as the previous example: there, the probability underflowed while the mean and log-likelihood remained representable; here, the mean itself is outside the numerical range.

Do not replace the infinite weight with an arbitrary cap and then treat the result as an ordinary fit. Flag the offending observation and predictor, check the initialized coefficients and predictor construction, and use a numerical method with an appropriate range if the model calculation must continue. Capping changes the likelihood and its derivatives.

A Practical Debugging Sequence

An extreme-predictor diagnosis is more informative when it records the calculation that first becomes unreliable. Apply the following sequence to each flagged observation, retaining the observation index and model family in any diagnostic output.

Check initialization and predictor construction.
Record the initialized coefficients, design row, offset if present, and resulting \(\eta_i\). A finite coefficient vector does not guarantee a moderate predictor; also check for overflow or accidental scaling in the dot product.

Check the mean and its components.
For a logit binomial model, compute both \(p_i\) and \(1-p_i\) safely. For a Poisson log-link model, check whether \(\exp(\eta_i)\) is finite before using it as a mean or weight.

Check the log-likelihood without exponentiating it.
Use the stable binomial softplus expression or the Poisson log-likelihood expression. Do not compute a tiny probability mass first and then take its logarithm.

Check score and working weight separately.
Compare the values with their mathematical expressions, and note whether a weight is genuinely tiny, has underflowed, or has overflowed. A score and its corresponding weight need not have similar magnitudes.

This sequence also helps locate whether the problem began before or after the likelihood evaluation. For example, a predictor that is unexpectedly extreme may come from an initialization or predictor-construction issue; a correct extreme predictor may instead expose an unstable probability calculation. Keep that distinction explicit when deciding what to investigate next.

Why the Distinction Matters

The likelihood determines what the model says about the data at the current coefficients, while working weights determine how strongly observations contribute to a local update. If a probability rounds to a boundary, or a positive weight is incorrectly set to zero, the numerical fit can behave as though an observation contains no curvature at all. If a Poisson mean overflows, the resulting infinite weight can dominate computations or contaminate an information matrix.

At the same time, not every extreme value is a bug. A binomial probability can be genuinely close to zero, with a genuinely small positive weight. A Poisson mean can be genuinely enormous. The aim is not to make such values look moderate; it is to compute representable quantities accurately, identify quantities outside the machine range, and avoid mistaking a failure of representation for a property of the statistical model.

Common Mistake

Wrong reasoning: “The software returned a probability of one, so the complementary probability is exactly zero; therefore a zero binomial weight and an infinite log-likelihood contribution are correct.” The conclusion does not follow. The stored probability may have rounded to one even though the model probability is strictly between zero and one.

Correction: For a logit model, compute the two log-probabilities using \(-s(-\eta_i)\) and \(-s(\eta_i)\), and compute the weight from separately evaluated \(p_i\) and \(1-p_i\) where representable. For a Poisson model, calculate the log-likelihood directly rather than taking the logarithm of an underflowed probability mass. Report overflow or underflow when the target quantity cannot be represented; do not silently substitute a boundary value.

Key takeaway: Trace each extreme-predictor observation from initialized coefficients to predictor, mean or probabilities, log-likelihood, and weight. Stable log-scale calculations prevent avoidable failures, while explicit range checks reveal when a mathematically valid quantity cannot be represented.

Check Your Understanding

For each question, identify the numerical calculation that should be checked, not only the fitted probability or mean.

In a binomial logit model with \(\eta=45\) and \(y=0\), why can computing \(\log(1-p)\) from a stored \(p\) give the wrong result? Give a stable expression for the log-likelihood contribution.
For a binomial observation with \(m=4\), explain why an exactly zero computed weight at a large positive predictor might not equal the mathematical weight.
In a Poisson log-link model, why can the probability of \(y=0\) underflow even when the log-likelihood contribution is finite?
At what point in the calculation should a Poisson fit with an infinite computed mean be flagged, and why is replacing that mean by a cap not an exact correction?
Explain why a tiny working weight does not by itself imply that the score contribution is tiny.