Code
rm(list=ls())
library(tidyverse)
options(dplyr.summarise.inform = FALSE)
library(scales)
library(rlang)
library(gt)
library(RColorBrewer)
library(forcats)
library(magrittr)
library(webshot2)
This conducts cost-effectiveness analysis and produces the main results tables charts used in the Results section of the report.
Open libraries
rm(list=ls())
library(tidyverse)
options(dplyr.summarise.inform = FALSE)
library(scales)
library(rlang)
library(gt)
library(RColorBrewer)
library(forcats)
library(magrittr)
library(webshot2)
Load functions and assumptions
#Useful monte carlo functions
source("a_functions.R")
#Take on board assupmtions for seabream, seabass and trout
source("b_fish_assumptions.R")
#Modelling assumptions
source("c_model_assumptions.R")
Write helper functions to call object/vectors
<- function(country,species) obtain(species,country,"tons")
consumption <- function(scenario,country) obtain("mshare",scenario,country)
mshare <- function(scenario) obtain("psuccess",scenario)
psuccess <- function(species) obtain("weight",species)
weight <- function(species) obtain("slaughter_minutes",species)
slaughter_minutes <- function(species) obtain(species,"stunned")
stun_share <- function(scenario) obtain("cost",scenario) cost
# Number of fish affected per year
<- function(country,species,scenario) {
no_fish_affected_per_year consumption(country,species)*1E3/
weight(species)*
1-stun_share(species))*
(mshare(scenario,country)*
psuccess(scenario)*
*
implementation_discount
fish_grocery
}
# Fish affected across all years of credit
<- function(country,species,scenario) {
no_fish_affected no_fish_affected_per_year(country,species,scenario)*years_credit
}
# Fish minutes affected
<- function(country,species,scenario) {
fish_hours_affected no_fish_affected(country,species,scenario)*slaughter_minutes(species)/60
}
I now produce a table that produces a possible combination of country and species
<- length(country_list)
no_of_countries <- expand.grid(country=country_list,species=species_list) %>%
combo mutate(country_species=str_c(country,"_",species))
I now write a function that:
Works out the total number of fish (or fish minutes) affected in a given scenario, country and species
Sums the result for all species and countries
Divides by the total program costs for a specific scenario
<- function(fish_function,program) {
output_per_dollar
<- map2( # Allows repeated calcs across different input vectors
temp_output$country, # Input vector 1: countries to run function over
combo$species, # Input vector 2: species to run function over
comboget(fish_function), # Calculation function (no of fish or minutes)
scenario=program # Scenario (pilot or scale)
%>%
) as.data.frame() %>% # Convert to data.frame
mutate(
total=rowSums(.))
return(
$total/(no_of_countries*cost(program))
temp_output
) }
Generate table for combinations of output and scenario
<- expand_grid(
results_combination measure=c("fish_hours_affected","no_fish_affected","no_fish_affected_per_year"),
program=c("pilot","scale")) %>%
mutate(description=str_c(program,"___",measure))
I now write a function that works out the number of fish minutes affected and total number of fish affected per dollar across both the pilot and scale scenarios. The results of the montecarlo simulation are in a single table called “impact_per_dollar”.
<- map2( # Function to iterate over vectors
impact_per_dollar $measure, # Vector of different measures
results_combination$program, # Vector of different scenarios
results_combination%>% # Function to repeat
output_per_dollar) as.data.frame()
# Assign names to data.frame
names(impact_per_dollar) <- results_combination$description
# Convert to long format
%<>%
impact_per_dollar pivot_longer(
cols=1:nrow(results_combination),
names_to="description",
values_to="number"
%>%
) separate_wider_delim(description,delim="___",names=c("scenario","measure"))
I first produce estimates of fish years affected per dollar in both scenarios.
<- output_per_dollar("fish_hours_affected","pilot")*hours_to_years
fypd_pilot <- output_per_dollar("fish_hours_affected","scale")*hours_to_years
fypd_scale <- function(scenario) obtain("fypd",scenario) fypd
<- function(scenario) {
wr_dollar_per_daly 1/(
fypd(scenario)*
*
salmon_wr*
duration_share/
fish_welfarerange_impact_stun
DALY_share_of_welfare_range) }
Calculate results and place into data frame
<-
results_dollar_per_daly tibble(
# core calculations
pilot=wr_dollar_per_daly("pilot"),
scale=wr_dollar_per_daly("scale")) %>%
# pivot to long format
pivot_longer(
cols=1:2,
values_to="dollars_per_daly",
names_to="scenario"
%>%
) # reformat for table production
mutate(scenario=str_to_sentence(scenario))
# Code to create impact per dollar density plot
<- function(desired_filter,description) {
impact_density_plot %>%
impact_per_dollar filter(measure==desired_filter) %>%
ggplot(
aes(
x=number,
color=str_to_sentence(scenario)
)+
) geom_density() +
theme_light() +
scale_x_log10(labels = scales::comma_format(drop0trailing = TRUE),n.breaks=8) +
labs(
title=description,
y="Probability density",
x=paste0(description," (log scale)"),
color="Scenario")
}
<- function(desired_filter,description) {
impact_table %>%
impact_per_dollar filter(measure==desired_filter) %>%
mutate(scenario=str_to_sentence(scenario)) %>%
group_by(scenario) %>%
summarise(
Mean=mean(number),
p5=quantile(number,0.05),
Median=median(number),
p95=quantile(number,0.95)
%>%
) pivot_longer(cols=2:5,names_to="Statistic",values_to="Value") %>%
pivot_wider(names_from = "scenario",values_from = "Value") %>%
gt() %>%
fmt_number(decimals = 1) %>%
tab_header(
title = description) %>%
tab_spanner(
label = "Scenario",
columns = c(Pilot,Scale))
}
<- impact_density_plot("no_fish_affected","Number of fish affected per dollar")
fig9_fish_pd <- impact_table("no_fish_affected","Number of fish affected per dollar")
tab8_fish_pd<- impact_table("no_fish_affected_per_year","Number of fish affected per dollar per year")
tab9_fish_pdpy <- impact_density_plot("fish_hours_affected","Fish hours affected per dollar")
fig10_fish_hours_pd <- impact_table("fish_hours_affected","Fish hours affected per dollar") tab10_fish_hours_pd
<-
fig11_dolperdaly_density %>%
results_dollar_per_daly ggplot(
aes(
x=dollars_per_daly,
color=scenario
)+
) geom_density() +
theme_light() +
scale_x_log10(
breaks=10^seq(0,15,1),
labels = scales::label_number(scale_cut = scales::cut_short_scale())) +
labs(
title = "$/DALY Density Plot",
y="Probability density",
x="$/DALY (log-scale)",
color="Scenario")
<-
tab11_dolperdaly_density %>%
results_dollar_per_daly group_by(scenario) %>%
summarise(
Mean=mean(dollars_per_daly),
p1=quantile(dollars_per_daly,0.01),
p5=quantile(dollars_per_daly,0.05),
p10=quantile(dollars_per_daly,0.10),
p25=quantile(dollars_per_daly,0.25),
Median=median(dollars_per_daly),
p75=quantile(dollars_per_daly,0.75),
p90=quantile(dollars_per_daly,0.90),
p95=quantile(dollars_per_daly,0.95),
p99=quantile(dollars_per_daly,0.99)
%>%
) pivot_longer(cols=2:11,names_to="Statistic",values_to="$/DALY") %>%
pivot_wider(names_from = "scenario",values_from = "$/DALY") %>%
gt() %>%
tab_spanner(
label = "Scenario",
columns = c(Pilot,Scale)) %>%
fmt_number(columns=2:3,suffixing=TRUE,n_sigfig =3)
<-
tab11_supplement %>%
results_dollar_per_daly mutate(
`$50/DALY`=dollars_per_daly<50,
`$1K/DALY`=dollars_per_daly<1000,
`$70K/DALY`=dollars_per_daly<70000) %>%
group_by(scenario) %>%
summarise_if(is.logical,sum,na.rm=TRUE) %>%
gt() %>%
tab_header(
title = "Number of simulations beating selected $/DALY benchmarks")
Impact per dollar charts and tables
fig9_fish_pd