为何Python恒等连接泊松GLM收敛但R版本不收敛？

Great question—this discrepancy boils down to key differences in how R's base glm() and Python's statsmodels GLM handle initialization, numerical stability, and convergence checks for Poisson models using an identity link. Let's break this down:

1. Initialization Strategy Differences #

R's glm() uses iterative reweighted least squares (IRLS) for fitting, and its default initial values for Poisson models with an identity link can be problematic for your dataset. Specifically, it starts with a naive initial guess (often based on a linear regression fit) that may lead to negative predicted means (mu) or situations where log(y/mu) produces NaNs (like when y=0 and mu is positive, leading to log(0) which is -Inf).

Python's statsmodels, by contrast, uses a more robust initialization approach for GLMs. For Poisson identity link models, it often first fits a model with a log link (which is numerically safer) and uses those coefficients as starting values for the identity link model. This avoids the unstable initial steps that trip up R's glm().

2. Handling of Zero Counts & Numerical Edge Cases #

Your dataset includes y=0 values, which are tricky for Poisson identity link models because the model enforces mu = β*M (where mu must be positive for the Poisson distribution to be valid). When R's glm() iterates, it may generate mu values that conflict with y (e.g., positive mu paired with y=0), leading to NaN in the log-likelihood calculation and triggering an early termination.

Statsmodels includes safeguards to handle these edge cases: it may clip extreme values, adjust the iteration step size, or use numerical smoothing to avoid NaNs during optimization, allowing the algorithm to keep converging even when temporary numerical issues arise.

3. Convergence Checks & Optimizer Tolerances #

R's glm() has strict default convergence criteria. If it encounters a single iteration with NaNs or invalid coefficients, it immediately stops with an error. Statsmodels, on the other hand, uses more lenient (or adaptive) convergence checks and may continue iterating even after minor numerical blips, as long as the coefficients are moving toward a stable solution.

Fixing the Convergence Issue in R #

You can resolve this in R by addressing the initialization problem directly:

Option 1: Use the Python-derived Coefficient as a Starting Value #

Since Python successfully converged to a coefficient of ~1578382.66, pass this as an initial value to R's glm():

init_coef <- 1578382.6635341388
model <- glm(b ~ 0 + M, family = poisson(identity), start = init_coef)
summary(model)

Option 2: Adjust Zero Values to Avoid Log(0) Errors #

Add a tiny positive value to your response variable to eliminate y=0 cases, which removes the source of NaNs in the log calculation:

# Add a minimal offset to zero values
b_adjusted <- pmax(b, 1e-6)
model <- glm(b_adjusted ~ 0 + M, family = poisson(identity))
summary(model)

Option 3: Use a More Robust GLM Implementation #

The glm2 package in R is designed to handle problematic convergence cases better than base glm():

install.packages("glm2")
library(glm2)
model <- glm2(b ~ 0 + M, family = poisson(identity))
summary(model)

内容的提问来源于stack exchange，提问作者librus