# =============================================================
# R code in C.1 Illustrative analysis of diabetic retinopathy
#
# December 11, 2018
# =============================================================


dcct <- dcct.ms[dcct.ms$cohort == 0,
                c("id","enum","visit","etime","state","statep","trx")]

# Data:
#   id enum visit etime state statep trx
#    1    1     0 0.000     1      1   0
#    1    2     2 0.671     1      1   0
#    1    3     4 1.180     1      1   0
#    1    4     6 1.665     1      1   0
#    1    5     8 2.166     1      1   0
#    1    6    10 2.680     1      1   0
#    1    7    12 3.181     2      2   0
#    1    8    14 3.680     2      2   0
#    1    9    16 4.227     3      3   0
#    1   10    18 4.646     2      3   0
#    1   11    20 5.095     3      3   0
#    1   12    22 5.626     2      3   0
#    1   13    24 6.182     2      3   0
#    1   14    26 6.642     2      3   0
#    1   15    28 7.121     3      3   0
#    1   16    30 7.619     3      3   0
#    1   17    32 8.159     3      3   0
#    1   18    34 8.676     2      3   0
#    1   19    41 9.161     3      3   0


library(msm)

# -------------------------------------------------------------
# C.1.1 Fitting the reversible markov model M1B with msm
# -------------------------------------------------------------

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

## Conventional
m1.conv <- msm(state ~ etime, subject=id, data=dcct[dcct$trx == 0,],
  qmatrix=mat.q, pci=c(3,6), center=FALSE, opt.method="optim")

qmatrix.msm(m1.conv, covariates=list('timeperiod[3,6)'=0, 'timeperiod[6,Inf)'=0), ci="delta")
qmatrix.msm(m1.conv, covariates=list('timeperiod[3,6)'=1, 'timeperiod[6,Inf)'=0), ci="delta")
qmatrix.msm(m1.conv, covariates=list('timeperiod[3,6)'=0, 'timeperiod[6,Inf)'=1), ci="delta")

## Intensive Therapy
m1.exp <- msm(state ~ etime, subject=id, data=dcct[dcct$trx == 1,],
  qmatrix=mat.q, pci=c(3,6), center=FALSE, opt.method="optim")

qmatrix.msm(m1.exp, covariates=list('timeperiod[3,6)'=0, 'timeperiod[6,Inf)'=0), ci="delta")
qmatrix.msm(m1.exp, covariates=list('timeperiod[3,6)'=1, 'timeperiod[6,Inf)'=0), ci="delta")
qmatrix.msm(m1.exp, covariates=list('timeperiod[3,6)'=0, 'timeperiod[6,Inf)'=1), ci="delta")


# -------------------------------------------------------------
# C.1.2 Fitting the progressive Markov model M2B with msm
# -------------------------------------------------------------

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

## Conventional
m2.conv <- msm(statep ~ etime, subject=id, data=dcct[dcct$trx == 0,],
  qmatrix=mat.q, pci=c(3,6), center=FALSE, opt.method="optim")

qmatrix.msm(m2.conv, covariates=list('timeperiod[3,6)'=0, 'timeperiod[6,Inf)'=0), ci="delta")
qmatrix.msm(m2.conv, covariates=list('timeperiod[3,6)'=1, 'timeperiod[6,Inf)'=0), ci="delta")
qmatrix.msm(m2.conv, covariates=list('timeperiod[3,6)'=0, 'timeperiod[6,Inf)'=1), ci="delta")

## Intensive Therapy
m2.exp <- msm(statep ~ etime, subject=id, data=dcct[dcct$trx == 1,],
  qmatrix=mat.q, pci=c(3,6), center=FALSE, opt.method="optim")

qmatrix.msm(m2.exp, covariates=list('timeperiod[3,6)'=0, 'timeperiod[6,Inf)'=0), ci="delta")
qmatrix.msm(m2.exp, covariates=list('timeperiod[3,6)'=1, 'timeperiod[6,Inf)'=0), ci="delta")
qmatrix.msm(m2.exp, covariates=list('timeperiod[3,6)'=0, 'timeperiod[6,Inf)'=1), ci="delta")


# -------------------------------------------------------------
# C.1.3  Fitting the hidden Markov model with msm
# -------------------------------------------------------------

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

emat  <- rbind(c(  0, 0.1, 0.1), 
               c(0.1,   0, 0.1), 
               c(0.1, 0.1,   0))

## Conventional
hmm.conv <- msm(state ~ etime, subject=id, data=dcct[dcct$trx == 0,],
  qmatrix=mat.q, ematrix=emat, pci=c(3,6), center=FALSE, opt.method="optim")

qmatrix.msm(hmm.conv, covariates=list('timeperiod[3,6)'=0, 'timeperiod[6,Inf)'=0), ci="delta")
qmatrix.msm(hmm.conv, covariates=list('timeperiod[3,6)'=1, 'timeperiod[6,Inf)'=0), ci="delta")
qmatrix.msm(hmm.conv, covariates=list('timeperiod[3,6)'=0, 'timeperiod[6,Inf)'=1), ci="delta")

ematrix.msm(hmm.conv, covariates=0, ci="delta")

## Intensive Therapy
hmm.exp <- msm(state ~ etime, subject=id, data=dcct[dcct$trx == 1,],
  qmatrix=mat.q, ematrix=emat, pci=c(3,6), center=FALSE, opt.method="optim")

qmatrix.msm(hmm.exp, covariates=list('timeperiod[3,6)'=0, 'timeperiod[6,Inf)'=0), ci="delta")
qmatrix.msm(hmm.exp, covariates=list('timeperiod[3,6)'=1, 'timeperiod[6,Inf)'=0), ci="delta")
qmatrix.msm(hmm.exp, covariates=list('timeperiod[3,6)'=0, 'timeperiod[6,Inf)'=1), ci="delta")

ematrix.msm(hmm.exp, covariates=0, ci="delta")