Shell Height - Meat Weight (SH-MW) Spatially Varying Coefficient Model for MAB

Spatial Allometry with sdmTMB

1 Overview

This analysis fits a Spatially Varying Coefficient (SVC) models to analyze the spatial, temporal, and environmental drivers of the shell height to meat weight relationship in sea scallops. We evaluate candidate models using AIC, perform residual diagnostics, conduct cross-validation, and visualize marginal effects and regional parameters.

2 Data Setup and Processing

In this section, we load required libraries, import the data, calculate station-level scallop density proxies, transform predictors, and project spatial coordinates to UTM Zone 19N (km) for spatial mesh construction.

# Load required libraries
library(sdmTMB)
library(dplyr)
library(ggplot2)
library(mgcv) # For s() smooth terms
library(sf)
library(viridis)
library(sdmTMBextra)
library(patchwork)
library(gstat)

# Read and prepare data
df <- read.csv("~/Scallops/Data/mashw0125.csv")

# Clean data, calculate density proxy, and scale predictors
df_clean <- df %>%
  filter(
    !is.na(mw), mw > 0,
    !is.na(sh), sh > 0,
    !is.na(Lat), !is.na(Lon),
    !is.na(date)
  ) %>%
  # Calculate station-level density (number of scallops sampled/caught per station)
  group_by(StationID) %>%
  mutate(station_density = n()) %>%
  ungroup() %>%
  mutate(
    log_sh = log(sh),                                 # Log shell height
    log_sh_c = log_sh - mean(log_sh, na.rm = TRUE),   # Centered log(sh)
    log_density = log(station_density),               # Log-transformed density proxy
    log_density_c = scale(log_density, center = TRUE, scale = TRUE)[, 1],
    year_f = as.factor(year),
    date_obj = as.Date(date)
  )

# Convert Lat/Lon (WGS84 EPSG:4326) to UTM Zone 19N (EPSG:32619) in kilometers
df_sf <- st_as_sf(df_clean, coords = c("Lon", "Lat"), crs = 4326) %>%
  st_transform(crs = 32619)

# Extract projected coordinates and scale to kilometers (km)
coords <- st_coordinates(df_sf) / 1000
df_clean$X <- coords[, 1]
df_clean$Y <- coords[, 2]

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

3 Data Explorations

Tows per year

tows_per_year <- df_clean %>%
  group_by(year) %>%
  summarize(
    num_tows = n_distinct(StationID),  # Count of unique station/tow IDs
    total_scallops = n(),              # Total individual scallops measured
    .groups = "drop"
  )

# 2. Plot number of sample tows per year
p_tows <- ggplot(tows_per_year, aes(x = factor(year), y = num_tows)) +
  geom_col(fill = "#2b5c8f", color = "white", alpha = 0.85, width = 0.7) +
  geom_text(aes(label = num_tows), vjust = -0.5, size = 3.5, color = "black") +
  scale_y_continuous(expand = expansion(mult = c(0, 0.12))) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
    panel.grid.major.x = element_blank()
  ) +
  labs(
    title = "Survey Effort: Number of Sample Tows Per Year",
    subtitle = "Count of unique stations sampled in each survey year",
    x = "Survey Year",
    y = "Number of Tows (Stations)"
  )

print(p_tows)

Spatial and seasonal coverage

# Station level aggregation for mapping
station_map_data <- df_clean %>%
  group_by(StationID, Lat, Lon) %>%
  summarize(
    mean_mw = mean(mw, na.rm = TRUE),
    mean_sh = mean(sh, na.rm = TRUE),
    mean_depth = mean(Depth, na.rm = TRUE),
    station_density = first(station_density),
    .groups = "drop"
  )

# Spatial Density Map
ggplot(station_map_data, aes(x = Lon, y = Lat)) +
  geom_point(aes(color = log(station_density), size = log(mean_mw)), alpha = 0.7) +
  scale_color_viridis_c(option = "viridis", name = "Log Station\nDensity") +
  coord_quickmap() +
  theme_minimal() +
  labs(
    title = "Spatial Survey Map: Station Density & Mean Meat Weight",
    x = "Longitude", y = "Latitude", size = "Mean MW (g)"
  )

# Seasonal Sampling Windows
ggplot(df_clean, aes(x = jday, y = year_f)) +
  geom_point(alpha = 0.2, color = "#2e7d32", size = 1) +
  theme_minimal() +
  labs(
    title = "Sampling Seasonality Window Across Survey Years",
    subtitle = "Checks for temporal shifts in survey timing (jday)",
    x = "Julian Day of Year (jday)", y = "Survey Year"
  )

Distributions

p1 <- ggplot(df_clean, aes(x = mw)) +
  geom_histogram(bins = 50, fill = "#2b5c8f", color = "white", alpha = 0.85) +
  theme_minimal() +
  labs(title = "Meat Weight (mw)", x = "Meat Weight (g)", y = "Count")

p2 <- ggplot(df_clean, aes(x = sh)) +
  geom_histogram(bins = 50, fill = "#2e7d32", color = "white", alpha = 0.85) +
  theme_minimal() +
  labs(title = "Shell Height (sh)", x = "Shell Height (mm)", y = "Count")

p3 <- ggplot(df_clean, aes(x = station_density)) +
  geom_histogram(bins = 50, fill = "#d95f02", color = "white", alpha = 0.85) +
  scale_x_log10() +
  theme_minimal() +
  labs(title = "Station Density (Sample Count/Station)", x = "Scallops per Station (log scale)", y = "Count")

p4 <- ggplot(df_clean, aes(x = Depth)) +
  geom_histogram(bins = 50, fill = "#7570b3", color = "white", alpha = 0.85) +
  theme_minimal() +
  labs(title = "Depth Distribution", x = "Depth (m)", y = "Count")

# Combine histograms into a 2x2 grid
(p1 + p2) / (p3 + p4)

Key Relationships

# Allometric Relationship
ggplot(df_clean, aes(x = log_sh_c, y = log(mw))) +
  geom_point(alpha = 0.1, size = 0.8, color = "#2b5c8f") +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 5), color = "firebrick", linewidth = 1.2) +
  theme_minimal() +
  labs(
    title = "Allometric Relationship: Centered Log Shell Height vs. Log Meat Weight",
    subtitle = "GAM smooth overlay testing for non-linear allometric shifts",
    x = "Centered Log Shell Height (log_sh_c)",
    y = "Log Meat Weight log(mw)"
  )

