glmer VS 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% ИРЧП включает 0: -0.55 to 2.13. Эффективный размер выборки намного превышает 10 000 для каждого b0 (в среднем около 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
1 4
glm

1 ответ:

Вы забыли включить среднее значение, поэтому ваш параметр перехвата фиксируется на нуле. Что-то вроде этого должно сработать: 

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.
5