3
$\begingroup$

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,
                   n.adapt = 100)

jags.samples(jags,
             c('Beta','T'),
             1000)
$\endgroup$

1 Answer 1

2
$\begingroup$

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

$\endgroup$
1
  • $\begingroup$ I had given up on this particular problem (I just wrote the gibbs sampler myself), but as I've found since then, this does indeed appear to be the solution to this kind of problem. Specifically, the problem is that JAGS has to conjure up its own initial values if you don't provide them, and sometimes they don't make sense in the context of the data. So the moral is: initialize your chains. Thanks for your reply! $\endgroup$ Oct 4, 2016 at 1:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.