# Chapter 1: CAV-history data analysis
# Using msm() for parameters estimation

# Ardo, 2016

# Prelim:
library(msm)
digits <- 3

# Rename data, and give info:
dta <- cav
subjects <- as.numeric(names(table(dta$PTNUM)))
N <- length(subjects)
cat("Sample size =",N,"\n")
cat("Frequencies observed state:"); print(table(dta$state))
cat("State table:"); print(statetable.msm(state,PTNUM,data=dta))

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

# Create history-state data:
Ostate <- dta$state
for(i in 1:N){
  dta.i <- dta[dta$PTNUM==subjects[i],]
  Hstate.i <- rep(NA, nrow(dta.i))
  for(j in 1:nrow(dta.i)){
    Hstate.i[j] <- max(dta.i$state[1:j])
  }
  dta$state[dta$PTNUM==subjects[i]] <- Hstate.i
}
dta <- cbind(dta,Ostate=Ostate)
# History data info:
cat("\nFrequencies observed history of state:"); print(table(dta$state))
cat("State table:"); print(statetable.msm(state,PTNUM,data=dta))

# Choose model:
Model <- 1
# Model formulation:
if(Model==1){
  # Generator matrix Q:
  q <- 0.01
  Q <- rbind(c(0,q,0,q), c(0,0,q,q),c(0,0,0,q),c(0,0,0,0))
  qnames <- c("q12","q14","q23","q24","q34")
  # Covariates:
  covariates <- as.formula("~years+bage+dage")
  constraint <- NULL
  fixedpars  <- c( 8:10,13:15)
  # Control:
  method <- "BFGS"
}

# Fit model using msm:
model <- msm(state~years, subject=PTNUM, data=dta, center=FALSE, 
            qmatrix=Q, death=TRUE, covariates=covariates, constraint=constraint,
            fixedpars=fixedpars, method=method, control=list(trace=0,REPORT=1))

# Generate output:
cat("\nModel",Model," with covariates: "); print(covariates)
cat("and constraints:\n"); print(constraint)
cat("and fixedpars:\n"); print(fixedpars)
cat("\n-2loglik =", model$minus2loglik,"\n")
conv <- model$opt$convergence; cat("Convergence code =", conv,"\n")
p    <- model$estimates; p.se <- sqrt(diag(model$covmat))
print(cbind(q=qnames,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)

# Transition probs for 1 year:
cat("\nTransition probs for 1 year\n")
bage0 <- median(dta$bage[dta$firstobs==1])
dage0 <- median(dta$dage[dta$firstobs==1])
pmat  <- pmatrix.msm(model, t=1, t1=0, covariates=list(years=0,
                    bage=bage0,dage=dage0), ci="normal")
print(pmat)