r - glmnet:相同的 lambda 但 glmnet() 和 cv.glmnet() 的系数不同

Question

即使我使用相同的 lambda，cv.glmnet() 生成的系数似乎与 glmnet() 生成的系数不同。为什么是这样？它们不应该相同吗？

library(glmnet)

# Data dimensions
num.samples <- 30
num.genes <- 17000

# Data objects - note that both X and Y are scaled
set.seed(123)
Y <- matrix(rnorm(num.samples), ncol=1)
set.seed(1234)
X <- matrix(rnorm(num.samples*num.genes), ncol=num.genes)

# Run cv.glmnet: get lambda.min and coef
fit.cv <- cv.glmnet(X, Y, nfolds=num.samples, intercept=FALSE)
fit.cv.lambda <- fit.cv$lambda.min
fit.cv.coef <- coef(fit.cv, s = fit.cv.lambda)[,1][2:(num.genes+1)]

# Run glmnet with lambda.min from cv.glmnet: get coef
second.lambda=fit.cv.lambda-0.0001 ## second.lambda included because glmnet manual recommends using >1 lambda for glmnet()
fit <- glmnet(X, Y, lambda=c(fit.cv.lambda,second.lambda), intercept=FALSE) 
fit.lambda <- fit$lambda[1]
fit.coef <- coef(fit, s = fit.cv.lambda)[,1][2:(num.genes+1)]

# Lambda is the same, but coefficients are not
fit.cv.lambda==fit.lambda ## TRUE
not.equal = which(fit.cv.coef != fit.coef)
length(not.equal) ## 18
mean(abs(fit.cv.coef[not.equal] - fit.coef[not.equal])) ## 0.0004038209

(我还注意到 glmnet() 和 cv.glmnet() 的系数在特定的 alpha 值上没有不同，但似乎没有明显的模式)

Answer 1

简答:这是一个数值准确性问题。您遇到的差异不是由于 cv.glmnet 和 glmnet 之间的差异造成的。相反，它们是由于以下原因的组合:

• 惩罚路径 lambda 在两个对象之间是不同的，我指的是整个惩罚路径，而不仅仅是感兴趣的惩罚是否在两个对象中
• 默认收敛阈值，thresh =

如果您希望从具有不同惩罚路径的两个 glmnet 或 cv.glmnet 对象获得的估计值相等(或至少非常接近)，请在两个函数中使用 thresh = 选项来减少收敛阈值。此外，我建议在 coef() 中设置 exact = TRUE。

扩展答案:下面我们将通过几个示例来说明这一点。在这样做之前，了解 coef() 函数的逻辑(它调用 predict.glmnet() 函数与 type = "coefficients ")。

• 如果您请求对已在对象的原始惩罚路径中计算的惩罚进行系数估计，coef() 将简单地返回原始对象的估计值。

• 如果您请求对不在原始路径中的惩罚进行系数估计并且exact = FALSE(这是默认设置)，则系数将根据来自原始路径中最近的惩罚。

• 如果您请求对原始路径中不存在的惩罚项进行系数估计并且 exact = TRUE，则新惩罚项将添加到原始路径并且整个模型将重新拟合为获得估计值

示例 1:相同的惩罚路径，相同的估计值

如果使用默认参数，glmnet() 和 cv.glmnet() 将为给定的数据集计算相同的惩罚路径(计算的惩罚次数尽管由于停止标准，它们之间的路径可能不同)。我们在下面显示:

library(glmnet)

# Data dimensions
num.samples <- 30
num.genes <- 17000

# Data objects - note that both X and Y are scaled
set.seed(123)
Y <- matrix(rnorm(num.samples), ncol=1)
X <- matrix(rnorm(num.samples*num.genes), ncol=num.genes)

# Run cv.glmnet and glmnet, obtain same penalty path up to min(num penalty)
fit <- glmnet(X, Y, intercept=FALSE)
cvfit <- cv.glmnet(X, Y, intercept= FALSE)
min_num_lambdas = min(length(fit$lambda), length(cvfit$lambda))
all.equal(fit$lambda[1:min_num_lambdas], cvfit$lambda[1:min_num_lambdas]) #TRUE

然后我们可以使用 coef() 从两个对象获得相同的估计，而不管新的惩罚是否在原始路径中。

# Requested penalty in original path
coef1 <- coef(fit, s = fit$lambda[10])
coef2 <- coef(cvfit, s = fit$lambda[10])
all.equal(coef1, coef2) #TRUE

# Requested penalty not in original path -- uses interpolation
coef1 <- coef(fit, s = 0.40)
coef2 <- coef(cvfit, s = 0.40)
all.equal(coef1, coef2) #TRUE

# Force glmnet to refit the model with s added to the penalty path
coef1 <- coef(fit, s = 0.40, exact = TRUE, x = X, y = Y)
coef2 <- coef(cvfit, s = 0.40, exact = TRUE, x = X, y = Y)
all.equal(coef1, coef2) #TRUE

示例 2:不同的惩罚路径，不同的估计值

我们通过在路径中请求 99 个惩罚而不是 100 个来对 glmnet() 拟合做一个小改动。这使得 fit 和 cvfit 之间的惩罚路径不同。现在，即使我们使用 exact = TRUE 选项，估计值也会不同。

fit <- glmnet(X, Y, intercept=FALSE, nlambda = 99)
coef1 <- coef(fit, s = fit$lambda[10], exact = TRUE, x = X, y = Y)
coef2 <- coef(cvfit, s = fit$lambda[10], exact = TRUE, x = X, y = Y)
all.equal(coef1, coef2) # Mean relative difference: 0.002006215

修复

为确保估计匹配，我们可以用降低的阈值重新拟合两个对象。

fit <- glmnet(X, Y, intercept=FALSE, nlambda = 99, thresh = 1e-20)
cvfit <- cv.glmnet(X, Y, intercept= FALSE, thresh = 1e-20)
coef1 <- coef(fit, s = fit$lambda[10], exact = TRUE, x = X, y = Y)
coef2 <- coef(cvfit, s = fit$lambda[10], exact = TRUE, x = X, y = Y)
all.equal(coef1, coef2) #TRUE

请注意，您还可以在 coef() 中使用 thresh = 选项，但这仅在新惩罚不在任何一个的原始路径中时才有效对象。

fit <- glmnet(X, Y, intercept=FALSE, nlambda = 99)
cvfit <- cv.glmnet(X, Y, intercept= FALSE)
coef1 <- coef(fit, s = 0.40, exact = TRUE, x = X, y = Y, thresh = 1e-20)
coef2 <- coef(cvfit, s = 0.40, exact = TRUE, x = X, y = Y, thresh = 1e-20)
all.equal(coef1, coef2) #TRUE

thresh = 需要降低到的值将完全取决于数据。