# Density Dependence
ggplot(df_clean, aes(x = log_density_c, y = log(mw))) +
  geom_point(alpha = 0.1, color = "darkslategrey", size = 0.8) +
  geom_smooth(method = "gam", color = "darkorange", linewidth = 1.2) +
  theme_minimal() +
  labs(
    title = "Effect of Scaled Station Density on Scallop Meat Weight",
    subtitle = "Testing for density-dependent growth or competition effects",
    x = "Scaled Log Station Density (log_density_c)",
    y = "Log Meat Weight log(mw)"
  )

# Temporal Variations
ggplot(df_clean, aes(x = year_f, y = mw)) +
  geom_boxplot(fill = "aliceblue", color = "#2b5c8f", outlier.alpha = 0.05, outlier.size = 0.5) +
  scale_y_log10() +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(
    title = "Meat Weight Variations Across Survey Years",
    x = "Year",
    y = "Meat Weight (g) [log scale]"
  )

# Raw vs Transformed Scale Side-by-Side
p_raw <- ggplot(df_clean, aes(x = sh, y = mw)) +
  geom_point(alpha = 0.1, size = 0.6, color = "#2b5c8f") +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 5), color = "firebrick") +
  theme_minimal() +
  labs(title = "A) Raw Scale (Non-Linear)", x = "Shell Height (mm)", y = "Meat Weight (g)")

p_log <- ggplot(df_clean, aes(x = log_sh_c, y = log(mw))) +
  geom_point(alpha = 0.1, size = 0.6, color = "#2e7d32") +
  geom_smooth(method = "lm", color = "firebrick") +
  theme_minimal() +
  labs(title = "B) Log-Transformed & Centered (Linearized)", x = "Centered Log Shell Height (log_sh_c)", y = "Log Meat Weight log(mw)")

p_raw + p_log

Spatial and Seasonal Coverage

# Station level aggregation for mapping
station_map_data <- df_clean %>%
  group_by(StationID, Lat, Lon) %>%
  summarize(
    mean_mw = mean(mw, na.rm = TRUE),
    mean_sh = mean(sh, na.rm = TRUE),
    mean_depth = mean(Depth, na.rm = TRUE),
    station_density = first(station_density),
    .groups = "drop"
  )

# Spatial Density Map
ggplot(station_map_data, aes(x = Lon, y = Lat)) +
  geom_point(aes(color = log(station_density), size = log(mean_mw)), alpha = 0.7) +
  scale_color_viridis_c(option = "viridis", name = "Log Station\nDensity") +
  coord_quickmap() +
  theme_minimal() +
  labs(
    title = "Spatial Survey Map: Station Density & Mean Meat Weight",
    x = "Longitude", y = "Latitude", size = "Mean MW (g)"
  )

# Seasonal Sampling Windows
ggplot(df_clean, aes(x = jday, y = year_f)) +
  geom_point(alpha = 0.2, color = "#2e7d32", size = 1) +
  theme_minimal() +
  labs(
    title = "Sampling Seasonality Window Across Survey Years",
    subtitle = "Checks for temporal shifts in survey timing (jday)",
    x = "Julian Day of Year (jday)", y = "Survey Year"
  )

4 Model Fitting and Selection

# Model 1: Base Allometric Model
fit_m1_base <- sdmTMB(
  formula = mw ~ log_sh_c,
  data = df_clean,
  mesh = mesh,
  family = Gamma(link = "log"),
  spatial = "on",
  time = "year",
  spatiotemporal = "iid"
)

# Model 2: Base + Depth
fit_m2_depth <- sdmTMB(
  formula = mw ~ log_sh_c + s(Depth),
  data = df_clean,
  mesh = mesh,
  family = Gamma(link = "log"),
  spatial = "on",
  time = "year",
  spatiotemporal = "iid"
)

# Model 3: Base + Depth + Seasonal Smooth (jday)
fit_m3_jday <- sdmTMB(
  formula = mw ~ log_sh_c + s(Depth) + s(jday),
  data = df_clean,
  mesh = mesh,
  family = Gamma(link = "log"),
  spatial = "on",
  time = "year",
  spatiotemporal = "iid"
)

# Model 4: Base + Depth + Seasonal Smooth (jday) + Density
fit_m4_density <- sdmTMB(
  formula = mw ~ log_sh_c + s(Depth) + s(jday) + log_density_c,
  data = df_clean,
  mesh = mesh,
  family = Gamma(link = "log"),
  spatial = "on",
  time = "year",
  spatiotemporal = "iid"
)

# Model 5: Full Model with Spatially Varying Coefficient (SVC for Shell Height)
fit_m5_svc <- sdmTMB(
  formula = mw ~ log_sh_c + s(Depth) + s(jday) + log_density_c,
  data = df_clean,
  mesh = mesh,
  family = Gamma(link = "log"),
  spatial = "on",
  spatial_varying = ~ log_sh_c,  # Spatial allometric slope
  time = "year",
  spatiotemporal = "iid"
)

AIC Table

model_list <- list(
  "M1: Base Allometry"               = fit_m1_base,
  "M2: + Depth"                      = fit_m2_depth,
  "M3: + Depth Julian Day"           = fit_m3_jday,
  "M4: + Depth Julian Day & Density" = fit_m4_density,
  "M5: Full Model + SVC"             = fit_m5_svc
)

aic_results <- data.frame(
  Model  = names(model_list),
  logLik = sapply(model_list, logLik),
  df     = sapply(model_list, function(m) attr(logLik(m), "df")),
  AIC    = sapply(model_list, AIC)
) %>%
  mutate(
    deltaAIC = AIC - min(AIC),
    weight   = exp(-0.5 * deltaAIC) / sum(exp(-0.5 * deltaAIC))
  ) %>%
  arrange(deltaAIC)

knitr::kable(aic_results, digits = 3, caption = "Model Selection via AIC Table")

Model Selection via AIC Table

	Model	logLik	df	AIC	deltaAIC	weight
