Spatiotemporal Modeling of Scallop Shell Height and Meat Weight

Comparing Stationary and Spatially Varying Slope Models with sdmTMB

Background

• Shell height/weight relationships allow the conversion of shell height data into biomass, which is important because shell height data are easier to collect.

• Scallop meat and gonad weight at a given shell height can vary with season and location.

• Hennen and Hart used GLMMs with data from 2001-2010.

• Research Track used a GAM.

• Yin et al. (2022) found reduced bias using spatiotemporal models

---

Available Data

• Cleaned NEFSC dredge data from 2001-2023

• Georges Bank

---

Objectives

1. Estimate scallop shell height to meat weight relationships across space and time.

2. Prepare spatial coordinate system (UTM Zone 19N) and log-transform covariates.

3. Fit baseline model assuming a constant relationship across space.

4. Construct a spatial footprint mask to avoid out of domain predictions.

5. Fit an advanced model with a Spatially Varying Coefficient (SVC) for shell height.

6. Evaluate diagnostics residuals using DHARMa.

---

Why sdmTMB

• Accounting for Spatial Autocorrelation

• Biological samples collected close together in space (e.g., scallops from neighboring survey tows within the same fishing bed) share environmental conditions like temperature, ocean currents, and food availability (phytoplankton). Standard models assume every observation is independent, which can dramatically underestimates standard errors.

• The spatial random field (\(\Omega_s\)) models the underlying spatial correlation using Gaussian random fields (via Matérn covariance functions). This accounts for patchiness in the ocean and ensures your parameter uncertainty estimates are statistically valid.

• Capturing Persistent Sweet Spots vs. Interannual Anomalies

• Persistent Spatial Field (\(\Omega_s\)): Identifies permanent geographic “sweet spots” (e.g., nutrient rich gravel beds where scallops consistently produce heavier meats relative to shell size year after year).

• Spatiotemporal Field (\(\mathcal{E}_{s,t}\)): Captures localized, time-varying anomalies (e.g., a warm water plume or localized algal bloom in 2021 that temporarily boosted meat weight in Zone A, but not Zone B).

• Localized Allometry: Growth Curves Vary Across Space

• Traditional length-weight equations assume a single global exponent (\(b\)) for the entire stock.

• By allowing \(\log(\text{SH})\) to vary smoothly across space, the model acknowledges that a 100 mm scallop on the shallow northern edge of the bank may have a vastly different meat yield than a 100 mm scallop in a deeper southern patch. This localized scaling exponent (\(b\)) is critical for spatial fisheries management.

• Correcting for Spatial Sampling Bias

• Survey effort in marine environments rarely covers every square kilometer uniformly every single year. If high density survey tows randomly cluster in high productivity areas during certain years, standard non-spatial averages will artificially skew stock biomass estimates upward.

• Model predictions are standardly evaluated over a consistent, unbiased, density weighted spatial prediction grid, effectively decoupling sampling effort bias from true population trends.

Data Prep

• Data log-transformed (MW, SH, Depth).

• Lat/Lon converted to UTM Projection (km) for accurate Euclidean spatial distance calculations.

## Setup & Data Preparation

# Load Required Libraries
library(sdmTMB)
library(dplyr)
library(sf)
library(ggplot2)
library(akima)
library(sdmTMBextra)
library(DHARMa)
library(patchwork)

# Load Data
df <- read.csv("~/Scallops/Data/gbshwv3.csv")

# Clean & Log-transform Variables
df_clean <- df %>% 
  filter(!is.na(mw), !is.na(sh), !is.na(lat), !is.na(lon), !is.na(year), !is.na(depth)) %>% 
  filter(mw > 0, sh > 0, depth > 0) %>% 
  mutate(
    log_mw = log(mw),
    log_sh = log(sh),
    log_depth = log(depth)
  )

# Project to UTM (Kilometers) - Zone 19N
df_sf <- st_as_sf(df_clean, coords = c("lon", "lat"), crs = 4326) %>% 
  st_transform(crs = 32619)

df_clean$X <- st_coordinates(df_sf)[,1] / 1000
df_clean$Y <- st_coordinates(df_sf)[,2] / 1000

Data Visualization

• Covariate Distributions & Log Transformations

# Raw vs. Log Shell Height
p1 <- ggplot(df_clean, aes(x = sh)) +
  geom_histogram(bins = 40, fill = "slategray", color = "white") +
  labs(title = "Raw Shell Height (mm)", x = "Shell Height", y = "Count") +
  theme_minimal()

