# Solved – 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,1times 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,1times 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,                    n.adapt = 100)  jags.samples(jags,              c('Beta','T'),              1000) ``
Contents

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