Population time series - trends, density dependence and drivers of growth (state-space models)

1. Variance components from one series. sd_process is biased (REML -17% at 25 yr; Bayes +31%) and sd_observation is biased in the Gompertz and Bayes models. Trend and abundance intervals stay valid. Pooling sites (marss) or an independent estimate of counting error (prior_sd_observation) fixes it.

2. Density dependence is not recoverable from a single count series (Knape 2008; Auger-Methe et al. 2016). Point estimates of b are overestimated by 30-140% depending on length and engine. The profile interval for b excludes 0 in 47% of series without density dependence.

   • The bootstrap LRT holds its size (0.04-0.07) but has almost no power at 15 years.
   • In bayes the U(-2, 0) prior means the interval for b never contains 0.
   • The interpretation text states all of this.
3. equilibrium_N (Gompertz ML) is undefined when b is near 0. The mean of the MLEs is dominated by huge values and coverage is 0.889 at 30 yr. The report suppresses K when its interval spans more than a 20-fold range.
4. Trend intervals are conservative. The false-positive rate is 0.025, and coverage is 0.97-1.0 below 20 years. Power is therefore lower than with plug-in intervals, which are invalid (coverage 0.81).
5. process_cov: equalvarcov (MARSS EM) converged in only 75% of replicates. Treat process_corr as exploratory.
6. loglinear with 3 sites x 10 years: 14.5% of fits have a singular Hessian. It also does not separate process from observation error.
7. Constant detection is assumed. Abundance is on the count scale, and a trend in detectability is indistinguishable from a population trend.
8. The importance sampler for ML abundance intervals can have low ESS (<200) for Gompertz fits with a strongly curved likelihood. The fit then adds a note.
9. Bayes DD fits mix slowly when b is poorly identified. In the 60-year cross-check, minimum ESS was 188 at 4 x 40,000 iterations; the report flags runs below ESS 400.

1. REML trend intervals are very conservative when counting error dominates. At 20 years with sd_process 0.06 and sd_observation 0.25, a +4%/yr trend was never detected (power 0.000, false-positive rate 0.000, coverage 1.000), whereas GLS with the true variances would detect it in 60% of series. The intervals are valid but can be about 3x wider than necessary; see "Independent verification". The report now says that "no clear trend" is not evidence of stability and points to a power analysis.
2. loglinear trend intervals are too narrow when the population itself wanders (random walk with drift plus counting error): coverage 0.825 at 4 sites x 25 years. The report warns when the fitted AR(1) correlation exceeds 0.9 and always says to compare with marss/exponential.

Every number below was produced by the scripts named next to it. Rerun the full recovery set with bash dev/timeseries/run_recover.sh <nsim> <cfg...> (configs in dev/timeseries/<cfg>.yaml, results in results/timeseries/recover_<cfg>/). bin/Rscript dev/timeseries/gate_md.R regenerates the tables in the appendix. All recoveries used the module's fast_settings:

• Bayes: 3 chains x 12,000 iterations, 4,000 burn-in, thin 4; the long design used 30,000 / 10,000 / 10.
• ML/REML: 300 trajectory draws, no bootstrap density-dependence test inside each replicate.

The gate is recovery_gate() with its fixed thresholds:

• bias	<= max(10% of truth, 3 MC SE);

• coverage >= 0.95 - 3 MC SE (0.885 at n = 100, 0.904 at n = 200);
• >= 90% of replicates converged.

Default parameters (module$defaults): growth_rate_mu = -0.03 (-2.96%/yr), sd_process = 0.10, sd_observation = 0.15, density_dependence_b = -0.3, equilibrium_N = 500, growth_beta = -0.15 per SD, start N = 500. Final run: 2026-10-20, after the SIR change to ML abundance intervals (below).

Bottom line for researchers:

