MAB_RiskAssess_2019update.Rmd

---
bibliography: riskassess.bib
csl: plos.csl
fontsize: 10pt
geometry: left=2cm, right=2cm, top=2cm, bottom=3cm, footskip = .5cm
link-citations: yes
output:
  html_document:
    df_print: paged
  pdf_document:
    includes:
      in_header: latex/header.tex
    keep_tex: yes
subparagraph: yes
---


```{r Directory and Data Set-up, echo = F, message = F , warning = F}

#Default Rmd options
knitr::opts_chunk$set(echo = FALSE,
                      message = FALSE,
                      warning = FALSE,
                      fig.align = 'center') #allows for inserting R code into captions

#Plotting and data libraries
library(ggplot2)
library(dplyr)
library(tidyr)
library(ecodata)
library(here)
library(kableExtra)
library(ggrepel)
library(stringr)
library(patchwork)
library(grid)
library(ggiraph)
library(vegan)
library(rpart)

#GIS libraries
library(sf)
library(rgdal)
library(raster)
library(rnaturalearth)

#Data directories
image.dir <- here("images")
gis.dir <- here("gis")

#GIS directory
#gis.dir <- here::here("inst","extdata","gridded")

#General inline text input for report
#Council
council <- "Mid-Atlantic Fishery Management Council"
council_abbr <- "MAFMC"

#Region identifiers
epu <- "Mid-Atlantic Bight"
epu_abbr <- "MAB"
region <- "Mid-Atlantic"
region_abbr <- "MA" #Some commercial data organized by "MA" or "NE" regions, not by EPU 

```

```{r plot constants, echo = F, message = F , warning = F}
#Time series constants
shade.alpha <- 0.3
shade.fill <- "lightgrey"
lwd <- 1
pcex <- 2
trend.alpha <- 0.5
trend.size <- 2
hline.size <- 1
hline.alpha <- 0.35
hline.lty <- "dashed"
label.size <- 5
hjust.label <- 1.5
letter_size <- 4
feeding.guilds <- c("Apex Predator","Piscivore","Planktivore","Benthivore","Benthos")
x.shade.min <- 2009
x.shade.max <- 2018
#Function for custom ggplot facet labels
label <- function(variable,value){
  return(facet_names[value])
}

```


