Discounted fundamental matrix of a time-homogeneous Markov chain at a constant discount rate

Signature

N_disc = (I - Q / (1 + r)^Delta)^(-1); t_disc = N_disc * c

Function

Discounted and episode-level occupancy of a time-homogeneous absorbing Markov chain

Computational function

• Computational function: expected discounted rewards of a time-homogeneous absorbing Markov chain by matrix inversion

  Returns the expected total reward per person of a time-homogeneous absorbing Markov chain from a starting distribution across the transient states and one reward per cycle for each transient state, using the discounted fundamental matrix of HE-FM-MKCH-001 in place of a cohort run. A reward of 1 in every transient state gives discounted life expectancy and utilities give quality-adjusted life expectancy. It applies Sonnenberg and Beck's half-cycle subtraction on request and checks that Q is substochastic and that the discounted series converges. The inputs differ from the formula's: a starting distribution, a reward vector and a half-cycle option.

  Inputs and outputs: Q: Transient block of the transition matrix, the same in every cycle; required. Unit: probability per cycle.; m0: Starting distribution across the transient states, summing to 1 for a cohort that starts alive; required. Unit: proportion.; u: Reward per cycle in each transient state, for example 1 for life years or a utility for QALYs; required. Unit: reward per cycle.; r: Annual discount rate, default 0. Unit: proportion per year.; Delta: Cycle length in years, default 1. Unit: years.; half_cycle: Whether to subtract half a cycle of the starting reward, default FALSE. Unit: logical.; N_disc: Discounted fundamental matrix. Unit: discounted cycles.; value: Expected discounted reward per person. Unit: discounted reward.; substochastic, converges: Checks that every entry of Q is non-negative with rows summing to at most 1, and that the largest absolute eigenvalue of the discounted Q is below 1. Unit: logical.

  Assumption: Q and the rewards are the same in every cycle, the discount rate is constant, membership is counted at the start of each cycle with the first cycle undiscounted, and the half-cycle subtraction removes half a cycle of the starting reward only, as Sonnenberg and Beck describe for a matrix solution. A model whose probabilities change with age needs the cohort trace of HE-CF-CSIM-001.

  Worked example (Life expectancy from Well, undiscounted): The article's chain gives 16 years, or 15.5 after the half-cycle subtraction. Q = [[0.85, 0.10], [0.20, 0.70]]; m0 = [1, 0]; u = [1, 1]; r = 0; Delta = 1; value = 16; value_half_cycle = 15.5

  Worked example (Quality-adjusted life expectancy from Well, undiscounted): With utilities of 0.85 for Well and 0.60 for Ill the value is 12.6, or 12.175 after subtracting half a cycle of the Well utility, as in the article. Q = [[0.85, 0.10], [0.20, 0.70]]; m0 = [1, 0]; u = [0.85, 0.60]; r = 0; Delta = 1; value = 12.6; value_half_cycle = 12.175

  Worked example (Discounted life expectancy from Well at 3.5%): The value is 10.7260 years, or 10.2260 with the half-cycle subtraction, as in the article. Q = [[0.85, 0.10], [0.20, 0.70]]; m0 = [1, 0]; u = [1, 1]; r = 0.035; Delta = 1; value = 10.7260; value_half_cycle = 10.2260

  Worked example (Discounted quality-adjusted life expectancy from Well at 3.5%): The value is 8.5007, or 8.0757 with the half-cycle subtraction (computed here for illustration). Q = [[0.85, 0.10], [0.20, 0.70]]; m0 = [1, 0]; u = [0.85, 0.60]; r = 0.035; Delta = 1; value = 8.5007; value_half_cycle = 8.0757

  Excel: With Q in TransientQ, the starting distribution as a row in StartRow, the rewards as a column in RewardCol, the rate in DiscRate, the cycle length in CycleYears and TRUE or FALSE in HalfCycle, Excel 365 returns the value: =INDEX(MMULT(MMULT(StartRow,MINVERSE(MUNIT(ROWS(TransientQ))-TransientQ/(1+DiscRate)^CycleYears)),RewardCol)-IF(HalfCycle,0.5*MMULT(StartRow,RewardCol),0),1,1)

  R: mc_value <- function(Q, m0, u, r = 0, Delta = 1, half_cycle = FALSE) { Q <- as.matrix(Q); d <- (1+r)^Delta; converges <- max(Mod(eigen(Q/d, only.values = TRUE)$values)) < 1; substochastic <- all(Q >= 0) && all(rowSums(Q) <= 1+1e-12); Nd <- solve(diag(nrow(Q))-Q/d); value <- drop(m0 %*% Nd %*% u); if (half_cycle) value <- value-0.5*sum(m0*u); list(N_disc = Nd, value = value, substochastic = substochastic, converges = converges) } Base R only; mc_value(rbind(c(0.85, 0.10), c(0.20, 0.70)), c(1, 0), c(0.85, 0.60), r = 0.035) returns 8.5007.

  Python: def mc_value(Q, m0, u, r=0.0, Delta=1.0, half_cycle=False): Q = np.asarray(Q, float); d = (1+r)**Delta; converges = bool(max(abs(np.linalg.eigvals(Q/d))) < 1); substochastic = bool((Q >= 0).all() and (Q.sum(axis=1) <= 1+1e-12).all()); Nd = np.linalg.inv(np.eye(len(Q))-Q/d); value = float(np.asarray(m0) @ Nd @ np.asarray(u)); value -= 0.5*float(np.dot(m0, u)) if half_cycle else 0.0; return {"N_disc": Nd, "value": value, "substochastic": substochastic, "converges": converges} Needs import numpy as np; returns the same values as the R function.

  Test (Matrix value matches a discounted cohort run): The discounted life expectancy from Well, 10.7260, equals a 3,000-cycle cohort run with HE-CF-CSIM-001 counted at the start of each cycle with discount factor 1.035 to the power minus k. Expected result: TRUE. Excel check: =ROUND(MatrixValue,4)=ROUND(TraceValue,4)

  Common error (Subtracting half a cycle from every state): Sonnenberg and Beck subtract half a cycle from the membership of the starting state only. Subtracting half a cycle of both the Well and the Ill utility from the undiscounted quality-adjusted value gives 11.875 in place of 12.175, because a person starting in Well spends no starting cycle in Ill.

  Source: Grinstead CM, Snell JL. Grinstead and Snell's Introduction to Probability. The CHANCE Project version of 4 July 2006, based on the 2nd edition published by the American Mathematical Society (full text read). Section 11.2, Theorems 11.4 and 11.5 (N = I + Q + Q^2 + ... and t = Nc); Sonnenberg FA, Beck JR. Markov models in medical decision making: a practical guide. Medical Decision Making. 1993;13(4):322-338. doi:10.1177/0272989X9301300409 (full text read). The Half-Cycle Correction (subtraction of one half cycle from the membership of each starting state) and Discounting.

  N_disc = (I - Q / (1 + r)^Delta)^(-1); value = m0 * N_disc * u - half_cycle * 0.5 * m0 * u