• Trends, the trajectory of true abundance and covariate effects on growth are recovered with honest intervals at realistic designs.
• Process and observation SDs from one series of 25 years are biased by about 15-40%. Do not over-interpret them.
• Density dependence from one count series is not recoverable without independent information on counting error. Even with it, b is overestimated (bayes_obsprior: -52%).

Single series (exponential growth with drift, REML):

Why REML is the default:

• Plain ML estimates sd_process at exactly 0 in many series, and the trend interval then collapses (coverage 0.675).
• REML integrates the variances out on a grid (flat prior on both SDs), which keeps trend intervals honest.
• The price is conservatism: coverage is 0.97-1.0 below 20 years, and the false-positive rate is 0.025 instead of 0.05 (section 4).

• dev/timeseries/prior_study.R (results/timeseries/prior_study.csv) compared four ways of handling the variances at 25 years (sd_p 0.1, sd_o 0.15). Each pair is coverage / power:

  • plug-in GLS: 0.81 / 0.42 (invalid);
  • log-uniform prior: 0.91 / 0.28;
  • sqrt prior: 0.95 / 0.20;
  • flat prior: 0.968 / 0.165.
• Flat was kept because it is the only option with coverage >= 0.90 in every scenario studied.

sd_process (single series) is the one default-design estimand that fails:

• REML puts it on the boundary (0) in about 14% of 25-year series.
• dev/timeseries/sd_estimator_study.R tried the grid-posterior median (+12%) and mean (+20%). No point estimator came within 10% at 25 years. At 50 years REML bias is -6.6%.
• Interval coverage is fine (0.98). The point estimate is not reliable, and the interpretation text says so when it is at the boundary.

Multi-site (marss, pooled REML, shared trend), nsim 200:

Density dependence, single series (gompertz, exact Kalman ML), nsim 200:

\* equilibrium_N passes the gate at 60 years only formally: the mean of the MLEs is dominated by a few series with b near 0, which makes 3 MC SE enormous. Report K only when its interval spans less than a 20-fold range. The interpretation enforces that rule.

Fix made in this round (ML abundance intervals). Abundance intervals for the ML engines previously drew parameters from the asymptotic normal (Hessian) distribution. That collapsed to zero width whenever sd_observation was estimated at 0, and exploded when b was drawn near -2. The coverage of Gompertz N_final was 0.573. They now use adaptive sampling-importance- resampling from the likelihood (kalman.R kf_sir):

• flat priors, with the SDs flat on the SD scale;
• a multivariate-t proposal, re-centred over 3 pilot rounds;
• typical importance ESS of 1,300-4,400 for exponential and 200-1,300 for Gompertz, out of 20,000 proposals; the fit notes when ESS < 200.

After the change, N_final coverage is 0.920 (Gompertz, 30 yr) and 0.960 (exponential ML, previously 0.79).

Bayesian state-space model (bayes, Nimble), nsim 100, 24 s per replicate:

