Range constraints for MVN nodes - a solution?

[Richard James]

Hello fellow JAGS users,

The model I am running (see below) has a multivariate normal node "S" which I require to be constrainted to a 0-100 interval (i.e. no negative values). I've tried using the "I(0,100) and "T(0,100)" terms to no avail, and I note that the manual it states that "I" cannot be used to constrain MVN nodes. Has anyone come across a way around this issue, as currently I'm generating negative S values that are illogical?

Thanks for any assistance given.

   model {
  # target data likelihood
  for(i in 1:N) {
    Y[i,1:J] ~ dmnorm(mu[i,],OmegaresZ)
    mu[i,1:J] <- p[i,1:K]%*%s[1:K,1:J,i]
  }

  # s
  for(i in 1:N) {
    for(k in 1:K) {
      s[k,1:J,i] ~ dmnorm(muS[,k],OmegaS[,,k])
    }
  }

  # p
  for(i in 1:N) {
    p[i,1:K] <- expphiV[i,1:K]/sum(expphiV[i,1:K])
    phiV[i,1:K] <- phi[i,1:(K-1)]%*%tV
    for(k in 1:K) {
      expphiV[i,k] <- exp(phiV[i,k])
    }
  }

  # phi
  for(i in 1:N) {
    for(k in 1:(K-1)) {
      phi[i,k] ~ dnorm(muphi[k],tauphi[k])
    }
    philast[i] <- -sum(phi[i,1:(K-1)])
  }

  # hyperparameters for phi prior
  for(k in 1:(K-1)) {
  muphi[k] ~ dnorm(0,1000)    
         tauphi[k] ~ dgamma(2,1)
    sigma2phi[k] <- 1/tauphi[k]
  }

  # Prior on OmegaresZ
  # combined precision matrix of residual + measurement error
  OmegaresZ ~ dwish(Rres,kres)
  SigmaresZ <- inverse(OmegaresZ) # for output

}

[Martyn Plummer]

The multivariate normal distribution can't be truncated. However, in your example it appears that the mean and variance of s are observed, because they do not appear on the left-hand side of any relation. If this is the case you can use the dinterval distribution to apply the constraint. Normally dinterval represents censoring, not truncation (see the manual) but when the parameters of a distribution are fixed they are the same.

To do this you would need to add a data step, or supply the equivalent data:

data {
   # s
   for(i in 1:N) {
      for(k in 1:K) {
         for (j in 1:J) {
            positive[i,j,k] <- 1
      }
   }
   zero <- 0
}

Then in your model

   # s
   for(i in 1:N) {
      for(k in 1:K) {
         for (j in 1:J) {
            positive[i,j,k] ~ dinterval(s[i,j,k], zero)
      }
   }

This encodes the information that every element of the array s is positive.

JAGS will fall back on a rather simple sampling method for s - the random walk Metropolis sampler so monitor the mixing of these parameters.

[mauricio]

Two question regarding this solution.
1) Would this work if I want to encode the information that every element of the array is negative?
2) Would this work if I want to encode the information that some elements of the array are positive and some are negative?

[Martyn Plummer]

Yes, the dinterval distribution is flexible enough to include these different constraints. For example, if you have a 2-dimensional multivariate normal that is unobserved:

Z[1:2] ~ dmnorm(mu[1:2], Tau[1:2,1:2])

and both mu and Tau are fixed quantities. Then you can censor Z[1] as positive and Z[2] as negative with:

~~~
for (i in 1:2) {
S[i] ~ dinterval(Z[i], 0)
}
~~~

where S is supplied as data:

~~~
S <- c(1, 0)
~~~

You will also need to supply starting values for Z such that Z[1] starts with a positive value and Z[2] starts with a negative value.