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

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

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

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

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

然后运行逻辑回归

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

我看到results就像

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

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

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

然后导入到Python中

  1. import pandas as pd
  2. census = pd.read_csv("census.csv")
  3. census_train = pd.read_csv("traincensus.csv")
  4. census_test = pd.read_csv("testcensus.csv")

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

如果我只是尝试

  1. import statsmodels.api as sm
  2. census_log_1 = sm.Logit.from_formula(f,census_train).fit()

我收到一个错误

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

猜你在找的Python相关文章