M5: Full Model + SVC	M5: Full Model + SVC	-170635.6	12	341295.2	0.000	1
M4: + Depth Julian Day & Density	M4: + Depth Julian Day & Density	-171237.2	11	342496.5	1201.312	0
M3: + Depth Julian Day	M3: + Depth Julian Day	-171252.5	10	342525.0	1229.801	0
M2: + Depth	M2: + Depth	-171261.7	8	342539.3	1244.172	0
M1: Base Allometry	M1: Base Allometry	-171364.9	6	342741.9	1446.706	0

5 Model Diagnostics

QQ-plot

df_clean$resids <- residuals(fit_m5_svc, type = "mle-mvn")

# Filter finite values
df_resids <- df_clean %>% filter(is.finite(resids))

# Set degrees of freedom for the t-distribution
df_param <- 5  # Replace with your desired degrees of freedom

# 1. Generate theoretical quantiles from the t-distribution
n <- length(df_resids$resids)
theoretical_quantiles <- qt(ppoints(n), df = df_param)

# 2. Plot sample residuals against theoretical t-quantiles
qqplot(
  theoretical_quantiles, df_resids$resids,
  xlab = paste0("Theoretical Quantiles (t-dist, df = ", df_param, ")"),
  ylab = "Sample Quantiles",
  main = "t-Distribution Q-Q Plot"
)

# 3. Add a reference line passing through the 1st and 3rd quartiles
qqline(
  df_resids$resids, 
  distribution = function(p) qt(p, df = df_param), 
  col = "red", 
  lwd = 2
)

Spatial Residuals

# Annual station-level residual aggregation
spatial_resids_by_year <- df_resids %>%
  group_by(year, StationID, X, Y) %>%
  summarize(
    mean_resid = mean(resids, na.rm = TRUE),
    n_obs = n(),
    .groups = "drop"
  )

# Plot Spatial Residuals Faceted by Year
ggplot(spatial_resids_by_year, aes(x = X, y = Y, color = mean_resid)) +
  geom_point(aes(size = n_obs), alpha = 0.8) +
  scale_color_gradient2(
    low = "#0571b0", mid = "white", high = "#ca0020",
    midpoint = 0, name = "Mean Residual"
  ) +
  scale_size_continuous(range = c(1, 3.5), guide = "none") +
  facet_wrap(~ year, ncol = 6) +
  coord_fixed() +
  theme_bw() +
  theme(panel.background = element_rect(fill = "white"), panel.grid = element_blank()) +
  labs(
    title = "Spatial Residuals by Year",
    subtitle = "Blue = Underpredicting | White = Zero | Red = Overpredicting",
    x = "UTM X (km)", y = "UTM Y (km)"
  )

6 Try to resolve residual patter

S-tail pattern in the QQ-plot which indicates that there are more extreme values than the assumed distribution. Let’s try some different model configurations to see if we can remove the pattern.

Try a finer mesh resolution

mesh_fine <- make_mesh(df_clean, xy_cols = c("X", "Y"), cutoff = 5)


fit_m6_svc <- sdmTMB(
  formula = mw ~ log_sh_c + s(Depth) + s(jday) + log_density_c,
  data = df_clean,
  mesh = mesh_fine,
  family = Gamma(link = "log"),
  spatial = "on",
  spatial_varying = ~ log_sh_c,  # Spatial allometric slope
  time = "year",
  spatiotemporal = "iid"
)

Compare residuals of finer mesh to old method

# Define degrees of freedom for the theoretical t-distribution
t_df <- 5 # Adjust this value as appropriate for your model

# ------------------------------------------------------------------------------
# Extract Randomized Quantile Residuals & Predictions
# ------------------------------------------------------------------------------
df_clean$resids_m5 <- residuals(fit_m5_svc, type = "mle-mvn")
df_clean$resids_m6 <- residuals(fit_m6_svc, type = "mle-mvn")

df_clean$fitted_m5 <- predict(fit_m5_svc)$est
df_clean$fitted_m6 <- predict(fit_m6_svc)$est

# Filter out non-finite values
df_resids <- df_clean %>% 
  filter(is.finite(resids_m5), is.finite(resids_m6))

# ------------------------------------------------------------------------------
# Q-Q Plot Comparison (t-Distribution)
# ------------------------------------------------------------------------------
qq_m5 <- ggplot(df_resids, aes(sample = resids_m5)) +
  stat_qq(distribution = qt, dparams = list(df = t_df), alpha = 0.2, color = "#2b5c8f", size = 0.7) +
  stat_qq_line(distribution = qt, dparams = list(df = t_df), color = "firebrick", linewidth = 0.8) +
  theme_minimal() +
  labs(title = "M1: Q-Q Plot (t-dist)", x = paste0("Theoretical Quantiles (t, df = ", t_df, ")"), y = "Sample Quantiles")

qq_m6 <- ggplot(df_resids, aes(sample = resids_m6)) +
  stat_qq(distribution = qt, dparams = list(df = t_df), alpha = 0.2, color = "#2e7d32", size = 0.7) +
  stat_qq_line(distribution = qt, dparams = list(df = t_df), color = "firebrick", linewidth = 0.8) +
  theme_minimal() +
  labs(title = "M5: Q-Q Plot (t-dist)", x = paste0("Theoretical Quantiles (t, df = ", t_df, ")"), y = "Sample Quantiles")

# ------------------------------------------------------------------------------
# Residuals vs. Fitted Values
# ------------------------------------------------------------------------------
fit_m5 <- ggplot(df_resids, aes(x = fitted_m5, y = resids_m5)) +
  geom_point(alpha = 0.15, size = 0.7, color = "#2b5c8f") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_smooth(method = "gam", color = "darkorange", se = TRUE) +
  theme_minimal() +
  labs(title = "M1: Residuals vs. Fitted", x = "Predicted Log MW", y = "Quantile Residuals")

fit_m6 <- ggplot(df_resids, aes(x = fitted_m6, y = resids_m6)) +
  geom_point(alpha = 0.15, size = 0.7, color = "#2e7d32") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_smooth(method = "gam", color = "darkorange", se = TRUE) +
  theme_minimal() +
  labs(title = "M5: Residuals vs. Fitted", x = "Predicted Log MW", y = "Quantile Residuals")

