Expected outcome per person from a cohort simulation trace with counting weights and discounting

Signature

Y = sum_(t=0)^T [w_t * delta_t * m_t * r]

• T Number of model cycles in the time horizon, so the trace has rows t = 0 to T (cycles)

Function

Cohort simulation of a state-transition model

Computational function

• Computational function: cohort trace run with rewards accumulated under a counting rule

  Runs a cohort simulation from the inputs a model holds, a starting distribution, a transition matrix, the reward per cycle in each state, a number of cycles, a counting rule and a discount rate per cycle, and returns the cohort trace and the expected total outcome per person. It first builds the trace row by row with the cohort update HE-FM-MM-001, then forms the counting weights and discount factors and applies HE-FM-CSIM-001. The inputs therefore differ from the formula's variables: the formula takes the trace rows, weights and discount factors as given, while the function generates them.

  Inputs and outputs: m_0: Starting distribution of the cohort across states, a row that sums to one; required. Unit: proportion of the cohort.; P: Transition matrix with rows that sum to one, the same in every cycle; required. Unit: probability per cycle.; r: Reward per cycle in each state; required. Unit: outcome per person per cycle.; T: Number of cycles; required, a positive whole number, and even for Simpson's rule. Unit: cycles.; rule: Counting rule, one of end, start, half or simpson; optional, default half. Unit: none.; d: Discount rate per cycle; optional, default 0. Unit: rate per cycle.; M: Cohort trace, one row per cycle from t = 0 to T. Unit: proportion of the cohort.; Y: Expected total outcome per person. Unit: the unit of r summed over the horizon.

  Assumption: The transition matrix is constant across cycles and the discount factor for cycle t is 1/(1+d)^t from t = 0. For probabilities that change with age or time, the R version takes a list of matrices by replacing P with P_list[[k]] inside Reduce. Each state is homogeneous, so a subgroup with different risks needs its own run.

  Worked example (Article QALYs over four cycles with half-cycle weights): The article's three-state model with utilities 0.9 and 0.6 returns 2.8445 QALYs per person, as in the article. m_0 = [1,0,0]; P = [[0.80,0.15,0.05],[0,0.70,0.30],[0,0,1]]; r = [0.9,0.6,0]; T = 4; d = 0; rule = half; Y = 2.8445

  Worked example (Article life years over four cycles): A reward of 1 in each alive state returns 3.4124 life years per person, the article's figure. m_0 = [1,0,0]; P = [[0.80,0.15,0.05],[0,0.70,0.30],[0,0,1]]; r = [1,1,0]; T = 4; d = 0; rule = half; Y = 3.4124

  Worked example (Run to absorption counted at the start of each cycle): Run for 400 cycles, by which time the cohort is absorbed, start-of-cycle counting returns 7.5 cycles alive, the 5 cycles in Healthy and 2.5 in Sick that the fundamental matrix HE-FM-ABS-001 gives. m_0 = [1,0,0]; P = [[0.80,0.15,0.05],[0,0.70,0.30],[0,0,1]]; r = [1,1,0]; T = 400; d = 0; rule = start; Y = 7.5

  Worked example (Half-cycle QALYs from a run to absorption): With half-cycle weights the same run returns the article's half-cycle-corrected 5.55 QALYs, the matrix result minus half a cycle in the starting state. m_0 = [1,0,0]; P = [[0.80,0.15,0.05],[0,0.70,0.30],[0,0,1]]; r = [0.9,0.6,0]; T = 400; d = 0; rule = half; Y = 5.55

  Excel: =MMULT(PreviousRow,TransitionMatrix) Copied down one row per cycle from the starting row, this builds the trace; =SUMPRODUCT(CycleWeights,DiscountFactors,MMULT(Trace,StateRewards)) then returns Y, as under HE-FM-CSIM-001.

  R: cohort_run <- function(m0, P, r, n_cycles, rule="half", disc=0) { if (rule == "simpson" && n_cycles %% 2 != 0) stop("Simpson rule needs an even number of cycles"); M <- do.call(rbind, Reduce(function(x, k) x %*% P, seq_len(n_cycles), matrix(m0, nrow=1), accumulate=TRUE)); w <- switch(rule, end=c(0, rep(1, n_cycles)), start=c(rep(1, n_cycles), 0), half=c(0.5, rep(1, n_cycles-1), 0.5), simpson=c(1, ifelse(seq_len(n_cycles-1) %% 2 == 1, 4, 2), 1)/3); delta <- 1/(1+disc)^(0:n_cycles); list(trace=M, Y=sum(w*delta*(M %*% r))) } Returns a list; the total is read with [["Y"]] and the trace with [["trace"]].

  Python: def cohort_run(m0, P, r, n_cycles, rule="half", disc=0.0): assert not (rule == "simpson" and n_cycles % 2), "Simpson rule needs an even number of cycles"; P = np.asarray(P, float); M = np.vstack([np.asarray(m0, float) @ np.linalg.matrix_power(P, t) for t in range(n_cycles+1)]); w = {"end": [0]+[1]*n_cycles, "start": [1]*n_cycles+[0], "half": [0.5]+[1]*(n_cycles-1)+[0.5], "simpson": [1/3]+[4/3 if t % 2 else 2/3 for t in range(1, n_cycles)]+[1/3]}[rule]; delta = (1+disc)**-np.arange(n_cycles+1.0); return M, float(np.sum(np.asarray(w)*delta*(M @ np.asarray(r, float)))) Requires numpy imported as np; returns the trace and the total.

  Test (Every trace row sums to one): Each row of the trace returned by the function adds to 1, because the states are exhaustive and every row of P sums to 1. Expected result: TRUE. Excel check: =ABS(SUM(TraceRow)-1)<1E-9

  Test (Long start-of-cycle run matches the fundamental matrix): A 400-cycle run counted at the start of each cycle with a reward of 1 in each alive state returns the total of HE-FM-ABS-001, 7.5 cycles in the article's example. Expected result: TRUE. Excel check: =ROUND(SUMPRODUCT(StartWeights,MMULT(LongTrace,AliveRewards)),4)=7.5

  Common error (One run at averaged transition probabilities): Outcomes are not linear in transition probabilities. If half a cohort survives each cycle with probability 0.9 and half with 0.5, the expected cycles alive counted at the start of each cycle are 10 and 2, a mean of 6, whereas one run at the average of 0.7 gives about 3.33. Subgroups that differ in risk need separate runs, with results weighted afterwards.

  Source: Sonnenberg FA, Beck JR. Markov models in medical decision making: a practical guide. Medical Decision Making. 1993;13(4):322-338. Sections on the Markov cohort simulation, the fundamental matrix solution and the half-cycle correction, which counts state membership at the start of each cycle in the matrix solution. Alarid-Escudero F, Krijkamp E, Enns EA, Yang A, Hunink MGM, Pechlivanoglou P, Jalal H. An introductory tutorial on cohort state-transition models in R using a cost-effectiveness analysis example. Medical Decision Making. 2023;43(1):3-20, equations 4 to 6.

  m_t = m_0 * P^t; Y = sum_(t=0)^T [w_t * delta_t * m_t * r]

