Learn how easy it is to sync an existing GitHub or Google Code repo to a SourceForge project! See Demo

Close

dmstate, msm module

Help
2013-01-21
2013-01-21
  • mark pahuta
    mark pahuta
    2013-01-21

    Hi, I'm trying to implement the JAGS msm module. I've created a simulated 3-state, progressive, time-homogenous multisate model dataset. It is panel observed.

    I've fit it using the R msm package, the MLE estimates agree with my simulation parameters. But when I fit it with JAGS, my transition intensity estimates are off by a factor of 10. I suspect I've set-up the model incorrectly. I was hoping for some help?

    Data and functional R code can be downloaded from: https://docs.google.com/folder/d/0B2_rKFnvrjMAdi1FVVU0dHZVQzA/edit

    model{
      for(i in 1:Npat){
        for (j in 2:nobs[i]){
          S[i,j] ~ dmstate (S[i,j-1], dtm[i,j], Q[,])
        }
      }
    
      Q[2,2]<- -1 * Q[2,3]
      Q[1,1]<- -1 * (Q[1,2] + Q[1,3])
      Q[2,1]<-0
      Q[3,1]<-0
      Q[3,2]<-0
      Q[3,3]<-0
      Q[1,2] <- a12
      Q[1,3] <- a13
      Q[2,3] <- a23
    
      a12<- (-1 * log(0.5))/ med12
      a13<- (-1 * log(0.5))/ med13
      a23<- (-1 * log(0.5))/ med23
    
      med12 ~ dgamma(2.4, 0.01)T(1,) 
      med13 ~ dgamma(2.4, 0.01)T(1,)
      med23 ~ dgamma(2.4, 0.01)T(1,) 
    }
    
     
  • mark pahuta
    mark pahuta
    2013-01-21

    So I found some documentation on dmstate, and it seems the first parameter, x, in dmstate(x,y,z) is state occupation at time 0. so I've made the appropriate change to the code: S[i,1] from S[i,j-1].

    I now get a sensible value for a12 (which is truly 0.002295). But I still get estimates for a13 and a23 that are 10 times to small....

    model{
      for(i in 1:Npat){
        for (j in 2:nobs[i]){
          S[i,j] ~ dmstate (S[i,1], dtm[i,j], Q[,])
        }
      }
    
      Q[2,2]<- -1 * Q[2,3]
      Q[1,1]<- -1 * (Q[1,2] + Q[1,3])
      Q[2,1]<-0
      Q[3,1]<-0
      Q[3,2]<-0
      Q[3,3]<-0
      Q[1,2] <- a12
      Q[1,3] <- a13
      Q[2,3] <- a23
    
      a12<- (-1 * log(0.5))/ med12
      a13<- (-1 * log(0.5))/ med13
      a23<- (-1 * log(0.5))/ med23
    
      med12 ~ dgamma(2.4, 0.01)T(1,) 
      med13 ~ dgamma(2.4, 0.01)T(1,)
      med23 ~ dgamma(2.4, 0.01)T(1,) 
    }
    
     
    Last edit: mark pahuta 2013-01-21