GLM / GLMM with biological interpretation

1. Odds/rate ratios for imprecise contrasts are biased upward in the mean. exp(beta-hat) is median-unbiased but not mean-unbiased: E[exp(b)] is about exp(beta + SE^2/2). When the SE of a log ratio exceeds about 0.35 (factor contrasts in Bernoulli GLMMs; small samples), the ratio estimand fails the mean-bias gate: +13% (binomial_glmm habitat), +15-90% in the small designs. The link-scale coefficient passes and the ratio's 95% CI keeps nominal coverage. Report ratios with their CIs; do not average point ratios across studies.
2. Few groups. With 10 sites x 4 visits (NB GLMM) Wald intervals for site-level effects under-cover (0.875-0.88) and the site SD is biased -42%. With 15 Bernoulli sites x 4 the site SD interval under-covers (0.845). The module flags fewer than 5 levels; between 5 and ~15 groups, treat site-level effects' intervals as somewhat optimistic and run a planted-truth check at the researcher's design. Power below ~20 sites is low for site-level effects in any case.
3. NB size theta needs data. With 100 observations and 30% structural zeros, theta is biased +29% (coverage 0.87). Overdispersion is then hard to separate from zero inflation; the coefficients and zero-inflation probability are still recovered.
4. Wald intervals. All intervals are Wald (log scale for SDs and dispersion). Profile or bootstrap intervals are not implemented. For random-effect SDs near zero the interval is meaningless and the report says the fit is singular.
5. mean_response_at_reference is for a typical group (random effect = 0), not the population average.
6. Not validated: random slopes (rejected by the parser), non-default links (probit, cloglog, inverse; accepted but their interpretation sentences are generic), zero-inflation formulas other than ~1 (supported, zi_beta_* reported, not recovery-tested), polynomial terms (supported in fitting and plots, not recovery-tested), the stats::glm fallback engine (cross-checked to agree, not separately recovery-tested), and DHARMa test calibration (the tests come from DHARMa, unmodified).
7. Spatial/temporal autocorrelation is not modelled or tested for.

Date: 2026-10-18. Code state: dev/glm/code_md5_final_suite.txt (all runs below). Later changes touch only plotting code in output.R and the engine description string in module.R (residual-panel subtitles, and a module-local grouped conditional-effects plot for interactions); recovery never calls it. Final hashes: dev/glm/code_md5_after_plot_tweak.txt. 2026-10-19: figures.R (publication figure suite, glm_figures) added and registered in module.R; it only reads the fit (no change to simulate/fit/estimates), so the results below still apply. Dev tests: dev/glm/figs/test_figs.R. Engine: glmmTMB 1.1.15.2, R 4.5.3. Gate: recovery_gate() with its fixed thresholds (|bias| <= max(10% of truth, 3 MC SE); coverage >= 0.95 - 3 MC SE, here 0.904 for 200 reps and 0.921 for 500; >= 90% converged). Recovery uses fast_settings (identical to the real-data settings except no DHARMa/R2), REML for mixed models.

Reproduce: bash dev/glm/run_validation.sh (configs in modules/glm/validation/; outputs in results/glm/validation/<run>/), bin/Rscript dev/glm/crosscheck.R, bin/Rscript dev/glm/nb_ll_check.R. Tables below come from dev/glm/validation_tables.R. Each fit takes 0.3-0.9 s, so all 3,700 recovery fits took about 25 min on 16 cores.

NB2 GLMM, 30 sites x 8 visits, site-level forest cover (normal 55 +/- 20) and habitat (3 levels), visit-level temperature, effort offset U(0.5, 2), theta = 2, site SD = 0.5. 200/200 converged.

(beta_habitat_wetland: 10.9% relative bias but |bias| 0.033 < 3 MC SE = 0.057, so it is within Monte Carlo noise.)

Full tables are in the appendix. Designs: 30 groups x 8 obs unless noted.

False-positive rate (effect = 0), nsim = 500 (null_effects: forest and temp effects 0, default design otherwise):

