Incorporating New Evidence and Uncertainty

Backwards Conversion, Solving for PSA Parameters and Copula-Based PSA Sampling

Learning Objectives

  1. Back-convert an existing transition probability matrix to incorporate new health states, strategies, and evidence.
  2. Solve for probabilisitic sensitivity analysis (PSA) distribution parameters.
  3. Sample correlated PSA distributions using copulas.

1. Backwards Conversion

  • Lectures 2 and 3 emphasized the importance of the transition rate matrix as the central “hub” of a Markov model.

Goal

What Rarely Happens

Adapting a Model Requires a Transition Rate Matrix

1. Backwards Conversion

1. Backwards Conversion

  • If starting a model from scratch, can simply define or estimate rates, and use them to construct the rate matrix.

1. Backwards Conversion

  • What if we have a model that is already defined in terms of a transition probability matrix?
  • How can we convert back to a rate matrix to add new health states, transition states, accumulators, etc. as needed?
  • Also useful if we want to keep everything the same, but change the model time step (e.g., see Chhatwal, Jayasuriya, and Elbasha (2016))

1. Backwards Conversion

  • We will next show you several ways to work backwards.
  • Boils down to solving for the continuous “generator matrix” for the observed transition probability matrix.

A Word of Warning

  • As you’ll see, it’s not always going to work perfectly.
  • If the original transition probability matrix was defined incorrectly (e.g., no jumpover probabilities), the generator matrix may not exist.
  • Identifiability is a more general issue, however.

Working Example: HIV Model

Transition probability matrix:

        Healthy LowCD4  AIDS Death
Healthy   0.721  0.202 0.067 0.010
LowCD4    0.000  0.581 0.407 0.012
AIDS      0.000  0.000 0.750 0.250
Death     0.000  0.000 0.000 1.000
  • From ‘Decision Modelling for Health Economic Evaluation’ by Briggs, Claxton, Sculpher.

Three Methods to Solve for the Generator

  • A. Compute the matrix logarithm of the transition probability matrix \(P\).

  • B. Maximum likelihood-based approach.

  • C. Convert the supplied transition probabilities to rates one-by-one.

  • (You can always try 2+ of these methods to see how they compare.)

A. Matrix Logarithm of P

  • In mathematical terms, the generator matrix is the matrix logarithm of the transition probability matrix.
  • A matrix has a logarithm if and only if if it is invertible.

A. Matrix Logarithm of P

  • In mathematical terms, the generator matrix is the matrix logarithm of the transition probability matrix.
  • A matrix has a logarithm if and only if if it is invertible.

\[ R = \log P \]

A. Matrix Logarithm of P

  • The \(\log\) can be found using spectral or eigenvalue decomposition.
  • If \(V\) is a matrix where each column is an eigenvector of \(P\), then

A. Matrix Logarithm of P

  • The \(\log\) can be found using spectral or eigenvalue decomposition.
  • If \(V\) is a matrix where each column is an eigenvector of \(P\), then

