Different GLM result for Python and R with binomial regression

Question

I have the following R code with binomial regression to fit the y and polynomial of x

res = glm(df.mat ~ poly(x, deg=degree), family=binomial(link="logit"))

and the result is

R code result

However, when I use the statsmodels GLM function in Python as

import statsmodels.api as sm

def poly(x, d):
    x = np.array(x)
    X = np.transpose(np.vstack((x**k for k in range(d+1))))
    return np.linalg.qr(X)[0][:,1:]

res_m = sm.GLM(np.array(df_mat), poly(x, 7), family = sm.families.Binomial())
result = res_m.fit()

with the polynomial function created by Python equivalent to R poly() function? then the result is

Python result

Either the difference arises in the poly calculation or the GLM fitting (or both). Have you compared the result of R's poly(x, deg = degree) with Python's poly(x, 7) to see if the issue is there? If you need more help than that, you should probably post some sample data. — Gregor Thomas
keeping zero power in polynomial would keep the constant. return np.linalg.qr(X)[0] — Josef

merv merv · Accepted Answer · 2020-10-06T16:18:23

According to the documentation, the statsmodels GLM does not automatically include an intercept term, which the R model does. Try including a column of ones in the matrix you pass as exog to the sm.GLM. The documentation notes a convenience function exists for this very purpose: statsmodels.tools.add_constant.

Different GLM result for Python and R with binomial regression

1 Answers