sum of multinomial distributions

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


    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, :
    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]) +

    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.



  • 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


  • 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, :
    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


Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:

No, thanks