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

Question

我对为什么我在 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 实际上不是等效的吗(我之前让它们工作以产生相同的结果)？是否需要对数据集进行一些额外的处理才能使其可供 statsmodels 使用？

Answer 1

错误是由于 patsy 将 LHS 变量扩展为完整的 Treatement 对比。 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