Reproducing Dias et al. (2025) from Dryad dataset bvq83bkk9¶
Paper. Dias AC, Lemoine NP, Yanco SW, Zimova M, Bay RA, Weeks BC (2025). Long-term increases in wing length occur independently of changes in climate and climate-driven shifts in body size. Proc. R. Soc. B 292: 20242556. doi:10.1098/rspb.2024.2556. The claims file was written from the version of record (Europe PMC, PMC11750374) and its electronic supplementary material (Tables S1 to S11, Fig. S1).
Dataset. DRYAD:bvq83bkk9, version id 333377 (Dryad's
version number 7, published 2024-12-09). The paper cites the dataset by its DOI, which always
resolves to the latest version; Dryad's API lists one version, 333377, so it is the one pinned. The
files are Data.zip (SHA-256 2185ac50e9a59c3a1dff59473216de2ab001d0d0044b1fa15442ae18a520d554)
and README.md (6c19de67b1c6c52802be42a6593c175a02bee0dab0ca2f6f958f08a0a3fa0482). Dryad needs a
token to download, so the workflow stages the files before the notebook runs and they are read in
place; the first cell checks their SHA-256.
Outside input. The paper's code, Zenodo
doi:10.5281/zenodo.14226724, Code.zip (SHA-256
cc9f58b9855aec4e3aca09be4138707d07fc34739ec5762e54e1f5bca48ddbc8), read in memory. It gives the
Stan model (lm_mod.stan), the column numbers the scripts select and the rows they drop. No data
come from outside the dataset: the phylogeny the preliminary models use, chitree.rds, is in it.
Notebooks. This notebook holds the data check, the sample sizes and every claim on the
preliminary phylogenetic mixed models (Tables S4 to S8). The hierarchical Stan models of Tables 1
to 3 are slow (see stan_models.ipynb, which times them, and autumn_ahy_main_a.ipynb, _c and _d, which fit
the smallest, once per filtering reading); the claims on them are blocked for time except where the attempt record says otherwise.
| Claim | Tested here | Status |
|---|---|---|
| sample-sizes | yes | tested |
| s4-s6-spring-migration-temperature-none | yes | tested; blocked, methods underspecified (seed readings differ) |
| s5-s6-ahy-prebreeding-wintering-temperature-none | yes | tested |
| s6-autumn-ahy-recent-breeding-temperature-none | yes | tested |
| s4-hy-autumn-migration-temperature-none | yes | tested |
| s6-autumn-ahy-recent-wintering-prebreeding-precipitation-none | yes | tested; blocked, methods underspecified (seed readings differ) |
| s6-autumn-ahy-recent-breeding-precipitation-none | yes | tested |
| s7-low-phylogenetic-signal | yes | tested |
| s4-s6-no-multicollinearity | yes | tested |
| s4-s6-climate-variable-selection | yes | tested |
| s4-s5-climate-size-interactions | yes | tested |
| s8-climate-trends | yes | tested |
| table3-autumn-ahy-prior-wintering-precipitation | preliminary part here, main model in autumn_ahy_main_a.ipynb, _c and _d |
see record |
| claims on the HY and spring AHY main models (Tables 1, 2, S1), climate-trends-cannot-explain-wing-increase, climate-explains-little, s9-s11-expanded-ranges | no | blocked, too_large |
| claims on the autumn AHY main model only (Table 3) | in autumn_ahy_main_a.ipynb, _c and _d |
see record |
In [1]:
import time
import warnings
import matplotlib.pyplot as plt
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)
pd.set_option("display.max_rows", 200)
started = time.time()
print("Dryad dataset", w.DRYAD_ID, "version", w.DRYAD_VERSION)
for name, digest in w.check_staged().items():
print(f"{name}: SHA-256 {digest} (matches Dryad's file list)")
with w.data_zip() as z:
print("\nData.zip holds:")
for info in z.infolist():
if "__MACOSX" not in info.filename and not info.filename.endswith(("/", ".DS_Store")):
print(f" {info.filename} {info.file_size} bytes")
Dryad dataset bvq83bkk9 version 333377 Data.zip: SHA-256 2185ac50e9a59c3a1dff59473216de2ab001d0d0044b1fa15442ae18a520d554 (matches Dryad's file list) README.md: SHA-256 6c19de67b1c6c52802be42a6593c175a02bee0dab0ca2f6f958f08a0a3fa0482 (matches Dryad's file list) Data.zip holds: Data/Datasets/ed_AHYfall.csv 223678 bytes Data/Datasets/ed_AHYspring.csv 221568 bytes Data/Datasets/ed_FHY.csv 223678 bytes Data/Datasets/morphdata.csv 2412347 bytes Data/Files/catnodata.csv 1422546 bytes Data/Files/chitree.rds 1687 bytes
Data check¶
Each check of the description against the values is stated under its output. The dataset has no time series, so there is no clock to align: every table is keyed by species and year.
In [2]:
raw = w.read_csv("morphdata.csv")
print("morphdata.csv:", raw.shape, "columns", list(raw.columns))
print("missing values per column:", raw.isna().sum().to_dict())
print("\nSpecimens by class (Season, NewAge):")
print(raw.groupby(["Season", "NewAge"]).size().to_string())
print("\nSpecimens by species and class:")
print(pd.crosstab(raw.Taxon, [raw.Season, raw.NewAge], margins=True).to_string())
print("\nSex values:", raw.Sex.value_counts().to_dict())
print("Years:", raw.Year.min(), "to", raw.Year.max(), "; species-years with specimens:",
raw.groupby(["Taxon", "Year"]).ngroups)
print("\nTarsus and wing (mm) by species:")
print(raw.groupby("Taxon")[["Tarsus", "Wing"]].agg(["min", "median", "max"]).to_string())
print("\nsize_idx and wing_idx (not used):", raw[["size_idx", "wing_idx"]].describe().loc[["min", "max"]].round(2).to_dict())
morphdata.csv: (27486, 10) columns ['Unnamed: 0', 'Taxon', 'Sex', 'NewAge', 'Year', 'Season', 'Tarsus', 'Wing', 'size_idx', 'wing_idx']
missing values per column: {'Unnamed: 0': 0, 'Taxon': 0, 'Sex': 0, 'NewAge': 0, 'Year': 0, 'Season': 0, 'Tarsus': 0, 'Wing': 0, 'size_idx': 0, 'wing_idx': 0}
Specimens by class (Season, NewAge):
Season NewAge
Fall AHY 3270
HY 11813
Spring AHY 12403
Specimens by species and class:
Season Fall Spring All
NewAge AHY HY AHY
Taxon
Junco_hyemalis 535 3172 1525 5232
Melospiza_georgiana 794 1391 2261 4446
Melospiza_melodia 190 447 4173 4810
Seiurus_aurocapilla 587 1710 1673 3970
Zonotrichia_albicollis 1164 5093 2771 9028
All 3270 11813 12403 27486
Sex values: {'m': 14157, 'f': 13329}
Years: 1981 to 2016 ; species-years with specimens: 180
Tarsus and wing (mm) by species:
Tarsus Wing
min median max min median max
Taxon
Junco_hyemalis 19.2 22.3 25.0 67.5 76.0 90.0
Melospiza_georgiana 19.4 22.8 26.2 52.0 60.0 69.0
Melospiza_melodia 20.0 23.1 25.7 55.5 65.5 76.0
Seiurus_aurocapilla 19.9 23.4 26.3 64.0 74.5 82.0
Zonotrichia_albicollis 21.3 24.9 28.0 61.0 72.0 82.0
size_idx and wing_idx (not used): {'size_idx': {'min': -4.63, 'max': 3.7}, 'wing_idx': {'min': -4.31, 'max': 4.54}}
The columns, species, sexes, ages, seasons and years (1981 to 2016) are those the README lists. No spring specimen is HY, as the Methods say. Tarsus (19.2 to 28.0 mm) and wing (52 to 90 mm) are in plausible ranges for these sparrows, junco and ovenbird. The file holds 27 486 specimens, 120 more than the 27 366 of the paper; the next cells show where they go.
In [3]:
cat = w.read_csv("catnodata.csv")
print("catnodata.csv:", cat.shape, "columns", list(cat.columns))
print("rows without a catalogue number:", int(cat.CatNo.isna().sum()))
known = cat.dropna(subset=["CatNo"])
print("catalogue numbers given more than once:", int(known.CatNo.duplicated(keep=False).sum()),
"rows,", int(known.CatNo.duplicated().sum()), "repeats; institution codes:", cat.InstCode.value_counts(dropna=False).to_dict())
dup = known[known.CatNo.duplicated(keep=False)].sort_values("CatNo")
print("repeated catalogue numbers that differ in species, sex, age, year or season:",
int(dup.groupby("CatNo")[["Taxon", "Sex", "Age", "Year", "Season"]].nunique().gt(1).any(axis=1).sum()))
counts = pd.concat({"morphdata": raw.groupby(["Season", "NewAge"]).size(),
"catnodata": cat.rename(columns={"Age": "NewAge"}).groupby(["Season", "NewAge"]).size()}, axis=1)
print("\nRows by class in the two files:")
print(counts.to_string())
print("\nmorphdata.csv has no catalogue number, so its rows cannot be matched to catnodata.csv one by one.")
same = raw.duplicated(subset=["Taxon", "Sex", "NewAge", "Year", "Season", "Tarsus", "Wing"], keep=False)
print("morphdata rows identical to another row in every recorded field:", int(same.sum()),
f"({same.mean():.1%}); measurements are rounded to 0.1 mm (tarsus) and 0.5 mm (wing), so such ties are expected")
print("wing values that are not a multiple of 0.5 mm:", int((raw.Wing * 2 % 1 != 0).sum()),
"; tarsus values with more than one decimal:", int(((raw.Tarsus * 10).round(6) % 1 != 0).sum()))
catnodata.csv: (26436, 8) columns ['Unnamed: 0', 'Taxon', 'Sex', 'Age', 'Year', 'Season', 'CatNo', 'InstCode']
rows without a catalogue number: 5061
catalogue numbers given more than once: 999 rows, 945 repeats; institution codes: {'FMNH': 26436}
repeated catalogue numbers that differ in species, sex, age, year or season: 44
Rows by class in the two files:
morphdata catnodata
Season NewAge
Fall AHY 3270 3166
HY 11813 11746
Spring AHY 12403 11524
morphdata.csv has no catalogue number, so its rows cannot be matched to catnodata.csv one by one.
morphdata rows identical to another row in every recorded field: 5514 (20.1%); measurements are rounded to 0.1 mm (tarsus) and 0.5 mm (wing), so such ties are expected
wing values that are not a multiple of 0.5 mm: 17 ; tarsus values with more than one decimal: 8
catnodata.csv, which the paper cites for catalogue numbers, has 26 436 rows, 1050 fewer than
morphdata.csv, and 5061 of them have no catalogue number. Catalogue numbers repeat, so the file
cannot show that each specimen was deposited once, and morphdata.csv has no identifier to link a
measurement to a specimen. Identical rows in morphdata.csv are expected from the rounding of the
measurements and are not evidence of duplicates.
In [4]:
ed = {cls: w.read_csv(w.ED_FILE[cls]) for cls in w.CLASSES}
for cls, df in ed.items():
print(f"{w.ED_FILE[cls]}: {df.shape}, y_id {df.y_id.min()} to {df.y_id.max()}, species {df.spec_id.nunique()}")
print("columns of ed_FHY.csv:", list(ed["HY"].columns))
print("columns of ed_AHYspring.csv not in ed_FHY.csv:", sorted(set(ed["spring AHY"].columns) - set(ed["HY"].columns)))
print("columns of ed_FHY.csv not in ed_AHYspring.csv:", sorted(set(ed["HY"].columns) - set(ed["spring AHY"].columns)))
print("\ned_FHY.csv and ed_AHYfall.csv are identical:", ed["HY"].equals(ed["autumn AHY"]))
print("ed_FHY.csv Year2 - y_id:", (ed["HY"].Year2 - ed["HY"].y_id).unique().tolist(),
"; ed_AHYspring.csv y_id - Year:", (ed["spring AHY"].y_id - ed["spring AHY"].Year).unique().tolist())
ed_hy, ed_sp = ed["HY"], ed["spring AHY"]
both = ed_hy.merge(ed_sp, left_on=["spec_id", "y_id"], right_on=["spec_id", "Year"], suffixes=("_f", "_s"))
climate_cols = [c for c in ed_sp.columns if c not in ("Unnamed: 0", "y_id", "spec_id", "Species2", "Year")]
diff = max(np.nanmax(np.abs(both[c + "_f"] - both[c + "_s"])) for c in climate_cols if c + "_f" in both)
print(f"ed_AHYspring.csv row for Year Y against ed_FHY.csv row for y_id Y, {len(both)} species-years: "
f"largest difference over the shared climate columns {diff}")
lag = ed_hy.merge(ed_hy, left_on=["spec_id", "y_id"], right_on=["spec_id", "Year2"], suffixes=("", "_prev"))
pairs = [("PrevTempWinMo", "TempWinMo"), ("PrevPrecipWinMo", "PrecipWinMo"), ("PrevTempBreedBase", "TempBreedMoBase"),
("PrevPrecipBreedBase", "PrecipBreedMoBase"), ("PrevMinTempFallMig", "MinTempFallMig"),
("PrevAvgTempPreBreedBase", "AvgTempPreBreedBase"), ("PrevMinTempSpringMig", "MinTempSpringMig")]
print("'Prev' columns against the same variable one year earlier (largest difference):",
{p: float(np.nanmax(np.abs(lag[p] - lag[c + "_prev"]))) for p, c in pairs})
values = ed_hy[[c for c in climate_cols if c in ed_hy.columns]]
per_species = ed_hy.groupby("spec_id")[["TempBreedMoBase", "PrecipBreedMoBase", "MinTempFallMig", "PrecipWinMo"]].agg(["mean", "std"])
print("\nClimate values (ed_FHY.csv), all columns: min", round(float(np.nanmin(values)), 2), "max", round(float(np.nanmax(values)), 2))
print("Mean and s.d. within species of four columns (README: scaled to mean 0 and s.d. 1 within species):")
print(per_species.round(3).to_string())
print("\nMissing values per climate column used by the models (ed_FHY.csv):")
print(ed_hy[["PrevTempWinMo", "PrevPrecipWinMo", "AvgTempPreBreedBase", "PrecipPreBreedBase", "TempBreedMoBase",
"PrecipBreedMoBase", "MinTempFallMig", "MinTempSpringMig", "PrevPrevPrecipWinMo", "PrevTempBreedBase"]].isna().sum().to_dict())
ed_FHY.csv: (180, 70), y_id 1981 to 2016, species 5
ed_AHYspring.csv: (180, 70), y_id 1981 to 2016, species 5
ed_AHYfall.csv: (180, 70), y_id 1981 to 2016, species 5
columns of ed_FHY.csv: ['Unnamed: 0', 'y_id', 'spec_id', 'Species2', 'MinTempSpringMig', 'MinTempFallMig', 'MinTempNextSpringMig', 'PrecipSpringMig', 'PrecipPreBreed', 'PrecipFallMig', 'PrecipNextSpringMig', 'TempBreedMo', 'TempWinMo', 'PrecipBreedMo', 'PrecipWinMo', 'TempBreedMoBase', 'TempBreedMoB1', 'TempBreedMoB2', 'TempBreedMoS1', 'TempBreedMoS2', 'TempBreedMoS3', 'PrecipBreedMoBase', 'PrecipBreedMoB1', 'PrecipBreedMoB2', 'PrecipBreedMoS1', 'PrecipBreedMoS2', 'PrecipBreedMoS3', 'PrecipPreBreedBase', 'PrecipPreBreedB1', 'PrecipPreBreedB2', 'PrecipPreBreedS1', 'PrecipPreBreedS2', 'PrecipPreBreedS3', 'AvgTempPreBreedBase', 'AvgTempPreBreedB1', 'AvgTempPreBreedB2', 'AvgTempPreBreedS1', 'AvgTempPreBreedS2', 'AvgTempPreBreedS3', 'Year2', 'PrevTempWinMo', 'PrevPrecipWinMo', 'PrevAvgTempPreBreedBase', 'PrevAvgTempPreBreedB1', 'PrevAvgTempPreBreedB2', 'PrevAvgTempPreBreedS1', 'PrevAvgTempPreBreedS2', 'PrevAvgTempPreBreedS3', 'PrevPrecipPreBreedBase', 'PrevPrecipPreBreedB1', 'PrevPrecipPreBreedB2', 'PrevPrecipPreBreedS1', 'PrevPrecipPreBreedS2', 'PrevPrecipPreBreedS3', 'PrevTempBreedBase', 'PrevTempBreedB1', 'PrevTempBreedB2', 'PrevTempBreedS1', 'PrevTempBreedS2', 'PrevTempBreedS3', 'PrevPrecipBreedBase', 'PrevPrecipBreedB1', 'PrevPrecipBreedB2', 'PrevPrecipBreedS1', 'PrevPrecipBreedS2', 'PrevPrecipBreedS3', 'PrevMinTempFallMig', 'PrevPrevTempWinMo', 'PrevPrevPrecipWinMo', 'PrevMinTempSpringMig']
columns of ed_AHYspring.csv not in ed_FHY.csv: ['Year']
columns of ed_FHY.csv not in ed_AHYspring.csv: ['Year2']
ed_FHY.csv and ed_AHYfall.csv are identical: True
ed_FHY.csv Year2 - y_id: [1] ; ed_AHYspring.csv y_id - Year: [1]
ed_AHYspring.csv row for Year Y against ed_FHY.csv row for y_id Y, 175 species-years: largest difference over the shared climate columns 0.0
'Prev' columns against the same variable one year earlier (largest difference): {'PrevTempWinMo': 0.0, 'PrevPrecipWinMo': 0.0, 'PrevTempBreedBase': 0.0, 'PrevPrecipBreedBase': 0.0, 'PrevMinTempFallMig': 0.0, 'PrevAvgTempPreBreedBase': 0.0, 'PrevMinTempSpringMig': 0.0}
Climate values (ed_FHY.csv), all columns: min -2.95 max 3.57
Mean and s.d. within species of four columns (README: scaled to mean 0 and s.d. 1 within species):
TempBreedMoBase PrecipBreedMoBase MinTempFallMig PrecipWinMo
mean std mean std mean std mean std
spec_id
Junco_hyemalis 0.013 1.023 0.040 1.012 -0.000 1.028 0.004 1.026
Melospiza_georgiana 0.012 1.023 0.038 1.013 -0.060 0.992 0.020 1.020
Melospiza_melodia 0.011 1.024 0.033 1.015 0.028 1.020 0.002 1.027
Seiurus_aurocapilla 0.004 1.026 0.029 1.020 -0.091 0.949 0.004 0.977
Zonotrichia_albicollis 0.014 1.023 0.038 1.012 -0.043 0.995 0.013 1.025
Missing values per climate column used by the models (ed_FHY.csv):
{'PrevTempWinMo': 0, 'PrevPrecipWinMo': 0, 'AvgTempPreBreedBase': 2, 'PrecipPreBreedBase': 0, 'TempBreedMoBase': 0, 'PrecipBreedMoBase': 0, 'MinTempFallMig': 1, 'MinTempSpringMig': 4, 'PrevPrevPrecipWinMo': 5, 'PrevTempBreedBase': 0}
The README says the climate is scaled to mean 0 and s.d. 1 within species. Over the years in the
files the means are near 0 and the s.d. near 1, but not exactly; this is consistent with scaling
over a longer span (the Methods give 1980 to 2017), not shown. The README says ed_FHY.csv and
ed_AHYfall.csv hold the climate of HY and autumn AHY birds; the two files are identical, which
fits, since both classes are collected in the autumn of y_id and the autumn AHY model uses the
"Prev" columns. The spring file's Year is y_id - 1, as the README says, and its rows equal the
ed_FHY.csv rows of the same year, so the three files are one climate table. The "Prev" columns
are the variable of the year before.
In [5]:
tree = w.read_tree()
print("chitree.rds: a phylo object with", len(tree["tip_label"]), "tips and", tree["Nnode"], "internal nodes")
print("the five species are tips:", all(s in tree["tip_label"] for s in w.SPECIES))
cov = w.tip_covariance(tree, w.SPECIES)
print("\nShared branch length from the root (vcv.phylo of the five tips, millions of years as in BirdTree):")
print(cov.round(2).to_string())
print("ultrametric over the five tips (equal root-to-tip lengths):", bool(np.allclose(np.diag(cov), cov.iloc[0, 0])))
chitree.rds: a phylo object with 46 tips and 45 internal nodes
the five species are tips: True
Shared branch length from the root (vcv.phylo of the five tips, millions of years as in BirdTree):
Junco_hyemalis Melospiza_georgiana Melospiza_melodia Seiurus_aurocapilla Zonotrichia_albicollis
Junco_hyemalis 78.76 69.62 69.62 59.90 73.70
Melospiza_georgiana 69.62 78.76 77.25 59.90 69.62
Melospiza_melodia 69.62 77.25 78.76 59.90 69.62
Seiurus_aurocapilla 59.90 59.90 59.90 78.76 59.90
Zonotrichia_albicollis 73.70 69.62 69.62 59.90 78.76
ultrametric over the five tips (equal root-to-tip lengths): True
Data preparation and sample sizes¶
The scripts select climate columns by their number in R's merged table and drop rows by
hard-coded number. wingallom.prepare rebuilds the merged table, selects by the same numbers,
finds the rows with a missing value and drops those; the cell below shows that they are the rows
the code drops, and which columns the numbers select.
In [6]:
READINGS = {"a": "as deposited, no filter", "b": "5 MAD per species, R mad()", "c": "5 MAD per species and sex",
"d": "5 MAD per species, unscaled MAD"}
removed = {r: w.outlier_mask(raw, r) for r in READINGS}
print("Specimens the outlier filter removes, by reading:")
for r, mask in removed.items():
by = raw[mask].groupby(["Season", "NewAge"]).size().to_dict()
print(f" ({r}) {READINGS[r]}: {int(mask.sum())} by class {by}")
print("readings (a) and (b) coincide:", bool((removed['a'] == removed['b']).all()))
RUN = ["a", "c", "d"] # (b) removes nothing, so it is (a)
morph = {r: w.morph_prepared(r) for r in RUN}
prep = {(r, which, cls): w.prepare(morph[r], cls, which) for r in RUN for which in ("main", "prelim") for cls in w.CLASSES}
for which in ("main", "prelim"):
for cls in w.CLASSES:
p = prep[("a", which, cls)]
print(f"\n{which} {cls}: {p['n_merged']} species-years with specimens; rows with a missing value {p['missing_rows']};"
f" rows the code drops {p['dropped_rows']}; same: {p['missing_rows'] == p['dropped_rows']}")
print(" selected columns:", p["selected"])
Specimens the outlier filter removes, by reading:
(a) as deposited, no filter: 0 by class {}
(b) 5 MAD per species, R mad(): 0 by class {}
(c) 5 MAD per species and sex: 17 by class {('Fall', 'AHY'): 1, ('Fall', 'HY'): 11, ('Spring', 'AHY'): 5}
(d) 5 MAD per species, unscaled MAD: 70 by class {('Fall', 'AHY'): 10, ('Fall', 'HY'): 35, ('Spring', 'AHY'): 25}
readings (a) and (b) coincide: True
main HY: 179 species-years with specimens; rows with a missing value [np.int64(3), np.int64(75)]; rows the code drops [3, 75]; same: True selected columns: ['PrevTempWinMo', 'PrevPrecipWinMo', 'AvgTempPreBreedBase', 'AvgTempPreBreedB1', 'AvgTempPreBreedB2', 'AvgTempPreBreedS1', 'AvgTempPreBreedS2', 'AvgTempPreBreedS3', 'PrecipPreBreedBase', 'PrecipPreBreedB1', 'PrecipPreBreedB2', 'PrecipPreBreedS1', 'PrecipPreBreedS2', 'PrecipPreBreedS3', 'TempBreedMoBase', 'TempBreedMoB1', 'TempBreedMoB2', 'TempBreedMoS1', 'TempBreedMoS2', 'TempBreedMoS3', 'PrecipBreedMoBase', 'PrecipBreedMoB1', 'PrecipBreedMoB2', 'PrecipBreedMoS1', 'PrecipBreedMoS2', 'PrecipBreedMoS3', 'spec_id', 'y_id', 'spyear'] main spring AHY: 180 species-years with specimens; rows with a missing value [np.int64(1), np.int64(37), np.int64(73), np.int64(109), np.int64(115), np.int64(145)]; rows the code drops [1, 37, 73, 109, 115, 145]; same: True selected columns: ['PrevPrecipWinMo', 'PrecipPreBreedBase', 'PrecipPreBreedB1', 'PrecipPreBreedB2', 'PrecipPreBreedS1', 'PrecipPreBreedS2', 'PrecipPreBreedS3', 'TempBreedMoBase', 'TempBreedMoB1', 'TempBreedMoB2', 'TempBreedMoS1', 'TempBreedMoS2', 'TempBreedMoS3', 'PrecipBreedMoBase', 'PrecipBreedMoB1', 'PrecipBreedMoB2', 'PrecipBreedMoS1', 'PrecipBreedMoS2', 'PrecipBreedMoS3', 'MinTempFallMig', 'PrecipWinMo', 'spec_id', 'y_id', 'spyear'] main autumn AHY: 177 species-years with specimens; rows with a missing value [np.int64(1), np.int64(36), np.int64(106), np.int64(111), np.int64(112), np.int64(142)]; rows the code drops [1, 36, 106, 111, 112, 142]; same: True selected columns: ['PrevPrevPrecipWinMo', 'PrevTempBreedBase', 'PrevTempBreedB1', 'PrevTempBreedB2', 'PrevTempBreedS1', 'PrevTempBreedS2', 'PrevTempBreedS3', 'PrevPrecipBreedBase', 'PrevPrecipBreedB1', 'PrevPrecipBreedB2', 'PrevPrecipBreedS1', 'PrevPrecipBreedS2', 'PrevPrecipBreedS3', 'PrevMinTempFallMig', 'MinTempFallMig', 'spec_id', 'y_id', 'spyear'] prelim HY: 179 species-years with specimens; rows with a missing value [np.int64(3), np.int64(6), np.int64(39), np.int64(75), np.int64(78), np.int64(113), np.int64(146)]; rows the code drops [3, 6, 39, 75, 78, 113, 146]; same: True selected columns: ['PrevTempWinMo', 'PrevPrecipWinMo', 'MinTempSpringMig', 'AvgTempPreBreedBase', 'AvgTempPreBreedB1', 'AvgTempPreBreedB2', 'AvgTempPreBreedS1', 'AvgTempPreBreedS2', 'AvgTempPreBreedS3', 'PrecipPreBreedBase', 'PrecipPreBreedB1', 'PrecipPreBreedB2', 'PrecipPreBreedS1', 'PrecipPreBreedS2', 'PrecipPreBreedS3', 'TempBreedMoBase', 'TempBreedMoB1', 'TempBreedMoB2', 'TempBreedMoS1', 'TempBreedMoS2', 'TempBreedMoS3', 'PrecipBreedMoBase', 'PrecipBreedMoB1', 'PrecipBreedMoB2', 'PrecipBreedMoS1', 'PrecipBreedMoS2', 'PrecipBreedMoS3', 'MinTempFallMig', 'spec_id', 'y_id', 'spyear'] prelim spring AHY: 180 species-years with specimens; rows with a missing value [np.int64(1), np.int64(4), np.int64(6), np.int64(7), np.int64(37), np.int64(39), np.int64(40), np.int64(73), np.int64(76), np.int64(78), np.int64(79), np.int64(109), np.int64(115), np.int64(145), np.int64(147), np.int64(148)]; rows the code drops [1, 4, 6, 7, 37, 39, 40, 73, 76, 78, 79, 109, 115, 145, 147, 148]; same: True selected columns: ['PrevTempWinMo', 'PrevPrecipWinMo', 'MinTempSpringMig', 'AvgTempPreBreedBase', 'AvgTempPreBreedB1', 'AvgTempPreBreedB2', 'AvgTempPreBreedS1', 'AvgTempPreBreedS2', 'AvgTempPreBreedS3', 'PrecipPreBreedBase', 'PrecipPreBreedB1', 'PrecipPreBreedB2', 'PrecipPreBreedS1', 'PrecipPreBreedS2', 'PrecipPreBreedS3', 'TempBreedMoBase', 'TempBreedMoB1', 'TempBreedMoB2', 'TempBreedMoS1', 'TempBreedMoS2', 'TempBreedMoS3', 'PrecipBreedMoBase', 'PrecipBreedMoB1', 'PrecipBreedMoB2', 'PrecipBreedMoS1', 'PrecipBreedMoS2', 'PrecipBreedMoS3', 'MinTempFallMig', 'TempWinMo', 'PrecipWinMo', 'MinTempNextSpringMig', 'spec_id', 'y_id', 'spyear'] prelim autumn AHY: 177 species-years with specimens; rows with a missing value [np.int64(1), np.int64(3), np.int64(4), np.int64(6), np.int64(7), np.int64(36), np.int64(38), np.int64(39), np.int64(73), np.int64(74), np.int64(76), np.int64(77), np.int64(106), np.int64(111), np.int64(112), np.int64(142), np.int64(144), np.int64(145)]; rows the code drops [1, 3, 4, 6, 7, 36, 38, 39, 73, 74, 76, 77, 106, 111, 112, 142, 144, 145]; same: True selected columns: ['PrevPrevTempWinMo', 'PrevPrevPrecipWinMo', 'PrevMinTempSpringMig', 'PrevAvgTempPreBreedBase', 'PrevAvgTempPreBreedB1', 'PrevAvgTempPreBreedB2', 'PrevAvgTempPreBreedS1', 'PrevAvgTempPreBreedS2', 'PrevAvgTempPreBreedS3', 'PrevPrecipPreBreedBase', 'PrevPrecipPreBreedB1', 'PrevPrecipPreBreedB2', 'PrevPrecipPreBreedS1', 'PrevPrecipPreBreedS2', 'PrevPrecipPreBreedS3', 'PrevTempBreedBase', 'PrevTempBreedB1', 'PrevTempBreedB2', 'PrevTempBreedS1', 'PrevTempBreedS2', 'PrevTempBreedS3', 'PrevPrecipBreedBase', 'PrevPrecipBreedB1', 'PrevPrecipBreedB2', 'PrevPrecipBreedS1', 'PrevPrecipBreedS2', 'PrevPrecipBreedS3', 'PrevMinTempFallMig', 'PrevTempWinMo', 'PrevPrecipWinMo', 'MinTempSpringMig', 'AvgTempPreBreedBase', 'AvgTempPreBreedB1', 'AvgTempPreBreedB2', 'AvgTempPreBreedS1', 'AvgTempPreBreedS2', 'AvgTempPreBreedS3', 'PrecipPreBreedBase', 'PrecipPreBreedB1', 'PrecipPreBreedB2', 'PrecipPreBreedS1', 'PrecipPreBreedS2', 'PrecipPreBreedS3', 'TempBreedMoBase', 'TempBreedMoB1', 'TempBreedMoB2', 'TempBreedMoS1', 'TempBreedMoS2', 'TempBreedMoS3', 'PrecipBreedMoBase', 'PrecipBreedMoB1', 'PrecipBreedMoB2', 'PrecipBreedMoS1', 'PrecipBreedMoS2', 'PrecipBreedMoS3', 'MinTempFallMig', 'spec_id', 'y_id', 'spyear']
In [7]:
expected = {("all", "specimens"): 27366,
("main", "HY", "species-years"): 177, ("main", "HY", "birds"): 11798,
("main", "spring AHY", "species-years"): 174, ("main", "spring AHY", "birds"): 12318,
("main", "autumn AHY", "species-years"): 171, ("main", "autumn AHY", "birds"): 3250,
("prelim", "HY", "birds"): 11731, ("prelim", "spring AHY", "birds"): 11481,
("prelim", "autumn AHY", "birds"): 3122}
rows = []
for key, value in expected.items():
row = {"quantity": " ".join(key), "paper": value, "tolerance": round(value * 0.1)}
for r in RUN:
if key[0] == "all":
found = len(morph[r])
else:
p = prep[(r, key[0], key[1])]
found = len(p["G"]) if key[2] == "species-years" else len(p["md"])
row[f"({r})"] = found
rows.append(row)
sizes = pd.DataFrame(rows)
for r in RUN:
sizes[f"within ({r})"] = (sizes[f"({r})"] - sizes.paper).abs() <= sizes.tolerance
print(sizes.to_string(index=False))
modelled = sum(len(prep[("a", "main", c)]["md"]) for c in w.CLASSES)
print("\nSpecimens in the three main models together, reading (a):", modelled,
"; specimens in species-years dropped for missing climate:", len(morph["a"]) - modelled)
print("Main-model specimens by species, reading (a):")
print(pd.DataFrame({c: prep[("a", "main", c)]["md"].Taxon.value_counts() for c in w.CLASSES}).to_string())
quantity paper tolerance (a) (c) (d) within (a) within (c) within (d)
all specimens 27366 2737 27486 27469 27416 True True True
main HY species-years 177 18 177 177 177 True True True
main HY birds 11798 1180 11798 11787 11763 True True True
main spring AHY species-years 174 17 174 174 174 True True True
main spring AHY birds 12318 1232 12318 12314 12293 True True True
main autumn AHY species-years 171 17 171 171 171 True True True
main autumn AHY birds 3250 325 3250 3249 3241 True True True
prelim HY birds 11731 1173 11731 11720 11696 True True True
prelim spring AHY birds 11481 1148 11481 11479 11457 True True True
prelim autumn AHY birds 3122 312 3122 3121 3114 True True True
Specimens in the three main models together, reading (a): 27366 ; specimens in species-years dropped for missing climate: 120
Main-model specimens by species, reading (a):
HY spring AHY autumn AHY
Taxon
Junco_hyemalis 3159 1505 529
Melospiza_georgiana 1391 2252 793
Melospiza_melodia 445 4151 190
Seiurus_aurocapilla 1710 1647 575
Zonotrichia_albicollis 5093 2763 1163
Result: reproduced. As deposited (reading a), every count equals the paper's: 177, 174 and 171 species-years and 11 798, 12 318 and 3250 birds in the main models, and 11 731, 11 481 and 3122 birds in the preliminary models. The three main models together hold 27 366 birds, the paper's total; the file's 27 486 include 120 birds in species-years dropped for a missing climate value. Readings (c) and (d) remove 17 and 70 specimens and every count stays within tolerance, but only reading (a) gives the paper's counts exactly. The rows the code drops by number are the rows with a missing value in the columns it selects, for all six tables.
Preliminary phylogenetic mixed models (Tables S4 to S7)¶
Script02 fits, per class, scaled wing on scaled tarsus, sex, every climate variable and year, each
with its tarsus interaction, with a random effect of species whose covariance is the phylogeny;
inverse-Wishart priors V = 1, nu = 0.02; 4 chains of 15 000 iterations with 5000 burn-in (thinning
10, MCMCglmm's default). MCMCglmm is an R package and cannot be installed in this environment;
wingallom.gibbs_glmm samples the same model with the same Gibbs scheme. The code sets no seed;
the seeds here are a standard choice, one generator per class, reading and chain. The models are
fitted under filtering readings (a), (c) and (d), as the claims require.
Interval readings, as the claims file names them: (i) the 2.5 and 97.5% quantiles of the four pooled chains; (ii) the 95% HPD interval of chain 4. pMCMC is from chain 4.
In [8]:
Ainv = w.scaled_tip_inverse(tree, w.SPECIES)
print("Inverse of the scaled tip covariance (inverseA, scale = TRUE), species in order", w.SPECIES)
print(np.round(Ainv, 3))
prelim = {}
t0 = time.time()
for ri, r in enumerate(RUN):
for ci, cls in enumerate(w.CLASSES):
X, names, y, groups = w.prelim_design(prep[(r, "prelim", cls)], w.PRELIM_TERMS[cls])
chains = [w.gibbs_glmm(X, y, groups, w.SPECIES, Ainv, np.random.default_rng([2025, ri, ci, k])) for k in range(4)]
lam = chains[3]["VCV"][:, 0] / chains[3]["VCV"].sum(axis=1)
prelim[(r, cls)] = {"summary": w.glmm_summary(chains, names), "vif": w.vif(X, names), "lambda": lam,
"n": len(y), "p": X.shape[1], "chains": chains, "names": names}
print(f"\n{len(prelim)} models of 4 chains in {time.time() - t0:.0f} s")
print(pd.DataFrame({k: {"birds": v["n"], "fixed effects": v["p"], "kept draws per chain": len(v["chains"][0]["Sol"])}
for k, v in prelim.items()}).T.to_string())
Inverse of the scaled tip covariance (inverseA, scale = TRUE), species in order ['Junco_hyemalis', 'Melospiza_georgiana', 'Melospiza_melodia', 'Seiurus_aurocapilla', 'Zonotrichia_albicollis'] [[ 9.319 -1.138 -1.138 -0.603 -6.251] [ -1.138 27.482 -24.581 -0.476 -1.138] [ -1.138 -24.581 27.482 -0.476 -1.138] [ -0.603 -0.476 -0.476 2.641 -0.603] [ -6.251 -1.138 -1.138 -0.603 9.319]]
9 models of 4 chains in 50 s
birds fixed effects kept draws per chain
a HY 11731 22 1000
spring AHY 11481 28 1000
autumn AHY 3122 38 1000
c HY 11720 22 1000
spring AHY 11479 28 1000
autumn AHY 3121 38 1000
d HY 11696 22 1000
spring AHY 11457 28 1000
autumn AHY 3114 38 1000
In [9]:
def gelman_rubin(chains, k):
x = np.array([c["Sol"][:, k] for c in chains])
n = x.shape[1]
wv = x.var(axis=1, ddof=1).mean()
b = n * x.mean(axis=1).var(ddof=1)
return np.sqrt(((n - 1) / n * wv + b / n) / wv)
for (r, cls), v in prelim.items():
psrf = [gelman_rubin(v["chains"], k) for k in range(v["p"])]
print(f"({r}) {cls}: largest Gelman-Rubin factor over the fixed effects {max(psrf):.3f}")
for cls in w.CLASSES:
print(f"\nTable S{4 + w.CLASSES.index(cls)}, {cls}, reading (a):")
print(prelim[("a", cls)]["summary"].round(3).to_string())
(a) HY: largest Gelman-Rubin factor over the fixed effects 1.001
(a) spring AHY: largest Gelman-Rubin factor over the fixed effects 1.001
(a) autumn AHY: largest Gelman-Rubin factor over the fixed effects 1.001
(c) HY: largest Gelman-Rubin factor over the fixed effects 1.001
(c) spring AHY: largest Gelman-Rubin factor over the fixed effects 1.001
(c) autumn AHY: largest Gelman-Rubin factor over the fixed effects 1.001
(d) HY: largest Gelman-Rubin factor over the fixed effects 1.002
(d) spring AHY: largest Gelman-Rubin factor over the fixed effects 1.001
(d) autumn AHY: largest Gelman-Rubin factor over the fixed effects 1.002
Table S4, HY, reading (a):
mean_pooled q2.5_pooled q97.5_pooled mean_c4 hpd_lo_c4 hpd_hi_c4 pMCMC_c4 excl0_i excl0_ii
term
(Intercept) 0.457 -0.052 0.964 0.464 -0.052 0.935 0.066 False False
scaleTarsus 0.113 0.095 0.131 0.112 0.094 0.131 0.001 True True
Sex21 -1.320 -1.345 -1.296 -1.320 -1.343 -1.294 0.001 True True
PrevTempWinMo 0.018 0.005 0.030 0.018 0.006 0.030 0.004 True True
PrevPrecipWinMo -0.019 -0.031 -0.007 -0.019 -0.032 -0.008 0.004 True True
MinTempSpringMig -0.011 -0.024 0.001 -0.011 -0.024 0.001 0.068 False False
AvgTempPreBreedBase 0.046 0.029 0.063 0.046 0.028 0.061 0.001 True True
PrecipPreBreedBase -0.027 -0.042 -0.012 -0.027 -0.042 -0.012 0.001 True True
TempBreedMoBase -0.065 -0.083 -0.048 -0.065 -0.083 -0.048 0.001 True True
PrecipBreedMoBase -0.046 -0.064 -0.029 -0.046 -0.062 -0.029 0.001 True True
MinTempFallMig -0.007 -0.020 0.006 -0.007 -0.020 0.005 0.294 False False
y_id 0.177 0.158 0.196 0.177 0.158 0.196 0.001 True True
scaleTarsus:Sex21 0.025 0.002 0.049 0.025 0.001 0.047 0.030 True True
scaleTarsus:PrevTempWinMo 0.016 0.003 0.028 0.016 0.005 0.028 0.010 True True
scaleTarsus:PrevPrecipWinMo 0.003 -0.010 0.015 0.003 -0.010 0.015 0.670 False False
scaleTarsus:MinTempSpringMig 0.005 -0.008 0.017 0.005 -0.008 0.016 0.406 False False
scaleTarsus:AvgTempPreBreedBase -0.016 -0.032 0.001 -0.016 -0.032 0.002 0.056 False False
scaleTarsus:PrecipPreBreedBase 0.005 -0.009 0.020 0.006 -0.009 0.020 0.426 False False
scaleTarsus:TempBreedMoBase 0.013 -0.004 0.031 0.013 -0.004 0.031 0.138 False False
scaleTarsus:PrecipBreedMoBase 0.008 -0.008 0.025 0.009 -0.009 0.023 0.290 False False
scaleTarsus:MinTempFallMig -0.000 -0.013 0.013 0.000 -0.013 0.013 0.998 False False
scaleTarsus:y_id -0.008 -0.026 0.010 -0.009 -0.027 0.009 0.350 False False
Table S5, spring AHY, reading (a):
mean_pooled q2.5_pooled q97.5_pooled mean_c4 hpd_lo_c4 hpd_hi_c4 pMCMC_c4 excl0_i excl0_ii
term
(Intercept) 0.650 0.111 1.176 0.650 0.098 1.128 0.022 True True
scaleTarsus 0.133 0.113 0.153 0.133 0.113 0.153 0.001 True True
Sex21 -1.385 -1.412 -1.358 -1.385 -1.412 -1.358 0.001 True True
PrevTempWinMo -0.015 -0.031 0.002 -0.015 -0.030 0.003 0.080 False False
PrevPrecipWinMo -0.023 -0.038 -0.008 -0.023 -0.037 -0.006 0.004 True True
MinTempSpringMig -0.006 -0.021 0.009 -0.006 -0.021 0.009 0.462 False False
AvgTempPreBreedBase 0.013 -0.002 0.029 0.013 -0.003 0.027 0.100 False False
PrecipPreBreedBase -0.023 -0.038 -0.008 -0.023 -0.038 -0.008 0.001 True True
TempBreedMoBase -0.029 -0.047 -0.012 -0.029 -0.046 -0.012 0.002 True True
PrecipBreedMoBase -0.029 -0.046 -0.012 -0.030 -0.046 -0.011 0.002 True True
MinTempFallMig -0.035 -0.049 -0.020 -0.035 -0.050 -0.021 0.001 True True
TempWinMo 0.009 -0.007 0.024 0.009 -0.006 0.025 0.274 False False
PrecipWinMo -0.021 -0.037 -0.005 -0.021 -0.037 -0.005 0.018 True True
MinTempNextSpringMig -0.013 -0.028 0.001 -0.014 -0.028 -0.000 0.056 False True
y_id 0.174 0.153 0.195 0.174 0.152 0.193 0.001 True True
scaleTarsus:Sex21 -0.010 -0.038 0.018 -0.010 -0.041 0.014 0.496 False False
scaleTarsus:PrevTempWinMo -0.001 -0.017 0.014 -0.001 -0.017 0.014 0.884 False False
scaleTarsus:PrevPrecipWinMo 0.015 0.000 0.030 0.015 0.001 0.030 0.038 True True
scaleTarsus:MinTempSpringMig 0.010 -0.005 0.025 0.010 -0.005 0.025 0.214 False False
scaleTarsus:AvgTempPreBreedBase -0.008 -0.023 0.007 -0.008 -0.024 0.007 0.290 False False
scaleTarsus:PrecipPreBreedBase -0.009 -0.024 0.005 -0.009 -0.024 0.005 0.210 False False
scaleTarsus:TempBreedMoBase -0.003 -0.021 0.015 -0.003 -0.020 0.015 0.690 False False
scaleTarsus:PrecipBreedMoBase -0.003 -0.020 0.014 -0.003 -0.019 0.013 0.740 False False
scaleTarsus:MinTempFallMig 0.011 -0.004 0.026 0.011 -0.004 0.025 0.166 False False
scaleTarsus:TempWinMo -0.003 -0.019 0.013 -0.003 -0.020 0.011 0.716 False False
scaleTarsus:PrecipWinMo -0.002 -0.018 0.014 -0.002 -0.017 0.014 0.794 False False
scaleTarsus:MinTempNextSpringMig 0.012 -0.003 0.026 0.012 -0.004 0.026 0.118 False False
scaleTarsus:y_id 0.006 -0.013 0.025 0.006 -0.014 0.025 0.562 False False
Table S6, autumn AHY, reading (a):
mean_pooled q2.5_pooled q97.5_pooled mean_c4 hpd_lo_c4 hpd_hi_c4 pMCMC_c4 excl0_i excl0_ii
term
(Intercept) 0.980 0.391 1.622 0.977 0.388 1.597 0.016 True True
scaleTarsus 0.112 0.073 0.149 0.112 0.077 0.149 0.001 True True
Sex21 -1.453 -1.501 -1.408 -1.453 -1.496 -1.405 0.001 True True
PrevPrevTempWinMo 0.028 -0.005 0.060 0.029 -0.004 0.059 0.086 False False
PrevPrevPrecipWinMo 0.030 0.003 0.058 0.030 0.002 0.058 0.036 True True
PrevMinTempSpringMig -0.010 -0.041 0.019 -0.010 -0.040 0.020 0.496 False False
PrevAvgTempPreBreedBase 0.020 -0.015 0.052 0.019 -0.019 0.047 0.262 False False
PrevPrecipPreBreedBase 0.027 -0.002 0.055 0.027 -0.003 0.052 0.068 False False
PrevTempBreedBase -0.118 -0.159 -0.076 -0.117 -0.156 -0.074 0.001 True True
PrevPrecipBreedBase -0.057 -0.094 -0.019 -0.057 -0.096 -0.022 0.006 True True
PrevMinTempFallMig -0.048 -0.074 -0.022 -0.048 -0.072 -0.020 0.001 True True
PrevTempWinMo 0.016 -0.016 0.048 0.016 -0.015 0.048 0.336 False False
PrevPrecipWinMo -0.020 -0.047 0.007 -0.020 -0.051 0.004 0.144 False False
MinTempSpringMig -0.002 -0.030 0.026 -0.002 -0.031 0.025 0.896 False False
AvgTempPreBreedBase 0.001 -0.036 0.037 0.000 -0.037 0.034 0.996 False False
PrecipPreBreedBase 0.003 -0.029 0.034 0.004 -0.028 0.035 0.804 False False
TempBreedMoBase 0.022 -0.017 0.061 0.023 -0.015 0.062 0.266 False False
PrecipBreedMoBase 0.020 -0.018 0.059 0.020 -0.020 0.059 0.314 False False
MinTempFallMig 0.034 0.004 0.062 0.034 0.010 0.064 0.020 True True
y_id 0.160 0.107 0.212 0.159 0.105 0.210 0.001 True True
scaleTarsus:Sex21 0.022 -0.028 0.070 0.022 -0.024 0.071 0.382 False False
scaleTarsus:PrevPrevTempWinMo 0.011 -0.020 0.044 0.012 -0.021 0.041 0.478 False False
scaleTarsus:PrevPrevPrecipWinMo 0.017 -0.011 0.046 0.017 -0.010 0.046 0.234 False False
scaleTarsus:PrevMinTempSpringMig -0.017 -0.047 0.013 -0.017 -0.047 0.011 0.284 False False
scaleTarsus:PrevAvgTempPreBreedBase 0.013 -0.020 0.046 0.014 -0.018 0.046 0.426 False False
scaleTarsus:PrevPrecipPreBreedBase 0.012 -0.017 0.041 0.013 -0.018 0.041 0.404 False False
scaleTarsus:PrevTempBreedBase -0.016 -0.058 0.028 -0.016 -0.055 0.027 0.470 False False
scaleTarsus:PrevPrecipBreedBase -0.010 -0.049 0.029 -0.009 -0.045 0.028 0.620 False False
scaleTarsus:PrevMinTempFallMig 0.008 -0.019 0.035 0.008 -0.020 0.033 0.588 False False
scaleTarsus:PrevTempWinMo -0.013 -0.046 0.021 -0.012 -0.044 0.022 0.464 False False
scaleTarsus:PrevPrecipWinMo -0.021 -0.048 0.007 -0.020 -0.044 0.007 0.146 False False
scaleTarsus:MinTempSpringMig 0.007 -0.021 0.034 0.006 -0.022 0.033 0.624 False False
scaleTarsus:AvgTempPreBreedBase -0.010 -0.046 0.025 -0.010 -0.044 0.026 0.594 False False
scaleTarsus:PrecipPreBreedBase 0.017 -0.016 0.050 0.017 -0.016 0.047 0.310 False False
scaleTarsus:TempBreedMoBase 0.020 -0.017 0.058 0.021 -0.017 0.058 0.274 False False
scaleTarsus:PrecipBreedMoBase 0.010 -0.030 0.048 0.010 -0.030 0.048 0.568 False False
scaleTarsus:MinTempFallMig -0.024 -0.052 0.005 -0.023 -0.050 0.006 0.118 False False
scaleTarsus:y_id -0.010 -0.060 0.041 -0.011 -0.060 0.040 0.684 False False
Helper for the claims of absence¶
For a list of climate terms, each claim below prints, per filtering reading, the estimate and both intervals of each main effect and tarsus interaction, and whether any interval excludes 0.
In [10]:
def absence(cls, terms, label):
rows = []
for r in RUN:
s = prelim[(r, cls)]["summary"]
for t in terms:
for name in (t, f"scaleTarsus:{t}"):
q = s.loc[name]
rows.append({"reading": r, "term": name, "mean": q.mean_pooled, "q2.5": q["q2.5_pooled"],
"q97.5": q["q97.5_pooled"], "hpd_lo_c4": q.hpd_lo_c4, "hpd_hi_c4": q.hpd_hi_c4,
"pMCMC_c4": q.pMCMC_c4, "excludes 0 (i)": q.excl0_i, "excludes 0 (ii)": q.excl0_ii})
out = pd.DataFrame(rows)
print(f"{label} ({cls} preliminary model):")
print(out.round(4).to_string(index=False))
print("any interval excluding 0, by reading:",
out.groupby("reading")[["excludes 0 (i)", "excludes 0 (ii)"]].any().to_dict("index"))
return out
Seed readings and the Monte Carlo check for intervals that end near 0¶
The code sets no seed, so a 95% bound that ends within a few thousandths of 0 can fall on either
side with another seed. The answer of 2026-10-27 to the attempt's question puts two seed readings
into the criteria of s4-s6-spring-migration-temperature-none and
s6-autumn-ahy-recent-wintering-prebreeding-precipitation-none: (1) one run with a recorded seed,
the fits above; (2) a set of seeds, a term counting as excluding 0 when its interval excludes 0 in
most of them. The set here is ten seed sets, one generator per seed set and chain, the same for
every model; "most" is more than 5 of 10. Each preliminary model is refitted with all ten, under
every filtering reading. seed_check prints the ten intervals of one term; seed_majority counts,
for a list of terms, the seed sets in which each interval excludes 0.
In [11]:
N_SETS = 10
t0 = time.time()
seed_sets = {}
for r in RUN:
for cls in w.CLASSES:
X, names, y, groups = w.prelim_design(prep[(r, "prelim", cls)], w.PRELIM_TERMS[cls])
seed_sets[(r, cls)] = [w.glmm_summary([w.gibbs_glmm(X, y, groups, w.SPECIES, Ainv, np.random.default_rng([99, seed, j]))
for j in range(4)], names) for seed in range(N_SETS)]
print(f"{sum(len(v) for v in seed_sets.values())} refits of 4 chains in {time.time() - t0:.0f} s")
def seed_check(r, cls, term):
out = pd.DataFrame([{"seed set": seed, "q2.5": s.loc[term, "q2.5_pooled"], "q97.5": s.loc[term, "q97.5_pooled"],
"hpd_lo_c4": s.loc[term, "hpd_lo_c4"], "hpd_hi_c4": s.loc[term, "hpd_hi_c4"],
"excludes 0 (i)": s.loc[term, "excl0_i"], "excludes 0 (ii)": s.loc[term, "excl0_ii"]}
for seed, s in enumerate(seed_sets[(r, cls)])])
print(f"{term}, {cls}, reading ({r}), {N_SETS} other seed sets:")
print(out.round(4).to_string(index=False))
print(f"interval (i) excludes 0 in {int(out['excludes 0 (i)'].sum())} of {N_SETS}; "
f"interval (ii) in {int(out['excludes 0 (ii)'].sum())} of {N_SETS}")
return out
def seed_majority(cls, terms, label):
rows = []
for r in RUN:
for t in terms:
for name in (t, f"scaleTarsus:{t}"):
ss = seed_sets[(r, cls)]
n_i = sum(bool(s.loc[name, "excl0_i"]) for s in ss)
n_ii = sum(bool(s.loc[name, "excl0_ii"]) for s in ss)
rows.append({"reading": r, "term": name,
"q97.5 range": f"{min(s.loc[name, 'q97.5_pooled'] for s in ss):.4f} to {max(s.loc[name, 'q97.5_pooled'] for s in ss):.4f}",
"q2.5 range": f"{min(s.loc[name, 'q2.5_pooled'] for s in ss):.4f} to {max(s.loc[name, 'q2.5_pooled'] for s in ss):.4f}",
"sets excluding 0 (i)": n_i, "sets excluding 0 (ii)": n_ii,
"excludes 0 in most (i)": n_i > N_SETS / 2, "excludes 0 in most (ii)": n_ii > N_SETS / 2})
out = pd.DataFrame(rows)
print(f"{label} ({cls} preliminary model), seed reading (2), {N_SETS} seed sets:")
print(out.to_string(index=False))
print("any term excluding 0 in most seed sets, by reading:",
out.groupby("reading")[["excludes 0 in most (i)", "excludes 0 in most (ii)"]].any().to_dict("index"))
return out
def seed_count(pairs, label):
"""For given (class, term) pairs, the seed sets in which each interval excludes 0, per filtering
reading, with the range of the bounds and of the estimate over the ten sets."""
rows = []
for r in RUN:
for cls, name in pairs:
ss = seed_sets[(r, cls)]
col = lambda k: [s.loc[name, k] for s in ss]
rows.append({"reading": r, "class": cls, "term": name,
"mean range": f"{min(col('mean_pooled')):.4f} to {max(col('mean_pooled')):.4f}",
"q2.5 range": f"{min(col('q2.5_pooled')):.4f} to {max(col('q2.5_pooled')):.4f}",
"q97.5 range": f"{min(col('q97.5_pooled')):.4f} to {max(col('q97.5_pooled')):.4f}",
"hpd_lo_c4 range": f"{min(col('hpd_lo_c4')):.4f} to {max(col('hpd_lo_c4')):.4f}",
"hpd_hi_c4 range": f"{min(col('hpd_hi_c4')):.4f} to {max(col('hpd_hi_c4')):.4f}",
"sets excluding 0 (i)": int(sum(col("excl0_i"))), "sets excluding 0 (ii)": int(sum(col("excl0_ii")))})
print(f"{label}, {N_SETS} other seed sets:")
print(pd.DataFrame(rows).to_string(index=False))
90 refits of 4 chains in 500 s
s4-s6-spring-migration-temperature-none¶
Spring migration temperature has no significant effect on wing length, or interaction with body size, in any of the three preliminary phylogenetic mixed models, and so is left out of the main models.
Terms: MinTempSpringMig (all three models), MinTempNextSpringMig (spring AHY) and PrevMinTempSpringMig (autumn AHY), with their tarsus interactions.
In [12]:
spring_mig = pd.concat([
absence("HY", ["MinTempSpringMig"], "Spring migration temperature"),
absence("spring AHY", ["MinTempSpringMig", "MinTempNextSpringMig"], "Spring migration temperature"),
absence("autumn AHY", ["PrevMinTempSpringMig", "MinTempSpringMig"], "Spring migration temperature")])
print("\nAny of the", len(spring_mig) // len(RUN), "intervals per reading excludes 0:",
spring_mig.groupby("reading")[["excludes 0 (i)", "excludes 0 (ii)"]].any().to_dict("index"))
Spring migration temperature (HY preliminary model):
reading term mean q2.5 q97.5 hpd_lo_c4 hpd_hi_c4 pMCMC_c4 excludes 0 (i) excludes 0 (ii)
a MinTempSpringMig -0.0110 -0.0238 0.0014 -0.0243 0.0007 0.068 False False
a scaleTarsus:MinTempSpringMig 0.0046 -0.0082 0.0167 -0.0079 0.0163 0.406 False False
c MinTempSpringMig -0.0102 -0.0228 0.0027 -0.0224 0.0020 0.120 False False
c scaleTarsus:MinTempSpringMig 0.0015 -0.0109 0.0138 -0.0108 0.0135 0.822 False False
d MinTempSpringMig -0.0103 -0.0230 0.0024 -0.0217 0.0022 0.102 False False
d scaleTarsus:MinTempSpringMig -0.0010 -0.0134 0.0115 -0.0123 0.0117 0.946 False False
any interval excluding 0, by reading: {'a': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}, 'c': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}, 'd': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}}
Spring migration temperature (spring AHY preliminary model):
reading term mean q2.5 q97.5 hpd_lo_c4 hpd_hi_c4 pMCMC_c4 excludes 0 (i) excludes 0 (ii)
a MinTempSpringMig -0.0061 -0.0214 0.0085 -0.0210 0.0086 0.462 False False
a scaleTarsus:MinTempSpringMig 0.0099 -0.0052 0.0250 -0.0049 0.0249 0.214 False False
a MinTempNextSpringMig -0.0134 -0.0275 0.0007 -0.0275 -0.0002 0.056 False True
a scaleTarsus:MinTempNextSpringMig 0.0117 -0.0027 0.0260 -0.0037 0.0255 0.118 False False
c MinTempSpringMig -0.0058 -0.0208 0.0093 -0.0206 0.0099 0.418 False False
c scaleTarsus:MinTempSpringMig 0.0090 -0.0057 0.0236 -0.0062 0.0231 0.252 False False
c MinTempNextSpringMig -0.0131 -0.0275 0.0015 -0.0292 -0.0005 0.068 False True
c scaleTarsus:MinTempNextSpringMig 0.0114 -0.0032 0.0253 -0.0026 0.0254 0.110 False False
d MinTempSpringMig -0.0062 -0.0212 0.0087 -0.0219 0.0078 0.452 False False
d scaleTarsus:MinTempSpringMig 0.0116 -0.0036 0.0261 -0.0044 0.0256 0.156 False False
d MinTempNextSpringMig -0.0136 -0.0279 0.0011 -0.0269 0.0019 0.084 False False
d scaleTarsus:MinTempNextSpringMig 0.0113 -0.0027 0.0256 -0.0030 0.0252 0.134 False False
any interval excluding 0, by reading: {'a': {'excludes 0 (i)': False, 'excludes 0 (ii)': True}, 'c': {'excludes 0 (i)': False, 'excludes 0 (ii)': True}, 'd': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}}
Spring migration temperature (autumn AHY preliminary model):
reading term mean q2.5 q97.5 hpd_lo_c4 hpd_hi_c4 pMCMC_c4 excludes 0 (i) excludes 0 (ii)
a PrevMinTempSpringMig -0.0103 -0.0410 0.0194 -0.0402 0.0204 0.496 False False
a scaleTarsus:PrevMinTempSpringMig -0.0167 -0.0474 0.0133 -0.0468 0.0111 0.284 False False
a MinTempSpringMig -0.0023 -0.0305 0.0260 -0.0307 0.0250 0.896 False False
a scaleTarsus:MinTempSpringMig 0.0065 -0.0213 0.0341 -0.0222 0.0333 0.624 False False
c PrevMinTempSpringMig -0.0094 -0.0400 0.0224 -0.0401 0.0227 0.538 False False
c scaleTarsus:PrevMinTempSpringMig -0.0168 -0.0466 0.0133 -0.0435 0.0149 0.290 False False
c MinTempSpringMig -0.0048 -0.0325 0.0231 -0.0320 0.0238 0.712 False False
c scaleTarsus:MinTempSpringMig 0.0098 -0.0177 0.0373 -0.0175 0.0348 0.454 False False
d PrevMinTempSpringMig -0.0083 -0.0391 0.0228 -0.0411 0.0214 0.646 False False
d scaleTarsus:PrevMinTempSpringMig -0.0154 -0.0455 0.0148 -0.0464 0.0120 0.306 False False
d MinTempSpringMig -0.0030 -0.0326 0.0251 -0.0288 0.0274 0.870 False False
d scaleTarsus:MinTempSpringMig 0.0087 -0.0187 0.0372 -0.0180 0.0343 0.494 False False
any interval excluding 0, by reading: {'a': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}, 'c': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}, 'd': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}}
Any of the 10 intervals per reading excludes 0: {'a': {'excludes 0 (i)': False, 'excludes 0 (ii)': True}, 'c': {'excludes 0 (i)': False, 'excludes 0 (ii)': True}, 'd': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}}
In [13]:
_ = seed_check("a", "spring AHY", "MinTempNextSpringMig")
_ = seed_check("c", "spring AHY", "MinTempNextSpringMig")
spring_mig_seeds = pd.concat([
seed_majority("HY", ["MinTempSpringMig"], "Spring migration temperature"),
seed_majority("spring AHY", ["MinTempSpringMig", "MinTempNextSpringMig"], "Spring migration temperature"),
seed_majority("autumn AHY", ["PrevMinTempSpringMig", "MinTempSpringMig"], "Spring migration temperature")])
print("\nSeed reading (2): any of the", len(spring_mig_seeds) // len(RUN), "terms per filtering reading excludes 0 in most seed sets:",
spring_mig_seeds.groupby("reading")[["excludes 0 in most (i)", "excludes 0 in most (ii)"]].any().to_dict("index"))
MinTempNextSpringMig, spring AHY, reading (a), 10 other seed sets:
seed set q2.5 q97.5 hpd_lo_c4 hpd_hi_c4 excludes 0 (i) excludes 0 (ii)
0 -0.0283 0.0009 -0.0284 0.0001 False False
1 -0.0279 0.0015 -0.0284 0.0013 False False
2 -0.0276 0.0010 -0.0267 0.0012 False False
3 -0.0275 0.0007 -0.0278 -0.0000 False True
4 -0.0281 0.0007 -0.0282 0.0003 False False
5 -0.0274 0.0009 -0.0271 0.0012 False False
6 -0.0275 0.0012 -0.0281 0.0010 False False
7 -0.0281 0.0011 -0.0261 0.0022 False False
8 -0.0276 0.0012 -0.0280 0.0003 False False
9 -0.0275 0.0009 -0.0275 0.0005 False False
interval (i) excludes 0 in 0 of 10; interval (ii) in 1 of 10
MinTempNextSpringMig, spring AHY, reading (c), 10 other seed sets:
seed set q2.5 q97.5 hpd_lo_c4 hpd_hi_c4 excludes 0 (i) excludes 0 (ii)
0 -0.0282 0.0012 -0.0283 0.0004 False False
1 -0.0278 0.0017 -0.0283 0.0015 False False
2 -0.0275 0.0013 -0.0265 0.0015 False False
3 -0.0274 0.0009 -0.0276 0.0002 False False
4 -0.0279 0.0009 -0.0281 0.0005 False False
5 -0.0273 0.0011 -0.0270 0.0014 False False
6 -0.0273 0.0014 -0.0280 0.0012 False False
7 -0.0280 0.0013 -0.0259 0.0024 False False
8 -0.0275 0.0014 -0.0278 0.0006 False False
9 -0.0274 0.0011 -0.0274 0.0008 False False
interval (i) excludes 0 in 0 of 10; interval (ii) in 0 of 10
Spring migration temperature (HY preliminary model), seed reading (2), 10 seed sets:
reading term q97.5 range q2.5 range sets excluding 0 (i) sets excluding 0 (ii) excludes 0 in most (i) excludes 0 in most (ii)
a MinTempSpringMig 0.0016 to 0.0020 -0.0240 to -0.0235 0 0 False False
a scaleTarsus:MinTempSpringMig 0.0166 to 0.0174 -0.0083 to -0.0078 0 0 False False
c MinTempSpringMig 0.0022 to 0.0026 -0.0232 to -0.0227 0 0 False False
c scaleTarsus:MinTempSpringMig 0.0132 to 0.0141 -0.0115 to -0.0109 0 0 False False
d MinTempSpringMig 0.0023 to 0.0028 -0.0231 to -0.0226 0 0 False False
d scaleTarsus:MinTempSpringMig 0.0109 to 0.0117 -0.0139 to -0.0134 0 0 False False
any term excluding 0 in most seed sets, by reading: {'a': {'excludes 0 in most (i)': False, 'excludes 0 in most (ii)': False}, 'c': {'excludes 0 in most (i)': False, 'excludes 0 in most (ii)': False}, 'd': {'excludes 0 in most (i)': False, 'excludes 0 in most (ii)': False}}
Spring migration temperature (spring AHY preliminary model), seed reading (2), 10 seed sets:
reading term q97.5 range q2.5 range sets excluding 0 (i) sets excluding 0 (ii) excludes 0 in most (i) excludes 0 in most (ii)
a MinTempSpringMig 0.0085 to 0.0095 -0.0220 to -0.0208 0 0 False False
a scaleTarsus:MinTempSpringMig 0.0240 to 0.0250 -0.0056 to -0.0046 0 0 False False
a MinTempNextSpringMig 0.0007 to 0.0015 -0.0283 to -0.0274 0 1 False False
a scaleTarsus:MinTempNextSpringMig 0.0257 to 0.0264 -0.0030 to -0.0021 0 0 False False
c MinTempSpringMig 0.0087 to 0.0097 -0.0218 to -0.0207 0 0 False False
c scaleTarsus:MinTempSpringMig 0.0233 to 0.0243 -0.0064 to -0.0054 0 0 False False
c MinTempNextSpringMig 0.0009 to 0.0017 -0.0282 to -0.0273 0 0 False False
c scaleTarsus:MinTempNextSpringMig 0.0254 to 0.0262 -0.0034 to -0.0024 0 0 False False
d MinTempSpringMig 0.0083 to 0.0093 -0.0222 to -0.0211 0 0 False False
d scaleTarsus:MinTempSpringMig 0.0259 to 0.0270 -0.0036 to -0.0025 0 0 False False
d MinTempNextSpringMig 0.0006 to 0.0014 -0.0285 to -0.0276 0 2 False False
d scaleTarsus:MinTempNextSpringMig 0.0256 to 0.0263 -0.0033 to -0.0023 0 0 False False
any term excluding 0 in most seed sets, by reading: {'a': {'excludes 0 in most (i)': False, 'excludes 0 in most (ii)': False}, 'c': {'excludes 0 in most (i)': False, 'excludes 0 in most (ii)': False}, 'd': {'excludes 0 in most (i)': False, 'excludes 0 in most (ii)': False}}
Spring migration temperature (autumn AHY preliminary model), seed reading (2), 10 seed sets:
reading term q97.5 range q2.5 range sets excluding 0 (i) sets excluding 0 (ii) excludes 0 in most (i) excludes 0 in most (ii)
a PrevMinTempSpringMig 0.0200 to 0.0221 -0.0415 to -0.0394 0 0 False False
a scaleTarsus:PrevMinTempSpringMig 0.0127 to 0.0150 -0.0472 to -0.0454 0 0 False False
a MinTempSpringMig 0.0249 to 0.0273 -0.0308 to -0.0289 0 0 False False
a scaleTarsus:MinTempSpringMig 0.0331 to 0.0350 -0.0224 to -0.0202 0 0 False False
c PrevMinTempSpringMig 0.0204 to 0.0226 -0.0406 to -0.0385 0 0 False False
c scaleTarsus:PrevMinTempSpringMig 0.0126 to 0.0148 -0.0468 to -0.0450 0 0 False False
c MinTempSpringMig 0.0221 to 0.0246 -0.0331 to -0.0313 0 0 False False
c scaleTarsus:MinTempSpringMig 0.0364 to 0.0384 -0.0186 to -0.0164 0 0 False False
d PrevMinTempSpringMig 0.0216 to 0.0239 -0.0400 to -0.0380 0 0 False False
d scaleTarsus:PrevMinTempSpringMig 0.0138 to 0.0159 -0.0463 to -0.0443 0 0 False False
d MinTempSpringMig 0.0235 to 0.0260 -0.0323 to -0.0304 0 0 False False
d scaleTarsus:MinTempSpringMig 0.0352 to 0.0374 -0.0193 to -0.0174 0 0 False False
any term excluding 0 in most seed sets, by reading: {'a': {'excludes 0 in most (i)': False, 'excludes 0 in most (ii)': False}, 'c': {'excludes 0 in most (i)': False, 'excludes 0 in most (ii)': False}, 'd': {'excludes 0 in most (i)': False, 'excludes 0 in most (ii)': False}}
Seed reading (2): any of the 10 terms per filtering reading excludes 0 in most seed sets: {'a': {'excludes 0 in most (i)': False, 'excludes 0 in most (ii)': False}, 'c': {'excludes 0 in most (i)': False, 'excludes 0 in most (ii)': False}, 'd': {'excludes 0 in most (i)': False, 'excludes 0 in most (ii)': False}}
Result: blocked, methods underspecified: the two seed readings give different outcomes. Seed reading (1), the recorded run: every interval of the spring migration temperatures and their tarsus interactions includes 0 under interval reading (i), in all three models and filtering readings, but under reading (ii) the chain-4 HPD interval of MinTempNextSpringMig in the spring AHY model ends at -0.0002 (filtering reading a) and -0.0005 (reading c), excluding 0, so this reading gives not reproduced; its pooled interval is -0.0275 to 0.0007, pMCMC 0.056. Seed reading (2), ten other seed sets: no term's interval excludes 0 in most of them under either interval reading or any filtering reading; the most is MinTempNextSpringMig's HPD interval in 2 of 10 sets under reading (d) and 1 of 10 under (a), and the pooled intervals exclude 0 in none. So this reading gives reproduced. The answer of 2026-10-27 blocks a claim whose seed readings differ, since the seed is a choice of the attempt.
s5-s6-ahy-prebreeding-wintering-temperature-none¶
Pre-breeding and wintering season temperatures have no significant association with AHY wing length in the preliminary models of spring and autumn AHY birds.
In [14]:
ahy_temp = pd.concat([
absence("spring AHY", ["TempWinMo", "PrevTempWinMo", "AvgTempPreBreedBase"], "Wintering and pre-breeding temperature"),
absence("autumn AHY", ["PrevTempWinMo", "PrevPrevTempWinMo", "AvgTempPreBreedBase", "PrevAvgTempPreBreedBase"],
"Wintering and pre-breeding temperature")])
print("\nAny interval excludes 0:", ahy_temp.groupby("reading")[["excludes 0 (i)", "excludes 0 (ii)"]].any().to_dict("index"))
print("Terms whose interval excludes 0:")
print(ahy_temp[ahy_temp["excludes 0 (i)"] | ahy_temp["excludes 0 (ii)"]][["reading", "term", "mean", "q2.5", "q97.5", "hpd_lo_c4", "hpd_hi_c4"]].round(4).to_string(index=False))
Wintering and pre-breeding temperature (spring AHY preliminary model):
reading term mean q2.5 q97.5 hpd_lo_c4 hpd_hi_c4 pMCMC_c4 excludes 0 (i) excludes 0 (ii)
a TempWinMo 0.0088 -0.0071 0.0243 -0.0064 0.0247 0.274 False False
a scaleTarsus:TempWinMo -0.0026 -0.0187 0.0129 -0.0195 0.0114 0.716 False False
a PrevTempWinMo -0.0146 -0.0306 0.0017 -0.0299 0.0025 0.080 False False
a scaleTarsus:PrevTempWinMo -0.0013 -0.0169 0.0144 -0.0171 0.0135 0.884 False False
a AvgTempPreBreedBase 0.0130 -0.0022 0.0285 -0.0028 0.0266 0.100 False False
a scaleTarsus:AvgTempPreBreedBase -0.0081 -0.0234 0.0071 -0.0237 0.0066 0.290 False False
c TempWinMo 0.0088 -0.0069 0.0250 -0.0062 0.0251 0.286 False False
c scaleTarsus:TempWinMo -0.0026 -0.0183 0.0131 -0.0193 0.0113 0.760 False False
c PrevTempWinMo -0.0145 -0.0306 0.0017 -0.0304 0.0010 0.072 False False
c scaleTarsus:PrevTempWinMo -0.0016 -0.0166 0.0142 -0.0141 0.0153 0.806 False False
c AvgTempPreBreedBase 0.0131 -0.0033 0.0289 -0.0032 0.0292 0.126 False False
c scaleTarsus:AvgTempPreBreedBase -0.0075 -0.0232 0.0080 -0.0241 0.0070 0.404 False False
d TempWinMo 0.0100 -0.0057 0.0258 -0.0056 0.0256 0.224 False False
d scaleTarsus:TempWinMo -0.0071 -0.0231 0.0094 -0.0231 0.0090 0.398 False False
d PrevTempWinMo -0.0136 -0.0296 0.0028 -0.0276 0.0033 0.104 False False
d scaleTarsus:PrevTempWinMo -0.0034 -0.0187 0.0122 -0.0181 0.0121 0.660 False False
d AvgTempPreBreedBase 0.0125 -0.0032 0.0282 -0.0032 0.0267 0.106 False False
d scaleTarsus:AvgTempPreBreedBase -0.0035 -0.0189 0.0119 -0.0187 0.0116 0.630 False False
any interval excluding 0, by reading: {'a': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}, 'c': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}, 'd': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}}
Wintering and pre-breeding temperature (autumn AHY preliminary model):
reading term mean q2.5 q97.5 hpd_lo_c4 hpd_hi_c4 pMCMC_c4 excludes 0 (i) excludes 0 (ii)
a PrevTempWinMo 0.0160 -0.0158 0.0484 -0.0147 0.0483 0.336 False False
a scaleTarsus:PrevTempWinMo -0.0127 -0.0460 0.0211 -0.0440 0.0217 0.464 False False
a PrevPrevTempWinMo 0.0277 -0.0049 0.0602 -0.0038 0.0595 0.086 False False
a scaleTarsus:PrevPrevTempWinMo 0.0115 -0.0204 0.0438 -0.0208 0.0415 0.478 False False
a AvgTempPreBreedBase 0.0005 -0.0357 0.0373 -0.0366 0.0340 0.996 False False
a scaleTarsus:AvgTempPreBreedBase -0.0098 -0.0457 0.0250 -0.0444 0.0261 0.594 False False
a PrevAvgTempPreBreedBase 0.0198 -0.0155 0.0522 -0.0190 0.0468 0.262 False False
a scaleTarsus:PrevAvgTempPreBreedBase 0.0133 -0.0200 0.0459 -0.0179 0.0459 0.426 False False
c PrevTempWinMo 0.0149 -0.0168 0.0482 -0.0165 0.0491 0.366 False False
c scaleTarsus:PrevTempWinMo -0.0115 -0.0454 0.0219 -0.0483 0.0185 0.474 False False
c PrevPrevTempWinMo 0.0256 -0.0072 0.0579 -0.0032 0.0617 0.096 False False
c scaleTarsus:PrevPrevTempWinMo 0.0128 -0.0189 0.0454 -0.0221 0.0431 0.486 False False
c AvgTempPreBreedBase 0.0028 -0.0339 0.0389 -0.0341 0.0375 0.830 False False
c scaleTarsus:AvgTempPreBreedBase -0.0135 -0.0479 0.0225 -0.0488 0.0202 0.448 False False
c PrevAvgTempPreBreedBase 0.0191 -0.0158 0.0536 -0.0119 0.0548 0.284 False False
c scaleTarsus:PrevAvgTempPreBreedBase 0.0130 -0.0199 0.0466 -0.0204 0.0455 0.460 False False
d PrevTempWinMo 0.0155 -0.0172 0.0487 -0.0174 0.0458 0.388 False False
d scaleTarsus:PrevTempWinMo -0.0077 -0.0406 0.0248 -0.0387 0.0262 0.672 False False
d PrevPrevTempWinMo 0.0264 -0.0059 0.0595 -0.0065 0.0601 0.114 False False
d scaleTarsus:PrevPrevTempWinMo 0.0123 -0.0198 0.0453 -0.0196 0.0449 0.404 False False
d AvgTempPreBreedBase 0.0057 -0.0312 0.0438 -0.0352 0.0401 0.752 False False
d scaleTarsus:AvgTempPreBreedBase -0.0034 -0.0404 0.0339 -0.0382 0.0347 0.866 False False
d PrevAvgTempPreBreedBase 0.0225 -0.0114 0.0578 -0.0131 0.0572 0.204 False False
d scaleTarsus:PrevAvgTempPreBreedBase 0.0127 -0.0207 0.0459 -0.0208 0.0453 0.466 False False
any interval excluding 0, by reading: {'a': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}, 'c': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}, 'd': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}}
Any interval excludes 0: {'a': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}, 'c': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}, 'd': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}}
Terms whose interval excludes 0:
Empty DataFrame
Columns: [reading, term, mean, q2.5, q97.5, hpd_lo_c4, hpd_hi_c4]
Index: []
Result: reproduced. In the spring and autumn AHY preliminary models, every interval of the wintering and pre-breeding temperatures (TempWinMo, PrevTempWinMo, PrevPrevTempWinMo, AvgTempPreBreedBase, PrevAvgTempPreBreedBase) and their tarsus interactions includes 0, under both interval readings and all three filtering readings.
s6-autumn-ahy-recent-breeding-temperature-none¶
In the autumn AHY preliminary model, temperature of the most recent breeding season is not significantly associated with wing length.
In [15]:
_ = absence("autumn AHY", ["TempBreedMoBase"], "Most recent breeding temperature")
Most recent breeding temperature (autumn AHY preliminary model):
reading term mean q2.5 q97.5 hpd_lo_c4 hpd_hi_c4 pMCMC_c4 excludes 0 (i) excludes 0 (ii)
a TempBreedMoBase 0.0224 -0.0170 0.0614 -0.0152 0.0624 0.266 False False
a scaleTarsus:TempBreedMoBase 0.0201 -0.0171 0.0579 -0.0168 0.0576 0.274 False False
c TempBreedMoBase 0.0205 -0.0179 0.0598 -0.0191 0.0576 0.304 False False
c scaleTarsus:TempBreedMoBase 0.0234 -0.0158 0.0615 -0.0160 0.0595 0.244 False False
d TempBreedMoBase 0.0163 -0.0236 0.0565 -0.0295 0.0522 0.396 False False
d scaleTarsus:TempBreedMoBase 0.0224 -0.0172 0.0622 -0.0153 0.0641 0.264 False False
any interval excluding 0, by reading: {'a': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}, 'c': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}, 'd': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}}
Result: reproduced. TempBreedMoBase (0.022, -0.017 to 0.061 under reading a) and its tarsus interaction include 0 under both interval readings and all filtering readings.
s4-hy-autumn-migration-temperature-none¶
In the HY preliminary model, autumn migration temperature has no significant association with wing length.
In [16]:
_ = absence("HY", ["MinTempFallMig"], "Autumn migration temperature")
Autumn migration temperature (HY preliminary model):
reading term mean q2.5 q97.5 hpd_lo_c4 hpd_hi_c4 pMCMC_c4 excludes 0 (i) excludes 0 (ii)
a MinTempFallMig -0.0069 -0.0197 0.0061 -0.0204 0.0055 0.294 False False
a scaleTarsus:MinTempFallMig -0.0001 -0.0128 0.0127 -0.0133 0.0128 0.998 False False
c MinTempFallMig -0.0069 -0.0193 0.0055 -0.0180 0.0067 0.282 False False
c scaleTarsus:MinTempFallMig 0.0006 -0.0118 0.0135 -0.0132 0.0132 0.930 False False
d MinTempFallMig -0.0068 -0.0200 0.0059 -0.0201 0.0059 0.312 False False
d scaleTarsus:MinTempFallMig -0.0008 -0.0137 0.0124 -0.0137 0.0138 0.932 False False
any interval excluding 0, by reading: {'a': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}, 'c': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}, 'd': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}}
Result: reproduced. In the HY model, MinTempFallMig (-0.007, -0.020 to 0.006 under reading a) and its tarsus interaction include 0 under both interval readings and all filtering readings.
s6-autumn-ahy-recent-wintering-prebreeding-precipitation-none¶
In the autumn AHY preliminary model, precipitation of the most recent wintering season and pre-breeding period is not associated with wing length.
Terms as the criterion gives them: PrevPrecipWinMo and PrecipPreBreedBase.
In [17]:
_ = absence("autumn AHY", ["PrevPrecipWinMo", "PrecipPreBreedBase"], "Most recent wintering and pre-breeding precipitation")
Most recent wintering and pre-breeding precipitation (autumn AHY preliminary model):
reading term mean q2.5 q97.5 hpd_lo_c4 hpd_hi_c4 pMCMC_c4 excludes 0 (i) excludes 0 (ii)
a PrevPrecipWinMo -0.0199 -0.0473 0.0074 -0.0507 0.0043 0.144 False False
a scaleTarsus:PrevPrecipWinMo -0.0211 -0.0477 0.0066 -0.0442 0.0067 0.146 False False
a PrecipPreBreedBase 0.0035 -0.0286 0.0343 -0.0279 0.0351 0.804 False False
a scaleTarsus:PrecipPreBreedBase 0.0166 -0.0161 0.0495 -0.0156 0.0467 0.310 False False
c PrevPrecipWinMo -0.0194 -0.0470 0.0073 -0.0459 0.0079 0.156 False False
c scaleTarsus:PrevPrecipWinMo -0.0208 -0.0485 0.0057 -0.0476 0.0076 0.144 False False
c PrecipPreBreedBase 0.0028 -0.0289 0.0338 -0.0277 0.0347 0.846 False False
c scaleTarsus:PrecipPreBreedBase 0.0176 -0.0159 0.0510 -0.0165 0.0477 0.292 False False
d PrevPrecipWinMo -0.0181 -0.0456 0.0100 -0.0453 0.0090 0.204 False False
d scaleTarsus:PrevPrecipWinMo -0.0262 -0.0532 0.0009 -0.0532 -0.0003 0.054 False True
d PrecipPreBreedBase 0.0004 -0.0307 0.0325 -0.0318 0.0317 0.992 False False
d scaleTarsus:PrecipPreBreedBase 0.0136 -0.0189 0.0463 -0.0190 0.0451 0.384 False False
any interval excluding 0, by reading: {'a': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}, 'c': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}, 'd': {'excludes 0 (i)': False, 'excludes 0 (ii)': True}}
In [18]:
_ = seed_check("d", "autumn AHY", "scaleTarsus:PrevPrecipWinMo")
_ = seed_majority("autumn AHY", ["PrevPrecipWinMo", "PrecipPreBreedBase"], "Most recent wintering and pre-breeding precipitation")
scaleTarsus:PrevPrecipWinMo, autumn AHY, reading (d), 10 other seed sets:
seed set q2.5 q97.5 hpd_lo_c4 hpd_hi_c4 excludes 0 (i) excludes 0 (ii)
0 -0.0521 0.0019 -0.0530 0.0021 False False
1 -0.0533 0.0012 -0.0532 0.0015 False False
2 -0.0524 0.0005 -0.0521 0.0005 False False
3 -0.0531 0.0018 -0.0513 0.0022 False False
4 -0.0532 0.0007 -0.0531 0.0040 False False
5 -0.0539 0.0008 -0.0536 0.0007 False False
6 -0.0527 0.0015 -0.0528 0.0012 False False
7 -0.0531 0.0011 -0.0555 -0.0013 False True
8 -0.0532 0.0008 -0.0505 0.0003 False False
9 -0.0528 0.0018 -0.0507 0.0012 False False
interval (i) excludes 0 in 0 of 10; interval (ii) in 1 of 10
Most recent wintering and pre-breeding precipitation (autumn AHY preliminary model), seed reading (2), 10 seed sets:
reading term q97.5 range q2.5 range sets excluding 0 (i) sets excluding 0 (ii) excludes 0 in most (i) excludes 0 in most (ii)
a PrevPrecipWinMo 0.0064 to 0.0089 -0.0484 to -0.0469 0 0 False False
a scaleTarsus:PrevPrecipWinMo 0.0054 to 0.0070 -0.0490 to -0.0472 0 0 False False
a PrecipPreBreedBase 0.0340 to 0.0359 -0.0296 to -0.0270 0 0 False False
a scaleTarsus:PrecipPreBreedBase 0.0472 to 0.0502 -0.0166 to -0.0147 0 0 False False
c PrevPrecipWinMo 0.0062 to 0.0087 -0.0482 to -0.0467 0 0 False False
c scaleTarsus:PrevPrecipWinMo 0.0050 to 0.0067 -0.0489 to -0.0471 0 0 False False
c PrecipPreBreedBase 0.0325 to 0.0343 -0.0307 to -0.0282 0 0 False False
c scaleTarsus:PrecipPreBreedBase 0.0483 to 0.0512 -0.0149 to -0.0130 0 0 False False
d PrevPrecipWinMo 0.0082 to 0.0107 -0.0466 to -0.0451 0 0 False False
d scaleTarsus:PrevPrecipWinMo 0.0005 to 0.0019 -0.0539 to -0.0521 0 1 False False
d PrecipPreBreedBase 0.0313 to 0.0331 -0.0322 to -0.0297 0 0 False False
d scaleTarsus:PrecipPreBreedBase 0.0446 to 0.0474 -0.0190 to -0.0170 0 0 False False
any term excluding 0 in most seed sets, by reading: {'a': {'excludes 0 in most (i)': False, 'excludes 0 in most (ii)': False}, 'c': {'excludes 0 in most (i)': False, 'excludes 0 in most (ii)': False}, 'd': {'excludes 0 in most (i)': False, 'excludes 0 in most (ii)': False}}
Result: blocked, methods underspecified: the two seed readings give different outcomes. Seed reading (1), the recorded run: PrevPrecipWinMo, PrecipPreBreedBase and their tarsus interactions include 0 under both interval readings in filtering readings (a) and (c), and under interval reading (i) in (d); under reading (d) the chain-4 HPD interval of the tarsus by PrevPrecipWinMo interaction is -0.0532 to -0.0003 and excludes 0, so this reading gives not reproduced. Seed reading (2), ten other seed sets: none of the four terms' intervals excludes 0 in most of them under either interval reading or any filtering reading; the most is that interaction's HPD interval in 1 of 10 sets under reading (d), and the pooled intervals exclude 0 in none. So this reading gives reproduced, and the claim is blocked under the answer of 2026-10-27.
s6-autumn-ahy-recent-breeding-precipitation-none¶
In the autumn AHY preliminary model, precipitation of the most recent breeding season is not significant.
In [19]:
_ = absence("autumn AHY", ["PrecipBreedMoBase"], "Most recent breeding precipitation")
Most recent breeding precipitation (autumn AHY preliminary model):
reading term mean q2.5 q97.5 hpd_lo_c4 hpd_hi_c4 pMCMC_c4 excludes 0 (i) excludes 0 (ii)
a PrecipBreedMoBase 0.0202 -0.0180 0.0594 -0.0198 0.0592 0.314 False False
a scaleTarsus:PrecipBreedMoBase 0.0096 -0.0300 0.0475 -0.0297 0.0481 0.568 False False
c PrecipBreedMoBase 0.0200 -0.0167 0.0591 -0.0170 0.0580 0.298 False False
c scaleTarsus:PrecipBreedMoBase 0.0113 -0.0260 0.0517 -0.0252 0.0507 0.562 False False
d PrecipBreedMoBase 0.0177 -0.0203 0.0561 -0.0201 0.0567 0.344 False False
d scaleTarsus:PrecipBreedMoBase 0.0136 -0.0249 0.0519 -0.0263 0.0502 0.428 False False
any interval excluding 0, by reading: {'a': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}, 'c': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}, 'd': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}}
Result: reproduced. PrecipBreedMoBase (0.020, -0.018 to 0.059 under reading a) and its tarsus interaction include 0 under both interval readings and all filtering readings.
table3-autumn-ahy-prior-wintering-precipitation, preliminary part¶
The claim also needs the autumn AHY main model (autumn_ahy_main_a.ipynb, _c and _d). Here: PrevPrecipPreBreedBase
and its tarsus interaction in the autumn AHY preliminary model (Table S6 reports -0.00 to 0.06).
In [20]:
_ = absence("autumn AHY", ["PrevPrecipPreBreedBase"], "Pre-breeding precipitation before development")
Pre-breeding precipitation before development (autumn AHY preliminary model):
reading term mean q2.5 q97.5 hpd_lo_c4 hpd_hi_c4 pMCMC_c4 excludes 0 (i) excludes 0 (ii)
a PrevPrecipPreBreedBase 0.0268 -0.0018 0.0549 -0.0034 0.0522 0.068 False False
a scaleTarsus:PrevPrecipPreBreedBase 0.0118 -0.0172 0.0414 -0.0178 0.0405 0.404 False False
c PrevPrecipPreBreedBase 0.0261 -0.0025 0.0552 -0.0017 0.0552 0.090 False False
c scaleTarsus:PrevPrecipPreBreedBase 0.0130 -0.0159 0.0434 -0.0180 0.0404 0.374 False False
d PrevPrecipPreBreedBase 0.0244 -0.0041 0.0527 -0.0029 0.0555 0.108 False False
d scaleTarsus:PrevPrecipPreBreedBase 0.0139 -0.0147 0.0425 -0.0138 0.0418 0.326 False False
any interval excluding 0, by reading: {'a': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}, 'c': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}, 'd': {'excludes 0 (i)': False, 'excludes 0 (ii)': False}}
Preliminary part: PrevPrecipPreBreedBase (0.027, -0.002 to 0.055 under reading a; the paper -0.00 to
0.06) and its tarsus interaction include 0 under both interval readings and all filtering readings.
The main-model part and the outcome are in autumn_ahy_main_a.ipynb, _c and _d.
s7-low-phylogenetic-signal¶
Phylogenetic signal is low in all three preliminary models (posterior mode of lambda 0.02, 0.02 and 0.04), so the main models do not control for phylogeny.
Lambda from chain 4, species variance over the sum of species and residual variances; posterior
mode (MCMCglmm's posterior.mode) and 95% HPD. Reported HPDs: HY 0.01 to 0.21, spring AHY 0.01 to
0.16, autumn AHY 0.01 to 0.21.
In [21]:
reported = {"HY": (0.02, 0.01, 0.21), "spring AHY": (0.02, 0.01, 0.16), "autumn AHY": (0.04, 0.01, 0.21)}
rows = []
for r in RUN:
for cls in w.CLASSES:
lam = prelim[(r, cls)]["lambda"]
lo, hi = w.hpd(lam)
mode = w.posterior_mode(lam)
pm, plo, phi = reported[cls]
rows.append({"reading": r, "class": cls, "mode": mode, "mean": lam.mean(), "hpd_lo": lo, "hpd_hi": hi,
"paper mode": pm, "paper HPD": f"{plo} to {phi}", "mode in paper HPD": plo <= mode <= phi})
lam_table = pd.DataFrame(rows)
print(lam_table.round(3).to_string(index=False))
fig, axes = plt.subplots(1, 3, figsize=(11, 3))
for ax, cls in zip(axes, w.CLASSES):
lam = prelim[("a", cls)]["lambda"]
ax.hist(lam, bins=50, range=(0, 1))
ax.set_title(f"{cls}, reading (a)")
ax.set_xlabel("lambda (chain 4)")
print(f"{cls}: lambda draws range {lam.min():.4f} to {lam.max():.4f}; all inside the 0 to 1 axis")
plt.tight_layout()
plt.show()
reading class mode mean hpd_lo hpd_hi paper mode paper HPD mode in paper HPD
a HY 0.056 0.136 0.013 0.385 0.02 0.01 to 0.21 True
a spring AHY 0.062 0.120 0.013 0.303 0.02 0.01 to 0.16 True
a autumn AHY 0.087 0.177 0.012 0.481 0.04 0.01 to 0.21 True
c HY 0.047 0.123 0.011 0.369 0.02 0.01 to 0.21 True
c spring AHY 0.051 0.121 0.013 0.351 0.02 0.01 to 0.16 True
c autumn AHY 0.067 0.170 0.013 0.439 0.04 0.01 to 0.21 True
d HY 0.038 0.128 0.012 0.383 0.02 0.01 to 0.21 True
d spring AHY 0.042 0.125 0.009 0.322 0.02 0.01 to 0.16 True
d autumn AHY 0.089 0.177 0.012 0.453 0.04 0.01 to 0.21 True
HY: lambda draws range 0.0088 to 0.8618; all inside the 0 to 1 axis spring AHY: lambda draws range 0.0113 to 0.8137; all inside the 0 to 1 axis autumn AHY: lambda draws range 0.0103 to 0.9166; all inside the 0 to 1 axis
The mode of a posterior that piles up against 0 depends on the chain. As a check of that spread, the lambda mode of each of the four chains, reading (a), and the mode and HPD interval of the four chains pooled, every reading:
In [22]:
for cls in w.CLASSES:
modes = [w.posterior_mode(c["VCV"][:, 0] / c["VCV"].sum(axis=1)) for c in prelim[("a", cls)]["chains"]]
print(cls, "modes of chains 1 to 4:", np.round(modes, 3))
print("\nLambda from the four chains pooled (the criterion uses chain 4, as Script02 does):")
rows = []
for r in RUN:
for cls in w.CLASSES:
lam = np.concatenate([ch["VCV"][:, 0] / ch["VCV"].sum(axis=1) for ch in prelim[(r, cls)]["chains"]])
lo, hi = w.hpd(lam)
pm, plo, phi = reported[cls]
rows.append({"reading": r, "class": cls, "draws": len(lam), "mode": w.posterior_mode(lam), "median": np.median(lam),
"hpd_lo": lo, "hpd_hi": hi, "paper mode": pm, "paper HPD": f"{plo} to {phi}"})
print(pd.DataFrame(rows).round(3).to_string(index=False))
HY modes of chains 1 to 4: [0.047 0.042 0.054 0.056] spring AHY modes of chains 1 to 4: [0.047 0.045 0.041 0.062] autumn AHY modes of chains 1 to 4: [0.096 0.095 0.076 0.087] Lambda from the four chains pooled (the criterion uses chain 4, as Script02 does):
reading class draws mode median hpd_lo hpd_hi paper mode paper HPD
a HY 4000 0.054 0.100 0.011 0.382 0.02 0.01 to 0.21
a spring AHY 4000 0.040 0.091 0.011 0.335 0.02 0.01 to 0.16
a autumn AHY 4000 0.096 0.133 0.009 0.475 0.04 0.01 to 0.21
c HY 4000 0.045 0.088 0.011 0.348 0.02 0.01 to 0.21
c spring AHY 4000 0.066 0.091 0.010 0.343 0.02 0.01 to 0.16
c autumn AHY 4000 0.069 0.137 0.013 0.469 0.04 0.01 to 0.21
d HY 4000 0.044 0.085 0.009 0.370 0.02 0.01 to 0.21
d spring AHY 4000 0.057 0.093 0.012 0.341 0.02 0.01 to 0.16
d autumn AHY 4000 0.068 0.138 0.012 0.484 0.04 0.01 to 0.21
Result: reproduced. The posterior modes of lambda (chain 4) are 0.056, 0.062 and 0.087 (HY, spring AHY, autumn AHY; reading a), and 0.038 to 0.089 over all readings, each inside the paper's HPD interval (0.01 to 0.21, 0.01 to 0.16, 0.01 to 0.21). The modes are larger than the paper's 0.02, 0.02 and 0.04, and the HPD intervals found here are wider (upper bounds 0.30 to 0.48). The four chains pooled give the same picture, modes 0.040 to 0.096 and HPD upper bounds 0.335 to 0.484, so the difference from the paper is not the noise of chain 4: the Gibbs sampler here gives a wider posterior of the species variance than the paper's MCMCglmm runs. Lambda is low by the paper's own measure, which names no level for "low".
s4-s6-no-multicollinearity¶
The predictors of the preliminary models are not collinear; every variance inflation factor is below 10 (largest 3.34 for HY, 2.43 for spring AHY, 6.08 for autumn AHY).
In [23]:
paper_vif = {"HY": 3.34, "spring AHY": 2.43, "autumn AHY": 6.08}
rows = []
for r in RUN:
for cls in w.CLASSES:
v = prelim[(r, cls)]["vif"]
rows.append({"reading": r, "class": cls, "terms": len(v), "largest VIF": v.max(), "term": v.idxmax(),
"all below 10": bool((v < 10).all()), "paper": paper_vif[cls],
"within tolerance": abs(v.max() - paper_vif[cls]) <= round(paper_vif[cls] * 0.1, 2)})
print(pd.DataFrame(rows).round(3).to_string(index=False))
print("\nVIFs, reading (a):")
print(pd.DataFrame({cls: prelim[("a", cls)]["vif"] for cls in w.CLASSES}).round(2).to_string())
reading class terms largest VIF term all below 10 paper within tolerance
a HY 21 3.339 scaleTarsus:y_id True 3.34 True
a spring AHY 27 2.427 scaleTarsus True 2.43 True
a autumn AHY 37 6.084 scaleTarsus:y_id True 6.08 True
c HY 21 3.312 scaleTarsus:y_id True 3.34 True
c spring AHY 27 2.426 scaleTarsus True 2.43 True
c autumn AHY 37 6.088 scaleTarsus:y_id True 6.08 True
d HY 21 3.379 scaleTarsus:y_id True 3.34 True
d spring AHY 27 2.400 scaleTarsus True 2.43 True
d autumn AHY 37 6.472 scaleTarsus:y_id True 6.08 True
VIFs, reading (a):
HY spring AHY autumn AHY
AvgTempPreBreedBase 2.10 1.57 2.37
MinTempFallMig 1.33 1.18 1.68
MinTempNextSpringMig NaN 1.18 NaN
MinTempSpringMig 1.23 1.16 1.75
PrecipBreedMoBase 1.94 1.72 2.67
PrecipPreBreedBase 1.29 1.31 1.55
PrecipWinMo NaN 1.35 NaN
PrevAvgTempPreBreedBase NaN NaN 2.08
PrevMinTempFallMig NaN NaN 1.36
PrevMinTempSpringMig NaN NaN 1.73
PrevPrecipBreedBase NaN NaN 2.68
PrevPrecipPreBreedBase NaN NaN 1.38
PrevPrecipWinMo 1.20 1.20 1.75
PrevPrevPrecipWinMo NaN NaN 1.67
PrevPrevTempWinMo NaN NaN 1.79
PrevTempBreedBase NaN NaN 2.97
PrevTempWinMo 1.17 1.38 2.04
Sex21 1.11 1.13 1.12
TempBreedMoBase 2.05 1.80 3.20
TempWinMo NaN 1.48 NaN
scaleTarsus 2.76 2.43 2.82
scaleTarsus:AvgTempPreBreedBase 2.13 1.47 2.04
scaleTarsus:MinTempFallMig 1.32 1.28 1.64
scaleTarsus:MinTempNextSpringMig NaN 1.19 NaN
scaleTarsus:MinTempSpringMig 1.30 1.19 1.67
scaleTarsus:PrecipBreedMoBase 2.14 1.70 2.81
scaleTarsus:PrecipPreBreedBase 1.31 1.36 1.49
scaleTarsus:PrecipWinMo NaN 1.35 NaN
scaleTarsus:PrevAvgTempPreBreedBase NaN NaN 2.04
scaleTarsus:PrevMinTempFallMig NaN NaN 1.35
scaleTarsus:PrevMinTempSpringMig NaN NaN 1.67
scaleTarsus:PrevPrecipBreedBase NaN NaN 2.66
scaleTarsus:PrevPrecipPreBreedBase NaN NaN 1.43
scaleTarsus:PrevPrecipWinMo 1.18 1.21 1.63
scaleTarsus:PrevPrevPrecipWinMo NaN NaN 1.64
scaleTarsus:PrevPrevTempWinMo NaN NaN 1.76
scaleTarsus:PrevTempBreedBase NaN NaN 3.52
scaleTarsus:PrevTempWinMo 1.19 1.43 2.16
scaleTarsus:Sex21 2.23 2.05 2.12
scaleTarsus:TempBreedMoBase 2.24 1.82 3.07
scaleTarsus:TempWinMo NaN 1.47 NaN
scaleTarsus:y_id 3.34 2.13 6.08
y_id 2.49 2.18 4.89
Result: reproduced. Every VIF is below 10. The largest are 3.339, 2.427 and 6.084 under reading (a) (paper 3.34, 2.43 and 6.08), and 3.312 to 3.379, 2.400 to 2.427 and 6.084 to 6.472 over the three filtering readings, all within tolerance.
s4-s6-climate-variable-selection¶
The climate variables with 95% intervals excluding 0 in each preliminary model are those included in that class's main model.
In [24]:
rows = []
for r in RUN:
for cls in w.CLASSES:
s = prelim[(r, cls)]["summary"]
climate = [t for t in w.PRELIM_TERMS[cls] if t != "y_id"]
main = set(w.MAIN_PREDICTORS[cls]) - {"y_id"}
for reading in ("i", "ii"):
sig = {t for t in climate if s.loc[t, f"excl0_{reading}"]}
rows.append({"reading": r, "interval": reading, "class": cls, "selected": sorted(sig),
"missing from selection": sorted(main - sig), "extra in selection": sorted(sig - main),
"differences": len(main ^ sig)})
sel = pd.DataFrame(rows)
print(sel[["reading", "interval", "class", "missing from selection", "extra in selection", "differences"]].to_string(index=False))
print("\nDifferences summed over the three classes:")
print(sel.groupby(["reading", "interval"]).differences.sum().to_string())
print("\nBounds near 0 that the claim names (unrounded), reading (a):")
near = [("HY", "MinTempSpringMig"), ("spring AHY", "PrecipWinMo"), ("autumn AHY", "PrevPrevPrecipWinMo"), ("autumn AHY", "MinTempFallMig")]
print(pd.DataFrame([{"class": c, "term": t, **prelim[("a", c)]["summary"].loc[t, ["mean_pooled", "q2.5_pooled", "q97.5_pooled", "hpd_lo_c4", "hpd_hi_c4"]].to_dict()}
for c, t in near]).round(4).to_string(index=False))
print()
seed_count([("HY", "MinTempSpringMig"), ("spring AHY", "PrecipWinMo"), ("autumn AHY", "PrevPrevPrecipWinMo"),
("autumn AHY", "MinTempFallMig"), ("spring AHY", "MinTempNextSpringMig")],
"Selection bounds near 0 (the four the claim names, and the one extra under reading ii)")
reading interval class missing from selection extra in selection differences
a i HY [] [] 0
a ii HY [] [] 0
a i spring AHY [] [] 0
a ii spring AHY [] [MinTempNextSpringMig] 1
a i autumn AHY [] [] 0
a ii autumn AHY [] [] 0
c i HY [] [] 0
c ii HY [] [] 0
c i spring AHY [] [] 0
c ii spring AHY [] [MinTempNextSpringMig] 1
c i autumn AHY [] [] 0
c ii autumn AHY [] [] 0
d i HY [] [] 0
d ii HY [] [] 0
d i spring AHY [] [] 0
d ii spring AHY [] [] 0
d i autumn AHY [] [] 0
d ii autumn AHY [] [] 0
Differences summed over the three classes:
reading interval
a i 0
ii 1
c i 0
ii 1
d i 0
ii 0
Bounds near 0 that the claim names (unrounded), reading (a):
class term mean_pooled q2.5_pooled q97.5_pooled hpd_lo_c4 hpd_hi_c4
HY MinTempSpringMig -0.0110 -0.0238 0.0014 -0.0243 0.0007
spring AHY PrecipWinMo -0.0207 -0.0371 -0.0048 -0.0366 -0.0045
autumn AHY PrevPrevPrecipWinMo 0.0305 0.0028 0.0581 0.0015 0.0580
autumn AHY MinTempFallMig 0.0338 0.0043 0.0621 0.0095 0.0638
Selection bounds near 0 (the four the claim names, and the one extra under reading ii), 10 other seed sets:
reading class term mean range q2.5 range q97.5 range hpd_lo_c4 range hpd_hi_c4 range sets excluding 0 (i) sets excluding 0 (ii)
a HY MinTempSpringMig -0.0112 to -0.0109 -0.0240 to -0.0235 0.0016 to 0.0020 -0.0247 to -0.0224 0.0006 to 0.0028 0 0
a spring AHY PrecipWinMo -0.0209 to -0.0205 -0.0371 to -0.0363 -0.0050 to -0.0042 -0.0374 to -0.0342 -0.0054 to -0.0028 10 10
a autumn AHY PrevPrevPrecipWinMo 0.0302 to 0.0312 0.0020 to 0.0035 0.0583 to 0.0597 0.0001 to 0.0061 0.0551 to 0.0625 10 10
a autumn AHY MinTempFallMig 0.0337 to 0.0343 0.0047 to 0.0068 0.0615 to 0.0630 0.0039 to 0.0097 0.0605 to 0.0646 10 10
a spring AHY MinTempNextSpringMig -0.0134 to -0.0132 -0.0283 to -0.0274 0.0007 to 0.0015 -0.0284 to -0.0261 -0.0000 to 0.0022 0 1
c HY MinTempSpringMig -0.0105 to -0.0102 -0.0232 to -0.0227 0.0022 to 0.0026 -0.0239 to -0.0216 0.0012 to 0.0030 0 0
c spring AHY PrecipWinMo -0.0209 to -0.0206 -0.0372 to -0.0364 -0.0050 to -0.0041 -0.0375 to -0.0343 -0.0053 to -0.0028 10 10
c autumn AHY PrevPrevPrecipWinMo 0.0280 to 0.0290 -0.0000 to 0.0015 0.0559 to 0.0572 -0.0019 to 0.0041 0.0527 to 0.0600 9 7
c autumn AHY MinTempFallMig 0.0328 to 0.0334 0.0039 to 0.0061 0.0604 to 0.0619 0.0032 to 0.0090 0.0594 to 0.0634 10 10
c spring AHY MinTempNextSpringMig -0.0132 to -0.0130 -0.0282 to -0.0273 0.0009 to 0.0017 -0.0283 to -0.0259 0.0002 to 0.0024 0 0
d HY MinTempSpringMig -0.0103 to -0.0101 -0.0231 to -0.0226 0.0023 to 0.0028 -0.0237 to -0.0216 0.0013 to 0.0031 0 0
d spring AHY PrecipWinMo -0.0207 to -0.0203 -0.0369 to -0.0362 -0.0048 to -0.0040 -0.0373 to -0.0341 -0.0052 to -0.0026 10 10
d autumn AHY PrevPrevPrecipWinMo 0.0316 to 0.0325 0.0034 to 0.0048 0.0596 to 0.0610 0.0010 to 0.0073 0.0564 to 0.0639 10 10
d autumn AHY MinTempFallMig 0.0304 to 0.0311 0.0014 to 0.0035 0.0583 to 0.0600 0.0005 to 0.0062 0.0572 to 0.0613 10 10
d spring AHY MinTempNextSpringMig -0.0136 to -0.0134 -0.0285 to -0.0276 0.0006 to 0.0014 -0.0287 to -0.0265 -0.0001 to 0.0018 0 2
Result: reproduced. Under interval reading (i) the climate main effects whose interval excludes 0 are exactly the main models' climate predictors in all three classes, under every filtering reading (0 differences). Under reading (ii) MinTempNextSpringMig is also selected in the spring AHY model under filtering readings (a) and (c) (its HPD bound is the one discussed above), one difference. Over the ten other seed sets the same selection holds: three of the four bounds the claim names exclude 0 in 10 of 10 sets, and the fourth, HY MinTempSpringMig, which is not a main-model predictor, in 0 of 10, except autumn AHY PrevPrevPrecipWinMo under filtering reading (c), 9 of 10 (i) and 7 of 10 (ii), and MinTempNextSpringMig excludes 0 in at most 2 of 10. The criterion has no seed rule, so the recorded run decides.
s4-s5-climate-size-interactions¶
In the preliminary models, previous wintering temperature interacts positively with body size in HY birds, previous wintering precipitation interacts positively with body size in spring AHY birds, and no climate variable interacts with body size in autumn AHY birds.
Tarsus by climate interactions only (not sex or year).
In [25]:
named = {"HY": "scaleTarsus:PrevTempWinMo", "spring AHY": "scaleTarsus:PrevPrecipWinMo"}
rows = []
for r in RUN:
for reading in ("i", "ii"):
for cls in w.CLASSES:
s = prelim[(r, cls)]["summary"]
inter = [f"scaleTarsus:{t}" for t in w.PRELIM_TERMS[cls] if t != "y_id"]
sig = [t for t in inter if s.loc[t, f"excl0_{reading}"]]
rows.append({"reading": r, "interval": reading, "class": cls,
"named": named.get(cls, "-"),
"named positive and excludes 0": (s.loc[named[cls], f"excl0_{reading}"] and s.loc[named[cls], "mean_pooled"] > 0) if cls in named else None,
"other interactions excluding 0": [t for t in sig if t != named.get(cls)]})
inter_table = pd.DataFrame(rows)
print(inter_table.to_string(index=False))
print("\nThe named interactions and the interactions near the line (unrounded), all readings:")
show = [("HY", "scaleTarsus:PrevTempWinMo"), ("spring AHY", "scaleTarsus:PrevPrecipWinMo"), ("HY", "scaleTarsus:AvgTempPreBreedBase"),
("autumn AHY", "scaleTarsus:MinTempFallMig")]
print(pd.DataFrame([{"reading": r, "class": c, "term": t, **prelim[(r, c)]["summary"].loc[t, ["mean_pooled", "q2.5_pooled", "q97.5_pooled", "hpd_lo_c4", "hpd_hi_c4", "pMCMC_c4"]].to_dict()}
for r in RUN for c, t in show]).round(4).to_string(index=False))
reading interval class named named positive and excludes 0 other interactions excluding 0
a i HY scaleTarsus:PrevTempWinMo True []
a i spring AHY scaleTarsus:PrevPrecipWinMo True []
a i autumn AHY - None []
a ii HY scaleTarsus:PrevTempWinMo True []
a ii spring AHY scaleTarsus:PrevPrecipWinMo True []
a ii autumn AHY - None []
c i HY scaleTarsus:PrevTempWinMo True []
c i spring AHY scaleTarsus:PrevPrecipWinMo True []
c i autumn AHY - None []
c ii HY scaleTarsus:PrevTempWinMo True []
c ii spring AHY scaleTarsus:PrevPrecipWinMo True []
c ii autumn AHY - None []
d i HY scaleTarsus:PrevTempWinMo True []
d i spring AHY scaleTarsus:PrevPrecipWinMo True []
d i autumn AHY - None []
d ii HY scaleTarsus:PrevTempWinMo True []
d ii spring AHY scaleTarsus:PrevPrecipWinMo True []
d ii autumn AHY - None [scaleTarsus:PrevPrecipWinMo]
The named interactions and the interactions near the line (unrounded), all readings:
reading class term mean_pooled q2.5_pooled q97.5_pooled hpd_lo_c4 hpd_hi_c4 pMCMC_c4
a HY scaleTarsus:PrevTempWinMo 0.0155 0.0030 0.0277 0.0048 0.0284 0.010
a spring AHY scaleTarsus:PrevPrecipWinMo 0.0151 0.0003 0.0295 0.0015 0.0303 0.038
a HY scaleTarsus:AvgTempPreBreedBase -0.0159 -0.0317 0.0008 -0.0318 0.0019 0.056
a autumn AHY scaleTarsus:MinTempFallMig -0.0238 -0.0523 0.0052 -0.0502 0.0062 0.118
c HY scaleTarsus:PrevTempWinMo 0.0152 0.0025 0.0278 0.0023 0.0277 0.018
c spring AHY scaleTarsus:PrevPrecipWinMo 0.0152 0.0006 0.0296 0.0007 0.0305 0.050
c HY scaleTarsus:AvgTempPreBreedBase -0.0163 -0.0325 0.0000 -0.0327 0.0010 0.068
c autumn AHY scaleTarsus:MinTempFallMig -0.0223 -0.0517 0.0062 -0.0537 0.0041 0.122
d HY scaleTarsus:PrevTempWinMo 0.0143 0.0015 0.0270 0.0013 0.0272 0.030
d spring AHY scaleTarsus:PrevPrecipWinMo 0.0167 0.0025 0.0316 0.0023 0.0326 0.028
d HY scaleTarsus:AvgTempPreBreedBase -0.0126 -0.0287 0.0039 -0.0284 0.0037 0.134
d autumn AHY scaleTarsus:MinTempFallMig -0.0235 -0.0530 0.0062 -0.0569 0.0033 0.124
The named interactions and the interactions near the line, over the ten other seed sets (the criterion has no seed rule; the recorded run decides):
In [26]:
_ = seed_check("a", "HY", "scaleTarsus:AvgTempPreBreedBase")
seed_count([("HY", "scaleTarsus:PrevTempWinMo"), ("spring AHY", "scaleTarsus:PrevPrecipWinMo"),
("HY", "scaleTarsus:AvgTempPreBreedBase"), ("autumn AHY", "scaleTarsus:PrevPrecipWinMo")],
"The named interactions and those near the line")
scaleTarsus:AvgTempPreBreedBase, HY, reading (a), 10 other seed sets:
seed set q2.5 q97.5 hpd_lo_c4 hpd_hi_c4 excludes 0 (i) excludes 0 (ii)
0 -0.0314 0.0005 -0.0308 0.0020 False False
1 -0.0321 0.0008 -0.0333 -0.0009 False True
2 -0.0323 0.0010 -0.0328 0.0005 False False
3 -0.0323 0.0009 -0.0326 -0.0000 False True
4 -0.0324 0.0005 -0.0301 0.0021 False False
5 -0.0317 0.0002 -0.0310 0.0003 False False
6 -0.0321 0.0007 -0.0320 0.0011 False False
7 -0.0321 0.0005 -0.0320 -0.0011 False True
8 -0.0319 0.0003 -0.0329 -0.0022 False True
9 -0.0320 0.0004 -0.0323 0.0001 False False
interval (i) excludes 0 in 0 of 10; interval (ii) in 4 of 10
The named interactions and those near the line, 10 other seed sets:
reading class term mean range q2.5 range q97.5 range hpd_lo_c4 range hpd_hi_c4 range sets excluding 0 (i) sets excluding 0 (ii)
a HY scaleTarsus:PrevTempWinMo 0.0153 to 0.0156 0.0023 to 0.0031 0.0278 to 0.0283 0.0017 to 0.0041 0.0270 to 0.0290 10 10
a spring AHY scaleTarsus:PrevPrecipWinMo 0.0151 to 0.0154 0.0002 to 0.0011 0.0292 to 0.0299 0.0001 to 0.0028 0.0291 to 0.0311 10 10
a HY scaleTarsus:AvgTempPreBreedBase -0.0160 to -0.0156 -0.0324 to -0.0314 0.0002 to 0.0010 -0.0333 to -0.0301 -0.0022 to 0.0021 0 4
a autumn AHY scaleTarsus:PrevPrecipWinMo -0.0212 to -0.0205 -0.0490 to -0.0472 0.0054 to 0.0070 -0.0513 to -0.0461 0.0027 to 0.0091 0 0
c HY scaleTarsus:PrevTempWinMo 0.0148 to 0.0151 0.0020 to 0.0028 0.0272 to 0.0277 0.0013 to 0.0037 0.0265 to 0.0285 10 10
c spring AHY scaleTarsus:PrevPrecipWinMo 0.0149 to 0.0152 0.0000 to 0.0009 0.0290 to 0.0297 0.0000 to 0.0026 0.0290 to 0.0309 10 10
c HY scaleTarsus:AvgTempPreBreedBase -0.0166 to -0.0162 -0.0327 to -0.0318 -0.0005 to 0.0003 -0.0336 to -0.0305 -0.0029 to 0.0014 6 8
c autumn AHY scaleTarsus:PrevPrecipWinMo -0.0213 to -0.0206 -0.0489 to -0.0471 0.0050 to 0.0067 -0.0512 to -0.0460 0.0024 to 0.0087 0 0
d HY scaleTarsus:PrevTempWinMo 0.0143 to 0.0146 0.0012 to 0.0020 0.0268 to 0.0275 0.0007 to 0.0031 0.0264 to 0.0284 10 10
d spring AHY scaleTarsus:PrevPrecipWinMo 0.0165 to 0.0169 0.0018 to 0.0027 0.0307 to 0.0314 0.0015 to 0.0043 0.0306 to 0.0328 10 10
d HY scaleTarsus:AvgTempPreBreedBase -0.0128 to -0.0125 -0.0290 to -0.0282 0.0032 to 0.0041 -0.0304 to -0.0273 0.0009 to 0.0048 0 0
d autumn AHY scaleTarsus:PrevPrecipWinMo -0.0262 to -0.0255 -0.0539 to -0.0521 0.0005 to 0.0019 -0.0555 to -0.0505 -0.0013 to 0.0040 0 1
Result: reproduced. Under interval reading (i), and every filtering reading, the HY tarsus by PrevTempWinMo (0.0155, 0.0030 to 0.0277 under reading a) and the spring AHY tarsus by PrevPrecipWinMo (0.0151, 0.0003 to 0.0295) interactions are positive and exclude 0, and no other tarsus by climate interaction excludes 0 in any of the three models. Under reading (ii) the same holds except under filtering reading (d), where the autumn AHY tarsus by PrevPrecipWinMo HPD also excludes 0. The criterion has no seed rule, and the outcome rests on the recorded run. Over the ten other seed sets the two named interactions exclude 0 in 10 of 10 under both interval readings and every filtering reading (spring AHY lower bounds 0.0000 to 0.0043 over both interval readings and the filtering readings), but the HY tarsus by pre-breeding temperature interaction, which the recorded run puts just inside the line, excludes 0 in most sets under filtering reading (c), 6 of 10 (i) and 8 of 10 (ii), and in 4 of 10 under reading (a) (ii). Seed reading (2), the majority of seeds, would give partially reproduced under reading (c), with one other interaction excluding 0. If the seed rule of the answer of 2026-10-27 applied to this claim, the recorded run (reproduced) and the majority of seeds (partially reproduced) would disagree, and the claim would be blocked as methods_underspecified; a question asks whether it applies.
s8-climate-trends¶
Over time, pre-breeding precipitation, breeding season temperature and precipitation, autumn migration temperature and wintering temperature increase, wintering precipitation decreases, and pre-breeding and spring migration temperatures show no trend.
Each variable modelled on Time = Year - 1980 with a random effect of species with the phylogenetic
covariance, priors and chains as in the preliminary models. The code's table (envdata3) is defined
nowhere; the claims file's readings are (a) one row per species and year of the deposited climate
data, and (b) only the species-years that have specimens. For (a), the same-year columns of
ed_FHY.csv give 1981 to 2016 and the spring file's rows (its Year) add 1980; the files hold no
same-year value for 2017. The check above showed the files agree where they overlap.
In [27]:
same_year = pd.concat([ed_hy.assign(Year=ed_hy.y_id)[["spec_id", "Year", *w.TREND_VARIABLES]],
ed_sp[~ed_sp.Year.isin(ed_hy.y_id)][["spec_id", "Year", *w.TREND_VARIABLES]]]).sort_values(["spec_id", "Year"])
with_specimens = raw.groupby(["Taxon", "Year"]).size().rename("n").reset_index()
tables = {"a": same_year,
"b": same_year.merge(with_specimens, left_on=["spec_id", "Year"], right_on=["Taxon", "Year"])}
for k, t in tables.items():
print(f"reading ({k}): {len(t)} species-years, years {t.Year.min()} to {t.Year.max()}")
expected_trend = {"AvgTempPreBreedBase": (0.0, -0.01, 0.01, "none"), "PrecipPreBreedBase": (0.02, 0.01, 0.04, "+"),
"TempBreedMoBase": (0.03, 0.02, 0.04, "+"), "PrecipBreedMoBase": (0.05, 0.04, 0.06, "+"),
"MinTempFallMig": (0.02, 0.00, 0.03, "+"), "TempWinMo": (0.02, 0.01, 0.03, "+"),
"PrecipWinMo": (-0.02, -0.03, -0.01, "-"), "MinTempSpringMig": (0.0, -0.01, 0.01, "none")}
tolerance = {"PrecipPreBreedBase": 0.015, "TempBreedMoBase": 0.01, "PrecipBreedMoBase": 0.01, "MinTempFallMig": 0.015,
"TempWinMo": 0.01, "PrecipWinMo": 0.01}
rows = []
for ti, (k, t) in enumerate(tables.items()):
for vi, var in enumerate(w.TREND_VARIABLES):
d = t.dropna(subset=[var])
X = np.column_stack([np.ones(len(d)), d.Year - 1980])
chains = [w.gibbs_glmm(X, d[var].to_numpy(float), d.spec_id.to_numpy(), w.SPECIES, Ainv,
np.random.default_rng([8, ti, vi, c])) for c in range(4)]
sm = w.glmm_summary(chains, ["(Intercept)", "Time"]).loc["Time"]
pm, plo, phi, sign = expected_trend[var]
rows.append({"reading": k, "variable": var, "n": len(d), "mean": sm.mean_pooled, "q2.5": sm["q2.5_pooled"],
"q97.5": sm["q97.5_pooled"], "mean_c4": sm.mean_c4, "hpd_lo_c4": sm.hpd_lo_c4, "hpd_hi_c4": sm.hpd_hi_c4,
"pMCMC_c4": sm.pMCMC_c4, "paper": f"{pm} ({plo}, {phi})", "paper trend": sign})
trends = pd.DataFrame(rows)
def call(lo, hi):
return "+" if lo > 0 else "-" if hi < 0 else "none"
trends["trend (i)"] = [call(a, b) for a, b in zip(trends["q2.5"], trends["q97.5"])]
trends["trend (ii)"] = [call(a, b) for a, b in zip(trends.hpd_lo_c4, trends.hpd_hi_c4)]
trends["mean within tolerance"] = [abs(m - expected_trend[v][0]) <= tolerance[v] if v in tolerance else None
for m, v in zip(trends["mean"], trends.variable)]
print(trends.round(4).to_string(index=False))
for k in tables:
for reading in ("i", "ii"):
sub = trends[trends.reading == k]
agree = (sub[f"trend ({reading})"] == sub["paper trend"]).sum()
means = sub["mean within tolerance"].dropna()
print(f"data ({k}), interval ({reading}): {agree} of 8 trend calls agree with the paper; "
f"{int(means.sum())} of {len(means)} means within tolerance")
reading (a): 185 species-years, years 1980 to 2016 reading (b): 180 species-years, years 1981 to 2016
reading variable n mean q2.5 q97.5 mean_c4 hpd_lo_c4 hpd_hi_c4 pMCMC_c4 paper paper trend trend (i) trend (ii) mean within tolerance
a AvgTempPreBreedBase 183 -0.0025 -0.0160 0.0112 -0.0029 -0.0167 0.0094 0.676 0.0 (-0.01, 0.01) none none none None
a PrecipPreBreedBase 185 0.0209 0.0075 0.0341 0.0207 0.0081 0.0354 0.004 0.02 (0.01, 0.04) + + + True
a TempBreedMoBase 185 0.0339 0.0212 0.0468 0.0340 0.0206 0.0458 0.001 0.03 (0.02, 0.04) + + + True
a PrecipBreedMoBase 185 0.0546 0.0434 0.0657 0.0549 0.0443 0.0659 0.001 0.05 (0.04, 0.06) + + + True
a MinTempFallMig 184 0.0155 0.0016 0.0294 0.0151 0.0018 0.0298 0.034 0.02 (0.0, 0.03) + + + True
a TempWinMo 185 0.0240 0.0109 0.0369 0.0240 0.0108 0.0361 0.001 0.02 (0.01, 0.03) + + + True
a PrecipWinMo 185 -0.0199 -0.0329 -0.0068 -0.0202 -0.0346 -0.0081 0.002 -0.02 (-0.03, -0.01) - - - True
a MinTempSpringMig 181 -0.0045 -0.0184 0.0090 -0.0044 -0.0179 0.0092 0.498 0.0 (-0.01, 0.01) none none none None
b AvgTempPreBreedBase 178 0.0007 -0.0133 0.0151 0.0008 -0.0122 0.0165 0.916 0.0 (-0.01, 0.01) none none none None
b PrecipPreBreedBase 180 0.0147 0.0020 0.0281 0.0147 0.0016 0.0285 0.038 0.02 (0.01, 0.04) + + + True
b TempBreedMoBase 180 0.0345 0.0207 0.0483 0.0344 0.0198 0.0475 0.001 0.03 (0.02, 0.04) + + + True
b PrecipBreedMoBase 180 0.0549 0.0429 0.0667 0.0551 0.0422 0.0663 0.001 0.05 (0.04, 0.06) + + + True
b MinTempFallMig 179 0.0204 0.0067 0.0338 0.0205 0.0079 0.0341 0.004 0.02 (0.0, 0.03) + + + True
b TempWinMo 180 0.0257 0.0121 0.0393 0.0259 0.0126 0.0399 0.001 0.02 (0.01, 0.03) + + + True
b PrecipWinMo 180 -0.0225 -0.0365 -0.0089 -0.0227 -0.0352 -0.0079 0.001 -0.02 (-0.03, -0.01) - - - True
b MinTempSpringMig 176 -0.0029 -0.0172 0.0113 -0.0030 -0.0164 0.0107 0.654 0.0 (-0.01, 0.01) none none none None
data (a), interval (i): 8 of 8 trend calls agree with the paper; 6 of 6 means within tolerance
data (a), interval (ii): 8 of 8 trend calls agree with the paper; 6 of 6 means within tolerance
data (b), interval (i): 8 of 8 trend calls agree with the paper; 6 of 6 means within tolerance
data (b), interval (ii): 8 of 8 trend calls agree with the paper; 6 of 6 means within tolerance
Result: reproduced. Under both data readings and both interval readings, pre-breeding precipitation, breeding temperature, breeding precipitation, autumn migration temperature and wintering temperature increase, wintering precipitation decreases, and pre-breeding and spring migration temperatures have intervals that include 0: 8 of 8 trend calls agree with Table S8. The six trends are 0.021, 0.034, 0.055, 0.016, 0.024 and -0.020 per year under data reading (a), all within tolerance.
Summary¶
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
| Claim | Outcome |
|---|---|
| sample-sizes | reproduced |
| s4-s6-spring-migration-temperature-none | blocked, methods underspecified (not reproduced with the recorded seed, on one HPD bound at -0.0002; reproduced over ten seed sets) |
| s5-s6-ahy-prebreeding-wintering-temperature-none | reproduced |
| s6-autumn-ahy-recent-breeding-temperature-none | reproduced |
| s4-hy-autumn-migration-temperature-none | reproduced |
| s6-autumn-ahy-recent-wintering-prebreeding-precipitation-none | blocked, methods underspecified (not reproduced with the recorded seed, on one HPD bound at -0.0003 under one filtering reading; reproduced over ten seed sets) |
| s6-autumn-ahy-recent-breeding-precipitation-none | reproduced |
| s7-low-phylogenetic-signal | reproduced |
| s4-s6-no-multicollinearity | reproduced |
| s4-s6-climate-variable-selection | reproduced |
| s4-s5-climate-size-interactions | reproduced |
| s8-climate-trends | reproduced |
| table3-autumn-ahy-prior-breeding-temperature | reproduced (autumn AHY notebooks) |
| table3-autumn-ahy-prior-autumn-temperature | reproduced (autumn AHY notebooks) |
| table3-autumn-ahy-recent-autumn-temperature | partially reproduced (autumn AHY notebooks) |
| table3-autumn-ahy-prior-wintering-precipitation | partially reproduced (autumn AHY notebooks and here) |
| the 16 claims on the HY and spring AHY Stan models, the expanded ranges and the claims that depend on them | blocked, too large (stan_models.ipynb) |
How the analysis differs from the paper's:
- The preliminary models are sampled by
wingallom.gibbs_glmm, a Gibbs sampler of the same model and priors as MCMCglmm (an R package that cannot be installed in this Python environment), with seeds chosen here; the paper's runs set no seed. Two claims of absence are judged under two seed readings, the recorded seed and the majority of ten seed sets, as the answer of 2026-10-27 to this attempt's question sets out. - The Stan model is the paper's
lm_mod.stanin current array syntax, run with CmdStan 2.40 in place of rstan; with the same seeds the draws still differ. - The outlier filter of the Methods is not in the code; readings (a) to (d) of the claims file are run, (b) being identical to (a).
- The Table S8 data table is not defined in the code; it is built here from the deposited climate files under the claims file's two readings.
What was found about the dataset: the deposited specimens are already filtered, or were never
filtered (reading a gives every count of the paper exactly), and the README's statement that the
filtering steps are in the code does not hold; the 120 specimens beyond the paper's 27 366 sit in
species-years without climate data. catnodata.csv has 1050 fewer rows than morphdata.csv, 5061
rows without a catalogue number and 945 repeated catalogue numbers, and morphdata.csv has no
identifier, so measurements cannot be traced to specimens. ed_FHY.csv and ed_AHYfall.csv are
identical, and the three climate files are one table shifted by a year, consistent with the README.
The code depends on R's column order and on row numbers fixed by hand, which this notebook confirms
for the deposited files.