Dear BUGS USERS,
KEYWORDS
ragged array
hierarchical centering
nested indexing
I have a problem with trying to implement a logistic regression with
random effects.
I have managed to get my model to run but I want to achieve a better
parameterisation as sampling is poor and convergence is very slow.
I wish to model a binomial response given the effect of 3 factors (with
all three 2-way and one 3-way interactions) and 2 random effects. The
data are population responses to one of three treatments, they can be
one of three growth forms and of three ordinal sizes, The random effects
are the species from which the population comes and the paper from which
the data come. Presently the random effects are simply additive terms
in a regression. Thus:
model {
for( i in 1 : 892) { sprt[i] ~ dbin(p[i],n[i])
logit(p[i]) <- alpha0
+ alpha1[lf[i]]
+ alpha2[dist[i]]
+ alpha3[sz4[i]]
+ alpha12[lf[i],dist[i]]
+ alpha13[lf[i],sz4[i]]
+ alpha23[dist[i],sz4[i]]
+ alpha123[lf[i],dist[i],sz4[i]]
+ pap[paper[i]]
+ sps[spp[i]] }
This samples very slowly and produces highly correlated chains.
I understand that centering covariates can be useful. However, I have
only seen this usedfor continuous covariates or binary covariates. Also,
I use the factors for indexing and to calculate cell means (like in the
ALLIGATORS example).
I have been trying to implement a hierarchical centering similar to that
used in the EPILEPSY example in BUGS 0.5 manual. One major problem seems
to be that I have very ragged arrays. So to estimate the contribution
for the species random effect I have to loop around j:
#loop around species, estimate random effect
for(j in 1:561){
sps[j] ~ dnorm(0.0,tau)}
tau~dgamma(0.001,0.001) }
I find it hard to work out how to associate a random effect for a
particular species [j] with the particular record [i] other than using
an additive term.
Is it possible to use hierarchical centering when the data are not
completely hierarchical?
Going right back to basics, removing all my factors I can get an
additive model like this:
model;
{
for( i in 1 : 892 ) {
sprt[i] ~ dbin(p[i],n[i])
logit(p[i]) <- mu[i]
mu[i] <- alpha0 + sps[spp[i]]
}
for( j in 1 : 561 ) {
sps[j] ~ dnorm( 0.0,tau)
}
tau ~ dgamma(0.001,0.001)
alpha0 ~ dnorm( 0.0,1.0E-6)
}
This runs, though chains for alpha0 and individual sps[j] are fairly
autocorrelated. What I would like to do is to take the sps additive term
out to introduce mu[i] ~ dnorm(alpha0, tau). But if I use code like that
below, the variation in mu[i] is completely random, as opposed to
associated with a species random effect. How can I associate the random
effect appropriately?
model;
{
for( i in 1 : 892 ) {
sprt[i] ~ dbin(p[i],n[i])
logit(p[i]) <- mu[i]
mu[i] ~ dnorm(theta[i],tau)
theta[i] <- alpha0
} for( j in 1 : 561 ) {
sps[j] ~ dnorm( 0.0,tau)
}
tau ~ dgamma(0.001,0.001)
alpha0 ~ dnorm( 0.0,1.0E-6)
}
Any help would be welcome,
Peter
--
Peter Vesk
Lecturer
School of Botany
The University of Melbourne, Victoria 3010
Australia
[log in to unmask]
+61 (0)3 8344 7480
-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask]
To mail the BUGS list, mail to [log in to unmask]
Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html
Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list
|