p2 <- ggplot(df_clean, aes(x = log_sh)) +
  geom_histogram(bins = 40, fill = "steelblue", color = "white") +
  labs(title = "Log(Shell Height)", x = "log(SH)", y = "Count") +
  theme_minimal()

# Raw vs. Log Meat Weight
p3 <- ggplot(df_clean, aes(x = mw)) +
  geom_histogram(bins = 40, fill = "slategray", color = "white") +
  labs(title = "Raw Meat Weight (g)", x = "Meat Weight", y = "Count") +
  theme_minimal()

p4 <- ggplot(df_clean, aes(x = log_mw)) +
  geom_histogram(bins = 40, fill = "steelblue", color = "white") +
  labs(title = "Log(Meat Weight)", x = "log(MW)", y = "Count") +
  theme_minimal()

# Display 2x2 layout
(p1 + p2) / (p3 + p4)

• Sampling Intensity Over Time (Annual Tow Counts)

  ggplot(df_clean, aes(x = factor(year))) +
    geom_bar(fill = "darkseagreen", color = "white") +
    geom_text(stat = "count", aes(label = after_stat(count)), vjust = -0.5, size = 3) +
    labs(
      title = "Number of Sample Observations per Year",
      x = "Survey Year",
      y = "Sample Size (N)"
    ) +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))

  

• Depth vs. Meat Weight Relationship (Checking non-Linearity)

  ggplot(df_clean, aes(x = depth, y = log_mw)) +
    geom_point(alpha = 0.15, size = 0.8, color = "slategray") +
    geom_smooth(method = "gam", formula = y ~ s(x, k = 5), color = "firebrick", size = 1.2) +
    labs(title = "Depth vs. log(Meat Weight)",
         subtitle = "GAM smoother (red) captures potential non-linear bathymetric effects",
         x = "Depth (m)", y = "log(Meat Weight)") +
    theme_minimal()

  

• Sampling Footprint & Spatial Tow Densities

  ggplot(df_clean, aes(x = X, y = Y)) +
    stat_bin_2d(bins = 50) +
    scale_fill_viridis_c(option = "magma", name = "Sample\nDensity") +
    labs(
      title = "Spatial Sampling Intensity across Survey Area",
      subtitle = "Hex/Grid counts highlight areas of high tow concentration",
      x = "UTM Easting (km)", 
      y = "UTM Northing (km)"
    ) +
    theme_minimal()

  

• Unadjusted Shell Height vs. Meat Weight by Year

  ggplot(df_clean, aes(x = sh, y = mw)) +
    geom_point(alpha = 0.1, size = 0.5, color = "gray30") +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = FALSE, color = "royalblue") +
    facet_wrap(~year, ncol = 5) +
    labs(
      title = "Raw Shell Height vs. Meat Weight Relationships by Year",
      subtitle = "Blue line shows basic quadratic fit per year",
      x = "Shell Height (mm)",
      y = "Meat Weight (g)"
    ) +
    theme_minimal() +
    theme(strip.background = element_rect(fill = "gray95", color = NA))

  

Spatial Mesh Construction

#| warning: false

# Build Spatial Mesh
mesh <- make_mesh(df_clean, xy_cols = c("X", "Y"), cutoff = 10)

plot(mesh)
points(df_clean$X, df_clean$Y, col = "red", pch = ".")

Fit Model 1: Baseline Spatiotemporal

• Fit a basic model with fixed year effects and depth splines

• \[ \\log(mw) = \\beta_0 + \\beta_1 \\log(sh) + \\factor(year) + s(\\log(depth)) + \\omega_s + \\epsilon\_{s,t} \]

fit_depth <- sdmTMB(
  formula = log_mw ~ log_sh + as.factor(year) + s(log_depth, k = 3),
  data = df_clean,
  mesh = mesh,
  time = "year",
  spatial = "on",
  spatiotemporal = "iid",
  family = gaussian(link = "identity")
)

Masking Grid to Catch Footprint

