使用 Iris 数据集使用 Python 在 R 中重现 LASSO / Logistic 回归结果

Posted 2023-03-12

【中文标题】使用 Iris 数据集使用 Python 在 R 中重现 LASSO / Logistic 回归结果

【英文标题】：Reproducing LASSO / Logistic Regression results in R with Python using the Iris Dataset

【发布时间】：2017-09-20 18:24:03

【问题描述】：

我正在尝试在 Python 中重现以下 R 结果。在这种特殊情况下，R 预测技能低于 Python 技能，但根据我的经验，这通常不是这种情况（因此希望在 Python 中重现结果），所以请忽略这里的细节。

目的是预测花卉种类（'versicolor' 0 或 'virginica' 1）。我们有 100 个标记样本，每个样本由 4 个花特征组成：萼片长度、萼片宽度、花瓣长度、花瓣宽度。我将数据分为训练集（60% 的数据）和测试集（40% 的数据）。对训练集进行 10 折交叉验证，寻找最优的 lambda（优化的参数在 scikit-learn 中为“C”）。

我在 R 中使用 glmnet，alpha 设置为 1（用于 LASSO 惩罚），对于 python，scikit-learn 的 LogisticRegressionCV 函数与“liblinear”求解器（唯一可以与L1 处罚）。交叉验证中使用的评分指标在两种语言之间是相同的。然而不知何故，模型结果是不同的（为每个特征找到的截距和系数差异很大）。

R 代码

library(glmnet)
library(datasets)
data(iris)

y <- as.numeric(iris[,5])
X <- iris[y!=1, 1:4]
y <- y[y!=1]-2

n_sample = NROW(X)

w = .6
X_train = X[0:(w * n_sample),]  # (60, 4)
y_train = y[0:(w * n_sample)]   # (60,)
X_test = X[((w * n_sample)+1):n_sample,]  # (40, 4)
y_test = y[((w * n_sample)+1):n_sample]   # (40,)

# set alpha=1 for LASSO and alpha=0 for ridge regression
# use class for logistic regression
set.seed(0)
model_lambda <- cv.glmnet(as.matrix(X_train), as.factor(y_train),
                        nfolds = 10, alpha=1, family="binomial", type.measure="class")

best_s  <- model_lambda$lambda.1se
pred <- as.numeric(predict(model_lambda, newx=as.matrix(X_test), type="class" , s=best_s))

# best lambda
print(best_s)
# 0.04136537

# fraction correct
print(sum(y_test==pred)/NROW(pred))   
# 0.75

# model coefficients
print(coef(model_lambda, s=best_s))
#(Intercept)  -14.680479
#Sepal.Length   0        
#Sepal.Width   0
#Petal.Length   1.181747
#Petal.Width    4.592025

Python 代码

from sklearn import datasets
from sklearn.linear_model import LogisticRegressionCV
from sklearn.preprocessing import StandardScaler
import numpy as np

iris = datasets.load_iris()
X = iris.data
y = iris.target
X = X[y != 0]  # four features. Disregard one of the 3 species.                                                                                                                 
y = y[y != 0]-1  # two species: 'versicolor' (0), 'virginica' (1). Disregard one of the 3 species.                                                                               

n_sample = len(X)

w = .6
X_train = X[:int(w * n_sample)]  # (60, 4)
y_train = y[:int(w * n_sample)]  # (60,)
X_test = X[int(w * n_sample):]  # (40, 4)
y_test = y[int(w * n_sample):]  # (40,)

X_train_fit = StandardScaler().fit(X_train)
X_train_transformed = X_train_fit.transform(X_train)

clf = LogisticRegressionCV(n_jobs=2, penalty='l1', solver='liblinear', cv=10, scoring = ‘accuracy’, random_state=0)
clf.fit(X_train_transformed, y_train)

print clf.score(X_train_fit.transform(X_test), y_test)  # score is 0.775
print clf.intercept_  #-1.83569557
print clf.coef_  # [ 0,  0, 0.65930981, 1.17808155] (sepal length, sepal width, petal length, petal width)
print clf.C_  # optimal lambda: 0.35938137

【参考方案1】：

上面的例子有几点不同：

系数的尺度