\[A' = V^{-1} A V\]

A. Matrix Logarithm of P

  • The \(\log\) can be found using spectral or eigenvalue decomposition.
  • If \(V\) is a matrix where each column is an eigenvector of \(P\), then

\[A' = V^{-1} A V\]

\[\log P = V (\log A') V^{-1}\]

A. Matrix Logarithm of P

V  <- eigen(mP)$vectors
iV <- solve(V)
Ap <- iV %*% mP %*% V
lAp <- diag(log(diag(Ap)), nrow(Ap), ncol(Ap))
R  <- V %*% lAp %*% iV
R[abs(R) < 1e-6 ] <- 0
dimnames(R) = dimnames(mP)
round(R,3)

A. Matrix Logarithm of P

V  <- eigen(mP)$vectors
iV <- solve(V)
Ap <- iV %*% mP %*% V
lAp <- diag(log(diag(Ap)), nrow(Ap), ncol(Ap))
R  <- V %*% lAp %*% iV
R[abs(R) < 1e-6 ] <- 0
dimnames(R) = dimnames(mP)
round(R,3)
        Healthy LowCD4   AIDS  Death
Healthy  -0.327  0.311  0.002  0.013
LowCD4    0.000 -0.543  0.615 -0.072
AIDS      0.000  0.000 -0.288  0.288
Death     0.000  0.000  0.000  0.000

A. Matrix Logarithm of P

Healthy LowCD4 AIDS Death
Healthy -0.327 0.311 0.002 0.013
LowCD4 0 -0.543 0.615 -0.072
AIDS 0 0 -0.288 0.288
Death 0 0 0 0
  • Note there is a negative rate from LowCD4 Death!
  • This is a transition from Death LowCD4 and should be moved to the other side of the matrix.

A. Matrix Logarithm of P

Healthy LowCD4 AIDS Death
Healthy -0.327 0.311 0.002 0.013
LowCD4 0 -0.543 0.615 -0.072
AIDS 0 0 -0.288 0.288
Death 0 0 0 0
  • The negative rate implies an identifiability issue.
  • Note rates from HealthyAIDS, HealthyDeath are relatively small, implying these are from statistical noise in observation.

A. Matrix Logarithm of P

  • expm::logm Higham (2008) method returns same as eigenvalue method.
  • Tweaking rates to get a proper model is ad-hoc and difficult.

A. Matrix Logarithm of P

  • Higham (2008) method returns same as eigenvalue method. - expm::logm()
  • Tweaking rates to get a proper model is ad-hoc and difficult.
R = expm::logm(mP)
        Healthy LowCD4   AIDS  Death
Healthy  -0.327  0.311  0.002  0.013
LowCD4    0.000 -0.543  0.615 -0.072
AIDS      0.000  0.000 -0.288  0.288
Death     0.000  0.000  0.000  0.000

B. Multi-State Model Based Approach

  • Imposition of structural assumptions and fitting to pseudo-data derived from original Markov model can result in reasonable rate estimates.
  • We could assume that a patient with HIV at this point in time only gets sicker and external causes of death are negligible. Healthy LowCD4 AIDS Death constrains model.
  • We use the reported data from the original model at year 1.

B. Multi-State Model Based Approach

Steps:

  1. Run the original model for 1+ cycles to obtain the Markov trace.
  2. Construct psuedo-data from the resulting trace based on a cohort with reasonable size (e.g., 1,000 patients)
  3. Estimate transition hazards in the pseudo-data based on a multistate model.1
  4. Use the estimated transition hazards to construct the rate matrix.

1. Run the original model

tr0 <- # Initial state occupancy.
  c("Healthy" = 1000,"LowCD4" = 0, "AIDS" = 0, "Death" = 0)
tr1 <- # State occupancy after one cycle.
  tr0 %*% mP 

tr <- # Bind together into a data frame. 
    rbind.data.frame(tr0, tr1) %>%
    mutate(cycle = c(0, 1)) %>% 
    select(cycle,everything())
tr
  cycle Healthy LowCD4 AIDS Death
1     0    1000      0    0     0
2     1     721    202   67    10

2. Construct psuedo-data

  • Idea: create data for 1,000 simulated “patients” followed for two cycles.
  • Counts in the Markov trace govern state occupancy in each cycle.
    • In first cycle, all 1,000 are in Healthy state.
    • In second cycle, 721 remain in Healthy while 202 are in LowCD4, etc.

2. Construct psuedo-data

tr <- 
    rbind.data.frame(tr0, tr1) %>%
    mutate(cycle = c(0, 1)) %>% 
    gather(state,count,-cycle) %>% 
    mutate(count = round(count,0)) %>% 
    arrange(cycle) %>% 
    lapply(.,rep,.$count) %>% 
    cbind.data.frame() %>% 
    as_tibble() %>% 
    mutate(state = factor(state,
      levels = c("Healthy","LowCD4","AIDS","Death"))) %>% 
    arrange(cycle,state) %>% 
    mutate(ptnum = rep(1:1000,2)) %>% 
    select(-count) %>% 
    arrange(ptnum,cycle) %>% 
    mutate(state = as.numeric(state)) %>% 
    select(ptnum,cycle,state)
# A tibble: 2,000 × 3
   ptnum cycle state
   <int> <dbl> <dbl>
 1     1     0     1
 2     1     1     1
 3     2     0     1
 4     2     1     1
 5     3     0     1
 6     3     1     1
 7     4     0     1
 8     4     1     1
 9     5     0     1
10     5     1     1
# ℹ 1,990 more rows

3. Fit Multistate Model

  • Next, fit a multistate model to the pseudo-data.
  • The multistate model will use likelihood-based methods to estimate transition intensities among health states in the data.
  • These transition intensities are the rates that you can use in the rate matrix!

3. Fit Multistate Model

  1. Define the state table for the multistate model
library(msm)
statetable.msm(state, ptnum, data=tr)
    to
from   1   2   3   4
   1 721 202  67  10
  1. Check that cell transition counts match the trace.
# Markov trace
rbind.data.frame(tr0, tr1) %>%
    mutate(cycle = c(0, 1)) %>% 
    select(cycle,everything())
  cycle Healthy LowCD4 AIDS Death
1     0    1000      0    0     0
2     1     721    202   67    10

3. Fit Multistate Model

  • Initial rate guesses are based on the eigenvalue method, i.e., result of logm(mP).
Q.init <- rbind(rbind(c(0, 0.3, 0,   .0001),    
                      c(0, 0,   0.6, 0.01),    
                      c(0, 0,   0,   0.3),  
                      c(0, 0,   0,   0)))
dimnames(Q.init) = dimnames(mP)
  • Fit the model
hiv.msm <- 
  msm(state ~ cycle, subject = ptnum, data = tr, qmatrix = Q.init)

3. Fit Multistate Model

hiv.msm

Call:
msm(formula = state ~ cycle, subject = ptnum, data = tr, qmatrix = Q.init)

Maximum likelihood estimates

Transition intensities
                  Baseline                          
Healthy - Healthy -0.32712 ( -3.680e-01, -2.907e-01)
Healthy - LowCD4   0.32702 (  1.354e-01,  7.898e-01)
Healthy - Death    0.00010 (  0.000e+00,        Inf)
LowCD4 - LowCD4   -0.64477 ( -1.139e+01, -3.650e-02)
LowCD4 - AIDS      0.63465 (  9.115e-06,  4.419e+04)
LowCD4 - Death     0.01012 ( 8.238e-299, 1.242e+294)
AIDS - AIDS       -0.34837 ( -3.805e+40, -3.189e-42)
AIDS - Death       0.34837 (  3.189e-42,  3.805e+40)

-2 * log-likelihood:  1572.2 

4. Construct the Rate Matrix

Rmsm <- msm::qmatrix.msm(hiv.msm, ci = "none")
round(Rmsm,3)
        Healthy LowCD4   AIDS Death
Healthy  -0.327  0.327  0.000 0.000
LowCD4    0.000 -0.645  0.635 0.010
AIDS      0.000  0.000 -0.348 0.348
Death     0.000  0.000  0.000 0.000

5. Embed the Transition Probability Matrix

Multistate-Model:

round(expm(Rmsm),3)
        Healthy LowCD4  AIDS Death
Healthy   0.721  0.202 0.067 0.010
LowCD4    0.000  0.525 0.388 0.088
AIDS      0.000  0.000 0.706 0.294
Death     0.000  0.000 0.000 1.000

Original Model Matrix:

round(mP,3)
        Healthy LowCD4  AIDS Death
Healthy   0.721  0.202 0.067 0.010
LowCD4    0.000  0.581 0.407 0.012
AIDS      0.000  0.000 0.750 0.250
Death     0.000  0.000 0.000 1.000

Markov Trace Comparison

# Original Model
c(1000, 0, 0, 0) %*% mP
     Healthy LowCD4 AIDS Death
[1,]     721    202   67    10
# MSM-Based Model
c(1000, 0, 0, 0) %*% pmatrix.msm(hiv.msm, t=1)
     Healthy LowCD4 AIDS Death
[1,]     721    202   67    10
  • Resulting expectation after 1 cycle is nearly identical.
  • The msm method has bits down at 4th decimal. Round off made them identical.
  • Only a tiny bit of error (721 vs. 721.0005) makes log numerically unstable converting from probability to rate.

Key Takeaways

  • By imposing some constraints on the underlying transitions, we were able to yield a generator matrix that makes sense!
  • Critically, our multistate model-based generator matrix closely approximates the original transition probability matrix, but does not imply negative rates that bring people back from death.

Key Takeaways

Multistate-Model:

round(logm(mP),3)
       [,1]   [,2]   [,3]   [,4]
[1,] -0.327  0.311  0.002  0.013
[2,]  0.000 -0.543  0.615 -0.072
[3,]  0.000  0.000 -0.288  0.288
[4,]  0.000  0.000  0.000  0.000

Original Model Matrix:

round(Rmsm,3)
        Healthy LowCD4   AIDS Death
Healthy  -0.327  0.327  0.000 0.000
LowCD4    0.000 -0.645  0.635 0.010
AIDS      0.000  0.000 -0.348 0.348
Death     0.000  0.000  0.000 0.000

Key Takeaways

  • With a reasonable generator matrix defined, we can now augment the model as we see fit:

    • Add additional health states.
    • Add new strategies with evidence on efficiacy drawn from the literature (e.g., hazard rates).
    • Add accumulators and transition states.
    • Change the cycle length.
  • Once done with the above, just embed the new matrix.

2. Solving for PSA Distributions

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

Parameter Type Distribution
Probability beta
Rate gamma
Utility weight beta
Right skew (e.g., cost) gamma, lognormal
Relative risks or hazard ratios lognormal
Odds Ratio logistic

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

  • 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.
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

  • A bit more involved as it involves finding the root of a function.
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

  • Root (i.e., point where the function crosses the zero line) seems to be between 2 and 4.
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

  • We’ll search in this range for our \(\alpha\) value.
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

  • General solution that can be used to solve for parameters of many common PSA distributions (e.g., beta, gamma, etc.).

\[\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

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

\[\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

  • 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.

\[\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

  • Optimization returns \(\alpha\) and \(1/\beta\) for the gamma distribution 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

  • 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.
Tf <- function(x, shape, rate)
    pgamma(x,shape,rate,log=TRUE) - 
    pgamma(x,shape,rate,log=TRUE, lower.tail = FALSE)

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.
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.

Copula PSA Sampling

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)

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

\[ \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

  • Requires drawing a random number between 0 and 1 (the quantile)
  • Requires an inverse cumulative distribution function \(F^{-1}(.)\)
  • You can do this easily in R or Excel!

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

  • IQR Community Care Costs: 750 to 1250
    • cost_cc \(\sim\) lognormal(mu_cc, sigma_cc).
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

  • IQR Direct Medical Care Costs: 1,000 - 3,000
    • cost_dc \(\sim\) gamma(shape_dc, rate_dc)
# 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.

  • To get multivariate uniform, first draw from multivariate normal (MASS::mvrnorm()).
  • Then, apply the values to the normal CDF (pnorm()) to get the quantiles.
# Define PSA Sample Size
PSA_sample_size = 1e6
U <- pnorm(MASS::mvrnorm(n = PSA_sample_size, Sigma = sigma, mu = c(0,0)))

4. Inverse Transform

# Inverse transform using parameters solved for earlier
cost_cc <- qlnorm(U[,1], mu_cc, sigma_cc)
cost_dc <- qgamma(U[,2], shape_dc, rate_dc)

# Make sure the marginal distributions check out
quantile(cost_cc,c(0.25,0.75))
    25%     75% 
 749.62 1249.57 
quantile(cost_dc,c(0.25,0.75))
   25%    75% 
1000.3 2997.7 
# Check that correlation is there
cor(cost_cc,cost_dc)
[1] 0.78337

Your Turn!

References

Chhatwal, Jagpreet, Suren Jayasuriya, and Elamin H. Elbasha. 2016. “Changing Cycle Lengths in State-Transition Models: Challenges and Solutions.” Medical Decision Making 36 (8): 952–64. https://doi.org/10.1177/0272989X16656165.
Cook, John D. 2010. “Determining Distribution Parameters from Quantiles.” Working paper.
Higham, Nicholas J. 2008. Functions of Matrices: Theory and Computation. SIAM.