• To prevent predicting into unsampled regions (e.g., middle of GB), we mask our interpolation grid to a 10 km buffer around positive catches.

  # Buffer around positive catches
  positive_catches <- df_clean %>% 
    filter(mw > 0) %>% 
    select(X, Y) %>% 
    distinct() %>% 
    st_as_sf(coords = c("X", "Y"))
  
  catch_buffers <- st_buffer(positive_catches, dist = 10) %>% 
    st_union()
  
  # Generate prediction grid
  grid_dense <- expand.grid(
    X = seq(min(df_clean$X), max(df_clean$X), length.out = 150),
    Y = seq(min(df_clean$Y), max(df_clean$Y), length.out = 150)
  )
  
  grid_sf <- st_as_sf(grid_dense, coords = c("X", "Y"))
  grid_masked_sf <- st_intersection(grid_sf, catch_buffers)
  
  pred_grid_spatial <- as.data.frame(st_coordinates(grid_masked_sf)) %>% 
    rename(X = X, Y = Y)
  
  pred_grid <- do.call(rbind, lapply(unique(df_clean$year), function(y) {
    pred_grid_spatial %>% mutate(year = y)
  }))
  
  # Interpolate depth
  interp_depth <- interp(x = df_clean$X, y = df_clean$Y, z = df_clean$log_depth, 
                         xo = unique(grid_dense$X), yo = unique(grid_dense$Y), duplicate = "mean")
  
  grid_depths <- expand.grid(X = interp_depth$x, Y = interp_depth$y)
  grid_depths$log_depth <- as.vector(interp_depth$z)
  
  pred_grid <- merge(pred_grid, grid_depths, by = c("X", "Y")) %>% 
    filter(!is.na(log_depth))
  
  pred_grid$log_sh <- mean(df_clean$log_sh)

Model 1 Predictions

preds <- predict(fit_depth, newdata = pred_grid)
preds_2023 <- preds %>% filter(year == 2023)

ggplot() +
  geom_raster(data = preds_2023, aes(X, Y, fill = est)) +
  scale_fill_viridis_c(name = "Predicted\nlog(Meat Wt)", na.value = "transparent") +
  geom_point(data = filter(df_clean, year == 2023), aes(X, Y), 
             color = "red", size = 0.5, alpha = 0.4) +
  labs(title = "Sampling Footprint Overlay (2023)",
       x = "UTM X (km)", y = "UTM Y (km)") +
  theme_minimal()

Look at All Years

ggplot(preds, aes(X, Y, fill = est)) +
  geom_raster() +
  scale_fill_viridis_c(name = "Predicted\nlog(Meat Wt)", option = "viridis") +
  facet_wrap(~year, ncol = 4) + 
  labs(title = "Spatio-temporal Meat Weight Index Across Years",
       x = "UTM X (km)", y = "UTM Y (km)") +
  theme_minimal() +
  theme(strip.background = element_rect(fill = "gray90", color = NA),
        strip.text = element_text(face = "bold"))

Random Fields Decomposition

• Persistent Spatial (\(\Omega_s\)) fields represent spatial variation in shell height that remains constant across time.

• It models spatial structure caused by fixed physical or biological conditions that aren’t explicitly included as covariates in your model (e.g., permanent substrate features, gravel vs. sand beds, long term bathymetry, persistent current eddies, or steady food availability zones).

  preds_fields <- predict(fit_depth, newdata = pred_grid)
  omega_plot <- preds_fields %>% filter(year == max(year))
  
  ggplot(omega_plot, aes(X, Y, fill = omega_s)) +
    geom_raster() +
    scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0, name = "Omega (ω_s)") +
    labs(title = "Persistent Spatial Random Field (Omega)",
         subtitle = "Time-invariant spatial variation",
         x = "UTM X (km)", y = "UTM Y (km)") +
    theme_minimal()

  

• Spatiotemporal Anomalies (\(\mathcal{E}_{s,t}\)) represent year to year deviations from the long-term baseline at specific locations.

  ggplot(preds_fields, aes(X, Y, fill = epsilon_st)) +
    geom_raster() +
    scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0, name = "Epsilon (ε_st)") +
    facet_wrap(~year, ncol = 4) +
    labs(title = "Spatio-Temporal Anomalies (Epsilon) by Year",
         x = "UTM X (km)", y = "UTM Y (km)") +
    theme_minimal() +
    theme(strip.background = element_rect(fill = "gray90", color = NA),
          strip.text = element_text(face = "bold"))

  

• Population Level Allometric Relationships refers to the average scaling relationship between shell height and meat weight

  sh_range <- seq(min(df_clean$sh), max(df_clean$sh), length.out = 200)
  
  overall_pred_grid <- data.frame(
    sh = sh_range,
    log_sh = log(sh_range),
    log_depth = mean(df_clean$log_depth),
    year = as.integer(median(df_clean$year))
  )
  
  overall_preds <- predict(fit_depth, newdata = overall_pred_grid, re_form = NA) %>%
    mutate(pred_mw = exp(est))
  
  ggplot() +
    geom_point(data = df_clean, aes(x = sh, y = mw), alpha = 0.15, color = "slategray") +
    geom_line(data = overall_preds, aes(x = sh, y = pred_mw), color = "firebrick", size = 1.2) +
    labs(title = "Overall Shell Height vs. Meat Weight Relationship",
         x = "Shell Height (mm)", y = "Meat Weight (g)") +
    theme_minimal()

  

