4. Annex on non-stunning interventions

Author

Sagar Shah, Rethink Priorities

Published

March 18, 2024

About this file

This file contains cost-effectiveness analysis of non-slaughter intervention contained within Annex 2 of the report.

Generic prep

Code
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

Code
# Load functions and assumptions
source("a_functions.R")
source("b_fish_assumptions.R") 
source("c_model_assumptions.R") 

Modelling assumptions

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
Code
# Share of fish life affected by welfare intervention
  fishlifeshare_values <- c(0.01,0.02,0.05,0.10,0.15,0.25,0.35,0.50)
  
# Moral value points
  mv_table_values <- c(0.01,0.05,0.1,0.25,0.5,0.75,1,5,10,25)
  
# Moral value points
  wr_table_values <- c(0.001,0.01,0.025,0.5,0.1,0.15,0.35,0.5,0.75,0.9,1)

Calculations

Helper functions

Write helper functions to call object/vectors

Code
consumption <- function(country,species)  obtain(species,country,"tons")
mshare <- function(country) obtain("mshare_welfare",country)
psuccess <- function() get("psuccess_welfare")
weight <- function(species) obtain("weight",species)
lifexp <- function(species) obtain("lifexp",species)
cost <- function() get("cost_welfare")
stun_share <- function(species) obtain(species,"stunned")

Number of fish affected/$ and associated lifespan

Core calculation functions

Calculate number of fish affected per dollar spent

Code
# Number of fish affected per dollar
no_fish_affected <- function(country,species) {
consumption(country,species)*1E3/
  weight(species)*
  mshare(country)*
  psuccess()*
  stun_share(species)*
  years_credit*
  implementation_discount*
  fish_grocery
}

Life span of affected fish

Code
# Fish years of affected fish}
lifespan_affected_fish <- function(country,species) {
no_fish_affected(country,species)*lifexp(species)/12
}

Execution functions

Produce a table that produces every possible combination of country and species, and a list of functions over which to calculate results.

Code
no_of_countries <- length(country_list)
country_species <- expand.grid(country=country_list,species=species_list) 
measure_list <- c("lifespan_affected_fish","no_fish_affected")

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.

Code
output_per_dollar <- function(fish_function) {

      pmap(                           # Allows repeated calcs across different input vectors
          country_species,                # Input vector 1: countries to run function over
          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

}

Run calculations

I then produce tables with montecarlo simulation results for each measure in both wide and long format.

Code
impact_per_dollar_wide <- map(measure_list,output_per_dollar) %>%                            
                              set_names(measure_list) %>%
                              as.data.frame() 
                      
impact_per_dollar_long <-  impact_per_dollar_wide %>% 
                            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.

Code
fypd <- output_per_dollar("lifespan_affected_fish")

$/DALY results

Core calculation functions

Code
results_dollar_per_daly <-
   1/(
    fypd*
    salmon_wr*
    interventionlifesshare*
    fish_welfarerange_impact/
    DALY_share_of_welfare_range) %>%
    tibble() %>%
    set_names("dollars_per_daly")

Share of simulations beating $$/DALY benchmark

Core calculation functions

Code
#Moral value approach
share_mv <- function(xvar,fishlifeshare,bar) {
  mean(bar>(1/(fypd*fishlifeshare*xvar)))
} 

# Welfare range approach
share_wr <- function(xvar,fishlifeshare,bar) {
  mean(bar>
      (1/
        (xvar*
        fypd*
        fishlifeshare*
        salmon_wr/
        DALY_share_of_welfare_range
        )))
}

Calculations

Define xvar (moral value or welfare range) values for calculating share of simulations beating a particular $/DALY bar and create combinations of xvar values, bar values and, lifeshare values over which to calculate results.

Code
 chart_values<- 10^seq(-4,4,0.01)
  xvar_values <- c(chart_values,mv_table_values,wr_table_values)
  
share_sims_calculation_grid <-   expand_grid(
                                    xvar=xvar_values,
                                    bar=bar_values,
                                    fishlifeshare=fishlifeshare_values)

Calculate results

Code
share_sims_results_table <- 
# Generate one row for unique combination of xvar_values and scenario lists
  share_sims_calculation_grid %>%
# Core calcualtions using share_mv and share_wr calculation functions
    rowwise() %>%
    mutate(
      moral_value=share_mv(xvar,fishlifeshare,bar),
      welfare_range=share_wr(xvar,fishlifeshare,bar)
    ) 

Format results table

Code
share_sims_results_table %<>%
# Convert pivot into longer format to make results easier to plot
  pivot_longer(
    cols=c(moral_value,welfare_range),
    names_to="approach",
    values_to="share_of_simulations") %>%
# Change descriptions of variables to make results easier to interpert
  mutate(
    approach=str_to_sentence(str_replace_all(approach,"_"," ")),
    bar_factor=factor(bar,levels=bar_values,labels=c("50","1K","70K")),
    fishlifeshare_factor=factor(fishlifeshare,levels=rev(fishlifeshare_values),labels=paste0(rev(fishlifeshare_values)*100,"%")))

Results tables/charts production functions

Impact per dollar density plots

Code
# Code to create impact per dollar density plot 

impact_density_plot <- function(desired_filter,named_description) { 
  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)"))
}

