publicising_species_discoveries.Rmd

---
title: 'Data analysis and graphics for: The costs and benefits of publicising species
  discoveries.'
author: "GE Ryan"
output:
  html_document:
    df_print: paged
  html_notebook: default
  #md_document
---


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```



```{r packages}
library(dplyr)
library(ggplot2)
library(tidyr)
library(stringr)
library(readxl)
library(viridis)
library(patchwork)
library(magrittr)
library(purrr)
library(extraDistr)
library(ggforce)
library(cowplot)
library(gridExtra)
```

Load functions
```{r source functions}
source("functions/calculate.zeta.R")
source("functions/cumulative.density.R")
source("functions/find.b.R")
source("functions/find.lambda.R")
source("functions/find.sigma.R")
source("functions/plot_cd_wide_zoom.R")
source("functions/plot_cd_zoom.R")
source("functions/sample.fun.R")
```


Set seed
```{r set seed}
set.seed(1984)
```


# Example scenarios

Load data. 
Data available from download at [DOI: 10.26188/11357126](https://dx.doi/10.26188/11357126)
```{r load example}
example_data <- readxl::read_xlsx(
  path = "data/species_publicity_data.xlsx",
  sheet = "example_case_studies"
)

example_data
```

Calculate decision score, zeta, for each option
```{r example_results}
example_results <- example_data %>%
  as_tibble %>%
  mutate(
    zeta = benefit - cost - probability_of_leaking*(cost_of_publicity - cost),
    involve = factor(
      involve,
      levels = c(
        "Government",
        "Local community group",
        "Public"
      )
    ),
    involve = recode(
      involve,
      "Local community group" = "Local community\ngroup",
    ),
    disclose = factor(
      disclose,
      levels = c(
        "Public",
        "Part",
        "Secret"
      ),
    )
  ) %>%
  group_by(species) %>%
  mutate(
    maxval = max(zeta)
  ) %>%
  ungroup %>%
  mutate(
    optimal = ifelse(
      zeta == maxval,
      "optimal",
      "not optimal"
    )
  ) %>%
  select(-maxval)

example_results
```

Plot results from examples scenarios in figure 1.
```{r figure_1}
figure_1 <- ggplot(example_results) +
  geom_hline(
    yintercept = 0,
    colour = "darkgrey"
  ) +
  geom_point(
    aes(
      x = disclose,
      y = zeta,
      colour = involve#,
      #shape = optimal,
      #size = optimal
    ),
    size = 4#,
    # position = position_dodge(
    #   width = 0.1
    # )
  ) +
  facet_wrap(
    facets = ~ species,
    nrow = 2,
    scales = "free_y"
  ) +
  theme_grey() +
  theme(
    axis.ticks.y = element_blank(),
    axis.text.y = element_blank(),
    axis.title.y = element_text(
      angle = 0
    ),
    axis.text.x = element_text(
      angle = 270
    )#,
    #legend.position = "bottom"
  ) +
  ylab("Decision\nscore") +
  xlab("Decision to disclose") +
  scale_colour_manual(
    values = plasma(11)[c(2,6,10)],
    guide = guide_legend(
      title = "Decision to\ninvolve",
      direction = "vertical"
    )
  ) #+
  # scale_shape_manual(
  #   values = c(19, 17),
  #   guide = guide_legend(
  #     title = "Decision is\noptimal?",
  #     direction = "vertical"
  #   )
  # ) +
  # scale_size_manual(
  #   values = c(4, 5),
  #   guide = guide_legend(
  #     title = "Decision is\noptimal?",
  #     direction = "vertical"
  #   )
  # ) 

figure_1
```

```{r figure_1 png}
png(
  filename = "plots/figure_1.png",
  width = 20,
  height = 12.5,
  units = "cm",
  res = 300
)
figure_1
dev.off()
```


```{r table_s2}
table_s2 <- example_results %>%
  select(species, disclose, involve, zeta) %>% 
  pivot_wider(
    names_from = disclose,
    values_from = zeta
  ) %>%
  arrange(species, involve)