• Variations by Depth

  depth_quantiles <- quantile(df_clean$depth, probs = c(0.1, 0.5, 0.9))
  
  grid_depth_vary <- expand.grid(
    sh = sh_range,
    depth = depth_quantiles,
    year = as.integer(median(df_clean$year))
  ) %>% mutate(log_sh = log(sh), log_depth = log(depth))
  
  preds_depth <- predict(fit_depth, newdata = grid_depth_vary, re_form = NA) %>%
    mutate(pred_mw = exp(est), Depth_Group = factor(paste0(round(depth), " m")))
  
  ggplot(preds_depth, aes(x = sh, y = pred_mw, color = Depth_Group)) +
    geom_line(size = 1.2) +
    scale_color_viridis_d(name = "Survey Depth") +
    labs(title = "Growth Curve Variation by Depth",
         x = "Shell Height (mm)", y = "Predicted Meat Weight (g)") +
    theme_minimal()

  

• Multi Year Shifts

  all_years <- as.integer(unique(df_clean$year))
  
  grid_all_years <- expand.grid(
    sh = sh_range,
    year = all_years,
    log_depth = mean(df_clean$log_depth)
  ) %>% mutate(log_sh = log(sh))
  
  preds_all_years <- predict(fit_depth, newdata = grid_all_years, re_form = NA) %>%
    mutate(pred_mw = exp(est))
  
  ggplot(preds_all_years, aes(x = sh, y = pred_mw, group = year, color = year)) +
    geom_line(size = 0.8, alpha = 0.85) +
    scale_color_viridis_c(name = "Survey Year", option = "plasma") +
    labs(title = "Long-Term Shifts in SH-MW Ratios",
         x = "Shell Height (mm)", y = "Predicted Meat Weight (g)") +
    theme_minimal()

  

Fit Model 2: Spatially Varying Coefficients (SVC)

• Allow the slope parameter \(\log(\text{SH})\) to vary smoothly across space

• \[ \\log(\\MW\_{s,t}) = (\\beta_1 + \\zeta_s)\\log(\\SH) + \\Factor(\\Year) + s(\\log(\\text{Depth})) + \\omega_s + \\epsilon\_{s,t} \]

  fit_svc <- sdmTMB(
    formula = log_mw ~ log_sh + as.factor(year) + s(log_depth, k = 3),
    data = df_clean,
    mesh = mesh,
    time = "year",
    spatial = "on",
    spatiotemporal = "iid",
    spatial_varying = ~ 0 + log_sh, 
    family = gaussian(link = "identity")
  )

• Localized Allometric Scaling Exponent (\(b\))

  preds_svc <- predict(fit_svc, newdata = pred_grid)
  
  global_slope <- tidy(fit_svc, effects = "fixed") %>% 
    filter(term == "log_sh") %>% 
    pull(estimate)
  
  preds_svc <- preds_svc %>% 
    mutate(local_slope = global_slope + zeta_s_log_sh)
  
  preds_slope_map <- preds_svc %>% filter(year == max(year))
  
  ggplot(preds_slope_map, aes(X, Y, fill = local_slope)) +
    geom_raster() +
    scale_fill_viridis_c(name = "Local Allometric\nExponent (b)", option = "magma") +
    labs(title = "Spatial Variation in Allometric Exponent (b)",
         x = "UTM X (km)", y = "UTM Y (km)") +
    theme_minimal()

  

• Depth smooth effect

  # Visually inspect the fitted GAM smoother for depth
  plot_smooth(fit_svc, select = 1)

  

• Calculate COG to evaluate whether the population or heavy yield zones have shifted over tim

  preds_cog <- predict(fit_svc, newdata = pred_grid, return_tmb_object = TRUE)
  cog <- get_cog(preds_cog, format = "wide")

  Bias correction is turned off.
  It is recommended to turn this on for final inference.
  Bias correction is turned off.
  It is recommended to turn this on for final inference.

  # Plot Spatial Shift of Population over Time
  ggplot(cog, aes(x = est_x, y = est_y, color = year)) +
    geom_path(arrow = arrow(length = unit(0.2, "cm")), color = "black", alpha = 0.5) +
    geom_point(size = 3) +
    scale_color_viridis_c(option = "plasma") +
    labs(
      title = "Center of Gravity (CoG) Shift Across Years",
      x = "UTM Easting (km)", y = "UTM Northing (km)"
    ) +
    theme_minimal()

  

