Fin Whale Group-Size GP Model - Diagnostics

Model Overview

This document contains diagnostic plots and summary statistics for the fin whale joint visual-acoustic density model with group-size likelihood structure.

Model features: - Hilbert-space Gaussian Process (HSGP) for spatial smoothing - Separate GP per season (Winter, Spring, Summer, Fall) - Group-based visual likelihood (models number of groups, not individuals) - Poisson model for group size distribution - Joint acoustic likelihood (20Hz and 40Hz calls)

---

Load Libraries and Data

library(rstan)
library(dplyr)
library(tidyr)
library(ggplot2)
library(stringr)
library(bayesplot)
library(posterior)
library(gridExtra)

# Set paths
code_dir <- "/Users/michaelaalksne/Desktop/PhD/Chapter2/CalCOFI/CalCOFI_integrated_SDM/code/stan/"
data_dir <- "/Users/michaelaalksne/Desktop/PhD/Chapter2/CalCOFI/CalCOFI_integrated_SDM/data/"
fig_dir  <- "/Users/michaelaalksne/Desktop/PhD/Chapter2/CalCOFI/CalCOFI_integrated_SDM/figures/modeled_outputs/fin_gp_group_size/"

season_labels <- c("1"="Winter", "2"="Spring", "3"="Summer", "4"="Fall")

# Load fitted model
fit <- readRDS(paste0(code_dir, "fitted_finwhale_gp_groupsize.rds"))
posterior_samples <- readRDS(paste0(code_dir, "posterior_samples_finwhale_gp_groupsize.rds"))

# Load original data for comparison
sighting_raw <- readr::read_csv(paste0(data_dir, "Bp_sighting_data.csv"), 
                                show_col_types = FALSE)
subsegment_raw <- readr::read_csv(paste0(data_dir, "Bp_groupsize_data.csv"),
                                  show_col_types = FALSE)

cat("Posterior samples:", nrow(posterior_samples), "\n")

## Posterior samples: 2000

cat("Total sightings:", nrow(sighting_raw), "\n")

## Total sightings: 288

cat("Total subsegments:", nrow(subsegment_raw), "\n")

## Total subsegments: 10694

---

1. Model Summary

cat("FIN WHALE GROUP-SIZE GP MODEL - SUMMARY\n")

## FIN WHALE GROUP-SIZE GP MODEL - SUMMARY

cat("DATA:\n")

## DATA:

cat("  Visual subsegments:", nrow(subsegment_raw), "\n")

##   Visual subsegments: 10694

cat("  Subsegments with groups:", sum(subsegment_raw$num_groups_observed > 0), "\n")

##   Subsegments with groups: 213

cat("  Total sightings:", nrow(sighting_raw), "\n")

##   Total sightings: 288

cat("  Mean observed group size:", round(mean(sighting_raw$group_size), 2), "\n\n")

##   Mean observed group size: 1.92

# Get convergence info early
summary_df <- as.data.frame(summary(fit)$summary)
rhat_issues <- rownames(summary_df)[summary_df$Rhat > 1.01 & !is.nan(summary_df$Rhat)]

cat("CONVERGENCE:\n")

## CONVERGENCE:

cat("  Rhat > 1.01:", length(rhat_issues), "parameters\n")

##   Rhat > 1.01: 0 parameters

cat("  Status:", ifelse(length(rhat_issues) == 0, "✓ CONVERGED", "⚠ CHECK"), "\n\n")

##   Status: ✓ CONVERGED

# Get key posteriors
exp_group_size_posterior <- posterior_samples$exp_group_size
seasonal_dens <- data.frame(
  Season = factor(season_labels, levels = c("Winter", "Spring", "Summer", "Fall")),
  mean_density = colMeans(as.matrix(posterior_samples[, paste0("mean_density_season.", 1:4)])),
  lower_95 = apply(as.matrix(posterior_samples[, paste0("mean_density_season.", 1:4)]), 
                   2, quantile, 0.025),
  upper_95 = apply(as.matrix(posterior_samples[, paste0("mean_density_season.", 1:4)]), 
                   2, quantile, 0.975)
)

rho_long <- bind_rows(lapply(1:4, function(m) {
  data.frame(
    rho = posterior_samples[[paste0("rho.", m)]],
    Season = season_labels[as.character(m)]
  )
})) %>%
  mutate(Season = factor(Season, levels = c("Winter", "Spring", "Summer", "Fall")))

cat("KEY PARAMETERS:\n")

## KEY PARAMETERS:

cat("  Expected group size:", round(mean(exp_group_size_posterior), 2),
    " (95% CI: ", round(quantile(exp_group_size_posterior, 0.025), 2), "-",
    round(quantile(exp_group_size_posterior, 0.975), 2), ")\n", sep = "")

##   Expected group size:1.93 (95% CI: 1.82-2.05)

cat("\n  Seasonal density (whales/km²):\n")

## 
##   Seasonal density (whales/km²):

for (i in 1:4) {
  cat("    ", season_labels[as.character(i)], ": ",
      sprintf("%.5f", seasonal_dens$mean_density[i]), "\n", sep = "")
}

