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) 

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

Similar Posts:

Rate this post

Leave a Comment