# ------------------------------------------------------------------------------
# Residuals vs. Depth (Covariate Check Example)
# ------------------------------------------------------------------------------
cov_m5 <- ggplot(df_resids, aes(x = Depth, y = resids_m5)) +
  geom_point(alpha = 0.15, size = 0.7, color = "#2b5c8f") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_smooth(method = "gam", color = "darkorange", se = TRUE) +
  theme_minimal() +
  labs(title = "M1: Residuals vs. Depth", x = "Depth (m)", y = "Quantile Residuals")

cov_m6 <- ggplot(df_resids, aes(x = Depth, y = resids_m6)) +
  geom_point(alpha = 0.15, size = 0.7, color = "#2e7d32") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_smooth(method = "gam", color = "darkorange", se = TRUE) +
  theme_minimal() +
  labs(title = "M5: Residuals vs. Depth", x = "Depth (m)", y = "Quantile Residuals")

# ------------------------------------------------------------------------------
# Spatial Empirical Variograms
# ------------------------------------------------------------------------------
v_m5 <- variogram(resids_m5 ~ 1, locations = ~ X + Y, data = df_resids, cutoff = 100)
v_m5$Model <- "M5: Full SVC Model"

v_m6 <- variogram(resids_m6 ~ 1, locations = ~ X + Y, data = df_resids, cutoff = 100)
v_m6$Model <- "M6: Full SVC Model Fine Mesh"

vario_df <- rbind(v_m5, v_m6)

p_vario <- ggplot(vario_df, aes(x = dist, y = gamma, color = Model)) +
  geom_point(size = 2) +
  geom_line(linewidth = 1) +
  scale_color_manual(values = c("M5: Full SVC Model" = "#2b5c8f", "M6: Full SVC Model Fine Mesh" = "#2e7d32")) +
  theme_minimal() +
  labs(
    title = "Empirical Residual Variogram Comparison",
    subtitle = "A flat line indicates spatial autocorrelation has been fully accounted for",
    x = "Distance (km)",
    y = "Semi-variance"
  )

# ------------------------------------------------------------------------------
# Display Dashboard Layout
# ------------------------------------------------------------------------------
# Arrange 3x2 Grid for Diagnostic Checks
(qq_m5 + qq_m6) / (fit_m5 + fit_m6) / (cov_m5 + cov_m6)

# Display Spatial Variogram
print(p_vario)

# ------------------------------------------------------------------------------
# Standalone Q-Q Comparison (t-Distribution)
# ------------------------------------------------------------------------------
df_qq_comp <- df_clean %>%
  select(resids_m5, resids_m6) %>%
  filter(is.finite(resids_m5), is.finite(resids_m6))

p_qq_m5 <- ggplot(df_qq_comp, aes(sample = resids_m5)) +
  stat_qq(distribution = qt, dparams = list(df = t_df), alpha = 0.25, color = "#2b5c8f", size = 0.8) +
  stat_qq_line(distribution = qt, dparams = list(df = t_df), color = "firebrick", linewidth = 1) +
  theme_minimal() +
  labs(
    title = "M1: Base Model",
    subtitle = "Formula: mw ~ log_sh_c",
    x = paste0("Theoretical Quantiles (t, df = ", t_df, ")"),
    y = "Sample Quantiles"
  )

p_qq_m6 <- ggplot(df_qq_comp, aes(sample = resids_m6)) +
  stat_qq(distribution = qt, dparams = list(df = t_df), alpha = 0.25, color = "#2e7d32", size = 0.8) +
  stat_qq_line(distribution = qt, dparams = list(df = t_df), color = "firebrick", linewidth = 1) +
  theme_minimal() +
  labs(
    title = "M5: Full SVC Model",
    subtitle = "Formula: mw ~ log_sh_c + s(Depth) + s(jday) + log_density_c",
    x = paste0("Theoretical Quantiles (t, df = ", t_df, ")"),
    y = "Sample Quantiles"
  )

# Combine plots side-by-side
p_qq_m5 + p_qq_m6

Try a different correlation structure

fit_m7_svc <- sdmTMB(
  formula = mw ~ log_sh_c + s(Depth) + s(jday) + log_density_c,
  data = df_clean,
  mesh = mesh,
  family = Gamma(link = "log"),
  spatial = "on",
  spatial_varying = ~ log_sh_c,  # Spatial allometric slope
  time = "year",
  spatiotemporal = "ar1",
  extra_time = c(2003)
)

fit_m8_svc <- sdmTMB(
  formula = mw ~ log_sh_c + s(Depth) + s(jday) + log_density_c,
  data = df_clean,
  mesh = mesh,
  family = Gamma(link = "log"),
  spatial = "on",
  spatial_varying = ~ log_sh_c,  # Spatial allometric slope
  time = "year",
  spatiotemporal = "rw",
  extra_time = c(2003)
)

# Define degrees of freedom for the theoretical t-distribution
t_df <- 5 # Adjust this value as appropriate for your model

# ------------------------------------------------------------------------------
# Calculate Randomized Quantile Residuals & Predictions
# ------------------------------------------------------------------------------
# Extract "mle-mvn" randomized quantile residuals for GLMMs
df_clean$resids_m5 <- residuals(fit_m5_svc, type = "mle-mvn")
df_clean$resids_m7 <- residuals(fit_m7_svc, type = "mle-mvn")
df_clean$resids_m8 <- residuals(fit_m8_svc, type = "mle-mvn")

df_clean$fitted_m5 <- predict(fit_m5_svc)$est
df_clean$fitted_m7 <- predict(fit_m7_svc)$est
df_clean$fitted_m8 <- predict(fit_m8_svc)$est

# Filter out non-finite values across all models
df_resids <- df_clean %>% 
  filter(
    is.finite(resids_m5), 
    is.finite(resids_m7), 
    is.finite(resids_m8)
  )

# Define a consistent color palette for the 3 models
model_colors <- c(
  "M5: SVC Base"    = "#2b5c8f", 
  "M7: SVC Model 7" = "#2e7d32", 
  "M8: SVC Model 8" = "#d95f02"
)

