What is a Bodymass Size Spectrum?

A community size spectrum is the relationship between abundance and body size. This shape and structure of this relationship has been shown to be temporally stable and reflect characteristics of the environments that support it (productivity levels).

This relationship is often plotted on a log body size x log abundance scale which reveals a linear relationship between the two. The slope and intercept of this relationship have been used to make comparisons across ecosystems regarding their productivity and the efficiency of energy transfer within a community.

Regions of Interest

# Load the shapefiles for the polygons
new_england <- ne_states("united states of america", returnclass = "sf")
canada <- ne_states("canada", returnclass = "sf")

# Load the strata
survey_strata <- read_sf(str_c(res_path, "Shapefiles/BottomTrawlStrata/BTS_Strata.shp"))  %>% 
  janitor::clean_names() %>% 
  filter(strata >= 01010 ,
         strata <= 01760,
         strata != 1310,
         strata != 1320,
         strata != 1330,
         strata != 1350,
         strata != 1410,
         strata != 1420,
         strata != 1490) 


# Key to which strata = which regions
strata_key <- list(
  "Georges Bank"          = as.character(13:23),
  "Gulf of Maine"         = as.character(24:40),
  "Southern New England"  = str_pad(as.character(1:12), width = 2, pad = "0", side = "left"),
  "Mid-Atlantic Bight"    = as.character(61:76))


# Assign Areas
survey_strata <- survey_strata %>% 
  mutate(
    strata = str_pad(strata, width = 5, pad = "0", side = "left"),
    strata_num = str_sub(strata, 3, 4),
    area = case_when(
      strata_num %in% strata_key$`Georges Bank` ~ "Georges Bank",
      strata_num %in% strata_key$`Gulf of Maine` ~ "Gulf of Maine",
      strata_num %in% strata_key$`Southern New England` ~ "Southern New England",
      strata_num %in% strata_key$`Mid-Atlantic Bight` ~ "Mid-Atlantic Bight",
    TRUE ~ "Outside Major Study Areas"),
    area = factor(area, levels = c("Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight"))) %>% 
  select(finstr_id, strata, strata_num, area, a2, str2, set, stratuma, str3, geometry)


# Map it
strata_plot <- ggplot() +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  geom_sf(data = survey_strata, aes(fill = area), alpha = 0.8) +
  coord_sf(xlim = c(-77, -65), ylim = c(34, 45.75), expand = FALSE) +
  guides(fill = guide_legend(nrow = 2)) +
  scale_fill_gmri() +
  scale_x_continuous(breaks = seq(-180, 180, by =5)) +
  theme(legend.position = c(0.7, 0.125), legend.title = element_blank()) +
  labs(title = "Study Area",
       fill = "Survey Regions")

strata_plot

Comparing Spectra Parameters

# Load the spectra coefficients
withr::with_dir(rprojroot::find_root('_targets.R'), 
                tar_load(size_spectrum_indices))

# Pull the ones that are regions and years
region_dat <- size_spectrum_indices %>% 
  filter(`group ID` == "single years * region") %>% mutate(Year = as.numeric(Year)) %>% 
  mutate(survey_area = factor(survey_area, levels = c("GoM", "GB", "MAB", "SNE")))

# Put them in list for lata
region_indices <- region_dat %>% 
  split(.$survey_area)



slopes <- region_dat %>% 
  ggplot(aes(Year, l10_slope_strat)) +
  geom_line(size = 0.5, color = "gray50") +
  geom_point(size = 1,color = "gray10") +
  facet_wrap(~survey_area, ncol = 1) +
  labs(title = "Size Spectrum Slope", y = "Slope")

ints <- region_dat %>% 
  ggplot(aes(Year, l10_int_strat)) +
  geom_line(size = 0.5, color = "gray50") +
  geom_point(size = 1,color = "gray10") +
  facet_wrap(~survey_area, ncol = 1) +
  labs(title = "Size Spectrum Intercept", y ="Intercept")

slopes | ints

Visualizing the Spectra

The body size spectra were estimated for each of these regions independently. Abundances were area-stratified, then grouped into log bodymass bins. 3d Surfaces below display how the shape of the spectra change over time.

These surfaces are followed by the resulting size spectrum parameters of the annual size spectra.

# Split the data into the groups we want to compare:
areas_binned <- nefsc_1g_binned %>% 
  split(.$survey_area)


# Get annual totals within each area
areas_aggregated <- areas_binned %>% 
  map(function(area_data){
    area_data %>% 
      split(.$Year) %>% 
      map_dfr(function(yr_data){
        aggregate_l10_bins(yr_data) %>%
          mutate(log10_bins = as.numeric(log10_bins),
                 norm_strat_abund = log10(norm_strat_abund),
                 norm_strat_abund = ifelse(is.na(norm_strat_abund), 0, norm_strat_abund))
          
      }, .id = "Year")
  })
####  Surface Plots
spectra_surface <- function(l10_totals, plot_title){
  
  # Pull data for matrix
  dat <- select(l10_totals, left_lim, Year, norm_strat_abund)
  
  # Make matrix to plot
  spectrum_mat <- as.matrix(xtabs(norm_strat_abund~left_lim + Year, data = dat))

  # Build plotly plot
  surf_plot <- plot_ly(
      x = as.numeric(colnames(spectrum_mat)), 
      y = as.numeric(rownames(spectrum_mat)), 
      z = spectrum_mat
    ) %>% 
  add_surface() %>%
  layout(
    title = str_c(plot_title, " - Bodymass Spectrum"),
    scene = list(
      xaxis = list(title = "Year", autorange = "reversed"),
      yaxis = list(title = "Log10( BodyMass (g) )"),
      zaxis = list(title = "Log10( Abundance) "),
      camera = list(eye = list(x = 1, y = 1, z = 1.25))
    ))
  
  #return plot
  return(surf_plot)
  
}