Try this function

Implementations

• Excel

  Discounted fundamental matrix and its row sums in Excel

  With the transient block in the named range TransientQ, the annual discount rate in DiscRate and the cycle length in years in CycleYears, Excel 365 returns N_disc as a dynamic array with the first formula and the column t_disc with the second. MUNIT builds the identity matrix and SEQUENCE the column of ones.

  =MINVERSE(MUNIT(ROWS(TransientQ))-TransientQ/(1+DiscRate)^CycleYears); =MMULT(MINVERSE(MUNIT(ROWS(TransientQ))-TransientQ/(1+DiscRate)^CycleYears),SEQUENCE(ROWS(TransientQ),1,1,0))

Assumptions

• Constant transition matrix and constant discount rate for the discounted matrix

  Q is the same in every cycle and the discount rate does not change over the horizon. Sonnenberg and Beck wrote that discounting cannot be used with the fundamental matrix solution because of its time variance; a constant geometric discount multiplies every cycle by the same factor and so folds into Q, but a discount schedule that changes over time, or a matrix that changes with age, does not, and the discounted total then comes from a cohort trace (HE-FM-CSIM-001).

• Start-of-cycle counting with an undiscounted first cycle in the discounted fundamental matrix

  Membership at the start of cycle k, for k from 0 upwards, is counted in full and discounted by (1 + r) to the power minus k Delta, so the starting cycle has a factor of 1. No half-cycle correction is applied. Sonnenberg and Beck's correction for a matrix solution subtracts half a cycle from the membership of each starting state, which at a factor of 1 subtracts 0.5 from t_disc.

