The Finding

GAM_size_structure_Clean.R’s selectivity models build dat_perch/dat_roach with one row per observed catch event, plus one row per zero-catch net with length_cm = NA. Because s(length_cm, ...) needs a length value, mgcv::gam()’s default na.action = na.omit silently drops every zero-catch net — those models compared catch counts only among mesh/net-type combinations where something was caught, never against true absence.

GAM_size_structure_ZeroFilled.R fixes this by expanding every net × mesh × net-type combination onto the 8 Size_Class categories already used elsewhere in this analysis, with catch_n as a real, explicit 0 where nothing of that size was caught. Refitting the same models on this zero-filled data reverses the headline conclusion:

tribble(
  ~Species, ~`Zone scope`,   ~`p, catch-only (original)`, ~`p, zero-filled (this analysis)`,
  "Perch",  "Both zones",    "1.8 × 10⁻¹³",                "0.929",
  "Roach",  "Both zones",    "< 2.2 × 10⁻¹⁶",               "0.464",
  "Perch",  "Benthic only",  "5.6 × 10⁻¹²",                "0.718",
  "Roach",  "Benthic only",  "< 2.2 × 10⁻¹⁶",               "0.905"
) |>
  knitr::kable(
    caption = paste(
      "Likelihood-ratio test, mesh×net-type selectivity smooth (Model 1) vs shared-smooth",
      "baseline (Model 0). 'catch-only' p-values are from the original report",
      "(Gillnet_Selectivity_BothZones_vs_Benthic.Rmd); every one of the four scopes flips",
      "from highly significant to non-significant once zero-catch nets are counted."
    )
  )
Likelihood-ratio test, mesh×net-type selectivity smooth (Model 1) vs shared-smooth baseline (Model 0). ‘catch-only’ p-values are from the original report (Gillnet_Selectivity_BothZones_vs_Benthic.Rmd); every one of the four scopes flips from highly significant to non-significant once zero-catch nets are counted.
Species Zone scope p, catch-only (original) p, zero-filled (this analysis)
Perch Both zones 1.8 × 10⁻¹³ 0.929
Roach Both zones < 2.2 × 10⁻¹⁶ 0.464
Perch Benthic only 5.6 × 10⁻¹² 0.718
Roach Benthic only < 2.2 × 10⁻¹⁶ 0.905

Caveat carried through this whole report: this could be a resolution artefact, not (only) a real reversal. A follow-up check refit all four zero-filled models with finer, fixed 5cm-wide bins instead of the 8 Size_Class categories. Two of the four moved back toward significance (Roach both-zones: p = 0.019; Perch Benthic: p = 0.075) while the other two did not (Perch both zones: p = 0.518; Roach Benthic: p = 0.698) — still nowhere near the original’s p < 10⁻¹¹ range for any of the four. The coarser the bins, the less power the model has to detect a real shape difference, so the exact non-significance reported here is somewhat resolution-dependent — treat “no shape difference detected” as the qualitative takeaway, not the precise p-values.

Methods (brief)

  • Zero-fill construction. Every net × mesh size × net-type combination is crossed with all 8 Size_Class categories (0-5, 5-10, ..., 40-50, >50 cm), each represented by its bin midpoint; catches are aggregated into (net, mesh, net-type, size class) counts, 0 where nothing was caught.
  • Why 8 coarse bins, not something finer. A first attempt used a 1cm length grid (~103k rows) fit with mgcv::gam(method = "REML"); a single interaction model ran over 2 hours without converging. Coarsening to Size_Class (~32k rows) and switching to mgcv::bam(method = "fREML", discrete = TRUE) fits the same formula in seconds — but at the cost of the resolution caveat above.
  • GAM structure. catch_n ~ offset(log(area_m2)) + s(length_cm, by = mesh_net, k = 8) + mesh_mm + net_type [+ zone] + s(lake_year, bs = "re"), family = mgcv::nb(). zone is included only in the both-zones models; k is capped at 8 (only 8 unique length positions per mesh).
  • A parametric aliasing issue. The standalone mesh_mm/net_type main effects are only weakly identified in Model 1 (very large SEs) — aliased with the per-level baseline that s(length_cm, by = mesh_net) already absorbs, much more visible under this sparse 8-position binning than with continuous length. The tables below show only net_type/zone rows for this reason; the anova(m0, m1) comparison is the reliable read, not individual Model 1 coefficients.
  • Curve color and range fix. Each net type’s curve is drawn only across the lengths it actually caught something at (previously the range was pooled across net types, which let a single stray catch from one net type stretch the other’s curve into totally unsupported territory — visible as an extreme confidence-interval blowup, since fixed).

