Alan Polansky
2012-03-28
I have been doing some research on constrained models and have recently read
the paper:
Gunn and Dunson (2005) "A Transformation Approach for Incorporating Monotone
or Unimodel Constraints", Biostatistics, 6, 434-449
In this paper they advocate fitting an unconstrained hierarchical model, and
then applying the constraint on the posterior distribution. They cite the fact
that the usual Gibbs sampling routines (Gelfand et. al., JASA, 1992, 523-532)
are difficult to apply to constrained parameter problems when fitting a
hierarchical model.
My question is whether JAGS requires this or not, or is it able to implement
the constraints in the prior (where I would like them implemented).
Suppose I have the following data:
X1 <- c(327,125,7,6,107,277,54)
X2 <- c(637,40,197,36,54,53,97,63,216,118)
N1 <- 7
N2 <- 10
and I want to fit an isotonic regression for set of X values with the
hierarchical model:
X1_ ~ exponential(theta1(i)), i = 1,...,N1
X2_ ~ exponential(theta2(i)), i = 1,...,N2
theta1_ ~ exponential(delta1)
theta2_ ~ exponential(delta2)
delta1 ~ exponential(lambda)
delta2 ~ exponential(lambda)
lambda us a specified constant and we constrain
theta1 > theta1 > ... > theta1
theta2 > theta2 > ... > theta2
I specify the JAGS model:
model {
for(i in 1:N1) {
X1_ ~ dexp(theta1_)
theta10_ ~ dexp(d1)
}
for(i in 1:N2) {
X2_ ~ dexp(theta2_)
theta20_ ~ dexp(d2)
}
d1 ~ dexp(d0)
d2 ~ dexp(d0)
d0 <- 0.01
theta11 <- sort(theta10)
theta21 <- sort(theta20)
for(i in 1:N1) {
theta1_ <- theta11
}
for(i in 1:N2) {
theta2_ <- theta21
}
}
JAGS compiles the model and seems to run fine, and the results seem okay.
But is it really fitting the model that I think it is fitting?
Thanks for any help or suggestions.__
Martyn Plummer
2012-03-29
Yes it should be the same model. You can double check comparing the results
with WinBUGS, which does allow you to define order constraints explicitly.
Alan Polansky
2012-03-30
Thanks for your comments - I am grateful. I am wondering if you know if there
is something different about the JAGS algorithm that allows it to sidestep the
issues referred to in Gunn and Dunson (2005), or is my model simple enough
that these issues do not apply?
Thanks again for your response and any additional help you may be able to give
Tim Handley
2012-04-12
I see that this model will produce theta1's and theta2's that meet your
constraints, but I'm concerned that the mechanism could be confusing to the
sampler. If one of you could help me to understand some additional details,
I'd appreciate it.
My specific concern is that there's a sort of aliasing going on. Sometimes,
theta1 corresponds to theta10, but sometimes it corresponds to theta10 or some
other of the theta10. During the course of sampling, the source of values for
theta1 will shift around among all of the theta10. While all of the theta10
have the same prior, I still feel like the shifting\aliasing could cause some
sort of problem. How can you sample from something when the source is
changing? Unfortunately, I don't understand the MCMC algorithms well enough to
know if there is a real problem or not.
Would anyone care to help illuminate this issue?
Alan, did you ever try this with WinBUGS for comparison?
Tim Handley
2012-04-16
After some thinking, I propose an alternative model construction that enforces
order on parameters in a more intuitive manner. Say we have a random effect
"alpha", for which we know alpha>alpha, essentially the same as Alan's
example, but with the reverse ordering. One way to code that effect is like
this:
alpha.delta[1] <- 0 alpha[1] <- dexp(.1) for(j in 2:N) { alpha.delta[j] ~ dexp(.1) alpha[j] <- alpha[j-1]+alpha.delta[j] }
Any comments and criticism would be appreciated.
Martyn Plummer
2012-04-17
If you want each of alpha ... alpha to have a marginal dexp(lambda)
distribution, but subject to the order constraint alpha < alpha < ... < alpha,
then you need different rate parameter for each component alpha.delta.
alpha[1] ~ dexp(lambda*N) for (j in 1:(N-1)) { alpha.delta[j] ~ dexp(lambda*(N-j)) alpha[j + 1] <- alpha[j] + alpha.delta[j] }
For example, alpha is the smallest of N exponential distributions with rate
parameter lambda, and is therefore exponential with rate N*lambda. This trick
works for exponentials, but not any other distribution.
Alan Polansky
2012-04-17
I cannot attempt a Winbugs comparison because I am running on OSX and do not
savor the idea or running windows on it.
I would really have the same questions there though, I am quite new to
Bayesian calculation and am trying to catch up, so I would probably have the
same concern with Winbugs due to my own ignorance.
Now, I did try a simpler model for the pooled data, and got similar looking
results. I also hard coded a similar model in R using a conjugate model and
also got similar result.