Introduction

You are interested in comparing the cost and outcomes over a 20 year time peroid of preventative strategies for postmenopausal women with osteoporosis. You plan to compare two options; Treatment A 10mg/day, and Treatment B 70mg once a week.

The main outcomes measured are incidence of bone fracture and the effect on LYs and QALYs. It is estimated that about 5% of the patients taking Treatment A will have a fracture and 5% of those taking Treatment B will experience a fracture. The cost of treating a fracture is $3,000 and a fracture will decrease life-years by 5. The utility reduction associated with having a fracture is 0.5. Treatment A causes GI complications in 5% of the patients, and Treatment B causes GI complications in 2% of the patients. The cost of treating GI complications is $200, and LYs are decreased by 1. The utility reduction associated with GI complicaations is 0.2.

Treatment A Treatment B
Cost $240 $300
Cost of bone fracture $3000 $3000
Probability of bone fracture 5% 5%
Cost of GI complication $200 $200
Probability of GI complication 5% 2%

Methods

Determining treatment costs

In order to determine the average cost of each treatment we must first determine the cost associated with each branch of the tree. Typically there are \(2^n\) paths or branches in each tree; where \(n\) is the number of side effects due to the treatment. In this example we have two side effects to consider, GI irritation and bone fracture; therefor we have \(n=2\) and as a result \(2^2=4\) branches or paths total.

In the next part of this assignment we will create a function that calculates the average cost of the treatment by first computing the cost associated with each of the four branches and then summing them together.

compute_cost <- function(c, c1, p1, c2, p2){
  # Compute the cost of path 1 SE + GI
  path1 <- (p1*p2)*(c + c1 + c2)
  # Compute the cost of path 2 SE + No GI
  path2 <- (p1*(1-p2))*(c + c1)
  # Compute the cost of path 3 No SE + GI
  path3 <- ((1-p1)*p2)*(c + c2)
  # Compute the cost of path 4 No SE + No GI
  path4 <- ((1-p1)*(1-p2))*(c)
  # Return the average cost for the treatment
  return(path1 + path2 + path3 + path4)
}

We see that for Treatment A where Drug cost = $240, Probability of bone fracture = 5%, Cost of bone fracture = $3000, Probability of GI irritation = 5%, and Cost of GI irritation = $200 we get a cost of $400 dollars. We compute this value in the following way

# Compute cost for Treatment A
drug_cost <- 240 * 20 # Multiply by 20 since we're looking for cost over entire 20yr period
bone_cost <- 3000
bone_prob <- 0.05
gi_cost <- 200
gi_prob <- 0.05

costA <- compute_cost(drug_cost, bone_cost, bone_prob, gi_cost, gi_prob)
costA
## [1] 4960

Now lets compute the cost for Treatment B

# Compute cost for Treatment B
drug_cost <- 300 * 20 # Multiply by 20 since we're looking for cost over entire 20yr period
bone_cost <- 3000
bone_prob <- 0.05
gi_cost <- 200
gi_prob <- 0.02

costB <- compute_cost(drug_cost, bone_cost, bone_prob, gi_cost, gi_prob)
costB
## [1] 6154

Determining lifeyear reductions

This same method can be used to determine the lifeyears associated with a particular treatment by using the lifeyear reductions for each event instead of costs. For example, bone fracture takes away 5 life years, GI upset takes away only 1 life year. We can recycle the compute_cost() function above by replacing the cost in dollars with associated life year reductions. This will give us the average lifeyear reduction for each of the treatments.

# Compute LY reductions for Treatment A
drug_ly_reduction <- 0
bone_ly_reduction <- 5
bone_prob <- 0.05
gi_ly_reduction <- 1
gi_prob <- 0.05

compute_cost(drug_ly_reduction, bone_ly_reduction, bone_prob, gi_ly_reduction, gi_prob)
## [1] 0.3

We see that the life year reduction associated with Treatment A is 0.3 years, since these individuals are expected to live an average of 20 years we can compute the remaining lifeyears associated with each treatment by subtracting the years taken by the treatment from the initial 20 years they had to start with. Thus, for Treatment A we see that patients on average have 19.7 life years. An example is given below

# Calculate LY reduction for Treatment A
ly_red_A <- compute_cost(drug_ly_reduction, bone_ly_reduction, bone_prob, gi_ly_reduction, gi_prob)

# LY remaining after treatment
ly_A <- 20 - ly_red_A
ly_A
## [1] 19.7

The same procedure can be done for Treatment B

