了解 glm 和 link 函数:如何生成数据?

机器算法验证 r 广义线性模型 随机生成 链接功能
2022-03-31 07:27:42

我试图通过尝试为它们生成数据并检查输出行为的方式来理解某些概念是如何工作的。目前,我因此意识到我不太了解 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它不是正态分布的,但我从来没有真正理解这一点,因为正态性假设是针对残差,而不是,否则永远不会与均匀采样一起使用。yy=mx+bx

所以我有两个问题:

  • 我没有得到什么?
  • 我将如何生成对上述 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)。只是我的输入数据,还是我错过了什么?μ=kθ=1/ηk=1/η/θx

编辑 2

中指出的那样,形状参数保持不变(据我所知,GLM 假设使用的分布来自自然指数族,基于这个答案,伽马是具有固定形状参数的分布)。所以我们有k

yΓ(k,1kη)

这是现在有效的更正代码:

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)
```
2个回答

exp()调用中是否包含错误术语很重要。这是此代码的大问题。考虑一下:

#  first, set the random seed, so that everything is reproducible
set.seed(5671)
N  = 10000
e  = rnorm(N,0,1)
x1 = runif(N,10,30)
y1 = exp(5*x1+ 10  + e)
y2 = exp(5*x1+ 10) + e
mod1.1 = glm(y1 ~ x1,family=gaussian(link="log"))
mod1.2 = lm(log(y1) ~ x1)

mod2.1 = glm(y2 ~ x1,family=gaussian(link="log"))
mod2.2 = lm(log(y2) ~ x1)

以下是从 glm 生成的方法(可以移动某些项目的顺序):

  1. 选择您的家庭和链接功能。

  2. 为要模拟的每个观察选择预测变量 (IV)。

  3. 选择你的系数。

  4. 评估每个观测值的线性预测值。

  5. 通过链接函数的倒数进行变换以获得每个观察的条件均值。

  6. 选择任何其他参数。

  7. 在每次观察时对分布进行采样,您现在拥有所有参数。


让我们看看如何使用反向链接模拟一个简单的 Gamma GLM,遵循以下步骤:

  1. 选择您的家庭和链接功能。伽玛,逆

  2. 为要模拟的每个观察选择预测变量 (IV)。( )x

  3. 选择你的系数。在这种情况下选择一个特定的 &β0β1

  4. 评估每个观测值的线性预测值。( )ηi=β0+β1xi,i=1,...,n

  5. 通过链接函数的倒数进行变换以获得每个观察的条件均值。的倒数是ηi=1/μiμi=1/ηi

  6. 选择任何其他参数。(选择形状参数)

  7. 在每次观察时对分布进行采样,您现在拥有所有参数。(例如y=rgamma(length(x),shape,scale=mu/shape)——注意scale这里是一个值向量)