Guilherme D. Garcia
Guilherme D. Garcia

Reputation: 179

glmer VS JAGS: different results in intercept-only hierarchical model

JAGS

I have an intercept-only logistic model in JAGS, defined as follows:

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)
}

When I plot the posterior distribution of b0 collapsing across all subjects (i.e., all b0[j]), my 95% HDI includes 0: -0.55 to 2.13. The Effective Sample Size is way above 10,000 for every b0 (around 18,000 on average). Diagnostics look good.

glmer()

Now, this is the equivalent glmer() model:

glmer(response ~ 1 + (1|subject), data = myData, family = "binomial")

The result of this model, however, is as follows:

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

So here it says my estimate is significantly above 0.

What the data look like

Here are the proportions of 0s and 1s by subject. You can see that, for the vast majority of subjects, the proportion of 1 is above 50%.

Any ideas why JAGS and glmer() are so different here?

      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

Upvotes: 4

Views: 282

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226182

You forgot to include a mean value, so your intercept parameter is fixed to zero. Something like this should work:

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)
}

Now the posterior density of mu0 should match the sampling distribution of the intercept parameter from glmer reasonably well.

Alternatively, if you use response ~ -1 + (1|subject) as your glmer formula, you should get results that match your current JAGS model.

Upvotes: 5

Related Questions