Hidden Markov Model

Hidden Markov Model for Disease Progression Measured with Error

A hidden Markov model joins two processes, a Markov chain of true disease states that nobody sees directly and a sequence of recorded results, such as screening grades, biomarker values or coded diagnoses, that each state produces with some probability. In health economics it is used to estimate how a chronic disease progresses when the only evidence is a series of imperfect measurements taken at clinic visits or screening rounds. Fitting one separates true progression from measurement error, so the transition probabilities passed to a Markov cohort model describe the disease and not the noise in the test. It also gives the probability that a given patient is in each state, which is what a screening or monitoring policy acts on.

Two layers in one model

Rabiner's tutorial calls the model a doubly embedded stochastic process, with an underlying process that is hidden and can only be observed through another set of stochastic processes that produce the observations. The hidden layer is an ordinary Markov chain over health states, with a starting distribution and a transition matrix. The observed layer is an emission distribution for each state, giving the probability of each possible result when the patient is truly in that state.

$$\Pr(X, Y) = \pi_{x_1} \, b_{x_1}(y_1) \prod_{t=2}^{T} a_{x_{t-1} x_t} \, b_{x_t}(y_t)$$

where $X = (x_1, \dots, x_T)$ is the sequence of hidden states, $Y = (y_1, \dots, y_T)$ is the sequence of observed results, $T$ is the number of observation times, $\pi_i$ is the probability of starting in state $i$, $a_{ij}$ is the probability that a person in state $i$ is in state $j$ at the next observation, and $b_i(y)$ is the emission probability of result $y$ in state $i$.

When the results are themselves states, the emission matrix is a misclassification matrix. The msm package manual describes this case for screening, where the probability that observed and true states agree is the sensitivity for a disease state and the specificity for a disease-free one. Emissions can also be continuous, as with FEV1 after lung transplantation, where the manual notes that short-term fluctuation makes the exact disease stage hard to determine at any visit. It also describes hidden Markov models as mixture models whose mixing distribution changes over time with the hidden chain. The model treats each recorded state as a possible misclassification of the true one.

Evaluation, decoding and learning

Rabiner organised the theory around three problems. Evaluation computes the probability of an observed sequence under a given model, the likelihood used to fit and compare models. Decoding chooses the hidden state sequence that best explains the observations, and learning adjusts the parameters to make the observations most probable.

Evaluation means summing over every path through the hidden states. Rabiner notes that direct calculation sums over $N^T$ paths and needs on the order of twice $T \cdot N^T$ calculations, against about $N^2 T$ for the forward procedure, which works through the observations one at a time. With 5 states and 100 observations the forward procedure needs about 3,000 calculations.

$$\alpha_1(j) = \pi_j \, b_j(y_1), \qquad \alpha_t(j) = \Big[ \sum_{i=1}^{N} \alpha_{t-1}(i) \, a_{ij} \Big] b_j(y_t), \qquad \Pr(Y) = \sum_{j=1}^{N} \alpha_T(j)$$

where $\alpha_t(j)$ is the joint probability of the first $t$ observations and of being in hidden state $j$ at time $t$, $N$ is the number of hidden states, and the other symbols are as above. Dividing $\alpha_t(j)$ by its sum over $j$ gives the probability that the patient is in state $j$ at time $t$ given the results so far. In a spreadsheet each step is a matrix multiplication by the transition matrix followed by element-wise multiplication by the emission probabilities. Rabiner warns that the terms shrink exponentially with long sequences, so practical code rescales them at each step, and Viterbi decoding can be run on logarithms.

For decoding, Rabiner points out that choosing the individually most likely state at each time can produce a sequence that is impossible when some transitions have zero probability. The Viterbi algorithm instead finds the single best path by dynamic programming. For learning, the Baum-Welch method is a form of the EM algorithm. Rabiner states that it reaches local maxima only and that the likelihood surface usually has many of them, so starting values matter.

Continuous time and irregular clinic visits

Most estimation theory was developed for equally spaced discrete time, but the msm manual notes that disease evolves in continuous time and patients are monitored at irregular and differing intervals. The hidden chain is then specified by a transition intensity matrix, and the transition probabilities over any gap between visits follow from it.