Both are consistent with 0.05. For comparison, with ML instead of REML the site-level rate was 0.084 (coverage 0.916; results/glm/validation_ml/null_effects), which is why REML is the default for mixed models.

Mis-specified fits must fail the gate:

dev/glm/crosscheck.R: 20 simulated datasets per case. The module's estimates are compared with lme4 (lmer, glmer, glmer.nb), stats::glm, MASS::glm.nb, and, where no package was available, hand-written likelihoods maximised with optim (beta, ZIP and ZINB GLMs; betareg/pscl are not installed). GLMM comparisons use the module with REML: false because glmer is ML only. Output: results/glm/validation/crosscheck_summary.csv.

Explained differences:

• SE differences of about 6% for glm.nb and glm(Gamma) arise because glm uses expected information and, for Gamma, the Pearson dispersion. glmmTMB uses the observed Hessian at the ML shape.
• NB2 GLMM vs glmer.nb: coefficients differ by up to 0.14 SE (median 0.11 SE). In dev/glm/nb_ll_check.R the glmmTMB objective at the glmmTMB MLE is higher than at glmer.nb's estimates in all 20 datasets, by 0.006 to 0.024 log-likelihood units (results/glm/validation/nb_glmm_loglik_check.csv). The module therefore reaches the better optimum. glmer.nb's own reported log-likelihood differs from glmmTMB's by -0.06 to +0.12, so the two packages evaluate this likelihood differently; glmer.nb fits theta in an outer loop. The difference is small next to the SEs.

core/run.R power glm --vary n_groups=10,20,30,50,80 --nsim 100 (default design, 8 visits per site; results/glm/validation/power_ngroups/).

About 40 sites give 80% power for a site-level effect of 0.3 per SD (rate ratio 1.35); a 3-level habitat contrast of 0.4 needs more than 80 sites. At 10 sites coverage is below nominal (0.88), consistent with hard_nb_small.

bin/Rscript core/run.R fit glm --config modules/glm/example/config.yaml --out results/glm/example (plus config_binomial.yaml, config_beta.yaml). The examples were simulated with known values (dev/glm/make_examples.R, dev/glm/example_truth.txt). All three converged, DHARMa tests all p > 0.4, and every generating value is inside its 95% CI:

• Wood Thrush NB2 GLMM: canopy rate ratio per 10 points 1.38 (1.17-1.62), truth exp(0.25) = 1.28; conifer vs deciduous 0.50 (0.26-0.96), truth 0.41; elevation per 100 m 0.95 (0.86-1.06), truth 0.89; station SD 0.56 (0.34-0.93), truth 0.45.
• Hatching binomial GLMM: OR per 1 m nest height 1.28 (1.09-1.50), truth 1.28; per 50 m from edge 1.65 (1.30-2.09), truth 1.49; predator guard 1.87 (1.17-3.0), truth 2.01.
• Seagrass beta GLMM: OR per 1 m depth 0.645 (0.576-0.721), truth 0.638; per 5 NTU 0.727 (0.60-0.88), truth 0.741; phi 12.1 (8.8-16.6), truth 12.

I read the HTML reports. The interpretation sentences match the estimates table: rate ratios, percentage changes, response-scale predictions across observed ranges, the random-effect spread and the dispersion statements. Conditional-effect plots render on the response scale with 95% bands and rugs. Factors are drawn as points with intervals, and interactions as several coloured curves with bands (dev/glm/inter_png.R; the core cond_effect_plot() loses its grouping under ggplot2 4.0, so the module draws grouped curves itself). The forest plot and DHARMa panels render as well.

An earlier full suite used ML (results/glm/validation_ml/, same configs). Under ML, re_sd_site was biased low in every mixed design: -11.9% default, -11.6% Gaussian, -17.9% Bernoulli, -9.6% beta, -8.3% Gamma, -10.1% ZIP, -17.3% ZINB. Group-level contrasts were under-covered in the ZIP/ZINB GLMMs (0.880-0.890 < 0.904), and the null false-positive rate for the site-level covariate was 0.084. Under REML (started from the ML solution) all of these pass: re_sd within +/-3%, ZI coverage 0.915-0.95, false-positive rate 0.058. With glmmTMB's default start, REML stopped at a spurious zero-SD optimum in 7/200 beta GLMM datasets. Restarting from the ML estimates reached a better REML optimum (lower objective, positive-definite Hessian) in all 7, so the module always starts REML from ML. If REML still fails to converge, the ML fit is reported with a note.

