From aca497d9ac1e2d7eca52ec1d1ce68ab55cc6ba08 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc-Andr=C3=A9=20Bombeck?= Date: Thu, 8 Oct 2026 08:32:15 +0200 Subject: [PATCH] fix(insights): correlate residuals, not the calendar Two daily series that both follow the season (temperature and blood pressure, daylight and sleep) correlated strongly with nothing between them. In a seeded simulation of such a pair the engine surfaced a finding in 58 % of 400-day runs and 30 % of 180-day runs, and a real next-day effect under the season came back with the wrong sign. Every daily pair now goes through seasonallyAdjustedPearson: a linear trend and one annual harmonic are removed from both series (trend only under 120 paired days), the residuals are correlated, and the p-value uses the Pyper-Peterman effective sample size. Shrinkage and the phrasing tier read the effective size too. BH-FDR and the effect floor are unchanged. Pairs too short to fit the basis are not tested. The weather channels are tested once over the mean of the day before and the day itself against the same day's outcome (lagDays 0), which pairs the overnight low with the night it shaped now that sleep is keyed on the wake day. Discovered pairs carry an additive nEff, the same-day templates are wired, and the matrix skips the weather read when the environment module is off. scripts/compare-correlation-findings.ts counts findings per account before and after, read-only, for the release-candidate dry run. --- docs/api/openapi.yaml | 17 +- scripts/compare-correlation-findings.ts | 338 ++++++++++++++++++ src/__tests__/i18n-reverse-coverage.test.ts | 2 - src/lib/ai/coach/tools/correlations-read.ts | 7 +- src/lib/environment/fields.ts | 26 +- .../correlation-channel-series.test.ts | 54 +++ .../__tests__/correlation-discovery.test.ts | 155 ++++++-- .../correlation-seasonal-simulation.test.ts | 90 +++++ .../insights/__tests__/correlations.test.ts | 26 +- .../discovery-matrix-modules.test.ts | 13 + .../__tests__/helpers/seeded-series.ts | 164 +++++++++ .../__tests__/seasonal-adjust.test.ts | 282 +++++++++++++++ .../insights/correlation-channel-series.ts | 38 +- src/lib/insights/correlation-discovery.ts | 94 ++++- src/lib/insights/correlations.ts | 158 +++++++- src/lib/insights/discovery-matrix.ts | 20 +- .../__tests__/period-narrative.test.ts | 13 +- src/lib/insights/seasonal-adjust.ts | 315 ++++++++++++++++ src/lib/openapi/routes/cycle.ts | 6 + src/lib/openapi/routes/insights/schemas.ts | 28 +- 20 files changed, 1751 insertions(+), 95 deletions(-) create mode 100644 scripts/compare-correlation-findings.ts create mode 100644 src/lib/insights/__tests__/correlation-seasonal-simulation.test.ts create mode 100644 src/lib/insights/__tests__/helpers/seeded-series.ts create mode 100644 src/lib/insights/__tests__/seasonal-adjust.test.ts create mode 100644 src/lib/insights/seasonal-adjust.ts 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()