table_s2
```

```{r write_table_s2}
write.csv(
  x = table_s2,
  file = "output/table_s2.csv",
  row.names = FALSE
)
```

# Questionnaire numeric estimates of cost and benefits
Read in numeric questionnaire data.
As above, data available from download at [DOI: 10.26188/11357126](https://dx.doi/10.26188/11357126)
```{r numeric_data}
numeric_data <- readxl::read_xlsx(
  path = "data/species_publicity_data.xlsx",
  sheet = "questionnaire_numeric"
)


numeric_data[ , which(substr(colnames(numeric_data), 1,1) == "q")] <- apply(
  X = numeric_data[ , which(substr(colnames(numeric_data), 1,1) == "q")],
  MARGIN = 2,
  FUN = as.numeric
)

numeric_data
```

Transform to long format and tidy
```{r numdat}
numdat <- numeric_data %>%
  select(
    -status,
    -currency,
    -consent_sp_id,
    -consent_sp_loc,
    -sp_loc
  ) %>%
  pivot_longer(
    cols = starts_with("q"),
    names_to = "question",
  ) %>%
  mutate(
    decision = substr(
      x = question,
      start = 2,
      stop = 2
    ) %>%
      as.integer,
    type = substr(
      x = question,
      start = 4,
      stop = 4
    ) %>%
      as.integer,
    estimate = substr(
      x = question,
      start = 6,
      stop = 6
    ) %>%
      as.integer
  ) %>%
  mutate(
    decision = case_when(
      decision == 1 ~ "Public",
      decision == 2 ~ "Part",
      decision == 3 ~ "Secret"
    ),
    decision = factor(
      x = decision,
      levels = c(
        "Public",
        "Part",
        "Secret"
      )
    ),
    type = case_when(
      type == 1 ~ "Benefit",
      type == 2 ~ "Cost",
      type == 3 ~ "pd"
    ),
    estimate = case_when(
      estimate == 1 ~ "lower",
      estimate == 2 ~ "upper",
      estimate == 3 ~ "most.likely",
      estimate == 4 ~ "pconf"
    )
  ) %>%
  select(
    species,
    decision,
    type,
    estimate,
    value
  ) %>%
  filter(
    species != "Species2"
  )

numdat
```


Separate out cost and benefit data, widen format, and determine type of distribution to be fit to estimates.

*Where values are not given for some parameters (in the case of the Wollemi pine), the most likely value is imputed as the arithmetic mean of the upper and lower values, and the confidence that the value between the upper and lower bounds is given as 0.95.*
```{r cost_benefit}
cost_benefit <- numdat %>%
  filter(type != "pd") %>%
  pivot_wider(
    names_from = estimate,
    values_from = value
  ) %>%
  mutate(
    pconf = pconf/100
  ) %>%
  filter(species != "Species2") %>%
  mutate(
    most.likely = ifelse(
      test = is.na(most.likely),
      yes = (upper - lower)/2,
      no = most.likely
    ),
    pconf = ifelse(
      test = is.na(pconf),
      yes = 0.95,
      no = pconf
    ),
    dist = case_when(
      pconf == 1 ~ "triangular",
      lower == most.likely & lower != 0 & pconf > 0.5 ~ "triangular",
      lower == most.likely ~ "exponential",
      TRUE ~ "lognormal"
    )
  )


cost_benefit
```

Estimate parameters for lognormal distributions
```{r cb_ln}
cost_benefit_lognormal <- cost_benefit %>%
  filter(dist == "lognormal") %>%
  mutate(
    value = pmap(
      .l = list(
        lower = lower,
        most.likely = most.likely,
        upper = upper,
        p = pconf
      ),
      .f = find.sigma,
      print.sigma = FALSE
    )
  ) %>%
  unnest(value)

cost_benefit_lognormal
```

Estimate parameters for exponential distributions
```{r cb_ex}
cost_benefit_exponential <- cost_benefit %>%
  filter(dist == "exponential") %>%
  mutate(
    value = pmap(
      .l = list(
        lower = lower,
        upper = upper,
        p = pconf
      ),
      .f = find.lambda,
    )
  ) %>%
  unnest(value)