Scripts: dev/glm/verify/ (verify_recovery.R = own simulators + driver, run_verify.sh, rerun_after_fix.sh, usability.sh, messy/). Outputs: results/glm/verify/. Log: logs/glm_verify.log. Final code hashes: dev/glm/verify/code_md5_verified.txt.

Literature / parameterisation check. Engine assumptions were checked against glmmTMB 1.1.15.2 ?sigma.glmmTMB and Brooks et al. (2017): NB2 Var = mu(1 + mu/theta) with betadisp = log(theta); beta Var = mu(1-mu)/(1+phi) (Ferrari & Cribari-Neto 2004); Gamma sigma = 1/sqrt(shape) = CV, so the module's log CV = -0.5 x betadisp is right; zero inflation is a logit-link mixture with E[y] = (1-pi) mu. The module's simulator shares canon_family, link_inv and coef_label with the fitter and standardises with the sample SD exactly as the fitter does, so I did not reuse it.

Independent planted truth. My simulators (no module code) set truth on the RAW covariate scale and derive the per-SD estimands analytically (beta_SD = beta_raw x sample SD; intercept at the sample mean; for x:factor interactions, factor contrast = b2 + b3 x mean(x)). The per-unit estimands (*_per_<u>) are therefore checked against raw-scale truth that does not depend on the scaling code. Fitted with module$fit() + fast_settings; summarised with core summarise_recovery() / recovery_gate() (thresholds unchanged). All replicates converged.

Controls (my designs).

• False-positive rate, Gamma design with both effects 0, nsim 500: obs-level depth 0.048, pond-level treatment 0.066 (MC SE 0.010). All 9 estimands PASS.
• ZIP data fitted as plain Poisson: FAIL as required (beta_intercept -43%, coverage 0.000; slopes coverage 0.745-0.815).
• Gamma GLMM fitted without the pond random effect: FAIL as required (pond-level treatment coverage 0.740, intercept 0.680, gamma_cv +15.5% / coverage 0.115). The obs-level depth effect still passes.

Builder's numbers reproduced. Rerunning recover for default (200), binomial_glmm (200) and null_effects (500) with the same seeds gave the tables in sections 1, 2 and 4 exactly (default PASS 12/12; binomial_glmm odds ratio FAILs +13.1% / +13.3%; FPR 0.058 / 0.042). After my engine fix (below), default, poisson_glm and gamma_glmm summaries are byte-identical to the builder's (results/glm/verify/after_fix/).

Problems found and fixed (module files only; core untouched).

1. Offset partial matching (engine.R, module.R). sp$offset used R's partial $ matching: if the spec had no offset element but had offset_per, the value of offset_per (e.g. 1 or 10) was taken as the offset column, so log(df[[1]]) (often the response) entered as an offset. In my first smoke test this made every Gamma GLMM fit collapse (slopes 0, CV 1e-7) and every binomial-trials fit fail. Through read_data() the same happened when a config had offset_per: without offset: (stopped with a confusing "column(s) 10 not found"). Now exact [["offset"]] lookups, plus a clear check_data message for offset_per without offset. The builder's suite was not affected (its simulator always sets offset).
2. Misspelled group ids (module.R). WT03 vs WT-03 was silently fitted as a 46th station. check_data now stops with "grouping column station has labels that look like the same station spelled differently: 'WT-03' (3 rows) vs 'WT03' (1 rows) ..." (case, punctuation and leading zeros ignored). fit notes also flag single-observation groups when the typical group has >= 3.
3. Figure export (output.R). The figure named "Effect of canopy cover (%)" failed to save in core export_figures() (% in the file name; ggsave treats it as a format string). The figure list names now spell out % and /. The underlying core issue (file stem = figure name, unsanitised) remains for other modules.
4. Interaction sentences (output.R). They read "Each 1-SD increase (or, for a factor, the named level instead of the reference) in forest_type_conifer multiplies the rate ratio for canopy_cover_pct by 0.922", with raw column names. Now: "when forest type = conifer rather than deciduous, the rate ratio per 1 SD of canopy cover (%) (15.6) is multiplied by 0.922 (95% CI 0.49 to 1.74)". Identity-link and continuous x continuous versions were checked too.
5. Detection caveat (output.R). For count and 0/1 families the interpretation now says that effects on recorded counts/detections combine abundance/occurrence with detectability, and points to occupancy / N-mixture modules. Not added for successes/trials, beta, Gamma or gaussian.
6. Missing-value message now says "data rows N, not counting the header line".

