Coding the Cost Effectiveness Acceptability Curve (CEAC) in R

Introduction

The cost effectiveness acceptability curve (CEAC) is a tool used to describe the output of a probabilistic sensitivity analysis conducted on a model used in economic evaluations of health technologies. It communicates the probability of cost effectiveness conditional on our willingness to pay threshold for each unit of effectiveness or QALY gained. while this blog post will not provide the best treatment of theory behind the CEAC it will provide you some information on how to compute them in R.

Preliminaries: Simulating Cost Effectiveness Data.

For the computing of the CEAC, I have chosen to simulate Data on 2000 treated and untreated patients (1000 in each group). If you are working with (real) patient level data using the net benefits regression framework would be more appropriate, but we can view this crude simulation as communicating the output from some probabilistic sensitivity analysis of a model. While setting the seed is not required, it does help with replicating the picture you see above.

I have also included code for the CE plane which is prodc

#Load Packages
library('ggplot2')
library('ggthemes')

#Set seed
set.seed(1234)

#Simulating Data
Control_QALYs<-rbeta(1000,8,2)
Treatment_QALYs<-rbeta(1000,8,2)
#Simulate Costs
Control_Costs<-rnorm(1000,1000,200)
Treatment_Costs<-rnorm(1000,1500,200)
#Create Incremental Costs
D_Costs<-Treatment_Costs-Control_Costs
D_QALYs<-Treatment_QALYs-Control_QALYs

#CE Data
CEdata<-data.frame(D_Costs,D_QALYs)

#Plot CE Plane
ggplot(CEdata,aes(y=D_Costs))+
  geom_point(aes(x=D_QALYs),color="blue")+
  geom_vline(xintercept=0,size=1)+
  geom_hline(yintercept=0,size=1)+
  xlim(-1,1)+
  xlab("Incremental QALYs")+
  ylab("Incremental Costs")+
  ggtitle("The CE Plane")

Coding the CEAC

To actually code your CEAC we need to convert the measures of incremental costs and effects into incremental net benefits. I do this by first creating a list of willingness to pay thresholds which we will be considering, creating an array for storing the data and then using a for loop for populating this array

#Create Threshold Vector
T_Names<-c("T0","T5000","T10000","T15000","T20000","T250000")
Threshold<-c(0,5000,10000,15000,20000,25000)

#Create Dataframe
id<-seq(1,1000)
df<-matrix(NA,nrow=1000,ncol=6,dimnames=list(NA,WTP=T_Names))

#For loop for Incremental Net benefits Calculation
for(i in 1:6){
  df[,i]<-Threshold[i]*D_QALYs-D_Costs  
}

#ProbCE Calculations
ProbCE_0<-sum(df[,1]>=0)/1000
Prob_CE5k<-sum(df[,2]>=0)/1000
Prob_CE10k<-sum(df[,3]>=0)/1000
Prob_CE15k<-sum(df[,4]>=0)/1000
Prob_CE20k<-sum(df[,5]>=0)/1000
Prob_CE25k<-sum(df[,6]>=0)/1000

#StoreValues
Prob_CE<-c(ProbCE_0,Prob_CE5k,Prob_CE10k,Prob_CE15k,Prob_CE20k,Prob_CE25k)

#CEAC Dataframe
CEAC_Data<-data.frame(Threshold,Prob_CE)
CEAC_Data

#Plot CEAC
ggplot(CEAC_Data, aes(y=Prob_CE))+
  geom_line(aes(x=Threshold),color="blue",size=1.0)+
  geom_point(aes(x=Threshold),color="black", size=2.0)+
  ylim(0,1)+
  xlab("Willingness to Pay per QALY")+
  ylab("Prob of CE")+
  ggtitle("The Cost Effectiveness Acceptability Curve")+
  theme_economist_white()

Concluding Remarks

This blog post has demonstrated a way to compute the cost effectiveness acceptability curve in R from scratch. There is likely a smarter way to compute each probability of cost effectiveness at each threshold in the form of a loop but the simplest solution is to calculate each probability directly.