glmnet (https://cran.r-project.org/web/packages/glmnet/glmnet.pdf) 对数据进行标准化，并且“系数始终以原始比例返回”。因此，您在调用 glmnet 之前没有扩展数据。 Python 代码将数据标准化，然后适合该标准化数据。在这种情况下，系数是标准化比例，而不是原始比例。这使得示例之间的系数不可比。

LogisticRegressionCV 默认使用分层折叠。 glmnet 使用 k-fold。

它们拟合不同的方程。请注意，scikit-learn 逻辑适合 (http://scikit-learn.org/stable/modules/linear_model.html#logistic-regression) 与逻辑方面的正则化。 glmnet 将正则化置于惩罚之上。

选择正则化强度来尝试 - glmnet 默认为 100 个 lambdas 来尝试。 scikit LogisticRegressionCV 默认为 10。由于 scikit 求解方程，范围在 1e-4 和 1e4 (http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegressionCV.html#sklearn.linear_model.LogisticRegressionCV) 之间。

容忍度不同。在我遇到的一些问题中，收紧公差会显着改变系数。

glmnet 默认 thresh 为 1e-7 Lo​​gisticRegressionCV 默认 tol 为 1e-4 即使使它们相同，它们也可能不会测量相同的东西。我不知道什么是 liblinear 措施。 glmnet - “每个内部坐标下降循环都会继续，直到任何系数更新后目标的最大变化小于 thresh 乘以零偏差。”

您可能想尝试打印正则化路径以查看它们是否非常相似，只是停止在不同的强度上。然后你可以研究原因。

即使更改了您可以更改的内容（并非上述所有内容），您也可能不会获得相同的系数或结果。尽管您在不同的软件中解决相同的问题，但软件解决问题的方式可能会有所不同。我们看到了不同的尺度、不同的方程、不同的默认值、不同的求解器等。

【参考方案2】：

您在这里遇到的问题是数据集的排序（注意我没有检查 R 代码，但我确定这是问题所在）。如果我运行你的代码然后运行它

print np.bincount(y_train) # [50 10]
print np.bincount(y_test) # [ 0 40]

您可以看到训练集不代表测试集。但是，如果我对您的 Python 代码进行一些更改，那么我的测试准确度为 0.9。

from sklearn import datasets
from sklearn import preprocessing
from sklearn import model_selection
from sklearn.linear_model import LogisticRegressionCV
from sklearn.preprocessing import StandardScaler
import numpy as np

iris = datasets.load_iris()
X = iris.data
y = iris.target
X = X[y != 0]  # four features. Disregard one of the 3 species.                                                                                                                 
y = y[y != 0]-1  # two species: 'versicolor' (0), 'virginica' (1). Disregard one of the 3 species.                                                                               

X_train, X_test, y_train, y_test = model_selection.train_test_split(X, y, 
                                                                    test_size=0.4,
                                                                    random_state=42,
                                                                    stratify=y)


X_train_fit = StandardScaler().fit(X_train)
X_train_transformed = X_train_fit.transform(X_train)

clf = LogisticRegressionCV(n_jobs=2, penalty='l1', solver='liblinear', cv=10, scoring = 'accuracy', random_state=0)
clf.fit(X_train_transformed, y_train)

print clf.score(X_train_fit.transform(X_test), y_test)  # score is 0.9
print clf.intercept_  #0.
print clf.coef_  # [ 0., 0. ,0., 0.30066888] (sepal length, sepal width, petal length, petal width)
print clf.C_ # [ 0.04641589]

【讨论】：

非常感谢。 train_test_split 函数似乎很方便，但是（请参阅我对 Grr 的回复）我不确定这是否是两种语言之间存在差异的原因。我将尝试在两者之间实现平衡拆分（在 R 和 python 中），然后更新我的初始帖子。

我建议创建两个文件，一个用于训练集，一个用于测试集，然后将它们读入 Python 和 R。这是确保正确拆分数据的最安全方法。

【参考方案3】：

我不得不对这里的几件事感到不满。

首先，“对于 python，scikit-learn 的 LogisticRegressionCV 函数带有“liblinear”求解器（唯一可以与 L1 惩罚一起使用的求解器）”。这显然是错误的，除非您打算以更明确的方式对其进行限定。只需查看sklearn.linear_model 类的描述，您就会看到少数特别提到 L1。我相信其他人也允许你实现它，但我真的不想数它们。

其次，您拆分数据的方法不太理想。查看拆分后的输入和输出，您会发现在您的拆分中，所有测试样本的目标值都是 1，而目标值 1 仅占训练样本的 1/6。这种不代表目标分布的不平衡将导致您的模型拟合不佳。例如，只使用开箱即用的sklearn.model_selection.train_test_split，然后完全按照您的方式重新安装LogisticRegressionCV 分类器，结果是.92 的准确度

现在说了这么多，有一个glmnet package for python，你可以使用这个包复制你的结果。该项目的作者有一篇博客讨论了尝试使用 sklearn 重新创建 glmnet 结果的一些限制。具体来说：

“Scikit-Learn 有几个求解器，类似于 glmnet、ElasticNetCV 和 LogisticRegressionCV，但它们有一些局限性。第一个只适用于线性回归，后者不处理弹性网络惩罚。” - 比尔拉特纳GLMNET FOR PYTHON

【讨论】：

感谢您的宝贵时间。我应该说“使用 LogisticRegressionCV 函数时唯一可以与 L1 惩罚一起使用的求解器”。文档列出了四个可以使用的求解器（'newton-cg'、'lbfgs'、'liblinear'、'sag'）；只有 liblinear 可以与 L1 一起使用。是的，分裂并不理想。我不会在操作中这样做；但是，由于我在 R 和 Python 之间以相同的方式拆分，我不确定这是导致不同结果的原因（我不确定如何在 R 中进行平衡拆分）。 python 的 glmnet 包可能是解决方案。谢谢。