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

Question

我正在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"是否已经隐含地假设了对数链接？

- RobertF

---

我认为这可能会很困难。如果你深入研究glmnet代码，你会发现glmnet(..., family="gaussian")调用了elnet，而elnet又调用了Fortran的spelnet函数。(泊松回归调用fishnet，它又调用了fishnet或者spfishnet（用于稠密和稀疏模型矩阵）。所以我怀疑需要有人从头开始编写一个变种的弹性网络来处理对数链接... - Ben Bolker

Answer 1

有点令人困惑，但是在glmnet和glm中的family参数是不同的。在glm中，您可以指定一个像"gaussian"这样的character，或者您可以指定一个带有某些参数的函数，如gaussian(link="log")。在glmnet中，您只能使用一个character来指定family，比如"gaussian"，并且没有办法通过该参数自动设置link。 gaussian的默认链接是identity函数，也就是根本没有转换。但是，请记住，链接函数只是您的y变量的一个转换；您可以自己指定它：

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

请注意，poisson系列的默认链接是log，但目标函数会改变。请参见第一段的?glmnet下的详细信息。

---

您的评论让我重新考虑了我的答案；我有证据表明它是不正确的。

正如您所指出的，E[log(Y)]和log(E[Y])之间存在差异。我认为上面的代码做的是拟合E[log(Y)]，这不是您想要的。以下是一些生成数据并确认您在评论中提到的内容的代码：

# 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 

- nograpes

---

说实话，在五分钟前我并没有意识到log[E(y)]与E[log(y)]并不完全等价。我运行了一些代码来测试glm(y~x,family=gaussian(link="log"))和glm(log(y)~x,family=gaussian)是否等价，结果它们确实是等价的，所以也许我在这里缺少一些基本的东西。 - nograpes

嗯，我在自己的数据上进行了双重检查，发现glm(log(y)x,family=gaussian)和glm(yx,family=gaussian(link="log"))的回归系数不同，而lm(log(y)x)和glm(log(y)x,family=gaussian)则产生相同的系数。 - RobertF

没关系 - 我会向 glmnet 作者确认一下。我无法相信他们会忽视这个问题，因为在医疗成本数据中使用具有 log 连接的高斯家族的广义线性模型都是非常常见的。 - RobertF

我确信他们忽视了这一点。这不是微不足道的改变，但可能需要重新思考底层代码/算法 - 希望不会太难，但也不是他们忘记包含的东西。 - Ben Bolker

@RNs_Ghost 不是 - 再次阅读nograpes的回答，唯一的可能性是在glmnet中使用Poisson链接通常被视为具有对数链接的高斯家族的“等效”。但是，如果您在nograpes的玩具数据集上运行以下代码：glmnet.model <- glmnet(x,y,family="poisson", thresh=1e-8, lambda=0); c(glmnet.model$a0, glmnet.model$beta[,1])，结果仍然不完全匹配glm()模型。 - RobertF

Answer 2

我知道这是一个老问题，但在当前版本的glmnet（4.0-2）中，可以使用glm系列函数作为“family”的参数，而不是一个字符字符串，因此您可以使用以下方式：

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

请注意，当您使用字符串参数时，该软件包的速度更快。

Reference: https://glmnet.stanford.edu/articles/glmnetFamily.html

- Andre Barbosa