##     Winter: 0.00138
##     Spring: 0.00470
##     Summer: 0.00863
##     Fall: 0.00412

cat("\n  GP length scales (km):\n")

## 
##   GP length scales (km):

rho_summary <- rho_long %>%
  group_by(Season) %>%
  summarise(mean_rho = mean(rho), .groups = "drop")
for (i in 1:4) {
  cat("rho summary", rho_summary$Season[i], ": ",
      round(rho_summary$mean_rho[i], 1), " km\n", sep = "")
}

## rho summary1: 151.4 km
## rho summary2: 124.3 km
## rho summary3: 131.6 km
## rho summary4: 91.6 km

measure of distance of spatial correlation rho = “How far apart can two locations be before their densities become uncorrelated?”

rho = 90 km means: - Two locations 20 km apart: HIGHLY correlated (similar density) - Two locations 90 km apart: MODERATELY correlated (the threshold) - Two locations 200 km apart: UNCORRELATED (independent)

Bigger rho = smoother surface (correlation extends further) Smaller rho = bumpier surface (density varies on shorter scales)

interpretation:

Winter (151 km) - Most Dispersed: Whales are spread out in small groups across the SCB Density varies on large spatial scales Fewer distinct hotspots, more gradual gradients

Spring (124 km) - Intermediate: Some aggregation but not as tight as Fall. Similar to winter/summer.

Summer (131 km) - Intermediate: Some aggregation but not as tight as Fall. Similar to winter/spring. However, we know that there are WAY more whales during this time of year. Interesting that spatial correlation remains similar. High densities are kind of all over the place cause there’s so many whales?

Fall (92 km) - Most Localized: Fin whales appear to be in more concentrated hotspots Density varies on smaller spatial scales local aggregations separated by empty areas Matches the acoustics (lots of 20Hz in specific areas!)

---

2. Convergence Diagnostics

2.1 Summary Statistics

# Print summary for key parameters
print(fit, pars = c("mu_intercept", "gp_alpha", "rho",
                    "beta_20", "beta_40", "mu_group_size", "exp_group_size",
                    "phi_v", "phi_20", "phi_40",
                    "mean_density_season", "spatial_sd"))

## Inference for Stan model: anon_model.
## 4 chains, each with iter=1000; warmup=500; thin=1; 
## post-warmup draws per chain=500, total post-warmup draws=2000.
## 
##                          mean se_mean    sd  2.5%    25%    50%    75%  97.5%
## mu_intercept[1]         -7.41    0.01  0.23 -7.86  -7.56  -7.40  -7.25  -6.96
## mu_intercept[2]         -7.40    0.01  0.24 -7.87  -7.55  -7.39  -7.24  -6.95
## mu_intercept[3]         -6.64    0.00  0.17 -6.99  -6.75  -6.64  -6.53  -6.32
## mu_intercept[4]         -6.98    0.01  0.23 -7.46  -7.14  -6.98  -6.83  -6.54
## gp_alpha[1]              1.31    0.02  0.43  0.75   1.02   1.23   1.50   2.35
## gp_alpha[2]              1.58    0.01  0.48  0.92   1.25   1.49   1.80   2.63
## gp_alpha[3]              1.89    0.02  0.61  1.14   1.49   1.77   2.14   3.41
## gp_alpha[4]              1.83    0.02  0.52  1.05   1.46   1.76   2.11   3.10
## rho[1]                 151.43    2.31 59.88 68.45 112.92 139.16 176.56 304.06
## rho[2]                 124.34    1.89 45.19 55.74  95.61 117.66 144.51 247.36
## rho[3]                 131.56    1.19 30.47 79.30 111.19 128.66 149.22 198.02
## rho[4]                  91.63    1.05 26.65 43.73  73.60  90.14 107.55 145.28
## beta_20[1]               0.30    0.00  0.11  0.15   0.23   0.28   0.36   0.56
## beta_20[2]               0.15    0.00  0.05  0.08   0.11   0.14   0.17   0.27
## beta_20[3]               0.01    0.00  0.00  0.01   0.01   0.01   0.02   0.02
## beta_20[4]               0.31    0.00  0.10  0.17   0.24   0.30   0.37   0.54
## beta_40[1]               0.00    0.00  0.00  0.00   0.00   0.00   0.00   0.01
## beta_40[2]               0.01    0.00  0.00  0.00   0.00   0.01   0.01   0.01
## beta_40[3]               0.01    0.00  0.00  0.00   0.01   0.01   0.01   0.02
## beta_40[4]               0.00    0.00  0.00  0.00   0.00   0.00   0.00   0.00
## mu_group_size            0.93    0.00  0.06  0.82   0.89   0.93   0.97   1.05
## exp_group_size           1.93    0.00  0.06  1.82   1.89   1.93   1.97   2.05
## phi_v                    0.11    0.00  0.02  0.08   0.09   0.11   0.12   0.15
## phi_20                   0.03    0.00  0.00  0.03   0.03   0.03   0.03   0.04
## phi_40                   0.02    0.00  0.00  0.02   0.02   0.02   0.02   0.02
## mean_density_season[1]   0.00    0.00  0.00  0.00   0.00   0.00   0.00   0.00
## mean_density_season[2]   0.00    0.00  0.13  0.00   0.00   0.00   0.00   0.00
## mean_density_season[3]   0.01    0.00  0.03  0.00   0.00   0.00   0.01   0.03
## mean_density_season[4]   0.00    0.00  0.01  0.00   0.00   0.00   0.00   0.01
## spatial_sd[1]            1.21    0.01  0.27  0.78   1.03   1.18   1.34   1.84
## spatial_sd[2]            1.47    0.01  0.34  0.98   1.24   1.42   1.64   2.24
## spatial_sd[3]            1.74    0.02  0.41  1.16   1.46   1.67   1.95   2.78
## spatial_sd[4]            1.64    0.01  0.38  1.01   1.36   1.60   1.86   2.50
##                        n_eff Rhat
## mu_intercept[1]         1715 1.00
## mu_intercept[2]         1249 1.00
## mu_intercept[3]         1831 1.00
## mu_intercept[4]         1284 1.00
## gp_alpha[1]              703 1.01
## gp_alpha[2]             1049 1.00
## gp_alpha[3]              591 1.01
## gp_alpha[4]              611 1.01
## rho[1]                   669 1.00
## rho[2]                   571 1.01
## rho[3]                   652 1.01
## rho[4]                   647 1.00
## beta_20[1]              1579 1.00
## beta_20[2]              1686 1.00
## beta_20[3]              2308 1.00
## beta_20[4]              1962 1.00
## beta_40[1]              2028 1.00
## beta_40[2]              1995 1.00
## beta_40[3]              2014 1.00
## beta_40[4]              1969 1.00
## mu_group_size           3142 1.00
## exp_group_size          3142 1.00
## phi_v                   2552 1.00
## phi_20                  3670 1.00
## phi_40                  3426 1.00
## mean_density_season[1]  1073 1.00
## mean_density_season[2]  2008 1.00
## mean_density_season[3]  1312 1.00
## mean_density_season[4]  1643 1.00
## spatial_sd[1]           1228 1.00
## spatial_sd[2]           1090 1.00
## spatial_sd[3]            710 1.01
## spatial_sd[4]            842 1.00
## 
## Samples were drawn using NUTS(diag_e) at Wed Mar 11 01:41:04 2026.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

