Constrained Hierarchical Models in JAGS

  • Alan Polansky

    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

    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

    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

    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

    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:[1] <- 0
        alpha[1] <- dexp(.1)
        for(j in 2:N)
[j] ~ dexp(.1)
          alpha[j] <- alpha[j-1][j]

    Any comments and criticism would be appreciated.

  • Martyn Plummer

    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[1] ~ dexp(lambda*N)
    for (j in 1:(N-1)) {[j] ~ dexp(lambda*(N-j))
       alpha[j + 1] <- alpha[j] +[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

    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.


Log in to post a comment.

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

Sign up for the SourceForge newsletter:

JavaScript is required for this form.

No, thanks