--- title: "Deep Red Sea Fishing" date: "2026-04-05" output: html_document editor_options: chunk_output_type: console --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` # libraries Load ```{r} library(tidyverse) library(compositions) library(zoo) library(ggrepel) library(cowplot) ``` # Prepare Data Load all LMEs data, bind, clean, and save ```{r} files <- list.files("Fisheries_Data", "\\.csv$", full.names = TRUE) df_list <- lapply(files, read.csv) df <- do.call(rbind, df_list) rm(files, df_list) ``` Check functional groups for uninformative ones ```{r} # Proportion of non-zero cells per functional group test <- df %>% group_by(functional_group) %>% summarise(prop_nonzero = mean(tonnes > 0), n_lmes = n_distinct(area_name[tonnes > 0]), first_nonzero = min(year[tonnes > 0])) ## drop krill (only appears in one LME) rm(test) ``` Aggregate to LME x year x functional group ```{r} # Drop Krill (1 LME only, appears 2001) # Sum across all fishing sectors and gear types catch_fg <- df %>% filter(functional_group != "Krill") %>% group_by(area_name, year, functional_group) %>% summarise(tonnes = sum(tonnes, na.rm = TRUE), .groups = "drop") ``` Artisanal pressure index ```{r} # Coastal small-scale = Artisanal + Subsistence # Trailing 5-year rolling mean (align = "right") # Loses 1950-1953, retains through 2019 pressure_index <- df %>% filter(fishing_sector %in% c("Artisanal", "Subsistence")) %>% group_by(area_name, year) %>% summarise(coastal_tonnes = sum(tonnes, na.rm = TRUE), .groups = "drop") %>% arrange(area_name, year) %>% group_by(area_name) %>% mutate( coastal_smooth = rollmean(coastal_tonnes, k = 5, fill = NA, align = "right"), peak_coastal = max(coastal_smooth, na.rm = TRUE), pressure_idx = coastal_smooth / peak_coastal, peak_year = year[which.max(coastal_smooth)] ) %>% ungroup() rm(df) ``` Wide format and CLR transformation ```{r} fg_cols <- catch_fg %>% distinct(functional_group) %>% pull(functional_group) # Pivot wide, fill structural zeros catch_wide <- catch_fg %>% pivot_wider(names_from = functional_group, values_from = tonnes, values_fill = 0) rm(catch_fg) # Proportions within each LME-year catch_prop <- catch_wide %>% mutate( total = rowSums(across(all_of(fg_cols))), across(all_of(fg_cols), ~ . / total) ) %>% select(-total) rm(catch_wide) # Pseudocount = half the minimum observed non-zero proportion pseudocount <- catch_prop %>% select(all_of(fg_cols)) %>% pivot_longer(everything()) %>% filter(value > 0) %>% pull(value) %>% min() / 2 # Apply pseudocount and CLR transform clr_matrix <- catch_prop %>% mutate(across(all_of(fg_cols), ~ ifelse(. == 0, pseudocount, .))) %>% select(all_of(fg_cols)) %>% as.matrix() %>% clr() # Trailing 5-year rolling mean on CLR values clr_smooth <- bind_cols( catch_prop %>% select(area_name, year), as_tibble(clr_matrix)) %>% arrange(area_name, year) %>% group_by(area_name) %>% mutate(across(all_of(fg_cols), ~ rollmean(., k = 5, fill = NA, align = "right"))) %>% ungroup() %>% drop_na() rm(catch_prop, clr_matrix, pseudocount) ``` # Analysis: PCA Run PCA on Centred Log-Ratio data ```{r} pca_result <- prcomp( clr_smooth %>% select(all_of(fg_cols)), center = TRUE, scale. = FALSE # CLR already on comparable scale ) var_explained <- summary(pca_result)$importance[2, ] # PC scores with metadata pca_scores <- bind_cols( clr_smooth %>% select(area_name, year), as_tibble(pca_result$x)) %>% mutate(is_red_sea = area_name == "Red Sea") # Loadings (keep PC1-PC4 for inspection) loadings <- as_tibble(pca_result$rotation, rownames = "functional_group") %>% select(functional_group, PC1, PC2, PC3, PC4) rm(clr_smooth, pca_result) ``` # Analysis: Red Sea global positioning Global averages across years for PCA1 and PCA2? ```{r} library(tidyverse) df_decade <- pca_scores %>% mutate(decade = (year %/% 10) * 10) %>% # assigns to decade (e.g. 1990, 2000) group_by(decade) %>% summarise( PC1_mean = mean(PC1, na.rm = TRUE), PC1_se = sd(PC1, na.rm = TRUE) / sqrt(sum(!is.na(PC1))), PC2_mean = mean(PC2, na.rm = TRUE), PC2_se = sd(PC2, na.rm = TRUE) / sqrt(sum(!is.na(PC2))), n = n(), .groups = "drop" ) %>% filter(decade %in% c(1950, 1970, 1990, 2010)) ``` Standard ellipse display in PC1-PC2 space ```{r} # Each LME = one ellipse computed across all years # Red Sea highlighted — expect tight ellipse at negative PC1 extreme lme_centroids <- pca_scores %>% group_by(area_name, is_red_sea) %>% summarise(pc1_mean = mean(PC1), pc2_mean = mean(PC2), .groups = "drop") ## Red sea as a percentile of the mean of all my_ecdf <- ecdf(lme_centroids$pc1_mean) percentile_rank <- my_ecdf(-36.1552381) * 100 # Output: 60 (meaning 35 is at the 60th percentile) ``` Combine ellipses with loadings ```{r} # Scaling factor for PCA loadings to make them visible scale_factor <- min( (max(pca_scores$PC1) - min(pca_scores$PC1)) / (max(loadings$PC1) - min(loadings$PC1)), (max(pca_scores$PC2) - min(pca_scores$PC2)) / (max(loadings$PC2) - min(loadings$PC2))) * 0.9 # Prepare loadings cutoff so only most important are labelled label_cutoff <- 1 / sqrt(29) loadings2 <- loadings %>% mutate( mag = sqrt(PC1^2 + PC2^2), strong = mag >= label_cutoff, label_clean = gsub("\\s*\\([^\\)]+\\)$", "", functional_group), arrow_col = ifelse(strong, "steelblue", "grey70"), arrow_alpha = ifelse(strong, 0.8, 0.25) ) # Plot p01 <- ggplot(pca_scores, aes(x = PC1, y = PC2, group = area_name, colour = is_red_sea, fill = is_red_sea)) + # Ellipses stat_ellipse(aes(alpha = is_red_sea), geom = "polygon", level = 0.99, type = "norm", colour = "transparent", fill = "grey75", alpha = .25) + stat_ellipse(data = pca_scores %>% filter(is_red_sea), geom = "polygon", level = 0.99, fill = "red", colour = "darkred", alpha = .5, linewidth = .5, type = "norm") + # Red Sea Centroid label geom_text(data = lme_centroids %>% filter(is_red_sea), aes(x = pc1_mean + 6.5, y = pc2_mean + 2.5, label = area_name), colour = "black", fontface = "bold", size = 6 / .pt, inherit.aes = FALSE) + # Loadings (WEAK arrows) geom_segment(data = loadings2 %>% filter(!strong), aes(x = 0, y = 0, xend = PC1 * scale_factor, yend = PC2 * scale_factor), arrow = arrow(length = unit(0.1, "cm")), linewidth = .25, colour = "steelblue", alpha = 0.4, inherit.aes = FALSE) + # Loadings (STRONG arrows) geom_segment(data = loadings2 %>% filter(strong), aes(x = 0, y = 0, xend = PC1 * scale_factor, yend = PC2 * scale_factor), arrow = arrow(length = unit(0.1, "cm")), linewidth = .5, colour = "steelblue", alpha = 0.8, inherit.aes = FALSE) + # Loadings labels (ONLY strong) geom_text_repel(data = loadings2 %>% filter(strong), aes(x = PC1 * scale_factor, y = PC2 * scale_factor, label = label_clean), box.padding = 0.1, size = 5 / .pt, max.overlaps = 5, colour = "steelblue", inherit.aes = FALSE) + # Zero lines geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.3, linewidth = .25) + geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.3, linewidth = .25) + # Scales scale_colour_manual(values = c("grey50", "red")) + scale_fill_manual(values = c("grey70", "red")) + scale_alpha_manual(values = c(0.08, 0.5)) + scale_y_continuous(limits = c(-39, 46)) + scale_x_continuous(limits = c(-49.5, 55.5)) + # Theme coord_equal() + labs(x = paste0("PC1 (", round(var_explained[1] * 100, 1), "% of Variation)"), y = paste0("PC2 (", round(var_explained[2] * 100, 1), "% of Variation)")) + guides(colour = "none", fill = "none", alpha = "none") + theme_classic() + theme(text = element_text(size = 6), axis.line = element_blank(), axis.ticks = element_line(linewidth = 0.2)) p01 ``` ## Analysis: Artisanal Fisheries Pressure in the Red Sea Show the change in artisanal fisheries pressure over time in the Red Sea. ```{r} red_sea_pressure <- pressure_index %>% filter(area_name == "Red Sea") red_sea_peak_year <- red_sea_pressure %>% filter(!is.na(coastal_smooth)) %>% filter(coastal_smooth == max(coastal_smooth, na.rm = TRUE)) %>% pull(year) red_sea_peak_val <- red_sea_pressure %>% filter(year == red_sea_peak_year) %>% pull(coastal_smooth) p02 <- red_sea_pressure %>% filter(!is.na(coastal_smooth)) %>% ggplot(aes(x = year, y = coastal_smooth)) + geom_ribbon(data = red_sea_pressure %>% filter(!is.na(coastal_smooth), year >= red_sea_peak_year), aes(ymin = 0, ymax = coastal_smooth), fill = "red", alpha = 0.15) + geom_line(colour = "darkred", linewidth = .5) + geom_vline(xintercept = red_sea_peak_year, linetype = "dashed", colour = "darkred", alpha = 0.7, linewidth = .2) + annotate("text", x = red_sea_peak_year - 1, y = red_sea_peak_val * 1.02, label = paste0("Estimated peak: ", red_sea_peak_year), hjust = 1.1, vjust = 1, colour = "black", size = 6 / .pt) + scale_x_continuous(breaks = c(1960, 1980, 2000, 2020)) + scale_y_continuous(expand = c(0, NA), limits = c(0,120000)) + labs(x = "Year", y = "Artisanal and subsistence catch (tonnes)") + theme_classic() + theme(text = element_text(size = 6), axis.line.x = element_blank(), axis.ticks = element_line(linewidth = 0.2)) p02 ``` # Analysis: Change since peak artisanal PC1 trajectory post-peak for LMEs compositionally similar to the Red Sea ```{r} all_lme_combined <- pca_scores %>% select(area_name, year, PC1, PC2, is_red_sea) %>% left_join( pressure_index %>% select(area_name, year, coastal_smooth, pressure_idx, peak_year), by = c("area_name", "year") ) # Red Sea PC1 value at its artisanal peak year red_sea_pc1_at_peak <- all_lme_combined %>% filter(area_name == "Red Sea", year == red_sea_peak_year) %>% pull(PC1) # Each LME's PC1 value at its own peak year pc1_at_peak <- all_lme_combined %>% filter(!is.na(PC1), !is.na(peak_year)) %>% group_by(area_name) %>% filter(year == first(peak_year)) %>% summarise(pc1_at_peak = mean(PC1), peak_year = first(peak_year), .groups = "drop") # Filter to LMEs within +-35 PC1 units of Red Sea at peak # and peak >= 5 years before end of series similar_lmes <- pc1_at_peak %>% filter( abs(pc1_at_peak - red_sea_pc1_at_peak) <= 35, peak_year <= 2014 ) %>% pull(area_name) # Post-peak trajectories for similar LMEs post_peak_similar <- all_lme_combined %>% filter( area_name %in% similar_lmes, !is.na(pressure_idx), !is.na(PC1) ) %>% group_by(area_name) %>% mutate(years_since_peak = year - peak_year) %>% filter(years_since_peak >= 0) %>% ungroup() rm(all_lme_combined, red_sea_pc1_at_peak, pc1_at_peak, similar_lmes, pressure_index, red_sea_peak_year, pca_scores, loadings, fg_cols) ``` PC1 spaghetti ```{r} p03 <- ggplot() + geom_line( data = post_peak_similar %>% filter(!is_red_sea), aes(x = years_since_peak, y = PC1, group = area_name), colour = "grey75", alpha = 0.4, linewidth = 0.25 ) + geom_smooth( data = post_peak_similar %>% filter(!is_red_sea), aes(x = years_since_peak, y = PC1, group = area_name), method = "gam", formula = y ~ s(x, k = 3), colour = "grey75", linewidth = .5, se = FALSE ) + geom_line( data = post_peak_similar %>% filter(is_red_sea), aes(x = years_since_peak, y = PC1), colour = "red", linewidth = .5 ) + geom_point( data = post_peak_similar %>% filter(is_red_sea) %>% filter(years_since_peak == max(years_since_peak)), aes(x = years_since_peak, y = PC1), colour = "darkred", size = .5 ) + annotate("text", x = post_peak_similar %>% filter(is_red_sea) %>% pull(years_since_peak) %>% max() + 0.5, y = post_peak_similar %>% filter(is_red_sea) %>% filter(years_since_peak == max(years_since_peak)) %>% pull(PC1), label = "Red Sea", colour = "black", hjust = -.1, size = 6 / .pt, fontface = "bold") + labs( x = "Years since artisanal and subsistence catch peak", y = paste0("PC1 (", round(var_explained[1] * 100, 1), "%)") ) + theme_classic() + theme(text = element_text(size = 6), axis.line.x = element_blank(), axis.ticks = element_line(linewidth = 0.2)) p03 ``` # Plot Build combined plot ```{r} pFC <- ggdraw() + draw_plot(p01, x = 0, y = 0, width = 2/3, height = 1) + draw_plot(plot_grid(p02, p03, ncol = 1, align = "hv"), x = 2/3, y = 0, width = 1/3, height = 1) + draw_plot_label(c("a", "b", "c"), x = c(.01, 2/3 - .01, 2/3 - .01), y = c(.99, .99, .49), size = 8) pFC ggsave("RS_Fisheries.png", pFC, width = 18, height = 9.7, units = "cm", dpi = 600, bg = "white") ```