4-level hierarchical model won't converge?

  • Anonymous - 2011-04-05


    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 with four hierarchical levels: taxa, species, study, population

    Data input:

    n.taxa: number of taxa

    n.species: number of species

    n.pops: number of populations

    n.studies: number of studies

    n.data: number of data points

    X: 0 if the locus isn't cross-species, 1 if it is

    C: 0 if the species isn't collapsed, 1 if it is

    N: 0 if the paper isn't a primer note, 1 if it is

    taxon.list: the taxon that each species is from

    species.list: the species that each study is from

    study.list: the study that each population is from

    pop.list: the population that each data point is from

    L: each data point

    Deterministic Nodes

    tau: precision for data points within taxa

    Shat: expected diversity to each species

    Yhat: expected diversity for each study

    Lhat: expected diversity for each data point

    Output nodes:

    lower.T: lower limit to taxon diversity means

    upper.T: upper limit to taxon diversity means

    alpha.L: the cross-species locus effect

    alpha.s: the collapse effect

    alpha.y: the primer note effect

    sigma.S: standard deviation for species within taxa

    sigma.Y: standard deviation for study within species

    sigma.P: standard deviation for populations within studies

    sigma.L: standard deviation for data points within populations

    T: mean diversity for each taxon

    S: mean diversity for each species

    Y: mean diversity for each study

    P: mean diversity for each population

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



    ], 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

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

    L_ ~ dnorm(Lhat_, tau.L)





  • Martyn Plummer

    Martyn Plummer - 2011-04-06

    Try flattening your model a bit by centring your random effects on zero.
    Instead of this:

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


  • 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

    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)

    3) Loaded the glm module

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



  • Martyn Plummer

    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


Log in to post a comment.

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:

No, thanks