cost_benefit_exponential
```

Estimate parameters for triangular distributions
```{r cb_tr}
cost_benefit_triangular <- cost_benefit %>%
  filter(dist == "triangular") %>%
  mutate(
    value = pmap(
      .l = list(
        lower = lower,
        most.likely = most.likely,
        upper = upper,
        p = pconf
      ),
      .f = find.b,
    )
  ) %>%
  unnest(value)

cost_benefit_triangular
```

Number of random samples to take from distributions
```{r nsamples}
nsamples <- 1000
```

Bind cost benefit data data back together and sample distributions for each variable
```{r cb_all}
cost_benefit_all <- bind_rows(
  cost_benefit_lognormal,
  cost_benefit_exponential,
  cost_benefit_triangular
) %>%
   mutate(
    samples = pmap(
      .l = list(
        dist = dist,
        most.likely = most.likely,
        value = value
      ),
      .f = sample.fun,
      n = nsamples
    )
  ) 

cost_benefit_all
```

Table s3. Generate table of statistical distributions and parameters used for sampling
```{r table_s3}
distribution_table <- cost_benefit_all %>%
  select(
    species,
    decision,
    type,
    lower,
    most.likely,
    everything(),
    -samples
  ) %>%
  mutate(
    sigma = ifelse(
      dist == "lognormal",
      value,
      NA
    ),
    mu = ifelse(
      dist == "lognormal",
      log(most.likely) + sigma^2,
      NA
    ),
    lambda = ifelse(
      dist == "exponential",
      value,
      NA
    ),
    a = ifelse(
      dist == "triangular",
      0,
      NA
    ),
    b = ifelse(
      dist == "triangular",
      value,
      NA
    ),
    c = ifelse(
      dist == "triangular",
      most.likely,
      NA
    )
  ) %>%
  select(
    -value
  ) %>%
  arrange(species, decision, type)

distribution_table
```

Write out table s3
```{r write_table_s3}
write.csv(
  x = distribution_table,
  file = "output/table_s3.csv",
  row.names = FALSE
)
```

Probability information becomes public given a decision, pd.
```{r pd}
pd_all <- numdat %>%
  filter(type == "pd") %>%
  mutate(
    pd = value/100
  ) %>%
  select(
    species,
    decision,
    pd
  )

pd_all
```

Cost if a decision is public, Cp
```{r cp}
cost_public <- cost_benefit_all %>%
  filter(decision == "Public") %>%
  filter(type == "Cost") %>%
  select(
    species,
    Cost_public = samples
  )

cost_public
```

Calculate estimates of decision score zeta from samples
```{r zeta_estimates}
zeta_estimates <- cost_benefit_all %>%
  pivot_wider(
    id_cols = c(
      species, decision
    ),
    names_from = type,
    values_from = samples
  ) %>%
  arrange(
    species, decision
  ) %>%
  full_join(
    cost_public
  ) %>%
  full_join(
    pd_all
  ) %>%
  mutate(
    pd = ifelse(
      test = is.na(pd),
      yes = 1,
      pd
    )
  ) %>%
  mutate(
    id = list(1:nsamples)
  ) %>%
  unnest(
    cols = c(
      Benefit,
      Cost,
      Cost_public,
      id
    )
  ) %>%
  mutate(
    zeta = pmap(
      .l = list(
        Bd = Benefit,
        Cd = Cost,
        pd = pd,
        Cp = Cost_public
      ),
      .f = calculate.zeta
    )
  ) %>%
  unnest(zeta) 

zeta_estimates
```

Minimum and maximum decision scores
```{r zeta_estimates_minmax}
zeta_estimates_minmax <- zeta_estimates %>%
  group_by(species) %>%
  summarise(
    min = min(zeta),
    max = max(zeta)
  )

zeta_estimates_minmax
```

Sequential values from min to max along each range
```{r zeta_ranges}
zeta_ranges <- zeta_estimates_minmax %$%
  mapply(
    FUN = seq,
    from = min,
    to = max,
    MoreArgs = list(
      length.out = 1000
    )
  )

colnames(zeta_ranges) <- zeta_estimates_minmax$species

zeta_ranges <- zeta_ranges %>%
  as_tibble %>%
  pivot_longer(
    cols = everything(),
    names_to = "species",
    values_to = "value"
  ) %>%
  arrange(species, value)

