---
title: "Mission: Iconic Reefs and National Coral Reef Monitoring Program Benthic Analyses"
subtitle: "Summary Report for 2022 and 2024"
date: "`r format(Sys.time(), '%d %B %Y')`"
bibliography: References.bib
reference-location: section
reference-section-title: Bibliography
nocite: "@NOAA2022a; @NOAA2022b; @NOAA2024a; @NOAA2024b"
csl: special.csl
link-citations: true
format:
html:
theme: lumen
toc: true
number-sections: true
code-fold: true
code-tools: true
fig-cap-location: bottom
fig-align: center
self-contained: true
title-block-banner: true
title-block-style: default
include-before-body: logos.html
include-in-header:
text: |
<style>
.csl-entry {
margin-bottom: .7em;
}
</style>
pdf:
fig-cap-location: bottom
author: "Nicole Krampitz, Dione Swanson, Julia Mateski, Lexie Sturm, and Shay Viehman"
mainfont: Arial
toc: true
toc-title: |
\clearpage
Table of Contents
number-sections: true
fig-numbering: false
table-caption: false
include-in-header:
text: |
\usepackage{graphicx}
\usepackage{titling}
\usepackage{float}
\pretitle{%
\begin{center}
\includegraphics[width=6cm]{mir_logo.png}\\[1em]
\includegraphics[width=4cm]{noaa_crcp.jpg}\\[1em]
\LARGE
}
\usepackage{caption}
\captionsetup[figure]{labelformat=empty}
\captionsetup[table]{labelformat=empty}
execute:
freeze: auto
echo: true
warning: false
message: false
---
\newpage
## Introduction
[Mission: Iconic Reefs](https://www.fisheries.noaa.gov/southeast/habitat-conservation/restoring-seven-iconic-reefs-mission-recover-coral-reefs-florida-keys) (M:IR) is a multi-institutional initiative designed to restore ecological function and biodiversity across key reef sites in the Florida Keys. Led by NOAA in collaboration with federal, state, institutional, and non-profit partners, M:IR targets nearly three million square feet of reef at seven key coral reefs in the Florida Keys National Marine Sanctuary (FKNMS): Carysfort (North and South), Cheeca Rocks, Eastern Dry Rocks, Horseshoe Reef, Looe Key, Newfound Harbor, and Sombrero Reef (Figure 1). Across these reefs, M:IR aims to restore coral cover to an average of 25% across the seven sites -- an inititaitve unparalleled in scope and scale that will require hundreds of thousands of coral restoration outplants and over a decade of phased intervention.
{fig-align="center"}
These Iconic reefs were selected for restoration due to their diversity in reef types, levels of anthropogenic disturbances, initial coral cover, geographic distribution, and cultural significance. Restoration progress at these reefs is informed by ongoing monitoring of coral reef benthic communities and coral populations at M:IR sites using a range of metrics, such as coral cover, size structure, mortality prevalence, and diversity. Such multi-metric monitoring is essential to evaluating whether current restoration practices are effective in establishing and maintaining a healthy coral reef, as restoration ‘success’ encompasses a range of ecological and management metrics rather than a single measurement (e.g., coral cover) [@Goergen2020].
In addition to monitoring within the M:IR reefs, it is also important to assess the surrounding reefs in the Florida Keys to interpret restoration progress. The National Coral Reef Monitoring Program (NCRMP), which conducts standardized reef assessments across the Florida Keys on a biennial basis, provides a solid framework for comparison. To align with this program, the M:IR team applied NCRMP survey methodology within its Iconic Sites in 2022 and 2024 to complement NCRMP surveys. These surveys focused on shallow (\<30 m), hard-bottom reef habitats, using a stratified-random, one-stage design within 50 m × 50 m grid cells to ensure representative sampling across depth and rugosity strata [@Ault2021]. The survey design ensures that survey sites are allocated by cross-shelf zone, rugosity type, and depth strata, with sites distributed around the Florida Keys from nearshore to offshore to a maximum depth of 30 m. The NCRMP sample grid was spatially joined with the seven distinct M:IR zones (Figure 1). Survey sites were proportionally allocated within NCRMP strata to maintain consistent representation across habitats.
```{r libraries, include = FALSE}
# LIBRARIES
library(ggplot2)
library(dplyr)
library(tidyr)
library(cowplot)
library(grid)
library(patchwork)
library(ggthemes)
library(extrafont)
library(grid)
library(ggthemes)
library(ggpattern)
library(extrafont)
library(plotly)
library(tidyverse)
library(ncrmp.benthics.analysis)
library(ggplot2)
library(cowplot)
library(gridExtra)
library(egg)
library(mgcv)
# source("R Scripts/total_dis_code.R")
source("../Project/R Scripts/MIR_make_size_bins.R")
source("../Project/R Scripts/theme_publication.R")
source("../Project/R Scripts/QuadPanel.R")
# source("R Scripts/bleaching_estimates.R")
# source("R Scripts/total_bleach_code.R")
source("../Project/R Scripts/dis_ble_prev_ratio_est.R")
##code for terminal rendering:
#cd ~/Downloads/MIR-Krampitz/docs
#quarto render index.qmd
```
```{r MIR species, echo = FALSE}
MIR_species <- c( "ACR CERV", "ACR PALM", "COL NATA", "DEN CYLI",
"DIC STOK", "DIP LABY", "EUS_ AST", "MEA MEAN",
"PSE CLIV", "PSE STRI", "MON CAVE" )
```
Corals and benthic communities were monitored using two different NCRMP surveys: the Benthic Community Assessment survey and the Coral Demographics survey (NOAA CRCP 2022a, 2022b, 2024a, 2024b). The Benthic Community Assessment survey includes: 1) benthic cover (%) estimates along a 15-m transect with a line point intercept method, 2) presence/absence of Endangered Species Act (ESA)-listed coral species [@NMFS2014], 3) abundance of key macroinvertebrates, and 4) reef rugosity measurements within a 15 m x 2 m belt-transect area (NOAA CRCP 2022a, 2024b). At the same site, coral demographics were surveyed within a 10 m x 1 m belt-transect area (NOAA CRCP 2022b, 2024b). NCRMP coral demographic survey data were combined with complementarily allocated survey data from Florida Fish and Wildlife Conservation Commission's Disturbance Response Monitoring ([DRM](https://coraldrm.org/)). In all coral demographics surveys, all live coral colonies ≧ 4 cm were counted, identified to species, measured to the nearest centimeter, and estimates were made of the proportion per colony of any present mortality (recent or old), disease (absent, present: slow, fast), and/or bleaching (none, total, partial, paling). Only live coral colonies were included in the survey; dead colonies with 100% mortality were not surveyed (e.g., colonies killed by coral disease). Juvenile corals (\< 4 cm) were reported for species richness only and were not included in counts, size measurements, or estimates of condition (NOAA CRCP 2022b); in 2024, juveniles of select coral species were counted (NOAA CRCP 2024b). In 2024, NCRMP benthic surveys also included large area imaging (LAI) collected at most survey sites; analyses of these LAI data are in progress and not presented here.
To allow for direct comparisons of benthic communities inside the M:IR restoration areas with control, non-restored areas across the Florida Keys reefs, the NCRMP dataset was restricted to strata types and depth zones (0-12 m) found within M:IR sites. Statistical comparisons were conducted using a two-tailed t-test between M:IR and NCRMP estimates for each survey year (2022 or 2024) and between years for M:IR and NCRMP. This report focuses on coral populations and benthic community metrics; a separate [report](https://rob-harper.github.io/MissionIconicReef) can be found for fish population assessments [@Blondeau2025].
NCRMP analyses scripts for corals and benthic communities are open source and available at [NCRMP Benthic R package](https://github.com/MSE-NCCOS-NOAA/NCRMP_benthics) [@Groves2025].
## Benthic Survey Effort
In 2022, a total of 584 surveys occurred (Table 1) across 403 sites [@Viehman2023]. In 2024, 771 surveys occurred across 542 sites [@GroveInPrep].
```{r Benthic Survey Summaries, results='asis', echo=FALSE}
###### HTML VERSION ############
if (knitr::is_html_output()) {
TableCaption <- "<div style='font-size:0.9em; color: #666666; margin-bottom:0.5em; text-align: center;'>
Table 1. Number of demographic (Demo) and benthic community assessment (BCA) surveys completed by M:IR, NCRMP, and DRM in 2022 and 2024. </div>"
cat(TableCaption)
#pulling the whole table from a different script because it's messy due to the format
source("../Project/R Scripts/html_summary_table.R")
} else if (knitr::is_latex_output()) {
#for pdf, using a photo of the table from excel instead because format is very messy
cat("\\begin{center}\n")
#changing the table caption to be more similar to html, with light coloring and cnetering
cat("\\textcolor{darkgray}{\\textbf{Table 1.} Number of demographic (Demo) and benthic community assessment (BCA) surveys completed by M:IR, NCRMP, and DRM in 2022 and 2024.}\n\n")
#this is the photo of the table
cat("\\includegraphics[width=\\linewidth]{../Project/images/Proper_Table.png}\n")
#this centers the photo
cat("\\end{center}\n")
}
```
### Coral Demography
For the purpose of this comparison, surveys were divided into two groups: inside M:IR and outside. Almost all of the inside (M:IR) surveys were completed by the M:IR team; however, if an NCRMP or DRM sample happened to fall within the boundaries of an M:IR iconic reef, it was added to the inside group. The inverse occurred if an M:IR survey happened to fall outside the M:IR boundary.
```{r prepare data for benthic cover table, include = FALSE}
##Let's just start with MIR/NCRMP only (benthic cover)
NCRMP_MIR <- read_csv("../Project/data/NCRMP_MIR_cora_2022_24_psu_grps.csv") %>%
group_by(YEAR, strat_cora, analysis_group) %>%
summarise(count = n()) %>%
filter(analysis_group != "NCRMP_NA") %>%
rename(STRAT = strat_cora,
PROT = analysis_group) %>%
mutate(description = case_when(
STRAT == "CFK01" ~ "Inshore reef, all relief types and depths",
STRAT == "CFK02" ~ "Mid-channel patch reef, all relief types and depths",
STRAT == "CFK03" ~ "Offshore patch reef, all relief types and depths",
STRAT == "CFK04" ~ "Forereef, low relief, shallow (0-6 m)",
STRAT == "CFK05" ~ "Forereef, high relief, shallow (0-6 m)",
STRAT == "CFK06" ~ "Forereef, low relief, mid-shallow (6-12 m)",
STRAT == "CFK07" ~ "Forereef, high relief, mid-shallow (6-12 m)"))
table <- NCRMP_MIR %>%
pivot_wider(id_cols = STRAT:description, names_from = YEAR, values_from = count) %>%
mutate(`Study area` = case_when(
PROT == "MIR_GRP"~ "Inside (M:IR)",
PROT == "NCRMP_GRP" ~ "Outside (NCRMP)"
)) %>%
select(`Study area`, STRAT, description, `2022`, `2024`) %>%
arrange(`Study area`, STRAT)
name <- c("<b>Study Area</b>", "<b>Strata Name</b>", "<b>Strat Description</b>", "<b>2022 Sample Number</b>", "<b>2024 Sample Number</b>")
```
```{r create coral demo table data, include = FALSE}
###Demographics (including DRM)
DRM <- read_csv("../Project/data/DRM_cora_2022_24_psu_grps.csv") %>%
group_by(YEAR, strat_cora, analysis_group) %>%
summarise(count = n()) %>%
filter(analysis_group != "NCRMP_NA") %>%
rename(STRAT = strat_cora,
PROT = analysis_group) %>%
mutate(description = case_when(
STRAT == "CFK01" ~ "Inshore reef, all relief types and depths",
STRAT == "CFK02" ~ "Mid-channel patch reef, all relief types and depths",
STRAT == "CFK03" ~ "Offshore patch reef, all relief types and depths",
STRAT == "CFK04" ~ "Forereef, low relief, shallow (0-6 m)",
STRAT == "CFK05" ~ "Forereef, high relief, shallow (0-6 m)",
STRAT == "CFK06" ~ "Forereef, low relief, mid-shallow (6-12 m)",
STRAT == "CFK07" ~ "Forereef, high relief, mid-shallow (6-12 m)")) %>%
rbind(., NCRMP_MIR) %>%
group_by(STRAT, PROT, YEAR, description) %>%
summarise(count = sum(count)) %>%
ungroup()
table_2 <- DRM %>%
pivot_wider(id_cols = STRAT:description, names_from = YEAR, values_from = count) %>%
mutate(`Study area` = case_when(
PROT == "MIR_GRP"~ "Inside (M:IR)",
PROT == "NCRMP_GRP" ~ "Outside <br> (NCRMP + DRM)"
)) %>%
select(`Study area`, STRAT, description, `2022`, `2024`) %>%
arrange(`Study area`, STRAT)
```
```{r print coral demo table, results = 'asis', echo = FALSE}
table_2 <- table_2 %>% ungroup() %>%
mutate(across(everything(), ~ ifelse(is.na(.), "--", .)))
####### HTML VERSION ##############
if (knitr::is_html_output()) {
sjPlot::tab_df(
table_2,
title = "<div style='text-align:center'> Table 2. The number of coral demographic sites sampled inside and outside M:IR areas in each stratum (includes DRM data).</div>",
escape = FALSE,
show.footnote = TRUE,
alternate.rows = TRUE,
col.header = name,
CSS = list(
css.firsttablecol = "width: 140px;", # widen first column
css.caption = "font-weight: normal;",
css.tdata = c(rep("", ncol(table) - 2), "text-align: center;", "text-align: center;"))
)
######### PDF VERSION ###########
} else if (knitr::is_latex_output()) {
table_2 <- table_2 %>% mutate(`Study area` = case_when(
`Study area` == "Outside <br> (NCRMP + DRM)" ~ "Outside (NCRMP + DRM)",
TRUE ~ `Study area`##PDF doesn't accept break
))
# Styled caption printed manually before the table
cat("\\textcolor{darkgray}{\\textbf{Table 2.} The number of coral demographic sites sampled inside and outside M:IR areas in each stratum (includes DRM data).}\n")
# Print the table with kableExtra
kableExtra::kbl(
table_2,
booktabs = TRUE,
caption = NULL,
col.names = c("Study Area", "Strata", "Strata Description", "2022", "2024")
) %>%
kableExtra::kable_styling(
latex_options = c("striped", "hold_position"),
font_size = 10,
position = "center"
) %>%
kableExtra::row_spec(0, bold = TRUE, align = "c") %>%
print()
}
```
### Benthic Composition
Disturbance Response Monitoring (DRM) only performs coral demographic surveys (Table 1) and therefore DRM surveys are not included in the summary of Benthic Community Assesement surveys (Table 3).
```{r create benthic cover table, echo = FALSE, results = 'asis'}
table <- table %>% ungroup() %>%
mutate(across(everything(), ~ ifelse(is.na(.), "--", .)))
if (knitr::is_html_output()) {
####### HTML VERSION ##############
sjPlot::tab_df(
table,
title = "<div style='text-align:center'>Table 3. Number of Benthic Commnity Assesement (BCA) sites sampled inside and outside M:IR areas in each stratum.</div>",
escape = FALSE,
show.footnote = TRUE,
alternate.rows = TRUE,
col.header = name,
CSS = list(
css.caption = "font-weight: normal;",
css.tdata = c(rep("", ncol(table) - 2), "text-align: center;", "text-align: center;")
)
)
######### PDF VERSION ###########
} else if (knitr::is_latex_output()) {
# Manually print a styled caption with label BEFORE the table
cat("\\textcolor{darkgray}{\\textbf{Table 3.} Number of Benthic Community Assessment (BCA) sites sampled inside and outside M:IR areas in each stratum.}\n")
# Print the table itself in kableExtra
kableExtra::kbl(
table,
booktabs = TRUE,
caption = NULL, # no caption here
align = c("l", "l", "l", "c", "c"),
col.names = c("Study Area", "Strat", "Strat Description", "2022", "2024")
) %>%
kableExtra::kable_styling(
latex_options = c("striped"),
font_size = 10,
position = "center"
) %>%
kableExtra::row_spec(0, bold = TRUE, align = "c") %>%
print()
}
```
\newpage
## Benthic Community Composition
Temporal trends in benthic cover of functional groups such as coral and macroalgae cover provide information on changes to the benthic community composition over time.
```{r fig 1: benthic cover data, fig.width=12, fig.asp=0.6, echo = FALSE}
#get cover categories
cat_codes <- read_csv("../Project/data/bcov_cat_codes.csv") %>%
dplyr::rename(cover_group = CODE)
#load and clean cover data (MIR and NCRMP, no DRM)
cvr <- read_csv("../Project/data/fk2022_24_NCRMP_MIR_bcov_catcd_dom_NK.csv")%>%
dplyr::rename(
cover_group = COVER_CAT_CD,
REGION = analysis_group,
avCvr = wavp,
SE = se_wp,
n_sites = n
) %>%
left_join(cat_codes)
```
```{r fig_benthic_cover_data (OG), echo = FALSE}
#wanted_cover_cats <- c("HCORA" , "MACAL")
#Formatting Data for Combined Macroalgal/Hard Coral plot
cvr_analysis <- cvr %>%
group_by(REGION, COV_CAT_A, YEAR) %>%
# filter(COV_CAT_A %in% wanted_cover_cats) %>%
summarise(
avCvr = sum(avCvr),
SE = sqrt(sum(SE^2))#,
# LCI_1.96 = avCvr - 1.96 * SE,
# UCI_1.96 = avCvr + 1.96 * SE,
# LCI = avCvr - SE,
# UCI = avCvr + SE
)
#change names --> more readable
new_dat <- cvr_analysis %>%
mutate(COV_CAT_A = recode(COV_CAT_A,
"HCORA" = "Hard Coral",
"MACAL" = "Macroalgae",
"RAMIC" = "Ramicrusta",
"TURFA" = "Turf Algae",
"SPONG" = "Sponge",
"SCORA" = "Octocoral",
"CCALG" = "Crustose Coralline Algae (CCA)",
"OTHER" = "Other"
))
##Other includes Cyanobacteria, Bare, Palythoa, Other organisms, Millepora, Seagrass
## read in significance
sig <- read_csv("../Project/data/significance.csv")
dodge_width <- 0.06
# Get unique categories
categories <- unique(new_dat$COV_CAT_A)
# Define the plot for multi-function use in tabset
make_plot <- function(df, category) {
df_filtered <- df %>% filter(COV_CAT_A == category)
plot <- df_filtered %>%
ggplot(aes(x = as.integer(YEAR), y = avCvr *100, group = interaction(COV_CAT_A, REGION), color = REGION)) +
geom_point(size = 2, position = position_dodge(width = 0.06)) +
geom_line(size = 1, position = position_dodge(width = 0.06)) +
geom_errorbar(aes(ymin = pmax(0, avCvr * 100 - SE * 100), ymax = avCvr * 100 + SE * 100),
width = 0.1,
position = position_dodge(width = 0.06)) +
scale_color_manual(
values = c("MIR_GRP" = "#0072B2", "NCRMP_GRP" = "#D55E00"),
labels = c("M:IR", "NCRMP"),
name = ""
) +
labs(x = "Year", y = "Cover (%)", title = category ) +
# facet_wrap(~ COV_CAT_A) + ##This was when we filtered for just hard corals/macroalgae. If we output as a pdf, might be necessary to bring back
scale_x_continuous(breaks = c(2022, 2024)) +
coord_cartesian(ylim = c(-1, NA)) +
#scale_y_continuous(, labels = scales::percent_format(accuracy = 2) ) +
theme_Publication(base_size = 20) +
theme(
strip.text = element_text(face = "bold", size = 16),
legend.position = "bottom",
legend.direction = "horizontal",
legend.key.size = unit(1.2, "cm"),
legend.margin = margin(0, 0, 0, 0, "cm"),
legend.text = element_text(size = 16),
legend.background = element_blank(),
legend.title = element_blank()
) + guides(color = guide_legend(override.aes = list(size = 4)))
##Conditional Formatting
if (category == "Ramicrusta"){
plot <- plot + coord_cartesian(ylim = c(-0, 3))
}
##Adding the significance labels
#This may take a fair amount of tweaking ##
sig_cover <- sig %>%
filter(metric == "cover") %>%
filter(variable == category)
# For between years find a midpoint to go on the plot
midpoints <- df_filtered %>%
group_by(REGION) %>%
summarise(
y_mid = 1.02 * mean(avCvr * 100, na.rm = TRUE))
# For between groups find the y min to go on the plot
ymin <- df_filtered %>%
group_by(COV_CAT_A) %>%
summarise(y_min = 85 * (min(avCvr)- max(SE))) %>% # this hopefully gets a nice plot minimum just below wherever it ends
dplyr::rename(variable = COV_CAT_A) #for joining purposes
# Add them to the significance table
sig_cover <- sig_cover %>%
left_join(midpoints, by = "REGION") %>%
left_join(ymin, by = "variable")
##Plot it with annotations
(plot +
### Significance BETWEEN GROUPS on X axis
geom_text(data = sig_cover %>% filter (is.na(REGION)),
aes(x = YEAR, y = -0.5), label = "+",
inherit.aes = FALSE, size = 12, color = "black") +
### Significance BETWEEN YEARS on graph
geom_text(data = sig_cover %>% filter (is.na(YEAR)),
aes(x = 2023, y = y_mid, color = REGION), label = "*",
inherit.aes = FALSE, size = 15))
}
```
```{r COVER HTML LOOP, fig.cap="Figure 2. Cover estimates for main benthic groups inside and outside M:IR areas for 2022 and 2024. A black (+) indicates a significant difference between groups in that survey year (i.e., 2022 or 2024). A colored (*) indicates significance difference between years for the given group (i.e., M:IR or NCRMP). The category ‘Other’ comprises Cyanobacteria, Palythoa, <i>Millepora</i>, Seagrass, bare substrate, and other miscellaneous organisms.", fig.width=9, results='asis', fig.asp=0.8, echo=FALSE}
# Define manual tab order
ordered_categories <- c("Hard Coral", "Macroalgae", "Octocoral", "Sponge", "Turf Algae", "Crustose Coralline Algae (CCA)", "Ramicrusta", "Other")
if (knitr::is_html_output()){
# Loop through categories and print plots
cat("::: {.panel-tabset}\n\n") # open tabset container
walk(ordered_categories, function(cat) {
cat("### ", cat, " {.tab-pane}\n\n", sep = "")
cat_plot <- make_plot(new_dat, cat)
print(cat_plot)
cat("\n\n")
})
cat(":::\n") # close tabset container
}
```
```{r COVER PDF, fig.cap = "Figures 2a-f. Cover estimates for main benthic groups inside and outside M:IR areas for 2022 and 2024. A black plus (+) indicates a significant difference between groups in that survey year (i.e., 2022 or 2024). A colored asterisk (*) indicates significance difference between years for the given group (i.e., M:IR or NCRMP).", fig.width= 15, fig.height = 10, fig.asp = 1.2, echo = FALSE, fig.pos='H', results = 'asis'}
if (knitr::is_latex_output()) {
##Put two graphs per page
for (i in seq(1, length(ordered_categories), by = 2)) {
plots_to_print <- purrr::map(ordered_categories[i:min(i+1, length(ordered_categories))],
~ make_plot(new_dat, .x))
combined <- patchwork::wrap_plots(plots_to_print, ncol = 1)
# Print with increased height
print(combined + patchwork::plot_layout(heights = rep(2, length(plots_to_print))))
cat("\n\n")
}
}
```
## Coral Demography
Although several metrics reported here encompass all observed coral species, ten focal species are highlighted as particularly important to M:IR's goal: *Acropora palmata*, *Acropora cervicornis*, *Colpophyllia natans*, *Dendrogyra cylindrus*, *Dichocoenia stokesii*, *Diploria labyrinthiformis*, *Meandrina meandrites*, *Montastraea cavernosa*, *Pseudodiploria clivosa*, and *Pseudodiploria strigosa*. These species have experienced drastic population declines and are targeted for restoration at the M:IR sites. Specifically, most were severely affected by White Band Disease (*Acropora* spp.) or Stony Coral Tissue Loss Disease (SCTLD).
The following table summarizes the number of observations across all strata and survey years of each focal species and identifies the most abundant species encountered.
```{r make data for top species, include = FALSE}
##Everything calculated for density is done so by region/year, so maybe just do total colonies observed?
df <- read.csv( "../Project/data/FK2024_NCRMP_DRM_MIR_corsz_AR_NK.csv") %>%
filter(N > 0,
analysis_group != "NCRMP_NA") %>%
group_by(SPECIES_CD, YEAR) %>%
summarise(observed = n()) %>%
arrange(desc(observed)) %>%
pivot_wider(id_cols = SPECIES_CD, names_from = YEAR, values_from = observed)
##Top 15 species
top <- df %>% ungroup() %>% dplyr::slice(1:15)
##MIR Species
target_species <- c("ACR PALM", "ACR CERV", "COL NATA", "DEN CYLI",
"DIC STOK", "DIP LABY", "EUS_ AST", "MEA MEAN",
"PSE CLIV", "PSE STRI", "MON CAVE")
MIR <- df %>% filter(
SPECIES_CD %in% target_species
)
##Together
table <- full_join(MIR, top)
##Add names from origial DRM data
names <- read_csv("../Project/data/DRM_FLK_2014_2024_2stage_coral_demographics.csv") %>% select(SPECIES_CD, SPECIES_NAME) %>% unique()
final <- left_join(table, names) %>%
select(SPECIES_NAME, SPECIES_CD, `2022`, `2024`) %>%
arrange(desc(`2022`)) %>%
mutate(SPECIES_NAME =
ifelse(
`SPECIES_CD` %in% target_species,
paste0("*", "<i>", SPECIES_NAME, "</i>"),
paste0("<i>", SPECIES_NAME, "</i>")
)) %>%
mutate(`2024` = replace_na(`2024`, 0)) %>%
mutate(across(c(`2022`, `2024`), ~ format(., big.mark = ",", scientific = FALSE))) %>%
ungroup() %>%
select(-SPECIES_CD)
```
```{r create top species table, echo = FALSE, results = 'asis'}
######### HTML OUTPUT ###############
if (knitr::is_html_output()) {
sjPlot::tab_df(
final,
title = "<div style='text-align:center'>Table 4. Number of colonies observed for focal M:IR species and the most commonly encountered non-focal species in 2022 and 2024. </div>",
alternate.rows = TRUE,
show.footnote = TRUE,
footnote = "* indicates M:IR focal species",
col.header = c(
"<b>Species</b>",
"<b>Colonies in 2022</b>",
"<b>Colonies in 2024</b>"
),
CSS = list(
css.caption = "font-weight: normal;",
css.table = "margin-left: auto; margin-right: auto; width: 90%;",
css.tdata = c(rep("", ncol(table) - 2), "text-align: center;", "text-align: center;")
))
######### PDF OUTPUT ###############
} else if (knitr::is_latex_output()) {
##remove the italics styling for HTML
final <- final %>% mutate(SPECIES_NAME = gsub("</?i>", "", SPECIES_NAME))
# Custom caption printed manually above the table
cat("\\textcolor{darkgray}{\\textbf{Table 4.} Number of colonies observed for focal M:IR species and the most commonly encountered non-focal species in 2022 and 2024.}\n\n")
# Render table with kableExtra
kableExtra::kbl(
final,
booktabs = TRUE,
col.names = c("Species", "Colonies in 2022", "Colonies in 2024"),
align = c("l", "c", "c")
) %>%
kableExtra::kable_styling(
latex_options = c("striped", "hold_position"),
font_size = 10,
position = "center"
) %>%
kableExtra::row_spec(0, bold = TRUE, align = "c") %>%
kableExtra::column_spec(1, italic = TRUE) %>%
print()
}
```
### Focal Species
The stratified random surveys are designed to optimize sampling efficiency and provide high precision estimates of coral population metrics at relatively low sample sizes [@Smith2011]. The effectiveness of restoration can be quantified by comparison of these metrics between M:IR sites and similar habitats outside these restoration areas.
The following figures for these metrics are grouped by each focal species. Significant statistical differences using a two-tailed t-test are indicated for differences between survey years (2022 vs. 2024) and between areas (M:IR vs. NCRMP) for density, old mortality, and mean colony size.
```{r prep density data, fig.width=12, fig.asp=0.6, echo = FALSE}
# Read in clean dens data
den <- read_csv("../Project/data/fk2022_24_NCRMP_MIR_corabd_dom_NK.csv") %>%
rename(
avDen = wavdns,
SE = se_wdns,
CV_perc = wcv_dns,
n_sites = n)
#load full species name
names <- read_csv("../Project/data/NCRMP_FKEYS2022_Benthic_Data02_CoralDemographics.csv") %>%
select(SPECIES_CD, SPECIES_NAME)
#Join w species names
den <- den %>% left_join(names, by = "SPECIES_CD")
```
```{r find top species, echo = FALSE}
# Load data to find top species
AR_ready_demo <- read_csv("../Project/data/FK2024_NCRMP_DRM_MIR_corsz_AR_NK.csv")
df <- AR_ready_demo %>%
filter(N == 1) %>%
group_by(SPECIES_CD) %>%
summarise(total = n(), .groups = "drop") %>%
arrange(desc(total)) %>%
#filter(SPECIES_CD %in% codes) %>%
slice(1:10)
topspecies <- df$SPECIES_CD
# Combine top species with MIR_species
target_species <- unique(c(topspecies, MIR_species))
```
```{r prep data, echo = FALSE}
tmp <- quad_panel_MIR(target_species)
# Prepare Density Data
a <- tmp$AR_density %>%
select(YEAR, SPECIES_CD, wavdns, se_wdns, analysis_group) %>%
rename(avDen = wavdns, SE_avDen = se_wdns)
# Prepare old mort data
b <- tmp$AR_mortality %>%
filter(MORT_TYPE == "Old") %>%
select(YEAR, SPECIES_CD, avMort, SE, analysis_group) %>%
rename(avMortOld = avMort, SE_avMortOld = SE)
# old_mortality_ready_for_stat_testing <- b %>%
# mutate(UCI = avMortOld + SE_avMortOld,
# LCI = avMortOld - SE_avMortOld)
#
# view(old_mortality_ready_for_stat_testing)
#
# readr::write_csv(old_mortality_ready_for_stat_testing, "old_mortality_ready_for_stat_testing.csv")
#
# dens_ready_for_stat_testing <- a%>%
# mutate(UCI = avDen + SE_avDen,
# LCI = avDen - SE_avDen)
#
#
# readr::write_csv(dens_ready_for_stat_testing, "dens_ready_for_stat_testing.csv")
#
#
# view(dens_ready_for_stat_testing)
# Prepare size data
# c <- tmp$AR_col_size %>%
# mutate(SE_maxdiam = sqrt(var_maxdiam)) %>%
# rename(avMaxdiam = avg_maxdiam) %>%
# select(REGION, YEAR, SPECIES_CD, avMaxdiam, SE_maxdiam, analysis_group)
#RATIO ESTIMATION
ratio_est_2022 <- read_csv("../Project/data/fk2022_NCRMP_MIR_sppLbar_dom.csv")
ratio_est_2024 <- read_csv("../Project/data/fk2024_NCRMP_MIR_sppLbar_domv2.csv")
c <- bind_rows(ratio_est_2022, ratio_est_2024) %>%
rename(avMaxdiam = Lbar,
SE_maxdiam = SE_L)
# Prep make size bins
tmp <- MIR_make_size_bins(years = c(2022, 2024),
analyzed_species = MIR_species)
#gather and clean domain Est
domain_estimates <- tmp$length_freq_domain_est %>%
dplyr::mutate(YEAR = as.factor(YEAR)) %>%
dplyr::ungroup()
#gather and clean length demos
demos <- tmp$length_demos %>%
dplyr::mutate(YEAR = as.factor(YEAR)) %>%
dplyr::ungroup()
#get domain mortality data
domain_mort <- tmp$domain_mort_spp
# write_csv(a, "../Project/denisty_data.csv")
# write_csv(b, "../Project/old_mort.csv")
# write_csv(c, "../Project/max_diam.csv")
```
```{r create_species_plots function, fig.width=9, fig.asp=1.2, echo = FALSE}
create_species_plots <- function(species) {
den_filtered <- a %>% filter(SPECIES_CD == species)
b_filtered <- b %>% filter(SPECIES_CD == species)
c_filtered <- c %>% filter(SPECIES_CD == species)
if (nrow(den_filtered) == 0 || nrow(b_filtered) == 0 || nrow(c_filtered) == 0) {
return(NULL)
}
dodge_width <- 0.05
### DENSITY PLOT ##########
plt1 <- den_filtered %>%
ggplot(aes(x = YEAR, y = avDen, group = analysis_group, color = analysis_group)) +
geom_point(position = position_dodge(dodge_width)) +
geom_line(position = position_dodge(dodge_width)) +
geom_errorbar(aes(ymin = avDen - SE_avDen, ymax = avDen + SE_avDen), width = 0.1, position = position_dodge(dodge_width)) +
labs(y = "Density (col/m2)", x = "Sample Year", title = "Coral Density") +
theme_Publication(base_size = 10) +
theme(strip.text = element_text(face = "italic"))+
theme(legend.position = "none") +
scale_x_continuous(breaks = c(2022, 2024)) +
scale_color_manual(values = c("#0072B2", "#D55E00"), labels = c("MIR", "Outside"), name = "") +
scale_linetype_manual(values = c("solid", "dashed"), guide = "none")
##Adding the significance labels
sig_cover <- sig %>%
filter(metric == "density") %>%
filter(variable == species)
# For between years find a midpoint to go on the plot
midpoints <- den_filtered %>%
group_by(analysis_group) %>%
summarise(
y_mid = 1.02 * mean(avDen, na.rm = TRUE)) %>%
dplyr::rename(REGION = analysis_group) #for joining purposes
##Special case for MMEA since there is overlap on the stars (they're in the same location). Easiest to rewrite
if (species == "MEA MEAN") {
midpoints <- tibble(
REGION = c("MIR_GRP", "NCRMP_GRP"),
y_mid = c(0.0011, 0.00187))}
## For between groups find the y min to go on the plot
ymin <- den_filtered %>%
group_by(SPECIES_CD) %>%
summarise(y_min = .85 * (min(avDen)- max(SE_avDen))) %>%#this hopefully gets a nice plot minimum just below wherever it ends
dplyr::rename(variable = SPECIES_CD) #for joining purposes
# Add them to the significance table
sig_cover <- sig_cover %>%
left_join(midpoints, by = "REGION") %>%
left_join(ymin, by = "variable")
##Plot it with annotations
plt1 <- plt1 +
### Significance BETWEEN GROUPS on X axis
geom_text(data = sig_cover %>% filter (is.na(REGION)),
aes(x = YEAR, y = y_min), label = "+",
inherit.aes = FALSE, size = 6, color = "black") +
### Significance BETWEEN YEARS on graph
geom_text(data = sig_cover %>% filter (is.na(YEAR)),
aes(x = 2023, y = y_mid, color = REGION), label = "*",
inherit.aes = FALSE, size = 10)
### MORTALITY PLOT ##########
#conditional APAL because only one year, add an NA for plotting aesthetics
if (species == "ACR PALM"){
b_filtered <- b_filtered %>%
bind_rows(
tibble(YEAR = 2024, SPECIES_CD = "ACR PALM", avMortOld = NA, SE_avMortOld = NA, analysis_group = "MIR_GRP"))
}
plt2 <- b_filtered %>%
ggplot(aes(x = YEAR, y = avMortOld, group = analysis_group, color = analysis_group)) +
geom_point(position = position_dodge(dodge_width)) +
geom_line(position = position_dodge(dodge_width)) +
geom_errorbar(aes(ymin = avMortOld - SE_avMortOld, ymax = avMortOld + SE_avMortOld), width = 0.1, position = position_dodge(dodge_width)) +
labs(y = "Old Mortality (%)", x = "Sample Year", title = "Old Mortality") +
theme_Publication(base_size = 10) +
theme(strip.text = element_text(face = "italic"))+
theme(legend.position = "none") +
scale_x_continuous(breaks = c(2022, 2024)) +
scale_color_manual(values = c("#0072B2", "#D55E00"), labels = c("MIR", "Outside"), name = "") +
scale_linetype_manual(values = c("solid", "dashed"), guide = "none")
#
##Adding the significance labels
#This may take a fair amount of tweaking ##
sig_cover <- sig %>%
filter(metric == "mortality") %>%
filter(variable == species)
# For between years find a midpoint to go on the plot
midpoints <- b_filtered %>%
group_by(analysis_group) %>%
summarise(
y_mid = 1.02 * mean(avMortOld, na.rm = TRUE)) %>%
dplyr::rename(REGION = analysis_group) #for joining purposes
## For between groups find the y min to go on the plot
ymin <- b_filtered %>%
mutate(YEARS = mean(YEAR)) %>%
group_by(SPECIES_CD, YEARS) %>%
summarise(y_min = .85 * (min(avMortOld)- max(SE_avMortOld))) %>%
#changed here because APAL only has 2022
dplyr::rename(variable = SPECIES_CD) #for joining purposes
# Add them to the significance table
sig_cover <- sig_cover %>%
left_join(midpoints, by = "REGION") %>%
left_join(ymin, by = "variable")
##Plot it with annotations
plt2 <- plt2 +
### Significance BETWEEN GROUPS on X axis
geom_text(data = sig_cover %>% filter (is.na(REGION)),
aes(x = YEAR, y = y_min), label = "+",
inherit.aes = FALSE, size = 6, color = "black") +
### Significance BETWEEN YEARS on graph
geom_text(data = sig_cover %>% filter (is.na(YEAR)),
aes(x = YEARS, y = y_mid, color = REGION), label = "*",
inherit.aes = FALSE, size = 10)
##### COLONY SIZE GRAPH ########
if (species == "ACR PALM"){
c_filtered <- c_filtered %>%
bind_rows(
tibble(YEAR = 2024, SPECIES_CD = "ACR PALM", avMortOld = NA, SE_avMortOld = NA, analysis_group = "MIR_GRP"))
}
plt3 <- c_filtered %>%
ggplot(aes(x = YEAR, y = avMaxdiam, group = analysis_group, color = analysis_group)) +
geom_point(position = position_dodge(dodge_width)) +
geom_line(position = position_dodge(dodge_width)) +
geom_errorbar(aes(ymin = avMaxdiam - SE_maxdiam, ymax = avMaxdiam + SE_maxdiam), width = 0.1, position = position_dodge(dodge_width)) +
labs(y = "Mean Max Diameter (cm)", x = "Sample Year", title = "Colony Size") +
theme_Publication(base_size = 10) +
theme(strip.text = element_text(face = "italic"))+
theme(legend.position = "none") +
scale_x_continuous(breaks = c(2022, 2024)) +
scale_color_manual(values = c("#0072B2", "#D55E00"), labels = c("MIR", "Outside"), name = "") +
scale_linetype_manual(values = c("solid", "dashed"), guide = "none")
##Adding the significance labels
#This may take a fair amount of tweaking ##
sig_cover <- sig %>%
filter(metric == "size") %>%
filter(variable == species)
# For between years find a midpoint to go on the plot
midpoints <- c_filtered %>%
group_by(analysis_group) %>%
summarise(
y_mid = 1.02 * mean(avMaxdiam), na.rm = TRUE) %>%
dplyr::rename(REGION = analysis_group) #for joining purposes
## For between groups find the y min to go on the plot
ymin <- c_filtered %>%
group_by(SPECIES_CD) %>%
summarise(y_min = .85 * (min(avMaxdiam)- max(SE_maxdiam))) %>%#this hopefully gets a nice plot minimum just below wherever it ends
dplyr::rename(variable = SPECIES_CD) #for joining purposes
# Add them to the significance table
sig_cover <- sig_cover %>%
left_join(midpoints, by = "REGION") %>%
left_join(ymin, by = "variable")
##Plot it with annotations
plt3 <- plt3 +
### Significance BETWEEN GROUPS on X axis
geom_text(data = sig_cover %>% filter (is.na(REGION)),
aes(x = YEAR, y = y_min), label = "+",
inherit.aes = FALSE, size = 6, color = "black") +
### Significance BETWEEN YEARS on graph
geom_text(data = sig_cover %>% filter (is.na(YEAR)),
aes(x = 2023, y = y_mid, color = REGION), label = "*",
inherit.aes = FALSE, size = 10)
domain_sub = subset(domain_estimates, SPECIES_CD == species)
demos_sub = subset(demos, SPECIES_CD == species)
mort_sub = domain_mort %>% dplyr::filter(SPECIES_CD == species)
domain_sub <- domain_sub %>%
# add 0s in for species not observed
tidyr::expand(., REGION, analysis_group, SPECIES_CD, YEAR, bin_num) %>%
#filter(YEAR %in% years) %>% #Filter back to years analyzed only
# connect back to demo_data to fill in with NAs
dplyr::full_join(., domain_sub) %>%
# make a presence/absence column
dplyr::mutate(length_freq_domain = ifelse(is.na(length_freq_domain), 0, length_freq_domain)) %>%
# add in mortality data
dplyr::full_join(., mort_sub)
n_bins = max(demos_sub$n_bins)
text_size = 14
angle = 45
hjust = dplyr::if_else(angle == 45, 1, 0.5)
min = min(demos_sub$min)
bin_width = min(demos_sub$bin_width)
small = dplyr::if_else(n_bins > 4, 0, 1)
num_vec = seq_len(max(domain_sub$bin_num))
lab_vec <- c("4- 10", " 11-15", " 16-20", " 21-25", " 26-30", " 31-35", " 36-45", " 46-65", " 66-85", " 86-105", "106+")
species_name <- names %>%
filter(SPECIES_CD == species) %>%
distinct(SPECIES_NAME) %>%
pull(SPECIES_NAME)
plt4 <- ggplot(data = domain_sub,
aes(x = as.integer(bin_num), y = length_freq_domain, fill = analysis_group)) +
geom_bar(stat = "identity", position = "dodge2", width = 0.9, color = "black", size = 0.5) +
geom_point(aes(y = oldmort_domain, color = analysis_group), size = 2) +
geom_line(aes(y = oldmort_domain, color = analysis_group), size = 2, alpha = 0.5) +
facet_wrap(~YEAR) +
labs(y = "Frequency (Size Class & Mortality)", x = "Length (cm)") +
scale_x_continuous(breaks = seq_along(lab_vec), labels = lab_vec) +
scale_fill_manual(values = c("#0072B2", "#D55E00"), labels = c("M:IR", "NCRMP + DRM")) +
scale_color_manual(values = c("#0072B2", "#D55E00"), guide = "none") +
theme_Publication(base_size = 10) +
theme(
axis.text.x = element_text(size = text_size, angle = angle, hjust = hjust),
legend.title = element_blank(),
strip.text = element_text(face = "italic"))
#what goes between the three plots and the size freq/mort plot
middle_annotation <- ggplot() +
geom_text(aes(x = 0.5, y = 0.5, label = "Size Frequency Distributions"), size = 5, fontface = "bold") +
theme_void() +
theme(plot.margin = margin(0, 0, 0, 0))
library(magick)
library(png)
library(jpeg)
species_images <- list(
"Acropora cervicornis" = list(path = "../Project/coral_image_files/staghorn.jpg", type = "jpg"),
"Colpophyllia natans" = list(path = "../Project/coral_image_files/col_nata.jpg", type = "jpg"),
"Dichocoenia stokesii" = list(path = "../Project/coral_image_files/dic_stok.jpg", type = "jpg"),
"Diploria labyrinthiformis" = list(path = "../Project/coral_image_files/dip_laby.jpg", type = "jpg"),
"Meandrina meandrites" = list(path = "../Project/coral_image_files/mea_mean.jpg", type = "jpg"),
"Pseudodiploria clivosa" = list(path = "../Project/coral_image_files/pse_cliv.png", type = "png"),
"Pseudodiploria strigosa" = list(path = "../Project/coral_image_files/pse_stri.png", type = "png"),
"Montastraea cavernosa" = list(path = "../Project/coral_image_files/mon_cave.jpg", type = "jpg"),
"Acropora palmata" = list(path = "../Project/coral_image_files/acr_palm.png", type = "png")
)
get_species_image_plot <- function(species_name) {
info <- species_images[[species_name]]
if (!is.null(info)) {
img <- switch(info$type,
png = png::readPNG(info$path),
jpg = jpeg::readJPEG(info$path))
cowplot::ggdraw() + cowplot::draw_image(img)
} else {
ggplot() + theme_void() }
}
image_plot <- get_species_image_plot(species_name)
final_plot <- (image_plot/
(plt1 + plt2 + plt3) /
middle_annotation /
plt4
) +
plot_layout(heights = c(1, 1, 0.1, 1.4)) +
plot_annotation(
title = paste0("<b><i>", species_name, "</i></b>"),
subtitle = "Colony density, mortality, size, and size class distribution (2022 vs. 2024)",
theme = theme(
plot.title = ggtext::element_markdown(hjust = 0.5, size = 18),
plot.subtitle = element_text(hjust = 0.5, size = 14) ))
return(final_plot)
}
plots_list <- map(MIR_species, create_species_plots)
names(plots_list) <- MIR_species
lst_clean <- compact(plots_list)
```
<!-- Tabs for HTML output -->
```{r coral-tabyss, echo=FALSE, fig.cap="Figure 3. Coral density, old mortality, colony size, and size frequency distribution for focal M:IR species across sample years (2022 and 2024) and survey types (M:IR and NCRMP + DRM). A black plus (+) indicates a significant difference in cover estimates between groups in that survey year (i.e., 2022 or 2024). A colored asterisk (*) indicates significance between years for the given group (i.e., M:IR or NCRMP). Note: Axis bounds vary by species.", results='asis', fig.width=9, fig.asp=1.1}
###
if (knitr::is_html_output()){
# Open the tabset container
cat("::: {.panel-tabset}\n\n")
# Loop through each species in lst_clean
purrr::walk(names(lst_clean), function(species_code) {
# # Skip if the plot is NULL
# if (is.null(lst_clean[[species_code]])) return()
# Get species name for the tab heading
species_name <- names %>%
dplyr::filter(SPECIES_CD == species_code) %>%
dplyr::distinct(SPECIES_NAME) %>%
dplyr::pull(SPECIES_NAME)
# Print the tab heading
cat("### *", species_name, "* {.tab-pane}\n\n", sep = "")
# Print the patchwork plot
print(lst_clean[[species_code]])
# Optional spacing between tabs
cat("\n\n")
})
# Close the tabset container
cat(":::\n")
} else if (knitr::is_latex_output()) {
# PDF: just print the plots, one per page
purrr::walk(lst_clean, function(p) {
if (is.null(p)) return()
print(p)
cat("\n\n")
})
}
```
### Coral Density
```{r overall Density, fig.cap = "Figure 4. Species-specific coral density (colones/m<sup>2</sup>) outside (NCRMP and DRM) and inside M:IR areas in 2022 and 2024. Focal M:IR species are denoted with a (*).", fig.width=9, fig.asp=1.1, echo = FALSE}
###Adding species names and starring MIR species
a <- left_join(a, names %>% distinct(), by = "SPECIES_CD") %>%
dplyr::mutate(SPECIES_NAME = case_when(
SPECIES_CD %in% MIR_species ~ paste0("*", SPECIES_NAME),
TRUE ~ as.character(SPECIES_NAME)))
#Factor re-order for aesthetics
a <- a %>%
mutate(SPECIES_NAME =
fct_reorder(SPECIES_NAME, if_else(is.na(avDen), 0, avDen), .desc = FALSE))
plt1 <- a %>%
ggplot(aes(x = avDen, y = SPECIES_NAME, group = interaction(SPECIES_CD, analysis_group), color = analysis_group)) +
geom_point( position = position_dodge(width = 0.3)) +
geom_errorbar(aes(xmin = avDen - SE_avDen, xmax = avDen + SE_avDen), width = 0.1, position = position_dodge(width = 0.3)) +
labs(y = "Species", x = expression(bold("Density (col/m"^2*")"))) +
theme_Publication(base_size = 10) +
facet_grid(~YEAR) +
scale_x_continuous(expand = expansion(mult = c(0, 0.1))) +
scale_y_discrete(labels = function(x) parse(text = paste0("italic('", x, "')")))+
theme(strip.text = element_text(face = "italic"))+
scale_color_manual(values = c("#0072B2", "#D55E00"), labels = c("MIR", "NCRMP + DRM"), name = "") +
scale_linetype_manual(values = c("solid", "dashed"), guide = "none") +
theme(
strip.text = element_text(size = 14, face = "bold"),
axis.text.y = element_text(size = 11),
axis.title = element_text(size = 13),
legend.position = "bottom",
legend.text = element_text(size = 10),
panel.spacing = unit(1, "lines")
)
print(plt1)
```
### Old Mortality
```{r overall Old_Mort, fig.cap = "Figure 5. Percent old mortality (± SEM) by coral species outside (NCRMP and DRM) and inside M:IR areas in 2022 and 2024. Focal M:IR species are denoted with a (*).", fig.width=9, fig.asp=1.1, echo = FALSE}
###Adding species names and starring MIR species
b <- left_join(b, names %>% distinct(), by = "SPECIES_CD") %>%
dplyr::mutate(SPECIES_NAME = case_when(
SPECIES_CD %in% MIR_species ~ paste0("*", SPECIES_NAME),
TRUE ~ as.character(SPECIES_NAME)))
#Factor re-order for aesthetics
b <- b %>%
mutate(SPECIES_NAME =
fct_reorder(SPECIES_NAME, if_else(is.na(avMortOld), 0, avMortOld), .desc = FALSE))
plt1 <- b %>%
ggplot(aes(x = avMortOld, y = SPECIES_NAME, group = interaction(SPECIES_CD, analysis_group), color = analysis_group)) +
geom_point( position = position_dodge(width = 0.3)) +
scale_y_discrete(labels = function(x) parse(text = paste0("italic('", x, "')")))+
geom_errorbar(aes(xmin = avMortOld - SE_avMortOld, xmax = avMortOld + SE_avMortOld), width = 0.1, position = position_dodge(width = 0.3)) +
labs(y = "Species", x = "Old Mortality (%)") +
theme_Publication(base_size = 10) +
facet_grid(~YEAR) +
scale_x_continuous(expand = expansion(mult = c(0, .08))) +
theme(strip.text = element_text(face = "italic"))+
scale_color_manual(values = c("#0072B2", "#D55E00"), labels = c("MIR", "NCRMP + DRM"), name = "") +
scale_linetype_manual(values = c("solid", "dashed"), guide = "none") +
theme(
strip.text = element_text(size = 14, face = "bold"),
axis.text.y = element_text(size = 11),
axis.title = element_text(size = 13),
legend.position = "bottom",
legend.text = element_text(size = 10),
panel.spacing = unit(1, "lines")
)
plt1
```
### Coral Size
```{r overall Max_Dim, fig.cap = "Figure 6. Maximum diameter (cm) ± SEM by coral species outside (NCRMP and DRM) and inside M:IR areas in 2022 and 2024. Focal M:IR species are denoted with a (*).", fig.width=9, fig.asp=1.1, echo = FALSE}
spec_unq <- unique(b$SPECIES_CD)
###Adding species names and starring MIR species
c <- left_join(c, names %>% distinct(), by = "SPECIES_CD") %>%
dplyr::mutate(SPECIES_NAME = case_when(
SPECIES_CD %in% MIR_species ~ paste0("*", SPECIES_NAME),
TRUE ~ as.character(SPECIES_NAME)))
#Factor re-order for aesthetics
c <- c %>%
mutate(SPECIES_NAME =
fct_reorder(SPECIES_NAME, if_else(is.na(avMaxdiam), 0, avMaxdiam), .desc = FALSE))
c <- c %>% filter(SPECIES_CD %in% spec_unq)
plt1 <- c %>%
ggplot(aes(x = avMaxdiam, y = SPECIES_NAME, group = interaction(SPECIES_CD, analysis_group), color = analysis_group)) +
geom_point( position = position_dodge(width = 0.3)) +
scale_y_discrete(labels = function(x) parse(text = paste0("italic('", x, "')")))+
geom_errorbar(aes(xmin = avMaxdiam - SE_maxdiam, xmax = avMaxdiam + SE_maxdiam), width = 0.1, position = position_dodge(width = 0.3)) +
labs(y = "Species", x = "Mean Colony Size (cm)") +
theme_Publication(base_size = 10) +
facet_grid(~YEAR) +
scale_x_continuous(limits = c(0, NA), breaks = scales::pretty_breaks(n = 6),
expand = expansion(mult = c(0, 0.08))) +
theme(strip.text = element_text(face = "italic"))+
scale_color_manual(values = c("#0072B2", "#D55E00"), labels = c("MIR", "Outside"), name = "") +
scale_linetype_manual(values = c("solid", "dashed"), guide = "none") +
theme(
strip.text = element_text(size = 14, face = "bold"),
axis.text.y = element_text(size = 11),
axis.title = element_text(size = 13),
legend.position = "bottom",
legend.text = element_text(size = 10),
panel.spacing = unit(1, "lines")
)
print(plt1)
```
### Bleaching Estimates
```{r bleaching setup, echo = FALSE}
MIR_species <- c("ACR PALM", "ACR CERV", "COL NATA", "DEN CYLI", "DIC STOK", "DIP LABY", "EUS FAST", "MEA MEAN", "PSE CLIV", "PSE STRI", "MON CAVE")
names_clean <- names %>%
select(SPECIES_CD, SPECIES_NAME) %>%
distinct(SPECIES_CD, .keep_all = TRUE)
##Old proportion estimates, will update when dione gives new ones
# tmp_ble_2022 <- read_csv("../Project/data/fk2022_NCRMP_MIR_sppPVbar_Ball_dom.csv") %>% left_join(., names_clean, by = "SPECIES_CD")
# tmp <- read_csv("../Project/data/fk2024_NCRMP_MIR_sppPVbar_Ball_dom.csv") %>%
# left_join(names_clean, by = "SPECIES_CD")
### Prevalence estimates given by Dione in place of above ones
tmp_ble_2022 <- read_csv("../Project/data/cor_prev_fk22_MIR_NCRMP.csv") %>% left_join(., names_clean, by = "SPECIES_CD") %>%
dplyr::rename(PVbar = percent_Bprev) ## change the column name for same code later
tmp <- read_csv("../Project/data/cor_prev_fk24_MIR_NCRMP.csv") %>% left_join(names_clean, by = "SPECIES_CD") %>%
dplyr::rename(PVbar = percent_Bprev) ## change the column name for same code later
tmp_ble <- full_join(tmp_ble_2022, tmp)
top_species <- tmp_ble %>%
group_by(SPECIES_CD) %>%
summarise(mean_pv = mean(PVbar, na.rm = TRUE)) %>%
arrange(desc(mean_pv)) %>%
slice_head(n = 10) %>%
pull(SPECIES_CD)
tmp_filtered <- tmp_ble %>%
filter(SPECIES_CD %in% union(top_species, MIR_species))
##Add a * for those MIR species
tmp_filtered <- tmp_filtered %>% mutate(
SPECIES_NAME = case_when(
SPECIES_CD %in% MIR_species ~ paste0("*", SPECIES_NAME),
TRUE ~ as.character(SPECIES_NAME)))
#Factor re-order for aesthetics
tmp_filtered <- tmp_filtered %>%
mutate(
SPECIES_NAME = fct_reorder(SPECIES_NAME, if_else(is.na(PVbar), 0, PVbar), .desc = FALSE),
YEAR = factor(YEAR, levels = c("2022", "2024")),
analysis_group = recode(analysis_group,
"MIR_GRP" = "MIR",
"NCRMP_GRP" = "NCRMP + DRM"))
plt_bleach_MIR <- tmp_filtered %>%
filter(SPECIES_NAME != "NA") %>%
filter(analysis_group == "MIR") %>%
ggplot(aes(x = PVbar, y = SPECIES_NAME, pattern = YEAR, group = rev(YEAR))) +
ggpattern::geom_col_pattern(
position = position_dodge(width = 0.8),
width = 0.8,
fill = "#0072B2", # background color (used if no pattern)
pattern_fill = "white", # color of the pattern lines (overlay)
pattern_spacing = 0.0075, # very tight pattern
pattern_size = 0.025, # thin lines
pattern_density = .4 # fully patterned
) +
scale_pattern_manual(values = c("2022" = "none", "2024" = "stripe")) +
# geom_errorbar(aes(xmin = pmax(PV_LCI, 0), xmax = PV_UCI), position = position_dodge(width = 0.8), width = 0.1) +
labs(x = "Bleaching Prevalence (%)", y = "Species", fill = "") +
scale_x_continuous(expand = expansion(mult = c(0, 0.1))) +
theme_Publication(base_size = 10) +
coord_cartesian(xlim = c(0, 80), clip = "off") +
scale_y_discrete(labels = function(x) parse(text = paste0("italic('", x, "')"))) +
theme(
strip.text = element_text(size = 14, face = "bold"),
axis.text.y = element_text(size = 11),
axis.title = element_text(size = 13),
legend.position = "bottom",
title = element_text(size = 14),
legend.text = element_text(size = 10),
panel.spacing = unit(1, "lines")
) +
ggtitle("Mission: Iconic Reef")
plt_bleach_NCRMP <- tmp_filtered %>%
filter(SPECIES_NAME != "NA") %>%
filter(analysis_group == "NCRMP + DRM") %>%
ggplot(aes(x = PVbar, y = SPECIES_NAME, pattern = YEAR, group = rev(YEAR))) +
ggpattern::geom_col_pattern(
position = position_dodge(width = 0.8),
width = 0.8,
fill = "#D55E00", # background color (used if no pattern)
pattern_fill = "white", # color of the pattern lines (overlay)
pattern_spacing = 0.0075, #75 # very tight pattern
pattern_size = 0.025, # thin lines
pattern_density = .4 # fully patterned
) +
scale_pattern_manual(values = c("2022" = "none", "2024" = "stripe")) +
# geom_errorbar(aes(xmin = pmax(PV_LCI, 0), xmax = PV_UCI), position = position_dodge(width = 0.8), width = 0.1) +
# facet_wrap(~analysis_group) +
labs(x = "Bleaching Prevalence (%)", y = "Species", fill = "") +
scale_x_continuous(expand = expansion(mult = c(0, 0.1))) +
theme_Publication(base_size = 10) +
coord_cartesian(xlim = c(0, 80), clip = "off") +
scale_y_discrete(labels = function(x) parse(text = paste0("italic('", x, "')"))) +
theme(
strip.text = element_text(size = 14, face = "bold"),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.title = element_text(size = 13),
title = element_text(size = 14),
legend.position = "bottom",
legend.text = element_text(size = 10),
panel.spacing = unit(1, "lines")
) +
ggtitle("NCRMP + DRM")
```
```{r bleaching, fig.cap = "Figure 7. Species-specific bleaching frequency observed outside (NCRMP + DRM) and inside M:IR areas. Focal M:IR species are denoted with a (*).", fig.width=10, fig.asp=1.1, echo = FALSE}
# Combine the two plots side-by-side
final_plot <- (plt_bleach_MIR + plt_bleach_NCRMP) +
plot_layout(ncol = 2) # Ensures side-by-side layout
print(final_plot)
```
### Disease Estimates
```{r disease, echo= FALSE}
# Set path to data files
MIR_species <- c("ACR PALM", "ACR CERV", "COL NATA", "DEN CYLI", "DIC STOK", "DIP LABY", "EUS FAST", "MEA MEAN", "PSE CLIV", "PSE STRI", "MON CAVE")
names_clean <- names %>%
select(SPECIES_CD, SPECIES_NAME) %>%
distinct(SPECIES_CD, .keep_all = TRUE)
#old proportional estimates, will update when new ones come in from Dione
# tmp_dis_2022 <- read_csv("../Project/data/fk2022_NCRMP_MIR_sppPVbar_Dall_dom.csv") %>%
# left_join(names_clean, by = "SPECIES_CD")
#
# tmp_dis_2024 <- read_csv("../Project/data/fk2024_NCRMP_MIR_sppPVbar_Dall_dom.csv") %>%
# left_join(names_clean, by = "SPECIES_CD")
#Use prevalence dataset and keep names the same for following code
tmp_dis <- tmp_ble %>% select(-PVbar) %>% dplyr::rename(PVbar = percent_Dprev)
# tmp_dis <- full_join(tmp_dis_2022, tmp_dis_2024)
top_species <- tmp_dis %>%
group_by(SPECIES_CD) %>%
summarise(mean_pv = mean(PVbar, na.rm = TRUE)) %>%
arrange(desc(mean_pv)) %>%
slice_head(n = 10) %>%
pull(SPECIES_CD)
tmp_filtered <- tmp_dis %>%
filter(SPECIES_CD %in% union(top_species, MIR_species))
##Add a * for those MIR species
tmp_filtered <- tmp_filtered %>% mutate(
SPECIES_NAME = case_when(
SPECIES_CD %in% MIR_species ~ paste0( "*", SPECIES_NAME),
TRUE ~ as.character(SPECIES_NAME))
)
#Factor re-order for aesthetics
tmp_filtered <- tmp_filtered %>%
mutate(
SPECIES_NAME = fct_reorder(SPECIES_NAME, PVbar, .fun = function(x) max(x, na.rm = TRUE), .desc = FALSE),
YEAR = factor(YEAR, levels = c("2022", "2024")),
analysis_group = recode(analysis_group,
"MIR_GRP" = "MIR",
"NCRMP_GRP" = "NCRMP + DRM"))
plt_dis_MIR <- tmp_filtered %>%
filter(SPECIES_NAME != "NA") %>%
filter(analysis_group == "MIR") %>%
ggplot(aes(x = PVbar, y = SPECIES_NAME, group = rev(YEAR), pattern = YEAR)) +
ggpattern::geom_col_pattern(
position = position_dodge(width = 0.8),
width = 0.8,
fill = "#0072B2",
pattern_fill = "white",
pattern_spacing = 0.0075,
pattern_size = 0.025,
pattern_density = .4
) +
scale_pattern_manual(values = c("2022" = "none", "2024" = "stripe")) +
# geom_errorbar(aes(xmin = pmax(PV_LCI, 0), xmax = PV_UCI),
# position = position_dodge(width = 0.8), width = 0.1) +
labs(x = "Disease Prevalence", y = "Coral Species", fill = "") +
scale_x_continuous(expand = expansion(mult = c(0, 0.1))) +
theme_Publication(base_size = 10) +
coord_cartesian(xlim = c(0, 10), clip = "off") +
scale_y_discrete(labels = function(x) parse(text = paste0("italic('", x, "')"))) +
theme(
strip.text = element_text(size = 14, face = "bold"),
axis.text.y = element_text(size = 11),
axis.title = element_text(size = 13),
legend.position = "bottom",
title = element_text(size = 14),
legend.text = element_text(size = 10),
panel.spacing = unit(1, "lines")
) +
ggtitle("Mission: Iconic Reef")
plt_dis_NCRMP <- tmp_filtered %>%
filter(SPECIES_NAME != "NA") %>%
filter(analysis_group == "NCRMP + DRM") %>%
ggplot(aes(x = PVbar, y = SPECIES_NAME, group = rev(YEAR), pattern = YEAR)) +
ggpattern::geom_col_pattern(
position = position_dodge(width = 0.8),
width = 0.8,
fill = "#D55E00",
pattern_fill = "white",
pattern_spacing = 0.0075,
pattern_size = 0.025,
pattern_density = .4
) +
scale_pattern_manual(values = c("2022" = "none", "2024" = "stripe")) +
# geom_errorbar(aes(xmin = pmax(PV_LCI, 0), xmax = PV_UCI),
# position = position_dodge(width = 0.8), width = 0.1) +
labs(x = "Disease Prevalence", y = "Coral Species", fill = "") +
scale_x_continuous(expand = expansion(mult = c(0, 0.1))) +
theme_Publication(base_size = 10) +
coord_cartesian(xlim = c(0, 10), clip = "off") +
scale_y_discrete(labels = function(x) parse(text = paste0("italic('", x, "')"))) +
theme(
strip.text = element_text(size = 14, face = "bold"),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.title = element_text(size = 13),
title = element_text(size = 14),
legend.position = "bottom",
legend.text = element_text(size = 10),
panel.spacing = unit(1, "lines")
) +
ggtitle("NCRMP + DRM")
```
```{r, echo = FALSE, fig.cap = "Figure 8. Species-specific disease frequency observed outside (NCRMP + DRM) and inside M:IR areas. Focal M:IR species are denoted with a (*).", fig.width=10, fig.asp=1.1}
# Combine the two plots side-by-side
(final_plot <- (plt_dis_MIR + plt_dis_NCRMP) +
plot_layout(ncol = 2))
```
## Point of Contact
*NCRMP Atlantic Benthic Team Lead:* **Shay Viehman** ([shay.viehman\@noaa.gov](mailto:shay.viehman@noaa.gov))
## References
::: {#refs}
:::
</br>
**Data Sources:**
NOAA National Centers for Coastal Ocean Science; NOAA Southeast Fisheries Science Center (2024). National Coral Reef Monitoring Program: Assessment of coral reef benthic communities in Florida \[2022, 2024\]. NOAA National Centers for Environmental Information. Dataset. doi: [10.7289/v5xw4h4z](https://doi.org/10.7289/v5xw4h4z)
Florida Fish and Wildlife Conservation Commission's Disturbance Response Monitoring ([DRM](https://coraldrm.org/))
