4-level hierarchical model won't converge? document.SUBSCRIPTION_OPTIONS = { "thing": "thread", "subscribed": false, "url": "subscribe", "icon": { "css": "fa fa-envelope-o" } };

Help
Anonymous
2011-04-05
2012-09-01
• Anonymous - 2011-04-05

Hi,

I'm curious if anyone has tips or advice on getting a large, 4-level
hierarchical model to converge. In my data, I have 9180 loci (data points)
from 1200 populations nested within 320 studies from 135 species in 17 taxa.
In other words, this doesn't seem like a simple problem. However, the model
has only three important coefficients: one for data points (alpha.L), one for
studies (alpha.y), and one for species (alpha.s). The full code is below.

After running 2 million iterations (a couple days of crunching), alpha.L still
isn't stationary, and two chains haven't converged. Any advice on how to get
this working better? Any help or tips would be much appreciated.

Background: I'm running this in JAGS under Mac OS X with a 2.53 GHz Intel Core
2 Duo. I've also tried the problem under OpenBugs on a Dell PC with Windows
XP, with similar results.

Model:

model {

hyperpriors on the taxon means and the coefficients

lower.T ~ dunif(0,100)

upper.T ~ dunif(0,100)

alpha.L ~ dunif(-1, 1)

alpha.s ~ dunif(-1, 1)

alpha.y ~ dunif(-1, 1)

sigma.S ~ dunif(0, 100)

tau.S <- pow(sigma.S, -2)

sigma.Y ~ dunif(0, 100)

tau.Y <- pow(sigma.Y, -2)

sigma.P ~ dunif(0, 100)

tau.P <- pow(sigma.P, -2)

sigma.L ~ dunif(0, 100)

tau.L <- pow(sigma.L,-2)

go through each taxon

for (t in 1:n.taxa){

T is mean taxon diversity

T ~ dunif(lower.T, upper.T)

}

go through each species

for (s in 1:n.species){

Code the collapse effect as a proportional impact on diversity

Shat <- (1 + alpha.sC)T[taxon.list]

S is mean species diversity

S ~ dnorm(Shat, tau.S)

}

go through each study

for (y in 1:n.studies){

Primer note effect is a proportional effect

Yhat <- (1 + alpha.y * N) * S[species.list]

Y is mean study diversity

Y ~ dnorm(Yhat, tau.Y)

}

go through each population

for (p in 1:n.pops){

P

~
dnorm(Y[study.list

], tau.P)

}

go through each data point

for (i in 1:n.data){

Actual data points are distributed around the expected value for the species

plus an adjustment for whether or not the data point is from a cross-species
locus

Lhat_ <- (1 + alpha.LX_)P[pop.list_]

L_ ~ dnorm(Lhat_, tau.L)

}

}

Regards,

Malin_

• Martyn Plummer - 2011-04-06

Try flattening your model a bit by centring your random effects on zero.

```# Code the collapse effect as a proportional impact on diversity
Shat[t] <- (1 + alpha.s*C[t])*T[taxon.list[t]]
# S is mean species diversity
S[t] ~ dnorm(Shat[t], tau.S)

# Code the collapse effect as a proportional impact on diversity
Shat[t] <- (1 + alpha.s*C[t])*T[taxon.list[t]]
Serror[t] ~ dnorm(0, tau.S)
# S is mean species diversity
S[t] <- Shat[it] + Serror[t]
```

and so on. NB I changed the index from s to t to stop it being interpreted as
BBCode strikethrough.

You can also try loading the glm model, which will give you block updating,
but the compiler might choke on the model.

Martyn

• Anonymous - 2011-04-06

Hi Martyn,

Thanks for the good suggestion. Strangely, centering the random effects seemed
to make it run about 3x slower on a test piece of data. Perhaps because the
graph size increases?

I did have success speeding up convergence dramatically when I changed 3
things:

1) Used empirical means and standard deviations for the initial values of the
groups means (T, S, Y, and P) and standard deviations (sigma.L, sigma.P,
sigma.Y, sigma.S). I still use random initial values for the coefficients I
care about (alpha.s, alpha.y, and alpha.L).

2) Narrowed the prior for T and recoded without upper.T and lower.T:

```# go through each taxon
for (t in 1:n.taxa){
# T is mean taxon diversity
T[t] ~ dunif(0, 10)
}
```

I'm pretty sure #1 is the bit that sped this up so much. Can you see anything
wrong with this approach?

Thanks,

Malin

• Martyn Plummer - 2011-04-07

The automatically generated initial values can't be right for every problem,
so yes you should choose your own starting values when you think it
appropriate.