zeta_ranges
```

Calculate the cumlative density for each decision scenario at each value for each species
```{r zeta_cd_wide}
zeta_cd_wide <- zeta_ranges %>%
  rowwise %>%
  mutate(
    Public = cumulative.density(
      species,
      value,
      data = zeta_estimates,
      nsamples = nsamples,
      decision = "Public"
    ),
    Part = cumulative.density(
      species,
      value,
      data = zeta_estimates,
      nsamples = nsamples,
      decision = "Part"
    ),
    Secret = cumulative.density(
      species,
      value,
      data = zeta_estimates,
      nsamples = nsamples,
      decision = "Secret"
    )
  ) 

zeta_cd_wide
```

```{r zeta_dominance}
zeta_dominance <- zeta_cd_wide %>%
  group_by(species) %>%
  summarise(
    Public_Part   = sum(Public <= Part)   == nsamples,
    Public_Secret = sum(Public <= Secret) == nsamples,
    Part_Public   = sum(Part <= Public)   == nsamples,
    Part_Secret   = sum(Part <= Secret)   == nsamples,
    Secret_Public = sum(Secret <= Public) == nsamples,
    Secret_Part   = sum(Secret <= Part)   == nsamples
  ) %>%
  pivot_longer(
    cols = -species,
    names_to = "dominance",
    values_to = "dominates"
  ) %>%
  mutate(
    dominant = sub(
      pattern = "_.*",
      "",
      dominance
    ),
    dominated = sub(
      pattern = ".*_",
      "",
      dominance
    )
  ) %>%
  select(species, dominant, dominates, dominated) %>%
  pivot_wider(
    id_cols = c(species, dominant),
    names_from = dominated,
    values_from = dominates
  )

zeta_dominance
```

```{r write zeta_dominace}
# write.csv(
#   x = zeta_dominance,
#   file = "output/table_s6.csv",
#   row.names = FALSE
# )
```


```{r zeta_cd_limits}
Zeta_cd_limits <- zeta_estimates %>%
  group_by(species, decision) %>%
  summarise(
    low = quantile(
      x = zeta,
      probs = 0.1
    ),
    high = quantile(
      x = zeta,
      probs = 0.9
    )
  ) %>%
  group_by(species) %>%
  summarise(
    lower = min(low),
    upper = max(high)
  )

Zeta_cd_limits
```

Transform cumulative densitiy to long format for plotting
```{r zeta_cd}
zeta_cd <- zeta_cd_wide %>%
  pivot_longer(
    cols = c(Public, Part, Secret),
    names_to = "decision",
    values_to = "cd"
  ) %>%
  mutate(
    decision = factor(
      decision,
      levels = c(
        "Public",
        "Part",
        "Secret"
      )
    )
  )

zeta_cd
```


## Figure 2

```{r fig2_list}
fig2_list <- mapply(
  FUN = plot_cd_zoom,
  species = unique(zeta_estimates$species),
   xlab = c(
    FALSE,
    FALSE,
    FALSE,
    FALSE,
    FALSE,
    FALSE,
    FALSE
  ),
  ylab = c(
    FALSE,
    FALSE,
    FALSE,
    FALSE,
    FALSE,
    FALSE,
    FALSE
  ),
  # xtext = c(
  #   FALSE,
  #   FALSE,
  #   FALSE,
  #   FALSE,
  #   FALSE,
  #   TRUE,
  #   TRUE
  # ),
  xtext = TRUE,
  MoreArgs = list(
    data = zeta_cd,
    limits = Zeta_cd_limits,
    textsize = 10
  ),
  SIMPLIFY = FALSE
)

fig2_list[[1]]
```

```{r figure_2}
f2legend <- get_legend(fig2_list[[1]] + guides(col = guide_legend(ncol = 1))) 