2.2 Rhat Diagnostics

cat("\nParameters with Rhat > 1.005:", length(rhat_issues), "\n")

## 
## Parameters with Rhat > 1.005: 0

if (length(rhat_issues) > 0) {
  cat("\nProblematic parameters:\n")
  print(summary_df[rhat_issues, c("mean", "sd", "Rhat", "n_eff")])
} else {
  cat("✓ All Rhat < 1.005 -  convergence!\n")
}

## ✓ All Rhat < 1.005 -  convergence!

# Rhat histogram
ggplot(data.frame(Rhat = summary_df$Rhat), aes(x = Rhat)) +
  geom_histogram(bins = 50, fill = "steelblue", alpha = 0.7) +
  geom_vline(xintercept = 1.005, color = "red", linetype = "dashed") +
  labs(title = "Rhat Distribution Across All Parameters",
       subtitle = "Red line at 1.005 threshold",
       x = "Rhat", y = "Count") +
  theme_minimal()

R_hat - a measure of chain equilibrium, must be within 0.05 of 1.0 looks good

2.3 Effective Sample Size

ggplot(data.frame(n_eff = summary_df$n_eff), aes(x = n_eff)) +
  geom_histogram(bins = 50, fill = "darkgreen", alpha = 0.7) +
  geom_vline(xintercept = 400, color = "red", linetype = "dashed") +
  labs(title = "Effective Sample Size Distribution",
       subtitle = "Red line at 400",
       x = "Effective Sample Size", y = "Count") +
  theme_minimal()

I ran: 4 chains × 500 post-warmup iterations = 2000 total samples

N_eff - effective sample size - the number of independent draws in the sample

n_eff ≈ 2000 means: - almost NO autocorrelation in the chains - Each iteration is providing independent information - good mixing and exploration - efficient sampling

---

3. Trace Plots

3.1 Density Intercepts and GP Amplitude

mcmc_trace(fit, pars = c("mu_intercept[1]", "mu_intercept[2]",
                          "mu_intercept[3]", "mu_intercept[4]",
                          "gp_alpha[1]", "gp_alpha[2]",
                          "gp_alpha[3]", "gp_alpha[4]")) +
  labs(title = "Trace Plots - Density Intercepts and GP Amplitude") +
  theme_minimal()

All chains mixed together, can’t tell them apart Looks like thick fuzzy horizontal bands No trends, no patterns, just random wiggling around a stable value. Looks good.

Example: ───────────────── ← all 4 chains overlapping ───────────────── ─────────────────