# ------------------------------------------------------------------------------
# Compute AIC Model Comparison Table
# ------------------------------------------------------------------------------
aic_table <- data.frame(
  Model = c("M5: SVC Base", "M7: SVC Model 7", "M8: SVC Model 8"),
  df    = c(attr(logLik(fit_m5_svc), "df"), 
            attr(logLik(fit_m7_svc), "df"), 
            attr(logLik(fit_m8_svc), "df")),
  AIC   = c(AIC(fit_m5_svc), AIC(fit_m7_svc), AIC(fit_m8_svc))
) %>%
  mutate(
    delta_AIC = AIC - min(AIC),
    weight    = exp(-0.5 * delta_AIC) / sum(exp(-0.5 * delta_AIC))
  ) %>%
  arrange(AIC)

# ------------------------------------------------------------------------------
# Q-Q Plot Comparisons (Student's t-Distribution)
# ------------------------------------------------------------------------------
qq_m5 <- ggplot(df_resids, aes(sample = resids_m5)) +
  stat_qq(distribution = qt, dparams = list(df = t_df), alpha = 0.2, color = model_colors["M5: SVC Base"], size = 0.7) +
  stat_qq_line(distribution = qt, dparams = list(df = t_df), color = "firebrick", linewidth = 0.8) +
  theme_minimal() +
  labs(title = "M5: Q-Q Plot (t-dist)", x = paste0("Theoretical Quantiles (t, df = ", t_df, ")"), y = "Sample Quantiles")

qq_m7 <- ggplot(df_resids, aes(sample = resids_m7)) +
  stat_qq(distribution = qt, dparams = list(df = t_df), alpha = 0.2, color = model_colors["M7: SVC Model 7"], size = 0.7) +
  stat_qq_line(distribution = qt, dparams = list(df = t_df), color = "firebrick", linewidth = 0.8) +
  theme_minimal() +
  labs(title = "M7: Q-Q Plot (t-dist)", x = paste0("Theoretical Quantiles (t, df = ", t_df, ")"), y = "Sample Quantiles")

qq_m8 <- ggplot(df_resids, aes(sample = resids_m8)) +
  stat_qq(distribution = qt, dparams = list(df = t_df), alpha = 0.2, color = model_colors["M8: SVC Model 8"], size = 0.7) +
  stat_qq_line(distribution = qt, dparams = list(df = t_df), color = "firebrick", linewidth = 0.8) +
  theme_minimal() +
  labs(title = "M8: Q-Q Plot (t-dist)", x = paste0("Theoretical Quantiles (t, df = ", t_df, ")"), y = "Sample Quantiles")

# ------------------------------------------------------------------------------
# Residuals vs. Fitted Values
# ------------------------------------------------------------------------------
fit_m5 <- ggplot(df_resids, aes(x = fitted_m5, y = resids_m5)) +
  geom_point(alpha = 0.15, size = 0.7, color = model_colors["M5: SVC Base"]) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_smooth(method = "gam", color = "black", se = TRUE, linewidth = 0.8) +
  theme_minimal() +
  labs(title = "M5: Resids vs Fitted", x = "Predicted Log MW", y = "Quantile Residuals")

fit_m7 <- ggplot(df_resids, aes(x = fitted_m7, y = resids_m7)) +
  geom_point(alpha = 0.15, size = 0.7, color = model_colors["M7: SVC Model 7"]) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_smooth(method = "gam", color = "black", se = TRUE, linewidth = 0.8) +
  theme_minimal() +
  labs(title = "M7: Resids vs Fitted", x = "Predicted Log MW", y = "Quantile Residuals")

fit_m8 <- ggplot(df_resids, aes(x = fitted_m8, y = resids_m8)) +
  geom_point(alpha = 0.15, size = 0.7, color = model_colors["M8: SVC Model 8"]) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_smooth(method = "gam", color = "black", se = TRUE, linewidth = 0.8) +
  theme_minimal() +
  labs(title = "M8: Resids vs Fitted", x = "Predicted Log MW", y = "Quantile Residuals")

# ------------------------------------------------------------------------------
# Residuals vs. Depth (Covariate Check Example)
# ------------------------------------------------------------------------------
cov_m5 <- ggplot(df_resids, aes(x = Depth, y = resids_m5)) +
  geom_point(alpha = 0.15, size = 0.7, color = model_colors["M5: SVC Base"]) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_smooth(method = "gam", color = "black", se = TRUE, linewidth = 0.8) +
  theme_minimal() +
  labs(title = "M5: Resids vs Depth", x = "Depth (m)", y = "Quantile Residuals")

cov_m7 <- ggplot(df_resids, aes(x = Depth, y = resids_m7)) +
  geom_point(alpha = 0.15, size = 0.7, color = model_colors["M7: SVC Model 7"]) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_smooth(method = "gam", color = "black", se = TRUE, linewidth = 0.8) +
  theme_minimal() +
  labs(title = "M7: Resids vs Depth", x = "Depth (m)", y = "Quantile Residuals")

cov_m8 <- ggplot(df_resids, aes(x = Depth, y = resids_m8)) +
  geom_point(alpha = 0.15, size = 0.7, color = model_colors["M8: SVC Model 8"]) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_smooth(method = "gam", color = "black", se = TRUE, linewidth = 0.8) +
  theme_minimal() +
  labs(title = "M8: Resids vs Depth", x = "Depth (m)", y = "Quantile Residuals")

# ------------------------------------------------------------------------------
# Spatial Empirical Variograms
# ------------------------------------------------------------------------------
v_m5 <- variogram(resids_m5 ~ 1, locations = ~ X + Y, data = df_resids, cutoff = 100)
v_m5$Model <- "M5: SVC Base"

v_m7 <- variogram(resids_m7 ~ 1, locations = ~ X + Y, data = df_resids, cutoff = 100)
v_m7$Model <- "M7: SVC Model 7"

v_m8 <- variogram(resids_m8 ~ 1, locations = ~ X + Y, data = df_resids, cutoff = 100)
v_m8$Model <- "M8: SVC Model 8"

vario_df <- rbind(v_m5, v_m7, v_m8)

