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

Question

我一直在尝试使用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

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

- j1897

---

glmnet 的标准化是由底层的 Fortran 代码完成的，因此很难判断它和 scale 是否实际上完全做了相同的事情。 - JAD

无论使用哪种编程语言，规模应标准化数据。这意味着通过相应的列均值减去每个列值并将列标准差除以单位方差和零均值来实现。不太明白为什么事情应该如此复杂，当它本不应该是这样的 :-( - j1897

对于c2到c3的比较：在?glmnet的standardize参数中，当为TRUE时... *"The coefficients are always returned on the original scale.*，但是当您手动转换时，这将不会发生。因此，您可以手动将其转换回原始比例：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) - 这些更接近。 - user20650

我还没有测试你所说的内容。你的评论对我来说很有道理，谢谢。但现在的问题是为什么model1和model3输出不同的系数估计值？model1没有标准化数据，而model3则进行了标准化，但根据文档，系数是以原始比例返回的。无法将这些发现与glmnet的文档相一致。 - j1897

Answer 1

实际上，在 scale(x) & standardize = FALSE 和 x & standardize = TRUE 之间有一点不同。我们需要乘以 (N-1)/N。

请参见此处。

如果我们使用高斯族，

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

这些结果几乎相同。希望这个答案不算太迟。

- ming gao