# Spectrum parameter plots
plot_params <- function(area_abbrev){
  slope_p <- ggplot(region_indices[[area_abbrev]], aes(Year, l10_slope_strat)) +
    geom_line() +
    geom_point() +
    labs(y = "Size Spectrum Slope", x = "")
  
  int_p <- ggplot(region_indices[[area_abbrev]], aes(Year, l10_int_strat)) +
    geom_line() +
    geom_point() +
    labs(y = "Size Spectrum Intercept")
  
  
  bin_data <- areas_aggregated[[area_abbrev]] %>% 
    mutate(Year = as.numeric(Year),
           l10_bins = factor(left_lim))
  bin_p <- ggplot(bin_data, aes(x = Year, y = norm_strat_abund, color = l10_bins)) +
    geom_line(aes(group = log10_bins)) +
    geom_point() + 
    scale_color_viridis_d(direction = -1) +
    guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) +
    labs(x = "Year", y = "Log( Abundance )", color = "Body Size Bin log10(g)")
  
  layout <- "
  AA##
  AACC
  BBCC
  BB##
  "
  
  # title
  tidy_name <- switch (area_abbrev,
    "GoM" = "Gulf of Maine",
    "GB" = "Georges Bank",
    "SNE" = "Southern New England",
    "MAB" = "Mid-Atlantic Bight"
  )
  
  plot_mosaic <- (slope_p / int_p)# / bin_p
  plot_mosaic <- plot_mosaic + plot_annotation(title = tidy_name)
  return( plot_mosaic)
}

Regional Comparisons

# All Data, Just comparing Regions
allyr_slopes <- region_dat %>% 
  ggplot(aes(survey_area, l10_slope_strat, color = survey_area)) +
  ggdist::stat_halfeye(
     aes(fill = survey_area),
     adjust = .5, 
     width = .6, 
     .width = 0, 
     justification = -.2, 
     point_colour = NA) + 
  geom_boxplot(
    width = .15, 
    outlier.shape = NA) +
  geom_point(
    size = 1.3,
    alpha = .3,
    position = position_jitter(
      seed = 1, width = .1)) +
  scale_color_gmri() +
  scale_fill_gmri() +
  scale_y_continuous(limits = c(-1.65, -0.6)) +
  facet_wrap(~survey_area, nrow = 1, scales = "free_x") +
  theme(legend.position = "none") +
  labs(x = "", y = "Size Spectrum Slope", color = "Time Period", 
       title =  "Regional Comparison: Size Spectra Slopes (all years)")


allyr_slopes

# Early year slopes

early_yr_slopes <-  region_dat %>% 
  # mutate(recent_yr = ifelse(Year >= 2000, "2000 +", "< 2000")) %>% 
  filter(Year < 2000) %>% 
  ggplot(aes(survey_area, l10_slope_strat, color = survey_area)) +
  ggdist::stat_halfeye(
     aes(fill = survey_area),
     adjust = .5, 
     width = .6, 
     .width = 0, 
     justification = -.2, 
     point_colour = NA) + 
  geom_boxplot(
    width = .15, 
    outlier.shape = NA) +
  geom_point(
    size = 1.3,
    alpha = .3,
    position = position_jitter(
      seed = 1, width = .1
    )
  ) +
  scale_color_gmri() +
  scale_fill_gmri() +
   scale_y_continuous(limits = c(-1.65, -0.6)) +
  facet_wrap(~survey_area, nrow = 1, scales = "free_x") +
  theme(legend.position = "none") +
  labs(x = "", y = "Size Spectrum Slope", color = "Time Period", 
       title =  "Regional Comparison: Size Spectra Slopes (before 2000)")


early_yr_slopes

# Recent Year Comparison
recent_yr_slopes <-  region_dat %>% 
  # mutate(recent_yr = ifelse(Year >= 2000, "2000 +", "< 2000")) %>% 
  filter(Year >= 2000) %>% 
  ggplot(aes(survey_area, l10_slope_strat, color = survey_area)) +
  ggdist::stat_halfeye(
     aes(fill = survey_area),
     adjust = .5, 
     width = .6, 
     .width = 0, 
     justification = -.2, 
     point_colour = NA) + 
  geom_boxplot(
    width = .15, 
    outlier.shape = NA) +
  geom_point(
    size = 1.3,
    alpha = .3,
    position = position_jitter(
      seed = 1, width = .1
    )
  ) +
  scale_color_gmri() +
  scale_fill_gmri() +
  scale_y_continuous(limits = c(-1.65, -0.6)) +
  facet_wrap(~survey_area, nrow = 1, scales = "free_x") +
  theme(legend.position = "none") +
  labs(x = "", y = "Size Spectrum Slope", color = "Time Period", 
       title =  "Regional Comparison: Size Spectra Slopes (2000+)")


recent_yr_slopes

allyr_slopes / early_yr_slopes / recent_yr_slopes

Gulf of Maine

spectra_surface(l10_totals = areas_aggregated$GoM, "Gulf of Maine")
plot_params("GoM")

Georges Bank

spectra_surface(l10_totals = areas_aggregated$GB, "Georges Bank")
plot_params("GB")

Southern New England

spectra_surface(l10_totals = areas_aggregated$SNE, "Southern New England")
plot_params("SNE")

Mid-Atlantic Bight

spectra_surface(l10_totals = areas_aggregated$MAB, "Mid-Atlantic Bight")
plot_params("MAB")

 

A work by Adam A. Kemberling

Akemberling@gmri.org