3.2 Length Scale (rho) and Call Rates (beta)

mcmc_trace(fit, pars = c("rho[1]", "rho[2]", "rho[3]", "rho[4]",
                          "beta_20[1]", "beta_20[2]", "beta_20[3]", "beta_20[4]",
                          "beta_40[1]", "beta_40[2]", "beta_40[3]", "beta_40[4]")) +
  labs(title = "Trace Plots - GP Length Scale and Call Rates") +
  theme_minimal()

All chains mixed together, can’t tell them apart

3.3 Group Size and Overdispersion

mcmc_trace(fit, pars = c("mu_group_size", "exp_group_size",
                          "phi_v", "phi_20", "phi_40")) +
  labs(title = "Trace Plots - Group Size and Overdispersion") +
  theme_minimal()

All chains mixed together, can’t tell them apart

---

4. Group Size Parameter

4.1 Prior vs Posterior

# Extract posterior
exp_group_size_posterior <- posterior_samples$exp_group_size

# Generate prior samples
set.seed(456)
n_prior_samples <- length(exp_group_size_posterior)
mu_group_size_prior <- rexp(n_prior_samples, rate = 1)
exp_group_size_prior <- mu_group_size_prior + 1

# Combine for plotting
group_size_comparison <- bind_rows(
  data.frame(group_size = exp_group_size_prior, distribution = "Prior"),
  data.frame(group_size = exp_group_size_posterior, distribution = "Posterior")
) %>%
  mutate(distribution = factor(distribution, levels = c("Prior", "Posterior")))

# Observed means
observed_mean <- mean(sighting_raw$group_size)


ggplot(group_size_comparison, aes(x = group_size, fill = distribution, 
                                   alpha = distribution)) +
  geom_density(color = "black", linewidth = 0.3) +
  geom_vline(xintercept = observed_mean, color = "red", 
             linetype = "dashed", linewidth = 1) +
  annotate("text", x = observed_mean, y = 0.5, 
           label = paste0("Observed mean: ", round(observed_mean, 2)),
           color = "red", hjust = -0.1, size = 3.5) +
  scale_fill_manual(values = c("Prior" = "gray50", "Posterior" = "steelblue")) +
  scale_alpha_manual(values = c("Prior" = 0.4, "Posterior" = 0.8)) +
  labs(title = "Expected Group Size: Prior vs Posterior",
       subtitle = "Prior: exponential(1) + 1 | Posterior learned from 288 sightings",
       x = "Expected group size (whales per group)",
       y = "Density",
       fill = NULL, alpha = NULL) +
  theme_minimal() +
  theme(plot.title = element_text(size = 14, face = "bold"),
        legend.position = "top")

This shows the comparison between the prior distribution we provided and the posterior that the model learned. Also with the observed mean from the data, which is basically what the model learned. Our prior was exp(1), to represent the long-tail distribution of the data. This plot shows that the model learned that the true underlying distribution of group sizes was not really that, rather, it was tightly centered around 2 whales per group, and the larger numbers are just outliers….

4.2 Group Size Summary Statistics

cat("OBSERVED GROUP SIZES (from data):\n")

## OBSERVED GROUP SIZES (from data):

cat("  N sightings:", nrow(sighting_raw), "\n")

##   N sightings: 288

cat("  Mean:", round(mean(sighting_raw$group_size), 2), "\n")

##   Mean: 1.92

cat("  Median:", median(sighting_raw$group_size), "\n")

##   Median: 1

cat("  Range:", min(sighting_raw$group_size), "-", max(sighting_raw$group_size), "\n")

##   Range: 1 - 30

cat("  SD:", round(sd(sighting_raw$group_size), 2), "\n\n")

##   SD: 2.6

cat("POSTERIOR exp_group_size:\n")

## POSTERIOR exp_group_size:

cat("  Mean:", round(mean(exp_group_size_posterior), 2), "\n")

##   Mean: 1.93

cat("  Median:", round(median(exp_group_size_posterior), 2), "\n")

##   Median: 1.93

cat("  95% CI: [", round(quantile(exp_group_size_posterior, 0.025), 2), ", ",
    round(quantile(exp_group_size_posterior, 0.975), 2), "]\n", sep = "")

##   95% CI: [1.82, 2.05]

cat("  SD:", round(sd(exp_group_size_posterior), 3), "\n")

##   SD: 0.057

4.3 Group Size by Number of Groups per Subsegment

sighting_raw %>%
  left_join(subsegment_raw %>% select(SubSegID, num_groups_observed),
            by = "SubSegID") %>%
  group_by(num_groups_observed) %>%
  summarise(
    n_sightings = n(),
    mean_group_size = mean(group_size),
    sd_group_size = sd(group_size),
    .groups = "drop"
  ) %>%
  ggplot(aes(x = factor(num_groups_observed), y = mean_group_size)) +
  geom_col(fill = "steelblue", alpha = 0.7) +
  geom_errorbar(aes(ymin = pmax(0, mean_group_size - sd_group_size),
                  ymax = mean_group_size + sd_group_size),
              width = 0.2)+
  geom_text(aes(label = paste0("n=", n_sightings)), 
            vjust = -0.5, size = 3) +
  geom_hline(yintercept = mean(exp_group_size_posterior), 
             color = "red", linetype = "dashed") +
  labs(title = "Mean Group Size by Subsegment Density",
       subtitle = "Red line = posterior mean exp_group_size",
       x = "Number of groups in subsegment",
       y = "Mean group size (whales/group)") +
  theme_minimal()

