如何在 glmnet 中指定日志链接？答案

【问题标题】：How to specify log link in glmnet?如何在 glmnet 中指定日志链接？
【发布时间】：2014-10-02 03:08:32
【问题描述】：

我正在使用 R 中的 glmnet 和 caret 包在广义线性模型上运行弹性网络。

我的响应变量是成本（其中成本 > 0 美元），因此我想为我的 GLM 指定一个带有日志链接的高斯族。但是 glmnet 似乎不允许我指定 (link="log") 如下：

> lasso_fit <- glmnet(x, y, alpha=1, family="gaussian"(link="log"), lambda.min.ratio=.001)

我尝试了不同的变体，有无引号，但没有运气。 glmnet 文档没有讨论如何包含日志链接。

我错过了什么吗？ family="gaussian" 是否已经隐式假设了一个日志链接？

【问题讨论】：

我认为这可能很困难。如果您深入研究glmnet 代码，您会看到glmnet(..., family="gaussian") 调用elnet，它调用Fortran 的spelnet 函数。（泊松回归调用fishnet，调用fishnet 或spfishnet（用于密集与稀疏模型矩阵。）所以我怀疑有人必须从头开始编写处理日志链接的弹性网络的变体.. .

标签： r glmnet

【解决方案1】：

这有点令人困惑，但是glmnet 和glm 中的family 参数是完全不同的。在glm 中，您可以指定character，例如"gaussian"，或者您可以指定带有一些参数的函数，例如gaussian(link="log")。在glmnet 中，只能用character 指定族，如"gaussian"，无法通过该参数自动设置链接。

gaussian 的默认链接是identity 函数，即根本没有转换。但是，请记住，链接函数只是您的 y 变量的转换；你可以自己指定：

glmnet(x, log(y), family="gaussian")

还要注意poisson 系列的默认链接是log，但目标函数会发生变化。请参阅前几段中?glmnet 下的详细信息。

你们的 cmets 让我重新思考我的答案；我有证据表明它不正确。

正如您所指出的，E[log(Y)] 和 log(E[Y]) 之间存在差异。我认为上面的代码所做的是适合 E[log(Y)]，这不是你想要的。下面是一些代码，用于生成数据并确认您在 cmets 中记录的内容：

# Generate data
set.seed(1)
x <- replicate(3,runif(1000))
y <- exp(2*x[,1] + 3*x[,2] + x[,3] + runif(1000))
df <- data.frame(y,x)

# Run the model you *want*
glm(y~., family=gaussian(link="log"), data=df)$coef
# (Intercept)          X1          X2          X3 
#   0.4977746   2.0449443   3.0812333   0.9451073 

# Run the model you *don't want* (in two ways)    
glm(log(y)~., family=gaussian(link='identity'), data=df)$coef
# (Intercept)          X1          X2          X3 
#   0.4726745   2.0395798   3.0167274   0.9957110 
lm(log(y)~.,data=df)$coef
# (Intercept)          X1          X2          X3 
#   0.4726745   2.0395798   3.0167274   0.9957110 

# Run the glmnet code that I suggested - getting what you *don't want*.
library(glmnet)
glmnet.model <- glmnet(x,log(y),family="gaussian", thresh=1e-8, lambda=0)
c(glmnet.model$a0, glmnet.model$beta[,1])
#        s0        V1        V2        V3 
# 0.4726745 2.0395798 3.0167274 0.9957110

【讨论】：

太好了，谢谢。如果我按照您的建议指定 glmnet(x, log(y), family="gaussian") ，则 glmnet 实际上将估计 log[E(y)] （我想要的）而不是 E[log(y)] （a我不想要的对数线性 OLS），对吗？
老实说，直到我在五分钟前查看它时，我才意识到 log[E(y)] 并不完全等同于 E[log(y)]。我运行了一些代码来测试 glm(y~x,family=gaussian(link="log")) 是否等同于 glm(log(y)~x,family=gaussian)，确实如此，所以也许我在这里遗漏了一些基本的东西。
嗯，我仔细检查了自己的数据，我看到 glm(log(y)~x,family=gaussian) 与 glm(y~x,family 的回归系数不同=gaussian(link="log"))，而 lm(log(y)~x) 和 glm(log(y)~x,family=gaussian) 产生相同的系数。
别担心 - 我会与 glmnet 作者核实。我不敢相信他们忽略了这一点，因为带有日志链接的高斯家庭的 GLM 在医疗保健成本数据中相当普遍。
@RNs_Ghost 否 - 再次阅读 nograpes 的响应，唯一的另一种可能性是在 glmnet 中使用泊松链接通常用作具有日志链接的高斯族的“等价物”。但是，如果您在 nograpes 玩具数据集上运行以下代码：glmnet.model <- glmnet(x,y,family="poisson", thresh=1e-8, lambda=0); c(glmnet.model$a0, glmnet.model$beta[,1])，结果仍然与glm() 模型不太匹配。

【解决方案2】：

我知道这是一个老问题，但在当前版本的 glmnet (4.0-2) 中，是可以使用 glm 系列函数作为“family”的参数而不是 a字符串，所以你可以使用：

glmnet(x, y, family=gaussian(link="log"))

请注意，使用字符串参数时，包会更快。

参考： https://glmnet.stanford.edu/articles/glmnetFamily.html

【讨论】：