I have a transition matrix, describing the probability of an entity moving from one state to another in a time-period. I use this transition matrix to generate a series of "flow matrices", $F$ for a population across multiple time-periods, where $F_{i,j}$ denotes the number of transitions from state i to state j in that particular time point. e.g. $F = \left[ {\begin{array}{cc} 75 & 25 \\ 50 & 50 \\ \end{array} } \right]$ would indicate, for a 2-state system, that in this generation, there were 75 instances of state 1 that remained at state 1, while 25 instances of a state 1 transitioning to state 2. Likewise there were 50 instances of a state 2 transitioning to state 1, and 50 instances of a state 2 remaining as state 2. I can also define a total abundance vector $\theta$, denoting the total number of each state at the start of this phase. In this case, $\theta = \left[ {\begin{array}{c} 100 \\ 100 \\ \end{array} } \right]$.

Now, using multiple generations of these flow matrices, I can use JAGS to infer the background transition matrix, $P$, since each row of $F$ can be modelled via a multinomial distribution:

for(g in 1:(Generations-1)){

for(i in 1:N){

F[i,1:N,g] ~ dmulti(P[,i], theta[i,g]) }}

where `Generations' is the number of timepoints (and thus the number of $F$ matrices).

The problem is that I'm now trying to do the same thing, but instead using the multivariate normal approximation of the multinomial distribution. With mean and covariance parameters as defined here.

Here is the JAGS model code I'm using to run this:

for(g in 1:(Generations-1)){
for(i in 1:N){
#F[i,1:N,g] ~ dmulti(P[,i], theta[i,g])

for(Omega_i in 1:N){
for(Omega_j in 1:N){
Omega[Omega_i,Omega_j,i,g] <- ifelse(Omega_i==Omega_j, (theta[i,g]*P[Omega_i,i]*(1 - P[Omega_i,i])), -(theta[i,g]*P[Omega_i,i]*P[Omega_j,i]))

F[i,1:N,g] ~ dmnorm.vcov(P[,i]*theta[i,g], Omega[1:N,1:N,i,g])

for(i in 1:N){
P[1:N,i] ~ ddirch(alpha)

You can see in line 3 where I comment out the multinomial model. I then define the covariance matrix for each row of each F. However upon running this I get the error Error in node F[1,1:3,1] Invalid parent values. (when I simulate with N = 3).

Looking at similar queries online, the issue often seems to be when the arguments of the distribution are ill-posed, such as a precision matrix that can be non invertible, however I don't believe there are any requirements upon my model, and everything seems well-defined. The data I provide is $F$, $\theta$, $N$, Generations, alpha. Note that $N$ is the numer of states simulated.

Here is how I call the model:

model_output <- run.jags(model = jags_model, burnin = 5000, adapt = 100, sample = 500, n.chains = 3, thin = 5, method = 'parallel',
                         monitor = c('P'),
                         data = list(F=F,
                                     theta = theta,
                                     Generations = Generations, 
                                     N = N,
                                     alpha = rep(1,N)

Many thanks to anyone who can shine a lot on this for me!

EDIT: Okay, I'm now realising the issue is that the calculated Omega covariance matrix is not positive-definite. As such, the question is instead how on earth to make this approximation suitable?


Your Answer

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

Browse other questions tagged or ask your own question.