Learn how easy it is to sync an existing GitHub or Google Code repo to a SourceForge project! See Demo

Close

Constrained Hierarchical Models in JAGS

2012-03-28
2012-09-01
  • 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:

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