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 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.s~~*C )*

# 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.L*X_)*P[pop.list_]

L_ ~ dnorm(Lhat_, tau.L)

}

}

Regards,

~~Malin~~**_**