Exact interval probabilities in a progressive illness-death model with constant rates

Signature

p_WW = exp(-(q_WS + q_WD) * t); p_WS = q_WS / (q_SD - q_WS - q_WD) * (exp(-(q_WS + q_WD) * t) - exp(-q_SD * t)); p_WD = 1 - p_WW - p_WS

Function

Continuous-time transition probabilities and expected state occupancy from transition intensities

Computational function

• Computational function: transition probability matrix from an intensity matrix by the matrix exponential

  Computes the transition probability matrix over any interval from an intensity matrix of any size, P(t) = exp(tQ), by scaling and squaring a truncated Taylor series, for models with recovery or more states where no closed form like HE-FM-CONT-001 is at hand. The inputs differ from the formula's: a whole intensity matrix instead of three named rates.

  Inputs and outputs: Q: Square intensity matrix with the rate from state i to state j off the diagonal and minus the row's other entries on the diagonal, so each row adds to 0; required. Unit: events per person-year.; t: Interval length; required, above 0. Unit: years.; m: Number of squarings, default 10. Unit: count.; terms: Number of Taylor terms, default 12. Unit: count.; P: Transition probability matrix over t, rows adding to 1. Unit: probability.

  Assumption: Constant intensities over t (time-homogeneous Markov process); the series of t Q divided by 2 to the power m converges quickly when that scaled matrix is small, which the defaults ensure for annual rates of the size used in health economic models.

  Worked example (Illness-death model, article example): With rates of 0.30 Well to Sick, 0.05 Well to Dead and 0.40 Sick to Dead, the one-year matrix has first row 0.704688, 0.206208 and 0.089104 and second row 0, 0.670320 and 0.329680, matching the closed form; over a month the first row is 0.971255, 0.024231 and 0.004515. Q = [[-0.35, 0.30, 0.05], [0, -0.40, 0.40], [0, 0, 0]]; t = 1; P = [[0.704688, 0.206208, 0.089104], [0, 0.670320, 0.329680], [0, 0, 1]]

  Worked example (Recovery from Sick at 0.10 a year): Adding recovery makes the model non-progressive; the one-year first row becomes 0.714774, 0.197298 and 0.087928 and the second row 0.065766, 0.616125 and 0.318109 (computed here for illustration). Q = [[-0.35, 0.30, 0.05], [0.10, -0.50, 0.40], [0, 0, 0]]; t = 1; P = [[0.714774, 0.197298, 0.087928], [0.065766, 0.616125, 0.318109], [0, 0, 1]]

  Worked example (Liver disease cohort over 34 months, Chhatwal and colleagues): With the appendix intensity matrix and 187 patients starting in decompensated cirrhosis, exp(34/12 Q) gives about 72 still in cirrhosis, 15 with carcinoma and 100 dead, the 72 matching the original study. Q = [[-0.3369, 0.0967, 0.2402], [0, -0.5573, 0.5573], [0, 0, 0]]; t = 2.833333; counts = [72.0, 14.7, 100.3]

  Excel: Excel has no matrix exponential. With Q in a three-by-three range IntensityQ and the step h = t/1024 in StepH, =MUNIT(3)+StepH*IntensityQ+MMULT(StepH*IntensityQ,StepH*IntensityQ)/2 gives a one-step matrix A; applying =MMULT(A,A) ten times in turn gives P(t) to within about 1E-7 for the examples here.

  R: expm_ss <- function(Q, t, m = 10, terms = 12) { A <- Q*t/2^m; P <- diag(nrow(Q)); Tk <- diag(nrow(Q)); for (k in 1:terms) { Tk <- Tk %*% A/k; P <- P+Tk }; for (j in 1:m) P <- P %*% P; P } Base R only; expm_ss(rbind(c(-0.35, 0.30, 0.05), c(0, -0.40, 0.40), c(0, 0, 0)), 1) returns the first example's matrix. The Matrix package's expm() gives the same result.

  Python: def expm_ss(Q, t, m=10, terms=12): n = len(Q); mul = lambda A, B: [[sum(A[i][k]*B[k][j] for k in range(n)) for j in range(n)] for i in range(n)]; A = [[q*t/2**m for q in row] for row in Q]; I = [[float(i == j) for j in range(n)] for i in range(n)]; T = list(itertools.accumulate(range(1, terms+1), lambda T, k: [[v/k for v in row] for row in mul(T, A)], initial=I)); P = [[sum(Tk[i][j] for Tk in T) for j in range(n)] for i in range(n)]; return functools.reduce(lambda P, _: mul(P, P), range(m), P) Needs import itertools, functools; returns the same matrices as the R function, with Q as a list of rows. scipy.linalg.expm(t*Q) gives the same result.

  Test (Rows add to one): Every row of P adds to 1 because every row of Q adds to 0. Expected result: TRUE. Excel check, with P in ProbMatrix: =MAX(ABS(MMULT(ProbMatrix,{1;1;1})-1))<1E-9

  Test (Two half intervals equal one interval): P(t / 2) multiplied by itself equals P(t). Expected result: TRUE. Excel check, with the half-interval matrix in HalfMatrix built by the same steps with t / 2: =MAX(ABS(MMULT(HalfMatrix,HalfMatrix)-ProbMatrix))<1E-6

  Common error (Exponentiating each entry of Q): EXP applied cell by cell to tQ is not the matrix exponential; for the article's Q it puts exp(0.30), about 1.35, in the Well to Sick cell, which is not a probability.

  Source: Jackson C. Multi-state modelling with R: the msm package. Version 1.8.2. Cambridge: MRC Biostatistics Unit; 7 November 2024. Section 1.4, equation 2 (P(t) = Exp(tQ), defined by the power series with matrix products, difficult to calculate reliably in general); Chhatwal J, Jayasuriya S, Elbasha EH. Changing cycle lengths in state-transition models: challenges and solutions. Medical Decision Making. 2016;36(8):952-964. doi:10.1177/0272989X16656165. Theoretical issues and Appendix A (intensity matrix, annual and monthly matrices and the 34-month check against 72 patients).

  A = Q * t / 2^m; P = (sum_(k=0)^terms A^k / k!)^(2^m)