figure_2 <- plot_grid(
  NULL,
  NULL,
  NULL,
  NULL,
  fig2_list[[1]] + scale_colour_manual(values = viridis(3)[3:1], guide = FALSE),
  fig2_list[[2]] + scale_colour_manual(values = viridis(3)[3:1], guide = FALSE),
  fig2_list[[3]] + scale_colour_manual(values = viridis(3)[3:1], guide = FALSE),
  fig2_list[[4]] + scale_colour_manual(values = viridis(3)[3:1], guide = FALSE),
  fig2_list[[5]] + scale_colour_manual(values = viridis(3)[3:1], guide = FALSE),
  fig2_list[[6]] + scale_colour_manual(values = viridis(3)[3:1], guide = FALSE),
  fig2_list[[7]] + scale_colour_manual(values = viridis(3)[3:1], guide = FALSE),
  f2legend,
  ncol = 4,
  align = "hv",
  labels = c(
    "",
    "",
    "",
    "",
    "Alpine tree-frog",
    "Giant Ibis",
    "Red Handfish",
    "Smooth Newt",
    "Species1",
    "Tiger",
    "Wollemi Pine"#,
    #""
  ),
  #label_fontface = "plain",
  label_size = 12,
  vjust = -0.01,
  hjust = 0,
  rel_heights = c(0.1, 1, 1),
  axis.titles = "collect"
) %>%
  grid.arrange(left = "Cumulative probability density", bottom = "Decision score")

figure_2
```

```{r print figure_2}
# png(
#   filename = "plots/figure_2_x.png",
#   width = 20,
#   height = 12.5,
#   units = "cm",
#   res = 300
# )
# figure_2
# dev.off() # some weird thing with gridExtra::grid.arrange makes this not work :shrug:

ggsave("plots/figure_2.png", plot = figure_2, width = 20, height = 12.5, units = "cm", dpi = 300)
```

### Figure s1

```{r figs1_list}
figs1_list <- mapply(
  FUN = plot_cd_wide_zoom,
  species = unique(zeta_estimates$species),
   xlab = c(
    FALSE,
    FALSE,
    FALSE,
    FALSE,
    FALSE,
    FALSE,
    FALSE
  ),
  ylab = c(
    FALSE,
    FALSE,
    FALSE,
    FALSE,
    FALSE,
    FALSE,
    FALSE
  ),
  # xtext = c(
  #   FALSE,
  #   FALSE,
  #   FALSE,
  #   FALSE,
  #   FALSE,
  #   TRUE,
  #   TRUE
  # ),
  xtext = TRUE,
  MoreArgs = list(
    data = zeta_cd,
    limits = Zeta_cd_limits,
    textsize = 10
  ),
  SIMPLIFY = FALSE
)

figs1_list[[1]]
```

```{r figure_s1}
fs1legend <- get_legend(figs1_list[[1]])

figure_s1 <- plot_grid(
  NULL,
  NULL,
  NULL,
  NULL,
  figs1_list[[1]] + scale_colour_manual(values = viridis(3)[3:1], guide = FALSE),
  figs1_list[[2]] + scale_colour_manual(values = viridis(3)[3:1], guide = FALSE),
  figs1_list[[3]] + scale_colour_manual(values = viridis(3)[3:1], guide = FALSE),
  figs1_list[[4]] + scale_colour_manual(values = viridis(3)[3:1], guide = FALSE),
  figs1_list[[5]] + scale_colour_manual(values = viridis(3)[3:1], guide = FALSE),
  figs1_list[[6]] + scale_colour_manual(values = viridis(3)[3:1], guide = FALSE),
  figs1_list[[7]] + scale_colour_manual(values = viridis(3)[3:1], guide = FALSE),
  #fdlegend,
  ncol = 4,
  align = "hv",
  labels = c(
    "",
    "",
    "",
    "",
    "Alpine tree-frog",
    "Giant Ibis",
    "Red Handfish",
    "Smooth Newt",
    "Species1",
    "Tiger",
    "Wollemi Pine"#,
    #""
  ),
  #label_fontface = "plain",
  label_size = 12,
  vjust = -0.01,
  hjust = 0,
  rel_heights = c(0.1, 1, 1)
)

