#!/usr/bin/env python3 # -*- coding: utf-8 -*- """longevity_models.py — Component 5 · competing-model comparison (AICc workbench) Part of the Seed Longevity Lab's downloadable Python lab (site v0.17.0, longevity.seedsofhawaii.org/scripts). Embeds Components 1 and 3's grouping; independent of the RCI component entirely. Standalone by design: the shared cores below are embedded verbatim in every script of the suite, so this file has NO imports from its siblings and no third-party dependencies. Python 3.9+. Input: a raw completed-test CSV (one row per completed germination test) OR the suite's normalized_assays.csv interchange file. Scientific contract: the Longevity Lab developer handoff (specification/02_ANALYSIS_CONTRACT). Numerically cross-validated against the site's JS engine, against R-equivalent references, and against the matching R script (this component ships in both languages) by the site's test suite. Each component's page card at /scripts spells out exactly where behavior differs from the archival workflow behind the published catalog. Counts rule, everywhere: actual number sown/tested + number germinated are required. Nothing in this suite ever assumes a sample size. """ import sys # --------------------------------------------------------------------------- # SHARED CORE · statistics (identical block in every Longevity Lab script) # Pure standard library. The IRLS binomial-probit GLM here is the same # algorithm as the site engine (src/lib/formal.js) and is cross-validated # against it, against R-equivalent statsmodels references, and against the # viabilitymetrics R package lineage. No sample size is ever assumed. # --------------------------------------------------------------------------- import math def normcdf(z): return 0.5 * (1.0 + math.erf(z / math.sqrt(2.0))) def dnorm(z): return math.exp(-0.5 * z * z) / math.sqrt(2.0 * math.pi) def probit(p): """Inverse normal CDF (Acklam) polished with one Newton step.""" if p <= 0.0 or p >= 1.0: raise ValueError("probit needs 0 0: dev += g * math.log(g / (n * m)) if n - g > 0: dev += (n - g) * math.log((n - g) / (n * (1 - m))) return 2 * dev def glm_probit(rows): """Binomial GLM, probit link, general design matrix (IRLS). rows: list of (x_vector, n, g). Returns dict or None. Separation/extreme fits are FLAGGED, never bias-reduced (the site's 'flagged, not faked' rule; see docs/DIFFERENCES_FROM_R_REFERENCE.md).""" usable = [(list(x), float(n), float(g)) for x, n, g in rows if n > 0 and 0 <= g <= n and all(math.isfinite(v) for v in x)] if not usable: return None p = len(usable[0][0]) totn = sum(r[1] for r in usable) totg = sum(r[2] for r in usable) beta = [0.0] * p beta[0] = probit(min(0.999, max(0.001, totg / max(1.0, totn)))) dev = float("inf") converged = False for _ in range(200): eta = [sum(xi * beta[j] for j, xi in enumerate(x)) for x, n, g in usable] mu = [min(1 - 1e-12, max(1e-12, normcdf(e))) for e in eta] XtWX = [[0.0] * p for _ in range(p)] XtWz = [0.0] * p for i, (x, n, g) in enumerate(usable): phi = dnorm(eta[i]) v = mu[i] * (1 - mu[i]) w = n * phi * phi / max(v, 1e-12) z = eta[i] + (g / n - mu[i]) / max(phi, 1e-12) for a2 in range(p): XtWz[a2] += w * x[a2] * z for b2 in range(p): XtWX[a2][b2] += w * x[a2] * x[b2] nxt = _solve([row[:] for row in XtWX], XtWz) if nxt is None or any(not math.isfinite(v) for v in nxt): break cand = nxt for _h in range(21): ceta = [sum(xi * cand[j] for j, xi in enumerate(x)) for x, n, g in usable] cmu = [min(1 - 1e-12, max(1e-12, normcdf(e))) for e in ceta] cdev = _binom_deviance(usable, cmu) if not (math.isfinite(dev) and cdev > dev + 1e-8): break cand = [(v + beta[j]) / 2 for j, v in enumerate(cand)] delta = max(abs(v - beta[j]) for j, v in enumerate(cand)) beta = cand rel = abs(cdev - dev) / (abs(cdev) + 0.1) dev = cdev if delta < 1e-10 or rel < 1e-13: converged = True break eta = [sum(xi * beta[j] for j, xi in enumerate(x)) for x, n, g in usable] mu = [min(1 - 1e-12, max(1e-12, normcdf(e))) for e in eta] extreme = any(abs(e) > 6.3 for e in eta) pinned = all(m > 0.999 or m < 0.001 for m in mu) separated = (not converged) or pinned or any(abs(b3) > 40 for b3 in beta) XtWX = [[0.0] * p for _ in range(p)] for i, (x, n, g) in enumerate(usable): phi = dnorm(eta[i]) v = mu[i] * (1 - mu[i]) w = n * phi * phi / max(v, 1e-12) for a2 in range(p): for b2 in range(p): XtWX[a2][b2] += w * x[a2] * x[b2] cov = _invert(XtWX) loglik = sum(g * math.log(mu[i]) + (n - g) * math.log(1 - mu[i]) for i, (x, n, g) in enumerate(usable)) return {"coef": beta, "cov": cov, "deviance": dev, "df_residual": len(usable) - p, "converged": converged, "extreme_prob": extreme, "separated": separated, "fitted": mu, "loglik": loglik} def crossing(intercept, slope, target_prop): if not (math.isfinite(intercept) and math.isfinite(slope)) or slope >= 0: return None if not (0 < target_prop < 1): return None return (probit(target_prop) - intercept) / slope def crossing_display(raw): if raw is None or not math.isfinite(raw): return "not interpretable" if raw <= 0: return "not interpretable (crossing ≤0 y)" return f"{raw:.2f} y" def chisq1_sf(x): if not math.isfinite(x) or x < 0: return float("nan") return 2 * (1 - normcdf(math.sqrt(x))) # --------------------------------------------------------------------------- # SHARED CORE · upload/normalization layer (identical block in every script) # Port of the site's src/lib/upload.js — the handoff contract: # · one row = one completed germination test (same age may repeat) # · actual tested + germinated counts required; percent recomputed from # counts; NOTHING here assumes a sample size for percentage-only data # · time in days / months / years / dates; original value + unit retained # alongside exact normalized days and years; no rounding, ever # · Initial is a baseline state (is_initial preferred; legacy # storage="Initial" recognized), not a storage treatment # · protocol_key = substrate + germination treatment, trivial case and # whitespace normalization ONLY — never synonym inference # · drying treatment optional; internal row IDs; ambiguous same-date # duplicates flagged for review, never merged or deleted # --------------------------------------------------------------------------- import csv, io, re, sys from datetime import date, timedelta DAYS_PER_YEAR = 365.25 DAYS_PER_MONTH = 365.25 / 12.0 ROLES = { "family": ["family"], "taxon": ["taxonname", "taxon", "species", "scientificname", "speciesname"], "accession": ["acquisitionnum", "accession", "accessionnum", "accessionnumber", "accno", "lot", "lotid", "collectionid"], "assay_id": ["assayid", "testid", "germinationtestid", "assay", "testnumber", "testno"], "source_institution": ["sourceinstitution", "institution", "seedbank", "source"], "is_initial": ["isinitial", "initial", "initialtest"], "storage": ["storedrycode", "storagecondition", "storagecode", "storecode", "storage", "storagetemp", "storagetemperature", "storagetempc"], "drying": ["dryingtreatment", "drying", "drytreatment", "dryingprotocol"], "substrate": ["substrate", "medium"], "treatment": ["germinationtreatment", "treatment", "pretreatment"], "tested": ["numbersown", "numsown", "sown", "seedsown", "seedssown", "tested", "numbertested", "ntested", "nsown", "n"], "germinated": ["numbergermin", "numbergerminated", "germinated", "numgerminated", "germin", "germ", "ngerminated"], "percent": ["%germin", "percentgermination", "percentgerminated", "germinationpercent", "germinationpct", "pctgermin", "percent", "pct", "viability", "viabilitypct"], "age_days": ["storageagedays", "storagedays", "days", "day", "ageindays"], "age_months": ["storageagemonths", "storagemonths", "months", "month", "ageinmonths"], "age_years": ["yearsdifference", "storageageyears", "storageyears", "years", "year", "yr", "yeardiff", "yearsdiff", "ageinyears"], "age_ambiguous": ["storageage", "age", "time", "interval", "storagetime"], "collection_date": ["collectiondate", "storagestartdate", "collecteddate", "datecollected", "collectdate"], "test_date": ["teststartdate", "germinationteststartdate", "testdate", "germinationdate", "startdate", "sowdate", "datesown"], "notes": ["notes", "comments", "note"], } TRUTHY = {"true", "t", "yes", "y", "1", "initial"} FALSY = {"false", "f", "no", "n", "0", ""} def _norm_header(h): return re.sub(r"[^a-z0-9%]", "", str(h or "").lower()) def detect_roles(headers): mapping = {} for i, h in enumerate(headers): n = _norm_header(h) for role, names in ROLES.items(): if role not in mapping and n in names: mapping[role] = i break return mapping def norm_protocol_text(x): s = re.sub(r"\s+", " ", str(x or "").strip()).lower() return s if s else "(missing)" def protocol_key(substrate, treatment): s, t = norm_protocol_text(substrate), norm_protocol_text(treatment) if s == "(missing)" and t == "(missing)": return "unspecified" return f"{s} || {t}" def normalize_storage(raw): s = str(raw or "").strip() l = re.sub(r"[°º]", "", re.sub(r"\s+", "", s.lower())) l = re.sub(r"[−–—]", "-", l) if re.match(r"^init", s, re.I): return {"code": "Initial", "original": s, "legacy_initial": True} if re.match(r"^-18c?$", l) or l in ("minus18", "minus18c"): return {"code": "-18C", "original": s, "legacy_initial": False} if re.match(r"^\+?5c?$", l): return {"code": "5C", "original": s, "legacy_initial": False} if l.startswith("amb"): return {"code": "Ambient", "original": s, "legacy_initial": False} return {"code": s, "original": s, "legacy_initial": False} def _parse_date_loose(v): s = str(v or "").strip() if not s: return None if re.match(r"^\d+(\.\d+)?$", s): num = float(s) # Excel serial (days since 1899-12-30) if 59 < num < 200000: return date(1899, 12, 30) + timedelta(days=num) return None m = re.match(r"^(\d{4})-(\d{1,2})-(\d{1,2})", s) if m: try: return date(int(m.group(1)), int(m.group(2)), int(m.group(3))) except ValueError: return None m = re.match(r"^(\d{1,2})[/\-](\d{1,2})[/\-](\d{2,4})$", s) if m: yr = int(m.group(3)) yr = 2000 + yr if yr < 100 else yr try: return date(yr, int(m.group(1)), int(m.group(2))) except ValueError: return None return None def _f(v): try: x = float(str(v).replace("%", "").strip()) return x except (ValueError, TypeError): return None NORMALIZED_SIGNATURE = {"row_id", "is_initial", "protocol_key", "storage_age_years", "tested", "germinated"} def parse_completed_tests(text, time_unit=None): """Raw completed-test CSV/TSV → canonical assay dicts. Returns dict with assays, issues, mapping, headers, needs_time_unit, percent_only, duplicate_groups, rejected.""" issues, assays, rejected = [], [], [] lines = [l for l in text.lstrip("").strip().splitlines() if l.strip()] out = {"assays": assays, "issues": issues, "mapping": {}, "headers": [], "needs_time_unit": False, "percent_only": False, "duplicate_groups": [], "rejected": rejected} if not lines: issues.append(("error", "No data.")) return out delim = "\t" if "\t" in lines[0] else (";" if lines[0].count(";") > lines[0].count(",") else ",") rows = list(csv.reader(io.StringIO("\n".join(lines)), delimiter=delim)) headers = [h.strip() for h in rows[0]] out["headers"] = headers mapping = detect_roles(headers) out["mapping"] = mapping for role in ("taxon", "accession"): if role not in mapping: issues.append(("error", f"Missing a {role} column. Headers seen: {', '.join(headers)}")) if "storage" not in mapping and "is_initial" not in mapping: issues.append(("error", "Missing a storage-condition column (post-storage rows need one; Initial rows may use is_initial).")) has_counts = "tested" in mapping and "germinated" in mapping has_pct = "percent" in mapping if not has_counts and not has_pct: issues.append(("error", 'No germination data found — need "tested" (number sown) + "germinated" counts.')) out["percent_only"] = (not has_counts) and has_pct if out["percent_only"]: issues.append(("block", "Counts required for formal modelling: this file has germination percentages but no number " "sown/tested and number germinated. A binomial longevity model requires those counts because " "80% based on 8/10 seeds carries different statistical information from 80/100. Add the " "original counts and upload again — the site will not assume a sample size.")) has_days = "age_days" in mapping has_months = "age_months" in mapping has_years = "age_years" in mapping has_amb = "age_ambiguous" in mapping has_dates = "collection_date" in mapping and "test_date" in mapping amb_unit = None if not (has_days or has_months or has_years) and has_amb: if time_unit in ("days", "months", "years"): amb_unit = time_unit issues.append(("info", f'Ambiguous time column read as {time_unit} (your --time-unit selection).')) else: out["needs_time_unit"] = True issues.append(("choice", f'The time column "{headers[mapping["age_ambiguous"]]}" doesn\'t say its unit. ' "Pass --time-unit days|months|years — this script won't guess.")) if not (has_days or has_months or has_years or has_amb or has_dates): issues.append(("error", "No storage-age column found — need days / months / years, or CollectionDate + TestStartDate.")) if any(l == "error" for l, _ in issues) or out["needs_time_unit"]: return out n_legacy_initial = n_date_derived = n_discrepant = n_pct_mismatch = 0 for li, cells in enumerate(rows[1:], start=1): if not any(str(c).strip() for c in cells): continue row_id = f"row-{li}" warnings = [] def get(role): i = mapping.get(role) return (cells[i].strip() if i is not None and i < len(cells) else "") taxon, accession = get("taxon"), get("accession") if not taxon or not accession: rejected.append({"row_id": row_id, "line": li + 1, "reason": "missing taxon or accession"}) continue st = normalize_storage(get("storage")) is_initial = st["legacy_initial"] if "is_initial" in mapping: v = get("is_initial").lower() if v in TRUTHY: is_initial = True elif v not in FALSY: warnings.append(f'unrecognized is_initial value "{get("is_initial")}" — treated as FALSE') if st["legacy_initial"]: n_legacy_initial += 1 tested = germinated = pct = None pct_supplied = _f(get("percent")) if has_pct else None if has_counts: tested, germinated = _f(get("tested")), _f(get("germinated")) if tested is None or germinated is None: rejected.append({"row_id": row_id, "line": li + 1, "reason": "missing counts"}); continue if tested <= 0: rejected.append({"row_id": row_id, "line": li + 1, "reason": f"number sown must be > 0 (got {tested:g})"}); continue if germinated < 0: rejected.append({"row_id": row_id, "line": li + 1, "reason": "germinated < 0"}); continue if germinated > tested: rejected.append({"row_id": row_id, "line": li + 1, "reason": f"germinated ({germinated:g}) > sown ({tested:g})"}); continue pct = 100.0 * germinated / tested if pct_supplied is not None and abs(pct_supplied - pct) > 2: warnings.append(f"supplied percent {pct_supplied:g}% ≠ recomputed {pct:.1f}% — counts are authoritative") n_pct_mismatch += 1 elif pct_supplied is not None: pct = min(100.0, max(0.0, pct_supplied)) # inspection only — no invented n else: rejected.append({"row_id": row_id, "line": li + 1, "reason": "no germination value"}); continue age_original = age_unit = None days = None if has_days and get("age_days") != "": age_original, age_unit = _f(get("age_days")), "days"; days = age_original elif has_years and get("age_years") != "": age_original, age_unit = _f(get("age_years")), "years" days = None if age_original is None else age_original * DAYS_PER_YEAR elif has_months and get("age_months") != "": age_original, age_unit = _f(get("age_months")), "months" days = None if age_original is None else age_original * DAYS_PER_MONTH elif amb_unit and get("age_ambiguous") != "": age_original = _f(get("age_ambiguous")); age_unit = amb_unit scale = 1.0 if amb_unit == "days" else (DAYS_PER_MONTH if amb_unit == "months" else DAYS_PER_YEAR) days = None if age_original is None else age_original * scale date_derived_years = None cd, td = _parse_date_loose(get("collection_date")), _parse_date_loose(get("test_date")) if cd and td: dd = (td - cd).days # exact calendar-day difference, never rounded if dd >= 0: date_derived_years = dd / DAYS_PER_YEAR if days is None: age_original, age_unit, days = dd, "dates", float(dd) n_date_derived += 1 elif abs(days / DAYS_PER_YEAR - date_derived_years) > 0.15: warnings.append(f"explicit age {days / DAYS_PER_YEAR:.2f} y disagrees with date-derived {date_derived_years:.2f} y — explicit age used") n_discrepant += 1 if days is None or days < 0: rejected.append({"row_id": row_id, "line": li + 1, "reason": "missing or invalid storage age"}); continue substrate, treatment = get("substrate"), get("treatment") s_norm, t_norm = norm_protocol_text(substrate), norm_protocol_text(treatment) assays.append({ "row_id": row_id, "assay_id": get("assay_id") or None, "family": get("family"), "taxon": taxon, "accession": accession, "is_initial": is_initial, "storage_code": "Initial" if is_initial else st["code"], "storage_original": st["original"], "drying": get("drying") or None, "substrate": substrate or None, "treatment": treatment or None, "substrate_norm": s_norm, "treatment_norm": t_norm, "protocol_key": protocol_key(substrate, treatment), "protocol_complete": s_norm != "(missing)" and t_norm != "(missing)", "tested": tested, "germinated": germinated, "pct": pct, "pct_supplied": pct_supplied, "age_original": age_original, "age_unit_original": age_unit, "days": days, "years": days / DAYS_PER_YEAR, "date_derived_years": date_derived_years, "collection_date": cd.isoformat() if cd else None, "test_date": td.isoformat() if td else None, "source_institution": get("source_institution") or None, "notes": get("notes") or None, "warnings": warnings, }) # ambiguous duplicate / multiple same-date tests — flag, never delete dup = {} for a in assays: date_key = a["test_date"] or f"age:{a['days']}" k = (a["taxon"], a["accession"], a["storage_code"], a["protocol_key"], date_key, a["tested"], a["germinated"]) dup.setdefault(k, []).append(a) for k, group in dup.items(): if len(group) < 2: continue aids = [a["assay_id"] for a in group] if all(x is not None for x in aids) and len(set(aids)) == len(aids): continue # distinct institutional IDs establish distinct tests out["duplicate_groups"].append({ "taxon": k[0], "accession": k[1], "storage": k[2], "protocol": k[3], "date_key": k[4], "row_ids": [a["row_id"] for a in group]}) if out["duplicate_groups"]: issues.append(("warn", f"Possible duplicate or multiple same-date tests: {len(out['duplicate_groups'])} group(s) of " "identical rows preserved — confirm whether they are distinct completed tests or duplicate records.")) if n_legacy_initial: issues.append(("info", f'{n_legacy_initial} row(s) used legacy storage="Initial" — mapped to is_initial=TRUE.')) if n_date_derived: issues.append(("info", f"{n_date_derived} row(s): storage age derived from dates as exact elapsed days (no rounding).")) if n_discrepant: issues.append(("warn", f"{n_discrepant} row(s): explicit age and date-derived age disagree by >0.15 y (explicit age used).")) if n_pct_mismatch: issues.append(("warn", f"{n_pct_mismatch} row(s): supplied percent disagrees with counts by >2 points — counts are authoritative.")) if rejected: reasons = "; ".join(sorted({r["reason"] for r in rejected})) issues.append(("warn", f"{len(rejected)} row(s) rejected: {reasons}.")) return out NORMALIZED_COLUMNS = [ "row_id", "assay_id", "family", "taxon", "accession", "is_initial", "storage_condition", "storage_condition_original", "drying_treatment", "substrate", "germination_treatment", "protocol_key", "tested", "germinated", "germination_percent_recomputed", "germination_percent_supplied", "storage_age_original", "storage_age_unit_original", "storage_age_days", "storage_age_years", "date_derived_years", "collection_date", "test_date", "source_institution", "notes", "warnings", ] def assays_to_normalized_rows(assays): rows = [] for a in assays: rows.append({ "row_id": a["row_id"], "assay_id": a["assay_id"] or "", "family": a["family"], "taxon": a["taxon"], "accession": a["accession"], "is_initial": "TRUE" if a["is_initial"] else "FALSE", "storage_condition": "" if a["is_initial"] else a["storage_code"], "storage_condition_original": a["storage_original"], "drying_treatment": a["drying"] or "", "substrate": a["substrate"] or "", "germination_treatment": a["treatment"] or "", "protocol_key": a["protocol_key"], "tested": "" if a["tested"] is None else f"{a['tested']:g}", "germinated": "" if a["germinated"] is None else f"{a['germinated']:g}", "germination_percent_recomputed": "" if a["pct"] is None else f"{a['pct']:.4f}", "germination_percent_supplied": "" if a["pct_supplied"] is None else f"{a['pct_supplied']:g}", "storage_age_original": "" if a["age_original"] is None else f"{a['age_original']:g}", "storage_age_unit_original": a["age_unit_original"] or "", "storage_age_days": f"{a['days']:.4f}", "storage_age_years": f"{a['years']:.6f}", "date_derived_years": "" if a["date_derived_years"] is None else f"{a['date_derived_years']:.6f}", "collection_date": a["collection_date"] or "", "test_date": a["test_date"] or "", "source_institution": a["source_institution"] or "", "notes": a["notes"] or "", "warnings": " | ".join(a["warnings"]), }) return rows def normalized_rows_to_assays(dict_rows): assays = [] for r in dict_rows: tested = _f(r.get("tested")) germinated = _f(r.get("germinated")) years = _f(r.get("storage_age_years")) days = _f(r.get("storage_age_days")) if years is None and days is not None: years = days / DAYS_PER_YEAR if years is None: continue is_initial = str(r.get("is_initial", "")).strip().lower() in TRUTHY substrate = r.get("substrate") or None treatment = r.get("germination_treatment") or None s_norm, t_norm = norm_protocol_text(substrate), norm_protocol_text(treatment) pct = _f(r.get("germination_percent_recomputed")) if pct is None and tested and germinated is not None and tested > 0: pct = 100.0 * germinated / tested assays.append({ "row_id": r.get("row_id") or f"row-{len(assays)+1}", "assay_id": r.get("assay_id") or None, "family": r.get("family") or "", "taxon": r.get("taxon") or "", "accession": r.get("accession") or "", "is_initial": is_initial, "storage_code": "Initial" if is_initial else (r.get("storage_condition") or r.get("storage_condition_original") or ""), "storage_original": r.get("storage_condition_original") or "", "drying": r.get("drying_treatment") or None, "substrate": substrate, "treatment": treatment, "substrate_norm": s_norm, "treatment_norm": t_norm, "protocol_key": r.get("protocol_key") or protocol_key(substrate, treatment), "protocol_complete": s_norm != "(missing)" and t_norm != "(missing)", "tested": tested, "germinated": germinated, "pct": pct, "pct_supplied": _f(r.get("germination_percent_supplied")), "age_original": _f(r.get("storage_age_original")), "age_unit_original": r.get("storage_age_unit_original") or None, "days": days if days is not None else years * DAYS_PER_YEAR, "years": years, "date_derived_years": _f(r.get("date_derived_years")), "collection_date": r.get("collection_date") or None, "test_date": r.get("test_date") or None, "source_institution": r.get("source_institution") or None, "notes": r.get("notes") or None, "warnings": [], }) return assays def load_input(path, time_unit=None): """Read EITHER a raw completed-test CSV or a normalized_assays.csv (this suite's interchange format — also what the website exports). Returns (assays, parse_result_or_None, was_normalized).""" with open(path, newline="", encoding="utf-8-sig") as fh: text = fh.read() first = text.lstrip("").splitlines()[0] if text.strip() else "" header_set = {h.strip() for h in first.split(",")} if NORMALIZED_SIGNATURE.issubset(header_set): with io.StringIO(text) as fh2: dict_rows = list(csv.DictReader(fh2)) return normalized_rows_to_assays(dict_rows), None, True parsed = parse_completed_tests(text, time_unit=time_unit) return parsed["assays"], parsed, False def write_csv(path, columns, dict_rows): with open(path, "w", newline="", encoding="utf-8") as fh: w = csv.DictWriter(fh, fieldnames=columns, extrasaction="ignore") w.writeheader() for r in dict_rows: w.writerow(r) def report_issues(parsed, strict=True): """Print parser issues to stderr. Exit 2 on blocking/percent-only when strict (analytical scripts); exit 1 on hard errors.""" if parsed is None: return for level, msg in parsed["issues"]: print(f"[{level.upper()}] {msg}", file=sys.stderr) if any(l == "error" for l, _ in parsed["issues"]) or parsed["needs_time_unit"]: sys.exit(1) if strict and parsed["percent_only"]: sys.exit(2) # --------------------------------------------------------------------------- # SHARED CORE · formal protocol-specific survival (identical block in # formal-dependent scripts). Port of the site's src/lib/formal.js, itself # built to the handoff's 13_protocol_specific_survival_v1_6.R + # 10_survival_plateau_qc semantics: # · trajectory = taxon × accession × storage × protocol_key (+ drying only # when supplied and actually variable) — never pooled across protocols # · matched same-protocol Initial (≤ 1 y) pooled once at t = 0, never # borrowed across protocols; older "Initials" excluded AND reported # · same-time rows pool by summed counts (identical binomial likelihood) # · plateau + early-rise QC; bridge category/tier; raw crossing kept, # ≤0 y displayed as "not interpretable (crossing ≤0 y)" # --------------------------------------------------------------------------- INITIAL_MAX_YEAR = 1.0 def group_formal_trajectories(assays): counted = [a for a in assays if a["tested"] is not None and a["germinated"] is not None] probe = {} for a in counted: if not a["drying"]: continue k = (a["taxon"], a["accession"], a["storage_code"], a["protocol_key"]) probe.setdefault(k, set()).add(a["drying"].strip().lower()) drying_keyed = any(len(s) > 1 for s in probe.values()) dry = (lambda a: (a["drying"] or "(none)").strip().lower()) if drying_keyed else (lambda a: "") post = [a for a in counted if not a["is_initial"]] initial = [a for a in counted if a["is_initial"] and a["years"] <= INITIAL_MAX_YEAR] excluded_late = [{"row_id": a["row_id"], "taxon": a["taxon"], "accession": a["accession"], "years": a["years"]} for a in counted if a["is_initial"] and a["years"] > INITIAL_MAX_YEAR] init_by = {} for a in initial: k = (a["taxon"], a["accession"], a["protocol_key"], dry(a)) e = init_by.setdefault(k, {"n_assays": 0, "sown": 0.0, "germ": 0.0}) e["n_assays"] += 1; e["sown"] += a["tested"]; e["germ"] += a["germinated"] traj_map = {} for a in post: k = (a["taxon"], a["accession"], a["storage_code"], a["protocol_key"], dry(a)) t = traj_map.setdefault(k, { "key": k, "taxon": a["taxon"], "accession": a["accession"], "family": a["family"], "storage_code": a["storage_code"], "protocol_key": a["protocol_key"], "substrate": a["substrate"], "treatment": a["treatment"], "protocol_complete": a["protocol_complete"], "drying": a["drying"] or None, "post_rows": []}) t["post_rows"].append(a) trajectories = [] for t in traj_map.values(): ik = (t["taxon"], t["accession"], t["protocol_key"], (t["drying"] or "(none)").strip().lower() if drying_keyed else "") init = init_by.get(ik) t["matched_initial"] = None if not init else { "n_assays": init["n_assays"], "sown": init["sown"], "germ": init["germ"], "pct": 100.0 * init["germ"] / init["sown"]} by_t = {} def add(tt, n, g, is_init): e = by_t.setdefault(tt, {"t": tt, "n": 0.0, "g": 0.0, "n_assays": 0, "includes_initial": False}) e["n"] += n; e["g"] += g; e["n_assays"] += 1 e["includes_initial"] = e["includes_initial"] or is_init if init: add(0.0, init["sown"], init["germ"], True) for a in t["post_rows"]: add(a["years"], a["tested"], a["germinated"], False) t["points"] = sorted(by_t.values(), key=lambda p: p["t"]) post_years = sorted({a["years"] for a in t["post_rows"]}) t["n_post_timepoints"] = len(post_years) t["max_post_years"] = post_years[-1] if post_years else None t["n_seeds"] = sum(p["n"] for p in t["points"]) early = [a["pct"] for a in t["post_rows"] if 0 < a["years"] <= 5] max_early = max(early) if early else None if t["matched_initial"] and max_early is not None: rise = max_early - t["matched_initial"]["pct"] t["early_rise_pp"] = rise t["early_rise_ge5"], t["early_rise_ge10"] = rise >= 5, rise >= 10 t["early_rise_class"] = ("substantial early rise (≥10 pp)" if rise >= 10 else "early rise (≥5 to <10 pp)" if rise >= 5 else "no ≥5 pp early rise detected") else: t["early_rise_pp"] = None t["early_rise_ge5"] = t["early_rise_ge10"] = False t["early_rise_class"] = ("not assessable: no same-protocol Initial" if not t["matched_initial"] else "not assessable: no post-storage test within 0–5 y") trajectories.append(t) trajectories.sort(key=lambda t: (t["taxon"], str(t["accession"]), t["storage_code"], t["protocol_key"])) return {"trajectories": trajectories, "drying_keyed": drying_keyed, "excluded_late_initials": excluded_late} def plateau_diagnostics(points, near100=99.0): s = sorted(points, key=lambda p: p["t"]) pct = [100.0 * p["g"] / p["n"] for p in s] lead = 0 for v in pct: if v >= near100: lead += 1 else: break trail = 0 for v in reversed(pct): if v <= 0: trail += 1 else: break note = ("leading ~100% plateau and trailing 0% plateau" if lead >= 2 and trail >= 2 else "leading ~100% plateau" if lead >= 2 else "trailing 0% plateau" if trail >= 2 else "none detected") return {"leading_near100": lead, "trailing_zero": trail, "plateau_candidate": lead >= 2 or trail >= 2, "plateau_note": note} def fit_formal_survival(traj, threshold_pct=None): pts = traj["points"] n_time = len({p["t"] for p in pts}) max_t = max((p["t"] for p in pts), default=None) plateau = plateau_diagnostics(pts) series_tier = ("insufficient (<2 distinct times)" if n_time < 2 else "exploratory (2 distinct times)" if n_time == 2 else "interpretable candidate (≥3 distinct times)") base = { "taxon": traj["taxon"], "accession": traj["accession"], "family": traj["family"], "storage_code": traj["storage_code"], "protocol_key": traj["protocol_key"], "substrate": traj["substrate"], "treatment": traj["treatment"], "protocol_complete": traj["protocol_complete"], "drying": traj["drying"], "points": pts, "fit_ok": False, "fit_method": None, "fit_reason": None, "fit_warning": None, "extreme_prob_warning": False, "converged": None, "series_tier": series_tier, "n_model_rows": len(pts), "n_timepoints": n_time, "n_post_timepoints": traj["n_post_timepoints"], "n_seeds": traj["n_seeds"], "max_time": max_t, "matched_initial_available": bool(traj["matched_initial"]), "matched_initial_pct": traj["matched_initial"]["pct"] if traj["matched_initial"] else None, "Ki_probit": None, "Ki_percent": None, "slope": None, "slope_se": None, "slope_p": None, "sigma_years": None, "p50_years": None, "p70_absolute_years": None, "threshold_pct": threshold_pct, "fitted_rci_p70_years_raw": None, "fitted_rci_p70_display": "not interpretable", "crossing_status": "model unavailable", "sparse_series": n_time < 3, **plateau, "early_rise_pp": traj["early_rise_pp"], "early_rise_ge5": traj["early_rise_ge5"], "early_rise_ge10": traj["early_rise_ge10"], "early_rise_class": traj["early_rise_class"], "bridge_category": "model unavailable", "bridge_evidence_tier": "do not interpret", "display_class": "NO FORMAL FIT", "qc_display_status": None, } def qc_display(res): if not res["protocol_complete"]: return "review: incomplete protocol metadata" if res["early_rise_ge10"]: return "review: substantial early rise; monotonic survival assumption questionable" if res["early_rise_ge5"]: return "review: early rise; interpret monotonic survival fit cautiously" if not res["matched_initial_available"] and res["bridge_evidence_tier"] == "primary model evidence": return "caution: within-range crossing but no same-protocol Initial" if not res["matched_initial_available"] and res["bridge_evidence_tier"] == "extrapolated model evidence": return "caution: extrapolated crossing and no same-protocol Initial" return res["bridge_evidence_tier"] if len(pts) < 2 or n_time < 2: base["fit_reason"] = "fewer than 2 distinct storage times" base["qc_display_status"] = qc_display(base) return base m = glm_probit([([1.0, p["t"]], p["n"], p["g"]) for p in pts]) if m is None: base["fit_reason"] = "probit fit failed" base["qc_display_status"] = qc_display(base) return base b0, b1 = m["coef"] ki_pct = normcdf(b0) * 100 declining = math.isfinite(b1) and b1 < 0 slope_se = math.sqrt(max(0.0, m["cov"][1][1])) if m["cov"] else None slope_p = 2 * (1 - normcdf(abs(b1 / slope_se))) if slope_se and slope_se > 0 else None sigma = -1 / b1 if declining else None p50 = crossing(b0, b1, 0.5) if declining and ki_pct >= 50 else None p70a = crossing(b0, b1, 0.7) if declining and ki_pct >= 70 else None ug = (crossing(b0, b1, threshold_pct / 100) if declining and threshold_pct is not None and 0 < threshold_pct < 100 else None) status = ("non-declining fitted slope; longevity crossing not interpreted" if not declining else "RCI-P70 threshold unavailable" if threshold_pct is None else "threshold crossing unavailable" if ug is None or not math.isfinite(ug) else "fitted viability already below RCI P70 at model time 0" if ug < 0 else "crossing within observed storage-time range" if ug <= max_t else "crossing extrapolated beyond observed storage-time range") cat = ("non-declining slope" if not math.isfinite(b1) or b1 >= 0 else "insufficient time points" if n_time < 2 else "crossing unavailable" if ug is None or not math.isfinite(ug) else "fitted viability below threshold at model time 0" if ug < 0 else "crossing within observed time range" if ug <= max_t else "crossing extrapolated beyond observed time range") tier = ("do not interpret" if not math.isfinite(b1) or b1 >= 0 or n_time < 2 or ug is None or not math.isfinite(ug) or ug < 0 else "exploratory only" if n_time == 2 else "primary model evidence" if ug <= max_t else "extrapolated model evidence") display = ("DO NOT INTERPRET" if tier == "do not interpret" else "EXPLORATORY" if tier == "exploratory only" else "FORMAL SURVIVAL") base.update({ "fit_ok": True, "fit_method": "irls-probit", "fit_reason": "ok" if declining else "non-declining slope", "fit_warning": ("possible separation: fit did not stabilize at finite estimates" if m["separated"] else "fitted probabilities numerically 0 or 1 occurred" if m["extreme_prob"] else None), "extreme_prob_warning": m["extreme_prob"], "converged": m["converged"], "Ki_probit": b0, "Ki_percent": ki_pct, "slope": b1, "slope_se": slope_se, "slope_p": slope_p, "sigma_years": sigma, "p50_years": p50, "p70_absolute_years": p70a, "fitted_rci_p70_years_raw": ug, "fitted_rci_p70_display": ("no RCI-P70 threshold (accession baseline unavailable)" if threshold_pct is None else crossing_display(ug)), "crossing_status": status, "bridge_category": cat, "bridge_evidence_tier": tier, "display_class": display, }) base["qc_display_status"] = qc_display(base) return base def thresholds_from_baselines(baselines): """(family, taxon, accession) → RCI-P70 threshold %, from the RCI branch's accession early baselines — the ONE genuine contingency of the formal workflow on the RCI stage.""" return {k: v["threshold"] for k, v in baselines.items()} # --------------------------------------------------------------------------- # SHARED CORE · competing-model comparison (identical block where used) # Port of the site's src/lib/models.js: four legitimate seed-longevity # models fit by maximum likelihood on the SAME binomial counts with one # optimizer, ranked by AICc with Akaike weights — the data picks the winner, # and a more complex model must EARN its parameters. # ellis-roberts v = Φ(Ki − t/σ) (k=2, field baseline) # control-viability v = C·Φ(Ki − t/σ) (k=3, Mead & Gray 1999) # resistant-fraction v = f + (1−f)·Φ(Ki − t/σ) (k=3, Hay 2014/Butler 2009) # avrami v = v0·exp(−(t/τ)^β) (k=3, Niedzielski 2009) # This asks which decay SHAPE explains one trajectory — a different question # from the shared-Ki temperature comparison and from the observed RCI. # --------------------------------------------------------------------------- def _binom_negll(points, predict): ll = 0.0 for t, n, g in points: if n <= 0: continue pr = min(1 - 1e-9, max(1e-9, predict(t))) ll += g * math.log(pr) + (n - g) * math.log(1 - pr) return -ll def _fit_mle(points, spec, factory): usable = [(t, n, g) for t, n, g in points if n > 0 and t >= 0 and 0 <= g <= n] if len(usable) < 2: return None p = [s["init"] for s in spec] step = [s["step"] for s in spec] best = _binom_negll(usable, factory(p)) if not math.isfinite(best): return None for _ in range(800): improved = False for i in range(len(p)): for d in (-1, 1): trial = p[:] trial[i] += d * step[i] if trial[i] < spec[i]["min"] or trial[i] > spec[i]["max"]: continue v = _binom_negll(usable, factory(trial)) if v < best - 1e-12: best, p = v, trial improved = True if not improved: step = [s * 0.5 for s in step] if all(s < spec[i]["tol"] for i, s in enumerate(step)): break return {"params": p, "negll": best, "n": len(usable), "seeds": sum(n for t, n, g in usable)} def _crossing_t(predict, thresh, max_t): lo, hi = 0.0, max_t if predict(lo) < thresh: return 0.0 if predict(hi) > thresh: hi = max_t * 8 if predict(hi) > thresh: return None for _ in range(60): mid = (lo + hi) / 2 if predict(mid) > thresh: lo = mid else: hi = mid return (lo + hi) / 2 def _seed_er(points): emp = [(t, probit(min(0.999, max(0.001, g / n)))) for t, n, g in points if n > 0] n = len(emp) sx = sum(e[0] for e in emp); sy = sum(e[1] for e in emp) sxx = sum(e[0] ** 2 for e in emp); sxy = sum(e[0] * e[1] for e in emp) den = n * sxx - sx * sx slope = (n * sxy - sx * sy) / (den if den else 1e-9) ki0 = (sy - slope * sx) / n if n else 1.5 sig0 = -1 / slope if slope < 0 else 50.0 return (ki0 if math.isfinite(ki0) else 1.5, sig0 if sig0 > 0 and math.isfinite(sig0) else 50.0) def compare_models(points): """points: [(t, n, g)] — time in any consistent unit. Returns ranked [{key, name, k, params{}, loglik, aicc, delta_aicc, akaike_weight, rank, p50, ...}], best first.""" ki0, sig0 = _seed_er(points) max_t = max((t for t, n, g in points), default=10) or 10 max_obs = max((g / n for t, n, g in points if n > 0), default=0.5) min_obs = min((g / n for t, n, g in points if n > 0), default=0.0) defs = [ ("ellis-roberts", "Ellis–Roberts", 2, [{"init": ki0, "min": -3, "max": 6, "step": 0.5, "tol": 1e-4}, {"init": sig0, "min": 1e-3, "max": 1e6, "step": max(1.0, sig0 * 0.5), "tol": 1e-3}], lambda v: (lambda t: normcdf(v[0] - t / v[1])), lambda v: {"Ki": v[0], "sigma": v[1]}), ("control-viability", "Control-viability", 3, [{"init": min(0.999, max_obs + 0.02), "min": 0.05, "max": 1.0, "step": 0.1, "tol": 1e-4}, {"init": 2.0, "min": -1, "max": 6, "step": 0.5, "tol": 1e-4}, {"init": 50.0, "min": 1e-3, "max": 1e6, "step": 25.0, "tol": 1e-3}], lambda v: (lambda t: v[0] * normcdf(v[1] - t / v[2])), lambda v: {"C": v[0], "Ki": v[1], "sigma": v[2]}), ("resistant-fraction", "Resistant-fraction", 3, [{"init": max(0.0, min(0.4, min_obs * 0.5)), "min": 0.0, "max": 0.9, "step": 0.1, "tol": 1e-4}, {"init": 2.0, "min": -1, "max": 6, "step": 0.5, "tol": 1e-4}, {"init": 50.0, "min": 1e-3, "max": 1e6, "step": 25.0, "tol": 1e-3}], lambda v: (lambda t: v[0] + (1 - v[0]) * normcdf(v[1] - t / v[2])), lambda v: {"f": v[0], "Ki": v[1], "sigma": v[2]}), ("avrami", "Avrami", 3, [{"init": min(0.999, max_obs), "min": 0.05, "max": 1.0, "step": 0.1, "tol": 1e-4}, {"init": max_t / 2, "min": 1e-3, "max": 1e7, "step": max_t / 4, "tol": 1e-3}, {"init": 1.0, "min": 0.2, "max": 6.0, "step": 0.4, "tol": 1e-4}], lambda v: (lambda t: v[0] * math.exp(-((t / v[1]) ** v[2]))), lambda v: {"v0": v[0], "tau": v[1], "beta": v[2]}), ] results = [] for key, name, k, spec, factory, named in defs: try: r = _fit_mle(points, spec, factory) except (ValueError, OverflowError): r = None if not r: continue predict = factory(r["params"]) loglik = -r["negll"] aic = 2 * k - 2 * loglik aicc = aic + (2 * k * (k + 1)) / max(1, r["n"] - k - 1) if not math.isfinite(aicc): continue results.append({"key": key, "name": name, "k": k, "params": named(r["params"]), "loglik": loglik, "aicc": aicc, "n": r["n"], "seeds": r["seeds"], "p50": _crossing_t(predict, 0.5, max(max_t, 10) * 10)}) if not results: return [] mn = min(r["aicc"] for r in results) rel = [math.exp(-0.5 * (r["aicc"] - mn)) for r in results] s = sum(rel) for i, r in enumerate(results): r["delta_aicc"] = r["aicc"] - mn r["akaike_weight"] = rel[i] / s results.sort(key=lambda r: r["delta_aicc"]) for i, r in enumerate(results): r["rank"] = i + 1 return results # ============================ command line ============================ MODEL_COLUMNS = [ "taxon", "accession", "storage_code", "protocol_key", "rank", "model", "k", "loglik", "aicc", "delta_aicc", "akaike_weight", "params", "p50_years", "n_points", "n_seeds", "verdict", ] def main(): import argparse, os, json ap = argparse.ArgumentParser(description=( "Component 5 — competing-model comparison (raw CSV or " "normalized_assays.csv in). Fits four published seed-longevity models " "(Ellis–Roberts, control-viability, resistant-fraction, Avrami) to " "EACH formal trajectory by maximum likelihood on the same binomial " "counts and ranks them by AICc with Akaike weights — the data picks " "the winner, and near-ties go to the simpler model. This asks which " "decay SHAPE fits one trajectory: a different question from the " "shared-Ki temperature comparison and from the observed RCI, and it " "overrides neither. Uses the formal component's trajectory grouping " "(embedded); independent of the RCI component.")) ap.add_argument("input") ap.add_argument("-o", "--outdir", default=".") ap.add_argument("--time-unit", choices=["days", "months", "years"]) args = ap.parse_args() assays, parsed, was_norm = load_input(args.input, time_unit=args.time_unit) report_issues(parsed) if not any(a["tested"] is not None for a in assays): print("[BLOCK] counts required — no sample size is ever assumed.", file=sys.stderr) sys.exit(2) g = group_formal_trajectories(assays) rows = [] for t in g["trajectories"]: pts = [(p["t"], p["n"], p["g"]) for p in t["points"]] ranked = compare_models(pts) if not ranked: continue tie = len(ranked) > 1 and ranked[1]["delta_aicc"] < 2 for r in ranked: verdict = "" if r["rank"] == 1: verdict = ("statistical tie with next model — parsimony favors the simpler form" if tie else "clear best fit by AICc") rows.append({ "taxon": t["taxon"], "accession": t["accession"], "storage_code": t["storage_code"], "protocol_key": t["protocol_key"], "rank": r["rank"], "model": r["key"], "k": r["k"], "loglik": f"{r['loglik']:.4f}", "aicc": f"{r['aicc']:.4f}", "delta_aicc": f"{r['delta_aicc']:.4f}", "akaike_weight": f"{r['akaike_weight']:.4f}", "params": json.dumps({k2: round(v, 6) for k2, v in r["params"].items()}), "p50_years": "" if r["p50"] is None else f"{r['p50']:.3f}", "n_points": r["n"], "n_seeds": int(r["seeds"]), "verdict": verdict, }) os.makedirs(args.outdir, exist_ok=True) p1 = os.path.join(args.outdir, "model_comparison.csv") write_csv(p1, MODEL_COLUMNS, rows) print(f"{p1}: {len(rows)} model fits across {len(g['trajectories'])} trajectories") if __name__ == "__main__": main()