```{r Generic Plot Function, echo = F}
#Plot figure new
soe.plot <- function(data, x.var, y.var, x.label = '', y.label = '', tol = 0.1,
                     x.start = NA, x.end = NA, end.start = 2006, bg.col = background,
                     end.col = recent, trend.neg = main.neg, trend.pos = main.pos,
                     end.trend.neg = end.neg, end.trend.pos = end.pos,

                     stacked = NA, x.line = 2.6, y.line = 3.5, scale.axis = 1,
                     rel.y.num = 1.5, rel.y.text = 1.5){
 
  #Select Data
  x <- data[Var == y.var, ]
  x <- x[order(x[, get(x.var)]), ]
  setnames(x, x.var, 'X')
 
  #Set common time step if necessary
  if(is.na(x.start)) x.start <- min(x[, X])
  if(is.na(x.end))   x.end   <- max(x[, X])
  x <- x[X >= x.start, ]
 
  #Set up plot parameters
  y.max <- max(x[, Value], na.rm = T) + tol * max(x[, Value], na.rm = T)
  y.min <- min(x[, Value], na.rm = T) - tol * abs(min(x[, Value], na.rm = T))
  y.mean <- mean(x[, Value], na.rm = T)
  y.sd <- sd(x[, Value], na.rm = T)
 
  #Plot blank plot
  plot(x[X >= x.start, list(X, Var)], xlim = c(x.start, x.end),
       ylim = c(y.min,y.max), xlab = '', ylab = '', axes = F, ty = 'n')
 
  #Add background
  u <- par('usr')
  rect(u[1], u[3], u[2], u[4], border = NA, col = bg.col)
 
  #Add end period shading
  rect(end.start - 0.5, u[3], u[2], u[4], border = NA, col = end.col)
  abline(h = u[3], lwd=3)
  abline(v = u[1], lwd=3)
 
  #Add grides
  abline(h = y.mean + y.sd, col = 'white', lwd = 3, lty = 2)
  abline(h = y.mean - y.sd, col = 'white', lwd = 3, lty = 2)
 
  #Add data points/lines
  points(x[, list(X, Value)], pch = 16, cex = 1.5)
  lines( x[, list(X, Value)], lwd = 2)
 
  #Add axis
  if(is.na(stacked)) axis(1, cex.axis = 1.5)
  if(!is.na(stacked)){
    if(stacked!= 'A') axis(3, cex.axis = 1.5, tck = 0.1, labels = F)
  }
  #Stacked axes with 0 overlap so need to remove
  if(scale.axis != 1){
        labels <- axTicks(2) / scale.axis
        if(labels[1] == 0) labels[1] <- ''
        axis(2, at = axTicks(2), labels = as.numeric(labels), cex.axis = rel.y.num,
             las = T)
    } else {axis(2, cex.axis = rel.y.num, las = T)}

  #Identify significant trends
  #Whole time series
  mksts <- zyp.trend.vector(x[, Value])
  mkstsp <- round(unlist(mksts[6]), 3)
 
  if(mkstsp < 0.05){
    lmod <- lm(x[, Value] ~ x[, X])
    lmod_c <- unlist(lmod[1])
    lmod_i <- lmod_c[1]
    lmod_s <- lmod_c[2]
    if(lmod_s > 0){
      lines(x[, X], lmod_s * x[, X] + lmod_i,
            col = trend.pos, lty = 1, lwd = 7)
    }
    if(lmod_s < 0){
      lines(x[, X], lmod_s * x[, X] + lmod_i,
            col = trend.neg, lty = 1, lwd = 7)
    }
  }
 
  #Final portion of time series
  mksld <- zyp.trend.vector(x[X > (end.start - 1), Value])
  mksldp <- round(unlist(mksld[6]), 3)
  if(mksldp < 0.05){
    l10_x <- x[X > (end.start - 1), X]
    l10_y <- x[X > (end.start - 1), Value]
    lmod <- lm(l10_y ~ l10_x)
    lmod_c <- unlist(lmod[1])
    lmod_i <- lmod_c[1]
    lmod_s <- lmod_c[2]
    if(lmod_s > 0){
      lines(l10_x, lmod_s * l10_x + lmod_i,
            col = end.trend.pos, lwd=7)
    }
    if(lmod_s < 0){
      lines(l10_x, lmod_s * l10_x + lmod_i,
            col = end.trend.neg, lwd=7)
    }
  }
 
  #Add axis labels
  if(!is.na(stacked)) text(u[1], u[4], labels = stacked, cex = 2, adj = c(-0.5, 1.5))
  if(is.na(stacked)){
    mtext(1, text = x.label, line = x.line, cex = 1.5)
    mtext(2, text = y.label, line = y.line, cex = rel.y.text)
  }
}

#Add axis labels for stacked plots
soe.stacked.axis <- function(x.label, y.label, x.line = 2.6,
                             y.line = 3.5, rel.y.text = 1.5){
  axis(1, cex.axis = 1.5)
  mtext(1, text = x.label, line = x.line, cex = 1.5, outer = T)
  mtext(2, text = y.label, line = y.line, cex = rel.y.text, outer = T)
}

```


```{r Plot options, echo = F}
#Background colors
background   <- '#F4F7F2'
recent       <- '#E6E6E6'
#trend lines
main.pos <- rgb(253/255, 184/255, 99/255,  alpha = 0.8)
main.neg <- rgb(178/255, 171/255, 210/255, alpha = 0.8)
end.pos  <- rgb(230/255, 97/255,  1/255,   alpha = 0.8)
end.neg  <- rgb(94/255,  60/255,  153/255, alpha = 0.8)

```

