[This article was first published on
R Works, and kindly contributed to
R-bloggers]. (You can report issue about the content on this page
here)
Want to share your content on R-bloggers?
click here if you have a blog, or
here if you don't.
In a previous post, we presented an rjags version of a Bayesian model from Chapter 2 of the textbook: Evidence Synthesis for Decision Making in Healthcare by Welton et al. (2012). In this post, we continue our project to present rjags versions of the WinBugs models in the text by presenting the economic decision model from Chapter 3. Note that to run the code below, you must first install the JAGS program which is available from http://mcmc-jags.sourceforge.net/, and the rjags package is available from https://cran.r-project.org/web/packages/rjags/index.html.
Chapter 3 of the text explores the health economics decision of whether or not a healthcare system should adopt the practice of prophylactically administering an antibiotic to women who are about to give birth by Cesarean section. Two versions of a decision model are presented: a deterministic model and a stochastic model. The deterministic model lays out the decision framework and calculates the expected costs of the two options: administering the antibiotic or not. The stochastic model adds uncertainty to the decision framework by assigning probability distributions to the parameters in the deterministic model. The stochastic model is then used to conduct a cost-effectiveness analysis of the two options.
Note that cost effectiveness decision frameworks in healthcare are typically presented as a Cost-Utility analysis, conducted from a systems perspective that measures utility of medical effects in Quality Adjusted Life Years using ICERs (Incremental Cost Effectiveness Ratio) = (Additional Costs) / (Additional Benefits).
The final decision often comes down to examining the equation: Net Benefit = Uλ – C where:
- U = Lifetime utility of intervention (usually measured in QALYS)
- λ = money value attributed to a unit of health gain: i.e., the exchange rate
- C = lifetime cost
The following figure depicts the decision as a tree. The decision of whether or not to prescribe prophylactic antibiotics is represented by the square node labelled Px. The two branches from this node represent the treatment and control groups. These groups are represented as circles because they are “chance nodes”; people in those groups may or may not go on to get a wound infection. The probability of infection is different for each group (p1 and p2). The four leaves of the tree are the four possible outcomes: getting a wound infection or not in either the treatment or the control. We can compute a cost for each outcome. For example, the most expensive outcome is for patients who received the antibiotic and got an infection anyway because they have both the costs of providing the prophylaxis and the costs of treating the wound. Though not shown on the diagram, we can also compute a utility (in QALYs) for each outcome (e.g., time a patient spends in the hospital is not as valuable as time spent otherwise). The expected cost and utility of each main branch are a probability-weighted average of the costs and utilities of the leaves of that branch.
Deterministic Decision Analysis
In this section, we will present a deterministic model of the problem that can be described by means of the following directed graph. Variables in blue on the right side if the diagram represent cost data (see the table in the section “Resource Use and Cost Data” below). The costs of the four leaves of our decision tree are computed from these values. A patient with an infected wound will have a different cost per day of care than a patient without a wound infection and presumably a different length of stay in the hospital. In addition, the treatment group has costs associated with buying and administering multiple doses of the prophylactic antibiotic.
The variables nc1 and rc1 are the total number of Caesarean sections and wound infections in the health system making the decision. Since they do not yet use prophylactic antibiotics for the operation, these values let us estimate the prevalence of wound infections in the untreated group (p1).
The green variables in the upper left (a through d) represent data from a clinical trial (see the table in the “Clinical Data” section below). Since the prevalence of wound infection in the hospitals where the study was done is different from that in the target health system, we need to compute the relative risk for treated vs. control patients in the study. This is just the proportion of patients getting an infection in the treatment group, divided by the proportion in the control group:
%7D%7Bc/(c+d)%7D)
where:
- a is the number of women who received antibiotics and subsequently developed infections
- b is the number of women who received antibiotics who did not develop infections
- c is the number of women who received the placebo and subsequently became infected
- d is the number of women who received the placebo who did not become infected
The deterministic model presented in this section may be found on page 56 in section 3.4 of the text.
The deterministic model simply computes the products and sums, and displays the results. There is nothing Bayesian here yet; this model just does some simple math in JAGS to lay the foundation for subsequent probabilistic refinements.
The Clinical Data
The data for the effectiveness of prophylactic antibiotics to reduce infection in women delivering by cesarean section comes from a study by Bibi et al., (1994).
Preparations for running the model.
This section of code loads the necessary libraries and contains a couple of helper functions that are used later in the code. The first function, coda_sample_2_df, extracts data from a coda sample into a dataframe. The second function, plot_densities_from_coda_df, makes comparative density plots for the given variables.
Show the code
# libraries
library(coda, verbose = FALSE)
library(rjags, verbose = FALSE) # coda.samples
library(jagsUI, verbose = FALSE) # wrapper for rjags
library(mcmcplots) # caterplot
library(tidyverse, verbose = FALSE)
library(gridExtra)
# A couple of helper functions
#| message: FALSE
#| warning: FALSE
#| echo: FALSE
# functions
coda_sample_2_df <- function(my_coda_sample) {
# Extract data from a coda sample into a dataframe
seq_along(my_coda_sample) %>%
lapply(function(chain) {
df <- my_coda_sample[[chain]] %>%
as.matrix() %>%
as.data.frame()
df["chain"] <- chain
df
}) %>%
bind_rows() %>%
mutate(chain = factor(chain))
}
plot_densities_from_coda_df <- function(vars, coda_df) {
# Make comparative density plots for the given variables
coda_df %>%
select(all_of(vars)) %>%
pivot_longer(cols = vars) %>%
ggplot(aes(x = value, col = name, fill = name)) +
geom_density(alpha = 0.6)
}
Show the code
effectiveness_df <- data.frame(
Description = c("Prophylactic antibiotics", "Placebo"),
Infection = c(4, 28),
No_infection = c(129, 108),
row.names = c("Treatment", "Control")
) %>% mutate(Total = Infection + No_infection)
effectiveness_df
Description Infection No_infection Total
Treatment Prophylactic antibiotics 4 129 133
Control Placebo 28 108 136
Resource Use and Cost Data
The resource use and cost data comes from Table 3.3, p55. Data is from Mugford et al., 1989, except for the cost of administering antibiotic which was estimated by Welton et al.
Show the code
cost_df <- data.frame(
Parameter = c(
"Length of stay: infection",
"Length of stay: no infection",
"Cost per day: infection",
"Cost per day: no infection",
"Cost per dose: cephaslosporin",
"Cost per administering antibiotic",
"Doses administered",
"Total number of Caesarians",
"Number infections: no antibiotics"
),
Estimate = c(8.80, 6.70, 163.03, 107.26, 5.67, 7.00, 3, 486, 41),
Units = c("days", "days", "£", "£", "£", "£", "count", "count", "count"),
Variable_name = c("loswd", "losnwd", "cstwd", "cstnwd", "cstPx", "cstadmin", "dose", "nc1", "rc1")
)
cost_df
Parameter Estimate Units Variable_name
1 Length of stay: infection 8.80 days loswd
2 Length of stay: no infection 6.70 days losnwd
3 Cost per day: infection 163.03 £ cstwd
4 Cost per day: no infection 107.26 £ cstnwd
5 Cost per dose: cephaslosporin 5.67 £ cstPx
6 Cost per administering antibiotic 7.00 £ cstadmin
7 Doses administered 3.00 count dose
8 Total number of Caesarians 486.00 count nc1
9 Number infections: no antibiotics 41.00 count rc1
The Deterministic Model
The following code builds a model in the BUGS language for specifying probability distributions that will be evaluated by the rjags interface to the JAGS MCMC engine.
Show the code
model_code <- "
model {
cost_nwdpx <- losnwd * cstnwd + dose * (cstPx + cstadmin) # Cost (No infection/Px)
cost_wdpx <- loswd * cstwd + dose * (cstPx + cstadmin) # Cost (Infection/Px)
cost_nwd <- losnwd * cstnwd # Cost (No infection/no Px)
cost_wd <- loswd * cstwd # Cost (Infection/no Px)
RR <- (a/(a+b))/(c/(c+d)) # Relative risk using
# data from table 3.2
p1 <- rc1/nc1
p2 <- RR * p1
costtrt <- ((1-p2) * cost_nwdpx) + p2 * cost_wdpx # Total cost (payoff) Px
costctl <- ((1-p1) * cost_nwd) + p1 * cost_wd # Total cost (payoff) No Px
}"
cost_data <- with(cost_df, setNames(as.list(Estimate), nm = Variable_name))
effectiveness_data <- list(
a = effectiveness_df["Treatment", "Infection"],
b = effectiveness_df["Treatment", "No_infection"],
c = effectiveness_df["Control", "Infection"],
d = effectiveness_df["Control", "No_infection"]
)
deterministic_data <- append(cost_data, effectiveness_data)
This next block of code executes the model and displays the results. Notice that the expected costs of the treatment and control branches match those in the text on p56.
We expect the rate of infection in the treated group (p2) to be lower than in the untreated group (p1) and the overall cost of treatment to be lower than the cost of not treating.
Show the code
results
JAGS output for model '5', generated by jagsUI.
Estimates based on 1 chains of 1 iterations,
adaptation = 100 iterations (sufficient),
burn-in = 0 iterations and thin rate = 1,
yielding 1 total samples from the joint posterior.
MCMC ran for 0.004 minutes at time 2025-05-19 10:23:22.849254.
mean sd 2.5% 50% 97.5% overlap0 f
p1 0.084 NA 0.084 0.084 0.084 FALSE 1
p2 0.012 NA 0.012 0.012 0.012 FALSE 1
costtrt 765.476 NA 765.476 765.476 765.476 FALSE 1
costctl 779.047 NA 779.047 779.047 779.047 FALSE 1
overlap0 checks if 0 falls in the parameter's 95% credible interval.
f is the proportion of the posterior with the same sign as the mean;
i.e., our confidence that the parameter is positive or negative.
Stochastic Decision Analysis
Having established a deterministic baseline model, the next step is to account for the uncertainty in the modeling assumptions by developing a Bayesian stochastic model that uses probability distributions rather than point estimates to represent uncertain parameters.
This diagram shows the relationships between variables in the stochastic model. The red nodes are stochastic variables, as is the orange outcome node (differential cost). These depend on parameterized random distributions. For example, the relative risk (RR) and prevalence of infection in the target hospital (p1) are no longer simple point estimates as they were in the deterministic model; instead, they are modeled by statistical distributions parameterized by the observed data (see equations below). The cost of administering the drug is now randomly taken from a range, and the length of stay estimates are described as Normal distributions with a given mean and precision. Variables computed from these stochastic variables are stochastic as well, so ‘p2’, ‘cst.trt’, and ‘cst.ctl’ are also red, and the outcome of all this ‘diff.cost’ is, in turn, estimated stochastically.
Model with Distributions
The next step towards working our way up to a full Bayesian model is to assign probability distributions to the clinical variables for which there is uncertainty. These are:
Distribution of ln(Relative Risk):
%5Csim%20N(%5Ctheta,%20prec))
%7D%7Bc/(c+d)%7D%5Cright)%20)
%20-%201%20/%20(a%20+%20b)%20%20+%20(1/c)%20-%201(/(c%20+%20d)%7D%20%5C%5C)
where as above:
- a is the number of women who received antibiotics and subsequently developed infections
- b is the number of women who received antibiotics who did not develop infections
- c is the number of women who received the placebo and subsequently became infected
- d is the number of women who received the placebo who did not become infected
Distribution of Prob[infection | no antibiotic]
)
where:
- rc1 = Number infections: no antibiotics
- nc1 = Total number of Carnelians
is the distribution for the probability of an infection with the antibiotic.
)
is the distribution of the length of a hospital stay with an infected wound where:
%5E2%7D)
)
is the distribution for length of hospital stay without infection, where:
%5E2%7D)
)
is the distribution of the antibiotic dose administered.
((cstPx%20+%20cstadmin)3%20+%20(losnwd)(cstnwd))%20+%20p_2((cstPx%20+%20cstadmin)3%20+%20(loswd)(cstwd)))
is the total cost for the antibiotic arm of the decision tree.
(losnwd)(cstnwd)%20+%20p_1(loswd)(cstwd))
is the total cost for the no antibiotic arm of the decision tree.
is the differential cost.
Code for the Stochastic Model
This next block of BUGS code implements the stochastic model described above. The code is similar to the deterministic model, but with the addition of the distributions for the parameters that are uncertain.
Show the code
stochastic_model_code <- "
model{
lnRR ~ dnorm(theta, prec) # Distribution for ln(Relative Risk)
theta <- log( (a/(a+b)) / (c/(c+d)) )
prec <- 1/( (1/a) - (1/(a+b)) + (1/c) - (1/(c+d)) )
p1 ~ dbeta(alpha, beta) # Distribution for Prob(Infection/NoPx)
alpha <- rc1
beta <- nc1 - rc1
p2 <- exp(lnRR) * p1 # Distribution for Prob(Infection/Px)
loswd ~ dnorm(mnloswd, precwd) # Distribution for length of stay with infection
precwd <- 1/pow(sdloswd/sqrt(numwd), 2)
losnwd ~ dnorm(mnlosnwd, precnwd) # Distribution for length of stay w/o infection
precnwd <- 1/pow(sdlosnwd/sqrt(numnwd), 2)
cstadmin ~ dunif(4, 10) # Px administration
cst.trt <- (1-p2)*((cstPx + cstadmin)*3 + (losnwd*cstnwd)) + p2*((cstPx + cstadmin)*3 + (loswd*cstwd)) # Total cost (payoff) Px
cst.ctl <- (1-p1)*(losnwd*cstnwd) + p1*(loswd*cstwd) # Total cost (payoff) No Rx
diff.cost <- cst.trt - cst.ctl # Difference in cost
}" %>% textConnection()
Data for the Stochastic Model
This block of R code provides parameter values for the stochastic model. The values are taken from the text, but some of them have been changed to reflect the fact that we are now using a stochastic model. The changes are indicated in the comments.
Show the code
stochastic_data <- list(
rc1 = 41, nc1 = 486, cstwd = 163.03, cstnwd = 107.26,
mnloswd = 8.8, sdloswd = 3.5,
mnlosnwd = 6.7, sdlosnwd = 7.1,
numwd = 41, numnwd = 445, # !!! numwd == rc1; numnwd == (nc1 - rc1)
cstPx = 5.67,
# rt=4, nt=133, rc=28, nc=136, # !!! Unused variables
a = 4, b = 129, c = 28, d = 108
)
# Here I code the changes from the data for the deterministic model to the
# data for the stochastic one, in case that makes the differences easier to see.
# Add these:
# setdiff(names(stochastic_data), names(deterministic_data))
# "mnloswd" "sdloswd" "mnlosnwd" "sdlosnwd"
# "numwd" == rc1;
# "numnwd" == (nc1 - rc1)
#
# Remove or rename these:
# setdiff(names(deterministic_data), names(stochastic_data))
# "loswd" -> "mnloswd"
# "losnwd" -> "mnlosnwd"
# "cstadmin" hard-coded uniform distribution
# "dose": hard-coded !!!
# uncertainty <- c(sdloswd=3.5, sdlosnwd=7.1) # standard deviations from Table 3.4
# rename_me <- c("loswd"="mnloswd", "losnwd"="mnlosnwd")
#
# stochastic_data <- append(deterministic_data, uncertainty)
# for (old_name in names(rename_me)){
# new_name <- rename_me[[old_name]]
# names(stochastic_data)[names(stochastic_data) == old_name] <- new_name
# }
# stochastic_data["cstadmin"] <- NULL # now coded as uniform random
# stochastic_data["dose"] <- NULL # the 3 is hard coded now
#
# stochastic_data["numwd"] <- stochastic_data["rc1"];
# stochastic_data["numnwd"] <- with(stochastic_data, (nc1 - rc1))
parameters_to_save <- c("cst.trt", "cst.ctl", "diff.cost")
Code to Run the Stochastic Model
This code runs the stochastic model. The code is similar to the deterministic model, but with the addition of the distributions for the parameters that are uncertain.
Show the code
stochastic_results <- jags(
data = stochastic_data,
parameters.to.save = parameters_to_save,
model.file = stochastic_model_code,
n.chains = 4,
n.adapt = 100,
n.iter = 50000,
n.burnin = 20000,
verbose = FALSE
)
Model Results
The following table shows a summary of the model results. These include the mean, standard deviations, and quantiles for the posterior distributions of the cost of treatment and control, the differential cost and the relative risk, along with information on some diagnostic statistics. The results are similar to those in the text on p 66.
Show the code
stochastic_results_summary <- summary(stochastic_results) |>
as.data.frame() |>
select(!c(3, 4))
round(stochastic_results_summary, 3)
mean sd 50% 75% 97.5% Rhat n.eff overlap0 f
cst.trt 766.945 36.390 767.034 791.444 838.007 1 120000 0 1.000
cst.ctl 779.238 35.095 779.353 802.934 848.135 1 120000 0 1.000
diff.cost -12.292 12.803 -11.779 -3.518 11.234 1 120000 1 0.835
The diagnostic statistics are interpreted as follows:
- Rhat, the Gelman-Rubin Statistic, is a diagnostic that compares the variance within chains to the variance between chains. Values close to 1 (typically less than 1.1) indicate that the chains have converged.
- n.eff: provides an estimate of how many independent samples the samples in the chain are equivalent to. A higher number suggests more reliable estimates.
- overlap0 = 0 indicates that the 95% credible interval does not include 0, suggesting a statistically significant effect.
- f is the proportion of the posterior with the same sign as the mean.
Diagnostic Plots
In this section, we will plot the results of the MCMC sampling. The first plot in each row shows the trace history for the parameters of interest, and the second shows the density. Multiple colors indicate the four independent MCMC chains. Since they are on top of each other and mixing well, we can see that the chains have converged to a common distribution.
Show the code
plot(stochastic_results)
Prior Distributions
Here, we plot the Prior distributions that were described above.
Show the code
## Prior Distribution for log Relative Risk
a <- stochastic_data$a
b <- stochastic_data$b
c <- stochastic_data$c
d <- stochastic_data$d
theta <- log((a / (a + b)) / (c / (c + d)))
# theta
prec <- 1 / ((1 / a) - (1 / (a + b)) + (1 / c) - (1 / (c + d)))
# prec
var <- 1 / prec
# var
# Generate a sequence of x values
x <- seq(-15, 15, length.out = 100)
# Compute the corresponding beta density values
lnRR <- dnorm(x, theta, prec)
# Create a data frame for ggplot2
norm_data <- data.frame(x = x, y = lnRR)
# Plot the beta distribution
q1 <- ggplot(norm_data, aes(x = x, y = lnRR)) +
geom_line(color = "brown", linewidth = 1) +
geom_area(fill = "brown", alpha = 0.2) +
labs(
title = "Prior: Log Relative Risk",
x = "Log(RR)",
y = "Density"
)
### Prior Distribution for p1
# Define the beta distribution parameters
alpha <- stochastic_data$rc1
beta <- stochastic_data$nc1 - stochastic_data$rc1
# Generate a sequence of x values
x <- seq(0, .2, length.out = 100)
# Compute the corresponding beta density values
y <- dbeta(x, alpha, beta)
# Create a data frame for ggplot2
beta_data <- data.frame(x = x, y = y)
# Plot the beta distribution
q2 <- ggplot(beta_data, aes(x = x, y = y)) +
geom_line(color = "blue", linewidth = 1) +
geom_area(fill = "blue", alpha = 0.2) +
labs(
title = "Prior: Infection | no antibiotic",
x = "p1",
y = "Density"
)
### Prior distribution for length of stay, wound infection
# Define the normal distribution parameters
mean <- stochastic_data$mnloswd
sd <- stochastic_data$sdlosnwd / sqrt(stochastic_data$numwd)
# Generate a sequence of x values
x <- seq(5, 15, length.out = 100)
# Compute the corresponding beta density values
y <- dnorm(x, mean, sd)
# Create a data frame for ggplot2
norm_data <- data.frame(x = x, y = y)
# Plot the beta distribution
q3 <- ggplot(norm_data, aes(x = x, y = y)) +
geom_line(color = "red", linewidth = 1) +
geom_area(fill = "red", alpha = 0.2) +
xlim(5, 13) +
labs(
title = "Prior: length of stay | infection",
x = "Length of Stay (days) ",
y = "Density"
)
## Prior distribution for length of stay, no wound infection
# Define the normal distribution parameters
mean <- stochastic_data$mnlosnwd
sd <- stochastic_data$sdlosnwd / sqrt(stochastic_data$numnwd)
# Generate a sequence of x values
x <- seq(5, 10, length.out = 100)
# Compute the corresponding beta density values
y <- dnorm(x, mean, sd)
# Create a data frame for ggplot2
norm_data <- data.frame(x = x, y = y)
# Plot the beta distribution
q4 <- ggplot(norm_data, aes(x = x, y = y)) +
geom_line(color = "darkgreen", linewidth = 1) +
geom_area(fill = "green", alpha = 0.2) +
xlim(5, 13) +
labs(
title = "Prior: length of stay | no infection",
x = "Length of Stay (days)",
y = "Density"
)
grid.arrange(q1, q3, q2, q4)
Posterior Distributions
Here we re-sample an existing model and plot densities of several variables using the R functions defined at the beginning of this document.
Show the code
my_vars <- c("p1", "p2", "cst.trt", "cst.ctl", "diff.cost")
coda_df <- coda.samples(stochastic_results$model,
variable.names = my_vars,
n.iter = 10000
) %>%
coda_sample_2_df() %>%
sample_frac(1L) # randomly scramble chain order
Plot the posterior densities for the cost of treatment
Show the code
w1 <- plot_densities_from_coda_df(c("cst.trt", "cst.ctl"), coda_df) + ggtitle("Posterior Densities for Cost of Treatment")
w1
Plot the posterior density for Differential Cost
Show the code
w2 <- coda_df %>% ggplot(aes(x = diff.cost)) +
geom_density(fill = "blue", alpha = 0.5) +
ggtitle("Posterior Density for Differential Cost")
w2
Plot the densities for P1 and P2
Show the code
w3 <- plot_densities_from_coda_df(c("p1", "p2"), coda_df) + ggtitle("Posterior Densities for P1 and P2")
w3
The Economic Model
This next section of code integrates the computations for cost-effectiveness into the model model presented above. You can find this code on page 65 of the text.
Show the code
economic_model_code <- "
model{
lnRR ~ dnorm(theta, prec) # Distribution for ln(Relative Risk)
theta <- log( (a/(a+b)) / (c/(c+d)) )
prec <- 1/( (1/a) - (1/(a+b)) + (1/c) - (1/(c+d)) )
p1 ~ dbeta(alpha, beta) # Distribution for Prob(Infection/NoPx)
alpha <- rc1
beta <- nc1 - rc1
p2 <- exp(lnRR) * p1 # Distribution for Prob(Infection/Px)
loswd ~ dnorm(mnloswd, precwd) # Distribution for length of stay with infection
precwd <- 1/pow(sdloswd/sqrt(numwd), 2)
losnwd ~ dnorm(mnlosnwd, precnwd) # Distribution for length of stay w/o infection
precnwd <- 1/pow(sdlosnwd/sqrt(numnwd), 2)
cstadmin ~ dunif(4, 10) # Px administration
cst.trt <- (1-p2)*((cstPx + cstadmin)*3 + (losnwd*cstnwd)) + p2*((cstPx + cstadmin)*3 + (loswd*cstwd)) # Total cost (payoff) Px
cst.ctl <- (1-p1)*(losnwd*cstnwd) + p1*(loswd*cstwd) # Total cost (payoff) No Rx
diff.cost <- cst.trt - cst.ctl # Difference in cost
# Economic evaluation code from pp 65-66
totQALYs.wd <- ((QALYwd/365)*loswd) + ((Fullhealth/365)*(fllwupdays-loswd)) # QALYs (infection)
totQALYs.nwd <- ((QALYnwd/365)*losnwd) + ((Fullhealth/365)*(fllwupdays-losnwd)) # QALYs (No infection)
QALYs.trt <- (1-p2) * totQALYs.nwd + p2 * totQALYs.wd # QALYs (Px)
QALYs.ctl <- (1-p1) * totQALYs.nwd + p1 * totQALYs.wd # QALYs (no px)
diff.QALYs <- (QALYs.trt - QALYs.ctl) # Difference in QALYs
for (k in 1:M)
{
lambda[k] <- (k-1) * 2000
INB[k] <- lambda[k] * diff.QALYs - diff.cost # !!! not `delta.Qalys` or `delta.cost`
ProbCE[k] <- step(INB[k])
}
}
" %>% textConnection()
# economic_data <- list(
# rc1=41, nc1=486, cstwd=163.03, cstnwd=107.26,
# mnloswd=8.8, sdloswd=3.5,
# mnlosnwd=6.7, sdlosnwd=7.1,
# numwd=41, numnwd=445,
# cstPx=5.67,
# a=4, b=129, c=28, d=108,
# # additional data for economic evaluation
# M=21, QALYwd=0.68, QALYnwd=0.88, Fullhealth=1, fllwupdays=20
# )
economic_data <- append(
stochastic_data,
# additional data for economic evaluation
list(M = 21, QALYwd = 0.68, QALYnwd = 0.88, Fullhealth = 1, fllwupdays = 20)
)
parameters_to_save <- c( # "ProbCE",
"QALYs.trt", "QALYs.ctl", "diff.QALYs"
)
Run the Economic Model
Show the code
economic_results <- jags(
data = economic_data,
parameters.to.save = parameters_to_save,
model.file = economic_model_code,
n.chains = 1,
n.adapt = 100,
n.iter = 50000,
n.burnin = 20000,
verbose = FALSE
)
And now, the economic results
Show the code
economic_results_summary <- summary(economic_results) %>% as.data.frame()
round(economic_results_summary, 5) |> select(!c(3, 4))
mean sd 50% 75% 97.5% Rhat n.eff overlap0 f
QALYs.trt 0.05251 0.00012 0.05252 0.05259 0.05274 NA 1 0 1.00000
QALYs.ctl 0.05213 0.00013 0.05213 0.05221 0.05238 NA 1 0 1.00000
diff.QALYs 0.00039 0.00008 0.00038 0.00044 0.00055 NA 1 0 0.99997
The Economic Decision: Cost-Effectiveness
The effectiveness of a healthcare intervention is often expressed QALYS, Quality-Adjusted Life Year, which combine mortality and morbidity from treatment outcomes into a single value. By putting diverse medical outcomes on the same scale, it makes it possible to allocate funds among very different programs (cancer care, prenatal care, vaccine programs, etc.).
The following figure plots the results of many simulations of the economic model on the cost-effectiveness plane. This is the major tool for facilitating the decision as to whether to adopt the practice of prophylactically administering antibiotics. The x-axis represents the incremental QALYs and the y-axis represents the incremental costs. The dashed lines represent different values for 
the cut-off cost.
Show the code
XLIM <- c(0, 7e-4) # use same limits as chapter 3
YLIM <- c(-75, 50)
LAMBDAS <- 10000 * (1:3)
CE_1 <- ce_df %>%
ggplot(aes(x = diff.QALYs, y = diff.cost, col = chain)) +
geom_point(size = 0.5, alpha = 0.1) +
geom_abline(slope = 0, intercept = 0, linetype = "solid") +
geom_abline(slope = LAMBDAS, intercept = 0, linetype = "dotted") +
theme(legend.position = "none") +
labs(x = "Incremental QALYs", y = "Incremental costs", title = "Cost-effectiveness plane for Caesarian section example") +
coord_cartesian(xlim = XLIM, ylim = YLIM) +
annotate("text",
x = XLIM[2], y = LAMBDAS * XLIM[2],
size = 3, hjust = 0.8, vjust = -0.25,
label = paste("lambda ==", LAMBDAS), parse = TRUE
)
CE_1
Our final plot shows the expected incremental net benefit and 95% credible interval for a range of lambda values. The x-axis represents the value of lambda and the y-axis represents the expected incremental net benefit. The dashed lines represent the 2.5% and 97.5% quantiles of the posterior distribution of the incremental net benefit.
Show the code
inb_df_long %>%
group_by(metric, lambda) %>%
filter(metric == "INB") %>%
summarize(
mean_INB = mean(value),
lo_end = quantile(value, probs = 0.025),
hi_end = quantile(value, probs = 0.975)
) %>%
ggplot(aes(x = lambda, group = metric)) +
geom_line(aes(y = mean_INB)) +
geom_line(aes(y = lo_end), linetype = "dashed") +
geom_line(aes(y = hi_end), linetype = "dashed")
Continue reading:
Cost-Effectiveness Analysis