This shows that for small numbers of sightings on a subsegment, as we increase # of sightings from 1 to 2, we increase the mean group size (with a lot of variability). But then as we go up past that, to 3-7 sightings on a subsegment, we don’t have a consistent increase in the number of groups. These also become more rare as we go on. We also see that the posterior mean is really close to the mean group size of the data. Which is 2! We usually see fin whales in group sizes of 2 whales.

---

5. GP Parameters

5.1 Length Scale (rho) by Season

ggplot(rho_long, aes(x = rho, fill = Season)) +
  geom_density(alpha = 0.6, color = NA) +
  geom_vline(xintercept = 100, linetype = "dashed", 
             color = "black", alpha = 0.5) +
  facet_wrap(~Season, scales = "free_y", ncol = 2) +
  scale_fill_manual(values = c("Winter" = "steelblue", "Spring" = "darkgreen",
                                "Summer" = "orange", "Fall" = "coral")) +
  labs(title = "GP Length Scale (rho) by Season",
       subtitle = "Spatial correlation range (km) | Dashed line = prior mean 100 km",
       x = "rho (km)", y = "Density") +
  theme_minimal() +
  theme(plot.title = element_text(size = 14, face = "bold"),
        legend.position = "none")

Basically shows the same thing as above but now as a distribution. The model clearly learned something different than the prior we gave it (100 km for all seasons). I find it interested/ a bit unexpected that summer and winter are similar although their densities are so different. Like same underlying distributed spatial structure, but WAY more whales in the summer. But then something different happens in the fall. Hypothesis: females leave but males stay in their territories to sing…? We know this is peak singing season and the summer/fall spatial distribution matches, density is just way lower in the fall.

5.2 GP Amplitude by Season

alpha_long <- bind_rows(lapply(1:4, function(m) {
  data.frame(
    alpha = posterior_samples[[paste0("gp_alpha.", m)]],
    Season = season_labels[as.character(m)]
  )
})) %>%
  mutate(Season = factor(Season, levels = c("Winter", "Spring", "Summer", "Fall")))

ggplot(alpha_long, aes(x = alpha, fill = Season)) +
  geom_density(alpha = 0.6, color = NA) +
  facet_wrap(~Season, scales = "free_y", ncol = 2) +
  scale_fill_manual(values = c("Winter" = "steelblue", "Spring" = "darkgreen",
                                "Summer" = "orange", "Fall" = "coral")) +
  labs(title = "GP Amplitude (alpha) by Season",
       subtitle = "Magnitude of spatial variation",
       x = "GP amplitude", y = "Density") +
  theme_minimal() +
  theme(plot.title = element_text(size = 14, face = "bold"),
        legend.position = "none")

alpha_long %>%
  group_by(Season) %>%
  summarise(mean_alpha = mean(alpha))

## # A tibble: 4 × 2
##   Season mean_alpha
##   <fct>       <dbl>
## 1 Winter       1.31
## 2 Spring       1.58
## 3 Summer       1.89
## 4 Fall         1.83

GP alpha: standard deviation of spatial variation

Winter: 1.31 ← LOWEST - most uniform distribution Spring: 1.58 ← moderate patchiness Summer: 1.89 ← HIGH - very patchy (the most variability) Fall: 1.83 ← HIGH - very patchy

Summer has: - High alpha (1.89) = BIGGEST differences in density across space - Long rho (131 km) = However, those differences change more GRADUALLY over longer distances - There’s high density aggregations all over the SCB?

Fall has: - High alpha (1.83) = BIG differences in density across space - Long rho (92 km) = Those differences change more RAPIDLY over shorter distances - There’s only a few high density aggregations? So still a lot of variance, but less spatial correlation between them.

Winter: - Lower alpha (1.31) = Smaller differences in density across space - Long rho (150 km) = Those differences change more gradually over longer distances - Makes sense, less whales overall, low density aggregations.

Spring: - Lower alpha (1.58) = Intermediate differences in density across space - Long rho (124 km) = Intermediate. (Kind of the build up period, influx of whales?)

---

6. Call Rate Parameters (Beta)

6.1 Seasonal Beta Posteriors

beta_long <- bind_rows(lapply(1:4, function(m) {
  bind_rows(
    data.frame(
      Value = posterior_samples[[paste0("beta_20.", m)]],
      call_type = "20Hz",
      Season = season_labels[as.character(m)]
    ),
    data.frame(
      Value = posterior_samples[[paste0("beta_40.", m)]],
      call_type = "40Hz",
      Season = season_labels[as.character(m)]
    )
  )
})) %>%
  mutate(Season = factor(Season, levels = c("Winter", "Spring", "Summer", "Fall")))

