# Jags Implementation of Multivariate Response Probit Model

I am trying to implement the latent variable interpretation of a probit model with vector response (described on wiki here), but am receiving an error.

In this model, we have a matrix $X$, $n \times p$, of $p$ factors for $n$ data points, a matrix $Y$, $n \times q$, of $q$ outcomes for $n$ data points, and are trying to estimate a coefficient matrix $B$, $p \times q$, of $p$ times $q$ coefficients, one for each combination of factor and outcome. In order to estimate $B$, we must also wrestle with the covariance matrix $T$, $p \times p$, and the latent parameter matrix $Y^*$, $n \times q$.

Here is the model: $$Y^*_{i,1 \times q} = X_{i,1 \times p} B_{p \times q} + E_{i,1\times q} \ \text{for all i = 1, \ldots, n}$$ $$y_{i,q} = \begin{cases} 0, & \text{ if  y^*_{i,q} < 0, }\\ 1, & \text{ if  y^*_{i,q} > 0 } \end{cases}$$ $$B_q \sim MVNormal(\vec{b_0},B_0) \\ E_i \sim MVNormal(\vec{0},T_{q \times q}) \iff Y^*_i \sim MVNormal(X_{i,1\times p} B_{p \times q}, T_{q \times q}) \\ T \sim Wishart(I_{q \times q},p)$$

Here is the JAGS code I have written:

model  {

#Declare likelihood for Y, relationship between Y and Y_s
for (i in 1:n)  {
for (q in 1:Q) {
Y[i,q] ~ dinterval(Y_s[i,q],0)
mu[i,q] <- X[i,] %*% Beta[,q]
}
Y_s[i,1:Q] ~ dmnorm(mu[i,],Tau)
}

#Prior on Betas
for (q in 1:Q) {
Beta[1:P,q] ~ dmnorm(b_0,B_0)
}

#Prior on covariance matrix
Tau ~ dwish(ID,P)
}


Here is the error I recieve:

Error in node Y[1,2]
Node inconsistent with parents


To my understanding, this a very general error which doesn't give a good idea of where things are going screwy.

Here is a reproducible example, calling JAGS from R. I describe this example in my question here, under the heading Toy Example. Here is the example in brief:

• In the FooBar Baseball League, there are $q=2$ halls of fame

• Players with higher skill are more likely to be inducted into a hall of fame.

• The decision makers at each hall of fame are friends, and will "unduly" influence the probability that they induct a player. If one hall of fame likes a certain player, it is more likely that the other will as well.

The example:

##Imports
library(rjags)

##Params
n <- 100
testProp = 0.2
P <- 2
Q <- 2

##Simulation
#Generate Baseball Players with random skill (bounded between 0 and 80)
skill <- rnorm(n,40,10)
skill[skill > 80] <- 80
skill[skill < 0 ] <- 0

#So that there is some collusion between halls aside from player skill
hallFactor <- rbeta(100,0.5,0.5)+0.25

#probability of being admitted in Hall of Fame
p <- hallFactor*skill/100

#Randomly admit players to Halls of Fame A and B
HallA <- sapply(p, function(x) rbinom(1,1,x))
HallB <- sapply(p, function(x) rbinom(1,1,x))

##JAGS model fitting
X = cbind(1,skill)
Y = cbind(HallA,HallB)

ID <- diag(P)

b_0 <- rep(0,P)
B_0 <- diag(P)*0.00001
jags <- jags.model('/path/to/your/model/file.bug',
data = list('n' = n,
'X' = X,
'Y' = Y,
'Q' = Q,
'P' = P,
'b_0' = b_0,
'B_0' = B_0),
n.chains = 4,

jags.samples(jags,
c('Beta','T'),
1000)


Except if I'm mistaken the problem comes about by the fact that for JAGS Y_s has no initial values. By adding to your jags.model() the line
'Y_s' = Y,

jags will run without an error (whether initializing Y_s with Y is another matter: you might draw the initial values for Y_s from a multivariate normal distribution).