Solved – Problems of specifying constrained parameters in a hierarchical model under JAGS

jagsmultilevel-analysisr

I'm trying to fit a model under JAGS where I have a response, y, that depends linearly on eight predictors (with no intersect : I'm actually trying to evaluate the share of each socioeconomic statuts population that voted for a candidate, wth aggregated data).

My problem is that naive hierarchical modelling, under lmer() in R or in JAGS, gives me some negative estimates of my interest parameters, which is logically impossible. So I would like to constrain these parameters to take values between 0 and 1. Should I just cut off the posteriors ? I actually would like to integrate it in the model, but the problem is that since there are eight predictors, I use multivariate normal distribution, and I dont know how to use the dinterval() (or alternatively T() ?) distribution on it.

Here is the model : 

model{
   for (i in 1:n) {
    y[i] ~ dnorm(y.hat[i], tau.y)
    y.hat[i] <- inprod(B[dpt[i],], X[i,])
   }

   tau.y <- pow(sigma.y, -2)
   sigma.y ~ dunif(0, 100)

   for (j in 1:J) {
    for (k in 1:K) {
      B[j,k] <- xi[k] * B.raw[j,k]
    }
    B.raw[j,1:K] ~ dmnorm(mu.raw[], Tau.B.raw[,])
   }

   for (k in 1:K) {
    is.interval[k] ~ dinterval(mu[k], lim)
    mu[k] <- xi[k]*mu.raw[k]
    mu.raw[k] ~ dnorm(0, 0.0001)
    xi[k] ~ dunif(0,100)
   }

   Tau.B.raw[1:K,1:K] ~ dwish(W[,], df)
   df <- K+1
   Sigma.B.raw[1:K,1:K] <- inverse(Tau.B.raw[,])
   for (k in 1:K) {
    for (k.prime in 1:K) {
      rho.B[k, k.prime] <- Sigma.B.raw[k, k.prime]/
        sqrt(Sigma.B.raw[k,k]*Sigma.B.raw[k.prime, k.prime])
    }
    sigma.B[k] <- abs(xi[k])*sqrt(Sigma.B.raw[k,k])
    }
    }

with n, J, K(=8), y, dpt, X, W, lim, and is.interval provided as data. My main parameter of interest is B.

Thanks you very much for any help here !

Best Answer

There are several solutions to your problem. One is, as you suggest, just to collect all the MCMC output then delete the draws for which the parameters of interest are outside of [0,1]. You can think of this as resampling from your posterior using importance sampling, where the importance sampling weight is 0 if the parameters are outside [0,1] and 1 if they are OK; a little thought should convince you that this works.

Another approach, as you have almost realized, is to use the I() construct (what you want in JAGS) on the prior. Here's a line from a JAGS program that is running even as I type:

  b.dni ~ dnorm(-0.00021, 100000)I(-0.005,0.005)

constraining b.dni to be between -0.005 and 0.005. This naturally works better if there's any substantial probability in the unconstrained case of the parameter falling outside the range, since you don't waste MCMC samples, albeit at a slightly higher computational cost per sample.

Another approach would be to abandon the Normal distribution for one with compact support, e.g., a Beta distribution:

b.dni.raw ~ dbeta(a,b)
b.dni <- lb.b.dni + (ub.b.dni - lb.b.dni) * b.dni.raw

Distributions on the hyperparameters a and b can be constructed more intuitively than working with a and b themselves by putting distributions on the mean and std. deviation of b.dni, then assigning the appropriate values to a and b based on the sampled mean and standard deviation. This is somewhat more of a pain programmatically, though it does gain you some flexibility in prior shape.

Related Question