---
title: "Advanced phase-type models in msmbayes"
author: "Christopher Jackson <chris.jackson@mrc-bsu.cam.ac.uk>"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_document:
    toc: true
    toc_float: true
    theme: simplex
    number-sections: true
resource_files:
  - ../man/figures/twostatephase.png
  - ../man/figures/threestatephase.png
  - ../man/figures/twostate.png
vignette: >
  %\VignetteIndexEntry{Advanced phase-type models}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---



To fit a model with non-exponential sojourn distributions (semi-Markov
model) in `msmbayes`, the easiest way is to use phase-type
approximations to Weibull and Gamma distributions, as explained in the
[vignette about this technique](semimarkov.html)  These represent the sojourn distribution with two parameters - which may be all that is possible to identify from data
where the state is only intermittently observed.

A more advanced approach to semi-Markov modelling in `msmbayes` is to
specify a phase-type distribution directly.



# Phase-type semi-Markov models

In a multi-state model, suppose we do not believe the sojourn time in a particular state follows an exponential distribution.  For example, in the infection model, suppose we do not believe that the time spent in the "test positive" state is exponential. 

```{r,echo=FALSE,out.width="500px"}
#| fig.alt: >
#|   A two-state Markov model
knitr::include_graphics("../man/figures/twostate.png")
```

We can relax this assumption by building a **"semi-Markov"** model, where the sojourn time in this state follows a more complex distribution.    One way to do this is to replace this state by a series of two (or more) latent or hidden states, known as "phases".  For example, with two phases: 

```{r,echo=FALSE,out.width="600px"}
#| fig.alt: >
#|   A phase-type model with two observable states, one of which has two phases
knitr::include_graphics("../man/figures/twostatephase.png")
```

> In a phase-type model, we allow progression from one phase to the next, but not transitions from later to earlier phases (or jumps between non-adjacent phases, if there are 3 phases or more).    Otherwise, we allow the same transitions out of each phase as were allowed in the original model (to state 1, in this example)

