python GLM Poisson regression diverges in hand-coded implementation

Question

I'm having trouble with convergence of this R-code translated to Python: statsmodels.GLM gives correct result

import numpy as np
import pandas as pd
from scipy.special import expit   # overflow encountered in exp
import statsmodels.api as sm
from matplotlib import pyplot as plt

# data
y= np.array([2,3,6,7,8,9,10,12,15],  dtype = np.float)
x= np.array([1,1,1,1,1,1,1,1,1,-1, -1, 0, 0, 0, 0, 1, 1, 1],  dtype = np.float).reshape(2,9).T    # with intercept in 1st col
data= pd.DataFrame({"y": y, "intercept": x[:,0], "x": x[:,1]})
print(data)

# for glm function
x1= np.array([-1, -1, 0, 0, 0, 0, 1, 1, 1],  dtype = np.float)
_x1 = sm.add_constant(x1, prepend=False)
_y = y[:, np.newaxis]

########################
## GLM
#######################
# implement-poisson-regression
from statsmodels.genmod.cov_struct import (Exchangeable, Independence,Autoregressive)
from statsmodels.genmod.families import Poisson

ind= ind = Independence()
fam= sm.families.links.log()    ## Log-Linear Regression, also known as Poisson Regression

##data.endog DV, data.exog IV
poisson_model = sm.GLM(  _y, _x1, cov_struct=ind, family = Poisson(link= fam))
poisson_results = poisson_model.fit( atol=1e-8)     # attach_wls=True,
print(poisson_results.summary())
##print(poisson_results.params)

###########
# https://www.statsmodels.org/devel/examples/notebooks/generated/glm.html
# Plot yhat vs y:
from statsmodels.graphics.api import abline_plot

yhat = poisson_results.mu
fig, ax = plt.subplots()
ax.scatter(yhat, _y)
line_fit = sm.OLS(_y, sm.add_constant(yhat, prepend=True)).fit()
abline_plot(model_results=line_fit, ax=ax)

ax.set_title('Model Fit Plot')
ax.set_ylabel('Observed values')
ax.set_xlabel('Fitted values');
plt.show()

sklearn gives different betas - is the problem in log_scaling_preprocessing because of the nature of lambda param ? still cannot achieve correct coef=0.6698 & intercept=1.8893

#######################
## SKLEARN
#######################
from sklearn import linear_model

regr = linear_model.PoissonRegressor(tol= 1e-8)
regr.fit(data[[ "x"]], data.y)

print("\nPoissonRegressor")
print(regr.score(data[[ "x"]], data.y))
print("coef_ ", regr.coef_)
print("intercept_ ", regr.intercept_)

but the main problem from the above-linked article is hand-coded Newton-Raphson algorithm , as it seems to diverge, not converge, though being coded like in linked R-code -- where is my mistake in translating the logics to Python ?? Whether w_matrix should be calculated in the other way ?

def poisson_Newton_Raphson(x,y,b_init,tol=1e-8):
  change = np.inf
  b_old = b_init
  while(change > tol) :
    eta = x @ b_old  # linear predictor
    w = np.diag(expit(eta))
    temp = x.T @ w @ x         # ?? np.linalg.inv(x.T @ w @ x)
    b_new = b_old + (temp @ x.T @ (y-expit(eta)))
    change= np.linalg.norm(b_old - b_new, 2)
##    print(change)
    b_old = b_new
    print(b_old)
  return b_new

res= poisson_Newton_Raphson(x, y, np.zeros(2))
print(res)

how to make all 3 codes give the same results ?? (correct in GLM_example)

P.S. alternative -- anyway, my w-matrix seems to be errorness - how to correct it ?

P.P.S. Poisson Regression to model either counts or frequencies, e.g. for a non-constant λ, Poisson distribution has its variance equal to its mean, as lambda grows large the Poisson looks more and more like a normal distribution, to interpret

P.P.P.S. glm-Course_001, n.b. choice of distribution

Answer 1

All the methods you listed give exactly the same results. Here is a display:

x = np.array([ 1.,-1,1,-1,1,0,1,0,1,0,1,0,1,1,1,1,1,1]).reshape(-1,2)
y = np.array([2.0, 3.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 15.0])
 
def poisson_Newton_Raphson(x, y, b_init, tol = 1e-8):
   change =  1
   b = np.array(b_init) + 0.0
   while change > tol:
     eta = x @ b  # linear predictor
     w = np.diag(np.exp(eta))
     direction = np.linalg.solve(x.T @ w @ x, x.T @ (y - np.exp(eta)))
     change = np.linalg.norm(direction)
     b += direction 
   return b

poisson_Newton_Raphson(x,y,[1,1])
array([1.889272 , 0.6697856])

---

from sklearn.linear_model import PoissonRegressor
regr = PoissonRegressor(fit_intercept = False, tol = 1e-8, alpha = 0)
regr.fit(x, y).coef_
array([1.88927199, 0.66978561])

---

from statsmodels.discrete.discrete_model import Poisson
Poisson(y, x).fit(disp = 0).params
array([1.889272 , 0.6697856])

---

def poisson_IRLS(x, y, b_init, tol=1e-8):
  change = np.Inf
  b = np.array(b_init) + 0.0
  while change > tol:
    eta = x @ b
    weight = np.exp(eta)
    w = np.diag(weight)
    z = np.linalg.solve(w, y - weight)
    direction = np.linalg.solve( x.T @ w @ x, x.T @ w @ z)
    change = np.linalg.norm(direction)
    b += direction
  return b

poisson_IRLS(x,y,[1,1])
array([1.889272 , 0.6697856])

Answer 2

another hand-coded Newton-Raphson implementation with scipy.optimize.minimize - was not working, but corrections advised above by Onyambu did it working as well

import numpy as np
from scipy.optimize import minimize
import scipy.special

# data
x = np.array([ 1.,-1,1,-1,1,0,1,0,1,0,1,0,1,1,1,1,1,1]).reshape(-1,2)
y = np.array([2.0, 3.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 15.0])

########################

# Define the log-likelihood function
def loglik(b, X, y):
    expb = np.exp(np.dot(X, b))
    # print(expb)
    ret = np.sum(y * np.log(expb) - expb - np.log(scipy.special.factorial(y)))
    return ret

# Define the gradient of the log-likelihood function
def gradloglik(b, X, y):
    expb = np.exp(np.dot(X, b))
    ret = np.dot(X.T, (y - expb))
    # print(r)
    return ret

# Define the optimization problem
def opt_prob(b, X, y):
    return -1* loglik(b, X, y), -1* gradloglik(b, X, y)

# Set the initial values for the coefficients
##print(X.shape[1])
b0 = np.zeros(x.shape[1])

# Optimize the log-likelihood function
result = minimize(opt_prob, b0, args=(x, y), method='BFGS', tol=1.4e-08, jac=True)


# Extract the optimized coefficients
coef = result.x
print(coef)