Spaces:
Runtime error
Runtime error
| """ | |
| make_report.py - assemble REPORT.md from the artifacts (nothing in the report is typed by hand). | |
| Inputs: artifacts/results/results.json, permutation.json (inside results.json), survival.json, | |
| artifacts/features/qc.csv, gemini_semantic.csv, gemini_retest.csv, | |
| artifacts/literature/verified_claims.json, models/config.json (optional) | |
| Outputs: REPORT.md, artifacts/results/gemini_reliability.json, artifacts/results/qc_by_batch.csv | |
| """ | |
| import json | |
| import os | |
| import sys | |
| import numpy as np | |
| import pandas as pd | |
| HERE = os.path.dirname(os.path.abspath(__file__)) | |
| sys.path.insert(0, HERE) | |
| import data as bd # noqa: E402 | |
| RES = os.path.join(HERE, "artifacts", "results") | |
| FEAT = os.path.join(HERE, "artifacts", "features") | |
| LIT = os.path.join(HERE, "artifacts", "literature", "verified_claims.json") | |
| def fmt_ci(v, ci, signed=False): | |
| f = "{:+.2f}" if signed else "{:.2f}" | |
| return f"{f.format(v)} [{f.format(ci[0])}, {f.format(ci[1])}]" | |
| def gemini_reliability(df): | |
| p1, p2 = os.path.join(FEAT, "gemini_semantic.csv"), os.path.join(FEAT, "gemini_retest.csv") | |
| if not (os.path.exists(p1) and os.path.exists(p2)): | |
| return None | |
| a = pd.read_csv(p1).set_index("file"); b = pd.read_csv(p2).set_index("file") | |
| common = a.index.intersection(b.index) | |
| a, b = a.loc[common], b.loc[common] | |
| out = dict(n=int(len(common))) | |
| for col, key in (("gem_is_bee", "agree_is_bee"), ("gem_activity", "agree_activity"), ("gem_quality", "agree_quality")): | |
| m = a[col].notna() & b[col].notna() | |
| out[key] = float((a.loc[m, col] == b.loc[m, col]).mean()) if m.any() else None | |
| # chance-corrected agreement (Cohen's kappa) | |
| if m.any(): | |
| x, y = a.loc[m, col].astype(int), b.loc[m, col].astype(int) | |
| po = float((x == y).mean()) | |
| cats = sorted(set(x) | set(y)) | |
| pe = float(sum((x == c).mean() * (y == c).mean() for c in cats)) | |
| out[key.replace("agree", "kappa")] = float((po - pe) / (1 - pe)) if pe < 1 else None | |
| out["is_bee_rate_pass1"] = float(a["gem_is_bee"].mean()) | |
| out["is_bee_rate_pass2"] = float(b["gem_is_bee"].mean()) | |
| out["mean_abs_conf_diff"] = float((a["gem_confidence"] - b["gem_confidence"]).abs().mean()) | |
| # by site | |
| meta = df.set_index(df["sound_file"].str.lower()) | |
| out["is_bee_rate_by_country"] = {c: float(a.loc[[f for f in common if meta.loc[f, "country"] == c], "gem_is_bee"].mean()) | |
| for c in sorted(set(meta["country"]))} | |
| if "gem_notes" in a.columns: | |
| out["examples_disagreeing"] = [dict(file=f, pass1=str(a.loc[f, "gem_dominant_sounds"]), pass2=None) | |
| for f in list(common[(a["gem_is_bee"] != b["gem_is_bee"]).values])[:5]] | |
| return out | |
| def qc_by_batch(df): | |
| p = os.path.join(FEAT, "qc.csv") | |
| if not os.path.exists(p): | |
| return None | |
| q = pd.read_csv(p); q["file"] = q["file"].str.lower() | |
| m = df.assign(file=df["sound_file"].str.lower()).merge(q, on="file") | |
| g = m.groupby("recording_batch")[["duration_s", "rms_db", "noise_floor_db", "level_range_db", "clip_frac"]].median().round(2) | |
| g["n"] = m.groupby("recording_batch").size() | |
| g["n_terminal_46"] = m.groupby("recording_batch")["y46"].apply(lambda s: int((s == "T").sum())) | |
| g.to_csv(os.path.join(RES, "qc_by_batch.csv")) | |
| return g | |
| def main(): | |
| df = bd.load_workbook(os.path.join(HERE, "AI_Data_Training.xlsx"), audio_dir=None) | |
| R = json.load(open(os.path.join(RES, "results.json"))) | |
| meta, res = R["meta"], R["results"] | |
| S = json.load(open(os.path.join(RES, "survival.json"))) if os.path.exists(os.path.join(RES, "survival.json")) else None | |
| lit = json.load(open(LIT)) if os.path.exists(LIT) else [] | |
| rel = gemini_reliability(df) | |
| if rel: | |
| json.dump(rel, open(os.path.join(RES, "gemini_reliability.json"), "w"), indent=1) | |
| qb = qc_by_batch(df) | |
| L = [] | |
| L.append("# What hive sound can and cannot predict: colony survival, viral loads and the missing Varroa data") | |
| L.append("") | |
| L.append("Generated by make_report.py from the pipeline artifacts. Every number below is computed; none is typed in.") | |
| L.append("") | |
| # ------------------------------------------------------------------ data | |
| n_rows = len(df); n_col = df["colony"].nunique() | |
| L.append("## 1. Data") | |
| L.append("") | |
| L.append(f"- Workbook rows: {n_rows}; recordings with an audio file: {meta['n_recordings']}; colonies: {n_col} " | |
| f"(USA {df[df.country=='USA'].colony.nunique()}, NZ {df[df.country=='NZ'].colony.nunique()}).") | |
| L.append(f"- Measured quantities per recording month: {', '.join(c for c in df.columns if c in ('CBPV','DWV','KBV'))} (viral loads) and the two survival labels. " | |
| "There is no Varroa mite count, colony size, queen status, brood, weight or nutrition column in the workbook, so none of those can be modelled or tested here.") | |
| pos = df.groupby("country")[["CBPV", "DWV", "KBV"]].agg(lambda s: (s > 0).mean()) | |
| L.append("- Fraction of recordings with a detectable load, by site: " + "; ".join( | |
| f"{c}: CBPV {pos.loc[c,'CBPV']:.2f}, DWV {pos.loc[c,'DWV']:.2f}, KBV {pos.loc[c,'KBV']:.2f}" for c in pos.index) + ".") | |
| t46 = df[df.y46 == "T"]; t03 = df[df.y03 == "T"] | |
| L.append(f"- Terminal labels: {len(t03)} recordings (0-3 months) and {len(t46)} recordings (4-6 months) from " | |
| f"{df[(df.y03=='T')|(df.y46=='T')].colony.nunique()} colonies; 12 of those colonies (F1-F12) were recorded once, in October 2023, and are all Terminal.") | |
| L.append("- Label structure: within a colony the labels move monotonically S/S -> S/T -> T/T towards the month of death, i.e. the labels encode months-to-death.") | |
| L.append("") | |
| # ------------------------------------------------------------------ method | |
| L.append("## 2. Method") | |
| L.append("") | |
| L.append(f"- Features per recording: handcrafted acoustic features ({meta['families'].get('hc','-')}), AST embedding ({meta['families'].get('ast','-')}: mean and SD of 768-d window embeddings), " | |
| f"CLAP embedding ({meta['families'].get('clap','-')}), Gemini structured listening ({meta['families'].get('gem','-')}), recording QC ({meta['families'].get('qc','-')}), metadata ({', '.join(meta.get('metadata_columns') or [])}).") | |
| L.append(f"- Validation: {meta['n_splits']}-fold colony-grouped cross-validation repeated {meta['n_repeats']} times (all recordings of a colony on one side), " | |
| f"metrics on out-of-fold predictions; colony-cluster bootstrap ({meta['n_boot']} resamples) 95% CI; month-aligned colony permutation test ({meta['n_perm']} permutations: " | |
| "colony identities permuted within country and each recording given the donor colony's label for the same calendar month, so site, season, batch and within-colony label structure are preserved and p-values test for acoustic signal beyond them).") | |
| L.append(f"- Models fixed in advance: PCA({meta['pca_k']}) + L2 logistic / ridge for embeddings, gradient boosting for handcrafted features, rank-average ensemble as the primary model. No per-target model selection.") | |
| L.append(f"- Confound controls: metadata-only baseline and paired bootstrap of the gain of audio over it; within-country evaluation with its own metadata baseline; cross-site transfer; every target re-run without the October 2023 batch; probes that ask whether audio can identify site, season and the October 2023 batch. Evaluations need at least {meta.get('min_pos_colonies', 8)} positive and {meta.get('min_pos_colonies', 8)} negative colonies.") | |
| L.append("- Verdict rule (fixed before the run): level 0 = the primary CI includes chance (the upper bound is stated); level 1 = CI excludes chance but the permutation p is at least 0.05 or the paired gain over metadata is not positive; level 2 = CI excludes chance, permutation p below 0.05 and paired gain over metadata positive. Holm-adjusted p-values across the non-probe targets are shown in brackets.") | |
| L.append("") | |
| # ------------------------------------------------------------------ results table | |
| L.append("## 3. Which targets are predictable from sound") | |
| L.append("") | |
| L.append("Metric: ROC-AUC (chance 0.50) for binary targets, Spearman correlation (chance 0) for levels; primary model = ensemble. " | |
| "'Gain' is the paired-bootstrap difference between metadata+audio and metadata alone. Verdicts follow the pre-registered rule in evaluate.py.") | |
| L.append("") | |
| L.append("| target | n (pos) [positive colonies] | metadata only | ensemble | gain over metadata | within-country gain of audio over metadata | perm. p (Holm) | without Oct 2023 batch (ensemble) | verdict |") | |
| L.append("|---|---|---|---|---|---|---|---|---|") | |
| order = sorted(res, key=lambda k: (res[k]["verdict"]["level"] * -1, k)) | |
| for name in order: | |
| r = res[name] | |
| if "families" not in r: | |
| L.append(f"| {name} | {r.get('n','')} | | | | | | {r['verdict']['label']} ({r['verdict'].get('reason','')}) |") | |
| continue | |
| key = "auc" if r["kind"] == "binary" else "spearman" | |
| f = r["families"] | |
| def c(k): | |
| return fmt_ci(f[k][key], f[k]["ci"], signed=(key == "spearman")) if k in f else "-" | |
| d = r.get("delta_over_meta", {}).get("meta+all_audio") | |
| gain = f"{d['mean']:+.2f} [{d['lo']:+.2f}, {d['hi']:+.2f}]" if d else "-" | |
| wcs = [] | |
| for cn, wc in (r.get("within_country") or {}).items(): | |
| if "skipped" in wc: | |
| wcs.append(f"{cn}: {wc['skipped'].replace('not estimable: ', '')}") | |
| else: | |
| dd = (wc.get("delta_over_meta") or {}).get("meta+all_audio") | |
| au = wc["families"].get("all_audio") | |
| wcs.append(f"{cn}: audio {fmt_ci(au[key], au['ci'], signed=(key=='spearman'))}" + (f", gain {dd['mean']:+.2f} [{dd['lo']:+.2f}, {dd['hi']:+.2f}]" if dd else "")) | |
| within = "; ".join(wcs) or "-" | |
| v = r["verdict"] | |
| p = v.get("p_perm"); ph = v.get("p_perm_holm") | |
| ps = "-" if p is None else f"{p:.3f}" + (f" ({ph:.3f})" if ph is not None else "") | |
| wo = v.get("without_oct2023") | |
| wos = fmt_ci(wo["value"], wo["ci"], signed=(key == "spearman")) if wo else "-" | |
| npos = f"{r['n']} ({r.get('n_pos','')}) [{r.get('n_pos_colonies','')}]" if r["kind"] == "binary" else str(r["n"]) | |
| lab = v["label"] + (f"; {v['upper_bound']}" if v.get("upper_bound") else "") | |
| L.append(f"| {name} | {npos} | {c('meta')} | {c('ensemble')} | {gain} | {within} | {ps} | {wos} | {lab} |") | |
| L.append("") | |
| # per-family detail | |
| L.append("### 3a. Per-family detail (point estimate [95% CI])") | |
| L.append("") | |
| fams_all = [k for k in ("meta", "qc", "gem", "hc", "ast", "clap", "ast+hc", "all_audio", "meta+all_audio", "lab+meta", "lab+meta+all_audio", "ensemble")] | |
| present = [k for k in fams_all if any(k in r.get("families", {}) for r in res.values())] | |
| L.append("| target | " + " | ".join(present) + " |") | |
| L.append("|---|" + "---|" * len(present)) | |
| for name in order: | |
| r = res[name] | |
| if "families" not in r: | |
| continue | |
| key = "auc" if r["kind"] == "binary" else "spearman" | |
| cells = [fmt_ci(r["families"][k][key], r["families"][k]["ci"], signed=(key == "spearman")) if k in r["families"] else "-" for k in present] | |
| L.append(f"| {name} | " + " | ".join(cells) + " |") | |
| L.append("") | |
| # transfer + no-oct | |
| L.append("### 3b. Cross-site transfer and the October 2023 batch") | |
| L.append("") | |
| L.append("| target | family | USA->NZ | NZ->USA | without Oct 2023 batch (ensemble members) |") | |
| L.append("|---|---|---|---|---|") | |
| for name in order: | |
| r = res[name] | |
| tr = r.get("transfer", {}) | |
| key = "auc" if r["kind"] == "binary" else "spearman" | |
| for famk in ("all_audio",): | |
| t = tr.get(famk, {}) | |
| def tc(k): | |
| v = t.get(k); return fmt_ci(v[key], v["ci"], signed=(key == "spearman")) if v else "-" | |
| wo = r.get("without_oct2023", {}).get("families", {}) | |
| wos = "; ".join(f"{k}: {fmt_ci(v[key], v['ci'], signed=(key=='spearman'))}" for k, v in wo.items() if k in ("meta", "all_audio", "ensemble")) or "-" | |
| if t or wo: | |
| L.append(f"| {name} | {famk} | {tc('USA->NZ')} | {tc('NZ->USA')} | {wos} |") | |
| L.append("") | |
| # ------------------------------------------------------------------ probes | |
| L.append("## 4. What the audio actually encodes: site, season and recording conditions") | |
| L.append("") | |
| for name in ("probe_country_is_nz", "probe_batch_oct2023", "probe_month_sin", "probe_month_cos"): | |
| r = res.get(name) | |
| if r and "families" in r: | |
| key = "auc" if r["kind"] == "binary" else "spearman" | |
| best = max((k for k in ("hc", "ast", "clap", "all_audio") if k in r["families"]), key=lambda k: r["families"][k][key]) | |
| L.append(f"- {r['description']}: best audio family {best} {fmt_ci(r['families'][best][key], r['families'][best]['ci'], signed=(key=='spearman'))}.") | |
| if qb is not None: | |
| L.append("") | |
| L.append("Recording-quality statistics by batch (median per site-year-month; full table in artifacts/results/qc_by_batch.csv):") | |
| L.append("") | |
| L.append("| batch | n | Terminal(4-6) rows | duration s | RMS dB | noise floor dB | level range dB |") | |
| L.append("|---|---|---|---|---|---|---|") | |
| for b, row in qb.iterrows(): | |
| L.append(f"| {b} | {int(row['n'])} | {int(row['n_terminal_46'])} | {row['duration_s']:.1f} | {row['rms_db']:.1f} | {row['noise_floor_db']:.1f} | {row['level_range_db']:.1f} |") | |
| L.append("") | |
| # ------------------------------------------------------------------ gemini | |
| L.append("## 5. Gemini") | |
| L.append("") | |
| gm = [k for k in res if "gem" in res[k].get("families", {})] | |
| if rel: | |
| L.append(f"- Structured listening (model {os.environ.get('GEMINI_MODELS', 'gemini-3.8-flash').split(',')[0]}, 30 s excerpt, temperature 0) was run twice on {rel['n']} recordings. " | |
| f"Agreement between the two passes: 'is honeybee colony audio' {100*rel['agree_is_bee']:.0f}% (kappa {rel.get('kappa_is_bee', float('nan')):.2f}), " | |
| f"activity level {100*rel['agree_activity']:.0f}% (kappa {rel.get('kappa_activity', float('nan')):.2f}), recording quality {100*rel['agree_quality']:.0f}% (kappa {rel.get('kappa_quality', float('nan')):.2f}). " | |
| f"Share of recordings called honeybee audio: pass 1 {100*rel['is_bee_rate_pass1']:.0f}%, pass 2 {100*rel['is_bee_rate_pass2']:.0f}%; by site: " | |
| + ", ".join(f"{c} {100*v:.0f}%" for c, v in rel["is_bee_rate_by_country"].items()) + ".") | |
| if gm: | |
| L.append("- Predictive value of the Gemini listening features (family 'gem', linear model, colony-grouped CV):") | |
| for name in order: | |
| r = res[name] | |
| if "gem" in r.get("families", {}): | |
| key = "auc" if r["kind"] == "binary" else "spearman" | |
| L.append(f" - {name}: {fmt_ci(r['families']['gem'][key], r['families']['gem']['ci'], signed=(key=='spearman'))}") | |
| L.append("- The app also lets Gemini write a summary, but only from a JSON evidence object of computed numbers; the numbers are displayed next to the text so every statement can be checked.") | |
| L.append("") | |
| # ------------------------------------------------------------------ survival | |
| if S: | |
| L.append("## 6. Factors associated with colony death (observational)") | |
| L.append("") | |
| L.append("Full tables: artifacts/results/survival_report.md. Key estimates (odds ratio per unit log(1+load) or per presence, 95% CI):") | |
| L.append("") | |
| for g in S["gee"]: | |
| if g["outcome"] == "y46" and g.get("terms") and all(np.isfinite(v["odds_ratio"]) for v in g["terms"].values()): | |
| L.append(f"- GEE, {g['label']} (n={g['n']}, events={g['n_events']} rows): " + "; ".join( | |
| f"{k} OR {v['odds_ratio']:.2f} [{v['ci'][0]:.2f}, {v['ci'][1]:.2f}]" for k, v in g["terms"].items() if not k.startswith("month")) + ".") | |
| for g in S.get("firth", []): | |
| if g["outcome"] == "y46" and g["label"].startswith("PRIMARY"): | |
| L.append(f"- Firth (Wald CI), {g['label']} (n={g['n']}, event colonies={g['n_event_colonies']}, separation={g['separation_detected']}): " + "; ".join( | |
| (f"{k} OR {v['odds_ratio']:.2f} [{v['ci'][0]:.2f}, {v['ci'][1]:.2f}]" if v.get("estimable", True) else f"{k} not estimable") | |
| for k, v in g["terms"].items() if not k.startswith("month")) + ".") | |
| ptc = S.get("positive_terminal_counts", {}) | |
| if ptc: | |
| L.append("- Identifiability (Terminal colonies positive / Terminal colonies, outcome y46): " + "; ".join( | |
| f"{k.split('|')[0]} {k.split('|')[2]} {v['terminal_colonies_positive']}/{v['terminal_colonies']}" for k, v in ptc.items() if k.split("|")[1] == "y46") + ".") | |
| ct = S["colony_tests"] | |
| L.append("- Colony level (died vs survived): " + "; ".join( | |
| f"{c} ever positive {t['ever_positive_died']} vs {t['ever_positive_survived']} (Fisher p {t['fisher_p']:.3f})" for c, t in ct.items() if c not in ("n_recordings", "well_monitored")) + ".") | |
| avl = S.get("audio_vs_lab") | |
| if avl: | |
| for t, r in avl.items(): | |
| d = r.get("delta_audio_over_lab_meta") | |
| if d and r.get("lab+meta") and r.get("lab+meta+all_audio"): | |
| L.append(f"- {t}: laboratory+metadata AUC {fmt_ci(r['lab+meta']['auc'], r['lab+meta']['ci'])}; adding audio {fmt_ci(r['lab+meta+all_audio']['auc'], r['lab+meta+all_audio']['ci'])}; gain {d['mean']:+.3f} [{d['lo']:+.3f}, {d['hi']:+.3f}].") | |
| L.append("") | |
| L.append("Interpretation: the twelve October 2023 colonies were sampled because they were failing (selection on the outcome), and they supply " | |
| "every DWV-positive Terminal row in the USA. The primary specification therefore excludes them; read the all-recordings odds ratios as a sensitivity run only. " | |
| "Varroa, the established vector of DWV, is unmeasured, so the causal pathway Varroa -> DWV -> death cannot be tested with this workbook.") | |
| L.append("") | |
| # ------------------------------------------------------------------ recommendations | |
| L.append("## 7. What to collect next") | |
| L.append("") | |
| L.append("1. Varroa mite counts (alcohol wash or sugar roll per 300 bees) at every recording visit. Varroa is the strongest documented predictor of colony loss and the vector of DWV; it is absent from the workbook, and the only acoustic-Varroa study with real scale (Yildiz and Karabag 2026) does not state whether its split was hive-independent.") | |
| L.append("2. Colony strength at each visit (frames of bees, brood area, adult population estimate) and queen status. These are the colony states with peer-reviewed evidence of acoustic detectability (queen loss, swarming) and the most plausible mediators of any sound-survival link.") | |
| L.append("3. A recording protocol that fixes the confounds seen here: same recorder and gain at both sites, microphone position inside the hive, time of day and ambient temperature logged, and a short silence sample per visit. The audio currently identifies site and season far better than any pathogen, so recording conditions must be equalised before acoustic pathogen claims are possible.") | |
| L.append("4. Pathogen sampling with matched timing and units: the workbook mixes zero-inflated loads whose units are not documented; record the assay, the unit and the detection limit, and sample the same bees that were recorded.") | |
| L.append("5. Longer follow-up of the same colonies (outcome at 12 months) and more dying colonies with monthly series; the current 12 monitored deaths are too few for a reliable survival model of any kind.") | |
| L.append("") | |
| # ------------------------------------------------------------------ literature | |
| if lit: | |
| L.append("## 8. Verified references") | |
| L.append("") | |
| L.append("Each entry was located and read by an agent and then independently re-checked against the source by a second agent; only entries that survived that check are listed. Claims are stated as the checker rewrote them.") | |
| L.append("") | |
| for topic in lit: | |
| L.append(f"### {topic['topic']}") | |
| L.append("") | |
| for cl in topic["verified"]: | |
| v = cl["verdict"] | |
| if v.get("supported"): | |
| L.append(f"- {v.get('corrected_claim') or cl['claim']} ") | |
| L.append(f" {v.get('corrected_citation') or cl['citation']}. {cl['url']}") | |
| notes = topic.get("notes", "").strip() | |
| if notes: | |
| L.append("") | |
| L.append(f"Search notes (what could not be found or verified): {notes}") | |
| L.append("") | |
| # ------------------------------------------------------------------ limitations | |
| L.append("## 9. Limitations") | |
| L.append("") | |
| L.append("- Observational data; two sites with different virus profiles, recorders and seasons; 12 of 24 deaths come from one single-visit batch.") | |
| L.append("- Viral load units and assay detection limits are not documented in the workbook.") | |
| L.append("- Gemini listening is a large language model's impression of a 30 s excerpt; its repeatability is reported above and it must not be read as a measurement.") | |
| L.append("- The deployed models are refit on all recordings; their expected performance on new colonies is the colony-grouped estimate in section 3, not their fit on the training recordings.") | |
| md = "\n".join(L) | |
| open(os.path.join(HERE, "REPORT.md"), "w").write(md) | |
| print(md[:3000]) | |
| print("...\nwrote REPORT.md") | |
| if __name__ == "__main__": | |
| main() | |