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)
Best Answer
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:
- Solved – error in JAGS beta-regression model R – value out of range in ‘gammafn’
- Solved – A problem using Bayesian modeling to impute a variable measured with error
- Solved – Code Problem in Lognormal Bayesian Model in JAGS/WinBUGS
- Solved – Code Problem in Lognormal Bayesian Model in JAGS/WinBUGS
- Solved – Correlation and group difference with unequal sample sizes