为什么 glmnet 的系数估计在具有相同输入参数的模型之间变化很大？答案

【问题标题】：Why the coefficient estimates of glmnet varies a lot between models with same input parameters?为什么 glmnet 的系数估计在具有相同输入参数的模型之间变化很大？
【发布时间】：2017-12-13 23:00:14
【问题描述】：

我一直在尝试使用cv.glmnet 拟合套索模型。我尝试基于标准化实现四种不同的模型（3 个使用cv.glmnet，1 个使用caret::train）。所有四个模型都给出了非常不同的系数估计值，我不知道为什么。

这是一个完全可重现的代码：

library("glmnet")
data(iris)
iris <- iris
dat <- iris[iris$Species %in% c("setosa","versicolor"),]
X <- as.matrix(dat[,1:4])
Y <- as.factor(as.character(dat$Species))

set.seed(123)
model1 <- cv.glmnet(x = X,
                    y = Y,
                    family = "binomial",
                    standardize = FALSE,
                    alpha = 1,
                    lambda = rev(seq(0,1,length=100)),
                    nfolds = 3)

set.seed(123)
model2 <- cv.glmnet(x = scale(X, center = T, scale = T),
                    y = Y,
                    family = "binomial",
                    standardize = FALSE,
                    alpha = 1,
                    lambda = rev(seq(0,1,length=100)),
                    nfolds = 3)
set.seed(123)
model3 <- cv.glmnet(x = X,
                    y = Y,
                    family = "binomial",
                    standardize = TRUE,
                    alpha = 1,
                    lambda = rev(seq(0,1,length=100)),
                    nfolds = 3)

##Using caret
library("caret")

lambda.grid <- rev(seq(0,1,length=100)) #set of lambda values for cross-validation
alpha.grid <- 1 #alpha
trainControl <- trainControl(method ="cv",
                             number=3) #3-fold cross-validation
tuneGrid <- expand.grid(.alpha=alpha.grid, .lambda=lambda.grid) #these are tuning parameters to be passed into the train function below

set.seed(123)
model4 <- train(x = X,
                y = Y,
                method="glmnet",
                family="binomial",
                standardize = FALSE,
                trControl = trainControl,                          
                tuneGrid = tuneGrid)

c1 <- coef(model1, s=model1$lambda.min)
c2 <- coef(model2, s=model2$lambda.min)
c3 <- coef(model3, s=model3$lambda.min)
c4 <- coef(model4$finalModel, s=model4$finalModel$lambdaOpt)
c1 <- as.matrix(c1)
c2 <- as.matrix(c2)
c3 <- as.matrix(c3)
c4 <- as.matrix(c4)

model2 预先缩放自变量（向量X），model3 通过设置standardize = TRUE 来做到这一点。所以至少这两个模型应该返回相同的结果 - 但事实并非如此。

四个模型得到的lambda.min分别是：

model1 = 0

model2 = 0

model3 = 0

model4 = 0.6565657

模型之间的系数估计也有很大差异。为什么会发生这种情况？

【问题讨论】：

glmnet 的标准化是由下面的 fortran 代码完成的，所以很难判断它和 scale 是否真的 100% 做同样的事情。
无论使用何种底层编程语言，规模都应该标准化数据。这意味着用相应的列均值减去每列值，并将列标准差除以具有单位方差和零均值。不太明白为什么事情不应该如此复杂:-(
用于比较 c2 和 c3：在 ?glmnet standardize 参数处；当 TRUE... “系数始终以原始比例返回。，当您手动转换时不会发生这种情况。因此您可以手动转换回原始比例：xs = scale(X) ; sx = attr(xs, "scaled:scale") ; ce = attr(xs, "scaled:center") ; co = as.numeric(c2) ; co[-1] / sx ; co[1] - sum((co[-1] / sx)*sx) - 这是更近了
我没有测试过你说的。你的评论对我来说很有意义。感谢那。但现在的问题是，为什么模型 1 和模型 3 会输出不同的系数估计值？模型 1 不标准化数据，模型 3 可以，但根据文档，系数以原始比例返回。无法将这些发现与 glmnet 文档所说的相关联

标签： r r-caret glmnet

【解决方案1】：

其实scale(x) & standardize = FALSE 和x & standardize = TRUE 有一点不同。我们需要多个 (N-1)/N。

见here。

如果我们使用高斯族，

library(glmnet)
X <- matrix(runif(100, 0, 1), ncol=2)
y <- 1 -2*X[,1] + X[,2]

enet <- glmnet(X, y, lambda=0.1,standardize = T,family="gaussian")
coefficients(enet)
coef <- coefficients(enet)
coef[2]*sd(X[,1])/sd(y) #standardized coef
#[1] -0.6895065

enet1 <- glmnet(scale(X)/99*100, y/(99/100*sd(y)),lambda=0.1/(99/100*sd(y)),standardize = F,family="gaussian")
coefficients(enet1)[2]
#[1] -0.6894995

如果我们使用二项式家庭，

data(iris)
iris <- iris
dat <- iris[iris$Species %in% c("setosa","versicolor"),]
X <- as.matrix(dat[,1:4])
Y <- as.factor(as.character(dat$Species))

set.seed(123)
model1 <- cv.glmnet(x = X,
                y = Y,
                family = "binomial",
                standardize = T,
                alpha = 1,
                lambda = rev(seq(0,1,length=100)),
                nfolds = 3)
coefficients(model1,s=0.03)[3]*sd(X[,2])
#[1] -0.3374946

set.seed(123)
model3 <- cv.glmnet(x = scale(X)/99*100,
                y = Y,
                family = "binomial",
                standardize = F,
                alpha = 1,
                lambda = rev(seq(0,1,length=100)),
                nfolds = 3)
coefficients(model3,s=0.03)[3]
#[1] -0.3355027

这些结果几乎相同。希望这个答案还为时不晚。

【讨论】：