A Common Issue

• We have defined our underlying model (either from bottom up or through backwards conversion) but need to define PSA distributions for model parameters.

• Model draws on literature-based parameters that are reported as means, standard deviations, interquartile ranges, 95% confidence intervals, etc.

A Common Issue

• Straightforward to obtain base case values from literature.

• But how can we define PSA distributions based on limited information?

PSA Distributions

Example

• Cost parameter reported in literature has interquartile range of $300-$750.
• What are the parameters of a PSA distribution such that the 25th percentile is $300 and the 75th percentile is $750?

Some Options

• Analytic formulas for some distributions (normal, gamma, beta).
• Formulas take as their inputs the two values (\(x_1\), \(x_2\)) and their associated quantiles (\(p_1\),\(p_2\)).
• Formulas return the PSA distribution parameters that match up to these values.

Some Options

• ParameterSolver implemented in Windows, Link

• See Cook (2010) and Vasco Grilo’s blog post for more.

• The next few slides provide you with nearly all the tools you’ll need, however.

Example: Cost PSA Distribution

• Suppose the interquartile range for a key cost variable is [300,750]
• We may have obtained this from the literature, from tabulating cost data, or from expert opinions.

Example: Cost PSA Distribution

x1 = 300
p1 = 0.25

x2 = 750
p2 = 0.75

Analytic Solution: Uniform PSA

For a uniform distribution with minimum \(a\) and maximum \(b\)

\[ a = \frac{p_2x_1 - p_1x_2}{p_2-p_1} \] \[ b = \frac{(1-p_1)x_2-(1-p_2)x_1}{p_2-p_1} \]

Analytic Solution: Uniform PSA

a = ((p2*x1) - (p1 * x2)) / (p2 - p1)
a

[1] 75

b = ((1 - p1)*x2 - (1 - p2)*x1) / (p2 - p1)
b

[1] 975

Analytic Solution: Uniform PSA

qunif(0.25, min = 75, max = 975)

[1] 300

qunif(0.75, min = 75, max = 975)

[1] 750

Analytic Solution: Normal PSA

\[ \sigma = \frac{x_2 - x_1}{\Phi^{-1}(p_2)-\Phi^{-1}(p_1)} \]

\[ \mu = \frac{x_1\Phi^{-1}(p_2)-x_2\Phi^{-1}(p_1)}{\Phi^{-1}(p_2)-\Phi^{-1}(p_1)} \]

Analytic Solution: Normal PSA

mu = (qnorm(p2)*x1 - qnorm(p1)*x2) / (qnorm(p2)-qnorm(p1))
mu

[1] 525

sigma = (x2 - x1) / (qnorm(p2) - qnorm(p1))
sigma

[1] 333.59

Analytic Solution: Normal PSA

qnorm(0.25, mean = 525, sd = 333.59)

[1] 300

qnorm(0.75, mean = 525, sd = 333.59)

[1] 750

Analytic Solution: Lognormal PSA

• Just take log of \(x_1\) and \(x_2\)

\[ \sigma = \frac{\ln(x_2) - \ln(x_1)}{\Phi^{-1}(p_2)-\Phi^{-1}(p_1)} \]

\[ \mu = \frac{\ln(x_1)\Phi^{-1}(p_2)-\ln(x_2)\Phi^{-1}(p_1)}{\Phi^{-1}(p_2)-\Phi^{-1}(p_1)} \]

Analytic Solution: Lognormal PSA

mu = (qnorm(p2)*log(x1) - qnorm(p1)*log(x2)) / (qnorm(p2)-qnorm(p1))
mu

[1] 6.1619

sigma = (log(x2) - log(x1)) / (qnorm(p2) - qnorm(p1))
sigma

[1] 0.67925

Analytic Solution: Lognormal PSA

qlnorm(0.25, mean = 6.1619, sd = 0.67925)

[1] 299.99

qlnorm(0.75, mean = 6.1619, sd = 0.67925)

[1] 749.98

Analytic Solution: Gamma PSA

• A bit more involved as it involves finding the root of a function.

Analytic Solution: Gamma PSA

gamma_fn <- function(alpha) {
    x1*qgamma(p2,shape = alpha, scale =1) - x2 * qgamma(p1, shape = alpha, scale = 1)
}
curve(gamma_fn, xlim = c(1,10), col = "blue", lwd = 1.5, lty=2)
abline(a=0,b=0)

Analytic Solution: Gamma PSA

gamma_fn <- function(alpha) {
    x1*qgamma(p2,shape = alpha, scale =1) - x2 * qgamma(p1, shape = alpha, scale = 1)
}
curve(gamma_fn, xlim = c(1,10), col = "blue", lwd = 1.5, lty=2)
abline(a=0,b=0)

Analytic Solution: Gamma PSA

gamma_fn <- function(alpha) {
    x1*qgamma(p2,shape = alpha, scale =1) - x2 * qgamma(p1, shape = alpha, scale = 1)
}
curve(gamma_fn, xlim = c(1,10), col = "blue", lwd = 1.5, lty=2)
abline(a=0,b=0)