Try this function

Implementations

• Excel

  Expected outcome from a trace with weights and discounting

  With the trace in a range named Trace (one row per cycle from t = 0, one column per state), the rewards in a column named StateRewards and the counting weights and discount factors in columns CycleWeights and DiscountFactors with one row per cycle, the cell returns Y.

  =SUMPRODUCT(CycleWeights,DiscountFactors,MMULT(Trace,StateRewards))

Assumptions

• Mutually exclusive, exhaustive and homogeneous cohort states

  Every member of the cohort is in exactly one state in every cycle, so each trace row sums to one, and everyone in a state faces the same transition probabilities and rewards. Under these conditions Y is an expected value per person and the cohort size is bookkeeping only.

• Rewards per cycle include the cycle length

  r is the reward for spending one whole cycle in a state. For QALYs it is the utility weight multiplied by the cycle length, so a six-month cycle in a state with utility 0.8 earns 0.4 QALYs.

• Counting rule stated with the cycle length

  The trace records the cohort at cycle boundaries while people move at any time within a cycle. Full credit at the start of a cycle overestimates expected values and no credit underestimates them, so the counting rule, the cycle length and any correction are reported and can be varied in sensitivity analysis.

Worked examples

• Expected years in Healthy counted at the end of each cycle

  In the article's illustrative three-state model, the Healthy column of the trace over four annual cycles is 1, 0.8, 0.64, 0.512 and 0.4096. With a reward of 1 in Healthy and no discounting, end-of-cycle weights give 2.3616 years.

  w_t = [0, 1, 1, 1, 1]; delta_t = [1, 1, 1, 1, 1]; m_t = [1, 0.8, 0.64, 0.512, 0.4096]; r = 1; T = 4; Y = 2.3616