Model Diagnostics

• Spatial Residuals

  df_clean$residuals <- residuals(fit_svc, type = "response")
  
  diagnostics_df <- df_clean %>% filter(is.finite(residuals))
  
  ggplot(diagnostics_df, aes(x = lon, y = lat, z = residuals)) +
    stat_summary_hex(fun = mean, bins = 35) + 
    scale_fill_gradient2(low = "blue", mid = "white", high = "red", name = "Mean Residual") +
    labs(title = "Binned Spatial Residuals", x = "Longitude", y = "Latitude") +
    theme_minimal()

  

• DHARMa Residuals

  sims <- simulate(fit_svc, nsim = 300, type = "mle-mvn")
  dharma_resids <- dharma_residuals(sims, fit_svc, return_DHARMa = TRUE)
  
  plot(dharma_resids)

  

• Map Prediction Uncertainty (SE)

  # Take 200 simulation draws from the joint precision matrix
  sim_draws <- predict(
    fit_svc, 
    newdata = pred_grid, 
    nsim = 200,          # Takes 200 simulation draws
    type = "response"    # Use "response" or "link" depending on your scale
  )
  
  # Returns a matrix:
  # Rows = Rows of newdata
  # Columns = Individual simulation draws (1 to 200)
  
  # Calculate point estimates (mean) and standard deviations (SE) across rows
  pred_grid$est_mean <- apply(sim_draws, 1, mean)
  pred_grid$est_se   <- apply(sim_draws, 1, sd)
  
  # Calculate 95% Confidence Intervals directly from percentiles
  pred_grid$ci_lower <- apply(sim_draws, 1, quantile, probs = 0.025)
  pred_grid$ci_upper <- apply(sim_draws, 1, quantile, probs = 0.975)
  
  ggplot(pred_grid %>% filter(year == max(year)), aes(X, Y, fill = est_se)) +
    geom_raster() +
    scale_fill_viridis_c(option = "magma", name = "SE (nsim)") +
    labs(
      title = "Spatial Prediction Standard Error (from Joint Precision Matrix)",
      subtitle = "Fast approximation using 200 simulations from the precision matrix",
      x = "UTM X (km)", 
      y = "UTM Y (km)"
    ) +
    theme_minimal()

  

• Visualize Model Hyperparameters (Range & Spatial SD)

• Showing estimated variance components (\(\sigma_O\), \(\sigma_E\)) and the Matérn spatial range (the distance in kilometers at which spatial points become independent) helps explore the spatial scale.

  ran_pars <- tidy(fit_svc, effects = "ran_pars", conf.int = TRUE)
  
  # Format with knitr::kable for slides
  knitr::kable(
    ran_pars %>% select(term, estimate, conf.low, conf.high), 
    digits = 3,
    caption = "Estimated Spatial Hyperparameters & Range Parameter"
  )

  Estimated Spatial Hyperparameters & Range Parameter

  term	estimate	conf.low	conf.high
  range	43.383	37.395	50.330
  phi	0.180	0.178	0.181
  sigma_O	0.112	0.083	0.150
  sigma_E	0.141	0.132	0.151
  sigma_Z	0.027	0.022	0.034
  sd__s(log_depth)	0.173	0.029	1.038

• Range: The distance at which two sampling locations become spatially independent (km).

• Phi: The standard deviation of the Guassian residual error. It represents unexplained individula level measurment noise.

• Sigma_O: Measures the magnitude of he variation in the persistent spatial field. A high value mean there are strong “sweet” and “cold” spots.

• Sigma_E: Spatiotemporal SD, measures the magnitude of variation in year to year spatial anomalies. A high value means the environmental anomalies fluctuate significantly in space from year to year

• Sigma_Z: Spatiotemporal SD, measures how much the slope varies across space. A high value represents that it scales substantially with space.

• sd_s(log_depth): Represents how strong the linear relationship is with depth. A higher values represents a non-linear relationship.

• Restricting model projections to populated survey buffers avoids misrepresentation outside of sampled domains. This can be further refined.

• Depths and annual shifts play a significant role in determining scallop meat yield.

• Localized model fits demonstrate that scallops scale differently in growth dependent on localized spatial habitats.

Next steps

• Get the full data set (GB, MA and current data)

• Compare models to current methods

• Explore ecovs (e.g., bottom temp, phytoplankton, bottom type)

• Produce SH/MW ratios for each SAMS area?