glmer VS JAGS:仅拦截分层模型的不同结果
glmer VS JAGS: different results in intercept-only hierarchical model
JAGS
我在 JAGS 中有一个仅拦截逻辑模型,定义如下:
model{
for(i in 1:Ny){
y[i] ~ dbern(mu[s[i]])
}
for(j in 1:Ns){
mu[j] <- ilogit(b0[j])
b0[j] ~ dnorm(0, sigma)
}
sigma ~ dunif(0, 100)
}
当我绘制 b0
的后验分布在所有对象(即所有 b0[j]
)中崩溃时,我的 95% HDI 包括 0
:-0.55 to 2.13
。每个 b0
的有效样本量远高于 10,000(平均约 18,000)。诊断看起来不错。
glmer()
现在,这是等效的 glmer()
模型:
glmer(response ~ 1 + (1|subject), data = myData, family = "binomial")
然而,这个模型的结果如下:
Random effects:
Groups Name Variance Std.Dev.
speaker (Intercept) 0.3317 0.576
Number of obs: 1544, groups: subject, 27
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.7401 0.1247 5.935 2.94e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
所以这里说我的估计明显高于 0
.
数据是什么样的
这里是 0s
和 1s
的比例(按主题)。可以看到,对于绝大多数科目来说,1
的比例都在50%以上。
知道为什么 JAGS 和 glmer()
在这里如此不同吗?
0 1
1 0.47 0.53
2 0.36 0.64
3 0.29 0.71
4 0.42 0.58
5 0.12 0.88
6 0.22 0.78
7 0.54 0.46
8 0.39 0.61
9 0.30 0.70
10 0.32 0.68
11 0.36 0.64
12 0.66 0.34
13 0.38 0.62
14 0.49 0.51
15 0.35 0.65
16 0.32 0.68
17 0.12 0.88
18 0.45 0.55
19 0.36 0.64
20 0.36 0.64
21 0.28 0.72
22 0.40 0.60
23 0.41 0.59
24 0.19 0.81
25 0.27 0.73
26 0.08 0.92
27 0.12 0.88
您忘记包含平均值,因此您的截距参数固定为零。这样的事情应该有效:
model{
for(i in 1:Ny){
y[i] ~ dbern(mu[s[i]])
}
for(j in 1:Ns){
mu[j] <- ilogit(b0[j])
b0[j] ~ dnorm(mu0, sigma)
}
mu0 ~ dnorm(0,0.001)
sigma ~ dunif(0, 100)
}
现在 mu0
的后验密度应该与 glmer
的截距参数的采样分布相当匹配。
或者,如果您使用 response ~ -1 + (1|subject)
作为您的 glmer
公式,您应该会得到与您当前的 JAGS 模型匹配的结果。
JAGS
我在 JAGS 中有一个仅拦截逻辑模型,定义如下:
model{
for(i in 1:Ny){
y[i] ~ dbern(mu[s[i]])
}
for(j in 1:Ns){
mu[j] <- ilogit(b0[j])
b0[j] ~ dnorm(0, sigma)
}
sigma ~ dunif(0, 100)
}
当我绘制 b0
的后验分布在所有对象(即所有 b0[j]
)中崩溃时,我的 95% HDI 包括 0
:-0.55 to 2.13
。每个 b0
的有效样本量远高于 10,000(平均约 18,000)。诊断看起来不错。
glmer()
现在,这是等效的 glmer()
模型:
glmer(response ~ 1 + (1|subject), data = myData, family = "binomial")
然而,这个模型的结果如下:
Random effects:
Groups Name Variance Std.Dev.
speaker (Intercept) 0.3317 0.576
Number of obs: 1544, groups: subject, 27
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.7401 0.1247 5.935 2.94e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
所以这里说我的估计明显高于 0
.
数据是什么样的
这里是 0s
和 1s
的比例(按主题)。可以看到,对于绝大多数科目来说,1
的比例都在50%以上。
知道为什么 JAGS 和 glmer()
在这里如此不同吗?
0 1
1 0.47 0.53
2 0.36 0.64
3 0.29 0.71
4 0.42 0.58
5 0.12 0.88
6 0.22 0.78
7 0.54 0.46
8 0.39 0.61
9 0.30 0.70
10 0.32 0.68
11 0.36 0.64
12 0.66 0.34
13 0.38 0.62
14 0.49 0.51
15 0.35 0.65
16 0.32 0.68
17 0.12 0.88
18 0.45 0.55
19 0.36 0.64
20 0.36 0.64
21 0.28 0.72
22 0.40 0.60
23 0.41 0.59
24 0.19 0.81
25 0.27 0.73
26 0.08 0.92
27 0.12 0.88
您忘记包含平均值,因此您的截距参数固定为零。这样的事情应该有效:
model{
for(i in 1:Ny){
y[i] ~ dbern(mu[s[i]])
}
for(j in 1:Ns){
mu[j] <- ilogit(b0[j])
b0[j] ~ dnorm(mu0, sigma)
}
mu0 ~ dnorm(0,0.001)
sigma ~ dunif(0, 100)
}
现在 mu0
的后验密度应该与 glmer
的截距参数的采样分布相当匹配。
或者,如果您使用 response ~ -1 + (1|subject)
作为您的 glmer
公式,您应该会得到与您当前的 JAGS 模型匹配的结果。