• Absorbing chain for the discounted fundamental matrix

  Every transient state can reach an absorbing state, so Q to the power k tends to zero, as Grinstead and Snell show for any absorbing chain, and the discounted powers of Q tend to zero faster still; I minus Q divided by the discount factor can then be inverted.

Worked examples

• Well and Ill chain with recovery discounted at 3.5 per cent a year

  In the article's chain, 85% of Well stay Well, 10% fall ill and 5% die each year, and of Ill 20% recover, 70% stay Ill and 10% die. At 3.5% a year with annual cycles, N_disc has rows 8.2603, 2.4658 and 4.9315, 4.5616, so discounted life expectancy is 10.7260 years from Well, matching a discounted cohort run, or 10.2260 after subtracting half a cycle, as in the article. From Ill it is 9.4932 years (computed here for illustration).

  Q = [[0.85,0.10],[0.20,0.70]]; I = [[1,0],[0,1]]; r = 0.035; Delta = 1; c = [1,1]; N_disc = [[8.2603,2.4658],[4.9315,4.5616]]; t_disc = [10.7260,9.4932]

• Zero discount rate returns the undiscounted Well and Ill fundamental matrix

  With r equal to zero the discount factor is 1 and N_disc is the fundamental matrix of the article, rows 12, 4 and 8, 6, with 16 expected years before death from Well and 14 from Ill, the result of HE-FM-ABS-001.

  Q = [[0.85,0.10],[0.20,0.70]]; I = [[1,0],[0,1]]; r = 0; Delta = 1; c = [1,1]; N_disc = [[12,4],[8,6]]; t_disc = [16,14]

Common errors

• Discounting undiscounted Markov chain life expectancy as one lump

  Discounting the 16 undiscounted years from Well as if they all fell at year 16 gives 9.2273, and treating them as 16 certain years from the start gives 12.5174. Both miss the 10.7260 of N_disc, because survival time is spread over many cycles and the discount factor is not linear in time.

• Annual discount factor applied to shorter Markov chain cycles

  With six-month cycles the per-cycle factor is 1.035 to the power 0.5, about 1.0173. Dividing Q by 1.035 discounts every half-year as if it were a full year and understates discounted life expectancy.

• Mixing counting conventions with the discounted fundamental matrix

  Counting membership from the end of the first cycle, with factor 1.035 to the power minus 1 for the first cycle counted, drops the undiscounted starting cycle and gives 9.7260 from Well, exactly one cycle less than t_disc. That total should not then also have half a cycle subtracted, which would count the same timing adjustment twice.

Sources

• Grinstead and Snell on the series form of the fundamental matrix

  Grinstead CM, Snell JL. Grinstead and Snell's Introduction to Probability. The CHANCE Project version of 4 July 2006, based on the 2nd edition published by the American Mathematical Society (full text read). Chapter 11 Markov Chains, section 11.2 Absorbing Markov Chains: Theorem 11.3 (Q to the power n tends to 0), Theorem 11.4 (I minus Q has an inverse N and N = I + Q + Q^2 + ..., proved from Q^n tending to 0) and Theorem 11.5 (t = Nc). Section 11.5 notes that the series converges because Q^n tends to 0, so it acts like a convergent geometric series. The discounted matrix follows by replacing Q with Q divided by the discount factor, whose powers also tend to 0; this step is the Markov Chain article's derivation and is not stated by Grinstead and Snell.

  View source →

• Sonnenberg and Beck on counting, half-cycle correction and discounting with the matrix

  Sonnenberg FA, Beck JR. Markov models in medical decision making: a practical guide. Medical Decision Making. 1993;13(4):322-338. doi:10.1177/0272989X9301300409 (full text read). The Half-Cycle Correction: the fundamental matrix representation is equivalent to counting state membership at the beginning of each cycle, and the correction for a matrix solution is subtraction of one half cycle from the membership of each starting state. Discounting: because of the time variance, discounting cannot be used when the fundamental matrix solution is used; the Markov Chain article qualifies this for a constant geometric discount.

  View source →

Canonical Identity