我试图通过尝试为它们生成数据并检查输出行为的方式来理解某些概念是如何工作的。目前,我因此意识到我不太了解 GLM-s 的情况。
这是我的小代码:
N = 10000
e = rnorm(N,0,1)
x1 = runif(N,10,30)
y = exp(5*x1+ 10 + e)
mod1 = glm (y ~ x1,family=gaussian(link="log"))
mod2 = lm(log(y) ~ x1)
调用模型的摘要很快就会发现,这mod2是一个很好的选择,但mod1很疯狂。我试着复习这个话题,很多页面都谈到了转换 的平均值y,因为y它不是正态分布的,但我从来没有真正理解这一点,因为正态性假设是针对残差,而不是,否则永远不会与均匀采样一起使用。
所以我有两个问题:
- 我没有得到什么?
- 我将如何生成对上述 GLM 有效的数据?
编辑
我重新编写了代码,以更密切地反映数学背景(基于 Glen_b 的回答,我意识到我添加错误的方式并不适用于所有情况)。
x = seq(from = 1,to = 15,by = 0.1)
N = length(x)
eta = 5*x + 10
# original
set.seed(5671)
y = exp(eta) + rnorm(N,0,1)
mod = glm(y ~ x,gaussian(link = "log"))
# new
set.seed(5671)
inverse_link = function(x){exp(x)}
means = sapply(eta,function(x){inverse_link(x)})
y = sapply(means,function(x){rnorm(1,mean=x,sd=1)})
mod = glm(y ~ x,gaussian(link = "log"))
可以比较两种情况下的结果是相同的。基于此,我的期望是以下代码可以正确拟合我的参数:
x = seq(from = 1,to = 15,by = 0.1)
N = length(x)
eta = 5*x + 10
set.seed(5671)
inverse_link = function(x){1/x}
scale = 1
shapes = sapply(eta,function(x){inverse_link(x)/scale})
y = sapply(shapes,function(x){rgamma(1,shape=x,scale=scale)})
mod = glm (y ~ x,family=Gamma(link="inverse"))
我的理由是所以我需要一个形状参数的伽马分布。我的问题是,系数非常偏离(约为 1.6 ,截距约为 6.4)。只是我的输入数据,还是我错过了什么?
编辑 2
中指出的那样,形状参数保持不变(据我所知,GLM 假设使用的分布来自自然指数族,基于这个答案,伽马是具有固定形状参数的分布)。所以我们有
这是现在有效的更正代码:
x = seq(from = 1,to = 15,by = 0.1)
N = length(x)
eta = 5*x + 10
set.seed(5671)
inverse_link = function(x){1/x}
shape = 3
scales = sapply(eta,function(x){inverse_link(x)/shape})
y = sapply(scales,function(x){rgamma(1,shape=shape,scale=x)})
mod = glm (y ~ x,family=Gamma(link="inverse"))
summary(mod)
```