为什么 statsmodels 无法重现我的 R 逻辑回归结果?

2024-03-25

我很困惑为什么 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)

I see results https://gist.github.com/orome/9840196 like

                                           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 using

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

trying

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

导致同样的错误。我可以避免错误的唯一方法是使用 useGLM

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

但这会产生results https://gist.github.com/orome/9839624与(我认为的)等效 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 使用?


该错误是由于 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
本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)

为什么 statsmodels 无法重现我的 R 逻辑回归结果? 的相关文章

  • Pandas dataframe:每批行的操作

    我有一个熊猫数据框df我想计算每批行的一些统计信息 例如 假设我有一个batch size 200000 对于每批batch sizerows 我想要一列的唯一值的数量ID我的数据框 我怎样才能做这样的事情呢 这是我想要的一个例子 prin
  • PyTorch 给出 cuda 运行时错误

    我对我的代码做了一些小小的修改 以便它不使用 DataParallel and DistributedDataParallel 代码如下 import argparse import os import shutil import time
  • 读取R中打开的Excel文件

    有没有办法将打开的Excel文件读入R 当Excel中打开一个excel文件时 Excel会对文件加锁 比如R中的read方法无法访问该文件 你能绕过这个锁吗 Thanks 编辑 这发生在带有原始 Excel 的 Windows 下 发生错
  • ValueError:不支持连续[重复]

    这个问题在这里已经有答案了 我正在使用 GridSearchCV 进行线性回归的交叉验证 不是分类器也不是逻辑回归 我还使用 StandardScaler 对 X 进行标准化 我的数据框有 17 个特征 X 和 5 个目标 y 观察 约11
  • 排序因素与水平

    有人能解释一下 R 中 ordered 参数的用途吗 R says ordered逻辑标志来确定级别是否应被视为有序 按给定的顺序 所以如果我有一个名为名称的因素并设置ordered TRUE names lt factor c fred
  • 如何使用 javascript/jquery/AJAX 调用 Django REST API?

    我想使用 Javascript jQuery AJAX 在前端调用 Django Rest API 请求方法是 POST 但当我看到 API 调用它的调用 OPTIONS 方法时 所以 我开始了解access control allow o
  • 在相同任务上,Keras 比 TensorFlow 慢

    我正在使用 Python 运行斩首 DCNN 本例中为 Inception V3 来获取图像特征 我使用的是 Anaconda Py3 6 和 Windows7 使用 TensorFlow 时 我将会话保存在变量中 感谢 jdehesa 并
  • 设置 verify_certs=False 但 elasticsearch.Elasticsearch 因证书验证失败而引发 SSL 错误

    self host KibanaProxy 自我端口 443 self user 测试 self password 测试 我需要禁止证书验证 使用选项时它与curl一起使用 k在命令行上 但是 在使用 Elasticsearch pytho
  • 揭秘sharedctypes性能

    在 python 中 可以在多个进程之间共享 ctypes 对象 然而我注意到分配这些对象似乎非常昂贵 考虑以下代码 from multiprocessing import sharedctypes as sct import ctypes
  • Werkzeug 中的线程和本地代理。用法

    首先 我想确保我正确理解了功能的分配 分配本地代理功能以通过线程内的模块 包 共享变量 对象 我对吗 其次 用法对我来说仍然不清楚 也许是因为我误解了作业 我用烧瓶 如果我有两个 或更多 模块 A B 我想将对象C从模块A导入到模块B 但我
  • 如何定义“f_n-chi-square”函数并使用“uniroot”求置信区间?

    I want to get a 95 confidence interval for the following question 我已经写了函数f n在我的 R 代码中 我首先使用 Normal 随机采样 100 个样本 然后定义函数h
  • 从 python 检测 macOS 中的暗模式

    我正在编写一个 PyQt 应用程序 我必须添加一个补丁 以便在启用暗模式的 Macos 上可以读取字体 app QApplication Fix for the font colours on macos when running dark
  • Python对象初始化性能

    我只是做了一些快速的性能测试 我注意到一般情况下初始化列表比显式初始化列表慢大约四到六倍 这些可能是错误的术语 我不确定这里的行话 例如 gt gt gt import timeit gt gt gt print timeit timeit
  • 当有很多列时,使用 readr::read_csv() 导入数据时覆盖列类型

    我正在尝试使用 R 中的 readr read csv 读取 csv 文件 我导入的 csv 文件大约有 150 列 我只包含示例的前几列 我希望将第二列从默认类型 我执行 read csv 时为日期 覆盖为字符或其他日期格式 GIS Jo
  • 从列表python的单个列表中删除子列表

    我已经经历过从列表列表中删除子列表 https stackoverflow com questions 47209786 removing sublists from a list of lists 但当我为我的数据集扩展它时 它不适用于我
  • 按特定样本前缀对列名称向量进行子集化

    假设我有一个如下所示的数据框 ca01 lt c 1 10 ca02 lt c 2 11 ca03 lt c 3 12 stuff 1 lt rep test 10 other lt rep 9 10 data lt data frame
  • 如何绘制堆积比例图?

    我有一个数据框 x lt data frame id letters 1 3 val0 1 3 val1 4 6 val2 7 9 id val0 val1 val2 1 a 1 4 7 2 b 2 5 8 3 c 3 6 9 我想绘制一个
  • 计算互相关函数?

    In R 我在用ccf or acf计算成对互相关函数 以便我可以找出哪个移位给我带来最大值 从它的外观来看 R给我一个标准化的值序列 Python 的 scipy 中是否有类似的东西 或者我应该使用fft模块 目前 我正在这样做 xcor
  • python 日志记录会刷新每个日志吗?

    当我使用标准模块将日志写入文件时logging 每个日志会分别刷新到磁盘吗 例如 下面的代码会将日志刷新 10 次吗 logging basicConfig level logging DEBUG filename debug log fo
  • 使用 Python 将对象列表转为 JSON

    我在转换时遇到问题Object实例到 JSON ob Object list name scaping myObj base url u number page for ob in list name json string json du

随机推荐