Work at SourceForge, help us to make it a better place! We have an immediate need for a Support Technician in our San Francisco or Denver office.

Close

sum of multinomial distributions

Help
2013-11-20
2013-11-22
  • Emilio A. Laca
    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
    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
    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
    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