python – 为什么statsmodels不能重现我的R逻辑回归结果?

前端之家收集整理的这篇文章主要介绍了python – 为什么statsmodels不能重现我的R逻辑回归结果?前端之家小编觉得挺不错的,现在分享给大家,也给大家做个参考。

我很困惑为什么R和statsmodels中的逻辑回归模型不一致.

如果我在R中准备一些数据

# From https://courses.edx.org/c4x/MITx/15.071x/asset/census.csv
library(caTools) # for sample.split
census = read.csv("census.csv")
set.seed(2000)
split = sample.split(census$over50k,SplitRatio = 0.6)
censusTrain = subset(census,split==TRUE)
censusTest = subset(census,split==FALSE)

然后运行逻辑回归

CensusLog1 = glm(over50k ~.,data=censusTrain,family=binomial)

我看到results就像

                                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)                              -8.658e+00  1.379e+00  -6.279 3.41e-10 ***
age                                       2.548e-02  2.139e-03  11.916  < 2e-16 ***
workclass Federal-gov                     1.105e+00  2.014e-01   5.489 4.03e-08 ***
workclass Local-gov                       3.675e-01  1.821e-01   2.018 0.043641 *  
workclass Never-worked                   -1.283e+01  8.453e+02  -0.015 0.987885    
workclass Private                         6.012e-01  1.626e-01   3.698 0.000218 ***
workclass Self-emp-inc                    7.575e-01  1.950e-01   3.884 0.000103 ***
workclass Self-emp-not-inc                1.855e-01  1.774e-01   1.046 0.295646    
workclass State-gov                       4.012e-01  1.961e-01   2.046 0.040728 *  
workclass Without-pay                    -1.395e+01  6.597e+02  -0.021 0.983134   
...

但我在Python中使用相同的数据,首先从R导出

write.csv(censusTrain,file="traincensus.csv")
write.csv(censusTest,file="testcensus.csv")

然后导入到Python中

import pandas as pd

census = pd.read_csv("census.csv")
census_train = pd.read_csv("traincensus.csv")
census_test = pd.read_csv("testcensus.csv")

我得到的错误和奇怪的结果与我在R中得到的结果没有任何关系.

如果我只是尝试

import statsmodels.api as sm

census_log_1 = sm.Logit.from_formula(f,census_train).fit()

我收到一个错误

ValueError: operands could not be broadcast together with shapes (19187,2) (19187,) 

即使用patsy准备数据也是如此

import patsy
f = 'over50k ~ ' + ' + '.join(list(census.columns)[:-1])
y,X = patsy.dmatrices(f,census_train,return_type='dataframe')

census_log_1 = sm.Logit(y,X).fit()

导致相同的错误.我可以避免错误的唯一方法是使用GLM

census_log_1 = sm.GLM(y,X,family=sm.families.Binomial()).fit()

但是这会产生results,与完全不同于(我认为是)等效的R API产生的那些:

                                                   coef    std err          t      P>|t|      [95.0% Conf. Int.]
----------------------------------------------------------------------------------------------------------------
Intercept                                       10.6766      5.985      1.784      0.074        -1.055    22.408
age                                             -0.0255      0.002    -11.916      0.000        -0.030    -0.021
workclass[T. Federal-gov]                       -0.9775      4.498     -0.217      0.828        -9.794     7.839
workclass[T. Local-gov]                         -0.2395      4.498     -0.053      0.958        -9.055     8.576
workclass[T. Never-worked]                       8.8346    114.394      0.077      0.938      -215.374   233.043
workclass[T. Private]                           -0.4732      4.497     -0.105      0.916        -9.288     8.341
workclass[T. Self-emp-inc]                      -0.6296      4.498     -0.140      0.889        -9.446     8.187
workclass[T. Self-emp-not-inc]                  -0.0576      4.498     -0.013      0.990        -8.873     8.758
workclass[T. State-gov]                         -0.2733      4.498     -0.061      0.952        -9.090     8.544
workclass[T. Without-pay]                       10.0745     85.048      0.118      0.906      -156.616   176.765
...

为什么Python中的逻辑回归会产生错误,并且会产生与R生成错误不同的结果?这些API实际上是不相同的(我之前已经让它们工作以产生相同的结果)?是否需要对数据集进行一些额外的处理才能使它们被statsmodel使用?

最佳答案
错误是由于patsy将LHS变量扩展为完全治疗对比的事实.如文档字符串中所示,Logit不处理此问题,但正如您所见,GLM与二项式族一样.

如果没有完整的输出,我无法说出结果的差异.很可能它是对分类变量的不同默认处理,或者您使用的是不同的变量.并非所有都列在您的输出中.

您可以通过执行以下预处理步骤来使用logit.

census = census.replace(to_replace={'over50k' : {' <=50K' : 0,' >50K' : 1}})

另请注意,logit的默认解算器似乎不能很好地解决此问题.它遇到了奇异的矩阵问题.实际上,这个问题的条件数是巨大的,你在R中得到的可能不是一个完全融合的模型.您可以尝试减少虚拟变量的数量.

[~/]
[73]: np.linalg.cond(mod.exog)
[73]: 4.5139498536894682e+17

我必须使用以下内容来获得解决方

mod = sm.formula.logit(f,data=census)
res = mod.fit(method='bfgs',maxiter=1000)    

你的一些细胞最终变得非常小.这与其他稀疏虚拟变量相结合.

[~/]
[81]: pd.Categorical(census.occupation).describe()
[81]: 
                    counts     freqs
levels                              
?                    1816  0.056789
Adm-clerical         3721  0.116361
Armed-Forces            9  0.000281
Craft-repair         4030  0.126024
Exec-managerial      3992  0.124836
Farming-fishing       989  0.030928
Handlers-cleaners    1350  0.042217
Machine-op-inspct    1966  0.061480
Other-service        3212  0.100444
Priv-house-serv       143  0.004472
Prof-specialty       4038  0.126274
Protective-serv       644  0.020139
Sales                3584  0.112077
Tech-support          912  0.028520
Transport-moving     1572  0.049159

猜你在找的Python相关文章