Emilio A. Laca
2013-11-20
I am trying to model multistate data where the identities of subjects are not available. The data consist of the number of subjects in each state at each time. I would like to estimate the matrix of transition rates and the sojourn times for each state.
I think the distribution of the counts in the states at time t is a sum of multinomials. Is it possible to specify a likelihood equation in JAGS where the rhs is a sum of multinomials? If not, is there a possible work-around?
Any guidance would be greatly appreciated.
Here's what I have been trying without success:
Warning to others who might be trawling for code examples: this is my second attempt to model something in JAGS, so the code may be garbage.
I have 8 sets or cohorts. All start with 24 individuals in state a at time 0.
Time intervals are 2 or 3 days, so I have to use two different transition matrices.
There are only 6 possible transitions as indicated in the definition of the elements of Q.
I have the data organized in an array (19 dates, 5 states, 8 sets).
model{ #priors qab ~ dunif(lo, hi) qbc ~ dunif(lo, hi) qbh ~ dunif(lo, hi) qce ~ dunif(lo, hi) qch ~ dunif(lo, hi) qeh ~ dunif(lo, hi) # define transition matrix Q[1,1] <- -qab Q[1,2] <- 0 Q[1,3] <- 0 Q[1,4] <- 0 Q[1,5] <- 0 Q[2,1] <- 0 Q[2,2] <- -(qbc+qbh) Q[2,3] <- qbc Q[2,4] <- 0 Q[2,5] <- qbh Q[3,1] <- 0 Q[3,2] <- 0 Q[3,3] <- -(qce+qch) Q[3,4] <- qce Q[3,5] <- qch Q[4,1] <- 0 Q[4,2] <- 0 Q[4,3] <- 0 Q[4,4] <- -qeh Q[4,5] <- qeh Q[5,1] <- 0 Q[5,2] <- 0 Q[5,3] <- 0 Q[5,4] <- 0 Q[5,5] <- 0 # Transition vectors pi[i,j,t[k]] for (dts.t in 1:ndts){ pi.mat[1:5,1:5,dts.t] <- Imat5 %*% mexp(my.dt[dts.t]*Q) } #likelihood of the data for (set.j in 1:nsets){ for (date.k in 2:ndates){ from.a ~ dmulti(pi.mat[1,1:5,d.days.index[date.k]],counts[date.k-1,1,set.j]) from.b ~ dmulti(pi.mat[2,1:5,d.days.index[date.k]],counts[date.k-1,2,set.j]) from.c ~ dmulti(pi.mat[3,1:5,d.days.index[date.k]],counts[date.k-1,3,set.j]) from.e ~ dmulti(pi.mat[4,1:5,d.days.index[date.k]],counts[date.k-1,4,set.j]) from.h[1] <- 0 from.h[2] <- 0 from.h[3] <- 0 from.h[4] <- 0 from.h[5] <- counts[date.k-1,5,set.j] counts[date.k,1:5,set.j] <- from.a + from.b + from.c + from.e + from.h } } # end likelihood loop } # end model loop
Of course, this code results in an error:
Error in jags.model("jagsHRmodelBasic.R", data = list(counts = counts, :
RUNTIME ERROR:
Compilation error on line 57.
counts[2,1:5,1] is a logical node and cannot be observed
When I tried the following for likelihood
counts[date.k,1:5,set.j] ~ dmulti(c(0.5,0.5,0,0,0),24) dmulti(pi[state.i,1:5,set.j,d.days.index[date.k]],counts[date.k-1,1,set.j]) + dmulti(pi[state.i,1:5,set.j,d.days.index[date.k]],counts[date.k-1,2,set.j]) + dmulti(pi[state.i,1:5,set.j,d.days.index[date.k]],counts[date.k-1,3,set.j]) + dmulti(pi[state.i,1:5,set.j,d.days.index[date.k]],counts[date.k-1,4,set.j]) + from.h
I received the following error:
Error in jags.model("jagsHRmodelBasic.R", data = list(counts = counts, :
Error parsing model file:
syntax error on line 47 near "+"
I can send more details if anyone sees that there is any hope.
Thanks,
Emilio
Emilio A. Laca
2013-11-20
Correction 1: (Most likely not relevant for the issue) the line
Q[1,2] <- 0
should be
Q[1,2] <- qab
eal
Martyn Plummer
2013-11-22
People often want to model the observed sum of two or more unobserved random variables, which is why I introduced the dsum distribution (see the JAGS User Manual) and associated sampler.
Unfortunately, your problem will still not work because you have too many constraints:
1) The sum from.a[i] + from.b[i] + from.c[i] + from.e[i] + from.h[i]
is fixed for i = 1...5
. The dsum sampler is designed to work with this constraint.
2) But sum(from.a[1:5])
... sum(from.h[1:5])
are also fixed. This makes too many constraints for the dsum sampler.
I'm not saying it is impossible to do, just that it is not implemented.
Emilio A. Laca
2013-11-22
Thanks for taking the time to look into this. I write this mostly to confirm that the model using dsum() compiles but fails during runtime with the error below. I will continue to work on this trying alternative formulations if anything comes to mind, and will report to the forum when I succeed (or give up).
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph Size: 2398
Initializing model
Deleting model
Error in jags.model("jagsHRmodelBasic.R", data = list(counts = sim.counts, :
RUNTIME ERROR:
Inconsistent arguments for logDensity in distribution dsum
In addition: Warning message:
In jags.model("jagsHRmodelBasic.R", data = list(counts = sim.counts, :
Unused variable "n.draws" in data