p_vario <- ggplot(vario_df, aes(x = dist, y = gamma, color = Model)) +
  geom_point(size = 2) +
  geom_line(linewidth = 1) +
  scale_color_manual(values = model_colors) +
  theme_minimal() +
  labs(
    title = "Empirical Residual Variogram Comparison across M5, M7, M8",
    subtitle = "Flatter curves indicate better removal of spatial autocorrelation",
    x = "Distance (km)",
    y = "Semi-variance"
  )

# ------------------------------------------------------------------------------
# Display Model Comparison & Diagnostic Plots
# ------------------------------------------------------------------------------
# Print AIC Table
cat("\n=== AIC Model Comparison Table ===\n")

=== AIC Model Comparison Table ===

print(aic_table)

            Model df      AIC delta_AIC        weight
1 M7: SVC Model 7 13 341244.9   0.00000  1.000000e+00
2    M5: SVC Base 12 341295.2  50.22563  1.240632e-11
3 M8: SVC Model 8 12 341993.4 748.41764 3.042104e-163

# Display 3x3 Diagnostic Grid
dashboard_grid <- (qq_m5 | qq_m7 | qq_m8) / 
                  (fit_m5 | fit_m7 | fit_m8) / 
                  (cov_m5 | cov_m7 | cov_m8)
dashboard_grid

# Print Spatial Variogram
print(p_vario)

More residuals

# Calculate Randomized Quantile Residuals
# Note: Using "mle-mvn" provides true randomized quantile residuals for GLMMs
df_clean$resids <- residuals(fit_m7_svc, type = "mle-mvn")

# Filter out non-finite residual values
df_resids <- df_clean %>% filter(is.finite(resids))

# ------------------------------------------------------------------------------
# Residuals vs. Continuous Covariates
# ------------------------------------------------------------------------------

# A. Centered Log Shell Height (log_sh_c)
p_res_sh <- ggplot(df_resids, aes(x = log_sh_c, y = resids)) +
  geom_point(alpha = 0.15, size = 0.8, color = "#2b5c8f") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red", linewidth = 1) +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 5), color = "darkorange", se = TRUE) +
  theme_minimal() +
  labs(
    title = "Residuals vs. Centered Log Shell Height",
    x = "Centered Log Shell Height (log_sh_c)",
    y = "Quantile Residuals"
  )

# Water Depth (Depth)
p_res_depth <- ggplot(df_resids, aes(x = Depth, y = resids)) +
  geom_point(alpha = 0.15, size = 0.8, color = "#7570b3") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red", linewidth = 1) +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 5), color = "darkorange", se = TRUE) +
  theme_minimal() +
  labs(
    title = "Residuals vs. Water Depth",
    x = "Depth (m)",
    y = "Quantile Residuals"
  )

# Seasonal Timing (Julian Day - jday)
p_res_jday <- ggplot(df_resids, aes(x = jday, y = resids)) +
  geom_point(alpha = 0.15, size = 0.8, color = "#2e7d32") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red", linewidth = 1) +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 5), color = "darkorange", se = TRUE) +
  theme_minimal() +
  labs(
    title = "Residuals vs. Julian Day",
    x = "Julian Day of Year (jday)",
    y = "Quantile Residuals"
  )

# Scaled Log Station Density (log_density_c)
p_res_density <- ggplot(df_resids, aes(x = log_density_c, y = resids)) +
  geom_point(alpha = 0.15, size = 0.8, color = "#d95f02") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red", linewidth = 1) +
  geom_smooth(method = "gam", formula = y ~ s(x, k = 5), color = "darkorange", se = TRUE) +
  theme_minimal() +
  labs(
    title = "Residuals vs. Scaled Station Density",
    x = "Scaled Log Density (log_density_c)",
    y = "Quantile Residuals"
  )

# ------------------------------------------------------------------------------
# Residuals Across Factor / Discrete Covariates
# ------------------------------------------------------------------------------

# Survey Year (year_f)
p_res_year <- ggplot(df_resids, aes(x = year_f, y = resids)) +
  geom_boxplot(fill = "aliceblue", color = "#2b5c8f", outlier.alpha = 0.1, outlier.size = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red", linewidth = 1) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(
    title = "Residuals Across Survey Years",
    x = "Year",
    y = "Quantile Residuals"
  )

# ------------------------------------------------------------------------------
# Display Arranged Diagnostic Grid
# ------------------------------------------------------------------------------

# 2x2 Grid for Continuous Covariates
(p_res_sh + p_res_depth) / (p_res_jday + p_res_density)

# Year Boxplot
p_res_year

7 Cross Validation

# Runs k-fold CV across 5 random folds
#cv_random <- sdmTMB_cv(
#  formula = mw ~ log_sh_c + s(Depth) + s(jday) + log_density_c,
#  data = df_clean,
#  mesh = mesh,
#  family = Gamma(link = "log"),
#  spatial = "on",
#  spatial_varying = ~ log_sh_c,
#  time = "year",
#  spatiotemporal = "ar1",
#  k_folds = 5
#)

#cat("Random CV Sum Log-Likelihood:", cv_random$sum_loglik, "\n")

#rmse_random <- sqrt(mean((cv_random$data$mw - exp(cv_random$data$cv_predicted))^2))
#cat("Random CV Out-of-Sample RMSE:", round(rmse_random, 4), "\n")

8 Marginal Effects

# 1. Depth Effect
depth_seq <- seq(min(df_clean$Depth), max(df_clean$Depth), length.out = 100)
sh_targets <- c(80, 100, 120, 140)

pred_grid_depth <- expand.grid(
  Depth = depth_seq, sh = sh_targets,
  jday = mean(df_clean$jday, na.rm = TRUE),
  log_density_c = 0, year = min(df_clean$year)
) %>%
  mutate(
    log_sh = log(sh),
    log_sh_c = log_sh - mean(log(df_clean$sh), na.rm = TRUE),
    sh_factor = paste0("SH = ", sh, " mm")
  )

preds_depth <- predict(fit_m7_svc, newdata = pred_grid_depth, re_form = NA, se_fit = TRUE)