ggplot(beta_long, aes(x = Value, fill = call_type)) +
  geom_density(alpha = 0.6, color = NA) +
  facet_grid(call_type ~ Season, scales = "free") +
  scale_fill_manual(values = c("20Hz" = "steelblue", "40Hz" = "coral")) +
  labs(title = "Seasonal Call Rate (beta) Posteriors",
       subtitle = "Calls per whale per hour per km²",
       x = "Beta", y = "Density") +
  theme_minimal() +
  theme(plot.title = element_text(size = 14, face = "bold"),
        legend.position = "none")

6.2 Beta Summary Table

beta_summary <- beta_long %>%
  group_by(Season, call_type) %>%
  summarise(
    median = median(Value),
    mean = mean(Value),
    lower_95 = quantile(Value, 0.025),
    upper_95 = quantile(Value, 0.975),
    .groups = "drop"
  ) %>%
  arrange(call_type, Season)

print(beta_summary, digits = 4)

## # A tibble: 8 × 6
##   Season call_type  median    mean lower_95 upper_95
##   <fct>  <chr>       <dbl>   <dbl>    <dbl>    <dbl>
## 1 Winter 20Hz      0.283   0.302    0.146    0.563  
## 2 Spring 20Hz      0.136   0.145    0.0752   0.269  
## 3 Summer 20Hz      0.0135  0.0139   0.00840  0.0222 
## 4 Fall   20Hz      0.297   0.312    0.167    0.537  
## 5 Winter 40Hz      0.00361 0.00391  0.00178  0.00766
## 6 Spring 40Hz      0.00632 0.00689  0.00313  0.0141 
## 7 Summer 40Hz      0.00834 0.00892  0.00441  0.0167 
## 8 Fall   40Hz      0.00225 0.00245  0.00113  0.00491

6.3 Beta Prior vs Posterior

# Prior parameters (from Stan model)
prior_params_fw <- tibble(
  Season = rep(c("Winter", "Spring", "Summer", "Fall"), each = 2),
  call_type = rep(c("20Hz", "40Hz"), 4),
  meanlog = log(c(
    0.1026, 0.00209,  # Winter
    0.0731, 0.00467,  # Spring
    0.0221, 0.00830,  # Summer
    0.2075, 0.00114   # Fall
  )),
  sdlog = 1.0
)

# Generate prior samples
set.seed(789)
n_samples <- nrow(posterior_samples)
beta_prior_samples <- prior_params_fw %>%
  rowwise() %>%
  mutate(Value = list(rlnorm(n_samples, meanlog, sdlog))) %>%
  unnest(Value) %>%
  select(Season, call_type, Value) %>%
  mutate(distribution = "Prior")

# Combine with posteriors
beta_comparison <- bind_rows(
  beta_long %>% mutate(distribution = "Posterior"),
  beta_prior_samples
) %>%
  mutate(
    Season = factor(Season, levels = c("Winter", "Spring", "Summer", "Fall")),
    distribution = factor(distribution, levels = c("Prior", "Posterior"))
  )

ggplot(beta_comparison, aes(x = Value, fill = distribution, 
                             alpha = distribution)) +
  geom_density(color = "black", linewidth = 0.3) +
  facet_grid(call_type ~ Season, scales = "free") +
  scale_fill_manual(values = c("Prior" = "gray50", "Posterior" = "steelblue")) +
  scale_alpha_manual(values = c("Prior" = 0.4, "Posterior" = 0.8)) +
  coord_cartesian(xlim = c(0, 0.5)) +
  labs(
    title = "Call Rate Beta: Prior vs Posterior",
    subtitle = "Gray = data-informed lognormal prior | Blue = posterior learned from acoustic data",
    x = "Beta (calls per whale per hour)",
    y = "Density"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, face = "bold"),
    strip.text = element_text(size = 8, face = "bold"),
    legend.position = "top"
  )

20 Hz: the model is really confident that we have low call rates in the summer. Higher than predicted call rates in winter, spring, and fall!

40 Hz: we basically didn’t learn very much from the data? Other than, we are really confident that call rates are really low, similar to our prior. Which makes sense. This is the noisest, most sparse dataset and I’m not sure it teaches us anything.

6.4 Seasonal Beta Summary Plot

# Create summary with CIs
beta_summary_plot <- beta_long %>%
  group_by(Season, call_type) %>%
  summarise(
    median_beta = median(Value),
    lower_95    = quantile(Value, 0.025),
    upper_95    = quantile(Value, 0.975),
    .groups = "drop"
  ) %>%
  mutate(Season = factor(Season, levels = c("Winter", "Spring", "Summer", "Fall")),
         call_type = factor(call_type, levels = c("20Hz", "40Hz")))