figure_s1
```

Print figure_s1
```{r print figure_s1}
png(
  filename = "plots/figure_s1.png",
  width = 20,
  height = 20,
  units = "cm",
  res = 300
)
figure_s1
dev.off()
```


Find best option for each simulation, i.e., the one that maximises the decision score, and find whether the best option, as well as each disclosure decision would result in a value above zero, i.e., and expected net benefit.
```{r best_option}
best_option <- zeta_estimates %>%
  pivot_wider(
    id_cols = c(
      species,
      id
    ),
    names_from = decision,
    values_from = zeta
  ) %>%
  select(
    species, id, Public, Part, Secret
  ) %>%
  rowwise %>%
  mutate(
    zetamax = max(Public, Part, Secret),
    best = case_when(
      zetamax == Public ~ "Public",
      zetamax == Part ~ "Part",
      zetamax == Secret ~ "Secret"
    ),
    zetamin = min(Public, Part, Secret),
    worst = case_when(
      zetamin == Public ~ "Public",
      zetamin == Part ~ "Part",
      zetamin == Secret ~ "Secret"
    )
  ) %>%
  ungroup %>%
  mutate(
    max_positive = ifelse(
      test = zetamax > 1,
      yes = 1,
      no = 0
    ),
    pub_positive = ifelse(
      test = Public > 1,
      yes = 1,
      no = 0
    ),
    par_positive = ifelse(
      test = Part > 1,
      yes = 1,
      no = 0
    ),
    sec_positive = ifelse(
      test = Secret > 1,
      yes = 1,
      no = 0
    )
  )
  

best_option
```

Summarise data into proportions for plotting.
```{r best_option_proportion}
best_option_proportion <- best_option %>%
  group_by(species, best) %>%
  summarise(n = n()) %>%
  mutate(
    prop.best = n/1000,
    best = factor(
      x = best,
      levels = c(
        "Public",
        "Part",
        "Secret"
      )
    )
  )

best_option_proportion
```


```{r worst_option_proportion}
worst_option_proportion <- best_option %>%
  group_by(species, worst) %>%
  summarise(n = n()) %>%
  mutate(
    prop.worst = n/1000,
    worst = factor(
      x = worst,
      levels = c(
        "Public",
        "Part",
        "Secret"
      )
    )
  )

worst_option_proportion
```


## Figure 3.
Best (a) ...
```{r fig3a}
figure_3a <- ggplot(best_option_proportion) +
  geom_bar(
    aes(
      x = species,
      y = prop.best,
      fill = best
    ),
    position = "stack",
    stat = "identity"
  ) +
  ylab("Proportion of times decision is best") +
  xlab("Species") +
  scale_fill_manual(
    values = viridis(3)[3:1],
    guide = guide_legend(
      title = "Decision\nto disclose",
      direction = "vertical"
    )
    #guide = FALSE
  ) +
  ylim(c(0,1)) +
  theme(
    axis.text.x = element_text(
      angle = 270
    )
  )

figure_3a
```

... and worst (b) options. 
```{r fig3b}
figure_3b <- ggplot(worst_option_proportion) +
  geom_bar(
    aes(
      x = species,
      y = prop.worst,
      fill = worst
    ),
    position = "stack",
    stat = "identity"
  ) +
  ylab("Proportion of times decision is worst") +
  xlab("Species") +
  scale_fill_manual(
    values = viridis(3)[3:1],
    guide = guide_legend(
      title = "Decision\nto disclose",
      direction = "vertical"
    )
    #guide = FALSE
  ) +
  ylim(c(0,1)) +
  theme(
    axis.text.x = element_text(
      angle = 270
    )
    
  )

figure_3b
```

Figure 3
```{r figure_3}
figure_3 <- figure_3a +
  figure_3b  +
  plot_layout(guides = "collect") +
  plot_annotation(tag_levels = "a")

figure_3
```

```{r print figure_3}
png(
  filename = "plots/figure_3.png",
  width = 20,
  height = 12.5,
  units = "cm",
  res = 300
)
figure_3
dev.off()
```

Table of proportion of times best option
```{r best_option_table}
best_option_table <-  best_option_proportion %>%
  pivot_wider(
    id_cols = species,
    names_from = best,
    values_from = prop.best
  ) %>%
  select(
    species,
    Public,
    Part,
    Secret
  )