$$P(t) = \exp(tQ)$$

where $Q$ is the matrix of transition intensities between hidden states, $t$ is the time between two consecutive observations, $\exp$ is the matrix exponential, and $P(t)$ is the matrix of transition probabilities over that interval. Each patient's likelihood is then a product of matrices, one per visit, each combining $P(t)$ for that gap with the emission probabilities of the result recorded at that visit, which is the forward procedure with a different transition matrix at every step. The same conversion gives the cycle probabilities of a cohort model's transition matrix, and transitions known only to fall between visits are interval censored.

Jackson and Sharples used this approach for FEV1 after lung transplantation, estimating disease progression rates together with the effects of mismeasurement and covariates by maximum likelihood.

Worked example with an imperfect annual test

All numbers in this example are illustrative. A progressive condition is screened once a year. True states are Healthy (H) and Diseased (D), with no recovery, an annual onset probability of 0.10, and a starting distribution of 0.90 Healthy and 0.10 Diseased. The test has sensitivity 0.80 and specificity 0.90.

1. Set out the two layers. The transition matrix has rows (0.90, 0.10) for Healthy and (0, 1) for Diseased. The emission probabilities are 0.90 for a negative and 0.10 for a positive when Healthy, and 0.20 for a negative and 0.80 for a positive when Diseased.

2. Compare true and recorded prevalence. In year 1 true prevalence is 0.10, but the share testing positive is 0.10 × 0.80 + 0.90 × 0.10 = 0.17. In year 2 true prevalence is 0.90 × 0.10 + 0.10 = 0.19 and the positive share is 0.19 × 0.80 + 0.81 × 0.10 = 0.233.

3. Estimate transitions naively from test results. The probability of a negative in year 1 is 0.90 × 0.90 + 0.10 × 0.20 = 0.83, and the probability of a negative followed by a positive is 0.1537, so the apparent onset probability is 0.1537 / 0.83 = 0.185. The probability of a positive in year 1 is 0.17 and of a positive followed by a negative 0.0907, giving an apparent annual recovery of 0.0907 / 0.17 = 0.534. A model built on recorded states would nearly double onset and invent recovery from a condition that cannot remit.

4. Run the forward procedure for one patient. The patient tests negative, then positive, then negative.

Each prediction multiplies the previous forward vector by the transition matrix, for example 0.81 × 0.10 + 0.02 × 1 = 0.101, and each forward value multiplies the prediction by the emission probability, for example 0.101 × 0.80 = 0.0808. The likelihood of the three results is 0.059049 + 0.017618 = 0.076667, and the probability of disease after year 3 is 0.017618 / 0.076667 = 0.230.

5. Check by listing the paths and decode. Only four paths are possible. HHH has probability 0.81 × 0.09 × 0.81 = 0.059049, HHD has 0.0729 × 0.10 × 0.20 = 0.001458, HDD has 0.81 × 0.10 × 0.80 × 0.20 = 0.01296, and DDD has 0.02 × 0.80 × 0.20 = 0.0032. They sum to 0.076667, matching the forward result. The Viterbi path is HHH, with posterior probability 0.059049 / 0.076667 = 0.770, so the year-2 positive is most probably a false positive. Using all three results, the probability of disease in year 2 is (0.01296 + 0.0032) / 0.076667 = 0.211, well below the 0.526 available at the end of year 2.

6. Interpret. A clinician reading the record literally would say the patient had the disease in year 2 and then recovered. The model gives a 77% chance that the patient was never diseased and a 23% chance of disease now, the figure a referral decision should use. For a cost-effectiveness analysis, the true onset probability of 0.10 and the test accuracy enter the model as separate parameters, so a new test or a new screening interval can be evaluated without changing the natural history.

Fitted transitions as inputs to cost-effectiveness models

The usual route into economic evaluation has two stages. A hidden Markov model fitted to longitudinal screening or monitoring data estimates progression between true states, and a decision model adds costs, utilities and the strategies compared. Scotland and colleagues fitted a continuous-time hidden Markov model to national longitudinal data from the Scottish diabetic retinopathy screening programme and used the transition probabilities in a 30-year decision model. Annual against biennial screening cost roughly GBP 74,000 to GBP 232,000 per QALY gained across groups with no retinopathy, and roughly GBP 22,000 to GBP 85,000 in type 1 diabetes, so biennial screening for people with no retinopathy was judged likely to save money for a very small increase in risk.

