Two-parameter Cholesky step for a correlated pair of normal PSA draws

Signature

L_21 = rho * s_2; L_22 = s_2 * sqrt(1 - rho^2); x_1 = mu_1 + s_1 * z_1; x_2 = mu_2 + L_21 * z_1 + L_22 * z_2

Function

Probabilistic sensitivity analysis run that turns joint parameter draws into stored decision outputs

Computational function

• Computational function for a PSA run with correlated log baseline rate and log hazard ratio through a constant-rate survival model

  Runs a complete probabilistic sensitivity analysis on a small fixed set of draws. For each draw it takes a pair of independent standard normal values, turns them into a correlated log baseline mortality rate and log hazard ratio with the Cholesky step HE-FM-PSA-001, runs a constant-rate survival model for both options with those same parameters, and stores the new option's incremental net health benefit. It then returns the expected incremental net health benefit (HE-FM-ENB-002), the probability that the new option is cost-effective (HE-FM-CEAC-001, with ties counted as not cost-effective), the two expected losses of HE-FM-PSA-002, the Monte Carlo standard error of the mean (the form of HE-FM-ANTV-002 with one observation per draw) and the number of draws for a target standard error (HE-FM-ABM-003). As in the article, mean survival under a constant mortality rate m is 1 / m years and each life-year is worth q QALYs without discounting, so a draw with log rate b0 and log hazard ratio b1 gives an incremental QALY gain of q times exp(minus b0 minus b1) less exp(minus b0). The Math / Code line writes out five draws; the R and Python versions draw any number.

  Inputs and outputs: mu_b0 and s_b0 are the mean and standard deviation of the log annual mortality rate on the comparator, mu_b1 and s_b1 those of the log hazard ratio, rho their correlation, above minus 1 and below 1; z1_1 to z1_5 and z2_1 to z2_5 are the five pairs of standard normal values; q is the QALYs per life-year; dC is the fixed incremental cost in currency; lambda is the threshold in currency per QALY; e_target is the target Monte Carlo standard error in QALYs. The outputs are the drawn b0_1 to b0_5 and b1_1 to b1_5, the stored INHB_1 to INHB_5 in QALYs per person, ENHB, P_CE as a proportion, EL_new and EL_comp in QALYs per person, the sample standard deviation SD_INHB, its Monte Carlo standard error MCSE and the number of draws N_req before rounding up.

  Assumption: The two log parameters are jointly normal, each draw uses a fresh pair of standard normal values, both options in a draw use the same drawn parameters, mortality is constant over time, the incremental cost is fixed and the draws are equally likely. The illustrative inputs below are not estimates for any real technology.

  Worked example (Article's five draws reproduced as a limiting case): With s_b0 equal to zero, rho equal to zero, s_b1 equal to 1 and the second standard normal values set to the logs of the article's rate ratios 0.4, 0.8, 0.8, 1 and 1, the five draws are the article's mortality rates 0.10, 0.20, 0.20, 0.25 and 0.25 against 0.25 on usual care. The function returns the article's incremental net health benefits, an expected value of 0.2 QALYs, a probability of 0.2 and expected losses of 0.28 and 0.48. The Monte Carlo standard error of about 0.561 QALYs is larger than the mean itself, and about 15,750 draws would be needed for a standard error of 0.01 if the spread of these five draws held. mu_b0 = -1.3862944; s_b0 = 0; mu_b1 = 0; s_b1 = 1; rho = 0; z1_1 = 0; z1_2 = 0; z1_3 = 0; z1_4 = 0; z1_5 = 0; z2_1 = -0.9162907; z2_2 = -0.2231436; z2_3 = -0.2231436; z2_4 = 0; z2_5 = 0; q = 0.5; dC = 15000; lambda = 25000; e_target = 0.01; INHB_1 = 2.4; INHB_2 = -0.1; INHB_3 = -0.1; INHB_4 = -0.6; INHB_5 = -0.6; ENHB = 0.2; P_CE = 0.2; EL_new = 0.28; EL_comp = 0.48; SD_INHB = 1.25499; MCSE = 0.56125; N_req = 15750

  Worked example (Five illustrative correlated draws): With the illustrative log rate and log hazard ratio of HE-EX-PSA-001 (means minus 1.3862944 and minus 0.2231436, standard deviations 0.10 and 0.20, correlation minus 0.5), an incremental cost of GBP 15,000, a threshold of GBP 25,000 per QALY and five illustrative pairs of standard normal values, the expected incremental net health benefit is about minus 0.1009 QALYs, two draws in five favour the new option and the expected losses are about 0.2426 from adopting and 0.1417 from keeping the comparator. Five draws give a Monte Carlo standard error of about 0.205, so the result says nothing reliable about the sign. mu_b0 = -1.3862944; s_b0 = 0.1; mu_b1 = -0.2231436; s_b1 = 0.2; rho = -0.5; z1_1 = 0.5; z1_2 = -1.0; z1_3 = 1.3; z1_4 = -0.2; z1_5 = 0.8; z2_1 = -1.2; z2_2 = 0.4; z2_3 = 0.9; z2_4 = -0.6; z2_5 = 1.5; q = 0.5; dC = 15000; lambda = 25000; e_target = 0.01; INHB_1 = 0.5751; INHB_2 = -0.47768; INHB_3 = -0.21705; INHB_4 = 0.13339; INHB_5 = -0.51823; ENHB = -0.10090; P_CE = 0.4; EL_new = 0.24259; EL_comp = 0.14170; SD_INHB = 0.45866; MCSE = 0.20512; N_req = 2103.72

  Excel: With frozen standard normal pairs in columns A and B (each =NORM.S.INV(RAND()) pasted as values) and the named inputs of HE-IM-PSA-001, =MeanLogRate+SDLogRate*A2 and =MeanLogHR+CholOffDiag*A2+CholDiagSecond*B2 filled down give columns C and D, named DrawLogRate and DrawLogHR, and =QALYPerLifeYear*(EXP(-C2-D2)-EXP(-C2))-IncCost/Threshold filled down in column E gives the column named DrawIncNHB. =AVERAGE(DrawIncNHB), =COUNTIF(DrawIncNHB,">0")/COUNT(DrawIncNHB) and =STDEV.S(DrawIncNHB)/SQRT(COUNT(DrawIncNHB)) return the expected value, the probability and the Monte Carlo standard error, and the formulas of HE-IM-PSA-002 applied to DrawIncNHB return the two expected losses.

  R: psa_run <- function(N, mu, s, rho, q, dC, lambda, seed) { set.seed(seed); S <- matrix(c(s[1]^2, rho*s[1]*s[2], rho*s[1]*s[2], s[2]^2), 2); L <- t(chol(S)); Z <- matrix(rnorm(2*N), 2); B <- mu + L %*% Z; INHB <- q*(exp(-B[1, ]-B[2, ])-exp(-B[1, ]))-dC/lambda; c(ENHB = mean(INHB), P_CE = mean(INHB > 0), EL_new = mean(pmax(-INHB, 0)), EL_comp = mean(pmax(INHB, 0)), MCSE = sd(INHB)/sqrt(N), r = cor(B[1, ], B[2, ])) } Run as psa_run(10000, c(log(0.25), log(0.8)), c(0.1, 0.2), -0.5, 0.5, 15000, 25000, 2026), which returns an expected incremental net health benefit of about minus 0.0761 QALYs, a probability of 0.4136, expected losses of about 0.2342 and 0.1581, a Monte Carlo standard error of about 0.0049 and a sample correlation of about minus 0.5005. The exact expected value for these inputs, from the log-normal mean of HE-FM-EV-004 (derived here), is about minus 0.0722 QALYs, about 0.79 Monte Carlo standard errors from the R result.

  Python: def psa_run(N, mu, s, rho, q, dC, lam, seed): import random, math; g = random.Random(seed); L21 = rho*s[1]; L22 = s[1]*math.sqrt(1-rho**2); Z = [(g.gauss(0, 1), g.gauss(0, 1)) for _ in range(N)]; B = [(mu[0]+s[0]*a, mu[1]+L21*a+L22*b) for a, b in Z]; I = [q*(math.exp(-b0-b1)-math.exp(-b0))-dC/lam for b0, b1 in B]; m = sum(I)/N; sd = math.sqrt(sum((x-m)**2 for x in I)/(N-1)); return m, sum(x > 0 for x in I)/N, sum(max(-x, 0) for x in I)/N, sum(max(x, 0) for x in I)/N, sd/math.sqrt(N) Run with the same inputs; Python's generator gives different draws for the same seed, an expected value of about minus 0.0811 QALYs, a probability of 0.4094 and a Monte Carlo standard error of about 0.0048.

  Test (Article's expected value and probability reproduced): With the inputs of the first worked example in the Excel columns, the mean of DrawIncNHB rounds to 0.2 and the share of positive draws is 0.2. Expected result TRUE. FALSE shows the hazard ratio applied by division, so that survival is 4 times the hazard ratio, which gives a mean of minus 1.0 and no positive draw. Excel check =AND(ROUND(AVERAGE(DrawIncNHB),4)=0.2,COUNTIF(DrawIncNHB,">0")/COUNT(DrawIncNHB)=0.2)

  Test (Five frozen correlated draws reproduced): With the inputs of HE-EX-PSA-001 in the named cells, QALYPerLifeYear at 0.5, IncCost at 15000, Threshold at 25000 and the five standard normal pairs of the second worked example frozen in columns A and B, the mean of DrawIncNHB rounds to minus 0.1009 and two of the five draws are positive. Expected result TRUE. FALSE shows the correlation left out of the draws, CholOffDiag at zero and CholDiagSecond equal to SDLogHR, which gives a mean of about minus 0.1526, or the upper triangular factor used without transposing, which gives about minus 0.1817. Excel check =AND(ROUND(AVERAGE(DrawIncNHB),4)=-0.1009,COUNTIF(DrawIncNHB,">0")/COUNT(DrawIncNHB)=0.4)

  Common error (Fresh baseline draw for the new option): Drawing a separate log baseline rate for the new option while its log hazard ratio is still drawn from the comparator's first standard normal value breaks the correlation of minus 0.5 between the new option's baseline and its hazard ratio, so the log of its treated rate has a variance of 0.05 instead of 0.03 (derived). From the log-normal moments (derived here, and matched by a 4 million draw simulation) the standard deviation of the incremental net health benefit rises from about 0.487 to about 0.693 QALYs, about 2.03 times the variance, so about twice as many draws are needed for the same Monte Carlo standard error (HE-FM-ABM-003), and the expected value moves from about minus 0.0722 to about minus 0.0467 QALYs. The general rule that both options use the same draw is HE-ER-ENB-003.

  Source: Steuten L, Thavorn K, Cipriano LE, Guertin JR, McQueen B, Schmidler G, Sullivan S, Wailoo A, Sculpher M. Methods to analyze and present uncertainty in economic evaluation to support HTA decisions. HEMA Report 2. Health Economics Methods Advisory, 2026. Section 4.3 (page 28), which states that decision models are non-linear in their inputs, that evaluating the model at the mean value of each input does not necessarily give the expected costs and health effects, and that probabilistic analysis averages model outcomes over the joint distribution of uncertain inputs, and the reporting list in section 6.6 (page 55), which asks for the number of simulations and evidence that the probabilistic results are stable. National Institute for Health and Care Excellence. NICE technology appraisal and highly specialised technologies guidance: the manual (PMG36). Published 31 January 2022, last updated 13 October 2026. Section 4.7.11, which asks for enough simulations to minimise the effect of Monte Carlo error and suggests reviewing the variance around net benefits or ICERs as the number of simulations increases. Owen AB. Monte Carlo theory, methods and examples. Stanford University, 2013. Chapter 5, Random vectors and objects. Section 5.2, equation 5.3 and the Cholesky decomposition on page 9.

  b0_1 = mu_b0 + s_b0 * z1_1; b1_1 = mu_b1 + s_b1 * (rho * z1_1 + sqrt(1 - rho^2) * z2_1); INHB_1 = q * (exp(-b0_1 - b1_1) - exp(-b0_1)) - dC / lambda; b0_2 = mu_b0 + s_b0 * z1_2; b1_2 = mu_b1 + s_b1 * (rho * z1_2 + sqrt(1 - rho^2) * z2_2); INHB_2 = q * (exp(-b0_2 - b1_2) - exp(-b0_2)) - dC / lambda; b0_3 = mu_b0 + s_b0 * z1_3; b1_3 = mu_b1 + s_b1 * (rho * z1_3 + sqrt(1 - rho^2) * z2_3); INHB_3 = q * (exp(-b0_3 - b1_3) - exp(-b0_3)) - dC / lambda; b0_4 = mu_b0 + s_b0 * z1_4; b1_4 = mu_b1 + s_b1 * (rho * z1_4 + sqrt(1 - rho^2) * z2_4); INHB_4 = q * (exp(-b0_4 - b1_4) - exp(-b0_4)) - dC / lambda; b0_5 = mu_b0 + s_b0 * z1_5; b1_5 = mu_b1 + s_b1 * (rho * z1_5 + sqrt(1 - rho^2) * z2_5); INHB_5 = q * (exp(-b0_5 - b1_5) - exp(-b0_5)) - dC / lambda; ENHB = (INHB_1 + INHB_2 + INHB_3 + INHB_4 + INHB_5) / 5; P_CE = ((INHB_1 > 0) + (INHB_2 > 0) + (INHB_3 > 0) + (INHB_4 > 0) + (INHB_5 > 0)) / 5; EL_new = (max(-INHB_1, 0) + max(-INHB_2, 0) + max(-INHB_3, 0) + max(-INHB_4, 0) + max(-INHB_5, 0)) / 5; EL_comp = (max(INHB_1, 0) + max(INHB_2, 0) + max(INHB_3, 0) + max(INHB_4, 0) + max(INHB_5, 0)) / 5; SD_INHB = sqrt(((INHB_1 - ENHB)^2 + (INHB_2 - ENHB)^2 + (INHB_3 - ENHB)^2 + (INHB_4 - ENHB)^2 + (INHB_5 - ENHB)^2) / 4); MCSE = SD_INHB / sqrt(5); N_req = (SD_INHB / e_target)^2

Try this function

Implementations

• Excel

  Correlated pair of PSA draws from named cells in Excel

  With the means in cells named MeanLogRate and MeanLogHR, the standard deviations in SDLogRate and SDLogHR and the correlation in CorrRateHR, the first two formulas give the Cholesky entries in cells named CholOffDiag and CholDiagSecond. StdNormalFirst and StdNormalSecond each hold NORM.S.INV(RAND()), pasted as values when the draw must be kept, and the last two formulas return the draws in cells named DrawLogRate and DrawLogHR. Filling the four draw cells down one row per PSA draw gives the columns used in the computational function below.

  =CorrRateHR*SDLogHR; =SDLogHR*SQRT(1-CorrRateHR^2); =MeanLogRate+SDLogRate*StdNormalFirst; =MeanLogHR+CholOffDiag*StdNormalFirst+CholDiagSecond*StdNormalSecond

Assumptions

• Joint normality of the correlated PSA pair on the scale of estimation

  The two parameters are jointly normal on the scale on which they were estimated, for example the log scale of a rate and a hazard ratio from one regression. The ISPOR-SMDM task force states that the covariance matrix defines the uncertainty in regression coefficients and that multivariate normality is appropriate for the regression's linear predictor. Values on the natural scale are found by back-transforming each draw, and their means then follow the log-normal mean in HE-FM-EV-004, not the exponential of mu_1 or mu_2.

• Correlation strictly inside minus 1 and 1 with a stated basis

  rho is above minus 1 and below 1, so L_22 is above zero and the covariance matrix is positive definite, the condition under which Owen states that a lower triangular Cholesky factor exists. At rho equal to 1 or minus 1 the second parameter is fixed by the first. rho comes from the joint estimation where it exists, and PMG36 section 4.7.19 asks for evidence on correlation to be reflected in the probabilistic analysis and for the assumptions made about correlations to be presented clearly.

• One fresh pair of standard normal values per PSA draw shared by all options

  z_1 and z_2 are independent standard normal values, drawn afresh for each PSA draw and never reused across draws. The pair x_1 and x_2 found from them enters the model for every option in that draw, as the article requires of each complete parameter set.

Worked examples

• Illustrative log baseline rate and log hazard ratio with a correlation of minus 0.5

  An illustrative log annual mortality rate with mean minus 1.3862944 (a rate of 0.25) and standard deviation 0.10 is drawn with an illustrative log hazard ratio with mean minus 0.2231436 (a hazard ratio of 0.8) and standard deviation 0.20, at an illustrative correlation of minus 0.5. For standard normal values of 1 and 0.5 the Cholesky entries are minus 0.1 and about 0.1732051, so the draw is a log rate of minus 1.2862944 (a rate of about 0.2763) and a log hazard ratio of about minus 0.2365411 (a hazard ratio of about 0.7894).

  rho = -0.5; s_2 = 0.20; mu_1 = -1.3862944; s_1 = 0.10; z_1 = 1; mu_2 = -0.2231436; z_2 = 0.5; L_21 = -0.1; L_22 = 0.1732051; x_1 = -1.2862944; x_2 = -0.2365411

• Neine and Curran's correlation of 0.5 on the standard normal scale

  For two standardised parameters with a correlation of 0.5, the second column of the Cholesky matrix printed by Neine and Curran has entries 0.5 and 0.866, the values of L_21 and L_22 here. Standard normal values of 1 and 1 give a first draw of 1 and a second of about 1.3660254.

  rho = 0.5; s_2 = 1; mu_1 = 0; s_1 = 1; z_1 = 1; mu_2 = 0; z_2 = 1; L_21 = 0.5; L_22 = 0.8660254; x_1 = 1; x_2 = 1.3660254

• Zero correlation as the limiting case of the Cholesky PSA step

  With the inputs of the first example but rho equal to zero, L_21 is zero and L_22 equals s_2, so the log hazard ratio draw of about minus 0.1231436 depends on z_2 alone. This limiting case is the usual independent draw.

  rho = 0; s_2 = 0.20; mu_1 = -1.3862944; s_1 = 0.10; z_1 = 1; mu_2 = -0.2231436; z_2 = 0.5; L_21 = 0; L_22 = 0.2; x_1 = -1.2862944; x_2 = -0.1231436

Common errors

• Using an upper triangular Cholesky factor without transposing it

  R's chol returns the upper triangular factor R with R transposed times R equal to the covariance matrix, and Owen warns that failing to transpose it is an error hard to detect from simulation output. For the illustrative inputs mu plus R z gives draws with standard deviations of about 0.1414 and 0.1732 instead of 0.10 and 0.20 and a correlation of about minus 0.71 instead of minus 0.5 (derived from R times R transposed). The factor used must be the lower triangular L with L times L transposed equal to the covariance matrix.

• Drawing jointly estimated parameters as if independent

  Setting rho to zero for the illustrative log rate and log hazard ratio gives the log of the treated rate, their sum, a variance of 0.05 instead of 0.03 (derived), a standard deviation about 29% too large. For the inputs of the computational function below the same omission moves the exact expected incremental net health benefit (derived from the log-normal mean of HE-FM-EV-004) from about minus 0.072 to about minus 0.047 QALYs per person, because the model is not linear in the parameters.

Sources

• Owen on sampling the multivariate normal through a Cholesky factor

  Owen AB. Monte Carlo theory, methods and examples. Stanford University, 2013. Chapter 5, Random vectors and objects. Section 5.2, Multivariate normal and t, which on page 8 samples X equal to mu plus C Z with Z standard normal and C C transposed equal to Sigma (equation 5.3) and lists the conditional normal distribution of one block given the other (property MVN-4), and on page 9 states that when Sigma is positive definite C can be the lower triangular L with L L transposed equal to Sigma, the Cholesky decomposition, that some implementations return the upper triangular R equal to L transposed and failing to transpose it would be an error that could be hard to detect from simulation output, and that it is good practice to compute Sigma minus C C transposed before using C in a simulation.

  View source →

• Neine and Curran on Cholesky decomposition for correlated PSA inputs

  Neine M, Curran D. An algorithm to generate correlated input-parameters to be used in probabilistic sensitivity analyses. Journal of Market Access & Health Policy. 2020;9(1):1857052. Methods, sections Correlation of parameters and Sampling algorithm, which give a correlation matrix with all correlations 0.5 for three variables and its Cholesky decomposition matrix, with rows (1, 0.5, 0.5), (0, 0.866, 0.2887) and (0, 0, 0.8165), state that multiplying independent variables by this matrix gives correlated variables with that correlation matrix, and describe simulating uncorrelated standard normal variables, multiplying them by the Cholesky matrix and transforming the result through the normal distribution function and the inverse beta, gamma or log-normal distribution function to give correlated non-normal PSA inputs.

  View source →

• ISPOR-SMDM task force on covariance and correlation in PSA inputs

  Briggs AH, Weinstein MC, Fenwick EA, Karnon J, Sculpher MJ, Paltiel AD. Model parameter estimation and uncertainty analysis, a report of the ISPOR-SMDM Modeling Good Research Practices Task Force Working Group-6. Medical Decision Making. 2012;32(5):722-732. Section Multivariate estimation and correlation, which states that the covariance matrix defines the uncertainty in regression coefficients and that multivariate normality is appropriate for the regression's linear predictor, and Best Practice VI-10, which states that correlation among parameters should be considered and that jointly estimated parameters, such as those from a regression analysis, have direct evidence on correlation that should be reflected in the analysis.

  View source →

• NICE PMG36 section 4.7.19 on correlation in the probabilistic analysis

  National Institute for Health and Care Excellence. NICE technology appraisal and highly specialised technologies guidance: the manual (PMG36). Published 31 January 2022, last updated 13 October 2026. Section 4.7.19, which asks for evidence about the extent of correlation between individual parameters to be considered and reflected in the probabilistic analysis and for assumptions made about the correlations to be presented clearly.

  View source →

Canonical Identity