# ===========================================================
# C.2 Code for the onset of arthritis mutilans in PsA
# R code in C.2.2 Marginal model for time to entry to the absorbing state
#
# December 11, 2018
# ===========================================================


library(msm)

# -----------------------------------------------------------
# Fit a marginal model
# -----------------------------------------------------------

mutilans$state6 <- ifelse(mutilans$state == 6, 2, 1)
mutilans[mutilans$id == 9, c("id","sex.female","age.psa","b27","c3",
                             "times","state","state6","status")]

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


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

fitM12 <- msm(state6 ~ times, subject=id, data=mutilans[mutilans$status == 1,],
  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))


# -----------------------------------------------------------
# Create mutilansTTE dataframe
# -----------------------------------------------------------

mutilans$estatus <- mutilans$status

mutilansTTE <- lapply(sort(unique(mutilans$id)), function(pid, indata) {
  datai <- indata[indata$id == pid,]
  datai <- datai[order(datai$times),]

  if ( sum(datai$state6 == 2) > 0 ) {
    estop   <- min( datai$times[datai$state6 == 2] )
    estart  <- max( datai$times[datai$times < estop] ) 
    estatus <- 3
  }
  else {
    estop   <- NA
    estart  <- max(datai$times)
    estatus <- 0
  }
  return(data.frame(id=pid, estart, estop, estatus))
}, indata=mutilans[mutilans$estatus == 1,])
mutilansTTE <- do.call("rbind", mutilansTTE)
mutilansTTE[mutilansTTE$id %in% c(1,5,11,21),]

# id estart estop estatus
#  1 31.403    NA       0
#  5  0.000 26.73       3
# 11 21.205    NA       0
# 21 18.768    NA       0

write.table(mutilansTTE, file="mutilansTTE.dat", 
            append=TRUE, row.names=FALSE, col.names=TRUE)

## Note that Turnbull estimates will be obtained using kaplanMeier function
## in TIBCO Spotfire S+