4
$\begingroup$

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?

$\endgroup$
11
  • 1
    $\begingroup$ 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. $\endgroup$ Apr 6, 2014 at 15:26
  • 3
    $\begingroup$ 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. $\endgroup$
    – Dason
    Apr 6, 2014 at 18:36
  • 1
    $\begingroup$ 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. $\endgroup$
    – Dason
    Apr 6, 2014 at 18:49
  • 1
    $\begingroup$ @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. $\endgroup$
    – orome
    Apr 6, 2014 at 18:55
  • 1
    $\begingroup$ You can force an intercept by including a column of 1s in your predictors. $\endgroup$
    – Dason
    Apr 6, 2014 at 19:03

2 Answers 2

7
$\begingroup$

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", 
                   header=TRUE)
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

cuse = pd.read_table("http://data.princeton.edu/wws509/datasets/cuse.dat",
                     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.

$\endgroup$
6
  • $\begingroup$ I'm working on getting permission to post the data; not looking good. $\endgroup$
    – orome
    Apr 6, 2014 at 18:35
  • $\begingroup$ "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. $\endgroup$
    – Josef
    Apr 6, 2014 at 18:40
  • $\begingroup$ 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? $\endgroup$
    – orome
    Apr 6, 2014 at 18:53
  • $\begingroup$ 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 $\endgroup$
    – orome
    Apr 6, 2014 at 19:26
  • $\begingroup$ 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? $\endgroup$
    – orome
    Apr 6, 2014 at 19:45
4
$\begingroup$

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.

$\endgroup$
9
  • $\begingroup$ Accuracy of 0.802816901 (R) vs 0.745042493 (Python). $\endgroup$
    – orome
    Apr 6, 2014 at 15:46
  • $\begingroup$ 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). $\endgroup$
    – orome
    Apr 6, 2014 at 15:48
  • $\begingroup$ 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. $\endgroup$ Apr 6, 2014 at 15:50
  • $\begingroup$ 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?). $\endgroup$
    – orome
    Apr 6, 2014 at 16:27
  • $\begingroup$ 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. $\endgroup$ Apr 6, 2014 at 16:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.