# Should statsmodels's GLM produce the same results as R's lm?

Should Python's

statsmodels.api.GLM(train_y, train_X, family=sm.families.Binomial()).fit().predict(test_X)


always produce the same results as R's

predict(glm(y ~ ., data=train_X, family=binomial), newdata=test)


where train_y is a pandas DataFrame containing the y column in the corresponding R data.frame, train; and where test_X and train_X are dataframes containing the remaining columns from the test and train dataframes respectively?

If not, are there parameters that I can supply to statsmodels's GLM to make it produce the same results as R's glm?

• This question strikes me as being more about the underlying theory in statistics, rather than about how to program R & Python. I believe it is better suited to Cross Validated (stats.SE) than here. – gung - Reinstate Monica Apr 6 '14 at 15:26
• One thing to double check is that they are indeed fitting the same model. I haven't used the python package you mention but I know some packages don't incorporate an intercept term by default. R does include an intercept term by default. I haven't used python's version so I can't say if that is what might be causing differences but it's something to double check. – Dason Apr 6 '14 at 18:36
• I don't know. I don't even know if that's an issue. I just know that some packages don't default to adding an intercept for you so it's always good to check if you are indeed fitting the same model. Hopefully somebody else that knows statsmodels will come along and comment either telling you to ignore me (because you are fitting the same models) or how to fit an intercept term with your model. – Dason Apr 6 '14 at 18:49
• @Dason: Whether it's the cause, I do think it's a factor (no pun): I don't see an intercept in the Python params. – orome Apr 6 '14 at 18:55
• You can force an intercept by including a column of 1s in your predictors. – Dason Apr 6 '14 at 19:03

Yes, they should give the same answers if you fit the same model. Compare

R code

cuse <- read.table("http://data.princeton.edu/wws509/datasets/cuse.dat",
attach(cuse)
mod <- glm(cbind(using, notUsing) ~ age + education + wantsMore , family= binomial)
summary(mod)


Python code

import pandas as pd
import statsmodels.api as sm

sep=" +")
res = sm.formula.glm("using + notUsing ~ C(age, Treatment('<25')) + "
"education + wantsMore",  family=sm.families.Binomial(),
data=cuse).fit()
res.summary()


If there's a convergence issue here, I wouldn't trust either answer without knowing why there are convergence issues. I'd be interested to have a look at some data that can reproduce these convergence failures in R.

• I'm working on getting permission to post the data; not looking good. – orome Apr 6 '14 at 18:35
• "Should" applies to regular cases, up to convergence tolerance. However, in "irregular" cases like near singular design matrix or perfect separation, the difference in implementation will show up. But then the statistical result is ambiguous anyway. – Josef Apr 6 '14 at 18:40
• Could it be that Python is "succeeding" where R is "failing": R gives me a useless model with every one of 300 variables highly significant, while Python gives me one, that matches the first branch of a CART model on the same data? – orome Apr 6 '14 at 18:53
• Also, I can't use formula: I have a column called 'break' that always causes File "<string>", line 1 break ^ SyntaxError: unexpected EOF while parsing – orome Apr 6 '14 at 19:26
• And sm.GLM(train.Y, train[columns], family=sm.families.Binomial()).fit() gives me different results than sm.formula.glm('Y ~ ' + '+'.join(columns), data=tweets_train, family=sm.families.Binomial()).fit(). Shouldn't they be the same? – orome Apr 6 '14 at 19:45

I would guess that they shouldn't always produce identical results, although the difference in the results would almost always be negligible. When trying to fit non-linear GLiMs (such as a logistic regression), there is no closed form solution for estimating the betas as there is in OLS regression. Instead, a search algorithm is used. Typically this is the Newton-Raphson gradient descent method. As far as I know, there is no absolute guarantee that this will yield identical results, and small differences in the way it's implemented in different languages behind the scenes could cause small differences in outputs. However, unless there is some strong ambiguity in the data, I would suspect any differences would not show up for several decimal places; that is, well after the point where people would have rounded anyway.

• Accuracy of 0.802816901 (R) vs 0.745042493 (Python). – orome Apr 6 '14 at 15:46
• FWIW, If I start from scratch and build and preprocess my data in Python (split the same as in R, but otherwise managed from the start in Python) I get 0.742209632. I'd expect some subtle differences here (I'm tokenizing text); but the big difference between what R and Python do with the same data (same tokenization) is surprising (to me). – orome Apr 6 '14 at 15:48
• Hmmm... that's a larger difference than I would have expected. If the sole issue is in the estimation of the betas, it would imply that R's estimated betas are better optimized. It should also be possible to input different starting points to the optimization in both, & to request different optimization algorithms in both to see if that helps. – gung - Reinstate Monica Apr 6 '14 at 15:50
• Actually, looking more closely, here's what I see: Nearly all of the R model's (more than 300!) p-values are highly significant (<2e-16), while only one of the Python model's is (at 0.025236), and that one corresponds to the first branch of the CART tree fitted to the same data. Also, R produces "does not converge" warnings. I'm not sure what all that means (I'm new to this) but it seems like Python is doing a lot better job of coming up with a usable model (is R overfitting?). – orome Apr 6 '14 at 16:27
• Is your Python model a CART? Your R model is logistic regression. You shouldn't expect that these would given identical responses b/c they are different models. – gung - Reinstate Monica Apr 6 '14 at 16:33