# ===========================================================
# C.2 Code for the onset of arthritis mutilans in PsA
# R code in C.2.1 Data frame and fit of intensity-based model
#
# December 11, 2018
# ===========================================================


library(msm)


mutilans[mutilans$id == 9, c("id","sex.female","age.psa","b27","c3",
                             "times","state","status")]

# id sex.female age.psa b27 c3  times state status
#  9          1      26   0  0  0.000     1      1
#  9          1      26   0  0  9.489     2      1
#  9          1      26   0  0 11.488     2      1
#  9          1      26   0  0 13.528     2      1
#  9          1      26   0  0 20.627     4      1
#  9          1      26   0  0 23.072     6      1
#  9          1      26   0  0 45.024   999      0


mat.q <- rbind(c(-0.1,  0.1,    0,    0,    0,   0), 
               c(   0, -0.1,  0.1,    0,    0,   0),
               c(   0,    0, -0.1,  0.1,    0,   0), 
               c(   0,    0,    0, -0.1,  0.1,   0),
               c(   0,    0,    0,    0, -0.1, 0.1), 
               c(   0,    0,    0,    0,    0,   0))  


fitM <- msm(state ~ times, subject=id, data=mutilans[mutilans$status == 1,],
  covariates = ~ sex.female + age.psa + b27 + c3,
  constraint=list(sex.female=rep(1,5), age.psa=rep(1,5), b27=rep(1,5), c3=rep(1,5)),
  qmatrix=mat.q, gen.inits=TRUE, pci=c(6,12,18), opt.method="optim", center=FALSE,
  control=list(trace=2, fnscale=100000, reltol=1e-10, maxit=10000))