为什么statsmodels和R的logistic回归结果不同？

Question

我试图比较python的statsmodels和R中的logistic回归实现

Python版本：

import statsmodels.api as sm
import pandas as pd
import pylab as pl
import numpy as np
df = pd.read_csv("http://www.ats.ucla.edu/stat/data/binary.csv")
df.columns = list(df.columns)[:3] + ["prestige"]
# df.hist()
# pl.show()
dummy_ranks = pd.get_dummies(df["prestige"], prefix="prestige")
cols_to_keep = ["admit", "gre", "gpa"]
data = df[cols_to_keep].join(dummy_ranks.ix[:, "prestige_2":])
data["intercept"] = 1.0
train_cols = data.columns[1:]
logit = sm.Logit(data["admit"], data[train_cols])
result = logit.fit()
result.summary2()

结果：

^{pr2}$

R版本：

data = read.csv("http://www.ats.ucla.edu/stat/data/binary.csv", head=T)
require(reshape2)
data1 = dcast(data, admit + gre + gpa ~ rank)
require(dplyr)
names(data1)[4:7] = paste("rank", 1:4, sep="")
data1 = data1[, -4]
summary(glm(admit ~ gre + gpa + rank2 + rank3 + rank4, family=binomial, data=data1))

结果：

Call:
glm(formula = admit ~ gre + gpa + rank2 + rank3 + rank4, family = binomial,
    data = data1)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-1.5133  -0.8661  -0.6573   1.1808   2.0629

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.184029   1.162421  -3.599 0.000319 ***
gre          0.002358   0.001112   2.121 0.033954 *
gpa          0.770591   0.343908   2.241 0.025046 *
rank2       -0.369711   0.310342  -1.191 0.233535
rank3       -1.015012   0.335147  -3.029 0.002457 **
rank4       -1.249251   0.414416  -3.014 0.002574 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 466.13  on 377  degrees of freedom
Residual deviance: 434.12  on 372  degrees of freedom
AIC: 446.12

Number of Fisher Scoring iterations: 4

结果相差很大，例如，秩2的p值分别为0.03和0.2。我想知道造成这种差异的原因是什么？请注意，我为这两个版本创建了虚拟变量，为python版本创建了一个常量列，它在R中自动处理

而且，python似乎快了2倍：

##################################################
# python timing
def test():
    for i in range(5000):
        logit = sm.Logit(data["admit"], data[train_cols])
        result = logit.fit(disp=0)
import time
start = time.time()
test()
print(time.time() - start)
10.099738836288452
##################################################
# R timing
> f = function() for(i in 1:5000) {mod = glm(admit ~ gre + gpa + rank2 + rank3 + rank4, family=binomial, data=data1)}
> system.time(f())
   user  system elapsed
 17.505   0.021  17.526

Answer 1

我把R部分改成这样：

makeDummy = function(x, x1) { ifelse(is.na(x), NA, ifelse(x == x1, 1, 0)) }
data = read.csv("http://www.ats.ucla.edu/stat/data/binary.csv", head=T)
data$rank2 = makeDummy(data$rank, 2)
data$rank3 = makeDummy(data$rank, 3)
data$rank4 = makeDummy(data$rank, 4)
summary(glm(admit ~ gre + gpa + rank2 + rank3 + rank4, family=binomial, data=data))

结果与

^{pr2}$

我想要么我用错了dplyr::dcast，要么dcast出了问题。在

Answer 2

不确定您的数据操作意图是什么，但它们似乎在R运行中丢失了信息。如果我保留了所有的排名信息，那么我会在原始数据对象上得到这些信息（结果在它们重叠的区域看起来非常相似）。（可能性仅估计为任意常数，因此只能比较对数似然的差异。即使有这样的警告，偏差也应该是负对数可能性的两倍，因此这些结果也具有可比性。）

> summary(glm(admit ~ gre + gpa +as.factor( rank), family=binomial,
       data=data))  # notice that I'm using your original data-object

Call:
glm(formula = admit ~ gre + gpa + as.factor(rank), family = binomial, 
    data = data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6268  -0.8662  -0.6388   1.1490   2.0790  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -3.989979   1.139951  -3.500 0.000465 ***
gre               0.002264   0.001094   2.070 0.038465 *  
gpa               0.804038   0.331819   2.423 0.015388 *  
as.factor(rank)2 -0.675443   0.316490  -2.134 0.032829 *  
as.factor(rank)3 -1.340204   0.345306  -3.881 0.000104 ***
as.factor(rank)4 -1.551464   0.417832  -3.713 0.000205 ***
 -
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 499.98  on 399  degrees of freedom
Residual deviance: 458.52  on 394  degrees of freedom
AIC: 470.52

Number of Fisher Scoring iterations: 4

Answer 3

我只能回答，而不能在已接受的答案中添加评论。在python中，通常需要删除其中一个伪类，使其成为引用类，但我不认为您需要为R这样做，因为glm将为您设置引用类。基本上，如果我能正确理解你的代码，你就不需要这行了。。。在

data1 = data1[, -4]

试着把声望放在原样上，但要用作为因素（）首先。在