best_option_table
```

And worst option
```{r worst_option_table}
worst_option_table <-  worst_option_proportion %>%
  pivot_wider(
    id_cols = species,
    names_from = worst,
    values_from = prop.worst
  ) %>%
  select(
    species,
    Public,
    Part,
    Secret
  )

worst_option_table
```

Write out table s4a
```{r write_table_s4a}
write.csv(
  x = best_option_table,
  file = "output/table_s4a.csv",
  row.names = FALSE
)
```

Write out table s4b
```{r write_table_s4b}
write.csv(
  x = worst_option_table,
  file = "output/table_s4b.csv",
  row.names = FALSE
)
```

Table of proportion of times decision score is >0
```{r net_benefit_table}
net_benefit_table <- best_option %>%
  group_by(
    species
  ) %>%
  summarise(
    Best = sum(max_positive)/1000,
    Public = sum(pub_positive)/1000,
    Part = sum(par_positive)/1000,
    Secret = sum(sec_positive)/1000
  )

net_benefit_table
```

Write out table s5
```{r write_table_s5}
write.csv(
  x = net_benefit_table,
  file = "output/table_s5.csv",
  row.names = FALSE
)
```

Transform for plotting
```{r net_benefit_plot}
net_benefit_plot <- net_benefit_table %>%
  pivot_longer(
    cols = -species,
    names_to = "Decision",
    values_to = "Proportion positive"
  ) %>%
  mutate(
    Decision = factor(
      x = Decision,
      levels = c(
        "Best",
        "Public",
        "Part",
        "Secret"
      )
    )
  )

net_benefit_plot
```


### Figure s2
```{r figs2}
figure_s2 <- ggplot(net_benefit_plot) +
  geom_bar(
    aes(
      x = Decision,
      y = `Proportion positive`,
      fill = Decision
    ),
    stat = "identity"
  ) + 
  facet_wrap(
    facets = ~ species,
    nrow = 2
  ) +
  ylab("Proportion of times decision results in net benefit") +
  xlab("Decision to disclose") +
  scale_fill_manual(
    values = c("darkgrey", viridis(3)[3:1]),
    guide = guide_legend(
      title = "Decision\nto disclose",
      direction = "vertical"
    ) 
  ) +
  theme(
    axis.text.x = element_text(
      angle = 270
    ),
    legend.position = c(0.85, 0.2)
  ) +
  ylim(c(0,1))

figure_s2
```


```{r save figs2}
png(
  filename = "plots/figure_s2.png",
  width = 20,
  height = 12.5,
  units = "cm",
  res = 300
)
figure_s2
dev.off()
```





# Questionnaire written responses on types of costs and benefits
```{r text_dat}
text_dat <- readxl::read_xlsx(
  path = "data/species_publicity_data.xlsx",
  sheet = "questionnaire_text"
)
```

```{r text_results}
q_numeric <- sapply(text_dat, is.numeric)
sum_resp <- sapply(
  X = text_dat,
  FUN = function(x){
    
    isnum <- is.numeric(x)
    
    if(isnum){
      result <- sum(x, na.rm = TRUE)
    } else {
      result <- NA
    }
    
  }
)

perc_resp <- sapply(
  X = text_dat,
  FUN = function(x){
    
    isnum <- is.numeric(x)
    
    lenx <- length(x)
    
    if(isnum){
      result <- 100*sum(x, na.rm = TRUE)/lenx
    } else {
      result <- NA
    }
    
  }
)

text_results <- tibble(
  question = colnames(text_dat),
  sum_responses = sum_resp,
  percentage_responses = perc_resp
)