ggplot(beta_summary_plot, aes(x = Season, y = median_beta, 
                               color = call_type, group = call_type)) +
  geom_ribbon(aes(ymin = lower_95, ymax = upper_95, fill = call_type), 
              alpha = 0.2, color = NA) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 4) +
  scale_color_manual(values = c("20Hz" = "steelblue", "40Hz" = "coral")) +
  scale_fill_manual(values  = c("20Hz" = "steelblue", "40Hz" = "coral")) +
  labs(
    title    = "Seasonal Fin Whale Call Rate (Beta)",
    subtitle = "Posterior median ± 95% CI",
    x = "Season", y = "Beta (calls per whale per hour)",
    color = "Call Type", fill = "Call Type"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, face = "bold"),
    legend.position = "top"
  )

This is really cool.

Higher per whale call rates in the winter and spring and fall! Highest in the fall. This is the average in the whole region. We don’t have enough data for this to vary spatially, unfortunately. Fall average: 0.4 20 Hz calls per whale per hour. In the high density pockets, it’s probably more like 120-150 calls per whale per hr. (based on, cue rate: 2-3 calls per min/ICI: 20-30 sec). Also, that has changed over time…

---

7. Seasonal Density Estimates

ggplot(seasonal_dens, aes(x = Season, y = mean_density, group = 1)) +
  geom_ribbon(aes(ymin = lower_95, ymax = upper_95), 
              alpha = 0.2, fill = "steelblue") +
  geom_line(color = "steelblue", linewidth = 1.2) +
  geom_point(color = "steelblue", size = 4) +
  geom_hline(yintercept = 0.00273, color = "red", 
             linetype = "dashed", linewidth = 0.8) +
  annotate("text", x = 4.4, y = 0.00273, 
           label = "Campbell et al. 2015",
           color = "red", size = 3, vjust = -0.5) +
  labs(
    title = "Seasonal Fin Whale Density",
    subtitle = "Mean ± 95% CI averaged over prediction grid | Campbell et al. reference",
    x = "Season",
    y = "Density (whales/km²)"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(size = 14, face = "bold"))

print(seasonal_dens, digits = 5)

##   Season mean_density   lower_95  upper_95
## 1 Winter    0.0013832 0.00076595 0.0027027
## 2 Spring    0.0047028 0.00095828 0.0039418
## 3 Summer    0.0086267 0.00259432 0.0346409
## 4   Fall    0.0041166 0.00158867 0.0136715

matches what we’d expect. Nice to add the seasonal context. Campbell et al. and Becker et al. don’t add that. ## 7.2 Density Heatmaps

# Load density surface
density_surface <- readr::read_csv(paste0(fig_dir, "Bp_density_gp_groupsize_surface.csv"),
                                   show_col_types = FALSE) %>%
  mutate(Season = factor(Season, levels = c("Winter", "Spring", "Summer", "Fall")))

# Load coastline
coastline <- map_data("world", region = c("USA", "Mexico"))

# Posterior median density surface
ggplot() +
  geom_tile(data = density_surface,
            aes(x = lon, y = lat, fill = median_density)) +
  geom_polygon(data = coastline,
               aes(x = long, y = lat, group = group),
               fill = "gray35", color = "gray55", linewidth = 0.3) +
  facet_wrap(~Season, ncol = 2) +
  scale_fill_viridis_c(
    option   = "plasma",
    name     = "Density\n(whales/km²)",
    trans    = "sqrt",
    na.value = "gray20"
  ) +
  coord_fixed(ratio = 1, xlim = c(-125, -116), ylim = c(29.5, 36)) +
  labs(title = "Fin Whale Density - GP Smooth Surface",
       subtitle = "Posterior median | sqrt color scale | Separate GP per season",
       x = "Longitude", y = "Latitude") +
  theme_dark() +
  theme(plot.title   = element_text(size = 14, face = "bold", color = "white"),
        plot.subtitle = element_text(size = 9, color = "gray70"),
        strip.text   = element_text(size = 11, face = "bold"),
        panel.grid   = element_line(color = "gray30", linewidth = 0.2))

7.3 Uncertainty Maps

# 95% CI width
density_surface %>%
  mutate(ci_width = upper_95 - lower_95) %>%
  ggplot() +
  geom_tile(aes(x = lon, y = lat, fill = ci_width)) +
  geom_polygon(data = coastline,
               aes(x = long, y = lat, group = group),
               fill = "gray35", color = "gray55", linewidth = 0.3) +
  facet_wrap(~Season, ncol = 2) +
  scale_fill_viridis_c(option = "magma", name = "95% CI\nwidth",
                       na.value = "gray20") +
  coord_fixed(ratio = 1, xlim = c(-125, -116), ylim = c(29.5, 36)) +
  labs(title = "Density Uncertainty - GP Model",
       subtitle = "Width of posterior 95% credible interval",
       x = "Longitude", y = "Latitude") +
  theme_dark() +
  theme(plot.title = element_text(size = 14, face = "bold", color = "white"),
        strip.text = element_text(size = 11, face = "bold"))

---

8. Overdispersion Parameters

