Dias et al. (2025): how long the hierarchical models take¶
Dataset DRYAD:bvq83bkk9, version id 333377 (the only
version; see the header of reproduce_dias_2025.ipynb for the files and their SHA-256). Paper
doi:10.1098/rspb.2024.2556.
Tables 1 to 3, S1 to S3 and S9 to S11 rest on the hierarchical model of lm_mod.stan (the paper's
code, Zenodo doi:10.5281/zenodo.14226724), which Script01 runs with NUTS, adapt_delta 0.85 and 4
chains of 4000 iterations (2000 warm-up). The claims require each model under three filtering
readings ((a), (c) and (d); (b) removes no specimen and equals (a)). This notebook runs each main
model briefly, 4 chains of 100 warm-up and 100 kept iterations with the code's seed, and scales the
time to the full 4000 iterations, to decide which fits fit in a notebook (CI stops a notebook at 30
minutes). It also shows the port of the Stan file to current syntax.
import tempfile
import time
import warnings
from pathlib import Path
import numpy as np
import pandas as pd
import wingallom as w
warnings.filterwarnings("ignore", module="rdata")
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
import logging
started = time.time()
for name, digest in w.check_staged().items():
print(f"{name}: SHA-256 {digest} (matches Dryad's file list, version {w.DRYAD_VERSION} of {w.DRYAD_ID})")
t0 = time.time()
print("CmdStan", w.CMDSTAN_VERSION, "at", w.install_cmdstan(Path(tempfile.gettempdir()) / "cmdstan"),
f"({time.time() - t0:.0f} s to install or find)")
logging.getLogger("cmdstanpy").setLevel(logging.ERROR) # sampler messages; divergences are counted below
code_zip = w.code_zip()
original = code_zip.read("Code/lm_mod.stan").decode()
model_code = w.port_stan(original)
workdir = Path(tempfile.mkdtemp())
t0 = time.time()
model = w.compile_stan(model_code, workdir)
print(f"lm_mod.stan compiled in {time.time() - t0:.0f} s")
Data.zip: SHA-256 2185ac50e9a59c3a1dff59473216de2ab001d0d0044b1fa15442ae18a520d554 (matches Dryad's file list, version 333377 of bvq83bkk9) README.md: SHA-256 6c19de67b1c6c52802be42a6593c175a02bee0dab0ca2f6f958f08a0a3fa0482 (matches Dryad's file list, version 333377 of bvq83bkk9)
/home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html from .autonotebook import tqdm as notebook_tqdm
CmdStan install directory: /tmp/cmdstan
Installing CmdStan version: 2.40.0 Downloading CmdStan version 2.40.0 Download successful, file: /tmp/tmpn7i_0ge4 Extracting distribution
Unpacked download as cmdstan-2.40.0 Building version cmdstan-2.40.0, may take several minutes, depending on your system.
Installed cmdstan-2.40.0 Test model compilation
CmdStan 2.40.0 at /tmp/cmdstan/cmdstan-2.40.0 (131 s to install or find)
lm_mod.stan compiled in 23 s
Data check¶
import difflib
print("Changes made to lm_mod.stan (array syntax of Stan 2.33 and later; nothing else):")
for line in difflib.unified_diff(original.splitlines(), model_code.splitlines(), lineterm="", n=0):
print(line)
morph = w.morph_prepared("a")
data = {}
for cls in w.CLASSES:
p = w.prepare(morph, cls, "main")
d = w.stan_data(p, w.MAIN_PREDICTORS[cls])
data[cls] = d
print(f"\n{cls}: N = {d['N']} birds, J = {d['J']} species-years, K = {d['K']} species, E = {d['E']} columns of G "
f"(intercept and {', '.join(w.MAIN_PREDICTORS[cls])})")
print(f" scaled wing {d['WL'].min():.2f} to {d['WL'].max():.2f}; scaled tarsus {d['X'][:, 1].min():.2f} to {d['X'][:, 1].max():.2f};"
f" females {int(d['Z'][:, 0].sum())}; year column {d['G'][:, -1].min():.2f} to {d['G'][:, -1].max():.2f} (decades)")
print(f" G: rows {d['G'].shape[0]}, climate values {d['G'][:, 1:-1].min():.2f} to {d['G'][:, 1:-1].max():.2f}; "
f"birds per species-year {np.bincount(d['spyear'])[1:].min()} to {np.bincount(d['spyear'])[1:].max()}")
Changes made to lm_mod.stan (array syntax of Stan 2.33 and later; nothing else): --- +++ @@ -6,2 +6,2 @@ - int<lower=0> spyear[N]; // Unique species year for each observation - int<lower=0> species_obs[N]; // Unique species for every observation + array[N] int<lower=0> spyear; // Unique species year for each observation + array[N] int<lower=0> species_obs; // Unique species for every observation @@ -11 +11 @@ - int<lower=0> species[J]; // Unique species for each species-year combo + array[J] int<lower=0> species; // Unique species for each species-year combo @@ -15,4 +15,4 @@ - vector[2] B[J]; // A unique set of intercept and TL parameters for every spyear - vector[2] A[K]; // A unique set of sex and TL:sex parameters for every species - vector[E] I[K]; // A unique set of env effects on species-year intercepts for every species - vector[E] S[K]; // A unique set of env effects on species-year slopes for every species + array[J] vector[2] B; // A unique set of intercept and TL parameters for every spyear + array[K] vector[2] A; // A unique set of sex and TL:sex parameters for every species + array[K] vector[E] I; // A unique set of env effects on species-year intercepts for every species + array[K] vector[E] S; // A unique set of env effects on species-year slopes for every species @@ -33 +33 @@ - vector[2] Bhat[J]; + array[J] vector[2] Bhat; HY: N = 11798 birds, J = 177 species-years, K = 5 species, E = 8 columns of G (intercept and PrevTempWinMo, PrevPrecipWinMo, AvgTempPreBreedBase, PrecipPreBreedBase, TempBreedMoBase, PrecipBreedMoBase, y_id) scaled wing -4.49 to 4.56; scaled tarsus -4.83 to 3.16; females 6027; year column -1.70 to 1.80 (decades) G: rows 177, climate values -2.74 to 2.85; birds per species-year 2 to 641 spring AHY: N = 12318 birds, J = 174 species-years, K = 5 species, E = 8 columns of G (intercept and PrevPrecipWinMo, PrecipPreBreedBase, TempBreedMoBase, PrecipBreedMoBase, MinTempFallMig, PrecipWinMo, y_id) scaled wing -4.25 to 3.59; scaled tarsus -4.71 to 4.22; females 5681; year column -1.70 to 1.70 (decades) G: rows 174, climate values -2.95 to 2.65; birds per species-year 3 to 428 autumn AHY: N = 3250 birds, J = 171 species-years, K = 5 species, E = 7 columns of G (intercept and PrevPrevPrecipWinMo, PrevTempBreedBase, PrevPrecipBreedBase, PrevMinTempFallMig, MinTempFallMig, y_id) scaled wing -4.01 to 3.04; scaled tarsus -4.28 to 3.79; females 1568; year column -1.60 to 1.80 (decades) G: rows 171, climate values -2.95 to 2.65; birds per species-year 1 to 114
Timing¶
A run of 200 iterations per chain is not a valid fit; only its time is used.
rows = []
for cls in w.CLASSES:
t0 = time.time()
fit = model.sample(data=data[cls], chains=4, parallel_chains=4, iter_warmup=100, iter_sampling=100,
adapt_delta=0.85, seed=w.STAN_SEEDS[cls], output_dir=str(workdir), show_progress=False)
seconds = time.time() - t0
depth = fit.method_variables()["treedepth__"]
rows.append({"class": cls, "birds": data[cls]["N"], "seconds for 200 iterations": round(seconds),
"mean tree depth (last 100)": round(float(depth.mean()), 2),
"estimate for 4000 iterations, minutes": round(seconds * 20 / 60)})
timing = pd.DataFrame(rows)
print(timing.to_string(index=False))
main_minutes = timing["estimate for 4000 iterations, minutes"]
print(f"\nThe three main models under one reading: about {main_minutes.sum()} minutes; under readings (a), (c) and (d): "
f"about {3 * main_minutes.sum()} minutes.")
print(f"The six sensitivity models (Tables S9 to S11, two per class) under one reading: about "
f"{2 * main_minutes.sum()} minutes; under three readings: about {6 * main_minutes.sum()} minutes.")
print("Notebook limit: 30 minutes. A main model fits in a notebook of its own when its estimate is under about 25 minutes.")
print(f"\nNotebook time: {(time.time() - started) / 60:.1f} minutes")
class birds seconds for 200 iterations mean tree depth (last 100) estimate for 4000 iterations, minutes
HY 11798 266 8.79 89
spring AHY 12318 223 8.31 74
autumn AHY 3250 66 7.89 22
The three main models under one reading: about 185 minutes; under readings (a), (c) and (d): about 555 minutes.
The six sensitivity models (Tables S9 to S11, two per class) under one reading: about 370 minutes; under three readings: about 1110 minutes.
Notebook limit: 30 minutes. A main model fits in a notebook of its own when its estimate is under about 25 minutes.
Notebook time: 11.8 minutes
Result: the Stan fits do not fit the limits. The estimates for 4000 iterations printed above for HY and spring AHY under one filtering reading each exceed the 30-minute notebook limit, and the 18 fits the main and sensitivity claims need for those classes over three readings (the printed totals for three readings) exceed the time of this run. Only the autumn AHY model, whose estimate is under the limit, is fitted, once per reading; the full fits print their own time. The claims that need the other fits are blocked as too large. The estimates are wall-clock times of a short run and change from run to run.
The autumn AHY model is fitted in autumn_ahy_main_a.ipynb, autumn_ahy_main_c.ipynb and autumn_ahy_main_d.ipynb, one notebook per filtering reading, which also record the time of the full fit. The estimate here scales a run of 200 iterations, all of them early warm-up; the autumn AHY full fits print their own time (cell 5 of each), and the attempt record gives their ratio to the autumn AHY estimate above and the HY and spring AHY estimates scaled by it, which stay over the 30-minute notebook limit. The other fits are not run; the claims that need them are blocked as too_large, with the estimates above.