# Chapter 8. MC models with CAV data. Fitted with msm() 

# Ardo, UCL 2017.

# Prelim:
library(msm)
digits <- 3

# Data:
cat("\nUsing CAV data and msm-software:\n")
cav <- cav[!is.na(cav$pdiag),]
subjects <- as.numeric(names(table(cav$PTNUM)))
N <- length(subjects)
cat("Sample size =",N,"\n")
cat("\nFrequencies observed history of state:")
print(table(cav$state))
cat("\nState table:")
print(statetable.msm(state,PTNUM,data=cav))

# Add baseline age:
bage <- rep(NA,nrow(cav))
for(i in 1:N){
  select.i <- cav$PTNUM==subjects[i]
  bage[select.i] <- cav$age[select.i][1]
}
cav <- cbind(cav,bage=bage)

#########################################
# Fit model using code in Jackson (2011):
# MC model:
ematrix <- rbind(c(0, 0.1, 0, 0), c(0.1, 0, 0.1, 0), 
                 c(0, 0.1, 0, 0), c(0, 0, 0, 0))
rownames(ematrix) <- c("Well", "Mild", "Severe","Death")
colnames(ematrix) <- c("Well", "Mild", "Severe","Death")

# Msm model:
oneway4.q <- rbind(c(0, 0.1, 0, 0.04), c(0, 0, 0.3, 0.05),
                   c(0, 0, 0, 0.3), c(0, 0, 0, 0))
rownames(oneway4.q) <- c("Well", "Mild", "Severe","Death")
colnames(oneway4.q) <- c("Well", "Mild", "Severe","Death")

# Fit model:
misc.msm <- msm(state ~ years, subject = PTNUM, data = cav,
    qmatrix = oneway4.q, ematrix = ematrix, obstrue = firstobs,
     death = TRUE, method = "BFGS")

# Generate output in user-written format:
model  <- misc.msm
qnames <- c("q12","q14","q23","q24","q34")
cat("\nModel in Jackson's 2011 paper:")
cat("\n-2loglik =", model$minus2loglik,"\n")
cat("AIC =", model$minus2loglik+2*length(model$opt$par),"\n")
conv <- model$opt$convergence
cat("Convergence code =", conv,"\n")
index <- c(1:5,9,14,16,21)
p    <- model$estimates[index]
p.se <- sqrt(diag(model$covmat))[index]
print(cbind(q=c(qnames,rep("c",4)),p=round(p,digits),
            se=round(p.se,digits),"Wald ChiSq"=round((p/p.se)^2,digits),
            "Pr>ChiSq"=round(1-pchisq((p/p.se)^2,df=1),digits)),quote=FALSE)
# Transformed intercepts:
cat("\nTransformed parameters as in paper:\n")
print(cbind(lambda=qnames,round(exp(p[1:5]),digits)),quote=FALSE)
# Misclassification matrix:
cat("\nMisclassification matrix as in paper:\n")
E.msm <- ematrix.msm(model, covariates="mean", ci="delta")
print(E.msm)

#################################################
# Extended model (as in Chapter 1, but with MC):
# Covariates:
covariates <- as.formula("~years+bage+dage")
fixedpars  <- c( 8:10,13:15)

# Fit model:
misc.msm2 <- msm(state ~ years, subject = PTNUM, data = cav, center=FALSE,
    qmatrix = oneway4.q, ematrix = ematrix, obstrue = firstobs,
    covariates=covariates, fixedpars=fixedpars,
    death = TRUE, method = "BFGS")

# Generate output in user-written format:
model <- misc.msm2
cat("\nExtended model:")
cat("\n-2loglik =", model$minus2loglik,"\n")
cat("AIC =", model$minus2loglik+2*length(model$opt$par),"\n")
conv <- model$opt$convergence
cat("Convergence code =", conv,"\n")
index <- c(1:20,24,29,31,36)
p    <- model$estimates[index]
p.se <- sqrt(diag(model$covmat))[index]
print(cbind(q=c(rep(qnames,4),rep("c",4)),p=round(p,digits),
            se=round(p.se,digits),"Wald ChiSq"=round((p/p.se)^2,digits),
            "Pr>ChiSq"=round(1-pchisq((p/p.se)^2,df=1),digits)),quote=FALSE)
# Misclassification matrix:
cat("\nMisclassification matrix:\n")
E.msm <- ematrix.msm(model, covariates="mean", ci="delta")
print(E.msm)
# Transition probs for 1 year:
cat("\nTransition probs for 1 year\n")
bage0 <- median(cav$bage[cav$firstobs==1])
dage0 <- median(cav$dage[cav$firstobs==1])
pmat <-  pmatrix.msm(model, t=1, t1=0, 
          covariates=list(years=0,bage=bage0,dage=dage0), ci="normal")
print(pmat)