This toy model is based on this tutorial: https://cran.r-project.org/web/packages/hesim/vignettes/cea.html#decision-analysis
Is offering more vegetarian options at lunchtime a cost-effective way to get college students to switch from meat to plant-based meat? Let’s imagine a philanthropist wants to know the cost of switching one meal containing meat to a vegetarian meal. The long-term welfare impact of switching one meal on farmed animals is unclear, given that changes in meat consumption in the cafeteria may not have a sustained effect on diet beyond the cafeteria. We set aside this issue here, settling for a proxy.
We’ll need the following packages:
library(hesim)
library(tidyverse)
library(data.table)
Based on Garnett et al. 2019 (https://www.pnas.org/content/116/42/20923.short), we evaluate the effect of increasing the percentage of vegetarian food options at a college cafeteria from 25% to 50% of total meal choices. Their intervention increased vegetarian meals/decreased meat meals by 7.8% points (from 19.1% vegetarian to 26.9% vegetarian). Although this is not what Garnett did, let’s imagine that the intervention increases vegetarian options to 50% by swapping out some meat choices for plant-based imitations of those meats.
Let’s assume plant-based meat costs 3 USD while meat and non-plant-based meat vegetarian meals (e.g. tofu) costs 0.75 USD: https://www.cnbc.com/2021/08/25/impossible-foods-beyond-meat-battle-price-parity-with-real-meat.html
control_mean<-100000*0.75 #costs of 25% vegetarian options, no plant-based meat
experimental_mean<-(75000*0.75) + (25000*3.00)# costs of 50% vegetarian option, including plant-based meat
Of course, there is uncertainty about the magnitude of both costs and benefits. By sampling from a probability distribution, we can estimate over the distribution of plausible values, rather than putting all of our confidence in just one value for each parameter. Let’s start with costs, assuming a standard deviation of 2,000 USD in the control group and 5,000 USD in the experimental group. Because prices should only take positive values, we use log-normal distributions.
# https://msalganik.wordpress.com/2017/01/21/making-sense-of-the-rlnorm-function-in-r/
#control group (25% vegetarian options)
m_c <- control_mean
s_c <- 2000
location_c <- log(m_c^2 / sqrt(s_c^2 + m_c^2))
shape_c <- sqrt(log(1 + (s_c^2 / m_c^2)))
print(paste("location:", location_c))
## [1] "location: 11.2248879633227"
print(paste("shape:", shape_c))
## [1] "shape: 0.0266619277511154"
#experimental group (50% vegetarian options)
m_e <- experimental_mean
s_e <- 5000
location_e <- log(m_e^2 / sqrt(s_e^2 + m_e^2))
shape_e <- sqrt(log(1 + (s_e^2 / m_e^2)))
print(paste("location:", location_e))
## [1] "location: 11.7841340828918"
print(paste("shape:", shape_e))
## [1] "shape: 0.0380814275479553"
set.seed(131) # for replicability
n_samples <- 10000 # 10,000 random draws
# cost
c <- vector(mode = "list", length = 2)
names(c) <- c("Control", "Experimental")
c[[1]] <- rlnorm(n_samples, location_c, shape_c)
c[[2]] <- rlnorm(n_samples, location_e, shape_e)
# We use normal distributions because negative values are possible
#we'll assume a sd of 4,000 vegetarian meals in each group for simplicity
e <- c
e[[1]] <- rnorm(n_samples, 19100, 4000) #19.1% vegetarian consumption in control group
e[[2]] <- rnorm(n_samples, 26900, 4000) #26.9% vegetarian consumption in experimental group
A table of the costs and benefits of the experimental and control groups is the input for the rest of the analysis. We assign the control group of 25% vegetarian meals “Meal Plan1,” and the experimental group with 50% vegetarian meals due to the inclusion of plant-based meat “Meal Plan2.”
ce <- data.table(sample = rep(seq(n_samples), length(e)),
strategy = rep(paste0("Meal Plan", seq(1, 2)),
each = n_samples),
cost = do.call("c", c), veg_meals = do.call("c", e))
ce
## sample strategy cost veg_meals
## 1: 1 Meal Plan1 73354.27 16992.29
## 2: 2 Meal Plan1 73894.08 20407.30
## 3: 3 Meal Plan1 77042.07 18832.51
## 4: 4 Meal Plan1 75066.76 20965.69
## 5: 5 Meal Plan1 74079.00 23154.48
## ---
## 19996: 9996 Meal Plan2 132703.00 25597.28
## 19997: 9997 Meal Plan2 129974.51 27671.78
## 19998: 9998 Meal Plan2 133778.52 26787.81
## 19999: 9999 Meal Plan2 129547.13 32506.11
## 20000: 10000 Meal Plan2 131546.01 26328.39
Let’s imagine that the most a donor would be willing to pay to avert one meal containing meat is 50 USD. By multiplying willingness to pay by numbers of meat meals averted, we are able to put benefits on the same scale as costs. Thus, we can compute Net Monetary Benefit, which is the monetized benefits of the intervention less its costs. In an analysis where we vary maximum willingness to pay from 0 USD to 50 USD in 5 USD intervals, we find that the experimental group has higher net monetary benefit when willingness to pay is at least 10 USD.
ktop <- 50
cea_out <- cea(ce, k = seq(0, ktop, 5), sample = "sample", strategy = "strategy",
e = "veg_meals", c = "cost")
cea_out
## $summary
## strategy grp e_mean e_lower e_upper c_mean c_lower c_upper
## 1: Meal Plan1 1 19132.69 11245.07 27140.64 75008.59 71191.34 78994.32
## 2: Meal Plan2 1 26915.01 19100.12 34763.73 131218.57 121404.55 141119.51
##
## $mce
## k strategy grp best prob
## 1: 0 Meal Plan1 1 1 1.0000
## 2: 0 Meal Plan2 1 0 0.0000
## 3: 5 Meal Plan1 1 1 0.7264
## 4: 5 Meal Plan2 1 0 0.2736
## 5: 10 Meal Plan1 1 0 0.3468
## 6: 10 Meal Plan2 1 1 0.6532
## 7: 15 Meal Plan1 1 0 0.2395
## 8: 15 Meal Plan2 1 1 0.7605
## 9: 20 Meal Plan1 1 0 0.1897
## 10: 20 Meal Plan2 1 1 0.8103
## 11: 25 Meal Plan1 1 0 0.1654
## 12: 25 Meal Plan2 1 1 0.8346
## 13: 30 Meal Plan1 1 0 0.1493
## 14: 30 Meal Plan2 1 1 0.8507
## 15: 35 Meal Plan1 1 0 0.1393
## 16: 35 Meal Plan2 1 1 0.8607
## 17: 40 Meal Plan1 1 0 0.1314
## 18: 40 Meal Plan2 1 1 0.8686
## 19: 45 Meal Plan1 1 0 0.1261
## 20: 45 Meal Plan2 1 1 0.8739
## 21: 50 Meal Plan1 1 0 0.1220
## 22: 50 Meal Plan2 1 1 0.8780
## k strategy grp best prob
##
## $evpi
## grp k best enmbci enmbpi evpi
## 1: 1 0 Meal Plan1 -75008.59 -75008.59 1.164153e-10
## 2: 1 5 Meal Plan1 20654.88 25546.08 4.891207e+03
## 3: 1 10 Meal Plan2 137931.57 151520.99 1.358941e+04
## 4: 1 15 Meal Plan2 272506.65 284637.37 1.213072e+04
## 5: 1 20 Meal Plan2 407081.72 419206.68 1.212496e+04
## 6: 1 25 Meal Plan2 541656.79 554322.74 1.266595e+04
## 7: 1 30 Meal Plan2 676231.86 689662.74 1.343088e+04
## 8: 1 35 Meal Plan2 810806.93 825113.90 1.430696e+04
## 9: 1 40 Meal Plan2 945382.00 960640.84 1.525883e+04
## 10: 1 45 Meal Plan2 1079957.07 1096214.37 1.625730e+04
## 11: 1 50 Meal Plan2 1214532.15 1231818.78 1.728664e+04
##
## $nmb
## strategy grp k enmb lnmb unmb
## 1: Meal Plan1 1 0 -75008.592 -78994.32 -71191.34
## 2: Meal Plan2 1 0 -131218.570 -141119.51 -121404.55
## 3: Meal Plan1 1 5 20654.876 -19032.71 60656.40
## 4: Meal Plan2 1 5 3356.502 -36580.14 43798.32
## 5: Meal Plan1 1 10 116318.344 37404.94 195715.36
## 6: Meal Plan2 1 10 137931.574 58431.25 217209.25
## 7: Meal Plan1 1 15 211981.812 93822.01 331349.53
## 8: Meal Plan2 1 15 272506.645 153993.77 391030.33
## 9: Meal Plan1 1 20 307645.280 149573.70 466749.53
## 10: Meal Plan2 1 20 407081.717 249067.15 564956.06
## 11: Meal Plan1 1 25 403308.748 205651.05 602481.61
## 12: Meal Plan2 1 25 541656.788 344784.75 738112.03
## 13: Meal Plan1 1 30 498972.216 261915.25 738065.34
## 14: Meal Plan2 1 30 676231.860 440632.03 912246.72
## 15: Meal Plan1 1 35 594635.684 318120.12 873641.59
## 16: Meal Plan2 1 35 810806.931 535952.64 1086164.55
## 17: Meal Plan1 1 40 690299.152 374425.55 1009065.84
## 18: Meal Plan2 1 40 945382.003 631602.03 1260640.88
## 19: Meal Plan1 1 45 785962.620 430538.03 1144803.99
## 20: Meal Plan2 1 45 1079957.074 727235.78 1434500.88
## 21: Meal Plan1 1 50 881626.088 486652.83 1280385.13
## 22: Meal Plan2 1 50 1214532.146 822586.28 1607671.88
## strategy grp k enmb lnmb unmb
##
## attr(,"class")
## [1] "cea"
## attr(,"strategy")
## [1] "strategy"
## attr(,"grp")
## [1] "grp"
We now directly compare the experimental group to the control group. We find that the Incremental Cost-Effectiveness Ratio, or the cost of averting one additional meal containing meat is 7.22 USD. This is just the ratio of the differences in costs and the differences in benefits.
cea_pw_out <- cea_pw(ce, k = seq(0, ktop, 5), comparator = "Meal Plan1",
sample = "sample", strategy = "strategy",
e = "veg_meals", c = "cost")
cea_pw_out
## $summary
## strategy grp ie_mean ie_lower ie_upper ic_mean ic_lower ic_upper
## 1: Meal Plan2 1 7782.321 -3523.513 19033.06 56209.98 45710.33 66935.85
## icer
## 1: 7.222778
##
## $delta
## sample strategy grp ie ic
## 1: 1 Meal Plan2 1 12711.7755 74722.81
## 2: 2 Meal Plan2 1 9363.9379 56011.69
## 3: 3 Meal Plan2 1 -875.8377 48660.72
## 4: 4 Meal Plan2 1 6567.7975 49602.51
## 5: 5 Meal Plan2 1 1414.8035 54472.24
## ---
## 9996: 9996 Meal Plan2 1 8807.3787 58103.94
## 9997: 9997 Meal Plan2 1 6243.5862 51921.54
## 9998: 9998 Meal Plan2 1 4046.4816 57025.26
## 9999: 9999 Meal Plan2 1 17936.4341 56381.36
## 10000: 10000 Meal Plan2 1 5016.9480 58384.69
##
## $ceac
## k strategy grp prob
## 1: 0 Meal Plan2 1 0.0000
## 2: 5 Meal Plan2 1 0.2736
## 3: 10 Meal Plan2 1 0.6532
## 4: 15 Meal Plan2 1 0.7605
## 5: 20 Meal Plan2 1 0.8103
## 6: 25 Meal Plan2 1 0.8346
## 7: 30 Meal Plan2 1 0.8507
## 8: 35 Meal Plan2 1 0.8607
## 9: 40 Meal Plan2 1 0.8686
## 10: 45 Meal Plan2 1 0.8739
## 11: 50 Meal Plan2 1 0.8780
##
## $inmb
## strategy grp k einmb linmb uinmb
## 1: Meal Plan2 1 0 -56209.98 -66935.85 -45710.33
## 2: Meal Plan2 1 5 -17298.37 -75105.79 39448.10
## 3: Meal Plan2 1 10 21613.23 -91822.42 133843.20
## 4: Meal Plan2 1 15 60524.83 -109794.76 229130.72
## 5: Meal Plan2 1 20 99436.44 -127264.76 324693.16
## 6: Meal Plan2 1 25 138348.04 -144158.75 419111.82
## 7: Meal Plan2 1 30 177259.64 -161663.57 514255.72
## 8: Meal Plan2 1 35 216171.25 -179549.70 609094.87
## 9: Meal Plan2 1 40 255082.85 -197641.68 703564.09
## 10: Meal Plan2 1 45 293994.45 -214825.64 798989.53
## 11: Meal Plan2 1 50 332906.06 -232009.59 893897.05
##
## attr(,"class")
## [1] "cea_pw"
## attr(,"strategy")
## [1] "strategy"
## attr(,"grp")
## [1] "grp"
## attr(,"comparator")
## [1] "Meal Plan1"
## attr(,"comparator_pos")
## [1] 1
Let’s make some plots to visualize the results. The cost-effectiveness plane shows in which samples the experimental group was cost-effective, assuming a willingness to pay of 50 USD. Dots below the line are sufficiently cost-effective, while dots above the line are insufficiently cost-effective.
p<-plot_ceplane(cea_pw_out, k = 50)
p + labs(x = "Meat Meals Averted")
Now, let’s plot the probability that the 50% vegetarian plan (“Meal Plan2”) is more cost-effective than the 25% vegetarian group (“Meal Plan1”).
plot_ceac(cea_out) # lines cross each other around $7
plot_ceac(cea_pw_out)#greater than 80% chance pb meat is cost-effective if willing to pay at least $20
Finally, let’s put a dollar value on our uncertainty. The expected value of perfect information is the most money one would pay to completely eliminate uncertainty about whether introducing plant-based meat is sufficiently cost-effective by collecting more data.
plot_evpi(cea_out) # should be willing to pay about 13k to get perfect data if willing to pay $10 per meat meal avoided.