This appendix is the technical record. For a gentler introduction to each example, see the worked examples.

Appendix S1. Worked examples for every CascadiaTools module

Supporting Information for: CascadiaTools: a validated, AI-guided analysis harness for ecology, wildlife, fisheries and conservation genetics. T. Levi, Cascadia Wildlife Lab.

Each example follows the same workflow an AI agent follows for a researcher. It gives the question, the prompt the researcher typed, the study design, and the true values used to simulate the example dataset. It then shows planted-truth recovery at that design (the acceptance gate), estimates from the example data set against the truth, convergence and goodness of fit, the report's interpretation, and the module's default publication figures. A verdict of DONE means every interpreted estimand passed the gate at this design. OPEN means at least one did not; the reason is stated, and the example shows how the harness reports it. All example data are simulated, and every result can be reproduced with the listed commands.

Regression and causal inference

S1.1 glm: Wood Thrush point counts: negative-binomial GLMM with survey-effort offset and station random intercepts (simulated example)

DONE All interpreted estimands (canopy, elevation, forest-type effects, mean count) pass planted truth at this design (|rel. bias| <= 7.2%, coverage 0.905-0.940); fit converged; DHARMa all p > 0.5; module cross-checked and negative controls in VALIDATION.md. Station SD (FAIL, -17.6% bias, 23% singular fits) and NB theta (unbounded in 7/199) are nuisance parameters and are not interpreted.

Question

Does the number of Wood Thrushes counted change with canopy cover, elevation and forest type?

What the researcher typed

I have Wood Thrush point counts in wood_thrush_counts.csv: 45 stations, each counted 4 times in May-June 2025. Most counts were 10 minutes, some were 5. I want to know whether the number of thrushes changes with canopy cover, elevation and forest type. Please analyse it and give me a report with figures I can use in a paper.

Study design

True values used to simulate the example data

Planted-truth recovery at this design

Estimates from the example data

Convergence and goodness of fit

Interpretation from the report

Per 10 percentage points of canopy cover (%) this is associated with a 39% higher expected Wood Thrush count (rate ratio 1.39, 95% CI 1.18 to 1.64).

Across the observed range of canopy cover (%) (42 to 98), with other continuous predictors at their mean; forest type = deciduous; typical station, the expected Wood Thrush count goes from 0.432 to 2.78.

Per 100 m of elevation (m) this is associated with a 5.1% lower expected Wood Thrush count (rate ratio 0.949, 95% CI 0.854 to 1.05). The 95% interval includes no effect: the data are compatible with effects from 31% lower to 13% higher, including none.

Compared with forest type = deciduous (reference), forest type = conifer is associated with a 52% lower expected Wood Thrush count (rate ratio 0.482, 95% CI 0.248 to 0.935): 0.528 vs 1.1.

Compared with forest type = deciduous (reference), forest type = mixed is associated with a 15% higher expected Wood Thrush count (rate ratio 1.15, 95% CI 0.651 to 2.01): 1.25 vs 1.1. The 95% interval includes no effect: the data are compatible with effects from 35% lower to 101% higher, including none.

Model checks: DHARMa simulated residuals showed no detectable problems with distribution, dispersion, outliers or zeros (all p >= 0.05).

These are effects on the expected count as recorded, which combines how many/whether animals are present with how likely they are to be detected or recorded.

Figures

Limits that matter at this design

• At 45 stations x 4 visits with ~1 bird per count, the station random-intercept SD is biased -18% (23% of simulated fits singular) and NB theta is sometimes unbounded: neither is interpreted (this run).
• Rate ratios for imprecise contrasts are biased upward in the mean (exp of a noisy log ratio): +7.2% conifer, +4.9% mixed here; log-scale coefficients unbiased (Known limit 1).
• Wald intervals only; recovery coverage 0.905-0.940 at this design (Known limit 4).
• Expected counts are for a typical station (random effect 0), not the landscape mean (Known limit 5).
• Counts confound abundance and detection; a canopy effect could partly be detectability (N-mixture module separates them).
• Spatial autocorrelation among stations is not modelled (Known limit 7).
• Associations, not causal effects.

recovery settings

Commands to reproduce

bin/ejob -c 1 -m 2 -t 0:10:00 -l logs/glm_si_truth.log -- bin/Rscript dev/glm/si_truth.R
bin/ejob -c 4 -m 8 -t 0:20:00 -l logs/glm_si_selftest.log -- bin/Rscript core/run.R selftest glm --config modules/glm/example/recovery.yaml --nsim 8 --cores 4 --out dev/glm/si_selftest
bin/ejob -c 16 -m 32 -t 3:00:00 -l logs/glm_si_recover.log -- bin/Rscript core/run.R recover glm --config modules/glm/example/recovery.yaml --nsim 200 --cores 16 --out results/glm/si/recovery
bin/ejob -c 1 -m 2 -t 0:10:00 -l logs/glm_si_reccheck.log -- bin/Rscript dev/glm/si_recovery_check.R
bin/ejob -c 4 -m 16 -t 2:00:00 -l logs/glm_si_fit.log -- bin/Rscript core/run.R fit glm --config modules/glm/example/config.yaml --out results/glm/si/fit
bin/ejob -c 2 -m 8 -t 0:20:00 -l logs/glm_si_figs.log -- bin/Rscript core/run.R figures glm --fit results/glm/si/fit/glm_fit.rds --config modules/glm/example/config.yaml
bin/ejob -c 1 -m 2 -t 0:10:00 -l logs/glm_si_summary.log -- bin/Rscript dev/glm/si_summary.R

Step-by-step walkthrough: modules/glm/example/WALKTHROUGH.md.

S1.2 causal: Does riparian canopy cover cause more juvenile trout? Causal effect from 140 observational stream reaches (DAG + g-computation, simulated data with known truth)

DONE

Question

If riparian canopy cover were increased by 10 percentage points, how many more juvenile trout per 100 m reach would we expect, and how much of that effect runs through cooler August water temperature?

What the researcher typed

I surveyed 140 stream reaches last summer, one electrofishing pass each. The file is stream_trout_reaches.csv. I want to know whether riparian canopy cover actually causes more juvenile trout, and by how much for a 10 % increase in canopy, not just whether they're correlated. Elevation, farmland in the catchment and stream width all affect both canopy and trout. I think shade works partly by keeping the water cool. I also measured periphyton, which depends on both shade and how many trout graze. Please do a proper causal analysis and give me a report with figures I can use in a paper.

Study design

True values used to simulate the example data

Planted-truth recovery at this design

Estimates from the example data

Convergence and goodness of fit

Interpretation from the report

Average causal effect: increasing riparian canopy cover by 10 percentage points at every site would change the expected number of juvenile trout per 100 m reach by 1.18 (95% CI 0.258 to 2.10), from 6.30 to 7.48 on average. As a ratio (rate ratio), 1.19 (1.05 to 1.34): a 19% increase.

Controlled direct effect (holding August water temperature fixed): 0.566 (95% CI -0.296 to 1.43) on the count scale. On the model scale about 50% of the total effect runs through the mediator(s) (1 - direct/total = 0.50), a rough split that ignores the uncertainty in both.

Strategies the DAG says are BIASED give: No adjustment (crude association) -> 1.78 (+51%); Omit elevation_m (a confounder) -> 2.14 (+81%); Omit wetted_width_m (a confounder) -> 0.530 (-55%); Also adjust for mediator water_temp_c (over-control) -> 0.552 (-53%); Also adjust for periphyton_chla (collider) -> 1.59 (+34%).

Sensitivity to unmeasured confounding (E-value): to fully explain away this effect, an unmeasured common cause would have to raise both riparian canopy cover and juvenile trout per 100 m reach by a factor of 1.66 (risk-ratio scale) beyond the adjusted variables; to make the 95% CI include no effect, a factor of 1.29.

