diff --git a/docs/api/openapi.yaml b/docs/api/openapi.yaml index 31aa35fe0..f7350a170 100644 --- a/docs/api/openapi.yaml +++ b/docs/api/openapi.yaml @@ -34740,6 +34740,9 @@ components: type: integer minimum: -9007199254740991 maximum: 9007199254740991 + nEff: + description: v1.42 — effective sample size after serial correlation, ≤ `n`; the p-value uses it. + type: number r: type: number pValue: @@ -62417,13 +62420,18 @@ components: type: integer minimum: -9007199254740991 maximum: 9007199254740991 - description: Paired-day count after the day+1 lag join (≥ 20). + description: Paired-day count after the lag join (≥ 20). + nEff: + description: v1.42 — effective sample size after serial correlation (Pyper–Peterman), ≤ `n`. The p-value, the shrinkage + and the tier use it. + type: number r: type: number - description: Pearson r over the lag-joined daily series. + description: "Pearson r over the lag-joined daily series, computed on residuals: a linear trend (and, from 120 paired + days, one annual harmonic) is removed from both series first (v1.42)." pValue: type: number - description: Two-sided exact Student-t p-value (< 0.05). + description: Two-sided exact Student-t p-value (< 0.05) at the effective sample size. qValue: type: number description: Benjamini-Hochberg FDR-adjusted q-value (≤ the surface threshold). @@ -62446,7 +62454,8 @@ components: type: integer minimum: -9007199254740991 maximum: 9007199254740991 - description: Lag in days applied (1). + description: "Lag in days applied: 1 = behaviour day, outcome the next day; 0 = the environmental channels, whose + exposure is the mean of the day before and the day itself, paired with the same day's outcome (v1.42)." window: description: v1.22 — which window surfaced the pair. Absent on the 180-day scan (retrospective default); `recent` for an emerging early-detection pair. diff --git a/scripts/compare-correlation-findings.ts b/scripts/compare-correlation-findings.ts new file mode 100644 index 000000000..886e5e0f4 --- /dev/null +++ b/scripts/compare-correlation-findings.ts @@ -0,0 +1,338 @@ +/** + * Read-only before/after count of correlation findings per account, for the + * release-candidate dry run of the v1.42 statistics change. + * + * v1.42 changes what the discovery engine accepts: trend and season are + * removed from both series before correlating, p-values use the effective + * sample size, shrinkage and the phrasing tier read that effective size, and + * the weather channels are tested once over the averaged day-before-and-same- + * day exposure instead of a next-day lag. Every one of those can make an + * existing finding disappear. This script answers "how many, and where", + * on a copy of production, before the release ships. + * + * For each account it builds the matrix exactly as the correlations route + * does (same assembler, same 180-day window, same module map, same tiered + * read) and scans it twice: + * + * - BEFORE: the pre-v1.42 rules, restated here from the engine's own + * exported pieces (lag join at 1 day, Pearson with df = n − 2, BH-FDR, + * shrinkage and tier by n), over the raw daily weather values; + * - AFTER: the engine as it now ships. + * + * It prints counts only. No correlation, p-value, metric value or account + * identifier leaves the database: accounts are listed by a short hash of + * their id, channels by their stable key family. Nothing is written — the + * pattern store is not synced, no cache is touched. + * + * Run (inside the app container, which has DATABASE_URL): + * pnpm dlx tsx scripts/compare-correlation-findings.ts [--limit N] + */ +import "dotenv/config"; + +import { createHash } from "node:crypto"; + +import { prisma } from "@/lib/db"; +import { resolveModuleMap } from "@/lib/modules/gate"; +import { assembleDiscoveryMatrix } from "@/lib/insights/discovery-matrix"; +import { MIN_PAIRED_N, pearson } from "@/lib/insights/correlations"; +import { + benjaminiHochberg, + confidenceTier, + discoverCorrelations, + discoverEmergingCorrelations, + EARLY_FDR_Q, + EARLY_MIN_PAIRS, + EARLY_WINDOW_DAYS, + FDR_Q, + filterSeriesToWindow, + lagJoin, + metricFamily, + shrinkEstimate, + type DailySeriesPoint, + type NamedSeries, +} from "@/lib/insights/correlation-discovery"; +import { ENVIRONMENT_FIELDS } from "@/lib/environment/fields"; +import { + isSurfaceVisible, + correlationChannelSurfaceId, +} from "@/lib/modules/surface"; +import { wallClockInTz } from "@/lib/tz/wall-clock"; +import { DEFAULT_TIMEZONE } from "@/lib/tz/format"; + +const MS_PER_DAY = 86_400_000; +/** The correlations route's window. */ +const WINDOW_DAYS = 180; + +interface Finding { + behaviour: string; + outcome: string; +} + +interface LegacyFinding extends Finding { + tier: ReturnType; +} + +/** The pre-v1.42 discovery rules, restated from the engine's parts. */ +function legacyDiscover( + series: NamedSeries[], + opts: { minPairs: number; fdrQ: number }, +): { found: LegacyFinding[]; pairsTested: number } { + const behaviours = series.filter((s) => s.role === "behaviour"); + const outcomes = series.filter((s) => s.role === "outcome"); + const tested: Array = []; + for (const b of behaviours) { + for (const o of outcomes) { + if (metricFamily(b.key) === metricFamily(o.key)) continue; + const { xs, ys } = lagJoin(b.points, o.points, 1); + if (xs.length < opts.minPairs) continue; + const result = pearson({ xs, ys, minPairs: opts.minPairs }); + if (result.status !== "ok") continue; + tested.push({ + behaviour: b.key, + outcome: o.key, + n: result.n, + r: result.r, + pValue: result.pValue, + }); + } + } + const q = benjaminiHochberg(tested.map((t) => t.pValue)); + const found = tested + .filter((t, i) => t.pValue < 0.05 && q[i] <= opts.fdrQ) + .map((t) => ({ + behaviour: t.behaviour, + outcome: t.outcome, + tier: confidenceTier(shrinkEstimate(t.r, t.n), t.n), + })) + .filter((t) => t.tier !== null); + return { found, pairsTested: tested.length }; +} + +/** The pre-v1.42 emerging pass: same rules, trailing window, no faint tier. */ +function legacyEmerging( + series: NamedSeries[], + retrospective: Finding[], + recentFromDayKey: string, +): Finding[] { + const recent = legacyDiscover( + filterSeriesToWindow(series, recentFromDayKey), + { minPairs: EARLY_MIN_PAIRS, fdrQ: EARLY_FDR_Q }, + ); + const established = new Set(retrospective.map(key)); + return recent.found + .filter((f) => f.tier !== "faint") + .filter((f) => !established.has(key(f))); +} + +function key(f: Finding): string { + return `${f.behaviour}␟${f.outcome}`; +} + +function tzDayKey(at: Date, tz: string): string { + const { year, month, day } = wallClockInTz(at, tz); + return `${year}-${String(month).padStart(2, "0")}-${String(day).padStart(2, "0")}`; +} + +/** + * The weather channels as the engine read them before v1.42: one raw daily + * value per field, no averaging. Read here directly because the shipped + * series now carries the averaged window. + */ +async function legacyEnvironmentSeries( + userId: string, + since: Date, +): Promise> { + // Same UTC lower bound the engine's own environment read uses. + const sinceKey = since.toISOString().slice(0, 10); + const rows = await prisma.environmentContext.findMany({ + where: { userId, date: { gte: sinceKey } }, + orderBy: { date: "asc" }, + take: 1000, + }); + const out = new Map(); + for (const field of ENVIRONMENT_FIELDS) { + const points: DailySeriesPoint[] = []; + for (const row of rows) { + const raw = row[field.column]; + if (raw == null || !Number.isFinite(raw)) continue; + const value = + field.column === "sunshineSec" || field.column === "daylightSec" + ? raw / 3600 + : raw; + points.push({ day: row.date, value }); + } + out.set(field.key, points); + } + return out; +} + +/** Channel family label for the aggregate table (no user data). */ +function familyLabel(channel: string): string { + if (channel.startsWith("ENV_")) return "ENVIRONMENT"; + if (channel.startsWith("CUSTOM_METRIC:")) return "CUSTOM_METRIC"; + if (channel.startsWith("SYMPTOM:")) return "SYMPTOM_EVENT"; + if (channel.startsWith("FACTOR:")) return "MOOD_FACTOR"; + return channel; +} + +interface Tally { + before: number; + after: number; +} + +function bump(map: Map, k: string, side: keyof Tally) { + const t = map.get(k) ?? { before: 0, after: 0 }; + t[side]++; + map.set(k, t); +} + +function parseLimit(argv: string[]): number | undefined { + const i = argv.indexOf("--limit"); + if (i === -1) return undefined; + const n = Number(argv[i + 1]); + return Number.isInteger(n) && n > 0 ? n : undefined; +} + +async function main(): Promise { + const limit = parseLimit(process.argv.slice(2)); + const users = await prisma.user.findMany({ + select: { id: true, timezone: true }, + orderBy: { createdAt: "asc" }, + ...(limit ? { take: limit } : {}), + }); + + const now = new Date(); + const since = new Date(now.getTime() - WINDOW_DAYS * MS_PER_DAY); + const byBehaviour = new Map(); + const byOutcome = new Map(); + const totals = { + accounts: 0, + scanned: 0, + retroBefore: 0, + retroAfter: 0, + recentBefore: 0, + recentAfter: 0, + kept: 0, + lost: 0, + gained: 0, + }; + + console.log( + "account tested(b/a) retro b→a recent b→a kept lost new env b→a", + ); + for (const user of users) { + totals.accounts++; + const tz = user.timezone ?? DEFAULT_TIMEZONE; + const modules = await resolveModuleMap(user.id); + const { series } = await assembleDiscoveryMatrix(user.id, { + tz, + since, + fetchMode: "tiered", + modules, + }); + const envVisible = isSurfaceVisible( + correlationChannelSurfaceId(ENVIRONMENT_FIELDS[0].key), + modules, + ); + const rawEnv = envVisible + ? await legacyEnvironmentSeries(user.id, since) + : new Map(); + const legacySeries: NamedSeries[] = series.map((s) => + s.key.startsWith("ENV_") + ? { key: s.key, role: s.role, points: rawEnv.get(s.key) ?? [] } + : { key: s.key, role: s.role, label: s.label, points: s.points }, + ); + + const recentFromDayKey = tzDayKey( + new Date(now.getTime() - EARLY_WINDOW_DAYS * MS_PER_DAY), + tz, + ); + + const legacy = legacyDiscover(legacySeries, { + minPairs: MIN_PAIRED_N, + fdrQ: FDR_Q, + }); + const before = legacy.found; + const beforeRecent = legacyEmerging(legacySeries, before, recentFromDayKey); + const after = discoverCorrelations(series, { locale: "en" }); + const afterRecent = discoverEmergingCorrelations(series, after, { + recentFromDayKey, + locale: "en", + }); + + if (legacy.pairsTested === 0 && after.pairsTested === 0) continue; + totals.scanned++; + + const beforeKeys = new Set(before.map(key)); + const afterKeys = new Set(after.discovered.map(key)); + const kept = [...afterKeys].filter((k) => beforeKeys.has(k)).length; + const lost = beforeKeys.size - kept; + const gained = afterKeys.size - kept; + const envBefore = before.filter((f) => f.behaviour.startsWith("ENV_")); + const envAfter = after.discovered.filter((f) => + f.behaviour.startsWith("ENV_"), + ); + + for (const f of before) { + bump(byBehaviour, familyLabel(f.behaviour), "before"); + bump(byOutcome, familyLabel(f.outcome), "before"); + } + for (const f of after.discovered) { + bump(byBehaviour, familyLabel(f.behaviour), "after"); + bump(byOutcome, familyLabel(f.outcome), "after"); + } + + totals.retroBefore += before.length; + totals.retroAfter += after.discovered.length; + totals.recentBefore += beforeRecent.length; + totals.recentAfter += afterRecent.emerging.length; + totals.kept += kept; + totals.lost += lost; + totals.gained += gained; + + const tag = createHash("sha256").update(user.id).digest("hex").slice(0, 8); + console.log( + [ + tag.padEnd(9), + `${legacy.pairsTested}/${after.pairsTested}`.padEnd(12), + `${before.length}→${after.discovered.length}`.padEnd(10), + `${beforeRecent.length}→${afterRecent.emerging.length}`.padEnd(11), + String(kept).padStart(4), + String(lost).padStart(4), + String(gained).padStart(3), + ` ${envBefore.length}→${envAfter.length}`, + ].join(" "), + ); + } + + console.log(""); + console.log( + `accounts ${totals.accounts}, with a testable matrix ${totals.scanned}`, + ); + console.log( + `retrospective findings ${totals.retroBefore} → ${totals.retroAfter} (kept ${totals.kept}, lost ${totals.lost}, new ${totals.gained})`, + ); + console.log( + `emerging findings ${totals.recentBefore} → ${totals.recentAfter}`, + ); + const table = (title: string, map: Map) => { + console.log(""); + console.log(`${title}: before → after`); + for (const [k, t] of [...map.entries()].sort((a, b) => + a[0].localeCompare(b[0]), + )) { + console.log(` ${k.padEnd(24)} ${t.before} → ${t.after}`); + } + }; + table("by behaviour family", byBehaviour); + table("by outcome family", byOutcome); +} + +main() + .catch((err: unknown) => { + console.error(err); + process.exitCode = 1; + }) + .finally(() => { + void prisma.$disconnect(); + }); diff --git a/src/__tests__/i18n-reverse-coverage.test.ts b/src/__tests__/i18n-reverse-coverage.test.ts index 14e8c6832..e3bda6a41 100644 --- a/src/__tests__/i18n-reverse-coverage.test.ts +++ b/src/__tests__/i18n-reverse-coverage.test.ts @@ -91,8 +91,6 @@ const RELEASE_PENDING_PREFIXES: Record = { "environment.fields.europeanAqi": "v1.42 air-quality field label (#615)", "environment.fields.uvIndexMax": "v1.42 air-quality field label (#615)", "environment.fields.dustMax": "v1.42 air-quality field label (#615)", - "insights.correlation.dailySameDay": - "v1.42 same-day environment correlation templates (#615)", "settings.appleHealth.freshness": "v1.42 per-type HealthKit freshness list (#1173)", "insights.workouts.manual.editTitle": "v1.42 manual workout edit (#1162)", diff --git a/src/lib/ai/coach/tools/correlations-read.ts b/src/lib/ai/coach/tools/correlations-read.ts index 7c49af712..3c050a535 100644 --- a/src/lib/ai/coach/tools/correlations-read.ts +++ b/src/lib/ai/coach/tools/correlations-read.ts @@ -75,9 +75,12 @@ const COINCIDENT_WINDOW_DAYS = 30; export interface CoachCorrelationDriver { behaviour: string; outcome: string; - /** "higher" / "lower" — sign of the next-day association. */ + /** "higher" / "lower" — sign of the association. */ direction: "higher" | "lower"; - /** Lag in days (always 1 today). */ + /** + * Lag in days: 1 = next day; 0 = an environmental exposure averaged over + * the day before and the day itself, against the same day's outcome. + */ lagDays: number; /** Paired-day sample count after the lag join. */ n: number; diff --git a/src/lib/environment/fields.ts b/src/lib/environment/fields.ts index bb313fbe3..2c08e7f32 100644 --- a/src/lib/environment/fields.ts +++ b/src/lib/environment/fields.ts @@ -7,7 +7,7 @@ * - the SIGNAL REGISTRY (`src/lib/signals/registry.ts`) registers one * `kind:"environment"` signal per field here, so each gets a stable signal * key + display metadata for free; and - * - the CORRELATION ENGINE reads {@link ENVIRONMENT_CHANNEL_KEYS} to fold the + * - the CORRELATION ENGINE reads {@link ENVIRONMENT_FIELDS} to fold the * same fields in as lagged BEHAVIOUR (exposure) channels against mood / * sleep / vitals. * @@ -20,6 +20,28 @@ * `column` is the `EnvironmentContext` numeric column the field reads. Leaf * module (no imports) so both the registry and the pure correlation engine can * depend on it without a cycle. + * + * Lag semantics (v1.42). Every field here is tested as ONE hypothesis per + * outcome: the exposure is the mean of the day before and the day itself + * ("lag 0–1"), paired with the outcome on the same day. The correlation read + * in `correlation-channel-series.ts` builds that mean, emits a point only for + * a day whose previous day is stored too, and tags the series `lagDays: 0`; a + * discovered pair therefore reports `lagDays: 0` and narrates "on the day and + * the day before". Two reasons, both about getting the night right: + * + * - sleep is keyed on the day you wake up, so last night's sleep sits on day + * D, and the overnight low that shaped it is the early morning of D (and the + * evening of D−1). The old fixed next-day lag paired the night's low with + * the FOLLOWING night; + * - testing lag 0 and lag 1 separately would double the environmental tests + * in the false-discovery family (and taking whichever looked better would + * be selection without correction). The averaged window covers both at the + * cost of one test. + * + * Season and trend are removed from both series before any of this is + * correlated (see `src/lib/insights/seasonal-adjust.ts`), so a field that only + * follows the calendar does not surface. A new field added below inherits all + * of it with no further wiring. */ /** The `EnvironmentContext` numeric columns a correlation field can read. */ @@ -115,5 +137,3 @@ export const ENVIRONMENT_FIELDS: readonly EnvironmentField[] = [ i18nLabelKey: "environment.fields.pressureDelta", }, ] as const; - -/** The correlation channel keys, in declaration order. */ diff --git a/src/lib/insights/__tests__/correlation-channel-series.test.ts b/src/lib/insights/__tests__/correlation-channel-series.test.ts index e2026f620..f481d058e 100644 --- a/src/lib/insights/__tests__/correlation-channel-series.test.ts +++ b/src/lib/insights/__tests__/correlation-channel-series.test.ts @@ -4,15 +4,21 @@ import type { SleepStage } from "@/generated/prisma/client"; const measurementFindMany = vi.fn(); const moodFindMany = vi.fn(); const customMetricFindMany = vi.fn(); +const environmentFindMany = vi.fn(); vi.mock("@/lib/db", () => ({ prisma: { measurement: { findMany: (a: unknown) => measurementFindMany(a) }, moodEntry: { findMany: (a: unknown) => moodFindMany(a) }, customMetric: { findMany: (a: unknown) => customMetricFindMany(a) }, + environmentContext: { + findMany: (a: unknown) => environmentFindMany(a), + }, }, })); import { + averageWithPreviousDay, + fetchEnvironmentSeries, buildMeasurementDailySeries, fetchMeasurementWindowSeries, fetchCustomMetricBehaviourSeries, @@ -345,3 +351,51 @@ describe("fetchCustomMetricBehaviourSeries", () => { ]); }); }); + +describe("environment exposure: the averaged lag 0–1 window (v1.42)", () => { + it("averages each day with the day before, and only where both exist", () => { + const out = averageWithPreviousDay([ + { day: "2026-05-03", value: 14 }, + { day: "2026-05-01", value: 10 }, + { day: "2026-05-02", value: 12 }, + // 05-04 missing: neither 05-04 nor 05-05 gets a point. + { day: "2026-05-05", value: 20 }, + { day: "2026-05-06", value: 22 }, + ]); + expect(out).toEqual([ + { day: "2026-05-02", value: 11 }, + { day: "2026-05-03", value: 13 }, + { day: "2026-05-06", value: 21 }, + ]); + }); + + it("crosses a month and a year boundary on calendar days", () => { + expect( + averageWithPreviousDay([ + { day: "2025-12-31", value: 2 }, + { day: "2026-01-01", value: 4 }, + ]), + ).toEqual([{ day: "2026-01-01", value: 3 }]); + }); + + it("tags every environment channel lagDays 0 over averaged points", async () => { + environmentFindMany.mockResolvedValueOnce([ + { date: "2026-05-01", tempMin: 8, tempMean: 12, sunshineSec: 3600 }, + { date: "2026-05-02", tempMin: 10, tempMean: 14, sunshineSec: 7200 }, + ]); + const series = await fetchEnvironmentSeries( + "u1", + new Date("2026-04-01T00:00:00Z"), + ); + expect(series.length).toBeGreaterThan(0); + for (const s of series) { + expect(s.lagDays, s.key).toBe(0); + expect(s.role).toBe("behaviour"); + } + const tMin = series.find((s) => s.key === "ENV_TEMP_MIN")!; + // The night before 05-02 and its morning: (8 + 10) / 2, keyed on 05-02. + expect(tMin.points).toEqual([{ day: "2026-05-02", value: 9 }]); + const sun = series.find((s) => s.key === "ENV_SUNSHINE")!; + expect(sun.points).toEqual([{ day: "2026-05-02", value: 1.5 }]); + }); +}); diff --git a/src/lib/insights/__tests__/correlation-discovery.test.ts b/src/lib/insights/__tests__/correlation-discovery.test.ts index 8b09fdc4c..e826fe004 100644 --- a/src/lib/insights/__tests__/correlation-discovery.test.ts +++ b/src/lib/insights/__tests__/correlation-discovery.test.ts @@ -7,6 +7,7 @@ import { type DailySeriesPoint, type NamedSeries, } from "../correlation-discovery"; +import { coupledPair } from "./helpers/seeded-series"; /** Build a contiguous daily series from day-1 of a month. */ function series(values: number[], startDay = 1): DailySeriesPoint[] { @@ -327,9 +328,10 @@ describe("D2-2 — effect-size floor on discovered drivers", () => { }); it("a strong deep pair is narrated with confident phrasing (high tier)", () => { - const n = 90; - const behaviour = Array.from({ length: n }, (_, i) => i + (i % 4)); - const outcome = [0, ...behaviour.slice(0, n - 1).map((v) => v * 2 + 5)]; + // Seeded noise with a real next-day coupling (r ≈ 0.9). A deterministic + // periodic fixture would not do here: its autocorrelation is perfect, so + // the effective sample size collapses and the pair is honestly thin. + const { behaviour, outcome } = coupledPair(11, 90, 2); const result = discoverCorrelations( [ { @@ -352,21 +354,15 @@ describe("D2-2 — effect-size floor on discovered drivers", () => { }); it("a faint-tier pair is hedged, never confident", () => { + // A seed whose sample effect lands mid-band (shrunk r ≈ 0.25). + const FAINT_SEED = 2; // Construct a pair whose SHRUNK r lands in [floor, confident): a moderate // raw r on a sample deep enough to clear significance but whose shrunk // magnitude stays under the confident threshold. - const n = 60; - // Target shrunk r ≈ 0.25 ⇒ raw r ≈ 0.25 * (n+10)/n ≈ 0.29. - const behaviour = Array.from( - { length: n }, - (_, i) => Math.sin(i * 0.9) * 10 + 50, - ); - const outcome = [ - 0, - ...behaviour - .slice(0, n - 1) - .map((v, i) => v * 0.35 + Math.cos(i * 2.1) * 9 + 100), - ]; + // Seeded noise with a weak next-day coupling (slope 0.3, r ≈ 0.3) over + // 200 days: significant, but the shrunk effect stays under the confident + // threshold. + const { behaviour, outcome } = coupledPair(FAINT_SEED, 200, 0.3); const result = discoverCorrelations( [ { @@ -382,7 +378,7 @@ describe("D2-2 — effect-size floor on discovered drivers", () => { (p) => p.behaviour === "TIME_IN_DAYLIGHT" && p.outcome === "SLEEP_DURATION", ); - // Deterministic: raw r=0.35, shrunk to ≈0.30 by n=59 → faint tier. + // Deterministic for the fixed seed: shrunk r in [0.2, 0.3) → faint tier. expect(pair).toBeDefined(); expect(pair!.tier).toBe("faint"); expect(Math.abs(pair!.shrunkR)).toBeGreaterThanOrEqual(EFFECT_SIZE_FLOOR); @@ -401,9 +397,7 @@ describe("D2-2 — effect-size floor on discovered drivers", () => { // pins the fix: with `locale: "de"` the interpretation is the German template // (localised metric label + correctly-negated "keine Ursache"), NOT English. it("narrates in the reader's locale — German, never the English string", () => { - const n = 90; - const behaviour = Array.from({ length: n }, (_, i) => i + (i % 4)); - const outcome = [0, ...behaviour.slice(0, n - 1).map((v) => v * 2 + 5)]; + const { behaviour, outcome } = coupledPair(11, 90, 2); const seriesInput: NamedSeries[] = [ { key: "TIME_IN_DAYLIGHT", @@ -474,21 +468,12 @@ describe("compliance + symptom channels (FDREXTEND)", () => { }); it("surfaces an adherence-dip → next-day symptom-flare link at a confident tier", () => { - // 70 contiguous days. Lower adherence on day D drives a higher symptom - // burden on day D+1 (the flagship cross-metric link). Adherence oscillates - // 100/70/40; symptom[d+1] = (100 - adherence[d]) scaled into 0..3. - const n = 70; - const adherence = Array.from({ length: n }, (_, i) => - i % 3 === 0 ? 100 : i % 3 === 1 ? 70 : 40, - ); - const symptom = [ - 0, - ...adherence - .slice(0, n - 1) - // (100 - a) maps 0..60 → ~0..3, plus a tiny deterministic jitter so the - // pair is not a perfect line (a real-world-shaped strong link). - .map((a, i) => ((100 - a) / 20) * (i % 2 === 0 ? 1.0 : 0.98)), - ]; + // 120 contiguous days. Lower adherence on day D drives a higher symptom + // burden on day D+1 (the flagship cross-metric link). Adherence varies day to day (seeded); the symptom burden the day after + // falls as adherence rises, with independent noise on top. + const pair0 = coupledPair(23, 120, -2); + const adherence = pair0.behaviour.map((v) => 70 + v * 10); + const symptom = pair0.outcome.map((v) => (v - 40) / 7); const result = discoverCorrelations( [ { @@ -717,8 +702,9 @@ describe("discoverEmergingCorrelations (early detection)", () => { // retrospective scan needs ≥ 20 lagged pairs, so it sees NOTHING; the early // window (floor 10) catches it — exactly the early-detection mandate. const len = 18; - const beh = Array.from({ length: len }, (_, i) => (i % 2 ? 2 : 9)); - const out = beh.map((b, i) => 40 + b * 3 + (i % 2 ? 0.3 : -0.3)); + const { behaviour: beh, outcome: out } = coupledPair(5, len, 3, { + noiseSd: 0.5, + }); const series: NamedSeries[] = [ { key: "MOOD", role: "behaviour", points: longSeries(beh) }, { @@ -878,3 +864,100 @@ describe("discoverLabOutcomeCorrelations (labs ↔ outcomes)", () => { expect(result.pairsTested).toBe(0); }); }); + +describe("v1.42 — per-series lag and the averaged environment window", () => { + it("pairs a lagDays-0 exposure with the same day and narrates it so", () => { + // Exposure (already the D−1/D mean upstream) drives the outcome on the + // SAME day; the default next-day join would see noise. + const { behaviour, outcome } = coupledPair(31, 120, 2); + const sameDayOutcome = outcome.slice(1); + const exposure = behaviour.slice(0, -1); + const result = discoverCorrelations( + [ + { + key: "ENV_TEMP_MIN", + role: "behaviour", + points: longSeries(exposure), + lagDays: 0, + }, + { + key: "SLEEP_DURATION", + role: "outcome", + points: longSeries(sameDayOutcome), + }, + ], + { locale: "en" }, + ); + expect(result.pairsTested).toBe(1); + const pair = result.discovered[0]; + expect(pair).toBeDefined(); + expect(pair.lagDays).toBe(0); + expect(pair.nEff).toBeGreaterThan(0); + expect(pair.nEff).toBeLessThanOrEqual(pair.n); + expect(pair.interpretation).toMatch(/on the day and the day before/); + expect(pair.interpretation).not.toMatch(/next-day/); + }); + + it("keeps the next-day lag and wording for every other channel", () => { + const { behaviour, outcome } = coupledPair(11, 90, 2); + const result = discoverCorrelations( + [ + { + key: "ACTIVITY_STEPS", + role: "behaviour", + points: longSeries(behaviour), + }, + { key: "SLEEP_DURATION", role: "outcome", points: longSeries(outcome) }, + ], + { locale: "en" }, + ); + expect(result.discovered[0].lagDays).toBe(1); + expect(result.discovered[0].interpretation).toMatch(/next-day/); + }); + + it("counts each exposure × outcome once in the false-discovery family", () => { + const { behaviour, outcome } = coupledPair(13, 90, 0); + const result = discoverCorrelations( + [ + { + key: "ENV_PRESSURE_MEAN", + role: "behaviour", + points: longSeries(behaviour), + lagDays: 0, + }, + { key: "MOOD", role: "outcome", points: longSeries(outcome) }, + { + key: "SLEEP_DURATION", + role: "outcome", + points: longSeries(outcome), + }, + ], + { locale: "en" }, + ); + expect(result.pairsTested).toBe(2); + }); + + it("carries a series' own lag into the emerging pass", () => { + const { behaviour, outcome } = coupledPair(17, 40, 3, { noiseSd: 0.5 }); + const series: NamedSeries[] = [ + { + key: "ENV_TEMP_MEAN", + role: "behaviour", + points: longSeries(behaviour.slice(0, -1)), + lagDays: 0, + }, + { + key: "HEART_RATE_VARIABILITY", + role: "outcome", + points: longSeries(outcome.slice(1)), + }, + ]; + const empty = { discovered: [], pairsTested: 0, fdrQ: 0.1, minPairs: 20 }; + const { emerging } = discoverEmergingCorrelations(series, empty, { + recentFromDayKey: longSeries(behaviour)[0].day, + locale: "en", + minPairs: 10, + }); + expect(emerging[0]?.lagDays).toBe(0); + }); +}); diff --git a/src/lib/insights/__tests__/correlation-seasonal-simulation.test.ts b/src/lib/insights/__tests__/correlation-seasonal-simulation.test.ts new file mode 100644 index 000000000..ed06e317f --- /dev/null +++ b/src/lib/insights/__tests__/correlation-seasonal-simulation.test.ts @@ -0,0 +1,90 @@ +/** + * The discovery engine against the seeded season simulation, end to end. + * + * Where `seasonal-adjust.test.ts` pins the statistic, this file pins what a + * person actually sees: does `discoverCorrelations` SURFACE a finding. It + * deliberately touches nothing but the engine's long-standing public surface + * (`discoverCorrelations` over two `NamedSeries`), so the same file runs + * against the engine before the adjustment existed — and fails there. That is + * the proof these tests can fail: before v1.42, two series that shared only + * the calendar surfaced as related in most runs. + * + * A finding here is a pair in `discovered`: p < 0.05, BH q ≤ 0.10, and a + * shrunk effect at or above the 0.2 floor. With one pair per run the BH step + * is the identity, so the rate is the engine's own false-positive rate. + */ +import { describe, expect, it } from "vitest"; + +import { + discoverCorrelations, + type DailySeriesPoint, +} from "../correlation-discovery"; +import { + seasonalScenario, + type SeasonalScenario, +} from "./helpers/seeded-series"; + +const RUNS = 200; + +function surfaces( + temperature: DailySeriesPoint[], + systolic: DailySeriesPoint[], +): boolean { + const result = discoverCorrelations( + [ + { key: "ACTIVITY_STEPS", role: "behaviour", points: temperature }, + { key: "BLOOD_PRESSURE_SYS", role: "outcome", points: systolic }, + ], + { locale: "en" }, + ); + return result.discovered.length > 0; +} + +/** Share of seeded runs in which the engine surfaces the pair. */ +function surfacedRate(scenario: SeasonalScenario): number { + let found = 0; + for (let seed = 1; seed <= RUNS; seed++) { + const { temperature, systolic } = seasonalScenario(5000 + seed, scenario); + if (surfaces(temperature, systolic)) found++; + } + return found / RUNS; +} + +describe("discovery engine: series that share only the season", () => { + it("surfaces no more than ~6 % over 400 days from 1 January", () => { + expect(surfacedRate({ days: 400, startDoy: 0 })).toBeLessThanOrEqual(0.06); + }); + + it("surfaces no more than ~6 % on the 180-day window, any time of year", () => { + expect(surfacedRate({ days: 180, startDoy: "random" })).toBeLessThanOrEqual( + 0.06, + ); + }); + + it("surfaces no more than ~6 % over a 180-day spring (rising heat, falling pressure)", () => { + expect(surfacedRate({ days: 180, startDoy: 32 })).toBeLessThanOrEqual(0.06); + }); + + it("surfaces no more than ~8 % on a short 60-day series (trend-only path)", () => { + // Under 120 days only the trend is removed. The effective-n estimate is + // noisier on 59 pairs: over 1 000 seeds this scenario surfaces 6.0 %, and + // 200 seeds sample that with about ±1.7 points, so the bound sits at the + // research plan's 8 % rather than at the long-series 6 %. + expect(surfacedRate({ days: 60, startDoy: 60 })).toBeLessThanOrEqual(0.08); + }); +}); + +describe("discovery engine: a real next-day effect under the season", () => { + it("finds it in most runs over 180 days", () => { + // beta 1.5 mmHg per K of anomaly: r ≈ 0.5 on the residuals. + expect( + surfacedRate({ days: 180, startDoy: 32, beta: 1.5 }), + ).toBeGreaterThan(0.8); + }); + + it("finds it in most runs over 400 days", () => { + expect(surfacedRate({ days: 400, startDoy: 0, beta: 1.5 })).toBeGreaterThan( + 0.9, + ); + }); +}); diff --git a/src/lib/insights/__tests__/correlations.test.ts b/src/lib/insights/__tests__/correlations.test.ts index 168feb44f..e7f3d12ed 100644 --- a/src/lib/insights/__tests__/correlations.test.ts +++ b/src/lib/insights/__tests__/correlations.test.ts @@ -8,6 +8,7 @@ import { MIN_PAIRED_N, MAX_P_VALUE, } from "../correlations"; +import { seededRng } from "./helpers/seeded-series"; /** * v1.4.20 phase B3 — correlation discovery. @@ -265,11 +266,14 @@ describe("correlateBpCompliance", () => { }); it("returns ok with negative r for strong compliance↔BP support", () => { + // Day-to-day: a higher-compliance day has a lower systolic (seeded noise, + // no shared trend), so the relation survives detrending. + const rng = seededRng(7); const daily = buildPairs( - Array.from( - { length: 20 }, - (_, i) => [50 + i * 2.5, 150 - i * 1.5] as [number, number], - ), + Array.from({ length: 20 }, () => { + const c = 50 + rng.n() * 20; + return [c, 190 - c * 0.8 + rng.n() * 2] as [number, number]; + }), ); const result = correlateBpCompliance({ daily }); expect(result.status).toBe("ok"); @@ -283,6 +287,20 @@ describe("correlateBpCompliance", () => { } }); + it("does not report two series that only share a trend (v1.42)", () => { + // Compliance creeping up while systolic creeps down over the same weeks is + // a shared drift, not a day-to-day relation. Before v1.42 this pair read + // r ≈ −1; the trend is now removed from both arms first. + const daily = buildPairs( + Array.from( + { length: 20 }, + (_, i) => [50 + i * 2.5, 150 - i * 1.5] as [number, number], + ), + ); + const result = correlateBpCompliance({ daily }); + expect(result.status).toBe("insufficient"); + }); + it("returns insufficient (not_significant) when r is near zero", () => { // Deterministic jitter, n large enough — but no relation. const daily = buildPairs([ diff --git a/src/lib/insights/__tests__/discovery-matrix-modules.test.ts b/src/lib/insights/__tests__/discovery-matrix-modules.test.ts index 930b52b27..8f27d2c34 100644 --- a/src/lib/insights/__tests__/discovery-matrix-modules.test.ts +++ b/src/lib/insights/__tests__/discovery-matrix-modules.test.ts @@ -54,6 +54,7 @@ import { assembleDiscoveryMatrix, maskSeriesByModules, } from "@/lib/insights/discovery-matrix"; +import { fetchEnvironmentSeries } from "@/lib/insights/correlation-channel-series"; const OPTS = { tz: "UTC", @@ -119,6 +120,18 @@ describe("assembleDiscoveryMatrix — switched-off modules", () => { expect(keys).toContain("CUSTOM_METRIC:abc"); }); + it("does not read the weather at all with the environment module off", async () => { + const read = vi.mocked(fetchEnvironmentSeries); + read.mockClear(); + await assembleDiscoveryMatrix("u1", { + ...OPTS, + modules: { environment: false }, + }); + expect(read).not.toHaveBeenCalled(); + await assembleDiscoveryMatrix("u1", OPTS); + expect(read).toHaveBeenCalledTimes(1); + }); + it("drops nothing for AI analysis off", async () => { const all = await assembleDiscoveryMatrix("u1", OPTS); const off = await assembleDiscoveryMatrix("u1", { diff --git a/src/lib/insights/__tests__/helpers/seeded-series.ts b/src/lib/insights/__tests__/helpers/seeded-series.ts new file mode 100644 index 000000000..7720e7672 --- /dev/null +++ b/src/lib/insights/__tests__/helpers/seeded-series.ts @@ -0,0 +1,164 @@ +/** + * Deterministic synthetic daily series for the correlation tests. + * + * A fixed-seed generator (mulberry32 + Box–Muller), never `Math.random`, so + * every simulation run is reproducible and a regression shows up as the same + * number every time. The fixtures mirror the shapes real data has: AR(1) + * noise (today resembles yesterday), an optional annual cycle, gaps. + */ +import type { DailySeriesPoint } from "../../correlation-discovery"; + +export interface Rng { + /** Uniform in [0, 1). */ + u: () => number; + /** Standard normal. */ + n: () => number; +} + +/** mulberry32 with a Box–Muller normal on top. */ +export function seededRng(seed: number): Rng { + let a = seed >>> 0; + const u = () => { + a = (a + 0x6d2b79f5) | 0; + let t = Math.imul(a ^ (a >>> 15), 1 | a); + t = (t + Math.imul(t ^ (t >>> 7), 61 | t)) ^ t; + return ((t ^ (t >>> 14)) >>> 0) / 4294967296; + }; + const n = () => { + let x = 0; + while (x === 0) x = u(); + return Math.sqrt(-2 * Math.log(x)) * Math.cos(2 * Math.PI * u()); + }; + return { u, n }; +} + +/** Stationary AR(1) noise with lag-1 autocorrelation `rho` and SD `sd`. */ +export function ar1(rng: Rng, length: number, rho: number, sd: number) { + const out: number[] = []; + let v = rng.n() * sd; + for (let i = 0; i < length; i++) { + v = rho * v + Math.sqrt(1 - rho * rho) * sd * rng.n(); + out.push(v); + } + return out; +} + +const MS_PER_DAY = 86_400_000; + +/** `YYYY-MM-DD` of `start` plus `offset` days. */ +export function dayKey(start: string, offset: number): string { + const [y, m, d] = start.split("-").map(Number); + return new Date(Date.UTC(y, m - 1, d) + offset * MS_PER_DAY) + .toISOString() + .slice(0, 10); +} + +/** + * Day-keyed points from an array starting at `start`; `null` entries are gaps + * (no point emitted). + */ +export function toSeries( + values: ReadonlyArray, + start = "2025-01-01", +): DailySeriesPoint[] { + const out: DailySeriesPoint[] = []; + values.forEach((value, i) => { + if (value != null) out.push({ day: dayKey(start, i), value }); + }); + return out; +} + +/** + * A behaviour/outcome pair with a real next-day coupling: the outcome on day + * D+1 is `slope` × the behaviour on day D plus independent noise. Both arms + * are mildly autocorrelated, like real daily vitals. + */ +export function coupledPair( + seed: number, + length: number, + slope: number, + opts: { noiseSd?: number; behaviourRho?: number; outcomeRho?: number } = {}, +): { behaviour: number[]; outcome: number[] } { + const rng = seededRng(seed); + const behaviour = ar1(rng, length, opts.behaviourRho ?? 0.2, 1); + const noise = ar1(rng, length, opts.outcomeRho ?? 0.2, opts.noiseSd ?? 1); + const outcome = noise.map( + (e, d) => 50 + e + (d > 0 ? slope * behaviour[d - 1] : 0), + ); + return { behaviour: behaviour.map((v) => 10 + v), outcome }; +} + +export interface SeasonalScenario { + /** Days in the scenario. */ + days: number; + /** + * Day of the year the scenario starts on (0 = 1 January). `"random"` draws + * a phase per run from the same seeded stream. + */ + startDoy: number | "random"; + /** + * True next-day effect of the temperature ANOMALY on systolic, mmHg per K. + * 0 is the null scenario: the two series share the season and nothing else. + */ + beta?: number; + /** Fraction of days dropped at random from each series. */ + missing?: number; + /** Length of one contiguous block missing from the systolic series. */ + block?: number; +} + +/** + * The simulation model behind the v1.42 statistics change (#615 research): + * + * temperature = 10 + 8·sin(2π·doy/365) + AR(1)(ρ 0.7, sd 3) + * systolic = 125 − 3·sin(2π·doy/365) + AR(1)(ρ 0.3, sd 8) + * + beta · temperatureAnomaly(day − 1) + * + * Temperature is warm in summer, blood pressure lower: with `beta = 0` the + * two share only the season, so any "finding" is a false positive. Returns + * day-keyed series in 2025 (+ following year as needed). + */ +export function seasonalScenario( + seed: number, + scenario: SeasonalScenario, +): { + temperature: DailySeriesPoint[]; + systolic: DailySeriesPoint[]; +} { + const rng = seededRng(seed); + const doy0 = + scenario.startDoy === "random" + ? Math.floor(rng.u() * 365) + : scenario.startDoy; + const { days } = scenario; + const beta = scenario.beta ?? 0; + const tempNoise = ar1(rng, days, 0.7, 3); + const sysNoise = ar1(rng, days, 0.3, 8); + const temp: Array = []; + const sys: Array = []; + for (let d = 0; d < days; d++) { + const phase = (2 * Math.PI * (doy0 + d)) / 365; + temp.push(10 + 8 * Math.sin(phase) + tempNoise[d]); + sys.push( + 125 - + 3 * Math.sin(phase) + + sysNoise[d] + + (d > 0 ? beta * tempNoise[d - 1] : 0), + ); + } + const missing = scenario.missing ?? 0; + const drop = (values: Array) => + values.map((v) => (rng.u() < missing ? null : v)); + const tempOut = drop(temp); + const sysOut = drop(sys); + const block = scenario.block ?? 0; + if (block > 0) { + const startAt = Math.floor(rng.u() * (days - block)); + for (let d = startAt; d < startAt + block; d++) sysOut[d] = null; + } + const start = dayKey("2025-01-01", doy0); + return { + temperature: toSeries(tempOut, start), + systolic: toSeries(sysOut, start), + }; +} diff --git a/src/lib/insights/__tests__/seasonal-adjust.test.ts b/src/lib/insights/__tests__/seasonal-adjust.test.ts new file mode 100644 index 000000000..f0ba59c00 --- /dev/null +++ b/src/lib/insights/__tests__/seasonal-adjust.test.ts @@ -0,0 +1,282 @@ +/** + * The seasonal/trend adjustment and the effective sample size, pinned at the + * level where the statistics live: one adjusted Pearson per pair. + * + * Every random draw is seeded (see `helpers/seeded-series.ts`); no test reads + * `Math.random`, so each rate below is the same number on every run. The + * calibration tests are written so that removing either half of the method + * (the adjustment or the effective n) turns them red; the "without n_eff" + * case proves the second half is load-bearing, not decoration. + */ +import { describe, expect, it } from "vitest"; + +import { + pearson, + seasonallyAdjustedPearson, + type DayPair, +} from "../correlations"; +import { lagJoin, type DailySeriesPoint } from "../correlation-discovery"; +import { + autocorrelation, + basisForSpan, + dayKeyToIndex, + effectiveSampleSize, + minFitRows, + residualise, + residualisePair, + SEASONAL_MIN_SPAN_DAYS, +} from "../seasonal-adjust"; +import { + ar1, + seasonalScenario, + seededRng, + type SeasonalScenario, +} from "./helpers/seeded-series"; + +/** Next-day pairs (temperature on D, systolic on D+1), keyed on D. */ +function lagPairs( + behaviour: DailySeriesPoint[], + outcome: DailySeriesPoint[], +): DayPair[] { + const byDay = new Map(outcome.map((p) => [p.day, p.value])); + const pairs: DayPair[] = []; + for (const b of behaviour) { + const next = dayKeyToIndex(b.day) + 1; + const key = new Date(next * 86_400_000).toISOString().slice(0, 10); + const y = byDay.get(key); + if (y !== undefined) pairs.push({ day: b.day, x: b.value, y }); + } + return pairs; +} + +const RUNS = 200; + +/** Share of seeded runs whose adjusted p-value falls below 0.05. */ +function rejectionRate( + scenario: SeasonalScenario, + test: (pairs: DayPair[]) => number | null, +): number { + let rejected = 0; + let tested = 0; + for (let seed = 1; seed <= RUNS; seed++) { + const { temperature, systolic } = seasonalScenario(1000 + seed, scenario); + const p = test(lagPairs(temperature, systolic)); + if (p === null) continue; + tested++; + if (p < 0.05) rejected++; + } + expect(tested).toBe(RUNS); + return rejected / tested; +} + +const adjustedP = (pairs: DayPair[]) => { + const r = seasonallyAdjustedPearson(pairs); + return r.status === "ok" ? r.pValue : null; +}; + +/** The method with its second half removed: residuals, but df = n − 2. */ +const adjustedWithoutEffectiveN = (pairs: DayPair[]) => { + const days = pairs.map((p) => dayKeyToIndex(p.day as string)); + const res = residualisePair( + days, + pairs.map((p) => p.x), + pairs.map((p) => p.y), + ); + if (!res) return null; + const r = pearson({ xs: res.xs, ys: res.ys }); + return r.status === "ok" ? r.pValue : null; +}; + +const rawP = (pairs: DayPair[]) => { + const r = pearson({ xs: pairs.map((p) => p.x), ys: pairs.map((p) => p.y) }); + return r.status === "ok" ? r.pValue : null; +}; + +describe("residualise", () => { + it("removes an exact annual wave plus trend to zero", () => { + const days = Array.from({ length: 200 }, (_, i) => 20_000 + i); + const values = days.map( + (d) => 3 + 0.01 * d + 5 * Math.sin((2 * Math.PI * d) / 365.25 + 0.4), + ); + const res = residualise(days, values, "harmonic"); + expect(res).not.toBeNull(); + for (const v of res!) expect(Math.abs(v)).toBeLessThan(1e-6); + }); + + it("removes only the line under the trend basis", () => { + const days = Array.from({ length: 60 }, (_, i) => 20_000 + i); + const values = days.map((d, i) => 2 * d + (i % 2 === 0 ? 1 : -1)); + const res = residualise(days, values, "trend")!; + // The alternating ±1 survives; the slope is gone. + expect(Math.abs(res[0] - 1)).toBeLessThan(0.05); + expect(Math.abs(res[1] + 1)).toBeLessThan(0.05); + }); + + it("refuses to fit fewer than basis size + 10 rows", () => { + const rows = minFitRows("harmonic") - 1; + const days = Array.from({ length: rows }, (_, i) => i); + expect(residualise(days, days.map(Number), "harmonic")).toBeNull(); + expect(minFitRows("trend")).toBe(12); + expect(minFitRows("harmonic")).toBe(14); + }); + + it("drops to the trend basis under 120 days of span", () => { + expect(SEASONAL_MIN_SPAN_DAYS).toBe(120); + expect(basisForSpan(119)).toBe("trend"); + expect(basisForSpan(120)).toBe("harmonic"); + }); +}); + +describe("effectiveSampleSize (Pyper–Peterman)", () => { + it("keeps about n for two white-noise series", () => { + const rng = seededRng(3); + const n = 300; + const days = Array.from({ length: n }, (_, i) => i); + const x = Array.from({ length: n }, () => rng.n()); + const y = Array.from({ length: n }, () => rng.n()); + const ne = effectiveSampleSize(days, x, y); + expect(ne).toBeGreaterThan(0.85 * n); + expect(ne).toBeLessThanOrEqual(n); + }); + + it("shrinks n well below 0.6·n for two persistent AR(0.7) series", () => { + const rng = seededRng(4); + const n = 300; + const days = Array.from({ length: n }, (_, i) => i); + const x = ar1(rng, n, 0.7, 1); + const y = ar1(rng, n, 0.7, 1); + expect(effectiveSampleSize(days, x, y)).toBeLessThan(0.6 * n); + }); + + it("reads autocorrelation over calendar lags, skipping gaps", () => { + // Days 0..9 and 20..29: lag-1 pairs exist only inside each block. + const days = [ + ...Array.from({ length: 10 }, (_, i) => i), + ...Array.from({ length: 10 }, (_, i) => 20 + i), + ]; + const values = days.map((d) => (d % 2 === 0 ? 1 : -1)); + expect(autocorrelation(days, values, 1)).toBeCloseTo(-1, 5); + expect(autocorrelation(days, values, 2)).toBeCloseTo(1, 5); + }); +}); + +describe("seasonallyAdjustedPearson: the seeded simulation", () => { + it("rejects the season, not the weather (seed 42, 400 days from 1 January)", () => { + const { temperature, systolic } = seasonalScenario(42, { + days: 400, + startDoy: 0, + }); + const pairs = lagPairs(temperature, systolic); + const raw = pearson({ + xs: pairs.map((p) => p.x), + ys: pairs.map((p) => p.y), + }); + const adjusted = seasonallyAdjustedPearson(pairs); + expect(raw.status === "ok" && raw.pValue).toBeLessThan(1e-6); + expect(adjusted.status).toBe("ok"); + if (adjusted.status !== "ok") return; + expect(adjusted.basis).toBe("harmonic"); + expect(adjusted.pValue).toBeGreaterThan(0.05); + expect(adjusted.effectiveN).toBeLessThan(adjusted.n); + }); + + it("holds the false-positive rate near 5 % over 200 seeds of shared season (400 days)", () => { + const scenario: SeasonalScenario = { days: 400, startDoy: 0 }; + expect(rejectionRate(scenario, rawP)).toBeGreaterThan(0.9); + expect(rejectionRate(scenario, adjustedP)).toBeLessThanOrEqual(0.06); + // Remove the effective n and the same residuals reject too often: the + // serial correlation of the residuals is half of the problem. + expect(rejectionRate(scenario, adjustedWithoutEffectiveN)).toBeGreaterThan( + 0.08, + ); + }); + + it("holds it on the engine's 180-day window at a random point in the year", () => { + const scenario: SeasonalScenario = { days: 180, startDoy: "random" }; + expect(rejectionRate(scenario, rawP)).toBeGreaterThan(0.3); + expect(rejectionRate(scenario, adjustedP)).toBeLessThanOrEqual(0.06); + }); + + it("holds it with 30 % missing days and a three-week gap", () => { + const scenario: SeasonalScenario = { + days: 400, + startDoy: 0, + missing: 0.3, + block: 21, + }; + expect(rejectionRate(scenario, adjustedP)).toBeLessThanOrEqual(0.06); + }); + + it("finds a real next-day effect, with the sign the season was hiding", () => { + const { temperature, systolic } = seasonalScenario(42, { + days: 400, + startDoy: 0, + beta: 0.7, + }); + const pairs = lagPairs(temperature, systolic); + const raw = pearson({ + xs: pairs.map((p) => p.x), + ys: pairs.map((p) => p.y), + }); + const adjusted = seasonallyAdjustedPearson(pairs); + // Raw: warmer days read as LOWER next-day pressure, the season's sign. + expect(raw.status === "ok" && raw.r).toBeLessThan(0); + expect(adjusted.status).toBe("ok"); + if (adjusted.status !== "ok") return; + expect(adjusted.r).toBeGreaterThan(0.15); + expect(adjusted.pValue).toBeLessThan(0.01); + }); + + it("has the power to find that effect in most runs", () => { + const power = rejectionRate( + { days: 400, startDoy: 0, beta: 0.7 }, + adjustedP, + ); + expect(power).toBeGreaterThan(0.85); + }); + + it("uses the trend-only path on a short series and stays calibrated", () => { + const scenario: SeasonalScenario = { days: 60, startDoy: 60 }; + const { temperature, systolic } = seasonalScenario(7, scenario); + const one = seasonallyAdjustedPearson(lagPairs(temperature, systolic)); + expect(one.status === "ok" && one.basis).toBe("trend"); + // 6.0 % over 1 000 seeds; 200 seeds sample it with about ±1.7 points. + expect(rejectionRate(scenario, adjustedP)).toBeLessThanOrEqual(0.08); + }); + + it("does not test a pair too short to fit, rather than testing it raw", () => { + const { temperature, systolic } = seasonalScenario(9, { + days: 12, + startDoy: 0, + }); + const pairs = lagPairs(temperature, systolic); + expect(pairs.length).toBe(11); + const result = seasonallyAdjustedPearson(pairs, { minPairs: 5 }); + expect(result).toEqual({ + status: "insufficient", + reason: "too_few_pairs", + n: 11, + }); + }); + + it("reports no variance for a series that is nothing but trend", () => { + const pairs: DayPair[] = Array.from({ length: 40 }, (_, i) => ({ + day: 20_000 + i, + x: 3 + i, + y: Math.sin(i), + })); + expect(seasonallyAdjustedPearson(pairs).status).toBe("insufficient"); + }); + + it("agrees with the engine's lag join on which days pair", () => { + const { temperature, systolic } = seasonalScenario(5, { + days: 60, + startDoy: 0, + missing: 0.2, + }); + const joined = lagJoin(temperature, systolic, 1); + const pairs = lagPairs(temperature, systolic); + expect(pairs.map((p) => p.x)).toEqual(joined.xs); + expect(pairs.map((p) => p.y)).toEqual(joined.ys); + }); +}); diff --git a/src/lib/insights/correlation-channel-series.ts b/src/lib/insights/correlation-channel-series.ts index ecbca30b3..76b0c9c22 100644 --- a/src/lib/insights/correlation-channel-series.ts +++ b/src/lib/insights/correlation-channel-series.ts @@ -54,7 +54,7 @@ import { pickMainNightAndNaps, type SleepStageRow, } from "@/lib/analytics/sleep-night"; -import { isNearUtc } from "@/lib/tz/format"; +import { isNearUtc, shiftDateKey } from "@/lib/tz/format"; import { probeRollupCoverage } from "@/lib/rollups/measurement-coverage"; import { loadUserSourcePriority } from "@/lib/rollups/measurement-read"; import type { @@ -258,7 +258,7 @@ export async function fetchEnvironmentSeries( }); return ENVIRONMENT_FIELDS.map((field) => { - const points: DailySeriesPoint[] = []; + const daily: DailySeriesPoint[] = []; for (const row of rows) { const raw = row[field.column]; if (raw == null || !Number.isFinite(raw)) continue; @@ -267,12 +267,42 @@ export async function fetchEnvironmentSeries( field.column === "sunshineSec" || field.column === "daylightSec" ? raw / 3600 : raw; - points.push({ day: row.date, value }); + daily.push({ day: row.date, value }); } - return { key: field.key, role: "behaviour" as const, points }; + // v1.42 — one averaged-lag hypothesis per outcome: the exposure on day D + // is the mean of D−1 and D, paired with the outcome on D (`lagDays: 0`). + // See the lag note in `src/lib/environment/fields.ts`. + return { + key: field.key, + role: "behaviour" as const, + points: averageWithPreviousDay(daily), + lagDays: 0, + }; }); } +/** + * Exposure over "the day before and the day itself": for each day D whose + * previous calendar day is also present, the mean of the two values, keyed + * on D. A day without its predecessor yields no point (no half-window mean). + * Pure; input order is irrelevant, output is ascending by day. + */ +export function averageWithPreviousDay( + points: readonly DailySeriesPoint[], +): DailySeriesPoint[] { + const byDay = new Map(); + for (const p of points) { + if (Number.isFinite(p.value)) byDay.set(p.day, p.value); + } + const out: DailySeriesPoint[] = []; + for (const day of [...byDay.keys()].sort()) { + const prev = byDay.get(shiftDateKey(day, -1)); + if (prev === undefined) continue; + out.push({ day, value: (prev + (byDay.get(day) as number)) / 2 }); + } + return out; +} + /** Explicit deterministic ceiling on opt-in custom behaviour channels. */ export const MAX_CUSTOM_CORRELATION_CHANNELS = 8; export const CUSTOM_METRIC_CHANNEL_PREFIX = "CUSTOM_METRIC:"; diff --git a/src/lib/insights/correlation-discovery.ts b/src/lib/insights/correlation-discovery.ts index 0fd87e959..bc8edd017 100644 --- a/src/lib/insights/correlation-discovery.ts +++ b/src/lib/insights/correlation-discovery.ts @@ -37,8 +37,21 @@ * n ≥ 20 / p < 0.05 / BH-FDR / effect-size-floor / shrinkage gates as every * other channel, so a thin series degrades to absent rather than to a spurious * link. The series builders live in `correlation-series-builders.ts`. + * + * v1.42 — every daily pair is tested on seasonally adjusted residuals: the + * annual cycle (from 120 paired days on) and a linear trend are removed from + * both series, and the p-value uses the Pyper–Peterman effective sample size + * instead of n. Two series that only share the calendar no longer surface as + * related (see `seasonal-adjust.ts` for the method and the simulation behind + * it). Shrinkage and the phrasing tier read the effective size too, so a + * strongly autocorrelated pair cannot be narrated as a deep, confident one. */ -import { pearson, MIN_PAIRED_N } from "@/lib/insights/correlations"; +import { + pearson, + seasonallyAdjustedPearson, + MIN_PAIRED_N, + type DayPair, +} from "@/lib/insights/correlations"; import { type Locale } from "@/lib/i18n/config"; import { getServerTranslator } from "@/lib/i18n/server-translator"; import { shiftDateKey } from "@/lib/tz/format"; @@ -173,6 +186,14 @@ export interface NamedSeries { /** User-facing label for dynamic channels such as custom metrics. */ label?: string; points: DailySeriesPoint[]; + /** + * Behaviour channels only: the lag this channel pairs at, overriding the + * scan's default (1). The environmental-exposure channels set 0 because + * their points are already the mean of the day before and the day itself + * (see the environment read in `correlation-channel-series.ts`), so "lag 0" + * there means the averaged 0–1 window, tested as one hypothesis. + */ + lagDays?: number; } /** One discovered, FDR-surviving correlation pair. */ @@ -180,12 +201,20 @@ export interface DiscoveredCorrelation { behaviour: string; outcome: string; behaviourLabel?: string; - /** Paired-day count after the day+1 lag join. */ outcomeLabel?: string; + /** Paired-day count after the lag join. */ n: number; - /** Pearson r (lag-joined). */ + /** + * Effective sample size after serial correlation (Pyper–Peterman). The + * p-value, the shrinkage and the tier use it; `n` stays the day count. + */ + nEff?: number; + /** + * Partial correlation of the lag-joined pair with trend (and, from 120 + * paired days, the annual cycle) removed from both series. + */ r: number; - /** Two-sided exact Student-t p-value. */ + /** Two-sided exact Student-t p-value at `nEff` degrees of freedom. */ pValue: number; /** Benjamini-Hochberg adjusted q-value. */ qValue: number; @@ -202,7 +231,10 @@ export interface DiscoveredCorrelation { tier: ConfidenceTier; /** Conservative, descriptive interpretation — never causal, tier-hedged. */ interpretation: string; - /** Lag in days applied (always 1 here). */ + /** + * Lag in days applied: 1 (behaviour day D, outcome D+1) by default; 0 for + * the environmental channels, whose exposure is the mean of D−1 and D. + */ lagDays: number; /** * v1.22 — which window surfaced the pair. `retrospective` = the standard @@ -264,12 +296,31 @@ export function lagJoin( return { xs, ys }; } +/** {@link lagJoin} keeping the behaviour day each pair is keyed on. */ +function lagJoinPairs( + behaviour: DailySeriesPoint[], + outcome: DailySeriesPoint[], + lagDays: number, +): DayPair[] { + const outcomeByDay = new Map(outcome.map((p) => [p.day, p.value])); + const pairs: DayPair[] = []; + for (const b of behaviour) { + const target = outcomeByDay.get(shiftDay(b.day, lagDays)); + if (target != null && Number.isFinite(b.value) && Number.isFinite(target)) { + pairs.push({ day: b.day, x: b.value, y: target }); + } + } + return pairs; +} + interface RawPair { behaviour: string; outcome: string; n: number; + nEff: number; r: number; pValue: number; + lagDays: number; behaviourLabel?: string; outcomeLabel?: string; } @@ -317,6 +368,7 @@ function interpret( outcome: string, r: number, tier: ConfidenceTier, + lagDays: number, t: Translate, behaviourLabel?: string, outcomeLabel?: string, @@ -333,7 +385,10 @@ function interpret( // "on the evidence so far" phrasing. Each tier keeps its own template. const tierKey = tier === "faint" ? "faint" : tier === "high" ? "high" : "moderate"; - return t(`insights.correlation.daily.${tierKey}`, params); + // Lag 0 is the averaged same-day window of the exposure channels ("on the + // day and the day before"); every other pair reads as next-day. + const template = lagDays === 0 ? "dailySameDay" : "daily"; + return t(`insights.correlation.${template}.${tierKey}`, params); } /** @@ -430,9 +485,11 @@ export function discoveryMeasurementTypes(keys: readonly string[]): string[] { /** * Run the FDR-controlled discovery over the behaviour × outcome matrix. * - * 1. For every (behaviour, outcome) pair, lag-join (D → D+1) and run - * Pearson with the exact p-value. Pairs below `minPairs` paired days - * are dropped (not tested — they cannot be defensibly assessed). + * 1. For every (behaviour, outcome) pair, lag-join (D → D+lag, the series' + * own lag or the scan default), remove trend and season from both arms + * and run Pearson on the residuals with the effective-n p-value. Pairs + * below `minPairs` paired days, or too few to fit the basis, are dropped + * (not tested — they cannot be defensibly assessed). * 2. Across the tested pairs, compute BH q-values. * 3. Surface pairs with p < 0.05 AND q ≤ `fdrQ`, ranked by q then |r|. * @@ -476,16 +533,19 @@ export function discoverCorrelations( // auto-correlation tautology, not a cross-domain insight; excluding it // before testing stops it crowding out genuine cross-metric links. if (metricFamily(b.key) === metricFamily(o.key)) continue; - const { xs, ys } = lagJoin(b.points, o.points, lagDays); - if (xs.length < minPairs) continue; - const result = pearson({ xs, ys, minPairs }); + const pairLag = b.lagDays ?? lagDays; + const pairs = lagJoinPairs(b.points, o.points, pairLag); + if (pairs.length < minPairs) continue; + const result = seasonallyAdjustedPearson(pairs, { minPairs }); if (result.status !== "ok") continue; tested.push({ behaviour: b.key, outcome: o.key, n: result.n, + nEff: result.effectiveN, r: result.r, pValue: result.pValue, + lagDays: pairLag, behaviourLabel: b.label, outcomeLabel: o.label, }); @@ -506,8 +566,10 @@ export function discoverCorrelations( // tiers to `null` and is dropped: it is statistically real but explains too // little variance to narrate as a driver. .map((t) => { - const shrunkR = shrinkEstimate(t.r, t.n); - const tier = confidenceTier(shrunkR, t.n); + // v1.42 — shrink and tier by the EFFECTIVE size: a persistent pair's + // 180 days may be worth 30 independent ones, and is narrated as such. + const shrunkR = shrinkEstimate(t.r, t.nEff); + const tier = confidenceTier(shrunkR, t.nEff); return { ...t, shrunkR, tier }; }) .filter( @@ -525,6 +587,7 @@ export function discoverCorrelations( behaviourLabel: t.behaviourLabel, outcomeLabel: t.outcomeLabel, n: t.n, + nEff: t.nEff, r: t.r, pValue: t.pValue, qValue: Math.round(t.qValue * 1000) / 1000, @@ -535,11 +598,12 @@ export function discoverCorrelations( t.outcome, t.r, t.tier, + t.lagDays, translate, t.behaviourLabel, t.outcomeLabel, ), - lagDays, + lagDays: t.lagDays, })) // RECON1 (D2-2 / D4) — rank by the SHRUNK effect magnitude (a deep, strong // pair leads a thin one even if the thin one's raw r is higher), with q as diff --git a/src/lib/insights/correlations.ts b/src/lib/insights/correlations.ts index 9ceffa0c6..bedb700a6 100644 --- a/src/lib/insights/correlations.ts +++ b/src/lib/insights/correlations.ts @@ -41,6 +41,11 @@ */ import { defaultLocale } from "@/lib/i18n/config"; import { getServerTranslator } from "@/lib/i18n/server-translator"; +import { + dayKeyToIndex, + residualisePair, + type SeasonalBasis, +} from "@/lib/insights/seasonal-adjust"; /** Locale-aware string builder — mirrors the server translator's `t`. */ type Translate = ( @@ -149,6 +154,12 @@ export interface PearsonOk { pValue: number; /** 95 % confidence interval on r via Fisher z-transform. */ confidenceInterval: [number, number]; + /** + * The effective sample size the p-value and interval used, when the caller + * supplied one (serially correlated input). Absent on the plain path, where + * the p-value uses `n`. + */ + effectiveN?: number; } export interface PearsonInsufficient { @@ -168,6 +179,21 @@ export function pearson(input: { xs: readonly number[]; ys: readonly number[]; minPairs?: number; + /** + * Effective sample size for serially correlated input (see + * `effectiveSampleSize` in `seasonal-adjust.ts`). When given, the t-test + * uses `min(effectiveN, n − covariates) − 2` degrees of freedom and the + * Fisher-z interval `effectiveN − 3`, so a persistent series cannot claim the + * evidence of `n` independent days. `minPairs` still gates on the observed + * pair count. + */ + effectiveN?: number; + /** + * Nuisance columns already partialled out of both arms (trend, harmonic). + * Each costs one degree of freedom, as in any partial correlation. Only + * read when `effectiveN` is given. + */ + covariates?: number; }): PearsonResult { const { xs, ys } = input; const minPairs = input.minPairs ?? MIN_PAIRED_N; @@ -210,15 +236,24 @@ export function pearson(input: { // grid past the three pre-defined cards here; the exact p-value // pre-empts a class of false positives that the normal-approx // overstated at low df. + const effectiveN = + input.effectiveN !== undefined && Number.isFinite(input.effectiveN) + ? Math.min(n, Math.max(3, input.effectiveN)) + : undefined; + const df = + effectiveN === undefined + ? n - 2 + : Math.min(effectiveN, n - (input.covariates ?? 0)) - 2; const tStat = Math.abs(clamped) >= 1 ? Number.POSITIVE_INFINITY - : (clamped * Math.sqrt(n - 2)) / Math.sqrt(1 - clamped * clamped); - const pValue = twoSidedPFromT(Math.abs(tStat), n - 2); + : (clamped * Math.sqrt(Math.max(0, df))) / + Math.sqrt(1 - clamped * clamped); + const pValue = twoSidedPFromT(Math.abs(tStat), df); // Fisher z-transform for 95 % CI on r. const z = atanh(clamped); - const se = 1 / Math.sqrt(Math.max(1, n - 3)); + const se = 1 / Math.sqrt(Math.max(1, (effectiveN ?? n) - 3)); const zLow = z - 1.96 * se; const zHigh = z + 1.96 * se; const confidenceInterval: [number, number] = [tanh(zLow), tanh(zHigh)]; @@ -229,9 +264,97 @@ export function pearson(input: { n, pValue, confidenceInterval, + ...(effectiveN === undefined + ? {} + : { effectiveN: Math.round(effectiveN * 10) / 10 }), }; } +/** One paired observation keyed on a calendar day. */ +export interface DayPair { + /** Day index (days since 1970-01-01) or a `YYYY-MM-DD` key. */ + day: number | string; + x: number; + y: number; +} + +export type AdjustedPearsonResult = + | (PearsonOk & { effectiveN: number; basis: SeasonalBasis }) + | PearsonInsufficient; + +/** + * Pearson for two daily series with the shared season and trend removed and + * the p-value computed from the effective sample size. The discovery engine + * (every surface that scans the matrix) and the dated fixed-hypothesis cards + * below go through it. + * + * The reported `r` is the correlation of the residuals, i.e. the partial + * correlation with trend (and, over 120 days or more, the annual cycle) held + * fixed. `n` stays the observed pair count; `effectiveN` is what the evidence + * is worth after serial correlation. Pairs too few to fit the basis come back + * `too_few_pairs` and are not tested at all. + */ +export function seasonallyAdjustedPearson( + pairs: readonly DayPair[], + opts: { minPairs?: number } = {}, +): AdjustedPearsonResult { + const minPairs = opts.minPairs ?? MIN_PAIRED_N; + const n = pairs.length; + if (n < minPairs) + return { status: "insufficient", reason: "too_few_pairs", n }; + const days = pairs.map((p) => + typeof p.day === "number" ? p.day : dayKeyToIndex(p.day), + ); + const adjusted = residualisePair( + days, + pairs.map((p) => p.x), + pairs.map((p) => p.y), + ); + if (!adjusted) return { status: "insufficient", reason: "too_few_pairs", n }; + // A series that is nothing but trend (and season) leaves residuals that are + // zero up to rounding. Pearson on those would correlate floating-point + // noise, so treat a residual spread below 1e-9 of the raw spread as none. + if ( + negligibleSpread( + pairs.map((p) => p.x), + adjusted.xs, + ) || + negligibleSpread( + pairs.map((p) => p.y), + adjusted.ys, + ) + ) { + return { status: "insufficient", reason: "no_variance", n }; + } + const result = pearson({ + xs: adjusted.xs, + ys: adjusted.ys, + minPairs, + effectiveN: adjusted.effectiveN, + covariates: adjusted.covariates, + }); + if (result.status !== "ok") return result; + return { + ...result, + effectiveN: result.effectiveN ?? n, + basis: adjusted.basis, + }; +} + +function sumSquaresAboutMean(values: readonly number[]): number { + const mean = values.reduce((s, v) => s + v, 0) / values.length; + return values.reduce((s, v) => s + (v - mean) ** 2, 0); +} + +function negligibleSpread( + raw: readonly number[], + residuals: readonly number[], +): boolean { + const rawSs = sumSquaresAboutMean(raw); + if (rawSs === 0) return true; + return sumSquaresAboutMean(residuals) <= rawSs * 1e-9; +} + function atanh(x: number): number { // Math.atanh in case the runtime polyfill misbehaves; both should agree. return 0.5 * Math.log((1 + x) / (1 - x)); @@ -631,6 +754,15 @@ function twoSidedPFromF( return 1 - normalCdf(z); } +/** + * Day index of a day-start `Date` (the fast path builds these from user-zone + * day keys at UTC midnight). Rounded so a DST-shifted midnight lands on the + * same day. + */ +function dayIndexOfDate(date: Date): number { + return Math.round(date.getTime() / 86_400_000); +} + // ── Hypothesis #1 — BP × medication compliance ───────────────── export interface BpComplianceInput { @@ -657,9 +789,13 @@ export function correlateBpCompliance( x: d.compliancePct, y: d.systolic, })); - const xs = points.map((p) => p.x); - const ys = points.map((p) => p.y); - const r = pearson({ xs, ys }); + const r = seasonallyAdjustedPearson( + input.daily.map((d, i) => ({ + day: dayIndexOfDate(d.date), + x: points[i].x, + y: points[i].y, + })), + ); if (r.status === "insufficient") { return { kind: "bp-compliance", @@ -718,9 +854,13 @@ export function correlateMoodPulse( t: Translate = defaultTranslate(), ): CorrelationResult { const points = input.daily.map((d) => ({ x: d.mood, y: d.restingPulse })); - const xs = points.map((p) => p.x); - const ys = points.map((p) => p.y); - const r = pearson({ xs, ys }); + const r = seasonallyAdjustedPearson( + input.daily.map((d, i) => ({ + day: dayIndexOfDate(d.date), + x: points[i].x, + y: points[i].y, + })), + ); if (r.status === "insufficient") { return { kind: "mood-pulse", diff --git a/src/lib/insights/discovery-matrix.ts b/src/lib/insights/discovery-matrix.ts index 56a207a9c..947684c16 100644 --- a/src/lib/insights/discovery-matrix.ts +++ b/src/lib/insights/discovery-matrix.ts @@ -58,6 +58,7 @@ import { type MoodFactorWindowFetch, } from "@/lib/insights/correlation-channel-series"; import { loadUserSourcePriority } from "@/lib/rollups/measurement-read"; +import { ENVIRONMENT_FIELDS } from "@/lib/environment/fields"; import { correlationChannelSurfaceId, isSurfaceVisible, @@ -303,7 +304,16 @@ export async function assembleDiscoveryMatrix( readMood(userId, tz, since, opts.includeMoodFactors === true), fetchComplianceSeries(userId, tz, since), fetchSymptomSeries(userId, tz, since), - fetchEnvironmentSeries(userId, since), + // The weather read is skipped outright when the environment module is off: + // the mask below would drop every ENV_* series anyway, so reading them + // was a wasted query on every scan. Probed through the first field's + // surface id, which the surface map owns (all ENV_* share one module). + isSurfaceVisible( + correlationChannelSurfaceId(ENVIRONMENT_FIELDS[0].key), + opts.modules, + ) + ? fetchEnvironmentSeries(userId, since) + : Promise.resolve([] as NamedSeries[]), fetchCustomMetricBehaviourSeries(userId, tz, since), fetchSymptomEventSeries(userId, tz, since), ]); @@ -332,9 +342,11 @@ export async function assembleDiscoveryMatrix( series.push({ key, role: "outcome", points: points(key) }); } } - // Environmental exposure and custom metrics are lag SOURCES only: they pair - // (D → D+1) against every outcome above, under the unchanged n ≥ 20 / FDR / - // effect-size gates, so a thin weather or custom series degrades to absent. + // Environmental exposure and custom metrics are lag SOURCES only. Custom + // metrics pair D → D+1; the environmental series carry `lagDays: 0` over an + // exposure already averaged across D−1 and D (see `fetchEnvironmentSeries`). + // Both run under the same n ≥ 20 / FDR / effect-size gates, so a thin + // weather or custom series degrades to absent. series.push(...environmentSeries, ...customMetricSeries); // The person's own symptoms are OUTCOMES: "short sleep, then an aura the // next day" is the question they log them to answer. Each is its own diff --git a/src/lib/insights/narrative/__tests__/period-narrative.test.ts b/src/lib/insights/narrative/__tests__/period-narrative.test.ts index 83abb0ed5..c3f393266 100644 --- a/src/lib/insights/narrative/__tests__/period-narrative.test.ts +++ b/src/lib/insights/narrative/__tests__/period-narrative.test.ts @@ -12,6 +12,7 @@ import type { import { DEFAULT_UNIT_PREFERENCES } from "@/lib/measurements/display-transform"; import { buildDeterministicNarrative } from "../period-narrative-deterministic"; import { buildNarrativeUserPrompt } from "../period-narrative-generate"; +import { coupledPair } from "../../__tests__/helpers/seeded-series"; /** * Build a contiguous daily series ending at `endDay` (YYYY-MM-DD), one point @@ -310,12 +311,12 @@ describe("assemblePeriodNarrativeContext — band transitions", () => { describe("assemblePeriodNarrativeContext — drivers (FDR-controlled)", () => { it("surfaces a strong lagged relationship and keeps it descriptive-only", () => { - // 31 contiguous days; outcome[D+1] = 2 * behaviour[D] + small noise so the - // lag-1 join yields a strong Pearson r over ≥ 20 paired days. - const bVals = Array.from({ length: 31 }, (_, i) => 10 + i); - const oVals = Array.from({ length: 31 }, (_, i) => - i === 0 ? 0 : 2 * (10 + (i - 1)) + (i % 2 === 0 ? 0.5 : -0.5), - ); + // 31 contiguous days; outcome[D+1] = 2 * behaviour[D] + noise (seeded), + // so the lag-1 join yields a strong day-to-day r over ≥ 20 paired days. + // Not two straight lines: a shared trend is removed before correlating. + const { behaviour: bVals, outcome: oVals } = coupledPair(3, 31, 2, { + noiseSd: 0.5, + }); const behaviour = seriesEndingAt(bVals, "2026-04-30"); const outcome = seriesEndingAt(oVals, "2026-04-30"); const discoverySeries: NamedSeries[] = [ diff --git a/src/lib/insights/seasonal-adjust.ts b/src/lib/insights/seasonal-adjust.ts new file mode 100644 index 000000000..4db5e1c67 --- /dev/null +++ b/src/lib/insights/seasonal-adjust.ts @@ -0,0 +1,315 @@ +/** + * Seasonal and trend adjustment for the correlation engine. + * + * Two daily health series that both follow the calendar (blood pressure is + * lower in summer, temperature is higher, sleep and HRV drift with daylight) + * correlate strongly even when neither has anything to do with the other. A + * raw Pearson over a year of such data reports the season, not a relation: + * in a seeded simulation of two series that share only an annual cycle, the + * raw lag-1 test rejected at p < 0.05 in 97 % of runs. + * + * The fix has two halves, and both are needed: + * + * 1. Remove the shared slow component from BOTH series before correlating. + * Each series is regressed on `[1, t, sin(2πd/365.25), cos(2πd/365.25)]` + * (one annual harmonic plus a linear trend) and the residuals are + * correlated. By Frisch–Waugh–Lovell this is the partial correlation with + * season and trend held fixed. Over a short span the annual wave is not + * identifiable and is practically a straight line, so below + * {@link SEASONAL_MIN_SPAN_DAYS} only `[1, t]` is removed. + * 2. Count the information honestly. Daily residuals are still serially + * correlated (weather anomalies persist for days), so `n − 2` degrees of + * freedom overstate the evidence. The p-value uses the effective sample + * size of Pyper & Peterman (1998), which shrinks n by the product of the + * two residual series' autocorrelations. Without it the adjusted test + * still rejects 10–12 % of null pairs instead of 5 %. + * + * Pure: no imports, no clock, no I/O. Days are integer day indices (days since + * 1970-01-01, see {@link dayKeyToIndex}), so a gap is simply an absent index; + * the regression runs over the days that exist and nothing is imputed. + */ + +/** Span (first to last paired day, inclusive) below which only the trend is removed. */ +export const SEASONAL_MIN_SPAN_DAYS = 120; + +/** Length of the annual cycle the harmonic models. */ +export const ANNUAL_PERIOD_DAYS = 365.25; + +/** + * Rows required beyond the basis size to fit a series. Below `p + 10` points + * a series is left unfitted and its pair is not tested at all, rather than + * tested raw. + */ +export const FIT_MARGIN = 10; + +/** Smallest effective sample size the p-value may assume. */ +export const MIN_EFFECTIVE_N = 3; + +/** Which nuisance basis a fit removed. */ +export type SeasonalBasis = "harmonic" | "trend"; + +const MS_PER_DAY = 86_400_000; + +/** `YYYY-MM-DD` → whole days since 1970-01-01 (calendar arithmetic, no zone). */ +export function dayKeyToIndex(day: string): number { + const y = Number(day.slice(0, 4)); + const m = Number(day.slice(5, 7)); + const d = Number(day.slice(8, 10)); + return Math.round(Date.UTC(y, m - 1, d) / MS_PER_DAY); +} + +/** The basis a span of `spanDays` calendar days gets. */ +export function basisForSpan(spanDays: number): SeasonalBasis { + return spanDays < SEASONAL_MIN_SPAN_DAYS ? "trend" : "harmonic"; +} + +/** Number of columns in a basis (intercept included). */ +export function basisSize(basis: SeasonalBasis): number { + return basis === "harmonic" ? 4 : 2; +} + +/** Fewest rows a fit under `basis` accepts. */ +export function minFitRows(basis: SeasonalBasis): number { + return basisSize(basis) + FIT_MARGIN; +} + +/** + * Solve the small symmetric system `A·x = b` by Gauss–Jordan with partial + * pivoting. Returns null when the system is singular (e.g. every day equal). + */ +function solve(A: number[][], b: number[]): number[] | null { + const n = A.length; + const M = A.map((row, i) => [...row, b[i]]); + for (let c = 0; c < n; c++) { + let p = c; + for (let r = c + 1; r < n; r++) { + if (Math.abs(M[r][c]) > Math.abs(M[p][c])) p = r; + } + [M[c], M[p]] = [M[p], M[c]]; + const pivot = M[c][c]; + if (Math.abs(pivot) < 1e-10) return null; + for (let j = c; j <= n; j++) M[c][j] /= pivot; + for (let r = 0; r < n; r++) { + if (r === c) continue; + const f = M[r][c]; + if (f === 0) continue; + for (let j = c; j <= n; j++) M[r][j] -= f * M[c][j]; + } + } + return M.map((row) => row[n]); +} + +/** + * Residuals of `values` after an ordinary least-squares fit of the chosen + * basis over the given days. The trend column is centred and scaled to the + * span so the normal equations stay well conditioned. Returns null when there + * are fewer than {@link minFitRows} rows or the fit is singular. + */ +export function residualise( + days: readonly number[], + values: readonly number[], + basis: SeasonalBasis, +): number[] | null { + const n = days.length; + if (values.length !== n) { + throw new Error("residualise: days and values must have equal length"); + } + if (n < minFitRows(basis)) return null; + + let lo = Infinity; + let hi = -Infinity; + for (const d of days) { + if (d < lo) lo = d; + if (d > hi) hi = d; + } + const mid = (lo + hi) / 2; + const scale = Math.max(1, hi - lo); + const w = (2 * Math.PI) / ANNUAL_PERIOD_DAYS; + const p = basisSize(basis); + + const rows: number[][] = new Array(n); + for (let i = 0; i < n; i++) { + const d = days[i]; + const t = (d - mid) / scale; + rows[i] = + basis === "harmonic" ? [1, t, Math.sin(w * d), Math.cos(w * d)] : [1, t]; + } + + const XtX = Array.from({ length: p }, () => new Array(p).fill(0)); + const Xty = new Array(p).fill(0); + for (let i = 0; i < n; i++) { + const f = rows[i]; + const y = values[i]; + for (let a = 0; a < p; a++) { + Xty[a] += f[a] * y; + for (let b = a; b < p; b++) XtX[a][b] += f[a] * f[b]; + } + } + for (let a = 0; a < p; a++) { + for (let b = 0; b < a; b++) XtX[a][b] = XtX[b][a]; + } + + const beta = solve(XtX, Xty); + if (!beta) return null; + + const out = new Array(n); + for (let i = 0; i < n; i++) { + let fit = 0; + for (let k = 0; k < p; k++) fit += rows[i][k] * beta[k]; + out[i] = values[i] - fit; + } + return out; +} + +/** + * Autocorrelation of a gappy daily series at a CALENDAR lag of `k` days: the + * covariance over the pairs of days `k` apart that both exist, divided by the + * variance over every day that exists. Returns 0 when fewer than three such + * pairs exist or the series is constant. + * + * `dense` holds the series on a contiguous day grid starting at the first day, + * with `NaN` for a missing day. + */ +function autocorrelationDense( + dense: Float64Array, + mean: number, + variance: number, + k: number, +): number { + if (variance === 0) return 0; + let c = 0; + let cnt = 0; + for (let i = 0; i + k < dense.length; i++) { + const a = dense[i]; + const b = dense[i + k]; + if (Number.isNaN(a) || Number.isNaN(b)) continue; + c += (a - mean) * (b - mean); + cnt++; + } + return cnt < 3 ? 0 : c / cnt / variance; +} + +function densify( + days: readonly number[], + values: readonly number[], +): { dense: Float64Array; mean: number; variance: number } { + let lo = Infinity; + let hi = -Infinity; + for (const d of days) { + if (d < lo) lo = d; + if (d > hi) hi = d; + } + const dense = new Float64Array(hi - lo + 1).fill(Number.NaN); + let sum = 0; + for (let i = 0; i < days.length; i++) { + dense[days[i] - lo] = values[i]; + sum += values[i]; + } + const mean = sum / values.length; + let ss = 0; + for (const v of values) ss += (v - mean) ** 2; + return { dense, mean, variance: ss / values.length }; +} + +/** + * Autocorrelation of a gappy daily series at calendar lag `k`. Exposed for the + * tests; the engine goes through {@link effectiveSampleSize}. + */ +export function autocorrelation( + days: readonly number[], + values: readonly number[], + k: number, +): number { + if (values.length < 3) return 0; + const { dense, mean, variance } = densify(days, values); + return autocorrelationDense(dense, mean, variance, k); +} + +/** + * Pyper & Peterman (1998) effective sample size for the correlation of two + * serially correlated series observed on the same `days`: + * + * 1/n_eff = 1/n + (2/n) · Σ_{k=1}^{⌊n/5⌋} ((n − k)/n) · ρx(k) · ρy(k) + * + * with ρ the calendar-lag autocorrelations of each residual series. Clamped to + * `[MIN_EFFECTIVE_N, n]`: two persistent series lose most of their apparent + * information, two white-noise series keep about all of it, and opposite + * autocorrelation structures never inflate n beyond the observed pairs. + */ +export function effectiveSampleSize( + days: readonly number[], + xs: readonly number[], + ys: readonly number[], +): number { + const n = days.length; + if (xs.length !== n || ys.length !== n) { + throw new Error("effectiveSampleSize: inputs must have equal length"); + } + if (n < 3) return n; + const x = densify(days, xs); + const y = densify(days, ys); + const maxLag = Math.max(1, Math.floor(n / 5)); + let s = 0; + for (let k = 1; k <= maxLag; k++) { + const rx = autocorrelationDense(x.dense, x.mean, x.variance, k); + if (rx === 0) continue; + const ry = autocorrelationDense(y.dense, y.mean, y.variance, k); + s += ((n - k) / n) * rx * ry; + } + const inv = 1 / n + (2 / n) * s; + if (!(inv > 0)) return n; + return Math.min(n, Math.max(MIN_EFFECTIVE_N, 1 / inv)); +} + +/** Both residual series of one paired sample, ready for Pearson. */ +export interface SeasonalResiduals { + /** Which nuisance basis was removed (decided by the paired span). */ + basis: SeasonalBasis; + /** Columns removed besides the intercept (partial-correlation df cost). */ + covariates: number; + xs: number[]; + ys: number[]; + /** Pyper–Peterman effective sample size of the residual pair. */ + effectiveN: number; +} + +/** + * Residualise both arms of a paired sample on the same nuisance basis and + * compute the effective sample size of the result. + * + * `days` is the day index each pair is keyed on (the behaviour day for a + * lagged pair). Keying the outcome arm on the same index is exact: a constant + * lag shifts the trend into the intercept and rotates the sine/cosine pair + * into itself, so the column space is identical. + * + * Returns null when the sample is too small to fit the chosen basis; the + * caller then leaves the pair untested rather than testing it raw. + */ +export function residualisePair( + days: readonly number[], + xs: readonly number[], + ys: readonly number[], +): SeasonalResiduals | null { + const n = days.length; + if (xs.length !== n || ys.length !== n) { + throw new Error("residualisePair: inputs must have equal length"); + } + if (n === 0) return null; + let lo = Infinity; + let hi = -Infinity; + for (const d of days) { + if (d < lo) lo = d; + if (d > hi) hi = d; + } + const basis = basisForSpan(hi - lo + 1); + const rx = residualise(days, xs, basis); + const ry = residualise(days, ys, basis); + if (!rx || !ry) return null; + return { + basis, + covariates: basisSize(basis) - 1, + xs: rx, + ys: ry, + effectiveN: effectiveSampleSize(days, rx, ry), + }; +} diff --git a/src/lib/openapi/routes/cycle.ts b/src/lib/openapi/routes/cycle.ts index 2410b7f6d..e7580035c 100644 --- a/src/lib/openapi/routes/cycle.ts +++ b/src/lib/openapi/routes/cycle.ts @@ -376,6 +376,12 @@ const cyclePhaseLaggedPair = z behaviour: z.string(), outcome: z.string(), n: z.number().int(), + nEff: z + .number() + .optional() + .describe( + "v1.42 — effective sample size after serial correlation, ≤ `n`; the p-value uses it.", + ), r: z.number(), pValue: z.number(), qValue: z.number(), diff --git a/src/lib/openapi/routes/insights/schemas.ts b/src/lib/openapi/routes/insights/schemas.ts index 97c8b955e..3ba65922d 100644 --- a/src/lib/openapi/routes/insights/schemas.ts +++ b/src/lib/openapi/routes/insights/schemas.ts @@ -473,12 +473,23 @@ export const discoveredCorrelation = z .string() .optional() .describe("Display label for a dynamic outcome channel."), - n: z + n: z.number().int().describe("Paired-day count after the lag join (≥ 20)."), + nEff: z .number() - .int() - .describe("Paired-day count after the day+1 lag join (≥ 20)."), - r: z.number().describe("Pearson r over the lag-joined daily series."), - pValue: z.number().describe("Two-sided exact Student-t p-value (< 0.05)."), + .optional() + .describe( + "v1.42 — effective sample size after serial correlation (Pyper–Peterman), ≤ `n`. The p-value, the shrinkage and the tier use it.", + ), + r: z + .number() + .describe( + "Pearson r over the lag-joined daily series, computed on residuals: a linear trend (and, from 120 paired days, one annual harmonic) is removed from both series first (v1.42).", + ), + pValue: z + .number() + .describe( + "Two-sided exact Student-t p-value (< 0.05) at the effective sample size.", + ), qValue: z .number() .describe( @@ -497,7 +508,12 @@ export const discoveredCorrelation = z interpretation: z .string() .describe("Conservative, descriptive interpretation — never causal."), - lagDays: z.number().int().describe("Lag in days applied (1)."), + lagDays: z + .number() + .int() + .describe( + "Lag in days applied: 1 = behaviour day, outcome the next day; 0 = the environmental channels, whose exposure is the mean of the day before and the day itself, paired with the same day's outcome (v1.42).", + ), window: z .enum(["retrospective", "recent"]) .optional()