ggplot(preds_depth, aes(x = Depth, y = exp(est), color = reorder(sh_factor, sh))) +
  geom_line(linewidth = 1.3) +
  geom_ribbon(aes(ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se), fill = reorder(sh_factor, sh)), alpha = 0.15, color = NA) +
  scale_color_viridis_d(option = "plasma", name = "Scallop Size") +
  scale_fill_viridis_d(option = "plasma", name = "Scallop Size") +
  theme_minimal() +
  labs(title = "Effect of Depth on Predicted Scallop Meat Weight", x = "Water Depth (m)", y = "Predicted Meat Weight (g)")

# 2. Julian Day Effect
jday_seq <- seq(min(df_clean$jday), max(df_clean$jday), length.out = 100)

pred_grid_jday <- expand.grid(
  jday = jday_seq, sh = sh_targets,
  Depth = mean(df_clean$Depth, na.rm = TRUE),
  log_density_c = 0, year = min(df_clean$year)
) %>%
  mutate(
    log_sh = log(sh),
    log_sh_c = log_sh - mean(log(df_clean$sh), na.rm = TRUE),
    sh_factor = paste0("SH = ", sh, " mm")
  )

preds_jday <- predict(fit_m7_svc, newdata = pred_grid_jday, re_form = NA, se_fit = TRUE)

ggplot(preds_jday, aes(x = jday, y = exp(est), color = reorder(sh_factor, sh))) +
  geom_line(linewidth = 1.3) +
  geom_ribbon(aes(ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se), fill = reorder(sh_factor, sh)), alpha = 0.15, color = NA) +
  scale_color_viridis_d(option = "plasma", name = "Scallop Size") +
  scale_fill_viridis_d(option = "plasma", name = "Scallop Size") +
  theme_minimal() +
  labs(title = "Seasonal Effect of Julian Day on Predicted Meat Weight", x = "Julian Day (jday)", y = "Predicted Meat Weight (g)")

# 3. Local Density Effect
density_seq <- seq(min(df_clean$log_density_c), max(df_clean$log_density_c), length.out = 100)

pred_grid_density <- expand.grid(
  log_density_c = density_seq, sh = sh_targets,
  Depth = mean(df_clean$Depth, na.rm = TRUE),
  jday = mean(df_clean$jday, na.rm = TRUE),
  year = min(df_clean$year)
) %>%
  mutate(
    log_sh = log(sh),
    log_sh_c = log_sh - mean(log(df_clean$sh), na.rm = TRUE),
    sh_factor = paste0("SH = ", sh, " mm"),
    raw_density = exp((log_density_c * sd(df_clean$log_density)) + mean(df_clean$log_density))
  )

preds_density <- predict(fit_m7_svc, newdata = pred_grid_density, re_form = NA, se_fit = TRUE)

ggplot(preds_density, aes(x = raw_density, y = exp(est), color = reorder(sh_factor, sh))) +
  geom_line(linewidth = 1.3) +
  geom_ribbon(aes(ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se), fill = reorder(sh_factor, sh)), alpha = 0.15, color = NA) +
  scale_color_viridis_d(option = "plasma", name = "Scallop Size") +
  scale_fill_viridis_d(option = "plasma", name = "Scallop Size") +
  theme_minimal() +
  labs(title = "Effect of Local Density on Predicted Scallop Meat Weight", x = "Station Scallop Count (Density Proxy)", y = "Predicted Meat Weight (g)")

Spatial and spatialtemporal variation

# ==============================================================================
# Spatial and Spatiotemporal Variation Analysis for fit_m7_svc
# ==============================================================================

library(sdmTMB)
library(dplyr)
library(ggplot2)
library(viridis)
library(patchwork)

# ------------------------------------------------------------------------------
# 1. Generate Model Predictions Including Spatial & Spatiotemporal Random Fields
# ------------------------------------------------------------------------------
# Predict on data locations (or pass a regular spatial grid in newdata)
preds_m7 <- predict(fit_m7_svc)

# ------------------------------------------------------------------------------
# 2. Visualize Fixed Spatial Random Fields (omega_s and zeta_s_log_sh_c)
# ------------------------------------------------------------------------------

# Spatial Intercept Variation (omega_s)
p_omega <- ggplot(preds_m7, aes(x = X, y = Y, color = omega_s)) +
  geom_point(size = 1.2, alpha = 0.8) +
  scale_color_viridis_c(option = "mako", name = "Omega (s)") +
  coord_fixed() +
  theme_minimal() +
  labs(
    title = "Spatial Intercept Field (omega_s)",
    subtitle = "Persistent spatial variation in baseline condition",
    x = "UTM X (km)", y = "UTM Y (km)"
  )

# Spatially Varying Slope (zeta_s_log_sh_c)
p_zeta <- ggplot(preds_m7, aes(x = X, y = Y, color = zeta_s_log_sh_c)) +
  geom_point(size = 1.2, alpha = 0.8) +
  scale_color_viridis_c(option = "plasma", name = "Zeta (s)") +
  coord_fixed() +
  theme_minimal() +
  labs(
    title = "Spatially Varying Slope Field (zeta_s_log_sh_c)",
    subtitle = "Spatial variation in shell height allometric slope",
    x = "UTM X (km)", y = "UTM Y (km)"
  )

# Display Spatial Fields Side-by-Side
p_omega + p_zeta

# ------------------------------------------------------------------------------
# 3. Visualize Spatiotemporal Variation Faceted by Year (epsilon_st)
# ------------------------------------------------------------------------------

p_epsilon <- ggplot(preds_m7, aes(x = X, y = Y, color = epsilon_st)) +
  geom_point(size = 1, alpha = 0.8) +
  scale_color_gradient2(
    low = "#0571b0", mid = "white", high = "#ca0020",
    midpoint = 0, name = "Epsilon (st)"
  ) +
  facet_wrap(~ year, ncol = 6) +
  coord_fixed() +
  theme_bw() +
  theme(
    panel.background = element_rect(fill = "white"),
    panel.grid = element_blank()
  ) +
  labs(
    title = "Spatiotemporal Random Fields (epsilon_st - AR1 Process)",
    subtitle = "Year-to-year dynamic spatial deviations in model fit",
    x = "UTM X (km)", y = "UTM Y (km)"
  )

# Print Spatiotemporal Plot
print(p_epsilon)