Checking the DAG against the data: the DAG implies 16 conditional independencies that can be tested. None is contradicted by the data (Fisher's C = 31.4, df = 32, p = 0.50). Passing these tests does not prove the DAG.

Model check: the negbin outcome model shows problems (DHARMa dispersion p < 0.05); the effect estimate may be off. Check the residual plot; a single p < 0.05 among several checks also happens by chance in 1 of 20 good models.

Figures

Limits that matter at this design

• The DAG is untestable where it matters most: local tests detect missing arrows between measured variables (95% power at 120 sites) but not an unmeasured confounder of canopy and trout; the E-value (1.66; 1.29 for the CI) is reported for that reason.
• A moderate collider (periphyton) can pass on the ratio scale (+5.1% bias) while failing on the count-scale ATE (+18%) and beta; the DAG, not the data, must exclude it.
• Small samples: the count-scale ATE is biased upward (+18%) at 40 sites; at this design (140 reaches) it passed (+1.8%).
• Only controlled direct effects are estimated; the ~50% share through water temperature is a rough log-scale split, not a natural indirect effect. The direct effect has only 51-60% power at 140 reaches.
• No exposure x covariate interactions; confounders act linearly on the log scale; simulation tested only linear-Gaussian non-outcome nodes.
• Detection is not modelled (one electrofishing pass); if capture efficiency changes with canopy, the effect mixes ecology and detection.
• DHARMa flagged dispersion (p = 0.016) on this heavy-tailed draw; bootstrap intervals are ~35% wider and only just exclude no effect (ATE 0.03 to 2.52).

recovery settings

estimates bootstrap sensitivity

Commands to reproduce

bin/ejob -c 1 -m 4 -t 0:20:00 -l logs/causal_make_examples.log -- bin/Rscript dev/causal/make_examples.R
bin/ejob -c 16 -m 32 -t 3:00:00 -l logs/causal_si_recover.log -- bin/Rscript core/run.R recover causal --config modules/causal/example/recovery.yaml --nsim 200 --cores 16 --out results/causal/si/recovery
bin/ejob -c 4 -m 16 -t 2:00:00 -l logs/causal_si_fit.log -- bin/Rscript core/run.R fit causal --config modules/causal/example/config.yaml --out results/causal/si/fit
bin/ejob -c 4 -m 16 -t 1:00:00 -l logs/causal_si_boot.log -- bin/Rscript core/run.R fit causal --config dev/causal/si_boot/config_bootstrap.yaml --out dev/causal/si_boot/fit

Step-by-step walkthrough: modules/causal/example/WALKTHROUGH.md.

S1.3 sem: Fire severity, vegetation and songbirds after a wildfire (piecewise SEM with a Poisson GLMM; simulated data with known truth)

OPEN

Question

How does wildfire severity affect breeding songbird detections, directly and indirectly through canopy loss and shrub regrowth, and does soil moisture affect birds through shrubs?

What the researcher typed

I have a CSV of 120 forest plots from 12 burn units after a wildfire (burned_forest_plots.csv). For each plot I have fire severity (dNBR), canopy cover, soil moisture, shrub cover and the number of birds detected in one point count. I think fire reduces canopy, less canopy means more shrubs, more shrubs means more birds, wetter soil means more shrubs, and fire might also reduce birds directly. Can you test this diagram, tell me how strong each link is, and how much of the fire effect on birds goes through the vegetation? Plots in the same burn unit are not independent. Please give me a report and figures for a paper.

Study design

True values used to simulate the example data

Planted-truth recovery at this design

Estimates from the example data

Convergence and goodness of fit

Interpretation from the report

Global fit: Fisher's C = 11.06 (df = 8), P = 0.198: no evidence that the model leaves out a needed link. This supports, but does not prove, the causal structure drawn.

Fire severity -> canopy cover: an increase of 100 dNBR is associated with a difference of -5.87 (95% CI -6.61 to -5.14) % canopy cover.

Canopy cover -> shrub cover: an increase of 10 % canopy is associated with a difference of -2.63 (95% CI -3.41 to -1.85) % shrub cover.

Soil moisture -> shrub cover: an increase of 5 % moisture is associated with a difference of 3.29 (95% CI 2.12 to 4.45) % shrub cover.

Shrub cover -> bird detections: an increase of 10 % shrub is associated with expected bird detections multiplied by 1.33 (95% CI 1.20 to 1.47; +33%).

Fire severity -> bird detections (direct): an increase of 100 dNBR is associated with expected bird detections multiplied by 0.91 (95% CI 0.87 to 0.95; -9%).

Per 1-SD increase in fire severity (212 dNBR), the indirect route alone multiplies expected bird detections by 1.10 (95% CI 1.05 to 1.15); the total association corresponds to x0.90 (95% CI 0.83 to 0.99). Direct and indirect associations have opposite signs and partly cancel; a weak total does not mean the driver has little influence.

Soil moisture affects birds only indirectly through shrubs: x1.10 (95% CI 1.05 to 1.16) per 1 SD (5.02 %) of moisture.

Bird counts mix abundance with detectability; the SEM cannot separate an effect on birds from an effect on detecting them.

Paths are associations unless the diagram is causally complete; passing the d-sep test does not prove it.

Figures

Limits that matter at this design

• Planted-truth gate FAILED 6/24 at this design: standardised coefficients involving the Poisson GLMM bird component are inflated by about +17% (coverage of std shrub -> birds 0.880). Cause (dev/sem/si/std_count_check.R): the observation-empirical scaling sd(eta)/cor(y, fitted) uses fitted values that include the 12 burn-unit BLUPs, inflating cor (0.664 vs 0.623); fixed-only fitted values give +2.2%. Matches piecewiseSEM::coefs(); VALIDATION.md Known limit 3 (reported there as ~+7%). Report raw paths, rate ratios and effects in original units; use standardised bird values only to rank routes.
• Omitted paths or unmeasured confounders bias indirect effects (~-42% in the negative control); d-sep detects an omitted path of standardised size ~0.3 only 73% of the time at n = 100 (Known limit 5).
• Counts mix abundance and detectability (README not_for); use N-mixture/occupancy with repeat visits to separate them.
• Interactions, nonlinear paths and feedback loops are not supported (Known limit 6).
• Total fire -> birds has only 9% power at this design because direct and indirect routes nearly cancel; in this data set its 95% CI (-0.00090 to -0.00007) misses the true -0.00012, as does canopy -> shrub (CI -0.341 to -0.185 vs truth -0.35): 7 of 9 raw intervals contain the truth.
• Single d-sep claims can be false alarms: moisture -> birds P = 0.023 although no such link exists in the truth.

Commands to reproduce

bin/ejob -c 1 -m 4 -t 0:10:00 -l logs/sem_si_regen.log -- bin/Rscript dev/sem/si/make_example_check.R
bin/ejob -c 16 -m 32 -t 3:00:00 -l logs/sem_si_recover.log -- bin/Rscript core/run.R recover sem --config modules/sem/example/recovery.yaml --nsim 200 --cores 16 --out results/sem/si/recovery
bin/ejob -c 8 -m 16 -t 0:30:00 -l logs/sem_si_stdcheck.log -- bin/Rscript dev/sem/si/std_count_check.R
bin/ejob -c 4 -m 16 -t 2:00:00 -l logs/sem_si_fit.log -- bin/Rscript core/run.R fit sem --config modules/sem/example/config.yaml --out results/sem/si/fit
bin/Rscript core/run.R figures sem --fit results/sem/si/fit/sem_fit.rds --config modules/sem/example/config.yaml --out results/sem/si/fit

Step-by-step walkthrough: modules/sem/example/WALKTHROUGH.md.

Occupancy

S1.4 occu-single: Pileated woodpecker occupancy and forest cover (simulated call-playback surveys, 120 stands)

DONE

Question

Is pileated woodpecker occupancy of forest stands associated with forest cover and elevation, allowing for birds that were present but missed (detection varying with survey date and observer)?

What the researcher typed

I have pileated woodpecker call-playback surveys from 120 forest stands, up to 4 visits each in May-June 2025. The files are detections.csv (one row per visit: stand, visit number, date, day of year, observer, detected 1/0) and sites.csv (forest cover %, elevation in m, stand age). Does forest cover or elevation affect whether woodpeckers occupy a stand? Please account for birds we missed, and allow detection to change with date and observer. I need results and figures for a paper.

Study design

True values used to simulate the example data

Planted-truth recovery at this design

Estimates from the example data

Convergence and goodness of fit

Interpretation from the report

The species was detected at least once at 34 of 120 sites (naive occupancy 0.28). Allowing for missed detections, the estimated probability that a site with average covariate values is occupied is 0.33 (95% CI 0.20 to 0.49).

Of the 120 surveyed sites, an estimated 47 (95% CI 37 to 58) were actually occupied (proportion 0.39 (95% CI 0.31 to 0.48)); this finite-sample estimate applies to these sites only, not to the wider landscape.

With the typical 4 surveys per site, an occupied site is detected at least once with probability 0.87 (95% CI 0.70 to 0.96), so about 13% of occupied sites go undetected with this survey effort.

Going from the 10th to the 90th percentile of forest cover (22 to 85.3%), occupancy changes from 0.09 (95% CI 0.03 to 0.28) to 0.62 (95% CI 0.36 to 0.81); each 10-unit increase multiplies the odds by 1.54.

Going from the 10th to the 90th percentile of elevation (361 to 1250 m), occupancy changes from 0.69 (95% CI 0.41 to 0.88) to 0.12 (95% CI 0.04 to 0.27).

Detection per survey by observer: AB 0.39 (95% CI 0.26 to 0.55); JK 0.13 (95% CI 0.05 to 0.26); MR 0.52 (95% CI 0.30 to 0.71).

These are associations from an observational design; they do not by themselves show that a covariate causes the change.

Goodness of fit (MacKenzie-Bailey test, 500 bootstrap replicates): c-hat = 1.31, p = 0.138. No evidence of lack of fit; this test has low power against unmodelled differences in detection among sites.

Figures

Limits that matter at this design

• Small samples inflate occupancy slopes: at 120 stands the forest effect was +6.0% (Run A) and +9.6% (Run B) too strong on average, just inside the 10% gate (100 sites fails at +11.5%). Trust the direction; treat the size as possibly a little high.
• Low power for the secondary effects at this design: an elevation effect of the simulated size was detected in 43% of replicates, observer differences in 16-32%.
• Unmodelled detection heterogeneity among sites biases occupancy low (about -23% at logit SD 1.5) and the MacKenzie-Bailey test flags it only about 1 time in 10.
• Few surveys with low detection fail (2 surveys, p = 0.2); not the case here (cumulative detection 0.87 over 4 visits), though 27 stands had only 2 visits.
• The module simulator has one occupancy and one detection covariate, so the standard recovery (Run A) omits elevation and observer; Run B (dev/occu-single/si_full_recovery.R) covers the full example model.
• No false positives, no spatial autocorrelation, closure assumed within the May-June season.

Commands to reproduce

bin/Rscript dev/occu-single/make_example.R
bin/ejob -c 1 -m 4 -t 0:10:00 -l logs/occu-single_si_truth.log -- bin/Rscript dev/occu-single/si_truth.R
bin/ejob -c 16 -m 32 -t 3:00:00 -l logs/occu-single_si_recover.log -- bin/Rscript core/run.R recover occu-single --config modules/occu-single/example/recovery.yaml --nsim 200 --cores 16 --out results/occu-single/si/recovery
bin/ejob -c 16 -m 32 -t 1:00:00 -l logs/occu-single_si_fullrecover.log -- bin/Rscript dev/occu-single/si_full_recovery.R 200 16
bin/ejob -c 4 -m 16 -t 2:00:00 -l logs/occu-single_si_fit.log -- bin/Rscript core/run.R fit occu-single --config modules/occu-single/example/config.yaml --out results/occu-single/si/fit

Step-by-step walkthrough: modules/occu-single/example/WALKTHROUGH.md.

S1.5 occu-rota: Do gray foxes avoid coyotes, and does housing density change that? Multispecies occupancy with species interactions (Rota et al. 2016)

DONE

Question

Are gray foxes less likely to occupy camera sites where coyotes are present, and does that avoidance weaken with housing density, after allowing for imperfect detection and cameras that were down?

What the researcher typed

I have camera-trap data from 220 stations, 6 weekly occasions in June-July 2024. The files are camera_detections.csv (one row per camera and week: days the camera worked, and 1/0 for coyote, gray fox and bobcat) and camera_sites.csv (forest cover within 500 m, housing units per km2, whether the camera is on a trail). Are gray foxes less likely to occupy sites where coyotes are present, and does that change with housing density? Please account for animals the cameras missed and for cameras that were down. Detection probably depends on days active and on whether the camera faces a trail. I need results and figures for a paper.

Study design

True values used to simulate the example data

Planted-truth recovery at this design

Estimates from the example data

Convergence and goodness of fit

Interpretation from the report

Allowing for missed detections, the estimated probability that each species occupies a site with average covariate values is: coyote 0.60 (95% CI 0.50 to 0.69); gray fox 0.25 (95% CI 0.18 to 0.34).

At average covariate values, coyote and gray fox occur together less often than expected if they occurred independently (negative co-occurrence; this pattern alone does not show that one species avoids the other).

Where coyote is present, gray fox occupies a site with probability 0.15 (95% CI 0.08 to 0.27); where coyote is absent, 0.41 (95% CI 0.27 to 0.57).

The co-occurrence odds ratio is 0.25 (95% CI 0.09 to 0.66); a value of 1 would mean no association.

The coyote-gray fox co-occurrence odds ratio is multiplied by 2.06 (95% CI 0.85 to 4.96) per 1 SD (0.984) increase in log(housing_per_km2), i.e. by 1.66 (95% CI 0.89 to 3.09) per doubling of housing_per_km2 (the interval includes 1, so the data do not show that the association changes along this gradient).

At low housing_per_km2 (10th percentile), P(gray fox | coyote present) = 0.04 (95% CI 0.01 to 0.14) vs coyote absent 0.27 (95% CI 0.13 to 0.49); at high housing_per_km2 (90th percentile), 0.42 (95% CI 0.25 to 0.61) vs 0.55 (95% CI 0.32 to 0.77).

Across the observed range of forest_pct (10th percentile 15.8 to 90th percentile 87), marginal occupancy changes: coyote from 0.71 (95% CI 0.57 to 0.83) to 0.46 (95% CI 0.34 to 0.59); gray fox from 0.12 (95% CI 0.06 to 0.22) to 0.46 (95% CI 0.32 to 0.59). This is an association, not necessarily a causal effect.

Per-survey detection probability where the species is present (average survey conditions): coyote 0.23 (95% CI 0.18 to 0.29); gray fox 0.30 (95% CI 0.22 to 0.39). With 6 surveys, the chance of detecting a species at least once at an occupied site is coyote 0.80, gray fox 0.88.

Caution: co-occurrence is not proof of a species interaction (Blanchet et al. 2020).

Goodness of fit: no evidence of lack of fit (parametric bootstrap, p = 0.25, c-hat 1.12).

Figures

Limits that matter at this design

• Interaction slopes are inflated by about 10-13% at a few hundred sites (VALIDATION known limit 1); +10.3% in Run A, +3.9% in Run B at this design. Trust the sign and the conditional-probability curves more than the slope's size.
• Low power at 220 cameras: the average-site interaction is detected in about 54% of replicates and its housing slope in only 26-27%; a non-significant slope is not evidence of no change.
• Co-occurrence is not interaction: a negative association can come from avoidance, predation or unmodelled habitat differences; the validation assumes the model is true.
• Rare species (marginal occupancy ~0.12) bias the interaction at 400 sites; gray fox (~0.25) is above that level, but recheck recovery if a rarer species is added.
• Boundary fits are flagged (|estimate| > 8 or SE > 10); none here (largest SE 0.50).
• Three-way interactions and factor covariates on interactions are not validated; not used here.
• Run A uses the module simulator with one site and one survey covariate (standing in for log housing and days active); Run B uses the full example generating process and passed 21/21 pooled over 400 replicates (one marginal fail, gray fox days-active detection slope +12.0%, in the first 200-replicate batch).

Commands to reproduce

bin/ejob -c 1 -m 4 -t 0:10:00 -l logs/occu-rota_si_truth.log -- bin/Rscript dev/occu-rota/si/si_truth.R
bin/ejob -c 16 -m 32 -t 3:00:00 -l logs/occu-rota_si_recover.log -- bin/Rscript core/run.R recover occu-rota --config modules/occu-rota/example/recovery.yaml --nsim 200 --cores 16 --out results/occu-rota/si/recovery
bin/ejob -c 16 -m 32 -t 1:00:00 -l logs/occu-rota_si_fullrecover.log -- bin/Rscript dev/occu-rota/si/si_full_recovery.R 200 16
bin/ejob -c 16 -m 32 -t 1:00:00 -l logs/occu-rota_si_fullrecover_s2027000.log -- bin/Rscript dev/occu-rota/si/si_full_recovery.R 200 16 2027000
bin/ejob -c 1 -m 4 -t 0:10:00 -l logs/occu-rota_si_fullrecover_pool.log -- bin/Rscript dev/occu-rota/si/pool_full.R
bin/ejob -c 4 -m 16 -t 2:00:00 -l logs/occu-rota_si_fit.log -- bin/Rscript core/run.R fit occu-rota --config modules/occu-rota/example/config.yaml --out results/occu-rota/si/fit
python3 dev/occu-rota/si/make_si_summary.py

Step-by-step walkthrough: modules/occu-rota/example/WALKTHROUGH.md.

S1.6 occu-waddle: Do invasive Cuban treefrogs push out native green treefrogs? Two-species conditional occupancy from PVC-pipe surveys (simulated data)

DONE

Question

Is green treefrog occupancy lower at sites where Cuban treefrogs occur, and is green treefrog detection lower on visits when Cuban treefrogs are found in the pipes, after allowing for both species being missed?

What the researcher typed

I have pipe-survey data for Cuban treefrogs and green treefrogs at 180 sites, 5 visits each, in pipe_surveys.csv, and site information (distance to road, canopy cover) in sites.csv. I want to know whether green treefrogs occupy fewer sites where Cuban treefrogs are present, and whether green treefrogs are harder to detect on visits when Cuban treefrogs are found. Cuban treefrog detection probably depends on air temperature. Please run the analysis and give me a report with figures I can use in a paper.

Study design

True values used to simulate the example data

Planted-truth recovery at this design

Estimates from the example data

Convergence and goodness of fit

Interpretation from the report

Where Cuban treefrog is present, green treefrog occupies 0.31 (95% CrI 0.21 to 0.43) of sites; where Cuban treefrog is absent, 0.67 (95% CrI 0.56 to 0.78).

The odds of green treefrog occupying a site are 0.21 (95% CrI 0.10 to 0.44) times as high where Cuban treefrog is present.

The species interaction factor is 0.64 (95% CrI 0.46 to 0.81): at average conditions the two species co-occur less often than expected if they were distributed independently (a segregation pattern).

Occupancy of Cuban treefrog: the odds are multiplied by 0.40 (95% CrI 0.26 to 0.57) per 1 SD of log(dist_road_m); it decreases from 0.96 at 29 m to 0.09 at 5410 m.

Occupancy of green treefrog where Cuban treefrog is absent: the odds are multiplied by 2.13 (95% CrI 1.35 to 3.63) per 1 SD of canopy cover (from 0.33 at 3% to 0.89 at 89%); where Cuban treefrog is present there is no clear change (1.14, 95% CrI 0.67 to 1.98).

Where both occur, green treefrog is detected on 0.23 (95% CrI 0.13 to 0.36) of surveys that detected Cuban treefrog and on 0.50 (95% CrI 0.36 to 0.64) of surveys that missed it; log-odds ratio -1.18 (95% CrI -2.01 to -0.40). This is a detection effect, not an occupancy effect.

A difference between psiBA and psiBa, or an SIF away from 1, is a co-occurrence pattern, consistent with competition or predation but equally with both species responding to habitat features not in the model; it does not by itself show that Cuban treefrog excludes green treefrog.

Figures

Limits that matter at this design

• Power for the detection interaction is modest: 0.61 for a log-odds ratio of -1.05 at 180 sites x 5 visits.
• Canopy effects on green treefrog occupancy have low power (0.39 and 0.53); the psiBa slope magnitude is uncertain by about +-15% (VALIDATION.md sections 1, 8). In this data set the psiBA canopy slope was estimated at 0.13 vs true 0.60 (interval covers truth).
• The dominant/subordinate direction is assumed (Cuban treefrog dominant); the model does not test it.
• Co-occurrence is not interaction: shared unmeasured habitat drivers produce psiBA < psiBa (Blanchet et al. 2020).
• A detection interaction (rBA < rBa) can come from shared survey conditions affecting both species, not only behaviour.
• Small-sample bias toward 'no interaction' (psiBA, P(both), SIF) is documented at 60 sites; at this design (180 x 5) all were within +-2.5%.
• Not validated: detection < 0.3, factor covariates, covariates on rBA/rBa (none used here).

truth realised

recovery meta

estimates cover truth

18/18

Commands to reproduce

bin/ejob -c 1 -m 4 -t 0:15:00 -l logs/occu-waddle_si_truth.log -- bin/Rscript dev/occu-waddle/si_truth.R
bin/ejob -c 4 -m 8 -t 0:30:00 -l logs/occu-waddle_si_selftest.log -- bin/Rscript core/run.R selftest occu-waddle --config modules/occu-waddle/example/recovery.yaml --nsim 4 --cores 4 --out dev/occu-waddle/si_selftest
bin/ejob -c 16 -m 32 -t 3:00:00 -l logs/occu-waddle_si_recover.log -- bin/Rscript core/run.R recover occu-waddle --config modules/occu-waddle/example/recovery.yaml --nsim 200 --cores 16 --out results/occu-waddle/si/recovery
bin/ejob -c 1 -m 4 -t 0:10:00 -l logs/occu-waddle_si_mdtable.log -- bin/Rscript dev/occu-waddle/mdtable.R results/occu-waddle/si/recovery
bin/ejob -c 4 -m 16 -t 2:00:00 -l logs/occu-waddle_si_fit.log -- bin/Rscript core/run.R fit occu-waddle --config modules/occu-waddle/example/config.yaml --out results/occu-waddle/si/fit
bin/Rscript core/run.R figures occu-waddle --fit results/occu-waddle/si/fit/occu-waddle_fit.rds --config modules/occu-waddle/example/config.yaml

Step-by-step walkthrough: modules/occu-waddle/example/WALKTHROUGH.md.

S1.7 occu-spatial: Fisher occupancy, canopy cover and elevation with residual spatial clustering (simulated camera-trap survey, 160 stations)

DONE DONE for the estimands that answer the question (covariate effects, occupancy, detection, mean occupancy, proportion occupied): all pass planted truth at this design (100/100 converged), the real fit converged and passed the PPC; spatial variance and range fail recovery on bias and are reported descriptively only, not interpreted.

Question

Is fisher occupancy associated with canopy cover and elevation, allowing for missed detections and for residual spatial clustering, and where in the 30 x 30 km block are fishers likely to occur?

What the researcher typed

I have fisher camera-trap data from 160 stations in the western Cascades, summer 2024, up to 6 weekly occasions per camera (some cameras died early). fisher_detections.csv has one row per station-week (station, week, start date, days the camera worked, detected 1/0). camera_stations.csv has UTM coordinates, canopy cover (%), elevation (m) and ownership. prediction_grid_1km.csv is a 1-km grid with the same covariates. Is fisher occupancy related to canopy cover and elevation? Please account for fishers we missed and for spatial clustering, and make a map of predicted occupancy. I need results and figures for a paper.

Study design

True values used to simulate the example data

Planted-truth recovery at this design

Estimates from the example data

Convergence and goodness of fit

Interpretation from the report

Occupancy probability at average covariate values and zero spatial effect is 0.50 (95% CrI 0.36 to 0.65).

Averaged over the 160 surveyed sites, mean occupancy probability is 0.51 (95% CrI 0.43 to 0.59), and the estimated proportion of surveyed sites that are occupied is 0.51 (95% CrI 0.46 to 0.57), compared with 0.44 of sites where the species was detected at least once (the difference is the correction for missed detections).

Occupancy and Canopy cover (%): each 1-SD increase (1 SD = 12.9 in original units) multiplies the odds of occupancy by 2.75 (95% CrI 1.65 to 5.26); predicted occupancy goes from 0.28 at 1 SD below the mean to 0.73 at 1 SD above it (other covariates at their means). Per 10-unit increase the odds of occupancy are multiplied by 2.19 (95% CrI 1.48 to 3.64).

Occupancy and Elevation (m): each 1-SD increase (1 SD = 246 in original units) multiplies the odds of occupancy by 0.58 (95% CrI 0.32 to 0.97); predicted occupancy goes from 0.63 at 1 SD below the mean to 0.37 at 1 SD above it. Per 200-unit increase the odds of occupancy are multiplied by 0.64 (95% CrI 0.40 to 0.98).

Detection: when the species is present, a single survey detects it with probability 0.30 (95% CrI 0.25 to 0.35) (at average survey conditions).

Detection and Days camera active in the week: per 1-unit increase the odds of detection are multiplied by 1.48 (95% CrI 1.12 to 2.09). This is a detection effect (how easily the species is found), not a habitat effect.

The residual spatial variance is 0.95 (95% CrI 0.17 to 3.89) on the logit scale; its prior alone gives a 95% interval of 0.18 to 4.13, so the data add little information: these numbers mostly restate the prior and do not show whether residual clustering is strong or weak.

Model check: posterior predictive p-values (Freeman-Tukey, by site) are 0.55; no sign of lack of fit.

All effects are associations from observational data, not causal effects.

Figures

Limits that matter at this design

• Spatial variance and effective range are not recoverable at this design (+48% and +21% bias in planted truth) and their posteriors restate the priors; read the spatial maps as residual smoothing, not a measured clustering scale.
• Occupancy slopes lean away from zero at modest sample sizes (+8% canopy, +11.5% elevation in planted truth; 10-28% at 100-120 sites in VALIDATION.md; ~10% inflation if no real clustering). Trust direction; report the non-spatial PGOcc fit alongside (canopy 0.87 vs 1.01 per SD).
• Low power for the elevation effect: an effect of -0.40 per SD was detected in only 33% of replicates; the real-data interval only just excludes zero.
• Spatial confounding: smooth covariates such as elevation (15-km range) can have part of their effect absorbed by the spatial term (VALIDATION.md: -20% for a 90%-spatial covariate).
• Recovery simulator differs from the example generator: independent covariates, missing weeks at random rather than camera failures, standard-normal days_active.
• Closure over the season, no false positives and no unmodelled detection heterogeneity are assumed; stations 1 km apart may share individuals, so occupancy is better read as use.
• Core gate caveat: the '>= 90% converged' clause cannot fail for non-convergence; convergence (100/100) was counted by hand.

recovery settings

estimates nonspatial PGOcc

Commands to reproduce

bin/ejob -c 1 -m 4 -t 0:20:00 -l logs/occu-spatial_make_example.log -- bin/Rscript dev/occu-spatial/make_example.R
bin/ejob -c 1 -m 4 -t 0:20:00 -l logs/occu-spatial_si_truth.log -- bin/Rscript dev/occu-spatial/si_truth.R
bin/ejob -c 16 -m 32 -t 3:00:00 -l logs/occu-spatial_si_recover.log -- bin/Rscript core/run.R recover occu-spatial --config modules/occu-spatial/example/recovery.yaml --nsim 100 --cores 16 --out results/occu-spatial/si/recovery
bin/ejob -c 4 -m 16 -t 2:00:00 -l logs/occu-spatial_si_fit.log -- bin/Rscript core/run.R fit occu-spatial --config modules/occu-spatial/example/config.yaml --out results/occu-spatial/si/fit
bin/ejob -c 4 -m 16 -t 1:00:00 -l logs/occu-spatial_si_fit_pgocc.log -- bin/Rscript core/run.R fit occu-spatial --config dev/occu-spatial/si_config_pgocc.yaml --out results/occu-spatial/si/fit_pgocc

Step-by-step walkthrough: modules/occu-spatial/example/WALKTHROUGH.md.

S1.8 abund-interactions: Abundance-mediated species interactions: pine marten abundance and red squirrel occupancy (simulated camera-trap example)

OPEN Planted truth at the example design (260 sites) fails for gamma0 (psi_beta_dominant_abundance): relative bias -11.0% (limit 10%; 3.2 MCSE) with coverage 0.95. Known limit 1 in VALIDATION.md. Convergence, GOF, cross-checks, negative control and reproducibility pass. Read the gamma0 interval, not its point estimate.

Question

Is red squirrel occupancy lower at sites used by more pine martens, and does that association depend on the share of native broadleaf woodland?

What the researcher typed

I have camera-trap data from 260 woodland sites in Northern Ireland: deployments.csv (camera start and end dates), records.csv (every photo, with species) and sites.csv (woodland cover and broadleaf share within 1 km). Is red squirrel occupancy lower where there are more pine martens, and does that depend on how much of the woodland is native broadleaf? Please give me a report with figures I can use in a paper.

Study design

True values used to simulate the example data

Planted-truth recovery at this design

Estimates from the example data

Convergence and goodness of fit

Interpretation from the report

Each 1-SD increase in woodland_pct (1 SD = 17.4 units) is associated with 1.70 (95% CI 1.47 to 1.96) times the expected Pine marten abundance (a +70% change) (per 10 units: rate ratio 1.36 (95% CI 1.25 to 1.47)).

Interaction at average habitat: each additional Pine marten individual at a site is associated with 0.59 (95% CI 0.39 to 0.85) times the odds that Red squirrel occupies the site (odds ratio per individual); on the logit scale the coefficient is -0.53 (95% CI -0.93 to -0.17) per individual. The posterior probability that the association is negative is > 0.99.

Predicted Red squirrel occupancy at average habitat: 0.62 (0.52-0.72) at a site with no Pine marten, 0.50 (0.42-0.57) with 1, 0.37 (0.25-0.49) with 2 and 0.26 (0.12-0.44) with 3 individuals (posterior medians, 95% CI).

Habitat dependence: the interaction changes by 0.69 (95% CI 0.30 to 1.14) logit units per individual for each 1-SD (20.8 units) increase in broadleaf_pct; the 95% interval excludes zero, so the data support an interaction whose strength depends on habitat.

Where broadleaf_pct = 12.9 (10th percentile): odds ratio per Pine marten individual 0.25 (95% CI 0.12 to 0.50) (negative).

Where broadleaf_pct = 65.2 (90th percentile): odds ratio per Pine marten individual 1.42 (95% CI 0.84 to 2.53) (not clearly different from zero).

Comparison with the common two-step shortcut (estimate Pine marten abundance first, then use it as a fixed covariate): the shortcut gives -0.48 (95% CI -0.81 to -0.15), with an interval 86% as wide as the joint model's. The joint model carries the uncertainty in each site's abundance into the interaction; the shortcut treats estimated abundance as if it were known, and in simulations it is biased towards zero with too-narrow intervals.

Caution: this is an association between the dominant's estimated abundance and the subordinate's occupancy, not proof of predation or competition. Unmodelled habitat preferences shared (or opposed) by the two species produce the same signal, so include the main habitat drivers of both species. Royle-Nichols abundance assumes individuals are detected independently and abundance is Poisson among sites; absolute abundance is less reliable than relative differences. Simulation checks for this module (VALIDATION.md) recovered the interaction without bias at 600 sites (4 weekly occasions of daily records), including when the dominant's per-individual daily detection was as low as 0.03; in an independent check with 700 sites and a weaker interaction its size was overstated by about 12%, and with 150-300 sites by roughly 10-20% (credible intervals kept near-nominal coverage). With fewer than about 1000 sites, read the interval rather than the point estimate.

Figures

Limits that matter at this design

• 260 sites is below the ~600 sites validated for a habitat-dependent interaction; at this design gamma0 is overstated by ~11% (coverage stays 0.95). Read intervals, not point estimates.
• The broadleaf main effect on red squirrel occupancy is weakly identified at this design (planted-truth power 0.34).
• Association, not causation: unmodelled habitat shared by both species loads onto the interaction term.
• Royle-Nichols assumptions (Poisson abundance, independent detection of individuals, closure over 4 weeks) were not stress-tested; group-living or overdispersed martens would bias abundance and the interaction. The posterior predictive check is the only guard.
• Abundance is the number of martens using a camera site, not a density; one marten can use several sites.
• The two-step plug-in shortcut fails the gate (attenuated, overconfident) and is shown for comparison only.

Commands to reproduce

bin/ejob -c 1 -m 4 -t 0:20:00 -l logs/abund-interactions_make_example.log -- bin/Rscript dev/abund-interactions/make_example.R
bin/ejob -c 16 -m 32 -t 3:00:00 -l logs/abund-interactions_si_recover.log -- bin/Rscript core/run.R recover abund-interactions --config modules/abund-interactions/example/recovery.yaml --nsim 100 --cores 16 --out results/abund-interactions/si/recovery
bin/ejob -c 4 -m 16 -t 2:00:00 -l logs/abund-interactions_si_fit.log -- bin/Rscript core/run.R fit abund-interactions --config modules/abund-interactions/example/config.yaml --out results/abund-interactions/si/fit
../.venv/bin/python dev/abund-interactions/si/make_si_summary.py

Step-by-step walkthrough: modules/abund-interactions/example/WALKTHROUGH.md.

Abundance

S1.9 nmixture: Hermit warbler abundance from repeated point counts (static binomial N-mixture model)

DONE

Question

How many Hermit warblers are at 120 Coast Range survey points, allowing for birds present but not counted, and is abundance associated with canopy cover and elevation?

What the researcher typed

I have Hermit warbler point counts from 120 points, each visited up to 3 times between mid-May and mid-June. The counts are in point_counts.csv (one row per visit, with wind, day of year and observer). The habitat for each point is in points.csv (canopy cover and elevation). Please estimate how many warblers there are, allowing for birds we missed, and tell me whether abundance is related to canopy cover and elevation.

Study design

True values used to simulate the example data

Planted-truth recovery at this design

Estimates from the example data

Convergence and goodness of fit

Interpretation from the report

At a site with average covariate values the model expects 2.00 (95% CI 1.69 to 2.38) Hermit warblers.

Summed over the 120 surveyed sites, estimated abundance is 272 (95% CI 244 to 308) (an average of 2.26 (95% CI 2.03 to 2.57) per site).

The sum of the largest count at each site was 217, so the model infers that about 20% of the animals present were missed even on the visit with the highest count.

Each increase of 1 SD (19.74 units) in Canopy cover is associated with expected abundance being multiplied by 1.59 (95% CI 1.37 to 1.84), i.e. +59% (95% CI +37% to +84%).

Across the observed range (16.00 to 98.00) expected abundance per site goes from 0.71 to 4.90, other covariates held at their mean.

Each increase of 1 SD (293.79 units) in Elevation is associated with expected abundance being multiplied by 1.19 (95% CI 1.05 to 1.36), i.e. +19% (95% CI +5% to +36%).

Detection: on a visit with average covariate values (Observer 'JDK') each animal present had a 0.51 (95% CI 0.41 to 0.61) probability of being counted; over 3 such visits the probability of counting an animal at least once is 0.89 (0.80 to 0.94).

Each increase of 1 SD (1.10 units) in Wind multiplies the odds of detecting an animal by 0.49 (95% CI 0.40 to 0.60).

Model fit: no significant lack of fit at the 5% level, but the p-value is borderline (c-hat 1.17, smaller of the chi-square and Freeman-Tukey p-values = 0.06).

Relative differences among sites and covariate effects (rate ratios) are more robust than the absolute totals. All covariate relationships are associations, not demonstrated causes.

Figures

Limits that matter at this design

• Absolute N assumes closure, no double counting and independent detection; double counting and closure violations were not validated and GOF rarely detects them (VALIDATION.md limit 4).
• Unmodelled detection heterogeneity biases N downward with overconfident intervals; strong site-level heterogeneity is flagged by GOF 95% of the time, individual-level heterogeneity was not tested. GOF p = 0.06 here is borderline (limit 2).
• Overdispersed abundance fitted as Poisson biases N by about -20%; here AIC favoured Poisson and NB/ZIP gave the same abundance (limit 3).
• Low detection (p ~0.2) with few visits biases absolute N; this design (p ~0.5, 3 visits) passed recovery (limit 1).
• Recovery used fast settings and one covariate per submodel; elevation, day-of-year and observer effects were checked only against this example's known truth (limit 8).
• In this one data set two detection intervals miss the truth (wind OR 0.49 [0.40, 0.60] vs 0.61; observer MLR OR 1.57 [1.03, 2.40] vs 1.00), each about 2 SE off; recovery coverage for wind was 0.97, so this is sampling chance, but single detection covariates should be read cautiously.

Commands to reproduce

bin/ejob -c 1 -m 4 -t 0:20:00 -l logs/nmixture_make_example.log -- bin/Rscript dev/nmixture/make_example.R
bin/ejob -c 1 -m 4 -t 0:20:00 -l logs/nmixture_si_truth.log -- bin/Rscript dev/nmixture/si_truth.R
bin/ejob -c 16 -m 32 -t 3:00:00 -l logs/nmixture_si_recover.log -- bin/Rscript core/run.R recover nmixture --config modules/nmixture/example/recovery.yaml --nsim 200 --cores 16 --out results/nmixture/si/recovery
bin/ejob -c 4 -m 16 -t 2:00:00 -l logs/nmixture_si_fit.log -- bin/Rscript core/run.R fit nmixture --config modules/nmixture/example/config.yaml --out results/nmixture/si/fit

Step-by-step walkthrough: modules/nmixture/example/WALKTHROUGH.md.

S1.10 distance: Mule deer density on a winter range from line-transect distance sampling, with shrub-cover and road-distance effects (simulated data)

OPEN All estimands that answer the question (density, abundance, per-unit density, detection, habitat rate ratio) pass planted truth at this design, the fit converged and both GOF checks pass; but the HDS density at an average transect (whose uncertainty also sets the bands on the conditional-effect curves and the map) under-covers (0.885 vs bound 0.904; VALIDATION.md Known limit 9, prepared fix not installed).

Question

How many mule deer are on the 360 km2 winter range (per km2 and in total, overall and in the north and south units), and is deer density higher where shrub cover is higher and farther from the highway?

What the researcher typed

I have a mule deer line-transect survey in modules/distance/example: a file of deer groups seen with the perpendicular distance from the line (m), group size, observer and date; a file of the 36 transects with their length, management unit (north or south, 180 km2 each), shrub cover and distance to the highway; and a 1 km2 grid of the whole 360 km2 winter range with the same two habitat values. Please estimate how many deer there are per km2 and in total, for each unit and overall, and tell me whether deer are more common where shrub cover is higher and farther from the road. Check that the method works for a survey like mine before you trust it, and give me a report and figures I can put in a paper.

Study design

True values used to simulate the example data

Planted-truth recovery at this design

Estimates from the example data

Convergence and goodness of fit

Interpretation from the report

A group within 400 m of the line had probability 0.44 (95% CI 0.38 to 0.50) of being detected. Equivalently, the effective strip half-width is 175 m.

Density: 12.0 (95% CI 9.24 to 15.5) individuals per km2, which over the 360 km2 study area is 4,314 (95% CI 3,328 to 5,593) individuals.

This is 3.59 (95% CI 2.82 to 4.58) groups per km2 times a mean group size of 3.34 (95% CI 3.12 to 3.56).

By stratum: north 7.97 (95% CI 5.68 to 11.2) per km2; south 16.0 (95% CI 11.6 to 22.0) per km2.

Each increase of 10 percentage points in shrub cover multiplies expected density by 1.89 (95% CI 1.56 to 2.30), i.e. +89%.

Each increase of 3 km in distance from road multiplies expected density by 1.40 (95% CI 1.12 to 1.74), i.e. +40%.

This is an association from an observational survey, not proof that shrub cover causes the difference.

The hierarchical model's density over the study area (12.2 (95% CI 9.70 to 15.3) per km2) agrees with the design-based estimate above.

Goodness of fit: Cramer-von Mises test of the distances, p = 0.72. No evidence of lack of fit.

Figures

Limits that matter at this design

• HDS intercept (density_at_average_site) ignores detection-key selection uncertainty: coverage 0.885 at this design (0.92 when half-normal chosen, as in this data set; 0.81 otherwise). Design-based density, abundance and slopes unaffected (Known limit 9).
• Design-based density intervals are conservative when density follows habitat across transects (coverage 0.975 here; Known limit 3).
• g(0) = 1, no evasive movement and accurate distances cannot be tested from these data; 75 m evasive movement biased density -27% in simulation, heaping at zero biased it upward (Known limit 4).
• Recovery simulated one habitat covariate and a single detection scale; the second covariate (road) and the observer difference were checked only on the example fit.
• The density map extrapolates to grid cells with shrub cover outside the surveyed 13.9-45.4 %.

recovery meta

Commands to reproduce

bin/Rscript -e 'source("core/R/env.R"); source("dev/distance/si/example_truth.R")'
bin/ejob -c 16 -m 32 -t 3:00:00 -l logs/distance_si_recover.log -- bin/Rscript core/run.R recover distance --config modules/distance/example/recovery.yaml --nsim 200 --cores 16 --out results/distance/si/recovery
bin/ejob -c 4 -m 16 -t 2:00:00 -l logs/distance_si_fit.log -- bin/Rscript core/run.R fit distance --config modules/distance/example/config.yaml --out results/distance/si/fit
python3 dev/distance/si/make_si_summary.py

Step-by-step walkthrough: modules/distance/example/WALKTHROUGH.md.

Spatial capture–recapture

S1.11 scr: Bobcat density from camera traps: closed spatial capture-recapture with a canopy-cover density covariate (simulated data)

DONE

Question

How many bobcats live in the study area (density per km2 and abundance), and is bobcat density higher where forest canopy cover is higher, allowing for cameras on trails detecting bobcats more often and for cameras that stopped working?

What the researcher typed

I have bobcat camera-trap data in modules/scr/example: a camera file with UTM coordinates, whether each camera was on a trail and which weeks it worked; a file of photo-detections of individually identified bobcats by week and camera; and a 400 m habitat grid with canopy cover and a flag for the reservoir. Please estimate bobcat density and how many bobcats live in the area, and tell me whether density is higher where canopy cover is higher. Check that the method works for a design like mine before you trust it, and give me a report and figures I can put in a paper.

Study design

True values used to simulate the example data

Planted-truth recovery at this design

Estimates from the example data

Convergence and goodness of fit

Interpretation from the report

Density was estimated at 0.24 (95% CrI 0.16 to 0.34) animals per km², i.e. 23.7 per 100 km2.

Over the 350 km² of habitat in the state space (the detector array plus a buffer), that is 82.7 (95% CrI 58.0 to 118) animals, of whom 33 (40%) were detected.

The spatial scale of detection, sigma, was 1326 (95% CrI 1165 to 1510) m; a home-range area of about 33 (95% CrI 26 to 43) km². Real home ranges are rarely circular, so treat this as an index of ranging scale.

'Camera on trail' (1 vs 0) multiplied the odds of detection by 3.2 (95% interval 1.9 to 5.4). Baseline detection was 0.30 (0.21 to 0.41) where Camera on trail = 1 versus 0.12 where it = 0. This is a detection effect (how easily animals are recorded), not a difference in where animals live.

Each 1 SD (15 units) higher Canopy cover (%) is associated with 51% higher density (rate ratio 1.5, 95% interval 1.1 to 2.1). Posterior probability that the association is positive: 0.99.

Expected density is 0.13 per km² (0.070 to 0.20) where Canopy cover (%) = 34.9 (10th percentile of the habitat) and 0.37 (0.19 to 0.62) where it = 72.6 (90th percentile).

This is a correlation between density and habitat across the state space, not proof that the habitat causes higher density.

Figures

Limits that matter at this design

• Small-sample upward bias: posterior means of density/abundance lean high with few detected animals (+12% with ~18 detected, +4-8% with ~42; VALIDATION.md limit 2); with 33 detected, quote the interval and compare with secr.
• Density-covariate false positives: with ~33-42 detected animals and no true effect ~10% of 95% intervals excluded 0 (VALIDATION.md limit 3); the canopy interval's lower end (0.06) is close to 0, so the size of the association is uncertain.
• Buffer must be >= 3-4 sigma (1 sigma inflated density +47%); here 3.8 sigma.
• The camera array sits in below-average canopy; density over the whole state space relies on the canopy model to extrapolate (unmodelled heterogeneity biased density +64% in VALIDATION.md 3b).
• Model scope: half-normal detection, Euclidean distance, closed population, correct photo IDs, no sex-specific ranging (sex not recorded).
• Density covariate effects are associations, not causal effects.

recovery meta

estimates secr

Commands to reproduce

bin/ejob -c 1 -m 8 -t 0:30:00 -- bin/Rscript dev/scr/make_example.R
bin/ejob -c 16 -m 32 -t 3:00:00 -l logs/scr_si_recover.log -- bin/Rscript core/run.R recover scr --config modules/scr/example/recovery.yaml --nsim 100 --cores 16 --out results/scr/si/recovery
bin/ejob -c 16 -m 32 -t 3:00:00 -l logs/scr_si_recover_matched.log -- bin/Rscript core/run.R recover scr --config modules/scr/example/recovery_matched_n.yaml --nsim 100 --cores 16 --out results/scr/si/recovery_matched_n
bin/ejob -c 4 -m 16 -t 2:00:00 -l logs/scr_si_fit.log -- bin/Rscript core/run.R fit scr --config modules/scr/example/config.yaml --out results/scr/si/fit
bin/ejob -c 1 -m 8 -t 0:30:00 -- bin/Rscript core/run.R figures scr --fit results/scr/si/fit/scr_fit.rds --config modules/scr/example/config.yaml
../.venv/bin/python dev/scr/make_si_summary.py

Step-by-step walkthrough: modules/scr/example/WALKTHROUGH.md.

S1.12 scr-smr: Fisher density from ear tags and trail cameras: generalized spatial mark-resight (simulated data)

DONE

Question

How many fishers live in the study area (density per km2 and abundance), does density rise with canopy cover, and do trail cameras photograph fishers more often, given that tagged fishers were caught only at box traps in the middle of the camera grid?

What the researcher typed

I've put my fisher data in a folder: box-trap captures and trap locations from May, trail camera locations and a photo log from June to August, GPS fixes from 9 collared fishers, and a habitat grid with canopy cover. Some photos show a readable ear tag, some a tag we couldn't read, and most no tag. Please estimate fisher density and abundance, tell me whether density rises with canopy cover, and whether trail cameras get more photos. The box traps were only in the middle of the camera grid. Please check the method works for a design like mine before you trust it, and give me a report and figures I can use in a paper.

Study design

True values used to simulate the example data

Planted-truth recovery at this design

Estimates from the example data

Convergence and goodness of fit

Interpretation from the report

Density was estimated at 0.30 (95% CrI 0.24 to 0.38) fishers per km² (30.1 per 100 km2).

Over the 289 km² state space (camera array plus buffer) that is 87.1 (95% CrI 68.0 to 109) fishers, of which the 23 marked animals are a fraction 0.27 (95% CrI 0.21 to 0.34).

Sigma, the spatial scale of space use, was 1010 (95% CrI 966 to 1057) m: circular home range of 19 (95% CrI 18 to 21) km².

Each 1 SD (12 units) higher canopy cover is associated with 64% higher density (rate ratio 1.6, 95% interval 1.2 to 2.2); expected density 0.15 per km² (0.082 to 0.24) at 36% canopy and 0.59 (0.36 to 0.87) at 69% canopy.

This is an association between density and habitat across the state space, not proof that the habitat causes higher density.

A camera on a trail multiplied the sighting rate by 2.3 (95% interval 1.8 to 3.0). This is a detection effect, not a difference in where animals live.

Because marking was modelled, marked animals are not assumed to be a random sample: animals living near the marking traps are more likely to carry marks, and the model accounts for that (generalized SMR).

Figures

Limits that matter at this design

• Conventional SMR after spatially concentrated marking is badly biased (-31% N, coverage 0.10); the box-trap captures must be modelled (gSMR).
• Upward small-sample bias of +4-7% in N/density (posterior skew); +4.6% at this design, inside the gate.
• secr's joint gSMR fit was ~26% low in validation; use secr only to check sigma and detection rates.
• Not modelled: tag loss, misread tag status, individual heterogeneity in photo rate (e.g. collar effects), non-independent photo events, open populations.
• GPS fixes are assumed independent; autocorrelated data must be thinned.
• With <= 5 marked animals N is biased ~+15% and the gate fails; this design has 23.
• Only the Poisson resight model was validated at nsim = 100.
• Power to detect the canopy effect (beta = 0.5 per SD) is 0.77 at this design.

recovery meta

Commands to reproduce

bin/ejob -c 1 -m 8 -t 0:30:00 -l logs/scr-smr_make_example.log -- bin/Rscript dev/scr-smr/make_example.R
bin/ejob -c 16 -m 32 -t 3:00:00 -l logs/scr-smr_si_recover.log -- bin/Rscript core/run.R recover scr-smr --config modules/scr-smr/example/recovery.yaml --nsim 100 --cores 16 --out results/scr-smr/si/recovery
bin/ejob -c 4 -m 16 -t 2:00:00 -l logs/scr-smr_si_fit.log -- bin/Rscript core/run.R fit scr-smr --config modules/scr-smr/example/config.yaml --out results/scr-smr/si/fit
bin/ejob -c 2 -m 8 -t 0:30:00 -l logs/scr-smr_si_figures.log -- bin/Rscript core/run.R figures scr-smr --config modules/scr-smr/example/config.yaml --out results/scr-smr/si/fit

Step-by-step walkthrough: modules/scr-smr/example/WALKTHROUGH.md.

S1.13 scr-search: Fisher density from a scat-detection-dog survey with GPS-logged search routes (search-encounter SCR, simulated data with known truth)

DONE

Question

How many fishers live in the study area and at what density, is density higher where there is more forest, and do experienced dog teams find more scats per km searched?

What the researcher typed

I ran a scat-detection-dog survey for fishers in the western Oregon Cascades in summer 2025: 36 search routes in three monthly sessions, each logged by GPS (search_tracks_gps.csv), with a field sheet per route saying which dog team searched and whether it was one of our experienced teams (track_field_sheet.csv). The scats and their genotypes are in fisher_scats.csv (blank fisher_id = failed genotype). I also have forest cover on a 700 m grid (habitat_grid.csv). Please estimate fisher density and abundance, tell me whether density is related to forest cover, and whether experienced teams find more scats per km. I'd like a report and figures for a paper.

Study design

True values used to simulate the example data

Planted-truth recovery at this design

Estimates from the example data

Convergence and goodness of fit

Interpretation from the report

Density was estimated at 0.070 (95% CrI 0.045 to 0.10) animals per km², i.e. 7.04 per 100 km2. Over the 784 km² state space (the searched area plus a buffer), that is 54.7 (95% CrI 38.0 to 76.0) animals, of which 22 (40%) were identified from samples.

The spatial scale of space use, sigma, was 1561 (95% CrI 1345 to 1828) m. If space use is a circular bivariate normal around each activity centre, 95% of an animal's activity falls within 3821 m of its centre, a home-range area of about 46 (95% CrI 34 to 63) km²; treat this as an index of ranging scale.

Encounter rate: a searcher covering 1 km right at an animal's activity centre finds on average 0.18 (95% CrI 0.11 to 0.27) of its samples (at the reference level of the detection covariates); a long straight search route passing through the activity centre would yield 0.70 (95% CrI 0.47 to 1.0) samples of that animal.

Detection: 'Experienced dog team' (1 vs 0) multiplied the number of samples found per km searched by 1.65 (95% interval 1.04 to 2.58); per-km encounter rate at the activity centre 0.30 with Experienced dog team = 1 versus 0.18 with Experienced dog team = 0. This is a detection (search efficiency) effect, not a difference in where animals live.

Density and Forest cover (%): each 1 SD (12 units) higher Forest cover (%) is associated with 155% higher density (rate ratio 2.55, 95% interval 1.60 to 4.20). Expected density is 0.015 per km² (0.0036 to 0.037) where Forest cover (%) = 39.6 (10th percentile of the state space) and 0.16 (0.10 to 0.23) where it = 71.6 (90th percentile).

This is a correlation between density and habitat across the state space, not proof that the habitat causes higher density. Because search effort was modelled explicitly (km searched per cell), more samples in heavily searched areas are not mistaken for higher density.

Figures

Limits that matter at this design

• Sparse data: with 20-35 animals identified (here 22) and density CV ~0.2, the posterior MEAN of density and N ran 9-13% high in independent simulations (coverage nominal); here recovery bias was +4.8%. Report the posterior median (0.0697, N 54) and the secr ML estimate alongside the mean.
• Abundance refers to the whole 784 km2 state space; away from the routes density is extrapolated from the forest model (detection_map). Quote density rather than N alone.
• Search effort must be recorded: ignoring it more than doubled the apparent forest effect in simulations with road-following routes through forest.
• 1 km cells (0.64 sigma) inflate sigma by about 3% relative to exact locations.
• Power at ~250 km of searching: 0.57 for the forest effect (rate ratio 1.65 per SD) and 0.66 for the team effect; ~400 km needed for ~0.75 power on forest.
• Untested assumptions: closure over June-August, one sigma for both sexes, independent scats (no latrines), correct genotypes.

recovery notes

nsim = 100 (one fit ~4 min); 100/100 ran and converged (max R-hat 1.008); gate: |bias| <= max(10%, 3 MC SE), coverage >= 0.885; abundance truth is the mean over replicates (78.6); N-hat/N median 1.02 (5-95%: 0.71-1.42).

estimates secr cells

Commands to reproduce

bin/ejob -c 16 -m 32 -t 3:00:00 -l logs/scr-search_si_recover.log -- bin/Rscript core/run.R recover scr-search --config modules/scr-search/example/recovery.yaml --nsim 100 --cores 16 --out results/scr-search/si/recovery
bin/ejob -c 4 -m 16 -t 2:00:00 -l logs/scr-search_si_fit.log -- bin/Rscript core/run.R fit scr-search --config modules/scr-search/example/config.yaml --out results/scr-search/si/fit
bin/ejob -c 2 -m 8 -t 0:20:00 -- bin/Rscript core/run.R figures scr-search --fit results/scr-search/si/fit/scr-search_fit.rds --config modules/scr-search/example/config.yaml

Step-by-step walkthrough: modules/scr-search/example/WALKTHROUGH.md.

S1.14 scr-open: Pacific marten abundance, survival, recruitment and canopy association from five years of DNA hair snares (open-population spatial capture-recapture, simulated data with known truth)

DONE Planted truth 28/28 pass at this design; R-hat <= 1.005, ESS >= 1035; all posterior predictive p-values 0.39-0.64; NumPyro agrees with Nimble; module negative controls fail as required (VALIDATION.md); reproducible commands. Caveat: total abundance depends on the canopy model because the snares sample thin canopy.

Question

How many Pacific martens live in the study area each summer (2020-2024), is the population declining, what are annual apparent survival and recruitment, and is marten density higher under denser canopy?

What the researcher typed

I have five summers (2020-2024) of marten hair-snare data from 49 snares on a 1-km grid, each checked 5 times a summer. hair_snares.csv has the snare locations (UTM) and which checks each snare was working; marten_genotypes.csv has one row per genotyped sample with year, check, snare and marten ID; and habitat_grid.csv has canopy cover on a 400 m grid around the snares. I want to know how many martens there are each year, whether they are declining, annual survival and recruitment, and whether density is higher under denser canopy. Please give me a report with figures I can use in a paper.

Study design

True values used to simulate the example data

Planted-truth recovery at this design

Estimates from the example data

Convergence and goodness of fit

Interpretation from the report

Abundance in the 144 km² state space (traps plus buffer) was 119 (95% CrI 85.0 to 161) in period 2020 and 186 (95% CrI 144 to 238) in period 2024 (density 0.83 and 1.3 animals per km²).

Over the study the population changed by a factor of 1.58 (95% CrI 1.15 to 2.16); posterior probability of a decline from the first to the last period: < 0.01.

Apparent survival between consecutive periods was 0.84 (95% CrI 0.77 to 0.91): an animal alive in one period had an estimated 84% chance of still being alive and resident in the next. 'Apparent' because death and permanent emigration from the study area cannot be told apart.

Per-capita recruitment rates (recruits per animal alive in the previous period) were 2021: 0.36 (95% CrI 0.14 to 0.68); 2022: 0.23 (95% CrI 0.085 to 0.46); 2023: 0.23 (95% CrI 0.088 to 0.44); 2024: 0.31 (95% CrI 0.15 to 0.52).

In total 277 (95% CrI 222 to 341) different animals were alive in the state space at some point during the study (superpopulation), of whom 95 (34%) were ever detected.

A trap at an animal's activity centre detects it on a given occasion with probability 0.32 (95% CrI 0.28 to 0.35); the spatial scale sigma was 803 (95% CrI 764 to 846) m, so 95% of activity falls within 1966 m of the centre.

Each 1 SD (12) higher Canopy cover (%) is associated with 273% higher density (rate ratio 3.7, 95% CrI 2.9 to 5.0). This describes where animals' activity centres are, not proof that Canopy cover (%) causes higher density.

WARNING: the traps sample habitat of unrepresentative density (state-space average 1.94 x the trapped area), so total abundance rests on extrapolating the density covariate effect to habitat the traps did not cover.

Figures

Limits that matter at this design

• Snares sample below-average canopy: 72% of martens are estimated to live > 2 sigma from any snare and state-space density is 1.94 x trapped-area density, so total abundance rests on extrapolating the canopy effect; a homogeneous-density fit (openCR) gives ~1.7-1.9x lower abundance.
• In this realisation the canopy-effect interval (1.06-1.61) misses the simulation value 1.0 (true-centre oracle 1.18); the 100-replicate recovery shows no bias at this design (+1.0%, coverage 0.98).
• Static activity centres: range shifts between years bias survival low.
• Constant survival, p0 and sigma; sex/age differences, trap response and time-varying detection are not modelled (PPCs flag but were not validated against such violations).
• Apparent survival: death and permanent emigration are confounded; recruitment includes immigration.
• First-year recruitment is not identifiable; last-interval survival/recruitment are less precise; equal intervals between years assumed.
• With ~40 detected animals abundance is 4-8% high and per-capita recruitment 12-20% high; this design (95 animals) recovers within ~5%.
• Recovery validated with Nimble marginal only; NumPyro agrees on this dataset.

recovery note

nsim = 100 (one fit takes ~7 min), module fast_settings; truth for realised quantities is the mean over replicates; 100/100 converged; gate 28/28 PASS

Commands to reproduce

bin/ejob -c 1 -m 4 -t 0:20:00 -l logs/scr-open_make_example.log -- bin/Rscript dev/scr-open/make_example.R
bin/ejob -c 16 -m 32 -t 3:00:00 -l logs/scr-open_si_recover.log -- bin/Rscript core/run.R recover scr-open --config modules/scr-open/example/recovery.yaml --nsim 100 --cores 16 --out results/scr-open/si/recovery
bin/ejob -c 4 -m 16 -t 2:00:00 -l logs/scr-open_si_fit.log -- bin/Rscript core/run.R fit scr-open --config modules/scr-open/example/config.yaml --out results/scr-open/si/fit
bin/ejob -c 1 -m 8 -t 0:30:00 -l logs/scr-open_si_figures.log -- bin/Rscript core/run.R figures scr-open --config modules/scr-open/example/config.yaml --fit results/scr-open/si/fit/scr-open_fit.rds
python3 dev/scr-open/make_si_summary.py

Step-by-step walkthrough: modules/scr-open/example/WALKTHROUGH.md.

Demography and population dynamics

S1.15 cjs: Adult song sparrow apparent survival from constant-effort banding (Cormack-Jolly-Seber, simulated data with known truth)

DONE

Question

What is annual apparent survival of adult song sparrows at a banding station, does it differ between sexes, and is it associated with winter cold and body mass, allowing for transient birds that never return?

What the researcher typed

I have nine years of song sparrow banding data from my station (2016 to 2024) in three files: sosp_captures.csv (every capture, with band number and date), sosp_birds.csv (sex and body mass of each bird when first banded) and sosp_years.csv (mean January minimum temperature for each winter). I want annual adult survival, whether it differs between males and females, and whether it is related to winter cold and body mass. Many birds are only caught once; I think some are migrants passing through. Please analyse this and give me a report with figures I can use in a paper.

Study design

True values used to simulate the example data

Planted-truth recovery at this design

Estimates from the example data

Convergence and goodness of fit

Interpretation from the report

Apparent annual survival at average mean January minimum temperature and body mass of residents: Female 0.54 (95% CI 0.43-0.65); Male 0.62 (95% CI 0.53-0.70).

Apparent survival is the probability of surviving AND staying in the study area; animals that emigrate permanently are indistinguishable from deaths, so true survival is at least this high.

Sex: Male vs Female (apparent survival): odds ratio 1.34 (95% CI 0.76-2.40). The interval includes 1: the data do not show a clear difference.

Mean January minimum temperature: each 1 SD (2.04 deg C) increase multiplies the odds of apparent survival by 1.30 (95% CI 0.91-1.87); female survival 0.45 at -6.0 deg C vs 0.61 at -1.0 deg C. The interval includes no effect. This is an association, not proof of a causal effect.

Body mass: each 1 SD (1.66 g) increase multiplies the odds of apparent survival by 1.17 (95% CI 0.90-1.52); female survival 0.49 at 19.8 g vs 0.59 at 23.9 g. The interval includes no effect.

Transients: an estimated 0.26 (95% CI 0.00-0.53) (Female); 0.32 (95% CI 0.11-0.48) (Male) of newly marked animals were transients.

Detection (recapture) probability per occasion: Female 0.41 (95% CI 0.28-0.55); Male 0.57 (95% CI 0.47-0.65).

Goodness of fit: overall c-hat = 0.92 (p = 0.624). Test 3.SR detected transients (p = 0.027): the fitted model accounts for them (tsm). Excluding that component, c-hat = 0.50 (p = 0.990). No evidence of remaining lack of fit.

Figures

Limits that matter at this design

• Low power at this design: 17% for the sex effect on survival, 51-56% for the winter and body mass effects (200-replicate recovery); intervals that include zero are inconclusive, not evidence of no effect.
• Transients must be modelled: ignoring them biases survival about 13% low with coverage about 0.5 (validated negative control); 3.SR flagged transience here and tsm was fitted.
• proportion_transient is an unconstrained ratio; the female interval extends below 0 (-0.12) and is truncated to 0 in the report text.
• False-positive rate for individual covariates such as body mass is 7.5-8% at nominal 5% (Wald tests).
• The winter effect rests on 8 winters (one value per year) and is an association, not a causal effect.
• Not modelled: individual heterogeneity in recapture beyond sex, band loss, temporary emigration; apparent survival = true survival x site fidelity.

recovery meta

Commands to reproduce

bin/ejob -c 1 -m 4 -t 0:10:00 -l logs/cjs_make_example.log -- bin/Rscript dev/cjs/make_example.R
bin/ejob -c 4 -m 8 -t 0:30:00 -l logs/cjs_si_selftest.log -- bin/Rscript core/run.R selftest cjs --out dev/cjs/si_selftest
bin/ejob -c 16 -m 32 -t 3:00:00 -l logs/cjs_si_recover.log -- bin/Rscript core/run.R recover cjs --config modules/cjs/example/recovery.yaml --nsim 200 --cores 16 --out results/cjs/si/recovery
bin/ejob -c 4 -m 16 -t 2:00:00 -l logs/cjs_si_fit.log -- bin/Rscript core/run.R fit cjs --config modules/cjs/example/config.yaml --out results/cjs/si/fit
python3 dev/cjs/make_si_summary.py

Step-by-step walkthrough: modules/cjs/example/WALKTHROUGH.md.

S1.16 multistate: Chytrid (Bd) infection, survival and clearance in a PIT-tagged frog population (multievent capture-recapture, simulated data with known truth)

DONE

Question

Do infected frogs survive less well than uninfected frogs, how often do frogs become infected or clear Bd infection each year, is water temperature associated with infection and clearance, and do larger frogs survive better, when some swabs were lost or inconclusive?

What the researcher typed

I have eight years (2015-2022) of capture data for a PIT-tagged frog population in three files: frog_captures.csv (each capture with the frog's id, year and Bd swab result: uninfected, infected, or U when the swab was lost or inconclusive), frog_individuals.csv (snout-vent length of each frog) and survey_years.csv (mean water temperature and survey effort in person-hours for each year). I want to know whether infected frogs survive less well, how often frogs become infected or clear the infection each year, whether water temperature is related to infection and clearance, and whether bigger frogs survive better. Please analyse this and give me a report with figures I can use in a paper.

Study design

True values used to simulate the example data

Planted-truth recovery at this design

Estimates from the example data

Convergence and goodness of fit

Interpretation from the report

Survival over one year for an animal with average covariate values was 0.84 (95% credible interval 0.75-0.96) in state 'uninfected'; 0.50 (95% credible interval 0.37-0.63) in state 'infected'.

The difference in survival (uninfected minus infected) was 0.34 (95% interval 0.18 to 0.50); posterior probability that survival is higher in 'uninfected' = 1.00.

Transitions are probabilities of being in a new state at the next occasion, for animals that survive the year (at average transition covariates): from 'uninfected' to 'infected' 0.33 (95% credible interval 0.19-0.52); from 'infected' to 'uninfected' 0.34 (95% credible interval 0.18-0.54).

Expected number of years a surviving animal stays in a state before switching (1 / switching probability, median and 95% interval): 'uninfected' 3.1 (1.9-5.4); 'infected' 3.0 (1.8-5.5).

Survival and svl_mm: odds ratio per 1 SD (= 6) 1.77 (95% interval 1.34-2.50); P(effect > 0) = 1.00. Going from the 10th to the 90th percentile of svl_mm (44.5 to 59), predicted survival changes in 'uninfected' from 0.72 to 0.91, in 'infected' from 0.33 to 0.65.

Transitions and water_temp_C (from 14.6 to 19.1, 10th to 90th percentile): 'uninfected' to 'infected': 0.48 at the 10th percentile vs 0.23 at the 90th (difference -0.24, 95% interval -0.41 to -0.06); 'infected' to 'uninfected': 0.16 at the 10th percentile vs 0.56 at the 90th (difference 0.40, 95% interval 0.10 to 0.69). Associations only.

Detection (the chance that a live animal is caught or resighted at an occasion, at average detection covariates) was 0.61 (95% credible interval 0.49-0.75) in 'uninfected'; 0.39 (95% credible interval 0.23-0.63) in 'infected'. Low detection in a state means fewer data on that state, not fewer animals.

State uncertainty (multievent model): when a detected animal was in a given state, the probability that its state was recorded (rather than coded unknown) was 0.74 (95% credible interval 0.67-0.81) in 'uninfected'; 0.86 (95% credible interval 0.73-0.98) in 'infected'.

Goodness of fit: posterior predictive checks showed no clear misfit (Bayesian p = 0.43 sightings per occasion, 0.64 state switches, 0.60 animals never seen again).

Figures

Limits that matter at this design

• Transition slopes on a yearly covariate rest on only 7 temperature values: power at this design was 0.76 for infection and 0.26 for clearance, and detected slopes tend to be inflated (example clearance slope 0.92 vs true 0.50). Associations only; confounded with any other year effect.
• False-positive rate for transition slopes was slightly above nominal (0.10 for one slope in the null validation run).
• Misclassification (e.g. false-negative qPCR recorded as uninfected) is not modelled; only lost/inconclusive swabs are treated as uncertain.
• Assignment (swab success) is assumed constant over time and the same at first capture and recaptures.
• Survival is apparent survival (death and permanent emigration are confounded).
• U-CARE tests are approximate for multievent data (unknown sightings treated as not seen).
• Covariate slopes on survival and detection are shared by both states; body length is fixed at first capture.
• The stationary distribution describes switching only and is not a prevalence forecast.

Commands to reproduce

bin/ejob -c 1 -m 4 -- bin/Rscript dev/multistate/make_example.R
bin/ejob -c 16 -m 32 -t 3:00:00 -l logs/multistate_si_recover.log -- bin/Rscript core/run.R recover multistate --config modules/multistate/example/recovery.yaml --nsim 100 --cores 16 --out results/multistate/si/recovery
bin/ejob -c 4 -m 16 -t 2:00:00 -l logs/multistate_si_fit.log -- bin/Rscript core/run.R fit multistate --config modules/multistate/example/config.yaml --out results/multistate/si/fit
bin/ejob -c 1 -m 8 -t 0:30:00 -- bin/Rscript core/run.R figures multistate --fit results/multistate/si/fit/multistate_fit.rds --config modules/multistate/example/config.yaml --out results/multistate/si/fit

Step-by-step walkthrough: modules/multistate/example/WALKTHROUGH.md.

S1.17 jolly-seber: Western pond turtle abundance, recruitment and survival from eight years of PIT-tag hoop-trapping (Jolly-Seber, simulated data with known truth)

DONE Every interpreted estimand passes planted truth at this design (42/45 overall; N_t rel. bias -2.5% to +2.8%, coverage 0.925-0.965; survival -0.6% / +0.1%, coverage 0.955 / 0.95); fit converged (R-hat 1.001, min ESS 2583); PPC p 0.52 / 0.62 and all R2ucare p > 0.5; module cross-checked against marked and latent Nimble, negative controls in VALIDATION.md. First-interval recruitment (entrants_2, entry_prob_2, recruitment_rate_1) fails the bias criterion and is excluded from interpretation.

Question

How many western pond turtles used the pond complex each year and over the whole study, how many new turtles joined each year, what share survive from year to year (and do males survive less well than females), and is the population declining?

What the researcher typed

I have eight summers (2015 to 2022) of hoop-trapping data for western pond turtles at one pond complex, in two files: turtle_captures.csv (every capture, with the turtle's ID, the season and its sex; U means we could not sex it) and turtle_seasons.csv (trap-nights of effort in each season). I want to know how many turtles were there each year and in total, how many new turtles joined each year, what share survive from year to year and whether males survive less well than females, and whether the population is declining. Please analyse this and give me a report with figures I can use in a paper.

Study design

True values used to simulate the example data

Planted-truth recovery at this design

Estimates from the example data

Convergence and goodness of fit

Interpretation from the report

175 western pond turtles were caught at least once over 8 occasions (2015 to 2022); 87 (50%) were caught on more than one occasion, and 40.2 were caught per occasion on average.

Estimated abundance was 108 (95% CI 76 to 148) at the first occasion (2015) and 112 (95% CI 87 to 145) at the last (2022). Across the study the average growth rate was 1.006 (95% CI 0.960 to 1.055) per year, i.e. a change of +0.6% per year; the posterior probability that the population declined (growth rate < 1) is 0.39.

In total 268 (95% CI 236 to 310) western pond turtles are estimated to have used the study population at some time during the study (the superpopulation), of which 175 (65%) were ever caught.

Apparent survival per year was 0.85 (95% CI 0.78 to 0.90) for female and 0.68 (95% CI 0.58 to 0.78) for male. The odds of surviving were 0.39 times as high for male as for female (95% CI 0.21 to 0.71); posterior probability that male animals survive less well than female animals > 0.99.

'Apparent' survival means that death and permanent emigration cannot be told apart.

Capture probability of an animal present was 0.37 (95% CI 0.31 to 0.44) at average sampling effort. Capture probability rose from 0.29 to 0.47 between trapping effort (trap-nights) = 90 and 150 (one SD below to one SD above the mean; odds ratio per SD 1.46, 95% CI 1.15 to 1.85). This is a sampling (detection) effect, not a change in the population.

Agent caveat: the report's average per-capita recruitment (0.23, 95% CI 0.16-0.30) includes the 2015-2016 interval, which failed planted truth; read recruitment from 2016-17 onward (11-35 entrants per year).

Figures

Limits that matter at this design

• First-interval recruitment is weakly identified: 2016 entrants and entry probability +26%, 2015-16 per-capita recruitment +47% in planted truth (coverage 0.97-0.975); not interpreted (VALIDATION Known limit 3).
• Low power for the sex effect on survival at this design: 0.27 for logit -0.4 (200-replicate recovery). The example estimate (-0.93, 95% CI -1.54 to -0.34) is about twice the truth, so quote the interval.
• Temporary emigration is not modelled: in the negative control it inflated N_super +22% (coverage 0.43) and biased capture probability -33%; a robust design is needed to test it (Known limit 4).
• No individual capture heterogeneity or trap-dependence model: heterogeneity biased N about 20% low in the verifier negative control; here the PPC and R2ucare 2.CT did not flag it (Known limit 5).
• Survival is apparent (death plus permanent emigration); abundance counts turtles using the trapped ponds, not a density.
• The marked cross-check cannot use unknown-sex turtles (sex effect dropped, survival pooled) and reports expected rather than realised N_t (Known limit 6).

recovery meta