# Introduction
The Council approved an EAFM Guidance Document in 2016 which outlined a path forward to more fully incorporate ecosystem considerations into marine fisheries management^[http://www.mafmc.org/s/EAFM_Guidance-Doc_2017-02-07.pdf], and revised the document in February 2019^[http://www.mafmc.org/s/EAFM-Doc-Revised-2019-02-08.pdf]. The Council’s stated goal for EAFM is "to manage for ecologically sustainable utilization of living marine resources while maintaining ecosystem productivity, structure, and function." Ecologically sustainable utilization is further defined as "utilization that accommodates the needs of present and future generations, while maintaining the integrity, health, and diversity of the marine ecosystem." Of particular interest to the Council was the development of tools to incorporate the effects of species, fleet, habitat and climate interactions into its management and science programs. To accomplish this, the Council agreed to adopt a structured framework to first prioritize ecosystem interactions, second to specify key questions regarding high priority interactions and third tailor appropriate analyses to address them [@gaichas_framework_2016]. Because there are so many possible ecosystem interactions to consider, a risk assessment was adopted as the first step to identify a subset of high priority interactions [@holsman_ecosystem-based_2017]. The risk elements included in the Council's initial assessment spanned biological, ecological, social and economic issues (Table \ref{riskel}) and risk criteria for the assessment were based on a range of indicators and expert knowledge (Table \ref{allcriteria}).

This document updates the Mid-Atlantic Council’s initial EAFM risk assessment with indicators from the 2019 State of the Ecosystem report. The risk assessment was designed to help the Council decide where to focus limited resources to address ecosystem considerations by first clarifying priorities. Overall, the purpose of the EAFM risk assessment is to provide the Council with a proactive strategic planning tool for the sustainable management of marine resources under its jurisdiction, while taking interactions within the ecosystem into account. 

Many risk rankings are unchanged based on the updated indicators for 2019 and the Council's risk criteria. Below, we highlight only the elements where updated information has changed the perception of risk. In addition, we present new indicators based on Council feedback on the original risk analysis that the Council may wish to include in future updates to the EAFM risk assessment.

\newpage
```{r riskel, echo=FALSE, message=FALSE, warnings=FALSE, results='asis'}

#tab.cap="Risk Elements, Definitions, and Indicators Used\\label{riskel}",

elem <-read.table("riskelements.txt", sep="|", header=F, strip.white = T, stringsAsFactors = F)
elem <- elem[,2:4]
names(elem) <- c("Element", "Definition", "Indicator")

# elem$Element <- factor(all$Element, levels=c("Assessment performance", "F status", "B status", "Food web  (Council Predator)", "Food web  (Council Prey)", "Food web  (Protected Species Prey)",
#                                             "Ecosystem productivity", "Climate", "Distribution shifts", "Estuarine habitat", "Offshore habitat", "Commercial Revenue",
#                                             "Recreational Angler Days/Trips", "Commercial Fishery Resilience (Revenue Diversity)", "Commercial Fishery Resilience (Shoreside Support)",
#                                             "Fleet Resilience", "Social-Cultural", "Commercial", "Recreational", "Control", "Interactions", "Other ocean uses", "Regulatory complexity", 
#                                             "Discards", "Allocation"))

kable(elem, format = "latex", booktabs = T, longtable=T, caption="Risk Elements, Definitions, and Indicators Used\\label{riskel}") %>%
  kable_styling(font_size=8, latex_options=c("repeat_header")) %>%
  column_spec(1, width="2.5cm") %>%
  column_spec(2:3, width="7cm") %>%
  group_rows("Ecological",1,11) %>%
  group_rows("Economic",12,15) %>%
  group_rows("Social",16,17) %>%
  group_rows("Food Production",18,19) %>%
  group_rows("Management",20,25)
  #landscape()


```
 

\newpage
\pagestyle{plain}
```{r allcriteria,  echo=FALSE, message=FALSE, warnings=FALSE, results='asis'}

#tab.cap="Risk Ranking Criteria used for each Risk Element\\label{allcriteria}",

all<-read.table("riskrankingcriteria.txt", sep="|", header=T, strip.white = T, stringsAsFactors = F)
names(all) <- c("Element", "Ranking", "Criteria")
all$Ranking <- factor(all$Ranking, levels=c("Low", "Low-Moderate",  "Moderate-High", "High"))
all$Element <- factor(all$Element, levels=c("Assessment performance", "F status", "B status", "Food web  (MAFMC Predator)", "Food web  (MAFMC Prey)", "Food web  (Protected Species Prey)",
                                            "Ecosystem productivity", "Climate", "Distribution shifts", "Estuarine habitat", "Offshore habitat", "Commercial Revenue",
                                            "Recreational Angler Days/Trips", "Commercial Fishery Resilience (Revenue Diversity)", "Commercial Fishery Resilience (Shoreside Support)",
                                            "Fleet Resilience", "Social-Cultural", "Commercial", "Recreational", "Control", "Interactions", "Other ocean uses", "Regulatory complexity", 
                                            "Discards", "Allocation"))

allwide <- all %>% 
  spread(Ranking, Criteria)

kable(allwide, format = "latex", booktabs = T, longtable=T, caption="Risk Ranking Criteria used for each Risk Element\\label{allcriteria}") %>%
  kable_styling(font_size=8, latex_options=c("repeat_header")) %>%
  column_spec(1, width="2cm") %>%
  column_spec(2:5, width="5cm") %>%
  landscape()


```
\clearpage
\pagestyle{fancy}

# Changes from 2018

## Decreased Risk
Summer flounder fishing mortality (F) status has improved from high risk (F>Fmsy) to low risk (F<Fmsy) based on the new benchmark assessment (Table \ref{sptable}).

Updated commercial fleet diversity (fleet count and fleet diversity) have no long term trends, thus improving from moderate-high risk to low risk according to risk criteria for this element (Table \ref{ecotable}). 

## Increased Risk
No indicators for individual elements have changed enough to warrant increased risk rankings according to the Council risk critiera.

However, we note that most management elements were not re-evaluated for 2019 (Table \ref{spsectable}). Quantitative evaluation of the risks posed by other ocean uses was delayed due to the government shutdown. In addition, the poorer condition of north Atlantic right whales relative to the 2018 report along with the continued increase in ocean temperature indicate that both protected species interactions and climate conditions continue to pose risks to Council-managed fisheries. 

## Re-evaluate Risk
Indicators for recreational opportunities based on updated Marine Recreational Information Program (MRIP) data show generally similar patterns of decreased angler days and trips over the past 10 years, but the declines are less pronounced than measured previously. A reduction from the highest risk ranking to a lower risk category may be warranted. 

# Potential New Indicators
All recreational indicators have been updated with new Marine Recreational Information Program (MRIP) data, and new indicators for recreational diversity are presented in this report at the request of the MAFMC.

## Recreational Diversity
Newly developed indicators for the diversity of recreational effort (i.e. access to recreational opportunities) by mode (party/charter boats, private boats, shore-based), and diversity of catch (NEFMC and MAFMC managed species only) show significant long-term downward trends. The downward effort diversity trend is driven by party/charter contraction (from a high of 24% of angler trips to 6% currently), with a shift towards shorebased angling. Effort in private boats remained stable between 36-37% of angler trips across the entire series. The long-term decrease in catch diversity in the Mid-Atlantic states contrasts with an increase in recreational catch diversity in New England states over the same time period; this trend requires further investigation as SAFMC managed species are not currently tracked separately (Fig. \ref{fig:rec-div})

Recreational diversity indices could be considered as additional risk element(s) to complement the existing Commercial fishery resilience (revenue diversity) element. 

We seek feedback whether the Council would like to include recreational diversity as an indicator for a new risk element, what risk criteria should be applied, and whether the recreational species diversity index should be modified to account for SAFMC managed species. 

```{r rec-div, fig.width = 4, fig.asp = 0.8, fig.cap = paste0("Recreational effort diversity and diversity of recreational catch in the ",region,".")}

recdat <- ecodata::recdat %>% 
  filter(EPU == region_abbr) %>% 
  group_by(Var) %>% 
  mutate(hline = mean(Value))

ylim_re <- c(2e7, 7e7)
ylim_rd <- c(1.75,2.75)
ylim_ra  <- c(1e6, 3.5e6)

series.col <- "black"

rec_div <- recdat %>% 
  filter(Var == "Recreational fleet effort diversity across modes") %>% 
  ggplot() + 
 #Highlight last ten years
  annotate("rect", fill = shade.fill, alpha = shade.alpha,
      xmin = x.shade.min , xmax = x.shade.max,
      ymin = -Inf, ymax = Inf) +
  # annotate("text", 
  #          x = label_loc[label_loc$Var == "Recreational fleet effort diversity across modes",]$xloc,
  #          y = label_loc[label_loc$Var == "Recreational fleet effort diversity across modes",]$yloc,
  #          label = label_loc[label_loc$Var == "Recreational fleet effort diversity across modes",]$labels,
  #          size = letter_size)+
  geom_gls(aes(x = Time, y = Value,
               group = Var),
             alpha = trend.alpha, size = trend.size) +
  geom_line(aes(x = Time, y = Value, color = Var), size = lwd) +
  geom_point(aes(x = Time, y = Value, color = Var), size = pcex) +
  ylim(ylim_rd)+
  scale_x_continuous(expand = c(0.01, 0.01)) +
  scale_color_manual(values = series.col, aesthetics = "color")+
  guides(color = FALSE) +
  ggtitle("Rec. fleet effort diversity")+
  ylab(expression("Effective Shannon")) +
  xlab("")+
  geom_hline(aes(yintercept = hline,
               color = Var),
           size = hline.size,
           alpha = hline.alpha,
           linetype = hline.lty) +
  theme_ts()

rec_div_catch <- recdat %>% 
  filter(Var == "Recreational Diversity of Catch") %>% 
  ggplot() + 
 #Highlight last ten years
  annotate("rect", fill = shade.fill, alpha = shade.alpha,
      xmin = x.shade.min , xmax = x.shade.max,
      ymin = -Inf, ymax = Inf) +
    # annotate("text", 
    #        x = label_loc[label_loc$Var == "Recreational anglers",]$xloc,
    #        y = label_loc[label_loc$Var == "Recreational anglers",]$yloc,
    #        label = label_loc[label_loc$Var == "Recreational anglers",]$labels,
    #        size = letter_size)+
  geom_gls(aes(x = Time, y = Value,
               group = Var),
             alpha = trend.alpha, size = trend.size) +
  geom_line(aes(x = Time, y = Value, color = Var), size = lwd) +
  geom_point(aes(x = Time, y = Value, color = Var), size = pcex) +

  scale_x_continuous(expand = c(0.01, 0.01)) +
  scale_color_manual(values = series.col, aesthetics = "color")+
  guides(color = FALSE) +
  ggtitle("Rec. diversity of catch")+
  ylab(expression("Effective Shannon")) +
  xlab("Time")+
  geom_hline(aes(yintercept = hline,
               color = Var),
           size = hline.size,
           alpha = hline.alpha,
           linetype = hline.lty) +
  theme_ts()


cowplot::plot_grid(#rec_effort, 
                   rec_div, 
                   #rec_anglers, 
                   rec_div_catch,
                   ncol = 1, 
                   align = "hv") +
    theme(plot.margin = unit(c(0.1, 0, 0, 0), "cm"))

```

## Chesapeake Bay Water Quality
Many important MAFMC managed species use estuarine habitats as nurseries or are considered estuarine and nearshore coastal-dependent (summer flounder, scup, black sea bass, and bluefish), and interact with other important estuarine-dependent species (e.g., striped bass and menhaden). An integrated measure of multiple water quality criteria shows a significantly increasing proportion of Chesapeake Bay waters meeting or exceeding EPA water quality standards over time (@zhang_chesapeake_2018; Fig. \ref{fig:cb-attainment}). This pattern was statistically linked to total nitrogen reduction, indicating responsiveness of water quality status to management actions implemented to reduce nutrients. 

This improvement in estuarine water quality could result in a future improvement in the estuarine habitat quality risk ranking for estuarine dependent species. This (currently high risk) ranking could change if other Mid-Atlantic estuaries have similar improvements in water quality and if this overall improvement in water quality moves the EPA assessment of estuarine condition from poor to fair. Estuarine water quality is just one component of estuarine condition. EPA ratings were based on 2003–2006 nearshore and estuarine summer sampling. Coastal waters in the Mid-Atlantic region rated fair to poor for water quality, fair for sediment quality, poor for benthic quality, good to fair for coastal habitat, and fair to poor for fish contamination. 

```{r cb-attainment,fig.width = 5, fig.asp = 0.45, fig.cap = "Estimated water quality standards attainment of Chesapeake Bay tidal waters for the combined assessment of dissolved oxygen, underwater bay grasses/water clarity and chlorophyll a using rolling three year assessment periods."}

minlab <- seq(1985,2015,5)
maxlab <- seq(1987,2017,5)



ches_bay_wq %>% 
  mutate(hline = mean(Value)) %>% 
  ggplot(aes(x = Time, y = Value)) +
       annotate("rect", fill = shade.fill, alpha = shade.alpha,
      xmin = x.shade.min , xmax = x.shade.max,
      ymin = -Inf, ymax = Inf) +
  geom_line() +
  geom_point() +
  geom_gls() +
      ylab("Estimated attainment (%)") +
  ggtitle("Chesapeake Bay Estimated Water Quality Standards Attainment") +
  scale_x_continuous(breaks = minlab,labels = paste0(minlab,"-",maxlab),expand = c(0.01, 0.01)) +
    geom_hline(aes(yintercept = hline),
           size = hline.size,
           alpha = hline.alpha,
           linetype = hline.lty) +
  theme_ts() +
  theme(axis.text.x = element_text(angle = 20, vjust = 0.65),
          plot.title = element_text(size = 10))
```

\newpage

```{r sptable,  echo=FALSE, message=FALSE, warnings=FALSE, results='asis'}

#tab.cap="Species level risk analysis results; l=low risk (green), lm= low-moderate risk (yellow), mh=moderate to high risk (orange), h=high risk (red)\\label{sptable}",

# spplist     oc,  sc,  flk, scp, bsb, mack, but, lsq, ssq, gtile,  btile,  blu, dog, monk
risk.species<-data.frame(
  Species = c("Ocean Quahog", "Surfclam", "Summer flounder", "Scup", "Black sea bass", "Atl. mackerel", "Butterfish", "Longfin squid", "Shortfin squid", "Golden tilefish", "Blueline tilefish", "Bluefish", "Spiny dogfish", "Monkfish", "Unmanaged forage", "Deepsea corals"),
  Assess  = c("l", "l", "l",  "l", "l", "l",  "l", "lm", "lm", "l", "h", "l", "lm", "h", "na", "na"),
  Fstatus = c("l", "l", "l", "l", "l", "h", "l", "lm", "lm", "l", "h", "l", "l", "lm", "na", "na"),
  Bstatus = c("l", "l", "lm", "l", "l", "h", "l", "lm", "lm", "lm", "mh", "lm", "lm", "lm", "na", "na"),
  FW1Pred = c("l", "l", "l", "l", "l", "l", "l", "l", "l", "l", "l", "l", "l", "l", "l", "l"),
  FW1Prey = c("l", "l", "l", "l", "l", "l", "l", "l", "l", "l", "l", "l", "l", "l", "lm", "l"),
  FW2Prey = c("l", "l", "l", "l", "l", "l", "l", "lm", "lm", "l", "l", "l", "l", "l", "lm", "l"),
  Climate = c("h", "mh", "lm", "lm", "mh", "lm", "l", "l", "l", "mh", "mh","l", "l", "l", "na", "na"),
  DistShift = c("mh", "mh", "mh", "mh", "mh", "mh", "h", "mh", "h", "l", "l", "mh", "h", "mh", "na", "na"),
  EstHabitat = c("l", "l", "h", "h", "h", "l", "l", "l", "l", "l", "l", "h", "l", "l", "na", "na")#,
#  OffHabitat = c("na", "na", "l", "l", "l", "l", "l", "l", "h", "na", "na", "na", "l", "l", "na", "na")#,
)

# these elements were removed by the council
#  PopDiv = c("na", "na", "na", "na", "na", "na", "na", "na", "na", "na", "na", "na", "na", "na"),
#  FoodSafe = c(),

# one column test
# risk.species %>%
#   mutate(Fstatus = 
#     cell_spec(Fstatus, format="latex", color = "black", align = "c", background =factor(Fstatus, c("na", "l", "lm", "mh", "h"),c("white", "green", "yellow", "orange", "red")))) %>%
#   kable(risk.species, format="latex", escape = F, booktabs = T, linesep = "")

#generalize to all
risk.species %>%
  mutate_at(vars(-Species), function(x){ 
    cell_spec(x, format="latex", color = "gray", align = "c", background =factor(x, c("na", "l", "lm", "mh", "h"),c("white", "green", "yellow", "orange", "red")))}) %>%
  kable(risk.species, format="latex", escape = F, booktabs = T, linesep = "", 
        caption="Species level risk analysis results; l=low risk (green), lm= low-moderate risk (yellow), mh=moderate to high risk (orange), h=high risk (red)\\label{sptable}") %>%
  kable_styling(latex_options = "scale_down") #%>%
  #kable_as_image()


```

```{r ecotable,  echo=FALSE, message=FALSE, warnings=FALSE, results='asis'}
  
#tab.cap="Ecosystem level risk analysis results; l=low risk (green), lm= low-moderate risk (yellow), mh=moderate to high risk (orange), h=high risk (red)\\label{sptable}",

risk.eco<-data.frame(
  System = c("Mid-Atlantic"),
  EcoProd = c("lm"),
  #EcoDiv = c("lm"),
  CommRev = c("mh"),
  RecVal = c("h"),
  FishRes1 = c("l"),
  FishRes4 = c("mh"),
  #CommJobs = c("mh"),
  #RecJobs = c("l"),
  FleetDiv = c("l"),
  Social = c("lm"),
  ComFood = c("h"),
  RecFood = c("mh")
)

#make table
risk.eco %>%
  mutate_at(vars(-System), function(x){ 
    cell_spec(x, format="latex", color = "gray", align = "c", background =factor(x, c("na", "l", "lm", "mh", "h"),c("white", "green", "yellow", "orange", "red")))}) %>%
  kable(risk.eco, format="latex", escape = F, booktabs = T, linesep = "",
        caption="Ecosystem level risk analysis results; l=low risk (green), lm= low-moderate risk (yellow), mh=moderate to high risk (orange), h=high risk (red)\\label{ecotable}") %>%
  kable_styling(latex_options = "scale_down") #%>%
  #kable_as_image()

```

```{r spsectable,  echo=FALSE, message=FALSE, warnings=FALSE, results='asis'}
#tab.cap="Species and sector level risk analysis results; l=low risk (green), lm= low-moderate risk (yellow), mh=moderate to high risk (orange), h=high risk (red)\\label{sptable}",

risk.sppsector<-data.frame(
  Species = c("Ocean Quahog-C", "Surfclam-C", "Summer flounder-R", "Summer flounder-C","Scup-R", "Scup-C","Black sea bass-R", "Black sea bass-C","Atl. mackerel-R", "Atl. mackerel-C","Butterfish-C", "Longfin squid-C", "Shortfin squid-C", "Golden tilefish-R", "Golden tilefish-C","Blueline tilefish-R","Blueline tilefish-C", "Bluefish-R", "Bluefish-C","Spiny dogfish-R", "Spiny dogfish-C", "Unmanaged forage", "Deepsea corals"),
  MgtControl = c(1,1,3,2,1,1,4,2,1,1,1,1,1,9,1,1,1,2,1,1,1,9,9),
  TecInteract = c(1,1,1,3,1,3,1,2,1,2,2,3,2,1,1,1,1,1,1,1,3,9,9),
  OceanUse = c(2,2,2,2,2,2,3,4,1,3,3,4,2,1,1,1,1,1,2,1,3,9,3),
  RegComplex = c(1,1,4,3,3,3,4,3,1,4,4,4,2,1,1,3,3,1,2,1,3,9,9),
  Discards = c(1,1,4,2,3,3,3,2,1,2,3,4,1,1,1,1,1,3,2,1,2,9,9),
  Allocation = c(1,1,4,4,1,1,4,4,4,4,1,4,1,1,1,4,4,4,4,1,4,9,9) 
)

#convert to text for consistency
risk.sppsector <- risk.sppsector %>%
     mutate_at(vars(-Species), function(x){
       recode(x,'1'="l",'2'="lm",'3'="mh",'4'="h",'9'="na")}) %>%
     as.data.frame()

#make table
risk.sppsector %>%
  mutate_at(vars(-Species), function(x){ 
    cell_spec(x, format="latex", color = "gray", align = "c", background =factor(x, c("na", "l", "lm", "mh", "h"),c("white", "green", "yellow", "orange", "red")))}) %>%
  kable(risk.sppsector, format="latex", escape = F, booktabs = T, linesep = "",
        caption="Species and sector level risk analysis results; l=low risk (green), lm= low-moderate risk (yellow), mh=moderate to high risk (orange), h=high risk (red)\\label{spsectable}") %>%
  kable_styling(font_size =  9) #%>%
  #kable_as_image()

```


# References