• Convergence (R-hat < 1.05 for every reported quantity) held in 96-100% of replicates.
• The SD biases come from a half-normal(0, 1) prior meeting the Knape/Dennis trade-off. With an informative observation-error prior (Knape's recommended cure) everything except b passes.
• The realised trend fails only at the DD default design. There the true realised trend is about 0 (a population at equilibrium), so the 10% relative tolerance is meaningless. Its absolute bias is 0.2 percentage points per year against an interval width of 2.8.

Log-linear GLMM (loglinear), nsim 200:

• 5 sites x 20 yr: trend bias -0.1%, coverage 0.938, power 0.933; covariate +1.6%, coverage 0.943. All PASS.
• 3 sites x 10 yr: estimates fine (coverage 0.92-0.94), but only 171/200 fits converged (glmmTMB pdHess false: the AR(1) variance collapses). FAIL (conv).

Ignoring observation error (observation_error: false, the classic diffusion/AR(1) regression on log counts):

Ignoring counting error makes density dependence look 2.6x stronger than the truth. This is the regression-to-the-mean artefact of Shenk et al. (1998), Dennis et al. (2006) and Knape (2008). Other controls that fail as they should:

• plain ML instead of REML (trend coverage 0.675);
• MARSS ML (trend coverage 0.877).

Source: dev/timeseries/ddtest_study.R (logs/timeseries_ddtest.log, logs/timeseries_ddtest30.log).

How to use these results:

• Test density dependence with the bootstrap LRT, never with the interval for b.
• With realistic noise (sd_p 0.1, sd_o 0.15) the test has essentially no power at 15 or even 30 years (rejection rate at b = -0.3 is no higher than at b = 0), as Knape (2008) predicts. A non-significant test therefore says nothing about whether density dependence is absent.

core/run.R power with dev/timeseries/power_exp5.yaml: single series, exponential REML, -5%/yr (growth_rate_mu = -0.0513), sd_process 0.1, sd_observation 0.15, nsim 200 per point. Results are in results/timeseries/power_5pct/. Power = share of 95% trend intervals excluding 0.

Decline size at 20 years (dev/timeseries/power_exp20.yaml, vary growth_rate_mu, nsim 200):

What this means for managers:

• With real year-to-year variation of 10% (sd_process 0.1), a single monitored population declining 5% a year is detected with 80% power only after more than 40 years.
• The drift of a random walk is inherently hard to estimate, and the module's integrated intervals are deliberately conservative.
• Monitoring several sub-populations with a shared trend helps a lot: marss at 4 sites x 25 years has power 0.735 for a 3%/yr decline, against 0.145 for one site.
• Plan with --vary n_sites as well as n_years.

Source: dev/timeseries/crosscheck.R, logs/timeseries_crosscheck3.log.

example/herd_counts.csv (31 winters, 3 not flown) was simulated with dev/timeseries/make_example.R. The true values are b = -0.3, K = 500, sd_p 0.10, sd_o 0.15, N_final 643 and a snow effect of -13.9%/SD.

The bayes fit (example/config.yaml) gives:

• N_final 614 (95% CrI 499-768);
• K 536 (478-599);
• snow effect -15.3% (-22.2 to -7.0);
• sd_p 0.10 and sd_o 0.18.

All of these cover the truth. b = -0.59 (-1.11 to -0.29) misses the true -0.3, which is the documented overestimation of density dependence.

The gompertz fit's bootstrap LRT gives p = 0.585, i.e. no detectable density dependence from this one series.

Everything below was produced by a verifier that did not build the module. Scripts: dev/timeseries/verify/. Results: results/timeseries/verify/<scenario>/ and results/timeseries/verify_rerun_<cfg>/. Logs: logs/timeseries_verify_*.log.

Literature and code check

• GSS/EGSS equations match Dennis et al. (2006) as restated by Dennis & Ponciano (2014): X_t = a + c X_{t-1} + E_t, Y_t = X_t + F_t. The module's c = 1 + b and a = -b log K give the stationary mean log K (= a/(1 - c)). EGSS is c = 1. Unlike Dennis et al. (2006), which start from the stationary distribution, the module estimates the initial state freely (the nonstationary variant; MARSS x0 at t = 0). It is stated in the methods text.
• REML for EGSS: the differenced log counts have mean mu k_i, variance k_i sd_p^2 + 2 sd_o^2 and lag-1 covariance -sd_o^2 (Staples et al. 2004; Humbert et al. 2009). Code (egss_S, egss_negreml) matches.
• Bayes: equilibrium log K = xbar - r0/b follows from setting the expected change r0 + b(x - xbar) to 0. Priors in code match README/methods.
• Knape (2008): b is weakly identified when both variances are unknown, and how badly depends on the parameter values; a known observation variance removes the problem. This matches the module's text and the new parameter-dependence result below.
• Fixed: the exponential (REML) methods text claimed a GLS trend with a t-based interval. The code reports the mean and quantiles of the grid-integrated mixture. The methods text and the reml.R header now describe what the code does.
• The simulators share no code with fit() beyond ts_drop_years, and the independent simulators below share none at all.

Rerun of the builder's recoveries (same configs, nsim as stated)

exp_default (200), marss_default (200), gomp_null (200), bayes_default (100) and exp_noobs (200) reproduced the appendix tables exactly. For example: exp_default trend bias -5.89%, coverage 0.975; marss_default trend coverage 0.955, power 0.735; gomp_null interval for b excludes 0 in 0.472; bayes_default covariate coverage 0.958 and b -62%; exp_noobs sd_process +133% and N_final coverage 0.

Independent planted truth (own simulators, indep_recovery.R; module fit() with fast_settings; core summarise_recovery + recovery_gate)

What this adds to the builder's results:

• Trends, abundance, equilibrium N and covariate effects pass at independent designs too. False-positive rates are 0.059-0.082, except REML single-series trends, which are conservative (0.000 here, 0.025 at the default).
• Density dependence is parameter-dependent, not uniformly unrecoverable. With b = -0.5 and process noise larger than counting error, 40 years of Gompertz ML recovered b without bias (coverage 0.96). The Bayes model (U(-2, 0) prior, 30 yr) still overstated it by 33%. The builder's FAIL for b holds at its own noise ratio (sd_o > sd_p) and should not be read as "never". Without independent evidence about the noise ratio, a researcher cannot tell which regime applies, so the report's caution stays.
• Power of REML single-series trends (limit 10). Oracle comparison (oracle_trend.R, 1000 series at the exp_v design): GLS with the true variances had coverage 0.936, power 0.595 and interval width 0.070 (log scale). Plug-in REML GLS had coverage 0.856 (invalid) and power 0.58. The module had coverage 1.000, power 0.000 and width 0.196. Validity is kept at a high price. The builder's prior study found the same trade-off at the default design. This is a design limit to report, not a bug.
• loglinear under random-walk dynamics (ll_rw_rho.R): the fitted AR(1) rho exceeded 0.9 in 66% of series (median 0.945). Coverage was 0.89 for those and 0.69 for the rest, so the rho warning catches only some cases. The fixed caveat sentence therefore states the 0.825 coverage for every loglinear fit.

Usability (as a non-coding user would meet it)

• fit on example/config.yaml reproduced the reported example: N_final 614 (499-768), K 536 (478-599), snow effect -15.3% (-22.2 to -7.0), b -0.59. Max R-hat 1.006, min ESS 684, PPC p = 0.61 / 0.75. Every number in the interpretation was checked by hand against the estimates (half-life log 0.5 / log(1 - 0.591) = 0.78 yr; compounding of trends; 10th/90th percentile growth).

• Messy data (dev/timeseries/verify/messy/):

  • wrong column name: clear message listing the columns present;
  • "n/a" in counts: clear message;
  • duplicated year: clear message;
  • covariates with a model that cannot use them: clear message;
  • misspelled site name ("Nrth Fork" in 2 of 25 years): previously passed check_data silently (MARSS would have fitted a fifth site or stopped with an unhelpful error).
• Figures: trend, variances, trace, coefficients, growth_snow_depth_cm and the MARSS residual panels were opened. No overlaps or clipping, and units are correct. The growth curve crosses 1 at about 53 cm, matching r0 ~ 0 at the covariate mean. Gompertz has no trend figure because it has no trend estimand.

Changes made by the verifier (text, checks and interpretation only; no change to simulate() or fit() numerics)

1. R/common.R ts_check_data: flags site names that differ only in case/punctuation, or by at most 2 edits when one has fewer than half the other's counted years, plus any site with fewer than 3 counted years. Non-numeric-value messages now give the spreadsheet line as well as the data row.
2. module.R methods text (exponential REML) and R/reml.R header: describe the integrated-mixture trend interval actually computed.

3. module.R interpretation:

   • the half-life sentence is given only when the interval for b excludes 0, and is flagged as likely too short;
   • when the Gompertz profile interval for b reaches the -2 search limit, the report says so;
   • the advice to use prior_sd_observation now says it needs model: bayes (the Gompertz ML fit ignored it);
   • "no clear trend" now adds the compounded range and a pointer to power analysis;
   • loglinear warns when AR(1) rho > 0.9 and states the 0.825 random-walk coverage.
4. Example reports regenerated with these texts (results/timeseries/example*/). No gate threshold was changed. Nothing in core/ was changed.

bayes_default (nsim 100, 0 crashed, 24.0 s/rep; config dev/timeseries/bayes_default.yaml)

bayes_hard (nsim 100, 0 crashed, 22.1 s/rep; config dev/timeseries/bayes_hard.yaml)

bayes_long (nsim 100, 0 crashed, 50.0 s/rep; config dev/timeseries/bayes_long.yaml)

bayes_missing (nsim 100, 0 crashed, 93.8 s/rep; config dev/timeseries/bayes_missing.yaml)

bayes_nodd (nsim 100, 0 crashed, 21.6 s/rep; config dev/timeseries/bayes_nodd.yaml)

bayes_null (nsim 100, 0 crashed, 22.1 s/rep; config dev/timeseries/bayes_null.yaml)

bayes_obsprior (nsim 100, 0 crashed, 24.3 s/rep; config dev/timeseries/bayes_obsprior.yaml)

exp_default (nsim 200, 0 crashed, 1.8 s/rep; config dev/timeseries/exp_default.yaml)

exp_hard (nsim 200, 0 crashed, 1.6 s/rep; config dev/timeseries/exp_hard.yaml)

exp_missing (nsim 200, 0 crashed, 1.7 s/rep; config dev/timeseries/exp_missing.yaml)

exp_ml (nsim 200, 0 crashed, 7.6 s/rep; config dev/timeseries/exp_ml.yaml)

exp_noobs (nsim 200, 0 crashed, 0.4 s/rep; config dev/timeseries/exp_noobs.yaml)

exp_null (nsim 200, 0 crashed, 1.8 s/rep; config dev/timeseries/exp_null.yaml)

gomp_default (nsim 200, 0 crashed, 22.6 s/rep; config dev/timeseries/gomp_default.yaml)

gomp_hard (nsim 200, 0 crashed, 10.1 s/rep; config dev/timeseries/gomp_hard.yaml)

gomp_long (nsim 200, 0 crashed, 45.2 s/rep; config dev/timeseries/gomp_long.yaml)

gomp_noobs (nsim 200, 0 crashed, 1.7 s/rep; config dev/timeseries/gomp_noobs.yaml)

gomp_null (nsim 200, 0 crashed, 22.2 s/rep; config dev/timeseries/gomp_null.yaml)

ll_default (nsim 200, 0 crashed, 1.1 s/rep; config dev/timeseries/ll_default.yaml)

ll_hard (nsim 200, 0 crashed, 0.8 s/rep; config dev/timeseries/ll_hard.yaml)

ll_null (nsim 200, 0 crashed, 1.1 s/rep; config dev/timeseries/ll_null.yaml)

marss_corr (nsim 200, 0 crashed, 9.9 s/rep; config dev/timeseries/marss_corr.yaml)

marss_default (nsim 200, 0 crashed, 4.7 s/rep; config dev/timeseries/marss_default.yaml)

marss_hard (nsim 200, 0 crashed, 2.5 s/rep; config dev/timeseries/marss_hard.yaml)

marss_ml (nsim 200, 0 crashed, 1.9 s/rep; config dev/timeseries/marss_ml.yaml)

marss_null (nsim 200, 0 crashed, 4.8 s/rep; config dev/timeseries/marss_null.yaml)

marss_sep (nsim 200, 0 crashed, 4.9 s/rep; config dev/timeseries/marss_sep.yaml)