Analytic Solution: Gamma PSA

alpha_ <- uniroot(gamma_fn,c(2,4))$root
alpha_

[1] 2.4558

calc_beta <- function(x1,p1,alpha) {
    x1 / qgamma(p1,alpha,1)
}

beta_ <- calc_beta(x1 = x1,  p1 = p1, alpha = alpha_)
beta_

[1] 230.16

Analytic Solution: Gamma PSA

qgamma(0.25,shape = 2.4558, scale = 230.16)

[1] 300

qgamma(0.75,shape = 2.4558, scale = 230.16)

[1] 750

Alternative: Optimization

\[\min_{\vec{\theta} \in R}\,\, (F(x_1|\vec{\theta})-p_1)^2+(F(x_2|\vec{\theta})-p_2)^2\]

Alternative: Optimization

• Requires cumulative distribution function (CDF), \(F\).
• See Hans W Borchers’ great slides for tips on optimization in R.

Some Advice

• We have found that while analytically correct, this optimization formula is not numerically stable.
• Transfinite scaling of the probabilities from CDFs \(F(.)\) stabilizes optimization.

Some Advice

\[\min_{\vec{\theta} \in R}\,\, (g(F(x_1|\vec{\theta}))-p_1)^2+(g(F(x_2|\vec{\theta}))-p_2)^2\]

Stable Optimization

• For probabilities, this is the \(\text{logit}(p)=\log \frac{p}{1-p}=\log p - \log (1-p)\).
• Tweaking parameters to converge can still happen.
• Example of theory versus practice.

Stable Optimization

\[\min_{\vec{\theta} \in R} \sum_{i \in 1,2}\,\, (\text{logit} (F(x_i|\vec{\theta}))-\text{logit} (p_i))^2\]

Transfinite Example, Part 1

# Transfinite scaling
Tf <- function(x, shape, rate)
    pgamma(x,shape,rate,log=TRUE) - 
    pgamma(x,shape,rate,log=TRUE, lower.tail = FALSE)

# Function to minimize
norm <- function(x1, p1, x2, p2, shape,rate)
    (Tf(x1, shape, rate)-(log(p1)-log(1-p1)) )^2 +
    (Tf(x2, shape, rate)-(log(p2)-log(1-p2)) )^2

# Bundle it all together into a single function.
fn <- function(x) norm(x1, p1, x2, p2, x[1], x[2])

Transfinite Example, Part 2

# Run general-purpose optimization on the function.
gamma_optim <- 
  optim(c(0.5, 0.1), # initial parameter guesses
    fn,
    gr = function(x) pracma::grad(fn, x),
    lower=c(-Inf, 1e-5),
    method = "L-BFGS-B",
    control=list(factr=1e-10, maxit=100))$par

# Note: gamma_optim$par[2] is 1/beta
gamma_optim

[1] 2.4558556 0.0043449

Transfinite Example

# Analytic formula
alpha_

[1] 2.4558

# Optimization
gamma_optim[[1]]

[1] 2.4559

# Analytic formula
beta_

[1] 230.16

# Optimization
1/gamma_optim[[2]]

[1] 230.16

Optimization for Other Distributions

• Just swap in the distribution function for the distribution you want to match to (e.g., pbeta rather than pgamma).
• Also make sure the supplied parameter names match the distribution you’re aiming for.

Optimization for Other Distributions

Tf <- function(x, shape, rate)
    pgamma(x,shape,rate,log=TRUE) - 
    pgamma(x,shape,rate,log=TRUE, lower.tail = FALSE)

Optimization for Other Distributions

Tf <- function(x, shape, rate)
    pgamma(x,shape,rate,log=TRUE) - 
    pgamma(x,shape,rate,log=TRUE, lower.tail = FALSE)

Summary on Distribution Fitting

• Use analytical formulas if they exist.
• Optimize on \(\text{logit}\) scale.
• Again, avoid \(\log\) on probabilities, use log=TRUE.
• Plot results for visual check.

Sensitivity Analysis for Decision Models

• It is standard in CEA analyses to explore sensitivity of model outputs and outcomes to variation in model inputs.
• Many decision problems are complicated by dependencies between model inputs (e.g., costs, utilities, etc.).
• These correlations may have meaningful impacts on model results and sensitivity.

HIV Model

In the HIV model discussed earlier, patients accrue two types of costs:

1. Direct medical costs associated with each health state (Healthy, LowCD4, AIDS)
2. Community medical costs associated with each health state.

HIV Model

• These costs are likely correlated with each other in meaningful ways.
• Similarly, other model parameters (e.g., state utility payoffs, hazard rates for treatment strategies, etc.) may also be correlated.

HIV Model

• How can we reflect this correlation in n-way sensitivity analyses, such as PSAs?

What You Were Taught

• Nearly all PSA exercises draw from independent distributions.
• To fully explore the space of correlated parameters, you need a lot of PSA draws!