text_results
```

```{r all_results}
all_results <- text_results %>%
  filter(substr(question, 1, 3) == "all") %>%
  mutate(
    c1 = sub(
      pattern = ".{0,4}_",
      replacement = "",
      question
    ),
    c2 = sub(
      pattern = ".{0,8}_",
      replacement = "",
      c1
    ),
    cost_benefit= sub(
      pattern = "_.*",
      replacement = "",
      c1
    ),
    group = sub(
      pattern = "_.*",
      replacement = "",
      c2
    ),
    response = sub(
      pattern = ".{0,8}_",
      replacement = "",
      c2
    )
  ) %>%
  select(cost_benefit, group, response, sum_responses, percentage_responses) %>%
  mutate(
    response = recode(
      response,
      "funding" = "Funding",
      "awareness-publicsupport" = "Awareness &\npublic support",
      "protection-recovery" = "Protection &\nrecovery",
      "species" = "Species",
      "individuals" = "Individual\npeople",
      "researchersandinstitutes" = "Researchers &\ninstitutes",
      "gov" = "Government",
      "business" = "Business",
      "ngo" = "NGOs",
      "genpublic" = "General\npublic",
      "genpub" = "General\npublic",
      "poaching" = "Poaching",
      "extinctionextirpation" = "Extinction &\nextirpation",
      "disease" = "Disease",
      "funding-financial" = "Funding",
      "reputation" = "Reputation"
    )
  )

all_results
```

Split each kind for plotting
```{r separate}
benefits <- all_results %>% filter(cost_benefit == "benefits")
costs <- all_results %>% filter(cost_benefit == "costs")

benefits_kind     <- benefits %>% filter(group == "kind")
benefits_entities <- benefits %>% filter(group == "entities")
benefits_source   <- benefits %>% filter(group == "source")

costs_kind     <- costs %>% filter(group == "kind")
costs_entities <- costs %>% filter(group == "entities")
```


```{r fig4a}
figure_4a <- ggplot(benefits_kind) +
  geom_bar(
    aes(
      x = response,
      y = percentage_responses
    ),
    fill = viridis(9)[8],
    stat = "identity",
    width = 0.9*3/7
  ) +
  #ylab("Percentage of informants") +
  xlab("Kinds of benefit") +
  ylim(0, 100) +
  coord_flip() +
  theme(
    axis.title.x = element_blank()
  )

figure_4a
```

```{r fig4c}
figure_4c <- ggplot(benefits_entities) +
  geom_bar(
    aes(
      x = response,
      y = percentage_responses
    ),
    fill = viridis(9)[8],
    stat = "identity",
    width = 0.9
  ) +
  #ylab("Percentage of informants") +
  xlab("Entities\nrecieving benefits") +
  ylim(0, 100) +
  theme(
    axis.text.y = element_text(
      size = 7,
      lineheight = 0.7
    ),
    axis.title.x = element_blank()
  )+
  coord_flip()

figure_4c
```

```{r fig4e}
figure_4e <- ggplot(benefits_source) +
  geom_bar(
    aes(
      x = response,
      y = percentage_responses
    ),
    fill = viridis(9)[8],
    stat = "identity",
    width = 0.9*3/7
  ) +
  ylab("Percentage of informants") +
  xlab("Source of benefits") +
  ylim(0, 100) +
  coord_flip()

figure_4e
```

```{r fig4b}
figure_4b <- ggplot(costs_kind) +
  geom_bar(
    aes(
      x = response,
      y = percentage_responses
    ),
    fill = viridis(9)[2],
    stat = "identity",
    width = 0.9*5/7
  ) +
  ylab("Percentage of informants") +
  xlab("Kinds of cost") +
  ylim(0, 100) +
  theme(
    axis.text.y = element_text(
      lineheight = 0.7
    ),
    axis.title.x = element_blank()
  ) +
  coord_flip()

figure_4b
```

```{r fig4d}
figure_4d <- ggplot(costs_entities) +
  geom_bar(
    aes(
      x = response,
      y = percentage_responses
    ),
    fill = viridis(9)[2],
    stat = "identity",
    width = 0.9*3/7
  ) +
  ylab("Percentage of informants") +
  xlab("Entities\nbearing costs") +
  ylim(0, 100) +
  coord_flip()

figure_4d
```

```{r fig4}
figure_4 <- figure_4a +
  figure_4b +
  figure_4c +
  figure_4d +
  figure_4e +
  plot_layout(
    ncol = 2
  ) +
  plot_annotation(
    tag_levels = "a"
  )

figure_4
```

```{r save fig4}
png(
  filename = "plots/figure_4.png",
  width = 20,
  height = 15,
  units = "cm",
  res = 300
)
figure_4
dev.off()
```