# Calculate LY reduction for Treatment B
drug_ly_reduction <- 0
bone_ly_reduction <- 5
bone_prob <- 0.05
gi_ly_reduction <- 1
gi_prob <- 0.02
ly_red_B <- compute_cost(drug_ly_reduction, bone_ly_reduction, bone_prob, gi_ly_reduction, gi_prob)
ly_B <- 20 - ly_red_B
ly_B
## [1] 19.73

Computing average utility for each treatment

Now that we’ve computed that average LY give per treatment we can apply the same process to compute the avereage utily for each treatment. Once we compute the average utily we can calculate the Quality Adjusted Life Years for each treatment and use this as a means of comparing the quality of life the patient can expect to experience while on the particular treatment.

We apply the same compute_cost() function as we did in the previous examples however, in this case we are looking at the utility reduction associated with each event. The reduction in utility associated with GI irritation is 0.2 and the reduction in utility associated with bone fracture is 0.5, there is not reduction in utility associatd with actually taking the medication.

# Compute the average reduction in utility for Treatment A
drug_utility_reduction <- ly_A
bone_utility_reduction <- -0.5
bone_prob <- 0.05
gi_utility_reduction <- -0.2
gi_prob <- 0.05

util_red_A <- compute_cost(drug_utility_reduction,
                           bone_utility_reduction,
                           bone_prob,
                           gi_utility_reduction,
                           gi_prob)
utility_A <- util_red_A # Subtract from 1 since everyone starts with 1 (full utility) prior to treatment
utility_A
## [1] 19.665

We do the same for Treatment B

# Compute the average reduction in utility for Treatment B
drug_utility_reduction <- ly_B
bone_utility_reduction <- -0.5
bone_prob <- 0.05
gi_utilty_reduction <- -0.2
gi_prob <- 0.02

util_red_B <- compute_cost(drug_utility_reduction,
                           bone_utility_reduction,
                           bone_prob,
                           gi_utility_reduction,
                           gi_prob)
utility_B <- util_red_B # Subtract from baseline utility
utility_B
## [1] 19.701

Computing Quality Adjusted Life Years (QALYs)

Now that we have the average utility associated with each of the treatments, we can compute the QALY for each treatment. This is done by multiplyting the expected LYs of each treatment with the average utility associated with each treatment (calculated above).

# Compute QALY for Treatment A
QALY_A <- utility_A
QALY_A
## [1] 19.665
# Compute QALY for Treatment B
QALY_B <- utility_B
QALY_B
## [1] 19.701

Now we are ready to compute incremental costs per LY gained and incremental cost per QALY gained for each treatment.

Computing incremental cost per LY gained

We compute the incremental cost per LY gained according to the formulat below \[ ICER = \frac{Cost_B - Cost_A}{Lifeyear_B - Lifeyear_A}\]

# Compute incremental cost per LY saved
ICER <- (costB - costA)/(ly_B - ly_A)
ICER
## [1] 39800

We now compute the incremental cost per QALY gained similarly \[ICUR = \frac{Cost_B - Cost_A}{QALY_B - QALY_B}\]

ICUR <- (costB - costA)/(QALY_B - QALY_A)
ICUR
## [1] 33166.67
Treatment A Treatment B
Cost 4960 6154
Lifeyears 19.7 19.73
Cost/LY — 39800.00
Cost/QALY — 33166.67

Sensitivity analysis

We are now interested in seeing how the costs and benefits associated with each treatment change with respect to small changes in parameter estimates. The variation in the parameters comes from the uncertainty associated with trying to measure them. For example, utility values associated with certain events are often determined using surveys or a variety of different methods which may give slightly different results Additionally, the probabilities of experiencing certain side effects may have some uncertaintly in the as well. The goal of sensitivity analysis is to determin how the uncertainties alter our finaly decision.

For this part of the excercise lets create another function to run our simulation. The function should manipulate each of the inputs and return the cost. In this example we will simply vary each variable by a percentage of its original value. In practice, you would want to specify the low and high ranges based on values from the literature.

sensitivity <- function(c, c1, p1, c2, p2, var=0.2, n=1e4){
  sampled_params <- data.frame(drug_cost=runif(n, c-(c*var), c+(c*var)),
                                   se1_cost=runif(n,c1-(c1*var), c1+(c1*var)),
                                   se1_prob=runif(n,p1-(p1*var), p1+(p1*var)),
                                   se2_cost=runif(n,c2-(c2*var), c2+(c2*var)),
                                   se2_prob=runif(n,p2-(p2*var), p2+(p2*var)))
  apply(sampled_params, 1, function(x) compute_cost(x[1], x[2], x[3], x[4], x[5]))
}

