为什么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)
我看到的结果像这样：
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()
但这产生的结果与我认为是等效的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中可用？
1 个回答
4
这个错误是因为patsy把左边的变量扩展成了一个完整的处理对比。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