• Expected years in Healthy counted at the start of each cycle

  Start-of-cycle weights credit the full cohort at t = 0 and nothing at t = 4, giving 2.952 years, as in the article.

  w_t = [1, 1, 1, 1, 0]; delta_t = [1, 1, 1, 1, 1]; m_t = [1, 0.8, 0.64, 0.512, 0.4096]; r = 1; T = 4; Y = 2.952

• Expected years in Healthy with half-cycle weights

  Half weights at t = 0 and t = 4 give 2.6568 years, the mean of the start and end counts and the article's figure.

  w_t = [0.5, 1, 1, 1, 0.5]; delta_t = [1, 1, 1, 1, 1]; m_t = [1, 0.8, 0.64, 0.512, 0.4096]; r = 1; T = 4; Y = 2.6568

• Expected years in Healthy with Simpson's 1/3 weights

  Simpson's weights of one third, four thirds, two thirds, four thirds and one third, entered to nine decimals, give 2.6459 years. The continuous-time area from HE-FM-BTH-001 is about 2.6458, computed here with unrounded ln 0.8; the article prints 2.6459. The two agree to within 0.0001.

  w_t = [0.333333333, 1.333333333, 0.666666667, 1.333333333, 0.333333333]; delta_t = [1, 1, 1, 1, 1]; m_t = [1, 0.8, 0.64, 0.512, 0.4096]; r = 1; T = 4; Y = 2.6459

• Discounted years in Healthy with half-cycle weights

  Computed here for illustration: discounting the half-cycle count at an illustrative 3.5% a year, with discount factors of 1 over 1.035 to the power t rounded to six decimals, gives 2.5107 discounted years in Healthy instead of 2.6568.

  w_t = [0.5, 1, 1, 1, 0.5]; delta_t = [1, 0.966184, 0.933511, 0.901943, 0.871442]; m_t = [1, 0.8, 0.64, 0.512, 0.4096]; r = 1; T = 4; Y = 2.5107

Common errors

• Counting rule left unstated in a cohort simulation

  End-of-cycle counting gives 2.3616 expected years in Healthy in the article's example, about 10.7% below the continuous-time area, and start-of-cycle counting gives 2.952, about 11.6% above it. The errors need not cancel in incremental results, so the rule has to be stated.

• Half-cycle weight dropped at the end of a truncated horizon

  Applying the half weight only at t = 0 when the horizon stops at four cycles gives 2.8616 years in Healthy instead of 2.6568, computed here for illustration. Over a horizon that is not lifetime, the final cycle also takes a half weight.

• Simpson weights with an odd number of cycles

  Simpson's 1/3 rule needs an even T. With three cycles the pattern of one third, four thirds, two thirds and one third sums to two and two thirds instead of 3, computed here for illustration, so the count is too low.

Sources

• Trace, reward vector, discount and within-cycle correction vectors

  Alarid-Escudero F, Krijkamp E, Enns EA, Yang A, Hunink MGM, Pechlivanoglou P, Jalal H. An introductory tutorial on cohort state-transition models in R using a cost-effectiveness analysis example. Medical Decision Making. 2023;43(1):3-20. Equations 4 to 6, which multiply the cohort trace by a vector of state rewards and weight the result by discount factors 1/(1+d)^t from t = 0 and a within-cycle correction vector, and the section on within-cycle correction giving the Simpson's 1/3 weights (numbering as in the arXiv preprint 2001.07824v4).

  View source →

• Within-cycle corrections as numerical integration methods

  Elbasha EH, Chhatwal J. Theoretical foundations and practical applications of within-cycle correction methods. Medical Decision Making. 2016;36(1):115-131. Abstract: seven methods from numerical integration, including trapezoids and Simpson's 1/3 rule; the standard half-cycle correction gives the same results as the trapezoidal rule, and errors need not cancel in incremental outcomes.

  View source →

• Half-cycle correction in the ISPOR-SMDM state-transition report

  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. Section on the half-cycle correction and best practice III-14, applying it to costs and effectiveness in the first cycle and in the final cycle if a lifetime horizon is not used.

  View source →

Canonical Identity