phi_data <- data.frame(
  parameter = c("phi_v (visual)", "phi_20 (20Hz)", "phi_40 (40Hz)"),
  mean = c(mean(posterior_samples$phi_v),
           mean(posterior_samples$phi_20),
           mean(posterior_samples$phi_40)),
  median = c(median(posterior_samples$phi_v),
             median(posterior_samples$phi_20),
             median(posterior_samples$phi_40)),
  lower_95 = c(quantile(posterior_samples$phi_v, 0.025),
               quantile(posterior_samples$phi_20, 0.025),
               quantile(posterior_samples$phi_40, 0.025)),
  upper_95 = c(quantile(posterior_samples$phi_v, 0.975),
               quantile(posterior_samples$phi_20, 0.975),
               quantile(posterior_samples$phi_40, 0.975))
)

ggplot(phi_data, aes(x = parameter, y = mean)) +
  geom_col(fill = "steelblue", alpha = 0.7) +
  geom_errorbar(aes(ymin = lower_95, ymax = upper_95), width = 0.2) +
  geom_text(aes(label = round(mean, 2)), vjust = -0.5) +
  labs(
    title = "Overdispersion Parameters (phi)",
    subtitle = "NegBinomial dispersion for visual and acoustic likelihoods",
    x = NULL,
    y = "phi value"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(size = 14, face = "bold"))

print(phi_data, digits = 3)

##        parameter   mean median lower_95 upper_95
## 1 phi_v (visual) 0.1069 0.1050   0.0774   0.1451
## 2  phi_20 (20Hz) 0.0313 0.0312   0.0278   0.0350
## 3  phi_40 (40Hz) 0.0187 0.0185   0.0153   0.0224

Negative Binomial distribution has two parameters: - μ (mean) = what the model predicts - φ (phi) = overdispersion parameter

In my model: obs_count ~ NegBinomial(predicted_count, phi)

HIGHER phi = LESS overdispersion = data fits model well LOWER phi = MORE overdispersion = lots of extra noise

So visual data fits better than acoustics, 20 Hz better than 40 Hz.

Visual: NegBinomial variance = μ + μ²/0.11

If model predicts 1 group: variance = 1 + 1/0.11 = 1 + 9 = 10 SD = sqrt(10)= 3.2 groups

There’s unexplained variability in how many groups we see per subsegment beyond what the spatial GP predicts. Might get better if we add environmental covariates?

Acoustics: more overdispersion than visual

If model predicts 5 calls: variance = 5 + 25/0.03 = 5 + 833 = 838 SD = sqrt(838) = 29 calls

Whale 20 Hz calling behavior is EXTREMELY variable.

40 Hz behavior is even more variable. What we expect!

---

Session Info

sessionInfo()

## R version 4.4.3 (2025-02-28)
## Platform: x86_64-apple-darwin20
## Running under: macOS Sequoia 15.7.4
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/Los_Angeles
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] gridExtra_2.3       posterior_1.6.1     bayesplot_1.15.0   
## [4] stringr_1.5.1       ggplot2_3.5.1       tidyr_1.3.1        
## [7] dplyr_1.1.4         rstan_2.32.7        StanHeaders_2.32.10
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6         tensorA_0.36.2.1     xfun_0.56           
##  [4] bslib_0.9.0          QuickJSR_1.6.0       inline_0.3.21       
##  [7] lattice_0.22-6       tzdb_0.5.0           vctrs_0.6.5         
## [10] tools_4.4.3          generics_0.1.3       stats4_4.4.3        
## [13] parallel_4.4.3       tibble_3.2.1         pkgconfig_2.0.3     
## [16] Matrix_1.7-3         checkmate_2.3.2      distributional_0.5.0
## [19] RcppParallel_5.1.10  lifecycle_1.0.4      compiler_4.4.3      
## [22] farver_2.1.2         munsell_0.5.1        codetools_0.2-20    
## [25] htmltools_0.5.8.1    maps_3.4.3           sass_0.4.9          
## [28] yaml_2.3.10          pillar_1.10.1        crayon_1.5.3        
## [31] jquerylib_0.1.4      cachem_1.1.0         abind_1.4-8         
## [34] tidyselect_1.2.1     digest_0.6.37        stringi_1.8.4       
## [37] reshape2_1.4.5       purrr_1.0.4          labeling_0.4.3      
## [40] fastmap_1.2.0        grid_4.4.3           colorspace_2.1-1    
## [43] cli_3.6.5            magrittr_2.0.3       loo_2.8.0           
## [46] pkgbuild_1.4.8       utf8_1.2.4           readr_2.1.5         
## [49] withr_3.0.2          scales_1.3.0         backports_1.5.0     
## [52] bit64_4.6.0-1        rmarkdown_2.29       matrixStats_1.5.0   
## [55] bit_4.6.0            hms_1.1.3            evaluate_1.0.3      
## [58] knitr_1.50           viridisLite_0.4.2    rlang_1.1.6         
## [61] Rcpp_1.0.14          glue_1.8.0           rstudioapi_0.17.1   
## [64] vroom_1.6.5          jsonlite_1.9.1       R6_2.6.1            
## [67] plyr_1.8.9