We then assume transitions between and out of the phases follow exponential distributions.  That is, we assume a Markov model on the latent state space.  The test-positive state then has a sojourn distribution known as the ["Coxian" phase-type](https://en.wikipedia.org/wiki/Phase-type_distribution#Coxian_distribution) distribution, instead of the exponential.

> **Note**: This is an example of a **hidden Markov model**, though one of a specific form where some states are observed correctly, and some states are latent.  Phase-type models are a convenient way to build semi-Markov models for intermittently-observed data, because there are standard algorithms for computing the likelihood of a hidden Markov model.



# Fitting phase-type models in msmbayes

To assume a phase-type sojourn distribution for one or more states in a `msmbayes`
model, set the `nphase` argument.  This is a vector of length equal to the number of states, giving the number of phases per state.  Any elements of `nphase` that are 1 correspond to states with the usual exponential sojourn distribution.


## Example: data simulated from a standard Markov model

Here we extend the 2-state model for the simulated infection data to give state 2 (infection) a two-phase sojourn distribution.   These data were originally simulated assuming an exponential sojourn distribution in state 2 (with rate 3).   In this situation, we would expect the estimated rates of transition out of each phase to be identical.

```{r}
library(msmbayes)
```

```{r,results="hide"}
Q <- rbind(c(0, 1), 
           c(1, 0))
draws <- msmbayes(infsim2, state="state", time="months", subject="subject",
                  Q=Q, nphase=c(1,2), fit_method = "optimize")
```

```{r}
summary(draws)
```

The phased states are labelled specially here, e.g. the first phase of state 2 is labelled `"2p1"`. The estimated transition rates from the two phases of state 2 to state 1 are not significantly different from 2 - but there is large uncertainty around the estimated rates, since this model is over-fitted to the data. Implementing the same model in `msm` would give extremely large confidence intervals around the estimated rates.


## Example: data simulated from a phase-type model 

The following code simulates a dataset from the following phase-type model structure 

```{r,echo=FALSE,out.width="600px"}
#| fig.alt: >
#|   A phase-type model with three observable states, where the first state has two phases
knitr::include_graphics("../man/figures/threestatephase.png")
```

with a transition intensity matrix (on the latent state space) given in the R object `Qg`.  This code uses the ability of the function `simmulti.msm()` from the `msm` package to simulate from hidden Markov models.
```{r}
Qg <- rbind(c(0, 0.18, 0.008, 0.012),
            c(0, 0,    0.016, 0.024),
            c(0, 0,    0,     0.2),
            c(0, 0,    0,     0))
E <- rbind(c(1,0,0,0), # hidden Markov model misclassification matrix
           c(1,0,0,0),
           c(0,1,0,0),
           c(0,0,1,0))
nsubj <- 10000; nobspt <- 10
set.seed(1)
sim.df <- data.frame(subject = rep(1:nsubj, each=nobspt),
                     time = seq(0, 100, length=nobspt))
library(msm)
sim.df <- simmulti.msm(sim.df[,1:2], qmatrix=Qg, ematrix=E)
```

We fit the phase-type model to the simulated data using both `msmbayes` and `msm`.

The fast optimization-based method is used to approximate the Bayesian posterior here,
for speed of demonstration, but in practice we would generally want to cross-check the results with MCMC.
The posterior distributions for these models are likely to have a skewed and multimodal 
shape which is not represented well by a normal approximation.

```{r,results="hide",warning=FALSE,error=FALSE,message=FALSE}
#| fig.alt: >
#|   Posterior density estimates
Q3 <- rbind(c(0,1,1),c(0,0,1),c(0,0,0))
draws <- msmbayes(data=sim.df, state="obs", time="time", subject="subject",
                  Q=Q3, nphase=c(2,1,1), fit_method="optimize")
bayesplot::mcmc_dens(draws, pars=sprintf("logq[%s]",1:6))
```


`msm` needs some tuning to converge, in particular the use of `fnscale` and 
explicit initial values for the transition rate.  In practice, sensitivity analysis
to these initial values would be needed to confirm that the reported estimates are
the global maximum likelihood estimates, rather than one of the local maxima or
saddle points.

```{r,results="hide"}
Q3 <- rbind(c(0,0.5,0.5),c(0,0,0.5),c(0,0,0))
s.msm <- msm(obs ~ time, subject=subject, data=sim.df, phase.states=1, qmatrix=Q3,
             phase.inits=list(list(trans=0.05,  # seems to need these
                                   exit=matrix(rep(0.05,4),nrow=2,byrow=TRUE))),
             control = list(trace=1,REPORT=1,fnscale=50000,maxit=10000))
```

Reassuringly, the Bayesian and frequentist methods give similar estimates of the transition intensities, which agree (within estimation error) with the values used for simulation. 

```{r}
qmatrix(draws)
qmatrix.msm(s.msm,ci="none")
Qg
```

The mean sojourn times in states of a phase-type model can either be calculated for the observable state space (`states="obs"`, the default), or the latent state space (`states="phase"`). 
```{r}
mean_sojourn(draws)
mean_sojourn(draws, states="phase") 
```

More practical experience of using general phase-type models for intermittently-observed multi-state data is needed.  In particular:

* How to specify substantive prior information.  In the examples above, the default
log-normal(-2,2) priors were used for transition intensities.  But we need a better 
way to choose these based on judgements about interpretable quantities, e.g. the mean sojourn time in the phased state.

* Limited experience suggests that the Bayesian method is more likely than `msm` to produce a plausible result without the need for tuning.   However, can the sampling algorithms always be relied upon?

* How to build and interpret models with covariates on intensities.  Covariates can be placed on intensities for transitions in the phased state-space, but this has not been tested in practice.  There is a risk of over-fitting, and strong constraints on covariate effects are expected to be necessary.

However, with intermittently-observed data, it is uncertain to what extent we will be able to identify a detailed shape for the sojourn distribution from the data.
Therefore using [phase-type approximations to shape-scale distributions](semimarkov.html) is expected to be a more useful approach to semi-Markov modelling.