Oke and colleagues built the first stage of this route for kidney function monitoring, where trials of monitoring strategies are lacking. Using electronic health records from 643 UK general practices in the Clinical Practice Research Datalink, they fitted hidden Markov models for eGFR stage progression and estimated that between 12.1% and 14.7% of tests misclassified the stage, so apparent change between consecutive results may be misclassification.

Assumptions and modelling choices

The likelihood rests on two assumptions that the msm manual states explicitly. Observations are conditionally independent given the underlying states, and the hidden process has the Markov property, so the next true state depends only on the current one. The second is the usual Markov assumption of a cohort model, applied to the true states. Results correlated for other reasons, such as repeated readings by one grader, break the first.

Identifiability is the main practical limit. The manual warns that over-complex models fitted to insufficient data will not be identifiable, recommends running models from several sets of initial values, and suggests constraints such as fixing misclassification between non-adjacent stages at zero or restricting state-specific means to plausible ranges. The number and type of hidden states are a further choice, for which Rabiner found no simple, theoretically correct rule. Each choice changes the transition probabilities, so a decision model should treat alternative specifications as structural uncertainty and carry parameter uncertainty into probabilistic sensitivity analysis.

Boundaries with neighbouring models

A Markov model for cost-effectiveness analysis assumes the modelled states are known at each cycle, and a hidden Markov model is often the step that estimates its inputs. A state space model has the same two layers with general, often continuous, states, and the Kalman filter handles its linear Gaussian case. A hidden Markov model has a finite set of discrete states. Latent transition analysis applies a similar structure to latent classes measured by several categorical indicators. When the decision itself depends on the current belief about the hidden state, the model becomes a partially observable Markov decision process, which Dorali and colleagues used to set personalised diabetic retinopathy screening recommendations.

Common misreadings

The most common error is treating the decoded path as observed data. The Viterbi path is the single most probable sequence, and in the example it still carries a 23% chance of being wrong, so state-specific costs and QALYs should be weighted by state probabilities. A second error is fitting an ordinary multi-state model to test results and reading apparent recoveries as real, which inflates both onset and remission. A third is reusing published sensitivity and specificity as the emission matrix without checking that the reference standard and population match. A fourth is reporting one fit without checking for other local maxima.

Sources

• Dorali P, Shahmoradi Z, Weng CY, Lee T. Cost-effectiveness analysis of a personalized, teleretinal-inclusive screening policy for diabetic retinopathy via Markov modeling. Ophthalmology Retina. 2023;7(6):532-542. https://doi.org/10.1016/j.oret.2023.01.001
• Jackson CH. Multi-state modelling with R: the msm package, version 1.8.2. MRC Biostatistics Unit; 2024. https://cran.r-project.org/web/packages/msm/vignettes/msm-manual.pdf
• Jackson CH, Sharples LD. Hidden Markov models for the onset and progression of bronchiolitis obliterans syndrome in lung transplant recipients. Statistics in Medicine. 2002;21(1):113-128. https://doi.org/10.1002/sim.886
• Oke JL, Feakins BG, Schlackow I, Mihaylova B, Simons C, O'Callaghan CA, et al. Statistical models for the deterioration of kidney function in a primary care population: a retrospective database analysis. F1000Research. 2019;8:1618. https://doi.org/10.12688/f1000research.20229.2
• Rabiner LR. A tutorial on hidden Markov models and selected applications in speech recognition. Proceedings of the IEEE. 1989;77(2):257-286. https://doi.org/10.1109/5.18626
• Scotland G, McKeigue P, Philip S, Leese GP, Olson JA, Looker HC, et al. Modelling the cost-effectiveness of adopting risk-stratified approaches to extended screening intervals in the national diabetic retinopathy screening programme in Scotland. Diabetic Medicine. 2016;33(7):886-895. https://doi.org/10.1111/dme.13129