// Analysis for the Vellora Research report "Smoking Prevalence and Lung Cancer
// Incidence Across Texas Counties". Pure and deterministic.
//
// Analysis plan, fixed before any association was computed:
//   Design: ecological, cross-sectional; unit = Texas county (canonical 5-digit FIPS).
//   Exposure: CDC PLACES 2025 release, current cigarette smoking among adults,
//     age-adjusted prevalence (CSMOKING, AgeAdjPrv; BRFSS 2023).
//   Outcome: State Cancer Profiles (Texas Cancer Registry), lung and bronchus
//     cancer incidence, age-adjusted rate per 100,000, 2019–2023 average.
//   Sample: counties with a published value for both; every other county is
//     listed with the source's reason. No imputation, no zeros.
//   Primary hypothesis (H1, the only confirmatory test): the two are positively
//     associated across counties. Spearman's ρ, two-sided p from the t
//     approximation, 95% CI from Fisher z with the Bonett–Wright standard
//     error, α = 0.05. One confirmatory test, so no multiplicity correction.
//   Exploratory, reported with unadjusted p and not used for the conclusion:
//     E1 Pearson r (Fisher z CI); E2 Spearman within metropolitan and within
//     nonmetropolitan counties (USDA RUCC 2023); E3 leave-one-out Spearman
//     range (influence); E4 crude smoking prevalence in place of age-adjusted;
//     E5 how excluded counties differ from included ones.
//   No county is removed as an outlier.
//
// Results are written to public/data/research/lung-smoking-results.v1.json by
// `npm run research:lung`; tests/lung-study.test.ts fails if that file is
// stale, and scripts/verify-lung-study.py recomputes it with SciPy.

import { describe, pearson, quantile } from "../stats";
import { spearmanCI } from "./cardio-study";

export const STUDY_ID = "lung-smoking-2026";

type Series = (number | null)[];
interface Cell {
  values: Series;
  lo: Series;
  hi: Series;
  reasons?: Record<string, string>;
}
export interface LungStudyInputs {
  countyOrder: string[];
  counties: { fips: string; name: string; metro: boolean; pop2020: number }[];
  smokingAA: Cell;
  smokingCrude: Cell;
  lung: Cell;
}

const round = (v: number, d = 4) => Number(v.toFixed(d));

function descriptives(v: number[], names: string[]) {
  const d = describe(v)!;
  const at = (x: number) => names.filter((_, i) => v[i] === x);
  return { n: d.n, mean: round(d.mean), sd: round(d.sd), median: round(d.median), q1: round(d.q1), q3: round(d.q3), min: d.min, max: d.max, minCounties: at(d.min), maxCounties: at(d.max) };
}

const spear = (x: number[], y: number[]) => {
  const c = spearmanCI(x, y);
  return { n: c.n, rho: round(c.rho), ci: [round(c.ci[0]), round(c.ci[1])] as [number, number], p: c.p! };
};

export function analyze(inp: LungStudyInputs) {
  const { countyOrder: fips, counties } = inp;
  if (counties.map((c) => c.fips).join() !== fips.join()) throw new Error("county order differs between inputs");
  const reasonOf = (cell: Cell, i: number) => cell.reasons?.[fips[i]] ?? "No value published";

  const included: number[] = [];
  const excluded: { fips: string; name: string; reasons: string[] }[] = [];
  fips.forEach((f, i) => {
    const reasons = [
      ...(inp.smokingAA.values[i] === null ? [`Smoking: ${reasonOf(inp.smokingAA, i)}`] : []),
      ...(inp.lung.values[i] === null ? [`Lung cancer: ${reasonOf(inp.lung, i)}`] : []),
    ];
    if (reasons.length) excluded.push({ fips: f, name: counties[i].name, reasons });
    else included.push(i);
  });

  const x = included.map((i) => inp.smokingAA.values[i]!);
  const y = included.map((i) => inp.lung.values[i]!);
  const names = included.map((i) => counties[i].name);
  const metro = included.map((i) => counties[i].metro);

  const p = pearson(x, y);
  const group = (m: boolean) => spear(x.filter((_, k) => metro[k] === m), y.filter((_, k) => metro[k] === m));
  const loo = included.map((_, k) => ({ county: names[k], rho: spearmanCI(x.filter((_, j) => j !== k), y.filter((_, j) => j !== k)).rho }));
  const looMin = loo.reduce((a, b) => (b.rho < a.rho ? b : a));
  const looMax = loo.reduce((a, b) => (b.rho > a.rho ? b : a));
  const crude = included.filter((i) => inp.smokingCrude.values[i] !== null);

  const exIdx = excluded.map((e) => fips.indexOf(e.fips));
  const popMedian = (idx: number[]) => quantile(idx.map((i) => counties[i].pop2020), 0.5);
  const smokingMedian = (idx: number[]) => quantile(idx.flatMap((i) => (inp.smokingAA.values[i] === null ? [] : [inp.smokingAA.values[i]!])), 0.5);

  return {
    study: STUDY_ID,
    parameters: { alpha: 0.05, primary: "Spearman rank correlation, two-sided, t approximation; 95% CI by Fisher z with Bonett–Wright SE", multipleTesting: "One confirmatory test; exploratory p-values unadjusted" },
    counties: { total: fips.length, included: included.length, excluded: excluded.length, includedMetro: metro.filter(Boolean).length, includedNonmetro: metro.filter((m) => !m).length },
    excluded,
    descriptives: { smoking: descriptives(x, names), lung: descriptives(y, names) },
    primary: spear(x, y),
    exploratory: {
      pearson: { n: p.n, r: round(p.r), ci: [round(p.ci![0]), round(p.ci![1])] as [number, number], p: p.p! },
      metro: group(true),
      nonmetro: group(false),
      leaveOneOut: { min: { county: looMin.county, rho: round(looMin.rho) }, max: { county: looMax.county, rho: round(looMax.rho) } },
      crudeSmoking: spear(crude.map((i) => inp.smokingCrude.values[i]!), crude.map((i) => inp.lung.values[i]!)),
      excludedProfile: {
        medianPop2020Included: popMedian(included),
        medianPop2020Excluded: popMedian(exIdx),
        maxPop2020Excluded: Math.max(...exIdx.map((i) => counties[i].pop2020)),
        metroExcluded: exIdx.filter((i) => counties[i].metro).length,
        medianSmokingIncluded: round(smokingMedian(included)),
        medianSmokingExcluded: round(smokingMedian(exIdx)),
      },
    },
    /** Per-county values behind the scatter plot (analytic sample only). */
    points: included.map((i, k) => ({
      fips: fips[i],
      name: names[k],
      metro: metro[k],
      smoking: x[k],
      smokingCI: [inp.smokingAA.lo[i], inp.smokingAA.hi[i]] as [number, number],
      lung: y[k],
      lungCI: [inp.lung.lo[i], inp.lung.hi[i]] as [number, number],
    })),
  };
}

export type LungStudyResults = ReturnType<typeof analyze>;