Code
rm(list=ls())
library(tidyverse)
library(scales)
library(rlang)
library(gt)
library(RColorBrewer)
library(forcats)
library(magrittr)
library(webshot2)
This file contains cost-effectiveness analysis of non-slaughter intervention contained within Annex 2 of the report.
rm(list=ls())
library(tidyverse)
library(scales)
library(rlang)
library(gt)
library(RColorBrewer)
library(forcats)
library(magrittr)
library(webshot2)
Run code to load useful montecarlo simulation functions and general seabream seabass assumptions
# Load functions and assumptions
source("a_functions.R")
source("b_fish_assumptions.R")
source("c_model_assumptions.R")
Define parameters over which to calculate share of simulation results. Definition of xvar provided in the table below.
Approach | Meaning of ‘xvar’ variable means in context of calculations |
---|---|
Moral Value | Moral value of improving fish welfare for 1 year through intervention relative to averting 1 DALY |
Welfare Range | Welfare gain from intervention - expressed as % of entire fish welfare range (negative to positive |
# Share of fish life affected by welfare intervention
<- c(0.01,0.02,0.05,0.10,0.15,0.25,0.35,0.50)
fishlifeshare_values
# Moral value points
<- c(0.01,0.05,0.1,0.25,0.5,0.75,1,5,10,25)
mv_table_values
# Moral value points
<- c(0.001,0.01,0.025,0.5,0.1,0.15,0.35,0.5,0.75,0.9,1) wr_table_values
Write helper functions to call object/vectors
<- function(country,species) obtain(species,country,"tons")
consumption <- function(country) obtain("mshare_welfare",country)
mshare <- function() get("psuccess_welfare")
psuccess <- function(species) obtain("weight",species)
weight <- function(species) obtain("lifexp",species)
lifexp <- function() get("cost_welfare")
cost <- function(species) obtain(species,"stunned") stun_share
Calculate number of fish affected per dollar spent
# Number of fish affected per dollar
<- function(country,species) {
no_fish_affected consumption(country,species)*1E3/
weight(species)*
mshare(country)*
psuccess()*
stun_share(species)*
*
years_credit*
implementation_discount
fish_grocery }
Life span of affected fish
# Fish years of affected fish}
<- function(country,species) {
lifespan_affected_fish no_fish_affected(country,species)*lifexp(species)/12
}
Produce a table that produces every possible combination of country and species, and a list of functions over which to calculate results.
<- length(country_list)
no_of_countries <- expand.grid(country=country_list,species=species_list)
country_species <- c("lifespan_affected_fish","no_fish_affected") measure_list
I now write a function that calculates number of fish / lifespan affected per country/species pair, sums across all the country/species pairs, and then divides by costs across all countries.
<- function(fish_function) {
output_per_dollar
pmap( # Allows repeated calcs across different input vectors
# Input vector 1: countries to run function over
country_species, get(fish_function) # Calculation function (no of fish or minutes)
%>%
) as.data.frame() %>% # Convert to data.frame
mutate(
total=rowSums(.)) %>% # Work out sum of all country and species combinations
pull(total)/(no_of_countries*cost()) # Divide by costs across all countries
}
I then produce tables with montecarlo simulation results for each measure in both wide and long format.
<- map(measure_list,output_per_dollar) %>%
impact_per_dollar_wide set_names(measure_list) %>%
as.data.frame()
<- impact_per_dollar_wide %>%
impact_per_dollar_long pivot_longer(
cols=1:length(measure_list),
names_to="description",
values_to="number")
Finally I produce a vector called fypd, that calculates the lifespan of affected fish, which is needed in the $/DALY and share of simulation results.
<- output_per_dollar("lifespan_affected_fish") fypd
<-
results_dollar_per_daly 1/(
*
fypd*
salmon_wr*
interventionlifesshare/
fish_welfarerange_impact%>%
DALY_share_of_welfare_range) tibble() %>%
set_names("dollars_per_daly")
# Code to create impact per dollar density plot
<- function(desired_filter,named_description) {
impact_density_plot %>%
impact_per_dollar_long filter(description==desired_filter) %>%
ggplot(
aes(
x=number
)+
) geom_density() +
theme_light() +
scale_x_log10(labels = scales::comma_format(drop0trailing = TRUE),n.breaks=8) +
labs(
title=named_description,
y="Probability density",
x=paste0(named_description," (log scale)"))
}
Density plot: number of fish affected
impact_density_plot("no_fish_affected","Number of fish affected per dollar")
Density plot: lifespan of fish affected
impact_density_plot("lifespan_affected_fish","Lifespan of affected fish per dollar")