What You Were Taught

• Independent PSA draws of two model parameters
  • \(c_{direct} \sim \text{gamma}(2.5, 1/230)\)
  • \(c_{community} \sim \text{gamma}(1.5, 1/150)\)

What We’ll Teach You

• \(\text{cor}(c_{direct},c_{community}) = 0.8\)

What We’ll Teach You

• Red points: independent PSA draws
• Black points: correlated PSA draws

Copula-Based PSA Sampling

• We can easily draw correlated PSA samples via copulas.

• Copulas allow us to sample from the joint distribution of correlated parameters.

Copulas

• Copula: a multivariate cumulative distribution function with uniform marginals.

• Intuition: the binding “glue” between random variables.

Sklar (1959)

• Sklar (1959) demonstrates that the joint distribution of two random variables can be represented as

Sklar (1959)

\[ \mathrm{F}_{X_{1}, X_{0}}\left(x_1, x_0\right)=C_{X_{1}, X_{0}}\left(\mathrm{~F}_{X_{1} }\left(x_1\right), \mathrm{F}_{X_{0} }\left(x_0\right)\right) \]

Copula-Based PSA Draws

• We’ll draw on copulas so we can sample our PSA as the joint distribution of uncertain parameters.

• You can always just set the correlations to zero if you want independent draws!

\[ F(x_1,...,x_d) = C \big ( F_1(x_1), ..., F_d(x_d) \big ) \]

How Do We Sample?

• Recall that the copula specifies multivariate CDF with uniform marginals.
• What if our PSA distributions are not uniform (e.g., costs \(\sim\) gamma, lognormal, etc.)?
• Easy! Use the inverse transform method.
• Idea: randomly sample from quantiles in the PSA distribution.

How Do We Draw PSA Values

Normal Distribution: PDF to CDF

Normal Distribution: CDF

Drawing from an Arbitrary PSA Distribution

1. Define PSA distribution parameters (e.g., \(\mu,\sigma\) if normal, etc.)
2. Draw a uniform random number between 0 and 1.
3. Feed this number through the quantile function for the specified distribution (e.g., qnorm(), qgamma(), etc.). This returns the value that maps to the randomly drawn quantile.
4. Do this for as many PSA samples as you need!

Copula-Based Sampling

• Sampling from a copula adds one more step to this process.
• We must also define a correlation matrix Sigma (\(\Sigma\)).
• We then draw uniform random numbers from a multivariate uniform distribution.
• Once drawn, we can feed each PSA variable back through its repective inverse CDF to obtain the final PSA sample.

Copula-Based Sampling

1. Define PSA distributions for each of \(n\) uncertain parameters.
2. Define an \(n \times n\) correlation matrix \(\Sigma\) for the joint distribution of the PSA parameters.
3. Draw an \(n\)-column multivariate uniform sample.
4. Apply the inverse transform method (as necessary) to each column.
5. Do this for as many PSA samples as you need!

Putting it All Together

• IQR Community Care Costs: 750 to 1250
  • cost_cc \(\sim\) lognormal(mu_cc, sigma_cc).
• IQR Direct Medical Care Costs: 1,000 - 3,000
  • cost_dc \(\sim\) gamma(shape_dc, rate_dc)
• COV(Community cost, Direct cost) = 0.8

1. Define PSA distributions

mu_cc = (qnorm(0.75)*log(750) - qnorm(0.25)*log(1250)) / (qnorm(0.75)-qnorm(.25)); mu_cc

[1] 6.8755

sigma_cc = (log(1250) - log(750)) / (qnorm(.75) - qnorm(.25)); sigma_cc

[1] 0.37868

1. Define PSA distributions

# Run general-purpose optimization on the function.
gamma_optim <- 
  optim(c(0.5, 0.1), # initial parameter guesses
    fn,
    gr = function(x) pracma::grad(fn, x),
    lower=c(-Inf, 1e-5),
    method = "L-BFGS-B",
    control=list(factr=1e-10, maxit=100))$par

shape_dc = gamma_optim[[1]]; shape_dc

[1] 1.7904

rate_dc = gamma_optim[[2]]; rate_dc

[1] 0.00080927

2. Define correlation matrix \(\Sigma\)

rho_cc_dc = 0.8
sigma <- matrix(c(1,rho_cc_dc,rho_cc_dc,1),nrow = 2, byrow=TRUE,
                dimnames=list(c("cost_cc","cost_dc"),c("cost_cc","cost_dc")))
sigma

        cost_cc cost_dc
cost_cc     1.0     0.8
cost_dc     0.8     1.0

3. Draw multivariate uniform sample.

• To get multivariate uniform, first draw from multivariate normal (MASS::mvrnorm()).
• Then, apply the values to the normal CDF (pnorm()) to get the quantiles.

3. Draw multivariate uniform sample.

# Define PSA Sample Size
PSA_sample_size = 1e6
U <- pnorm(MASS::mvrnorm(n = PSA_sample_size, Sigma = sigma, mu = c(0,0)))