# ------------------------------------------------------------------------------
# 4. Extract Variance Components & Spatial Parameters
# ------------------------------------------------------------------------------
cat("\n=== Spatial & Spatiotemporal Variance Parameters ===\n")

=== Spatial & Spatiotemporal Variance Parameters ===

tidy(fit_m7_svc, effects = "ran_pars", conf.int = TRUE)

# A tibble: 8 × 5
  term         estimate std.error conf.low conf.high
  <chr>           <dbl>     <dbl>    <dbl>     <dbl>
1 range         52.9      2.81     47.7      58.7   
2 phi           29.1      0.178    28.7      29.4   
3 sigma_O        0.0740   0.00759   0.0605    0.0905
4 sigma_E        0.151    0.00475   0.142     0.161 
5 sigma_Z        0.417    0.0355    0.353     0.493 
6 rho            0.294   NA         0.216     0.368 
7 sd__s(Depth)   0.552   NA         0.311     0.978 
8 sd__s(jday)    0.106   NA         0.0350    0.324 

Regional Parameter Deviations (SAMS Areas)

# ==============================================================================
# SAMS Area Shell Height - Meat Weight (SH-MW) Allometric Relationships
# Model: fit_m7_svc
# ==============================================================================

# ------------------------------------------------------------------------------
# 1. Generate Model Predictions & Extract Spatial Random Fields
# ------------------------------------------------------------------------------
# Predict on dataset to get spatial random fields for fit_m7_svc
preds_m7 <- predict(fit_m7_svc)

# ------------------------------------------------------------------------------
# 2. Extract Global Model Parameters
# ------------------------------------------------------------------------------
# Mean log shell height used for centering
mean_log_sh <- mean(log(df_clean$sh), na.rm = TRUE)

# Global fixed-effect slope for centered log shell height (log_sh_c)
beta_slope <- tidy(fit_m7_svc) %>% 
  filter(term == "log_sh_c") %>% 
  pull(estimate)

# Global predicted log(MW) at baseline covariate levels (log_sh_c = 0)
standard_data <- data.frame(
  log_sh_c      = 0,
  Depth         = mean(df_clean$Depth, na.rm = TRUE),
  jday          = mean(df_clean$jday, na.rm = TRUE),
  log_density_c = 0,
  year          = min(df_clean$year)
)

C_standard <- predict(fit_m7_svc, newdata = standard_data, re_form = NA)$est

# ------------------------------------------------------------------------------
# 3. Calculate Regional SAMS Parameters (a and b in MW = a * SH^b)
# ------------------------------------------------------------------------------
sams_params_m7 <- preds_m7 %>%
  filter(!is.na(Sams)) %>%
  group_by(Sams) %>%
  summarize(
    mean_omega = mean(omega_s, na.rm = TRUE),            # Spatial intercept deviation
    mean_zeta  = mean(zeta_s_log_sh_c, na.rm = TRUE),   # Spatial slope deviation
    n_obs      = n(),
    .groups    = "drop"
  ) %>%
  filter(n_obs > 30) %>% # Filter out areas with insufficient sample size
  mutate(
    b = beta_slope + mean_zeta,                          # Allometric slope exponent (b)
    log_a = C_standard + mean_omega - (b * mean_log_sh), # Intercept back-transform
    a = exp(log_a),                                      # Multiplicative coefficient (a)
    equation_label = sprintf("MW = %.2e * SH^%.3f", a, b)
  )

# Print SAMS Parameter Table
cat("\n=== SAMS Area Regional SH-MW Parameters (fit_m7_svc) ===\n")

=== SAMS Area Regional SH-MW Parameters (fit_m7_svc) ===

knitr::kable(
  sams_params_m7 %>% select(Sams, n_obs, a, b, equation_label), 
  digits = 4, 
  caption = "Model-Derived SAMS Regional Parameters (MW = a * SH^b)"
)

Model-Derived SAMS Regional Parameters (MW = a * SH^b)

Sams	n_obs	a	b	equation_label
	39	1e-04	2.7314	MW = 6.15e-05 * SH^2.731
BI	1312	1e-04	2.6325	MW = 1.00e-04 * SH^2.632
DMV	2382	2e-04	2.4747	MW = 2.02e-04 * SH^2.475
ET	13681	1e-04	2.6951	MW = 7.61e-05 * SH^2.695
HCS	10872	1e-04	2.7307	MW = 6.37e-05 * SH^2.731
Inshore	1709	1e-04	2.7527	MW = 6.18e-05 * SH^2.753
LI	18014	0e+00	2.8333	MW = 3.99e-05 * SH^2.833
NYB	8199	0e+00	2.8412	MW = 3.92e-05 * SH^2.841
VIR	217	0e+00	2.9969	MW = 1.71e-05 * SH^2.997

# ------------------------------------------------------------------------------
# 4. Plot Regional Growth Curves
# ------------------------------------------------------------------------------
# Create continuous prediction grid across SAMS areas
grid_sams_m7 <- expand.grid(
  Sams = sams_params_m7$Sams, 
  sh   = seq(40, 180, length.out = 200)
) %>%
  left_join(sams_params_m7, by = "Sams") %>%
  mutate(mw_pred = a * (sh ^ b))

# Global reference curve
grid_global <- data.frame(sh = seq(40, 180, length.out = 200)) %>%
  mutate(
    b_global = beta_slope,
    a_global = exp(C_standard - (beta_slope * mean_log_sh)),
    mw_pred  = a_global * (sh ^ b_global)
  )

# Plot SAMS curves overlaid with global average
ggplot(grid_sams_m7, aes(x = sh, y = mw_pred, color = Sams)) +
  geom_line(linewidth = 1.2) +
  geom_line(
    data = grid_global, aes(x = sh, y = mw_pred), 
    color = "black", linewidth = 1.5, linetype = "dashed", inherit.aes = FALSE
  ) +
  scale_color_viridis_d(option = "turbo") +
  theme_minimal() +
  labs(
    title = "SAMS Regional Allometric Growth Curves (fit_m7_svc)",
    subtitle = "Dashed black line indicates global model average across regions",
    x = "Shell Height (mm)",
    y = "Predicted Meat Weight (g)",
    color = "SAMS Region"
  )