Now we can call sensitivity() just as we did compute_cost() however it will sample each parameter from a distribution and compute the cost and return the sample, the var argument specifies the range or the variance associated with each parameter, while the n argument specifies the total number of samples to run.

Let’s go ahead and run the simulation and plot the cost distribution for the two treatments

library(ggplot2)
library(devtools)
library(digest)
source_url("https://raw.github.com/low-decarie/FAAV/master/r/stat-ellipse.R")
## SHA-1 hash of file is 5fa2e55de033ed5a45b422ebaa230556bd7e236e
## Loading required package: proto
# Compute cost distributions for the two treatments
n <- 1e3
costA <- sensitivity(240, 3000, 0.05, 200, 0.05, n=n)
costB <- sensitivity(300, 3000, 0.05, 200, 0.02, n=n)
labs <- rep("A",2*n)
labs[n:(2*n)] <- "B"
cost_dist <- data.frame(cost=c(costA, costB), treatment=labs)

ggplot(data=cost_dist, aes(x=cost, fill=treatment)) + 
  geom_histogram(alpha=0.5, position="identity") +
  labs(title="Simulated cost distribution")
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.

Now we can do the same for the expected lifeyears for each treatment

n <- 1e3
lyA <- (20 - sensitivity(0, 5, 0.05, 1, 0.05, n=n))
lyB <- (20 - sensitivity(0, 5, 0.05, 1, 0.02, n=n))
ly_red_dist <- data.frame(Lifeyears=c(lyA, lyB), Treatment=labs)

ggplot(data=ly_red_dist, aes(x=Lifeyears, fill=Treatment)) + geom_histogram(alpha=0.5, position="identity") + 
  labs(title="Simulated lifeyear distributions")
## stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.

Now we can use the simulated data and plot the cost versus benefit (lifeyears)

# Plot distributions with 95% confidence ellipse
df <- cbind(cost_dist[,-2], ly_red_dist)
names(df) <- c("Cost", "Lifeyears", "Treatment")
ggplot(data=df, aes(x=Lifeyears, y=Cost, color=Treatment)) + geom_point(alpha=.7) + 
  stat_ellipse() + 
  labs(title="Cost vs Lifeyears")
## Loading required package: MASS

qaly_A <- 20*(1 - sensitivity(0, 0.5, 0.05, 0.2, 0.05, n=n))
qaly_B <- 20*(1 - sensitivity(0, 0.5, 0.05, 0.2, 0.02, n=n))
df <- data.frame(qaly=c(qaly_A, qaly_B), treatment=labs)

sens_analysis <- cbind(cost_dist[,-2], df)
names(sens_analysis) <- c("Cost", "QALY", "Treatment")

# Plot distributions with 95% confidence ellipse
ggplot(data=sens_analysis, aes(x=QALY, y=Cost, colour=Treatment)) + 
  geom_point(alpha=0.7) +
  stat_ellipse() + 
  labs(title="Cost vs QALY")

Now lets plot the incremental cost agains the incremental benefits

incremental_costs <- costB - costA
incremental_lys <- lyB - lyA
ICER <- data.frame(Cost=incremental_costs, Lifeyears=incremental_lys)
ggplot(data=ICER, aes(Lifeyears, Cost)) +
  geom_hline(yintercept=0) +
  geom_vline(xintercept=0) + 
  geom_point(alpha=0.5) + 
  stat_ellipse(col="green") +
  geom_point(aes(x=mean(Lifeyears), y=mean(Cost), col="orange", size=3)) + 
  theme(legend.position="none") +
  labs(title="Incremental Cost vs Lifeyears")

In the plot above, we see that the incremental cost associated with Treatment B has a mean of about $50 with an incremental lifeyears of about 0.03 years. Also note that the 95% confidence ellipse extends beyond the cost and lifeyear axes.

We can perform the same procedure using incremental Cost and QALYs.

incremental_qaly <- (qaly_B - qaly_A)
ICUR <- data.frame(Cost=incremental_costs, QALYs=incremental_qaly)
ggplot(data=ICUR, aes(QALYs, Cost)) +
  geom_hline(yintercept=0) +
  geom_vline(xintercept=0) + 
  geom_point(alpha=0.5) + 
  stat_ellipse(col="green") +
  geom_point(aes(x=mean(QALYs), y=mean(Cost), col="orange", size=3)) + 
  theme(legend.position="none") +
  labs(title="Incremental Cost vs QALYs")

From the plot above it appears that we get an incremental increase in QALY of about 0.0125 years however this comes at with an incremental cost of about $50. Again, the 95% confidence ellipse crosses both the Cost and QALY axis.