Usability. On the Wood Thrush example every sentence was checked by hand against the estimates table (for example RR per 10 points = exp(0.500 x 10/15.6) = 1.378; station spread exp(+/-1.96 x 0.56) = 0.334-3.0; NB variance/mean at the mean count 1 + 1.13/5.89 = 1.19; conifer 1.156 x 0.498 = 0.576). The canopy curve (0.47 to 2.84) and forest-type points (1.16 / 0.576 / 1.32) match. Wording is associational and states that the response is a count per 10 minutes. The NA, wrong-column-name (lists the columns present, explains canopy cover pct became canopy.cover.pct), typo and offset_per messages are actionable. Binomial and beta example sentences were checked the same way.

Not fixed / still open. The module has no figures(fit, data, opts) per the newer MODULE_SPEC / FIGURES.md; core falls back to plots(). Core's journal styling moves the interaction legend to the right, which squeezes the panel. Core export_figures() does not sanitise file names. The known limits above (ratio retransformation bias, few groups, theta with 100 obs) are confirmed, not changed.

default: 200 replicates, 200 converged (100.0%), 0 errors, coverage bound 0.904; mean fit time 0.87 s

gaussian_lmm: 200 replicates, 200 converged (100.0%), 0 errors, coverage bound 0.904; mean fit time 0.53 s

binomial_glmm: 200 replicates, 200 converged (100.0%), 0 errors, coverage bound 0.904; mean fit time 0.57 s

binomial_trials_glm: 200 replicates, 200 converged (100.0%), 0 errors, coverage bound 0.904; mean fit time 0.32 s

poisson_glm: 200 replicates, 200 converged (100.0%), 0 errors, coverage bound 0.904; mean fit time 0.34 s

beta_glmm: 200 replicates, 200 converged (100.0%), 0 errors, coverage bound 0.904; mean fit time 0.91 s

gamma_glmm: 200 replicates, 200 converged (100.0%), 0 errors, coverage bound 0.904; mean fit time 0.65 s

zip_glmm: 200 replicates, 200 converged (100.0%), 0 errors, coverage bound 0.904; mean fit time 0.61 s

zinb_glmm: 200 replicates, 200 converged (100.0%), 0 errors, coverage bound 0.904; mean fit time 0.88 s

interaction_glmm: 200 replicates, 200 converged (100.0%), 0 errors, coverage bound 0.904; mean fit time 0.62 s

hard_nb_small: 200 replicates, 200 converged (100.0%), 0 errors, coverage bound 0.904; mean fit time 0.56 s

hard_binomial_small: 200 replicates, 200 converged (100.0%), 0 errors, coverage bound 0.904; mean fit time 0.51 s

hard_zinb_small: 200 replicates, 200 converged (100.0%), 0 errors, coverage bound 0.904; mean fit time 0.36 s

null_effects: 500 replicates, 500 converged (100.0%), 0 errors, coverage bound 0.921; mean fit time 0.67 s

misspec_poisson: 200 replicates, 200 converged (100.0%), 0 errors, coverage bound 0.904; mean fit time 0.55 s

misspec_norandom: 200 replicates, 200 converged (100.0%), 0 errors, coverage bound 0.904; mean fit time 0.35 s