Perch

Model comparison and key coefficients

anova(m0_Perch, m1_Perch, test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: catch_n ~ offset(log(area_m2)) + s(length_cm, by = mesh_mm, k = 8) + 
##     mesh_mm + net_type + zone + s(lake_year, bs = "re")
## Model 2: catch_n ~ offset(log(area_m2)) + s(length_cm, by = mesh_net, 
##     k = 8) + mesh_mm + net_type + zone + s(lake_year, bs = "re")
##   Resid. Df Resid. Dev     Df Deviance Pr(>Chi)
## 1     32089      14896                         
## 2     32040      14861 48.215   34.646    0.929
bind_rows(
  `Model 0 (shared smooth)` = extract_coefs(m0_Perch),
  `Model 1 (interaction)` = extract_coefs(m1_Perch),
  .id = "Model"
) |> knitr::kable(digits = 3, caption = "Perch, both zones")
Perch, both zones
Model term Estimate Std. Error Pr(>|t|)
net_typeMOD…1 Model 0 (shared smooth) net_typeMOD -0.101 0.068 0.137
zonePelagic…2 Model 0 (shared smooth) zonePelagic -2.984 0.095 0.000
net_typeMOD…3 Model 1 (interaction) net_typeMOD -0.859 0.642 0.181
zonePelagic…4 Model 1 (interaction) zonePelagic -2.977 0.095 0.000 
anova(m0_Perch_Benthic, m1_Perch_Benthic, test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: catch_n ~ offset(log(area_m2)) + s(length_cm, by = mesh_mm, k = 8) + 
##     mesh_mm + net_type + s(lake_year, bs = "re")
## Model 2: catch_n ~ offset(log(area_m2)) + s(length_cm, by = mesh_net, 
##     k = 8) + mesh_mm + net_type + s(lake_year, bs = "re")
##   Resid. Df Resid. Dev     Df Deviance Pr(>Chi)
## 1     23210      12997                         
## 2     23162      12954 48.216   42.151    0.718
bind_rows(
  `Model 0 (shared smooth)` = extract_coefs(m0_Perch_Benthic),
  `Model 1 (interaction)` = extract_coefs(m1_Perch_Benthic),
  .id = "Model"
) |> knitr::kable(digits = 3, caption = "Perch, Benthic zone only")
Perch, Benthic zone only
Model term Estimate Std. Error Pr(>|t|)
net_typeMOD…1 Model 0 (shared smooth) net_typeMOD -0.191 0.070 0.007
net_typeMOD…2 Model 1 (interaction) net_typeMOD -1.069 0.692 0.122

Selectivity Curves — Both Zones vs Benthic Only

Both zones

Perch_Selectivity_Plot

Benthic zone

Perch_Benthic_Selectivity_Plot

Roach

Model comparison and key coefficients

anova(m0_roach, m1_roach, test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: catch_n ~ offset(log(area_m2)) + s(length_cm, by = mesh_mm, k = 8) + 
##     mesh_mm + net_type + zone + s(lake_year, bs = "re")
## Model 2: catch_n ~ offset(log(area_m2)) + s(length_cm, by = mesh_net, 
##     k = 8) + mesh_mm + net_type + zone + s(lake_year, bs = "re")
##   Resid. Df Resid. Dev     Df Deviance Pr(>Chi)
## 1     32093      10494                         
## 2     32051      10452 42.554   42.731   0.4636
bind_rows(
  `Model 0 (shared smooth)` = extract_coefs(m0_roach),
  `Model 1 (interaction)` = extract_coefs(m1_roach),
  .id = "Model"
) |> knitr::kable(digits = 3, caption = "Roach, both zones")
Roach, both zones
Model term Estimate Std. Error Pr(>|t|)
net_typeMOD…1 Model 0 (shared smooth) net_typeMOD -0.109 0.074 0.141
zonePelagic…2 Model 0 (shared smooth) zonePelagic -1.215 0.080 0.000
net_typeMOD…3 Model 1 (interaction) net_typeMOD -0.674 1.332 0.613
zonePelagic…4 Model 1 (interaction) zonePelagic -1.223 0.080 0.000 
anova(m0_roach_Benthic, m1_roach_Benthic, test = "Chisq")
## Analysis of Deviance Table
## 
## Model 1: catch_n ~ offset(log(area_m2)) + s(length_cm, by = mesh_mm, k = 8) + 
##     mesh_mm + net_type + s(lake_year, bs = "re")
## Model 2: catch_n ~ offset(log(area_m2)) + s(length_cm, by = mesh_net, 
##     k = 8) + mesh_mm + net_type + s(lake_year, bs = "re")
##   Resid. Df Resid. Dev     Df Deviance Pr(>Chi)
## 1     23218     7509.7                         
## 2     23178     7480.9 39.921   28.762   0.9053
bind_rows(
  `Model 0 (shared smooth)` = extract_coefs(m0_roach_Benthic),
  `Model 1 (interaction)` = extract_coefs(m1_roach_Benthic),
  .id = "Model"
) |> knitr::kable(digits = 3, caption = "Roach, Benthic zone only")
Roach, Benthic zone only
Model term Estimate Std. Error Pr(>|t|)
net_typeMOD…1 Model 0 (shared smooth) net_typeMOD -0.091 0.083 0.276
net_typeMOD…2 Model 1 (interaction) net_typeMOD 1.188 1.796 0.508

Selectivity Curves — Both Zones vs Benthic Only

Both zones

Roach_Selectivity_Plot

Benthic zone

Roach_Benthic_Selectivity_Plot

Conclusions

  • No statistically supported CEN vs MOD selectivity-shape difference, for either species, in either zone scope, once zero-catch nets are counted as real zero observations (all four p-values in the finding table above). The original catch-only analysis’s strong interaction signal for both species appears to have been substantially inflated by silently excluding zero-catch nets.
  • This could still be partly a resolution artefact — see the caveat above. The 5cm-bin follow-up pulled two of the four comparisons back toward significance, so don’t read the exact p-values as final; read “no shape difference detected at this resolution” as the finding.
  • Zone remains the dominant, robustly-estimated predictor for both species — zonePelagic ≈ −2.98 (Perch) / −1.22 (Roach) in Model 1, both with tight SEs (≈ 0.08–0.09) and p < 2×10⁻¹⁶, unlike any net_type estimate.
  • One notable, opposite-signed exception: Perch, Benthic-only, Model 0 (no interaction): net_typeMOD = −0.191, SE = 0.070, p = 0.0067 — MOD catching significantly fewer Perch than CEN, where the original catch-only analysis found MOD catching more. This appears in only one of eight fitted models here and reverses sign from the original result, so it should be treated as a hypothesis to test further (e.g. at finer resolution), not a settled finding.
  • net_type main effects are otherwise small and non-significant throughout (Model 0 estimates ranging p = 0.14–0.28 elsewhere), for both species, in both zone scopes.

Session Information

sessionInfo()
## R version 4.6.1 (2026-06-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 8 x64 (build 9200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=Portuguese_Portugal.utf8  LC_CTYPE=Portuguese_Portugal.utf8   
## [3] LC_MONETARY=Portuguese_Portugal.utf8 LC_NUMERIC=C                        
## [5] LC_TIME=Portuguese_Portugal.utf8    
## 
## time zone: Europe/Berlin
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] gratia_0.11.2   mgcv_1.9-4      nlme_3.1-169    lubridate_1.9.5
##  [5] forcats_1.0.1   stringr_1.6.0   dplyr_1.2.1     purrr_1.2.2    
##  [9] readr_2.2.0     tidyr_1.3.2     tibble_3.3.1    ggplot2_4.0.3  
## [13] tidyverse_2.0.0
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.10        generics_0.1.4     stringi_1.8.7      lattice_0.22-9    
##  [5] hms_1.1.4          digest_0.6.39      magrittr_2.0.5     evaluate_1.0.5    
##  [9] grid_4.6.1         timechange_0.4.0   RColorBrewer_1.1-3 fastmap_1.2.0     
## [13] jsonlite_2.0.0     Matrix_1.7-5       ggokabeito_0.1.0   scales_1.4.0      
## [17] jquerylib_0.1.4    cli_3.6.6          rlang_1.3.0        splines_4.6.1     
## [21] withr_3.0.3        cachem_1.1.0       yaml_2.3.12        otel_0.2.0        
## [25] tools_4.6.1        tzdb_0.5.0         nanonext_1.10.1    vctrs_0.7.3       
## [29] R6_2.6.1           lifecycle_1.0.5    tweedie_3.1.0      pkgconfig_2.0.3   
## [33] pillar_1.11.1      bslib_0.11.0       gtable_0.3.6       Rcpp_1.1.2        
## [37] glue_1.8.1         statmod_1.5.2      xfun_0.60          tidyselect_1.2.1  
## [41] knitr_1.51         farver_2.1.2       patchwork_1.3.2    htmltools_0.5.9   
## [45] mirai_2.7.1        labeling_0.4.3     rmarkdown_2.31     compiler_4.6.1    
## [49] S7_0.2.2           mvnfast_0.2.8