diff --git a/documents/ovando-etal-assessing-global-fisheries-long.Rmd b/documents/ovando-etal-assessing-global-fisheries-long.Rmd index 16096d5..f300b19 100644 --- a/documents/ovando-etal-assessing-global-fisheries-long.Rmd +++ b/documents/ovando-etal-assessing-global-fisheries-long.Rmd @@ -17,7 +17,7 @@ output: reference_docx: word-template.docx bookdown::pdf_document2: default params: - results_name: ["v0.5"] + results_name: ["v1.0"] min_years_catch: [25] abstract: | Assessments of the global state of fisheries play an important role in shaping the public narrative around ocean health, motivating future directions of research and funding, formulating evidence-based policy and tracking the implementation of the United Nations Sustainable Development Goals. While we have reliable estimates of stock status for fisheries accounting for 50% of global catch, our knowledge of the state of the remaining 50%, the worlds 'unassessed' fisheries, is poor. Numerous high-profile publications featuring a range of statistical methods have produced estimates of the global status of these unassessed fisheries, but limited quantity and quality of data along with methodological differences have produced counterintuitive and conflicting results. This is especially true in areas such as Southeast Asia and Africa which are also regions that have high dependence on fishery resources. How can we effectively estimate the state of global fishery sustainability and track progress towards the Sustainable Development Goals targets? We developed a flexible assessment model to assess the value of different kinds, quantities, and quality of data in improving estimates of fishery stock status. We then explore avenues for obtaining potentially impactful data, including through the use of local expert opinion through Fisheries Management Index scores, and increasingly available but historically underutilized data such as trawl footprints and effort data. These data are then used to illustrate how different types of information paint starkly different pictures of the state of fisheries around the world, and to identify priority data types for future collection. Our results provide further evidence that relying on catch data alone results in inaccurate classification of fishery status, often performing only as well or worse than a random guess. Obtaining accurate estimates of stock status for the world's unassessed fisheries depends on prioritizing the collection of high-priority data at a global scale, not on the development of new modeling methods alone. diff --git a/documents/ovando-etal-assessing-global-fisheries.Rmd b/documents/ovando-etal-assessing-global-fisheries.Rmd index 39a23ec..211165b 100644 --- a/documents/ovando-etal-assessing-global-fisheries.Rmd +++ b/documents/ovando-etal-assessing-global-fisheries.Rmd @@ -1,5 +1,7 @@ --- -title: "Global Assessments of Fishery Status need Better Data more than Better Models" +title: | + | Global Assessments of Fishery Status need + | Better Data more than Better Models author: - Daniel Ovando - Ray Hilborn @@ -13,14 +15,14 @@ date: "`r Sys.Date()`" bibliography: ["../references.bib"] csl: nature.csl output: + bookdown::pdf_document2: default bookdown::word_document2: reference_docx: word-template.docx - bookdown::pdf_document2: default params: - results_name: ["v0.5"] + results_name: ["v1.0"] min_years_catch: [25] abstract: | - Assessments of the global state of fisheries play an important role in tracking the implementation of the United Nations Sustainable Development Goals. While we have reliable estimates of stock status for fisheries accounting for 49% of global catch, our knowledge of the state of the remaining 51%, the worlds 'unassessed' fisheries, is poor. Numerous high-profile publications have produced estimates of the global status of these unassessed fisheries, but limited quantity and quality of data along with methodological differences have produced counterintuitive and conflicting results. Here, we show that despite numerous efforts, our understanding of the status of global fisheries remains poor, even when new sources of broadly available data are added. Obtaining accurate estimates of stock status for the world's fisheries depends on prioritizing the collection of high-priority data at a global scale, not on the development of new modeling methods alone. + Assessments of the global state of fisheries play an important role in tracking the implementation of the United Nations Sustainable Development Goals. While we have reliable estimates of stock status for fisheries accounting for 49% of global catch, our knowledge of the state of the remaining 51%, the world's 'unassessed' fisheries, is poor. Numerous high-profile publications have produced estimates of the global status of these unassessed fisheries, but limited quantity and quality of data along with methodological differences have produced counterintuitive and conflicting results. Here, we show that despite numerous efforts, our understanding of the status of global fisheries remains poor, even when new sources of broadly available data are added. Obtaining accurate estimates of stock status for the world's fisheries depends on prioritizing the collection of high-priority data at around the globe, rather than the development of new modeling methods alone. always_allow_html: true linkcolor: blue --- @@ -109,22 +111,38 @@ breaks <- c(0, breaks, Inf) labels <- c("over", "fully", "under") -rand_fits <- assess_ram_fits %>% - filter(data == "ram-data") %>% - mutate(mean = sample(c(0.1, 1, 2), n(), replace = TRUE), - data = "guess") +# rand_fits <- assess_ram_fits %>% +# filter(data == "ram-data") %>% +# mutate(mean = sample(c(0.1, 1, 2), n(), replace = TRUE), +# data = "guess") + perf <- assess_ram_fits %>% - bind_rows(rand_fits) %>% + # bind_rows(rand_fits) %>% mutate(sraplus_bin = cut(mean, breaks = breaks, labels = labels)) %>% mutate(ram_bin = cut(ram_b_v_bmsy, breaks = breaks, labels = labels)) %>% + mutate(data = as.factor(data)) %>% + mutate(data = forcats::fct_recode(data, + "RAM Index" = "ram-data", + "SAR" = "sar", + "FMI" = "fmi", + "Effective CPUE" = "cpue", + "Effective CPUE+" = "cpue-plus", + "Nominal CPUE" = "nominal-cpue", + "Nominal CPUE+" ="nominal-cpue-plus", + "CMSY" = "cmsy", + "Guess" = "guess", + "RAM U/Umsy" = "u_umsy" + + )) %>% group_by(data) %>% summarise( mpe = median((mean - ram_b_v_bmsy ) / ram_b_v_bmsy), mape = median(abs(mean - ram_b_v_bmsy ) / ram_b_v_bmsy), accuracy = mean(ram_bin == sraplus_bin, na.rm = TRUE) ) %>% - arrange((mape)) + arrange((mape)) %>% + rename("Data Used" = "data") @@ -192,27 +210,27 @@ theme_set(pub_theme) -# Publication Requirements + -Target Journal: [Nature Sustainability](https://www.nature.com/natsustain/about/content) + -3500 words main text + -50 references (main text) + -3000 words (max) methods + # Introduction -The United Nations Sustainable Development Goal 14 (SDG 14), related to "Life under water", calls for the global community to "Conserve and sustainably use the oceans, seas and marine resources for sustainable development". Meeting these targets depends in part on our ability to effectively track the status of global marine fish stocks. The Food and Agriculture Organization of the United Nations' (FAO) State of World Fisheries and Aquaculture (SOFIA) report is the standard source for tracking the global state of fisheries. As of the most recent report, the FAO estimates that as of 2017 59.6% of marine fish stocks are maximally sustainably fished, 6.2% are underfished, and 34.2% are overfished [@fao2020]. While foundational, the SOFIA assessment was designed in the 1970s based on the then available data and methods, With the surge in data availability and models designed for data-limited stocks, new global assessment methods are needed to meet the new demand for estimation and tracking progress towards the SDG goals. +The United Nations Sustainable Development Goal 14 (SDG 14), focusing on "Life under water", calls for the global community to "Conserve and sustainably use the oceans, seas and marine resources for sustainable development". Meeting these targets depends in part on our ability to effectively track the status of global marine fish stocks. The Food and Agriculture Organization of the United Nations' (FAO) State of World Fisheries and Aquaculture (SOFIA) report is the most widely used source for tracking the global state of fisheries. As of the most recent report, the FAO estimates that as of 2017 59.6% of marine fish stocks are maximally sustainably fished, 6.2% are underfished, and 34.2% are overfished [@fao2020]. The SOFIA assessment was designed in the 1970s based on the then available data and methods, With the surge in data availability, such as global assessments of management strength, trawl footprints, and fishing effort, and modeling methods, fresh approaches are needed to meet the new demand for tracking progress towards the SDG goals and supporting the research targets set out in the UNited Nations Decade on Oceans and Science plan. -Forty nine percent of landings of marine fishes reported by the FAO are represented in the RAM Legacy Stock Assessment Database, a repository of "gold standard" estimates of fishery status [@hilborn2020]. The SOFIA process uses these formal assessments for their determinations of stock status wherever possible. However, that leaves 51% of global fisheries landings lacking in formal stock assessments. While these "unassessed" stocks are generally individually smaller than the typically larger and more valuable stocks in the assessed category, collectively they are a vital source of food, employment, cultural value, and ecosystem services around the world. +Forty nine percent of landings of marine fishes reported by the FAO are represented in the RAM Legacy Stock Assessment Database [@ricard2012], a repository of best-available estimates of fishery status [@hilborn2020]. The SOFIA process uses these formal assessments for their determinations of stock status wherever possible. However, that leaves 51% of global fisheries landings lacking in formal stock assessments. While these "unassessed" stocks are generally individually smaller than the typically larger and more valuable stocks in the assessed category, collectively they are a vital source of food, employment, cultural value, and ecosystem services around the world. -The SOFIA report bases its estimates of included unassessed stocks mostly on expert opinion. While local experts can be well informed as to the status of their fish stocks, a more quantitative and reproducible process would be desirable. Numerous studies in recent years have put forward versions of "data-limited" models that have attempted to provide estimates for the global status of these unassessed stocks [@costello2012b; @costello2016a; @pauly2007; @rosenberg2018; @thorson2012a]. Due to data limitations, all of these global assessment efforts use forms of "catch-only" stock assessment models (Free et al. 2020 [@free2020] and references therein). These models seek to infer stock status, for example in terms of biomass *B* relative to the biomass at maximum sustainable yield *B~MSY~*, from characteristics of a fishery's catch history, for example the ratio of catch to maximum catch [@martell2013]. +The SOFIA report bases its estimates of included unassessed stocks mostly on expert opinion. While local experts can be well informed as to the status of their fish stocks, a more quantitative and reproducible process would be desirable, in order to increase reproducibility, transparency, and potential reliability. Numerous studies in recent years have put forward versions of "data-limited" models that have attempted to provide estimates for the global status of these unassessed stocks [@costello2012b; @costello2016a; @pauly2007; @rosenberg2018; @thorson2012a]. Due to data limitations, all of these global assessment efforts use forms of "catch-only" stock assessment models (Free et al. 2020 [@free2020] and references therein). These models seek to infer stock status, for example in terms of biomass *B* relative to the biomass at maximum sustainable yield *B~MSY~*, from characteristics of a fishery's catch history, for example the ratio of catch to maximum catch [@martell2013]. -However, Free et al. 2020 [@free2020] demonstrated that these catch-only models can often produce both imprecise and biased estimates of current stock status in terms of B/B~MSY~. These issues become apparent when we consider some of the macro-level predictions made by these models. The RAM Legacy Stock Assessment Database [@ricard2012] contains the best available estimates of B/B~MSY~ and other fishery reference points for hundreds of fisheries [@hilborn2020]. While it must noted that the estimates in RAM are themselves model outputs subject to their own non-trivial errors and biases, a simple benchmark is to compare the best available estimates of fishery status from RAM to those predicted by potentially less reliable methods intended for use when insufficient data are available for a full stock assessment model. +However, Free et al. 2020 [@free2020] demonstrated that these catch-only models can often produce both imprecise and biased estimates of current stock status in terms of B/B~MSY~. These issues become apparent when we consider some of the macro-level predictions made by these models. While the estimates B/B~MSY~ and other reference points reported in RAM are themselves model outputs subject to their own non-trivial errors and biases, a simple benchmark is to compare the estimates of fishery status from RAM to those predicted by potentially less reliable methods intended for use when insufficient data are available for a full stock assessment model. -Costello et al. 2016 [@costello2016a] finds similar rankings of regions in terms of stock status as RAM, but their estimate of state of fisheries in the Mediterranean/Black Sea regions and Southeast Asia seem to be over-optimistic, and the Northeast Pacific should be better by comparison. Rosenberg et al. 2018 [@rosenberg2018] demonstrates the same problem, with stocks in Southeast Asia estimated as doing better than the Northeast Atlantic or Northeast Pacific. The Pauly 2007 [@pauly2007] catch based approach finds the stocks of Southeast Asia in much better condition than the Northeast Pacific or Northeast Atlantic (Table.\@ref(tab:status-tab)). The methods besides RAM in Table.\@ref(tab:status-tab) include both formally assessed and unassessed stocks, and as such we would expect them to differ broadly in their estimates of regional stock status, particularly since we might expect unassessed stocks to be less rigorously managed and by extension have poorer stock status. However, the lack of consistency across heavily assessed regions such as the Northeast Pacific, and the lack of contrast in stock status between heavily and lightly managed regions is concerning. +Costello et al. 2016 [@costello2016a] finds similar rankings of regions in terms of stock status as RAM, but their estimate of state of fisheries in the Mediterranean/Black Sea regions and Southeast Asia seem to be over-optimistic, and the Northeast Pacific should be better by comparison. The super-ensemble method [@anderson2017a] used in Rosenberg et al. 2018 [@rosenberg2018] demonstrates the same problem, with stocks in Southeast Asia estimated as doing better than the Northeast Atlantic or Northeast Pacific. The Pauly 2007 [@pauly2007] catch based approach finds the stocks of Southeast Asia in much better condition than the Northeast Pacific or Northeast Atlantic (Table.\@ref(tab:status-tab)). The methods besides RAM in Table.\@ref(tab:status-tab) include both formally assessed and unassessed stocks, and as such we would expect them to differ broadly in their estimates of regional stock status, particularly since we might expect unassessed stocks to be less rigorously managed and by extension have poorer stock status. However, the lack of consistency across heavily assessed regions such as the Northeast Pacific, and the lack of contrast in stock status between heavily and lightly managed regions is concerning. ```{r status-tab} @@ -257,7 +275,7 @@ readr::read_csv(here("data","status-table.csv")) %>% left_join(state_space_results, by = c("FAO Area" = "fao_area")) %>% arrange(`FAO % Overfished`) %>% select(`FAO Area`,`FAO % Overfished`, everything()) %>% - knitr::kable(caption = "Estimates of B/B~MSY~ by FAO, RAM through Hilborn et al. 2020, Costello et al. 2016, + knitr::kable(caption = "Estimates of B/B~MSY~ by FAO, RAM Legacy Stock Assessment Database through Hilborn et al. 2020, Costello et al. 2016, Rosenberg et al. 2018, and Pauly 2007", booktabs = TRUE, format = "pipe") %>% @@ -266,7 +284,7 @@ readr::read_csv(here("data","status-table.csv")) %>% ``` -In this paper we ask, can combining the FAO's catch statistics with other broadly available data improve our understanding of the state of global fisheries? We use a flexible surplus-production based stock assessment package, [`sraplus`](www.github.com/danovando/sraplus) to demonstrate how different sources of data can be used to augment catch-only models at a global scale, and to evaluate how our perception of global stock status would vary depending on which sources of data we include. We show that our understanding of global fishery status is poor, and that improvements depend on an redoubled effort at global data collection. +In this paper we ask, can combining the FAO's catch statistics with other broadly available data improve our understanding of the state of global fisheries? We use a flexible surplus-production based stock assessment software package, [`sraplus`](www.github.com/danovando/sraplus) to demonstrate how different data sources can be used to augment catch-only models at a global scale, and to evaluate how our perception of global stock status would vary depending on which sources of data we include. We show that our understanding of global fishery status is poor, and that improvements depend on a redoubled effort at global data collection and synthesis. # Results @@ -288,11 +306,11 @@ We first present a case study demonstrating how different kinds of data can lead -For our case study, we selected `r n_distinct(exs$stockid)` stocks for which we have stock specific FMI and SAR scores. We then paired effort data at the resolution of year, country, and FAO statistical area from Rousseau et al. 2019 [@rousseau2019] to each stock. We first used the catch history heuristics internal to CMSY [@froese2017] to estimate stock status. We then used stock-specific data on SAR and FMI to generate priors on F/F~MSY~ for each of the stocks, which were then passed to `sraplus`. Lastly, we used the reconstructed effort data from [@rousseau2019] to create an index of abundance for each stock, and estimated stock status by fitting to this index while using priors on fishing mortality rates informed by each stock's FMI and SAR values. While CMSY systemically overestimated fishing mortality rates and underestimated stock status, use of the SAR, FMI, and effort data produced substantially more accurate results (Fig.\@ref(fig:cs-plot)). +For our case study, we selected `r n_distinct(exs$stockid)` stocks for which we have stock specific FMI and SAR scores. We then paired effort data at the resolution of year, country, and FAO statistical area from Rousseau et al. 2019 [@rousseau2019] to each stock. As a benchmark, we first estimated stock status for these case study fisheries using the CMSY [@froese2017] method, as this has become one of the most widely used catch-only models currently available. We then used stock-specific data on SAR and FMI to generate priors on F/F~MSY~ for each of the stocks, which were then passed to sraplus. Lastly, we used the reconstructed effort data from [@rousseau2019] to create an index of abundance for each stock, and estimated stock status by fitting to this index while using priors on fishing mortality rates informed by each stock's FMI and SAR values. While CMSY systemically overestimated fishing mortality rates and underestimated stock status, use of the SAR, FMI, and effort data produced substantially more accurate results (Fig.\@ref(fig:cs-plot)). -```{r cs-plot, fig.cap="RAM values of B/B~MSY~ and F/F~MSY~ (x-axes) for case study fisheries plotted against estimated values (y-axes) using CMSY [@froese2017], priors informed by stock-specific Fisheries Management Index (FMI) and swept area ratio (SAR) scores, and an abundance index based on reconstructed effort trends assuming a rate of technological increase of 2.6%. Each point is a stock, point size is a function of stock size. Black dashed line shows the 1:1 relationship."} +```{r cs-plot, fig.cap="RAM values of B/B~MSY~ and F/F~MSY~ (x-axes) for case study fisheries plotted against estimated values (y-axes) using CMSY [@froese2017], priors informed by stock-specific Fisheries Management Index (FMI) and swept area ratio (SAR) scores, and an abundance index based on reconstructed effort trends assuming a rate of technological increase of 2.6%. Each point is a stock. Black dashed line shows the 1:1 relationship."} ex_scatter_plot ``` @@ -300,37 +318,46 @@ ex_scatter_plot We next assessed the ability of FMI, SAR, and effort data to improve estimates of global stock status. We based this test on `r ram_comp_data$stockid %>% n_distinct()` fisheries from RAM, covering `r n_distinct(ram_comp_data$isscaap_group)` broad taxonomic groups, with estimates of B/B~MSY~ and greater than 25 years of continuous catch history. B/B~MSY~ values from RAM are themselves estimates, not data, but they are the best available information on global stock status. We then paired the catch histories for these RAM stocks with regional-level SAR, FMI, and effort data. This process approximates a global-level assessment exercise, where data are available at regional levels, but not for specific fisheries. -As a proof of concept, we also estimated B/B~MSY~ of our candidate RAM stocks by using `sraplus` to fit to an abundance index drawn directly from RAM. We then fit a range of models utilizing different combinations FMI, SAR, and effort data, along with the CMSY catch-only method described in Froese et al. 2017 [@froese2017] (See Table.S1). We assessed performance using three metrics: median percent error (MPE, a measure of bias), median absolute percent error (MAPE, a measure of accuracy), and classification accuracy. Classification accuracy is calculated as the proportion of times that use of a given combination of data resulted in a stock being classified into the correct FAO status classification (one of underfished, maximally sustainably fished, and overfished). +We also estimated B/B~MSY~ values of our candidate RAM stocks by using sraplus to fit to an abundance index drawn directly from RAM. We then fit a range of models utilizing different combinations FMI, SAR, and effort data, along with the CMSY method described in Froese et al. 2017 [@froese2017] (See Table.\@ref(tab:dat-desc)). We assessed performance using three metrics: median percent error (MPE, a measure of bias), median absolute percent error (MAPE, a measure of accuracy), and classification accuracy. Classification accuracy was calculated as the proportion of times that use of a given combination of data resulted in a stock being classified into the correct FAO status classification (one of underfished, maximally sustainably fished, and overfished). -Overall the `sraplus` estimates of B/B~MSY~ resulting from using the RAM data are reasonably good (median absolute percent error `r percent(perf$mape[perf$data == "ram-data"])`, accuracy = `r percent(perf$accuracy[perf$data == "ram-data"])`, Table.\@ref(tab:perf-tab), Fig.\@ref(fig:mpe-map)-\@ref(fig:acc-map)). This exercise tells us that given sufficiently high quality data, a surplus production model such as `sraplus` is reasonably capable of reproducing the global state of fisheries as understood from formally assessed fisheries. +Overall the sraplus estimates of B/B~MSY~ resulting from using the RAM data were relatively accurate and unbiased at a macro level (MPE `r percent(perf$mpe[perf$"Data Used" == "ram-data"])`,MAPE `r percent(perf$mape[perf$"Data Used" == "ram-data"])`, accuracy = `r percent(perf$accuracy[perf$"Data Used" == "ram-data"])`, Table.\@ref(tab:perf-tab), Fig.\@ref(fig:mpe-map)-\@ref(fig:acc-map)). This exercise tells us that given sufficiently high quality data, a surplus production model such as sraplus is reasonably capable of reproducing the global state of fisheries as understood from formally assessed fisheries. -Performance limitations then are likely to arise less from model misspecification than from the quality of the data themselves. These becomes clearer once we consider the performance of `sraplus` models fit to combinations of our broadly available datasets. Many of the datasets used produced similar levels of bias as the RAM data (Table.\@ref(tab:perf-tab)). However, this is somewhat an artifact of the data. The status of most stocks in RAM is also relatively good, with recent B/B~MSY~ values generally near one. This means that a model that more or less reproduces the global average of stock status will be relatively unbiased on average, but imprecise. Focusing on MAPE instead, the error of the models jumps dramatically as soon as data other than RAM are used, to a minimum value of `r percent(perf$mape[2])` and a maximum of `r percent(max(perf$mape))`. The mean accuracy across all non-RAM data fits was only `r percent(mean(perf$accuracy[!perf$data %in% c("ram-data","guess")]))`. Note that there are only three bins in the FAO stock status classifications, and a "model" that randomly assigns a stock to a status category has a mean accuracy of `r percent(perf$accuracy[perf$data == "guess"])`. +Performance limitations then are likely to arise less from model misspecification than from the quality of the data themselves. This becomes clearer once we consider the performance of sraplus models fit to combinations of our broadly available datasets. Many of the datasets used produced similar levels of bias as the RAM data (Table.\@ref(tab:perf-tab)). However, this is somewhat an artifact of the data. The status of most stocks in RAM is also relatively good, with recent B/B~MSY~ values generally near one. This means that a model that more or less reproduces the global average of stock status will be relatively unbiased on average, but imprecise. Focusing on MAPE instead, the error of the models jumps dramatically as soon as data other than RAM are used, to a minimum value of `r percent(perf$mape[2])` and a maximum of `r percent(max(perf$mape))`. The mean accuracy across all non-RAM data fits was only `r percent(mean(perf$accuracy[!perf$"Data Used" %in% c("ram-data","guess")]))`. Note that there are only three bins in the FAO stock status classifications, and a "model" that randomly assigns a stock to a status category has a mean accuracy of `r percent(perf$accuracy[perf$"Data Used" == "guess"])`. Looking geographically we see a similar pattern of a rapid decrease in performance for models using non-RAM data intended to simulate a global assessment process. Across the models, performance was not consistent in space: use of different data performed best or worst for different FAO regions. For example, models fit to nominal CPUE data substantially overestimate stock status in the Mediterranean, while models based on data using effective CPUE perform better in that region (but worse in others) Fig.\@ref(fig:mape-map). We find similarly inconsistent performance for both bias (Fig.\@ref(fig:mpe-map)) and accuracy (Fig.\@ref(fig:acc-map)). Overall, while some data sources performed slightly better than others by some metrics in some places, no models using any non-RAM data were able to capture the overall state or geographic distribution of stock status represented in RAM in a consistently satisfactory manner. ```{r perf-tab} knitr::kable(perf, digits = 2, - caption = "Global performance statistics in the most recent year available of models using different sources of data. mpe = median percent error (bias), mape = median absolute percent error (error), accuracy = percent of times that stocks were classified to the correct FAO status bin (underfished, maximally sustainably fished, overfished). Performance is judged relative to reported values in RAM Legacy Stock Assessment Database.", + caption = "Global performance statistics in the most recent year available of models using different sources of data. mpe = median percent error (bias), mape = median absolute percent error (error), accuracy = percent of times that stocks were classified to the correct FAO status bin (underfished, maximally sustainably fished, overfished). Performance is judged relative to reported values in RAM Legacy Stock Assessment Database. See Table.5 for details of data types.", format = "pipe") ``` -```{r mpe-map, fig.width=8, fig.height=8, fig.cap="Median percent error in most recent B/B~MSY~ by FAO statistical area from different data sources. ram-data refers to catch and abundance index drawn from RAM. CPUE refers to an index of abundance based on reconstructed effort data. cpue-plus uses CPUE along with Fisheries Management Index (FMI) and/or wept area ratio (SAR) data. For both CPUE series 'nominal' assumes a 0% technology rate, otherwise a 2.6% technology rate is assumed. ram_u_umsy assumes all fisheries in the region share a common U/Umsy series with formally assessed fisheries in the region. fmi uses fmi scores to develop a prior on recent fishing mortality rates, sar does the same but based on swept area ratio. CMSY uses the methods from Froese et al. 2017 [@froese2017]."} +```{r mpe-map, fig.width=8, fig.height=8, fig.cap="Median percent error in most recent B/B~MSY~ by FAO statistical area from different data sources. RAM Index refers to catch and abundance index drawn from RAM. Effective CPUE refers to an index of abundance based on reconstructed effort data. Effective CPUE+ uses CPUE along with Fisheries Management Index (FMI) and/or wept area ratio (SAR) data. For both CPUE series 'nominal' assumes a 0% technology creep, for effective a 2.6% technology creep is assumed. RAM U/U~MSY~ assumes all fisheries in the region share a common U/U~MSY~ series with formally assessed fisheries in the region. FMI uses FMI scores to develop a prior on recent fishing mortality rates, SAR does the same but based on swept area ratio. CMSY uses the methods from Froese et al. 2017 [@froese2017]. Guess assignes a random recent B/B~MSY~ of 0.4,1, or 1.6."} ram_mpe_map_plot ``` -```{r mape-map, fig.cap = "Median absolute percent error in most recent B/B~MSY~ by FAO statistical area from different data sources. ram-data refers to catch and abundance index drawn from RAM. CPUE refers to an index of abundance based on reconstructed effort data. cpue-plus uses CPUE along with Fisheries Management Index (FMI) and/or wept area ratio (SAR) data. For both CPUE series 'nominal' assumes a 0% technology rate, otherwise a 2.6% technology rate is assumed. ram_u_umsy assumes all fisheries in the region share a common U/Umsy series with formally assessed fisheries in the region. fmi uses fmi scores to develop a prior on recent fishing mortality rates, sar does the same but based on swept area ratio. CMSY uses the methods from Froese et al. 2017 [@froese2017]."} +```{r mape-map, fig.cap = "Median absolute percent error in most recent B/B~MSY~ by FAO statistical area from different data sources. RAM Index refers to catch and abundance index drawn from RAM. Effective CPUE refers to an index of abundance based on reconstructed effort data. Effective CPUE+ uses CPUE along with Fisheries Management Index (FMI) and/or wept area ratio (SAR) data. For both CPUE series 'nominal' assumes a 0% technology creep, for effective a 2.6% technology creep is assumed. RAM U/U~MSY~ assumes all fisheries in the region share a common U/U~MSY~ series with formally assessed fisheries in the region. FMI uses FMI scores to develop a prior on recent fishing mortality rates, SAR does the same but based on swept area ratio. CMSY uses the methods from Froese et al. 2017 [@froese2017]. Guess assignes a random recent B/B~MSY~ of 0.4,1, or 1.6."} ram_mape_map_plot ``` -```{r acc-map, fig.cap="Mean classification accuracy (assignment to FAO stock status category) by FAO statistical area arising from different data sources. ram-data refers to catch and abundance index drawn from RAM. CPUE refers to an index of abundance based on reconstructed effort data. cpue-plus uses CPUE along with Fisheries Management Index (FMI) and/or wept area ratio (SAR) data. For both CPUE series 'nominal' assumes a 0% technology rate, otherwise a 2.6% technology rate is assumed. ram_u_umsy assumes all fisheries in the region share a common U/Umsy series with formally assessed fisheries in the region. fmi uses fmi scores to develop a prior on recent fishing mortality rates, sar does the same but based on swept area ratio. CMSY uses the methods from Froese et al. 2017 [@froese2017]."} +```{r acc-map, fig.cap="Mean classification accuracy (assignment to FAO stock status category) by FAO statistical area arising from different data sources. RAM Index refers to catch and abundance index drawn from RAM. Effective CPUE refers to an index of abundance based on reconstructed effort data. Effective CPUE+ uses CPUE along with Fisheries Management Index (FMI) and/or wept area ratio (SAR) data. For both CPUE series 'nominal' assumes a 0% technology creep, for effective a 2.6% technology creep is assumed. RAM U/U~MSY~ assumes all fisheries in the region share a common U/U~MSY~ series with formally assessed fisheries in the region. FMI uses FMI scores to develop a prior on recent fishing mortality rates, SAR does the same but based on swept area ratio. CMSY uses the methods from Froese et al. 2017 [@froese2017]. Guess assignes a random recent B/B~MSY~ of 0.4,1, or 1.6."} ram_acc_map_plot ``` +\newpage + +The poor global performance of models fit to currently accessible broadly available datasets, relative to models fit to data from RAM, indicates that we must prioritize collection and aggregation of new or underutilized data sources if we are to improve our estimates of global fishery status. What sources of data might provide the greatest value in improving our estimates of global stock status? We used sraplus together with the RAM database to estimate the average reduction in error resulting from having access to different kinds of data (Fig.\@ref(fig:voi-plot)). While having access to complete index of abundance, such as a fishery independent survey, was on average able to reduce error relative to a baseline catch-only heuristic, using only the most recent quarter of the available abundance index actually increased error on average. We may have to wait many years for new surveys to provide substantial improvements in status estimates, or work to expand access to long-running existing surveys that have yet to be fully utilized in fisheries assessment [@maureaud2020]. + +```{r voi-plot, fig.cap="Posterior probability distributions of estimated effect of different data types on root mean squared error of B/B~MSY~ in the most recent 5 years of data available for each model fit. Distribution is full posterior probability distribution. Point is median, thicker black section inner 66th quantile of the posterior, the thinner black line the 95th. Change is relative to the mean performance of a catch-only heuristic model."} +b_voi_plot + + geom_vline(aes(xintercept = 0), alpha = 0.75) +``` + # Discussion @@ -341,7 +368,10 @@ ram_acc_map_plot Global-level assessment are critical for guiding management agendas for the world's oceans, and tracking critical indicators such as the United Nations Sustainable Development Goals. Despite this need, and despite advances in stock assessment methods and available data, we show that our understanding of the world's fisheries remains murky in many parts of the world. While in some cases addition of globally available data, such as quality of fisheries management or effort reconstructions, provided value above and beyond catch histories alone (Fig.\@ref(fig:cs-plot)), at the global level models fit using each of the available datasets, besides the RAM-derived indices, produced biased and imprecise estimates of stock status, frequently performing worse than a simple guess (Table.\@ref(tab:perf-tab)). -What quality of assessment is needed and what constitutes a meaningful improvement in assessment quality depends on the needs of those using the assessment outputs. It may be that for particular regions, species, or uses the results presented here or in other past global analyses are sufficient. In some instances using the data presented here did provide some improvement over use of catch-only style methods; the difficulty comes in attempting to apply data types uniformly across the globe. While it is unreasonable to expect any global-scale data to be able to perform as well as data pulled from RAM assessments themselves, or that data-limited methods would perform well for every individual stock, our hope would be that a data-limited approach based on globally available data sources would be able to correctly capture general trends in stock status in time and space. That none of the datasets collected here can achieve that, and that our test on the RAM data suggest that model misspecification is not the primary culprit, tells us that improvements in estimates of global stock status must come from improvements in the data themselves. +What quality of assessment is needed and what constitutes a meaningful improvement in assessment quality depends on the needs of those using the assessment outputs. It may be that for particular regions, species, or uses the results presented here or in other past global analyses are sufficient. In some instances using the data presented here did provide some improvement over use of catch-only style methods; the difficulty comes in attempting to apply data types uniformly across the globe. While it is unreasonable to expect any global-scale data to be able to perform as well as detailed stock assessments reported in RAM, or that data-limited methods would perform well for every individual stock, our hope would be that a data-limited approach based on globally available data sources would be able to correctly capture general patterns in stock status in time and space. That none of the datasets collected here can achieve that, and that our test on the RAM data suggest that model misspecification is not the primary culprit, tells us that improvements in estimates of global stock status must come from improvements in the quality and use of data themselves. + +We chose to test the performance of methods against estimated values reported in RAM. A reasonable critique of this choice is that unassessed stocks, on which these methods would actually be used, are likely to have vastly different dynamics than the heavily managed stocks in RAM. The clear and systemically poor performance of the methods tested here when applied to RAM stocks would require though that these methods have massively lower bias and higher accuracy for unassessed fisheries than RAM stocks if they are to be expected to provide reasonable picture of global stock status. Free et al. 2020's [@free2020] simulation testing of the types of methods included here suggests that this is unlikely to be true. z + @@ -349,40 +379,35 @@ Our results do not imply that the kinds of broadly available data presented here But, we must simultaneously consider data quality and resolution: applying one SAR value to all stocks in a region, even if that SAR value can provide valuable information for a subset of fisheries, causes inaccurate estimates of stock status when applied too broadly. Our analysis does not show that the data considered here are without value, but that attempting to indiscriminately apply these data to all areas results in meaningfully incorrect estimates of stock status for regions whose nature does not match the assumptions needed to apply these data sources. -What sources of data might provide the greatest value in improving our estimates of global stock status? We used `sraplus` together with the RAM database to estimate the average reduction in error resulting from having access to different kinds of data (Fig.\@ref(fig:voi-plot)). While having access to complete index of abundance, such as a fishery independent survey, was on average able to reduce error relative to a baseline catch-only heuristic, using only the most recent quarter of the available abundance index actually increased error on average. We may have to wait many years for new surveys to provide substantial improvements in status estimates, or work to expand access to long-running existing surveys that have yet to be fully utilized in fisheries assessment [@maureaud2020]. -Our value-of-information analysis also shows though the high utility of having access to even a recent snapshot of F/F~MSY~ (Fig.\@ref(fig:voi-plot)). Swept area ratios, Fisheries Management Index scores, or other similar metrics can be used to construct priors on fishing mortality rates, though care must be taken in applying them at the appropriate spatial resolution. Another avenue would be to prioritize the development of a global repository for length and age composition data. Given appropriate conditions, these length measurements can be used to estimate local fishing mortality rates [@hordyk2016; @prince2019; @rudd2017]. While length-based assessments come with a host of assumptions and pitfalls, properly implemented in some fisheries this may provide an overlooked source of fisheries data at a global scale, at least as an improvement over relying on catch-alone or regional proxies. Such a database could be used to construct stock or stock complex specific priors on fishing mortality for particular regions around the globe, which when paired with catch data and where possible indicies of abundance could meaningfully improve our understanding of global fisheries [@thorson2015f]. +Our value-of-information analysis also shows though the high utility of having access to even a recent snapshot of F/F~MSY~ (Fig.\@ref(fig:voi-plot)). Swept area ratios, Fisheries Management Index scores, or other similar metrics can be used to construct priors on fishing mortality rates, though care must be taken in applying them at the appropriate spatial resolution. Another avenue would be to prioritize the development of a global repository for length and age composition data. Given appropriate conditions, these length measurements can be used to estimate local fishing mortality rates [@hordyk2016; @prince2019; @rudd2017]. While length-based assessments come with a host of assumptions and pitfalls, properly implemented in some fisheries with appropriate life histories these methods may provide an overlooked source of information on fisheries at a global scale, at least as an improvement over relying on catch-only or regional proxies alone. Such a database could be used to construct stock or stock complex specific priors on fishing mortality for particular regions around the globe, which could meaningfully improve our understanding of global fisheries, particularly when paired with catch data and where possible indices of abundance [@thorson2015f;@rudd2017]. -```{r voi-plot, fig.cap="Posterior probability distributions of estimated effect of different data types on root mean squared error of B/B~MSY~ in the most recent 5 years of data available for each model fit. Distribution is full posterior probability distribution. Point is median, thicker black section inner 66th quantile of the posterior, the thinner black line the 95th. Change is relative to the mean performance of a catch-only heuristic model."} -b_voi_plot + - geom_vline(aes(xintercept = 0), alpha = 0.75) -``` We must also prioritize collection and curation of fish population survey data worldwide. Repositories of fishery-independent survey data would be immensely beneficial, such as those maintained by [FishStat](https://james-thorson.github.io//). Recent research confirms that there are bottom trawl data to support analysis of biomass-trends since 2001 and potentially earlier in many regions [@maureaud2020], and survey data are available for more stocks than have previously had stock assessments. Effort reconstructions such as those utilized here may help create fishery-dependent abundance indices in some instances, and going forward datasets such as those compiled by [Global Fishing Watch](https://globalfishingwatch.org/) in combination with the reconstruction approaches of [@rousseau2019] might allow us to construct and use timeseries of fishing effort specific to particular areas, fleets, and species complexes -Expanded training of fisheries scientists around the globe is another critical need. Even were we to dramatically expand the amount and types of data available for global assessment, individual fisheries and regions will need to make informed decisions about which sources of data may be applicable and which not, and to critically evaluate the results of any model based on local expertise. This is why stock assessments even in the data-rich world are not an automated process; the real challenge is often not in fitting a model to data but in understanding how best to use the data and the quality and limitations of the model used. Empowering a global network of fisheries scientists through training and peer-support would both improve the health and management of local stocks and provide a means for ensuring that global estimates of stock status are as accurate as they can be. +Expanded training of fisheries scientists around the globe is another critical need. Even were we to dramatically expand the amount and types of data available for global assessment, individual fisheries and regions will need to make informed decisions about which sources of data may be applicable and which not, and to critically evaluate the results of any model based on local expertise. This is why stock assessments even in the data-rich world are not an automated process; the real challenge is often not in fitting a model to data but in understanding how best to use the data and the quality and limitations of the model used. Empowering a global network of fisheries scientists through training and peer-support would help local experts make the most of available data, ensure the reliability of newly collected data, and improve the interpretation of assessment results. The coming decades are a critical time for the future of fisheries and ocean health. Achieving the United Nations Sustainable Development Goal 14 for the conservation and sustainable use of the world's oceans depends on our ability to effectively assess the status of fish stocks around the world. The RAM Legacy Stock Assessment Database combined with the FAO's expert elicitation of status for select stocks have dramatically improved our understanding of global fisheries in recent years. However, this process still leaves a substantial number of fisheries and global catch unassessed. Numerous catch-based data-limited approaches have attempted to fill that gap, and while these efforts have advanced our knowledge and interest in unassessed fisheries, none have yet been able to provide a clear solution to this problem. Improving estimates of global stock status depends on investing in an improved and expanded global network of fishery data and fisheries scientists. # Methods -All analysis were conducted in the R programming language [@rcoreteam2019]. Model fitting was conducted using Rcpp [eddelbuettel2011] and stan [@carpenter2017], implemented through Template Model Builder [@kristensen2016] by the tmbstan package [@monnahan2018]. The `sraplus` package is publicly available at [github.com/danovando/sraplus](https://github.com/DanOvando/sraplus), and all materials needed to fully reproduce this manuscript are available at [github.com/DanOvando/assessing-global-fisheries](https://github.com/DanOvando/assessing-global-fisheries). +All analysis were conducted in the R programming language [@rcoreteam2019]. Model fitting was conducted using Rcpp [@eddelbuettel2011] and stan [@carpenter2017], implemented through Template Model Builder [@kristensen2016] by the tmbstan package [@monnahan2018]. The sraplus package is publicly available at [github.com/danovando/sraplus](https://github.com/DanOvando/sraplus), and all materials needed to fully reproduce this manuscript are available at [github.com/DanOvando/assessing-global-fisheries](https://github.com/DanOvando/assessing-global-fisheries). -`sraplus` contains too many options to cover within this methods section. We encourage readers to explore the documentation available at the package website at www.github.com/danovando/sraplus. Below we describe the structure of the population model underpinning `sraplus`, the estimation models used, and the construction of priors used in this paper. +sraplus contains too many options to cover within this methods section. We encourage readers to explore the documentation available at the package website at www.github.com/danovando/sraplus. Below we describe the structure of the population model underpinning sraplus, the estimation models used, and the construction of priors used in this paper. ## Population Model -The core of `sraplus` is a Pella-Tomlinson [@pella1969] production model constructed in the manner of [@winker2018]. While models of these kinds abstract away many important details of fish biology and fleet behavior, they are the highest resolution model that the potential data evaluated here will support. The purpose of `sraplus` is not to make substantial improvements in the fitting of surplus production models, but to provide a flexible tool for exploring the impacts of adding different kinds of data and priors on estimates of fishery status +The core of sraplus is a Pella-Tomlinson [@pella1969] production model constructed in the manner of [@winker2018]. While models of these kinds abstract away many important details of fish biology and fleet behavior, they are the highest resolution model that the potential data evaluated here will support. The purpose of sraplus is not to make substantial improvements in the fitting of surplus production models, but to provide a flexible tool for exploring the impacts of adding different kinds of data and priors on estimates of fishery status The population growth equation is \begin{equation} f(x)=\begin{cases} -B_{t + 1} = \left(B_{t} + B_{t}\frac{r}{m - 1}\left(1 - \left(\frac{B_t}{K}\right)^{m- 1}\right) - \hat{c_t}\right)p_t +B_{t + 1} = \left(B_{t} + B_{t}\frac{r}{m - 1}\left(1 - \left(\frac{B_t}{K}\right)^{m- 1}\right) - c_t\right)p_t , & \text{if $B_t>0.25 \times K$}.\\ -B_{t + 1} = \left(B_{t} + \frac{B_{t}}{0.25 \times K}\left(B_{t}\frac{r}{m - 1}\left(1 - \left(\frac{B_t}{K}\right)^{m- 1}\right) - \hat{c_t}\right)\right)p_t, & \text{otherwise}. +B_{t + 1} = \left(B_{t} + \frac{B_{t}}{0.25 \times K}\left(B_{t}\frac{r}{m - 1}\left(1 - \left(\frac{B_t}{K}\right)^{m- 1}\right) - c_t\right)\right)p_t, & \text{otherwise}. \end{cases} (\#eq:pt) \end{equation} @@ -401,15 +426,15 @@ Process error *p* is assumed to be log-normally distributed, such that ## Estimation Model -All of our estimates are Bayesian in nature. We can break the use of `sraplus` into two distinct categories: with data and without. By "data", we refer to measurements which are used to confront model estimates within a likelihood function. In our context, these include fishery-independent survey data, or a CPUE index. When there are no data, the model amounts to filtering priors through the model (the combination of the Pella-Tomlinson operating model and the catches for the stock in question, along with any fixed parameters). Under this mode, the model is essentially a stock-reduction analysis model, in the manner of [@walters2006], in which we ask, which combinations of prior probability distributions of parameters do not crash the population (i.e. results in biomass less than catches in any time step in the fishery's history), given the constraints of the population model and the catches. This step updates the prior distribution of population parameters by eliminating combinations of priors that are impossible for a given catch history and a specified functional form. +All of our estimates are Bayesian in nature. We can break the use of sraplus into two distinct categories: with data and without. By "data", we refer to measurements which are used to confront model estimates within a likelihood function. In our context, these include fishery-independent survey data, or a CPUE index. When there are no data, the model amounts to filtering priors through the model (the combination of the Pella-Tomlinson operating model and the catches for the stock in question, along with any fixed parameters). Under this mode, the model is essentially a stock-reduction analysis model, in the manner of [@walters2006], in which we ask, which combinations of prior probability distributions of parameters do not crash the population (i.e. results in biomass less than catches in any time step in the fishery's history), given the constraints of the population model and the catches. This step updates the prior distribution of population parameters by eliminating combinations of priors that are impossible for a given catch history and a specified functional form. The full list of estimable parameters are listed in Table.\@ref(tab:param-tab). *r* and *K* are the only two parameters that are always estimated. Estimation of every other parameter can be turned on or off. When estimation is turned off estimable parameters are fixed at their initial values, which can either be set to model defaults or specified by the user. For our main sets of results (everything excluding the value of information analysis), the estimated parameters are *r*, *K*, $\sigma_{obs}$, $\gamma$, and *B0*. *q* is also estimated when needed. ```{r param-tab} -cd <- data_frame(Parameter = "Carrying Capacity", Abbreviation = "K", `Default Prior` = "Prior predictive tuning") %>% - rbind(c("Growth rate", "r", "Thorson, 2020 [@thorson2020] updated by prior predictive tuning")) %>% - rbind(c("Shape parameter", "m", "Drawn from Thorson et al. 2012 (@thorson2012)")) %>% - rbind(c("Catchability", "q", "$logn(1e^{-3}, 1)$")) %>% +cd <- data_frame(Parameter = "Carrying Capacity", Abbreviation = "$K$", `Default Prior` = "Prior predictive tuning") %>% + rbind(c("Growth rate", "$r$", "Thorson, 2020 [@thorson2020] updated by prior predictive tuning")) %>% + rbind(c("Shape parameter", "$m$", "Drawn from Thorson et al. 2012 (@thorson2012)")) %>% + rbind(c("Catchability", "$q$", "$logn(1e^{-3}, 1)$")) %>% rbind(c("Observation Error", "$\\sigma_{obs}$", "$logn(.05,1)$")) %>% rbind(c("Ratio of process to observation error", "$\\gamma$", "$logn(.5,0.25)$")) %>% rbind(c("Initial State", "$B0$", "Posterior probability dist. of catch-based regressions")) @@ -427,29 +452,31 @@ knitr::kable( kableExtra::row_spec(0, bold = T) ``` -`sraplus` can be run in two forms: either as a stock reduction analysis [SRA, @walters2006], or fit to an index of abundance (fishery dependent or independent). Unless there is an abundance index to fit to the model runs as a stock reduction analysis. A stock reduction analysis works by specifying prior distributions on population parameters and critically the recent state of the fishery. `sraplus` allows users to specify the most recent status in units of depletion, B/B~MSY~, F, or F/F~MSY~. We then sample from the prior distributions of the population model parameters and apply those to the production model, along with the catch history. Any run that results in the collapse of the population (catch greater than biomass in any time step) is immediately rejected. The remaining viable draws from the prior distributions are sampled in proportion to the supplied prior on recent stock status. Any model results based on data sources listed in Table.\@ref(tab:dat-desc) that do not contain "cpue" or "RAM data" are estimated through stock reduction analysis. All SRA style runs in our paper sampled 2,000 draws of the prior-predictive distribution from a total of 1e6 candidate draws. +sraplus can be run in two forms: either as a stock reduction analysis [SRA, @walters2006], or fit to an index of abundance (fishery dependent or independent). Unless there is an abundance index to fit to the model runs as a stock reduction analysis. A stock reduction analysis works by specifying prior distributions on population parameters and critically the recent state of the fishery. sraplus allows users to specify the most recent status in units of depletion, B/B~MSY~, F, or F/F~MSY~. We then sample from the prior distributions of the population model parameters and apply those to the production model, along with the catch history. Any run that results in the collapse of the population (catch greater than biomass in any time step) is immediately rejected. The remaining viable draws from the prior distributions are sampled in proportion to the supplied prior on recent stock status. Any model results based on data sources listed in Table.\@ref(tab:dat-desc) that do not contain "cpue" or "RAM data" are estimated through stock reduction analysis. All SRA style runs in our paper sampled 2,000 draws of the prior-predictive distribution from a total of 1e6 candidate draws. When an index of abundance is available the model estimates the posterior probability distributions of the estimated and transformed parameters using Hamiltonian Monte Carlo implemented in stan [@standevelopmentteam2018] accessed through the `tmbstan` interface [@monnahan2018]. By default the model uses 2000 draws with a 1000 step warm-up and one chain. Any detailed fit for a particular fishery would likely use more draws and chains, but we verified that this sampling routine produced an acceptable tradeoff of speed and convergence criteria. The model fits to a direct estimate of abundance (e.g. a fishery independent survey or a standardized catch-per-unit-effort index), the likelihood is $$log(a_t) \sim normal(f(r,K,m,B0,\pmb{p},\pmb{c}) \times q, sigma_{obs})$$ -where *f* is the Pella-Tomlinson production model (Equation.\@ref(eq:pt)). When an effort index is available, `sraplus` constructs an index of abundance based on the catch and effort data. [@rousseau2019] measure an index of abundance as catch divided by their effort index, either nominal or effective (assuming the 2.6% annualized technology rate). This rate of technology creep assumes that every unit increase in effort is log-linearly greater than the unit of effort before it. When effort increases dramatically above historic levels, this can create a CPUE index that decreases faster than the true population. This is due to the fact that the marginal fishing mortality produced by increasing unit of efforts increases decreases as effort approaches infinity (since fishing mortality is bounded between 0 and 1). To accommodate this, we generate a catch per effective harvest rate index of abundance, as +where *f* is the Pella-Tomlinson production model (Equation.\@ref(eq:pt)). When an effort index is available, sraplus constructs an index of abundance based on the catch and effort data. [@rousseau2019] measure an index of abundance as catch divided by their effort index, either nominal or effective (assuming the 2.6% annualized technology creep). This rate of technology creep assumes that every unit increase in effort is log-linearly greater than the unit of effort before it. When effort increases dramatically above historic levels, this can create a CPUE index that decreases faster than the true population. This is due to the fact that the marginal fishing mortality produced by increasing unit of efforts increases decreases as effort approaches infinity (since fishing mortality is bounded between 0 and 1). To accommodate this, we generate a catch per effective harvest rate index of abundance, as -$$cpue_t = \frac{catch_t}{(1 - e^{-f_t})}$$ +$$cpue_t = \frac{c_t}{(1 - e^{-f_t})}$$ $$f_t = q_tE_t$$ -Where $q_t$ can has a technology rate component $\tau$ +Where $q_t$ can has a technology creep component $\tau$ $$q_t = q_{t-1} \times (1 + \tau)$$ We then fit to the index of abundance per -$$log(cpue_t) \sim normal(f(r,K,m,B0,\pmb{p},\pmb{C}), \sigma_{obs})$$ +$$log(cpue_t) \sim normal(f(r,K,m,B0,\pmb{p},\pmb{c}), \sigma_{obs})$$ ### CMSY -In addition to the results from `sraplus`, we include a set of results produced by the CMSY method [@froese2017]. For computational efficiency, we used a ported version of the CMSY model available at . The only modification made is to convert the underlying population model to C++ for faster computation. For each stock we used all the default options provided by CMSY, except for resilience, which was pulled from the vulnerability scores from FishBase accessed through `rfishbase` [boettiger2012]. Vulnerability scores greater than 66 were scored as low resilience, between 33 and 66 medium resilience, and lower than 33 high resilience. +In addition to the results from sraplus, we include a set of results produced by the CMSY method [@froese2017]. For computational efficiency, we used a ported version of the CMSY model available at . The only modification made is to convert the underlying population model to C++ for faster computation. For each stock we used all the default options provided by CMSY, except for resilience, which was pulled from the vulnerability scores from FishBase accessed through `rfishbase` [@boettiger2012]. Vulnerability scores greater than 66 were scored as low resilience, between 33 and 66 medium resilience, and lower than 33 high resilience. + +In the absence of an index of abundance CMSY and sraplus are conceptually the same: Both sample from prior distributions of life history and stock status conditional on the stock not collapsing given a catch history and a functional form for the population growth equation. The core different between CMSY and sraplus as utilized here then is that when we refer to CMSY we mean CMSY run using the default heuristics used by the algorithm to construct priors on initial, intermediate, and terminal stock status. In contrast, sraplus fit using for example FMI data uses the empirical models described here to generate priors on stock status. In addition, all sraplus runs use a prior-predictive tuning procedure to resolve issues related to Borel's paradox, described below. ## Priors @@ -459,9 +486,9 @@ We address two critical features of prior use in sraplus below: tuning of the pr ### Prior Predictive Tuning -Suppose that the only thing we observe from a fishery is a catch history. Assuming Pella-Tomlinson population dynamics, the only thing we can learn from this catch history alone is the set of model parameters that ensure that the population still exists and never collapsed in the past. We can think of this as a binomial process in which we fit a model of the population, conditional on catches, to the fact that we know that population existed in each time step of the catch history. Beyond that though, baring additional information or model assumptions we have no way of knowing whether these catches represent a substantial proportion of a small population or a minuscule fraction of a massive population; all we know is that the current population must be greater than 0. +Suppose that the only thing we observe from a fishery is a catch history. Assuming Pella-Tomlinson population dynamics, the only thing we can learn from this catch history alone is the set of model parameters that ensure that the population still exists and never collapsed in the past. We can think of this as a binomial process in which we fit a model of the population, conditional on catches, to the fact that we know that population existed in each time step of the catch history. Beyond that though, baring additional information or model assumptions we have no way of knowing whether these catches represent a substantial proportion of a small population or a minuscule fraction of a massive population; all we know is that the current population must be greater than zero. -In the absence of any data to fit to, the SRA algorithm works by assuming that we know current stock status, and then finds feasible parameters to satisfy that belief. This creates a problem for the Bayesian nature of our analysis though. Consider a production model with two parameters, a growth rate *r* and a carrying capacity *K*. Once we specify prior distributions on *r* and *K*, and then apply these distributions to our model (the shape of the production function along with the catch histories), we have implicitly provided a prior on the status of the stock in all time periods, since each unique combination of *r* and *K* together with the model and the catch history produces a deterministic stock status in each time step. Doing so places essentially two priors on recent stock status: one implicit prior through the population parameter priors, and one explicit through the users perception of recent stock status, creating a problem termed Borel's Paradox. This may seem like an academic concern, and indeed in our experience when the data are sufficiently informative the Bayesian version of our model subject to Borel's paradox produces effectively identical results to those produce by the same model fit by maximum likelihood. However, particularly for the SRA version of `sraplus`, Borel's Paradox poses a particular problem. +In the absence of any data to fit to, the SRA algorithm works by assuming that we know current stock status, and then finds feasible parameters to satisfy that belief. This creates a problem for the Bayesian nature of our analysis though. Consider a production model with two parameters, a growth rate *r* and a carrying capacity *K*. Once we specify prior distributions on *r* and *K*, and then apply these distributions to our model (the shape of the production function along with the catch histories), we have implicitly provided a prior on the status of the stock in all time periods, since each unique combination of *r* and *K* together with the model and the catch history produces a deterministic stock status in each time step. Doing so places essentially two priors on recent stock status: one implicit prior through the population parameter priors, and one explicit through the users perception of recent stock status, creating a problem termed Borel's Paradox (See Poole and Raftery 2000 [@poole2000] and references therein for a discussion of Borel's Paradox in a fisheries context). This may seem like an academic concern, and indeed in our experience when the data are sufficiently informative the Bayesian version of our model subject to Borel's paradox produces effectively identical results to those produce by the same model fit by maximum likelihood. However, particularly for the SRA version of sraplus, Borel's Paradox poses a particular problem. The SRA algorithm works in two steps. First, the algorithm rejects any draws that resulted in the collapse of the population (biomass less than catch in a given timestep). From there a standard SRA would sample from the priors in proportion to the stated prior on recent stock status. If the bulk of the prior on terminal stock status was concentrated at 50% of *K*, combinations of *r* and *K* that produce terminal stock status near 50% of *K* are sampled proportionally more frequently. However, lower values of terminal stock status have fewer candidate values of *r* and *K*, since it becomes harder and harder to find viable pairs that come close to but do not crash the population at any time step. Conversely, in the absence of constraints higher values of stock status have infinite combinations of plausible *r* and *K* combinations: since under this model the population cannot be greater than carrying capacity, as for example *K* approaches infinity terminal stock status asymptotes at close to 100% of *K*. The net result of this is that even though individual combinations of *r* and *K* that produce higher stock status than the mean of the prior on recent stock status individually have lower probability of being sampled, there are many more opportunities for the lower-probability events that produce higher stock status to be sampled. As a result, the post-model-pre-data prior on terminal depletion will always be higher under this method than the supplied prior on stock status. @@ -533,9 +560,9 @@ comps %>% ### Priors Informed by Outside Data -Along with allowing users to supply their own priors, the `sraplus` package contains three built-in methods used in this manuscript for converting information on stock status from additional outside data into a form usable as a prior by `sraplus`. We paired data on catch histories, swept area ratio (SAR), and fisheries management index (FMI) with estimates of stock status from the RAM legacy stock assessment database. We then trained a model of the general form $log(status) \sim normal(variable,\sigma)$ for each of these three data types. Given values of these variables for a new fishery then, `sraplus` uses the fitted model to generate posterior predictive distributions of stock status based on these data, which can then be used as priors on stock status by `sraplus`. +Along with allowing users to supply their own priors, the sraplus package contains three built-in methods used in this manuscript for converting information on stock status from additional outside data into a form usable as a prior by sraplus. We paired data on catch histories, swept area ratio (SAR), and fisheries management index (FMI) with estimates of stock status from the RAM legacy stock assessment database. We then trained a model of the general form $log(status) \sim normal(variable,\sigma)$ for each of these three data types. Given values of these variables for a new fishery then, sraplus uses the fitted model to generate posterior predictive distributions of stock status based on these data, which can then be used as priors on stock status by sraplus. -All prior regression models where tested by out-of-sample predictive power, and where competing models were considered the final model was chosen by leave-on-out validation [@vehtari2017]. The final models are intended as a reasonably robust means of translating available data (catch histories, FMI, and SAR values) into a form usable by `sraplus`. Given the scope of this analysis, we do not claim that the presented regressions are the best possible model relating these data with the fishery status indicators of interest. Rather, each regression was tested to ensure that it is unlikely that, given the same data, an alternative model would perform substantially better than those presented here. +All prior regression models where tested by out-of-sample predictive power, and where competing models were considered the final model was chosen by leave-on-out validation [@vehtari2017]. The final models are intended as a reasonably robust means of translating available data (catch histories, FMI, and SAR values) into a form usable by sraplus. Given the scope of this analysis, we do not claim that the presented regressions are the best possible model relating these data with the fishery status indicators of interest. Rather, each regression was tested to ensure that it is unlikely that, given the same data, an alternative model would perform substantially better than those presented here. #### Catch-Only Priors @@ -572,9 +599,9 @@ $$log(value_i) \sim normal(s(SAR_i) + s(\frac{catch_i}{max(catch)_i)}) + 1,\sigm ## Value of Information Calculations -What sources of data might provide the greatest value in improving our estimates of global stock status? We fit 3,000 `sraplus` models to randomly sampled fisheries from RAM, each time varying the kind and quality of data made available to the model, and what parameters the model estimated. We then calculated the root-mean-squared-error (RMSE) between the observed and predicted B/B~MSY~ over the most recent five years of each fit, and then fit a regression to these data to estimate the posterior probability of the effect of different data types and model states on RMSE (Fig.\@ref(fig:voi-plot)). +What sources of data might provide the greatest value in improving our estimates of global stock status? We fit 3,000 sraplus models to randomly sampled fisheries from RAM, each time varying the kind and quality of data made available to the model, and what parameters the model estimated. We then calculated the root-mean-squared-error (RMSE) between the observed and predicted B/B~MSY~ over the most recent five years of each fit, and then fit a regression to these data to estimate the posterior probability of the effect of different data types and model states on RMSE (Fig.\@ref(fig:voi-plot)). -The value-of-information (VOI) calculations presented in Fig.\@ref(fig:voi-plot) help illustrate what types of data may be most beneficial to acquire at a global scale if we are to improve our knowledge of the state of global fisheries. The VOI analysis was performed by using `sraplus` to generate estimates of stock status (B/B~MSY~) for stocks in the RAM legacy stock assessment, and comparing the estimated values to the values reported in RAM. We generate fits for 3000 combinations of a RAM stock and available data. For any one draw, we randomly sample a RAM stock and a list of available data and data quality. For example, we might sample stock *A* with information on recent fishing mortality rates for the first iteration, and stock *A* again for the second iteration but now with information on recent fishing mortality rates and a recent index of abundance. The result is a set of model performance estimates where the characteristics of the stock and the data made available to the model are randomized. +The value-of-information (VOI) calculations presented in Fig.\@ref(fig:voi-plot) help illustrate what types of data may be most beneficial to acquire at a global scale if we are to improve our knowledge of the state of global fisheries. The VOI analysis was performed by using sraplus to generate estimates of stock status (B/B~MSY~) for stocks in the RAM legacy stock assessment, and comparing the estimated values to the values reported in RAM. We generate fits for 3000 combinations of a RAM stock and available data. For any one draw, we randomly sample a RAM stock and a list of available data and data quality. For example, we might sample stock *A* with information on recent fishing mortality rates for the first iteration, and stock *A* again for the second iteration but now with information on recent fishing mortality rates and a recent index of abundance. The result is a set of model performance estimates where the characteristics of the stock and the data made available to the model are randomized. Using this set of fits, we assess performance as the root-mean-squared-error of B/B~MSY~ over the most recent 5 years of the fishery, in order to evaluate the ability of the model to capture the recent trends in stock status and not just the most recent year. We evaluate the contributing of each data type to RMSE using a Gamma GLM with a log link of the form @@ -587,19 +614,30 @@ Where $\pmb{\beta}$ is the vector of coefficients associated with the matrix of ::: {#refs} ::: -# Supplemental Material +# Supplementary Information + +\renewcommand{\thefigure}{S\arabic{figure}} + +\setcounter{figure}{0} + + +\renewcommand{\thetable}{S\arabic{table}} + +\setcounter{table}{0} + + +## Data Sources Used ```{r dat-desc} dat_description <- tribble(~`Data Name`, ~Description, - "RAM-Data", "Fit to abundance index from RAM", - "sar", "Prior on terminal F/Fmsy set by regional swept area ratio", - "fmi", "Prior on terminal F/Fmsy set by regional fisheries management index scores", - "cpue", "Fit to CPUE index created from RAM catch and regional effort index. 2.5% tech. creep", - "cpue-plus", "Fit to CPUE index created from RAM catch and regional effort index with priors informed by SAR and FMI. 2.5% tech. creep", - "nominal-cpue", "Fit to CPUE index created from RAM catch and regional effort index. 0% tech. creep", - "nominal-cpue-plus", "Fit to CPUE index created from RAM catch and regional effort index with priors informed by SAR and FMI. 0% tech. creep", - "catch-only", "priors on terminal B/Bmsy informed by catch history", - "heuristic", "priors on initial and terminal B/Bmsy informed by CMSY heuristc") + "RAM Index", "Fit to abundance index from RAM", + "SAR", "Prior on terminal F/Fmsy set by regional swept area ratio", + "FMI", "Prior on terminal F/Fmsy set by regional fisheries management index scores", + "EFfective CPUE", "Fit to CPUE index created from RAM catch and regional effort index. 2.6% technology creep", + "Effective CPUE+", "Fit to CPUE index created from RAM catch and regional effort index with priors informed by SAR and FMI. 2.6% tech. creep", + "Nominal CPUE", "Fit to CPUE index created from RAM catch and regional effort index. 0% tech. creep", + "Nominal CPUE+", "Fit to CPUE index created from RAM catch and regional effort index with priors informed by SAR and FMI. 0% tech. creep", + "Guess", "Priors on terminal B/Bmsy randomly sampled from .4,1,1.6") knitr::kable(dat_description, caption = "Data sources used for terminal stock status estimate", @@ -609,8 +647,10 @@ knitr::kable(dat_description, # kableExtra::row_spec(0, bold = T) ``` +## Sample Prior Posterior Plots + -```{r ppplot, fig.cap="Prior-posterior plots of fits for case study fishery in Fig.6"} +```{r ppplot, fig.cap="Prior posterior plots of fits for case study fishery in Fig.6"} a = get_prior_posterior(co_fit, co_driors) @@ -626,3 +666,123 @@ a$prior_posterior %>% ``` + +## Population Model + +The core of our model is a Pella-Tomlinson [@pella1969] production model in the manner of @winker2018. While models of these kinds abstract away many important details of fish biology and fleet behavior, they are the highest resolution model that the potential data evaluated here will support. + +```{r si-prior-tab} +cd <- data_frame(Parameter = "Carrying Capacity", Abbreviation = "K", `Default Prior` = "$logn(10\\times{max(catch)},4)$") %>% + rbind(c("Growth rate", "r", "$logn(r_{fishlife}, \\sigma_{r,fishlife})$")) %>% + rbind(c("Shape parameter", "m", "$logn(1.01, 0.25)$")) %>% + rbind(c("Catchability", "q", "$logn(1e^{-3}, 0.3)$")) %>% + rbind(c("Observation Error", "$\\sigma_{obs}$", "$logn(.05,2)$")) %>% + rbind(c("Process Error", "$\\sigma_{proc}$", "$logn(.05,0.5)$")) %>% + rbind(c("Tech Creep", "$creep$", "$logn(0.025,0.2)$")) + + +knitr::kable( + cd, + "latex", + align = "l", + booktabs = TRUE, + escape = F, + row.names = F, + caption = "Potential parameters included in sraplus, abbreviations, and default priors" +) %>% + kableExtra::kable_styling(full_width = T, font_size = 9) %>% + kableExtra::row_spec(0, bold = T) +``` + + +The population growth equation is + +\begin{equation} + f(x)=\begin{cases} + B_{t + 1} = \left(B_{t} + B_{t}\frac{r}{m - 1}\left(1 - \left(\frac{B_t}{K}\right)^{m- 1}\right) - \hat{c_t}\right)p_t +, & \text{if $B_t>0.25 \times K$}.\\ + B_{t + 1} = \left(B_{t} + \frac{B_{t}}{0.25 \times K}\left(B_{t}\frac{r}{m - 1}\left(1 - \left(\frac{B_t}{K}\right)^{m- 1}\right) - \hat{c_t}\right)\right)p_t, & \text{otherwise}. + \end{cases} + (\#eq:sihockey) +\end{equation} + +## Prior Generating Regressions + +### Swept Area Ratio + +```{r sar-fit, fig.cap = "Observed (x-axis) vs posterior predictive (y-axis) F/F~MSY~ for regression of swept area ratio (SAR) on F/F~MSY~"} +fit <- sraplus::sar_models$fit[sraplus::sar_models$metric == "mean_uumsy"][[1]] + +fit_data <- fit$data + +pp <- posterior_predict(fit, draws = 1000) + +r2_plot <- tibble(r2 = bayes_R2(fit)) %>% + ggplot(aes(r2)) + + geom_histogram() + + labs(subtitle = "B) Bayesian R2") + +pi_plot <- ppc_intervals(y = fit_data$log_value, yrep = pp, x = fit_data$log_value) + + scale_x_continuous(name = "Observed Log F/Fmsy") + + scale_y_continuous(name = "Posterior Predictive F/Fmsy") + + labs(subtitle = "A) Observed vs. Posterior Predictive") + + pub_theme + +pi_plot + r2_plot + +``` + +### Fisheries Management Index + +```{r fmi-fit, fig.cap = "Observed (x-axis) vs posterior predictive (y-axis) F/F~MSY~ for regression of fisheries management index (FMI) on F/F~MSY~"} + +fit <- sraplus::fmi_models$fit[sraplus::fmi_models$metric == "mean_uumsy"][[1]] + +fit_data <- fit$data + +pp <- posterior_predict(fit, draws = 1000) + +r2_plot <- tibble(r2 = bayes_R2(fit)) %>% + ggplot(aes(r2)) + + geom_histogram() + + labs(subtitle = "B) Bayesian R2") + +pi_plot <- ppc_intervals(y = fit_data$log_value, yrep = pp, x = fit_data$log_value) + + scale_x_continuous(name = "Observed Log F/Fmsy") + + scale_y_continuous(name = "Posterior Predictive F/Fmsy") + + labs(subtitle = "A) Observed vs. Posterior Predictive") + + pub_theme + +pi_plot + r2_plot +``` + + + +## Catch History + +```{r catch-reg,fig.cap = "Observed (x-axis) vs posterior predictive (y-axis) B/B~MSY~ for regression of catch on B/B~MSY~"} + +fit <- sraplus::catch_b_model + +fit_data <- fit$data %>% + group_by(stockid) %>% + filter(year == max(year)) + +pp <- posterior_predict(fit, draws = 1000, newdata = fit_data) + +r2_plot <- tibble(r2 = bayes_R2(fit)) %>% + ggplot(aes(r2)) + + geom_histogram() + + labs(subtitle = "B) Bayesian R2") + +pi_plot <- ppc_scatter_avg(y = fit_data$log_value, yrep = pp) + + coord_flip() + + scale_y_continuous(name = "Observed Log B/Bmsy") + + scale_x_continuous(name = "Posterior Predictive B/Bmsy") + + labs(subtitle = "A) Observed vs. Posterior Predictive") + + pub_theme + + +pi_plot + r2_plot +``` + diff --git a/documents/si-ovando-etal-assessing-global-fisheries.Rmd b/documents/si-ovando-etal-assessing-global-fisheries.Rmd index bdd5e47..c59ec8a 100644 --- a/documents/si-ovando-etal-assessing-global-fisheries.Rmd +++ b/documents/si-ovando-etal-assessing-global-fisheries.Rmd @@ -17,7 +17,7 @@ output: bookdown::word_document2: reference_docx: word-template.docx params: - results_name: ["v0.5"] + results_name: ["v1.0"] min_years_catch: [25] abstract: | Assessments of the global state of fisheries play an important role in shaping the public narrative around ocean health, motivating future directions of research and funding, formulating evidence-based policy and tracking the implementation of the United Nations Sustainable Development Goals. While we have reliable estimates of stock status for fisheries accounting for 50% of global catch, our knowledge of the state of the remaining 50%, the worlds 'unassessed' fisheries, is poor. Numerous high-profile publications featuring a range of statistical methods have produced estimates of the global status of these unassessed fisheries, but limited quantity and quality of data along with methdological differences have produced counterintuitive and conflicting results. This is especially true in areas such as Southeast Asia and Africa which are also regions that have high dependence on fishery resources. How can we effectively estimate the state of global fishery sustainability and track progress towards the Sustainable Development Goals targets? We developed a flexible assessment model to assess the value of different kinds, quantities, and quality of data in improving estimates of fishery stock status. We then explore avenues for obtaining potentially impactful data, including through the use of local expert opinion through Fisheries Management Index scores, and increasingly available but historically underutilized data such as trawl footprints and effort data. These data are then used to illustrate how different types of information paint starkly different pictures of the state of fisheries around the world, and to identify priority data types for future collection. Our results provide further evidence that catch data by themselves cannot provide reliable estimates of fishery health. Our ability to obtain reliable estimates of stock status for the world's unassessed fisheries depends on prioritizing the collection of key data at a global scale, not on the development of new modeling methods alone. diff --git a/make-assessing-global-fisheries.R b/make-assessing-global-fisheries.R index 81ebe32..cb7e5bc 100644 --- a/make-assessing-global-fisheries.R +++ b/make-assessing-global-fisheries.R @@ -52,7 +52,7 @@ n_cores <- 2 # options(mc.cores = 1) -results_name <- "v0.5" +results_name <- "v1.0" results_description <- @@ -894,13 +894,13 @@ tmp2 <- tmp2 %>% ) ex_scatter_plot <- tmp2 %>% - ggplot(aes(truth, value, size = lifetime_catch)) + + ggplot(aes(truth, value)) + geom_vline(aes(xintercept = 0)) + geom_hline(aes(yintercept = 0)) + geom_abline(aes(slope = 1, intercept = 0),linetype = 2) + geom_point(alpha = 0.75) + - facet_grid(variable ~ data, scales = "free_y") + - scale_size(trans = "sqrt", name = "Lifetime Catch") + + facet_grid(variable ~ data, scales = "free_x") + + # scale_size(trans = "sqrt", name = "Lifetime Catch") + scale_x_continuous(name = "RAM Value", expand = expansion(add = c(0, .1))) + scale_y_continuous("Estimated Value", expand = expansion(add = c(0, .1))) @@ -1938,7 +1938,7 @@ assess_ram_fits <- ram_status_fits %>% null_model <- assess_ram_fits %>% filter(data == "cmsy") %>% mutate(mean = sample(c(.4,1,1.6), n(), replace = TRUE)) %>% - mutate(data = "Guess") %>% + mutate(data = "guess") %>% mutate(resid = ram_b_v_bmsy - mean, ae = abs(resid)) @@ -2028,6 +2028,19 @@ fao_area_ram_status %>% arrange(mape) +ram_labeller <- c( + "ram-data" = "RAM Index", + "sar" = "SAR", + "fmi" = "FMI", + "cpue" = "Effective CPUE", + "cpue-plus" = "Effective CPUE+", + "nominal-cpue" = "Nominal CPUE", + "nominal-cpue-plus" = "Nominal CPUE+", + "cmsy" = "CMSY", + "guess" = "Guess", + "u_umsy" = "RAM U/Umsy" +) + ram_mpe_map_plot <- fao_area_ram_status %>% filter(!is.na(data)) %>% mutate(data = fct_reorder(data, abs(mpe), .fun = median)) %>% @@ -2039,7 +2052,7 @@ ram_mpe_map_plot <- fao_area_ram_status %>% color = "black", size = 0.01 ) + - facet_wrap(~ data) + + facet_wrap(~ data, labeller = labeller(data = ram_labeller)) + scale_fill_gradient2( low = "steelblue", high = "tomato", @@ -2072,7 +2085,7 @@ ram_mape_map_plot <- fao_area_ram_status %>% color = "black", size = 0.01 ) + - facet_wrap(~ data) + + facet_wrap(~ data, labeller = labeller(data = ram_labeller)) + scale_fill_gradient( low = "white", # mid = "orange", @@ -2106,7 +2119,7 @@ ram_acc_map_plot <- fao_area_ram_status %>% color = "black", size = 0.01 ) + - facet_wrap(~ data) + + facet_wrap(~ data, labeller = labeller(data = ram_labeller)) + scale_fill_viridis( limits = c(0,1), breaks = c(0,.33, .66,1), @@ -2146,7 +2159,7 @@ ram_b_map_plot <- combo %>% color = "black", size = 0.01 ) + - facet_wrap( ~ data) + + facet_wrap(~ data, labeller = labeller(data = ram_labeller)) + scale_fill_viridis( limits = c(0, NA), name = "Median B/Bmsy",