# ===============================================================
# Generating multistate data with intermittent observation
# See PDF file for details
#
# December 20, 2018
# ===============================================================


library(msm)


generatedata.f <- function(m, CA, init.pi, qmat, rho) {
  diag(qmat) <- (-1)*apply(qmat, 1, sum)
  nstates <- ncol(qmat)

  outdata <- lapply(1:m, function(id, nstates, init.pi, qmat, rho, CA) {
    ai <- 0
    whichstate <- as.vector(rmultinom(1, size=1, prob=init.pi))
    Zai <- c(1:nstates)[whichstate == 1]

# Alternative approach:
#    Zai <- sample(c(1:nstates), size=1, replace=FALSE, prob=init.pi)

    r <- 1;  air <- ai[1]
    while (air < CA) {
      r <- r + 1
      Delta.air <- rexp(1, rate=rho)
      air <- Delta.air + ai[r-1]

      if (air < CA) {
        ai <- c(ai, air)
        PP <- MatrixExp(qmat, t=(ai[r] - ai[r-1]), n=50, k=3, method="series")

        whichstate <- as.vector(rmultinom(1, size=1, prob=PP[Zai[r-1],]))
        Zai <- c(Zai, c(1:nstates)[whichstate == 1])

# Alternative approach:
#       Zai <- c(Zai, sample(c(1:nstates), size=1, replace=FALSE, prob=PP[Zai[r-1],]))
      }
      else { break }
    }    

    ai <- c(ai, CA)
    Zai <- c(Zai, 999)    

    return( data.frame(id=rep(id, length(ai)), times=ai, states=Zai) )
  }, nstates=nstates, init.pi=init.pi, qmat=qmat, rho=rho, CA=CA)
  outdata <- do.call("rbind", outdata)
  outdata <- outdata[order(outdata$id, outdata$times),]
  return(outdata)   
}


set.seed(1000)
simdata <- generatedata.f(m=1000,
                          CA=1,
                          init.pi=c(0.4,0.3,0.3),
                          qmat=rbind(c(0,0.6,0.4), c(0.3,0,0.2), c(0.2,0.1,0)),
                          rho=5)
    
print(statetable.msm(state=states, subject=id, data=simdata))

fit <- msm(states ~ times, subject=id, data=simdata,
           qmatrix=rbind(c(0,0.6,0.4), c(0.3,0,0.2), c(0.2,0.1,0)),
           censor=999, censor.states=c(1,2,3),
           center=FALSE)
print(fit)