The most probable RECIST-like path when visits go missing
Oncology trials record tumour response on a visit calendar: Baseline \(\rightarrow\) Week 4 \(\rightarrow\) Week 8 \(\rightarrow\) Week 12 \(\rightarrow\) Week 24. At each post-baseline assessment the category is one of CR, PR, SD, or PD.1 Scans are missed, windows are missed, and a visit is coded NE.2 The table then looks like
\[ \text{SD} \;\rightarrow\; \text{PR} \;\rightarrow\; {?} \;\rightarrow\; \text{PR} \;\rightarrow\; \text{CR}. \]
Last observation carried forward (LOCF) fills the hole with PR. Copying the neighbour on either side does the same. Neither method asks whether the completed chain is a plausible disease course.
A hidden Markov model (HMM) treats the recorded category as a noisy or incomplete view of a latent status \(s_t\) (Baum and Petrie, 1966; Rabiner, 1989). The Viterbi algorithm (Viterbi, 1967; Forney, 1973) returns the single most probable state sequence given the recordings — the MAP3 path, not a visit-wise majority vote.
Restore the chain \(s_1 \rightarrow \cdots \rightarrow s_T\) that maximises \(P(s_{1:T} \mid y_{1:T})\), subject to clinical transition structure (PD absorbing, improvement SD \(\to\) PR \(\to\) CR allowed) and to whatever was actually seen.
This article implements that decoder from scratch in R, on a five-visit RECIST-like grid, including missing visits, unequal spacing, and a short misclassification extension.
A RECIST caveat, stated first. Under RECIST 1.1, Baseline is the reference measurement, not a response category (Therasse et al., 2000; Eisenhauer et al., 2009). CR/PR/SD/PD are defined post-baseline against baseline and nadir; confirmation rules and the 20% + 5 mm PD criterion are operational, not Markov. The five-point chain below is a latent status process on the visit calendar — a reconstruction and QC tool. It is not a best overall response (BOR) engine and must not replace the protocol estimand (ICH E9(R1), 2019).
library(dplyr, warn.conflicts = FALSE)
library(tidyr)
library(purrr)
library(ggplot2)
library(expm)Loading required package: Matrix
Attaching package: 'Matrix'The following objects are masked from 'package:tidyr':
expand, pack, unpack
Attaching package: 'expm'The following object is masked from 'package:Matrix':
expmSuppose Week 8 is missing and the recorded sequence is SD, PR, NA, PR, CR. In practice the chain gets completed three ways:
| Method | Completed chain |
|---|---|
| LOCF | SD \(\to\) PR \(\to\) PR \(\to\) PR \(\to\) CR |
| Copy next | SD \(\to\) PR \(\to\) PR \(\to\) PR \(\to\) CR |
| Independent draw | whatever the marginal at Week 8 suggests |
LOCF and copy-next agree here only because the recordings flanking the hole are both PR; they part ways as soon as the next observation differs from the last.
LOCF is a poor primary analysis for longitudinal means (Siddiqui, Hung, and O’Neill, 2009; National Research Council, 2010; Little et al., 2012). For a categorical disease course it fails a second way: it ignores transition constraints. A sequence PD, NA, CR becomes PD, PD, CR under LOCF, which contradicts absorbing progression. Independent imputation of each visit can invent the same impossibility.
Rabiner (1989) separates three HMM jobs:
This post is job 2. Job 3 is out of scope: the transition matrix here is protocol-structured, not estimated from the same small trial we then “complete”.
The decoding estimand is a path. The forward–backward algorithm gives the visit-wise posteriors \[P(S_t = s \mid y_{1:T}),\] yet the visit-wise modes can disagree with the Viterbi path (Rabiner, 1989; Cappé, Moulines, and Rydén, 2005). For “what was the most plausible course of disease?” the path is the right object.
Four latent statuses, in clinical order of worsening. The plate below is the same five-visit calendar: a Markov chain of hidden statuses \(s_t\) and a recording \(y_t\) that may be missing.
%%{init: {"theme": "base", "themeVariables": {"darkMode": false, "background": "#ffffff", "primaryColor": "#d6eaf8", "primaryTextColor": "#1a252f", "primaryBorderColor": "#1f6f8b", "lineColor": "#2c3e50", "mainBkg": "#d6eaf8", "clusterBkg": "#f7f9fb", "clusterBorder": "#5d6d7e", "titleColor": "#1a252f"}, "flowchart": {"padding": 12, "nodeSpacing": 18, "rankSpacing": 40, "curve": "linear"}}}%%
flowchart LR
subgraph hidden [Latent status]
S0["s1"]
S1["s2"]
S2["s3"]
S3["s4"]
S4["s5"]
S0 --> S1 --> S2 --> S3 --> S4
end
subgraph rec [Recorded RECIST]
Y0["y1"]
Y1["y2"]
Y2["y3 missing"]
Y3["y4"]
Y4["y5"]
end
S0 --> Y0
S1 --> Y1
S2 --> Y2
S3 --> Y3
S4 --> Y4
classDef hid fill:#d6eaf8,stroke:#1f6f8b,color:#1a252f,stroke-width:1.5px
classDef seen fill:#d5f5e3,stroke:#1e8449,color:#1a252f,stroke-width:1.5px
classDef miss fill:#fdebd0,stroke:#b9770e,color:#1a252f,stroke-width:1.8px
class S0,S1,S2,S3,S4 hid
class Y0,Y1,Y3,Y4 seen
class Y2 miss
states <- c("CR", "PR", "SD", "PD")Visits \(t = 1,\ldots,5\) sit at Baseline, Week 4, Week 8, Week 12, and Week 24. The latent process is a Markov chain
\[ P(s_1) = \pi(s_1), \qquad P(s_t \mid s_{1:t-1}) = P(s_t \mid s_{t-1}), \]
and the recording \(y_t\) depends only on \(s_t\) (memoryless emission). That is the classical HMM (Baum and Petrie, 1966).
At Baseline, eligibility usually forbids overt PD, CR is rare, and mass sits on SD (with some PR). A teaching prior:
pi0 <- c(CR = 0.05, PR = 0.20, SD = 0.70, PD = 0.05)
pi0 CR PR SD PD
0.05 0.20 0.70 0.05 These numbers are expert structure, not MLE.
Rows are origin \(s'\), columns destination \(s\), so the entry is \(P(s \mid s')\). PD is absorbing: once progressed, stay PD. Improvement SD \(\to\) PR \(\to\) CR is allowed; CR \(\to\) PR is the confirmation / relapse leak, not a typical recovery from PD.
P <- matrix(
c(
0.70, 0.20, 0.08, 0.02,
0.15, 0.60, 0.20, 0.05,
0.02, 0.18, 0.70, 0.10,
0.00, 0.00, 0.00, 1.00
),
nrow = 4, byrow = TRUE,
dimnames = list(from = states, to = states)
)
P to
from CR PR SD PD
CR 0.70 0.20 0.08 0.02
PR 0.15 0.60 0.20 0.05
SD 0.02 0.18 0.70 0.10
PD 0.00 0.00 0.00 1.00rowSums(P)CR PR SD PD
1 1 1 1 | From \(\backslash\) To | CR | PR | SD | PD |
|---|---|---|---|---|
| CR | 0.70 | 0.20 | 0.08 | 0.02 |
| PR | 0.15 | 0.60 | 0.20 | 0.05 |
| SD | 0.02 | 0.18 | 0.70 | 0.10 |
| PD | 0.00 | 0.00 | 0.00 | 1.00 |
The last hop on the calendar is 12 weeks, not 4. The matrix above is a 4-week step. Week 12 \(\to\) Week 24 is handled later with a generator \(Q\) and \(P(\Delta t) = \exp(Q\,\Delta t)\) (Kalbfleisch and Lawless, 1985). Until then, treat \(P\) as the one-step kernel on the discrete visit index, which is already enough to see Viterbi work.
The teaching model makes the recording exact:
Schwartz et al. (2016) clarify NE / unevaluable under RECIST 1.1; here NE is simply “no emission constraint”. The Noisy emissions section below replaces the indicator with a confusion matrix \(C_{s,k} = P(y=k \mid s)\) (Jackson et al., 2003).
emission_hard <- function(y) {
e <- setNames(rep(0, length(states)), states)
if (is.na(y) || y %in% c("NE", "MISSING")) {
e[] <- 1
} else if (!y %in% states) {
stop("unknown category: ", y)
} else {
e[y] <- 1
}
e
}
emission_hard("PR")CR PR SD PD
0 1 0 0 emission_hard(NA_character_)CR PR SD PD
1 1 1 1 There are \(4^5 = 1024\) chains on five visits — enumerable, but the point of the algorithm is the recurrence that stays \(O(T |S|^2)\) when \(T\) or \(|S|\) grows (Durbin, Eddy, Krogh, and Mitchison, 1998).
Let \(dp_t(s)\) be the probability of the most probable chain that ends in \(s\) at visit \(t\) and that accounts for \(y_{1:t}\). Then (Viterbi, 1967; Forney, 1973; Rabiner, 1989)
\[ dp_1(s) = \pi(s)\, e_1(s), \]
\[ dp_t(s) = e_t(s) \cdot \max_{s'} \bigl[ dp_{t-1}(s') \, P(s \mid s') \bigr], \]
\[ \psi_t(s) = \arg\max_{s'} \bigl[ dp_{t-1}(s') \, P(s \mid s') \bigr]. \]
At \(T\), \(s_T^\star = \arg\max_s dp_T(s)\). Backpointers \(\psi\) reconstruct \(s_{T-1}^\star, \ldots, s_1^\star\).
%%{init: {"theme": "base", "themeVariables": {"darkMode": false, "primaryColor": "#d6eaf8", "primaryTextColor": "#1a252f", "primaryBorderColor": "#1f6f8b", "lineColor": "#2c3e50", "mainBkg": "#d6eaf8", "titleColor": "#1a252f"}, "flowchart": {"padding": 12, "nodeSpacing": 18, "rankSpacing": 40, "curve": "linear"}}}%%
flowchart LR
prev["dp at t-1 over s'"]
trans["times P(s given s')"]
mx["max over s'"]
em["times e_t(s)"]
now["dp_t(s)"]
prev --> trans --> mx --> em --> now
classDef step fill:#d6eaf8,stroke:#1f6f8b,color:#1a252f,stroke-width:1.5px
class prev,trans,mx,em,now step
Products of probabilities underflow. The implementation uses logarithms (Rabiner, 1989):
\[ \log dp_t(s) = \log e_t(s) + \max_{s'} \bigl[ \log dp_{t-1}(s') + \log P(s \mid s') \bigr]. \]
If every candidate is impossible, \(\log dp_T = -\infty\): the recorded sequence is incompatible with \(P\) (typically PD then CR under an absorbing kernel).
Take \(y = (\mathrm{SD},\; \mathrm{NA},\; \mathrm{PR})\) and the matrices above. Visit 1 is pinned to SD:
\[ \log dp_1(\mathrm{SD}) = \log 0.70, \qquad \log dp_1(\mathrm{CR}) = \log dp_1(\mathrm{PR}) = \log dp_1(\mathrm{PD}) = -\infty. \]
Visit 2 is missing, so \(e_2 \equiv 1\) and the only surviving predecessor is SD:
\[ \log dp_2(s) = \log 0.70 + \log P(s \mid \mathrm{SD}). \]
Visit 3 is pinned to PR, so only \(s=\mathrm{PR}\) is legal:
\[ \log dp_3(\mathrm{PR}) = \max_{s'} \bigl[ \log dp_2(s') + \log P(\mathrm{PR} \mid s') \bigr]. \]
log0 <- function(x) {
out <- log(x)
out[!is.finite(out) | x <= 0] <- -Inf
out
}
logP <- log0(P)
dp1 <- log0(pi0) + log0(emission_hard("SD"))
dp2 <- sapply(states, function(s) {
max(dp1 + logP[, s])
})
dp3 <- log0(emission_hard("PR")) + sapply(states, function(s) {
max(dp2 + logP[, s])
})
round(cbind(dp1 = dp1, dp2 = dp2, dp3 = dp3), 3) dp1 dp2 dp3
CR -Inf -4.269 -Inf
PR -Inf -2.071 -2.428
SD -0.357 -0.713 -Inf
PD -Inf -2.659 -InfThe surviving path is forced at the ends (SD, ?, PR). The missing middle is the \(s'\) that attains the max into PR. Staying in SD at the hole wins: \(P(\mathrm{SD} \mid \mathrm{SD}) = 0.70\) so \(\log dp_2(\mathrm{SD})\) is far above \(\log dp_2(\mathrm{PR})\), and that occupancy beats \(P(\mathrm{PR} \mid \mathrm{PR}) = 0.60\) versus \(P(\mathrm{PR} \mid \mathrm{SD}) = 0.18\). The Viterbi chain is SD \(\to\) SD \(\to\) PR — the same fill LOCF would have written, for a different reason (the kernel, not a copy rule). A more sticky \(P(\mathrm{PR} \mid \mathrm{PR})\) would flip the backpointer.
psi3_into_PR <- states[which.max(dp2 + logP[, "PR"])]
psi3_into_PR[1] "SD"viterbi_tv <- function(y, P_list, pi0, emission) {
Tm <- length(y)
stopifnot(length(P_list) == Tm - 1L)
st <- colnames(P_list[[1]])
K <- length(st)
dp <- matrix(-Inf, Tm, K, dimnames = list(NULL, st))
psi <- matrix(NA_integer_, Tm, K)
dp[1, ] <- log0(pi0) + log0(emission(y[1]))
if (Tm > 1L) {
for (t in 2:Tm) {
logP <- log0(P_list[[t - 1L]])
loge <- log0(emission(y[t]))
for (s in seq_len(K)) {
scores <- dp[t - 1L, ] + logP[, s]
psi[t, s] <- which.max(scores)
dp[t, s] <- loge[s] + max(scores)
}
}
}
log_prob <- max(dp[Tm, ])
if (!is.finite(log_prob)) {
return(list(
path = rep(NA_character_, Tm),
log_prob = -Inf,
dp = dp,
psi = psi,
feasible = FALSE
))
}
path_idx <- integer(Tm)
path_idx[Tm] <- which.max(dp[Tm, ])
if (Tm > 1L) {
for (t in Tm:2) {
path_idx[t - 1L] <- psi[t, path_idx[t]]
}
}
list(
path = st[path_idx],
log_prob = log_prob,
dp = dp,
psi = psi,
feasible = TRUE
)
}
viterbi <- function(y, P, pi0, emission) {
viterbi_tv(y, rep(list(P), length(y) - 1L), pi0, emission)
}The recurrence is written once, against a list of one-step kernels — one per hop, P_list[[t - 1]] carrying hop \(t - 1 \to t\). The homogeneous decoder is the special case of feeding the same \(P\) to every hop. We will cash in that generality in the Unequal spacing section.
feasible = FALSE is the honest answer when the recordings contradict the kernel (hard PD then CR). which.max() on a row of \(-\infty\) would otherwise invent a traceback.
ADaM-style long data: USUBJID, VISIT, TIME, AVALC in {CR, PR, SD, PD, NA}. Subject SUBJ001 is the teaching sequence with a hole at Week 8. A few neighbours are simulated from the same \((\pi, P)\) with MAR-like missingness after Baseline.
set.seed(20260814)
visits <- tibble(
VISIT = factor(
c("Baseline", "Week 4", "Week 8", "Week 12", "Week 24"),
levels = c("Baseline", "Week 4", "Week 8", "Week 12", "Week 24")
),
TIME = c(0, 4, 8, 12, 24)
)
draw_chain <- function() {
s <- character(nrow(visits))
s[1] <- sample(states, 1L, prob = pi0)
for (t in 2:length(s)) {
s[t] <- sample(states, 1L, prob = P[s[t - 1L], ])
}
s
}
n_subj <- 8L
sim <- map_dfr(seq_len(n_subj), function(i) {
tibble(
USUBJID = sprintf("SUBJ%03d", i),
visits,
TRUE_S = draw_chain()
)
}) %>%
group_by(USUBJID) %>%
mutate(
miss = TIME > 0 & runif(n()) < 0.18,
AVALC = ifelse(miss, NA_character_, TRUE_S)
) %>%
ungroup()
# Teaching subject: SD → PR → (missing) → PR → CR
teach <- tibble(
USUBJID = "SUBJ001",
VISIT = visits$VISIT,
AVALC = c("SD", "PR", NA_character_, "PR", "CR"),
TRUE_S = c("SD", "PR", "PR", "PR", "CR")
)
sim <- sim %>%
rows_update(teach, by = c("USUBJID", "VISIT"))
sim %>%
filter(USUBJID == "SUBJ001") %>%
select(USUBJID, VISIT, TIME, AVALC)# A tibble: 5 × 4
USUBJID VISIT TIME AVALC
<chr> <fct> <dbl> <chr>
1 SUBJ001 Baseline 0 SD
2 SUBJ001 Week 4 4 PR
3 SUBJ001 Week 8 8 <NA>
4 SUBJ001 Week 12 12 PR
5 SUBJ001 Week 24 24 CR y001 <- sim %>%
filter(USUBJID == "SUBJ001") %>%
arrange(TIME) %>%
pull(AVALC)
vt001 <- viterbi(y001, P, pi0, emission_hard)
vt001$path[1] "SD" "PR" "PR" "PR" "CR"vt001$log_prob[1] -4.990245vt001$feasible[1] TRUEThe DP table (log-domain). \(-\)Inf marks states ruled out by a hard emission. At Week 8 every column is finite: the hole.
dp_tbl <- as.data.frame(round(vt001$dp, 3))
dp_tbl$VISIT <- levels(visits$VISIT)
dp_tbl$AVALC <- y001
dp_tbl %>%
select(VISIT, AVALC, all_of(states)) VISIT AVALC CR PR SD PD
1 Baseline SD -Inf -Inf -0.357 -Inf
2 Week 4 PR -Inf -2.071 -Inf -Inf
3 Week 8 <NA> -3.969 -2.582 -3.681 -5.067
4 Week 12 PR -Inf -3.093 -Inf -Inf
5 Week 24 CR -4.990 -Inf -Inf -InfWeek 8 is restored as PR: the most probable bridge from the Week-4 PR into the Week-12 PR, under this \(P\).
locf_fill <- function(y) {
last <- NA_character_
out <- y
for (i in seq_along(y)) {
if (!is.na(y[i])) last <- y[i]
out[i] <- last
}
out
}
compare_one <- function(y, P, pi0, emission) {
vt <- viterbi(y, P, pi0, emission)
tibble(
VISIT = levels(visits$VISIT),
TIME = visits$TIME,
observed = y,
viterbi = vt$path,
locf = locf_fill(y),
feasible = vt$feasible,
log_prob = vt$log_prob
)
}
cmp001 <- compare_one(y001, P, pi0, emission_hard)
cmp001 %>%
select(VISIT, observed, viterbi, locf)# A tibble: 5 × 4
VISIT observed viterbi locf
<chr> <chr> <chr> <chr>
1 Baseline SD SD SD
2 Week 4 PR PR PR
3 Week 8 <NA> PR PR
4 Week 12 PR PR PR
5 Week 24 CR CR CR On this particular hole LOCF and Viterbi agree. They need not. Replace the Week-12 recording with SD and Viterbi prefers an SD middle; LOCF still copies PR.
y_alt <- c("SD", "PR", NA, "SD", "SD")
cmp002 <- compare_one(y_alt, P, pi0, emission_hard)
cmp002 %>% select(VISIT, observed, viterbi, locf)# A tibble: 5 × 4
VISIT observed viterbi locf
<chr> <chr> <chr> <chr>
1 Baseline SD SD SD
2 Week 4 PR PR PR
3 Week 8 <NA> SD PR
4 Week 12 SD SD SD
5 Week 24 SD SD SD rank_of <- function(x) {
match(x, c("CR", "PR", "SD", "PD"))
}
plot_df <- bind_rows(
cmp002 %>%
transmute(VISIT, TIME, method = "Viterbi", state = viterbi),
cmp002 %>%
transmute(VISIT, TIME, method = "LOCF", state = locf),
cmp002 %>%
filter(!is.na(observed)) %>%
transmute(VISIT, TIME, method = "Observed", state = observed)
) %>%
mutate(
method = factor(method, levels = c("Observed", "Viterbi", "LOCF")),
y = rank_of(state)
)
ggplot(plot_df, aes(TIME, y, colour = method, shape = method)) +
geom_line(data = filter(plot_df, method != "Observed"), linewidth = 0.7) +
geom_point(size = 3) +
scale_x_continuous(
breaks = visits$TIME,
labels = as.character(visits$VISIT)
) +
scale_y_reverse(
breaks = 1:4,
labels = c("CR", "PR", "SD", "PD")
) +
scale_colour_manual(
values = c(Observed = "#222222", Viterbi = "#2C7FB8", LOCF = "#D95F0E")
) +
labs(
x = NULL, y = NULL, colour = NULL, shape = NULL,
title = "SUBJ001: observed recordings, Viterbi path, LOCF"
) +
theme_minimal(base_size = 11) +
theme(panel.grid.minor = element_blank())
Week 12 \(\to\) Week 24 is three times a 4-week step. Treating that hop like any other understates the chance the disease moved during the long gap. What exactly “treating it like any other” means comes in two flavors, and we will compare both against the honest one.
A continuous-time Markov chain (CTMC) describes the process with a generator \(Q\) instead of a step matrix: off-diagonal \(q_{ij} \ge 0\) is the instantaneous rate of jumping from state \(i\) to state \(j\) “per week”, the diagonal is \(-q_i\), the rate of leaving state \(i\) altogether, so rows sum to 0, and an absorbing PD is the all-zero row (\(Q_{\mathrm{PD},\cdot} = 0\)). The transition kernel over an interval \(\Delta t\) is then the matrix exponential
\[ P(\Delta t) = \exp(Q\,\Delta t) \]
(Kalbfleisch and Lawless, 1985; Jackson, 2011). The longer the visit gap, the further \(P(\Delta t)\) moves from the identity: more time to leave every state. Compute the exponential with a proper algorithm, not a truncated series (Moler and Van Loan, 2003); expm::expm() implements Pade with scaling and squaring.
A natural first move is to recover \(Q\) from the protocol \(P\) as a matrix logarithm, \(Q = \log(P)/4\), so that \(P(4) = P\) exactly, and then let \(P(12) = \exp(3 \cdot 4\,Q)\). This does not work here:
Ql <- as.matrix(logm(P)) / 4
Ql [,1] [,2] [,3] [,4]
[1,] -0.0979035783 0.07699224 0.01787038 0.003040954
[2,] 0.0598251562 -0.14913036 0.07772671 0.011578499
[3,] -0.0007348382 0.07203501 -0.09973973 0.028439562
[4,] 0.0000000000 0.00000000 0.00000000 0.000000000min(Ql[!diag(4)])[1] -0.0007348382The off-diagonal entry \(q_{\mathrm{SD},\mathrm{CR}} < 0\) disqualifies \(\log(P)/4\) as a generator: this \(P\) is a step kernel that is not embeddable in a CTMC — a 4-week step SD \(\to\) CR is allowed, but not as a rate of crossing SD to CR directly. (Embeddability is a real constraint, not a curiosity: a matrix can be a fine one-step kernel and still have no valid generator.) So \(Q\) below is its own expert structure, chosen so that \(P(4)\) resembles the protocol \(P\) — within \(\pm 0.05\) per entry — while rows sum to 0 and PD stays absorbing:
Q <- matrix(
c(
-0.10, 0.070, 0.020, 0.010,
0.040, -0.120, 0.060, 0.020,
0.005, 0.050, -0.080, 0.025,
0.000, 0.000, 0.000, 0.000
),
nrow = 4, byrow = TRUE,
dimnames = list(from = states, to = states)
)
rowSums(Q) CR PR SD PD
1.734723e-18 1.040834e-17 0.000000e+00 0.000000e+00 max(abs(P - as.matrix(expm(Q * 4))))[1] 0.04936662P4 <- as.matrix(expm(Q * 4))
P12 <- as.matrix(expm(Q * 12))
dimnames(P4) <- dimnames(P12) <- dimnames(Q)
round(P4, 3) to
from CR PR SD PD
CR 0.686 0.189 0.079 0.046
PR 0.106 0.649 0.168 0.076
SD 0.025 0.138 0.744 0.093
PD 0.000 0.000 0.000 1.000round(P12, 3) to
from CR PR SD PD
CR 0.370 0.284 0.191 0.155
PR 0.156 0.363 0.270 0.211
SD 0.070 0.218 0.467 0.245
PD 0.000 0.000 0.000 1.000Over 12 weeks the chain has more time to leave CR/PR/SD for PD; the PD row stays \((0,0,0,1)\). The cash-in from the Implementation section: viterbi_tv was built for a list of kernels, so the time-inhomogeneous decoder is just the call — \(P_t = P(\Delta t_t)\) at each hop instead of a single \(P\):
P_list <- list(P4, P4, P4, P12)
vt_tv <- viterbi_tv(y001, P_list, pi0, emission_hard)
vt_tv$path[1] "SD" "PR" "PR" "PR" "CR"Two shortcuts to the same 12-week hop are tempting, and both misplace probability relative to \(P(\Delta t)\):
For SUBJ001 all three decode the same hole as PR, so we add one chain where the choices disagree — SD, PR, --, --, PD — two holes, one long, and a terminal PD:
P3 <- P %*% P %*% P
dimnames(P3) <- dimnames(P)
spacing_table <- function(y) {
runs <- list(
`P once each hop` = list(P, P, P, P),
`P^3 on the long hop` = list(P, P, P, P3),
`exp(Q dt)` = list(P4, P4, P4, P12)
)
out <- map(runs, \(rl) viterbi_tv(y, rl, pi0, emission_hard))
tibble(
kernel = names(runs),
path = map_chr(out, \(o) paste0(o$path, collapse = " ")),
feasible = map_lgl(out, \(o) o$feasible),
log_prob = map_dbl(out, \(o) o$log_prob)
)
}
spacing_table(y001)# A tibble: 3 × 4
kernel path feasible log_prob
<chr> <chr> <lgl> <dbl>
1 P once each hop SD PR PR PR CR TRUE -4.99
2 P^3 on the long hop SD PR PR PR CR TRUE -4.66
3 exp(Q dt) SD PR PR PR CR TRUE -5.06spacing_table(c("SD", "PR", NA_character_, NA_character_, "PD"))# A tibble: 3 × 4
kernel path feasible log_prob
<chr> <chr> <lgl> <dbl>
1 P once each hop SD PR PD PD PD TRUE -5.07
2 P^3 on the long hop SD PR PR PR PD TRUE -4.95
3 exp(Q dt) SD PR PR PR PD TRUE -4.75Two things fall out of the second chain. \(P\)-once says progression already by Week 8 (SD PR PD PD PD); \(P^3\) and \(P(\Delta t)\) say stable, then PD at 24 — the naive kernel, which cannot move much per hop, pushes the move earlier and calls it done. And on the first chain the log-likelihoods separate: the same path is \(-4.99\) under \(P\)-once but \(-4.66\) under \(P^3\), so the shortcuts disagree about plausibility even when they agree about the path — which matters as soon as you score or compare chains.
The starkest disagreement is on CR, --, --, --, PD: \(P\)-once decodes it as CR PD PD PD PD (log-prob \(-6.91\)) — progression forced all the way by Week 4 — while \(P(\Delta t)\) decodes it as CR CR CR CR PD (\(-5.99\)): hold CR across the short hops and let the long gap do the moving. The single-4-week kernel, marched across a 12-week gap hop by hop, cannot say “the state is unknown because the interval was long”; it has to over-explain the movement somewhere. On irregular real calendars \(P(\Delta t)\) is the honest discrete skeleton of a continuous-time chain; msm fits that model from panel data (Jackson, 2011).
Independent review and local investigators disagree. A recorded PR can be a true CR or SD. Jackson et al. (2003) put that misclassification in the emission; Jackson and Sharples (2002) is the clinical HMM cousin (staged decline, noisy marker).
Let \(C_{s,k} = P(y = k \mid s)\), rows = true status, columns = recorded category:
C <- matrix(
c(
0.88, 0.09, 0.02, 0.01,
0.06, 0.82, 0.10, 0.02,
0.02, 0.08, 0.82, 0.08,
0.01, 0.03, 0.08, 0.88
),
nrow = 4, byrow = TRUE,
dimnames = list(true = states, obs = states)
)
rowSums(C)CR PR SD PD
1 1 1 1 emission_noisy <- function(y) {
e <- setNames(rep(0, length(states)), states)
if (is.na(y) || y %in% c("NE", "MISSING")) {
e[] <- 1
} else if (!y %in% states) {
stop("unknown category: ", y)
} else {
e[] <- C[, y]
}
e
}
emission_noisy("PR") CR PR SD PD
0.09 0.82 0.08 0.03 SD, PD, CR, CR, CR cannot occur if PD is absorbing and emissions are indicators: \(P(\mathrm{CR} \mid \mathrm{PD}) = 0\) and \(e(\mathrm{CR})\) kills every other state.
y_bad <- c("SD", "PD", "CR", "CR", "CR")
vt_hard <- viterbi(y_bad, P, pi0, emission_hard)
vt_hard$feasible[1] FALSEvt_hard$log_prob[1] -Infvt_hard$path[1] NA NA NA NA NAvt_soft <- viterbi(y_bad, P, pi0, emission_noisy)
vt_soft$feasible[1] TRUEvt_soft$log_prob[1] -8.446402vt_soft$path[1] "SD" "SD" "CR" "CR" "CR"The noisy decoder does not walk through absorbing PD into CR. It treats the isolated PD as a likely misrecord of SD (or PR) and keeps a chain that can legally reach the later CRs. The opposite reading — trust PD and treat the CRs as false — is possible only because \(C_{\mathrm{PD},\mathrm{CR}} > 0\); here that path loses.
A small \(\varepsilon\) floor on a hard emission, \(e_t(s) = 1\{s=y_t\}(1-\varepsilon) + \varepsilon/(|S|-1)\), is the same idea with a less interpretable \(C\).
The decoder above is the teaching object. Production HMM fitting lives elsewhere:
| Tool | Role |
|---|---|
depmixS4 (Visser and Speekenbrink, 2010) |
discrete-time HMM / dependent mixtures, EM |
seqHMM (Helske and Helske, 2019) |
categorical sequences, multi-channel, Viterbi helpers |
msm (Jackson, 2011) |
continuous-time Markov and HMM on irregular visits |
| Zucchini, MacDonald, and Langrock (2016) | textbook discrete-time HMM in R |
None of these is a RECIST engine. They estimate and decode; the clinical kernel is still yours.
| Method | What it maximises / assumes | Typical failure |
|---|---|---|
| LOCF / last non-missing | persistence of the last seen label | PD then CR; ignores \(P\) |
| Independent imputation | visit-wise margins | incoherent chains |
| Forward–backward modes | \(\arg\max_s P(S_t=s \mid y)\) at each \(t\) | modes need not form a legal path |
| Multiple imputation | posterior of a scientific estimand | right for means / ORR; not a single chain |
| Viterbi | \(\arg\max_{s_{1:T}} P(s_{1:T} \mid y)\) | wrong \(P\) or wrong \(e_t\); not a BOR |
Use Viterbi when the scientific sentence is “the most probable course of statuses, given this kernel and these recordings.” Use RECIST machinery for the endpoint; use MMRM / MI / ICH E9(R1) for the treatment-effect estimand; use a confusion-matrix emission when the recording itself is on trial.