Share of simulations

Chart labels

Code
xlab_wr="
    Welfare gain from intervention as a % of total fish welfare range"

xlab_mv="Moral value of improving a fish life year via intervention
relative to averting a human DALY"

Chart function

Code
share_sims_chart <- function (
                        chart_approach,
                        xvar_lower_lim=0,
                        xvar_upper_lim=max(xvar_values),
                        xlab_des,
                        bar_value=50
                        ) {
  
share_sims_results_table %>%
    filter(
    approach==chart_approach,
    xvar<=xvar_upper_lim,
    xvar>=xvar_lower_lim,
    bar==bar_value
    ) %>%
  ggplot(
    aes(
      x=xvar,
      y=share_of_simulations,
      color=fishlifeshare_factor))+
  geom_line()+
  scale_color_brewer(palette = "Dark2") +
  theme_light()+
  labs(
    subtitle = str_to_title(paste0(chart_approach," approach")),
    title=paste0("Share of simulations where intervention beats ",bar_value,"$/DALY benchmark"),
    y="Share of simulations",
    x=xlab_des,
    linetype="$/DALY Benchmark",
    color="Share of fish \nlifespan affected \nby intervention",
    caption= "Note log scale on x-axis")
}

Table function

Code
share_beating_bar_table <- function(values,chart_approach,bar_value) {

share_sims_results_table %>%
  filter(
      xvar %in% values,
      fishlifeshare %in% fishlifeshare_values,
      bar==bar_value,
      approach==chart_approach) %>%
      select(xvar,share_of_simulations,fishlifeshare_factor) %>%
      distinct() %>%
      pivot_wider(id_cols=xvar,names_from=fishlifeshare_factor,values_from=share_of_simulations) %>%
      arrange(xvar) %>%
      gt() %>%
      tab_header(title = paste0("Share of simulations beating $",bar_value,"/DALY bar")) %>%
      tab_spanner(label="Share of fish lifespan affected by intervention",columns=paste0(fishlifeshare_values*100,"%"))
}

wr_bar_table <- function(bar_value) {
  share_beating_bar_table(wr_table_values,"Welfare range",bar_value) %>%
  fmt_percent(decimals=1) %>%  
  fmt_percent(columns=xvar,decimals=0) %>%
  cols_label(xvar=xlab_wr) 
}

mv_bar_table <- function(bar_value) {
  share_beating_bar_table(mv_table_values,"Moral value",bar_value) %>%
  fmt_percent(decimals=1) %>%  
  fmt_number(columns=xvar,drop_trailing_zeros=TRUE) %>%
  cols_label(xvar=xlab_mv) 
}

Outputs

Impact per dollar

Density plot: number of fish affected

Code
impact_density_plot("no_fish_affected","Number of fish affected per dollar")

Density plot: lifespan of fish affected

Code
impact_density_plot("lifespan_affected_fish","Lifespan of affected fish per dollar")