Try this function

Implementations

• Excel

  Illness-death interval probabilities from named rates

  With RateWS, RateWD, RateSD and CycleLength named, the three formulas return the probabilities of staying Well, being Sick and being dead, held in ProbWW, ProbWS and ProbWD.

  =EXP(-(RateWS+RateWD)*CycleLength); =RateWS/(RateSD-RateWS-RateWD)*(EXP(-(RateWS+RateWD)*CycleLength)-EXP(-RateSD*CycleLength)); =1-ProbWW-ProbWS

Assumptions

• Constant transition rates over the illness-death interval

  The three intensities do not change within the interval (a time-homogeneous Markov process); with age-dependent rates the formula applies piece by piece over intervals short enough for the rates to be treated as constant.

• Progressive illness-death structure without recovery

  Nobody returns from Sick to Well, and the Sick exit rate differs from the Well exit rate; when they are equal the Sick probability is q_WS times t times exp(minus q_SD t). Models with recovery or more states need the matrix exponential (HE-CF-CONT-001).

Worked examples

• One-year probabilities in the article's illness-death example

  With rates of 0.30, 0.05 and 0.40 a year, staying Well is exp(minus 0.35), about 0.7047; being Sick is 6 times (0.704688 minus 0.670320), about 0.2062; and death by either route takes the remaining 0.0891, as in the article.

  q_WS = 0.3; q_WD = 0.05; q_SD = 0.4; t = 1; p_WW = 0.7047; p_WS = 0.2062; p_WD = 0.0891

• One-month probabilities in the article's illness-death example

  Over a month the same rates give about 0.9713, 0.0242 and 0.0045 (computed here for illustration).

  q_WS = 0.3; q_WD = 0.05; q_SD = 0.4; t = 0.083333; p_WW = 0.9713; p_WS = 0.0242; p_WD = 0.0045

• Decompensated cirrhosis, liver cancer and death

  Chhatwal and colleagues estimated rates of 0.0967 from decompensated cirrhosis to hepatocellular carcinoma, 0.2402 from cirrhosis to death and 0.5573 from carcinoma to death; the annual probabilities are about 0.7140, 0.0620 and 0.2241 (the source reports 0.0619 for the middle value).

  q_WS = 0.0967; q_WD = 0.2402; q_SD = 0.5573; t = 1; p_WW = 0.714; p_WS = 0.062; p_WD = 0.2241

Common errors

• Converting each rate separately into a probability

  Treating each rate as the only risk gives 1 minus exp(minus 0.30), about 0.2592, for Sick and 0.0488 for death, and a lifetime model with half-cycle counting that overstates life expectancy by 0.300 years, about 6 per cent, in the article's example; this traditional conversion is applicable only to two-state models.

• Using the competing-risk split for an annual illness-death matrix

  The split gets staying Well right but counts only the 0.0422 direct deaths, missing the 0.0469 who fall Sick and die within the year; the annual model with half-cycle counting then gives 5.486 life-years against the true 5.000.

• Entering a rate as a probability

  Using the Well to Sick rate of 0.30 as an annual probability overstates it against the single-risk 0.2592 and the exact 0.2062; conversions between time units go through rates, and probabilities should never be called rates.

Sources

• Matrix exponential and its closed form for the three-state illness-death model

  Jackson C. Multi-state modelling with R: the msm package. Version 1.8.2. Cambridge: MRC Biostatistics Unit; 7 November 2024. Section 1.4, equation 2: for a time-homogeneous process P(t) is the matrix exponential of t times the transition intensity matrix Q; for the three-state illness-death model with no recovery, p11(t) = exp(minus (q12 plus q13) t), p12(t) = q12 / (q12 plus q13 minus q23) times (exp(minus q23 t) minus exp(minus (q12 plus q13) t)) when q12 plus q13 differs from q23, q12 t exp(minus (q12 plus q13) t) when they are equal, and p13(t) the remainder.

  View source →

• Separate probability conversion fails with competing risks; worked liver disease intensity matrix

  Chhatwal J, Jayasuriya S, Elbasha EH. Changing cycle lengths in state-transition models: challenges and solutions. Medical Decision Making. 2016;36(8):952-964. doi:10.1177/0272989X16656165. Background (Transition probabilities): the traditional approach of converting each probability through its rate is applicable only to two-state models and can result in significant error; Theoretical issues: P(t) = exp(tQ) for an intensity matrix Q; Appendix A: rates of 0.0967 (decompensated cirrhosis to hepatocellular carcinoma), 0.2402 (death from cirrhosis) and 0.5573 (death from carcinoma) give the annual matrix with first row 0.7140, 0.0619, 0.2241 by the matrix exponential.

  View source →

• Probabilities converted through rates and never called rates

  Siebert U, Alagoz O, Bayoumi AM, Jahn B, Owens DK, Cohen DJ, Kuntz KM. State-transition modeling: a report of the ISPOR-SMDM Modeling Good Research Practices Task Force-3. Value in Health. 2012;15(6):812-820. doi:10.1016/j.jval.2012.06.014. Parameter derivation: the conversion of transition probabilities from one time unit to another should be done through rates, and to avoid confusion probabilities should never be called rates; cycle length should be short enough that an event occurs at most once per cycle.

  View source →

Canonical Identity