From cfd7c18d48b0a09d15fac0669e659f91f18503c3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marc-Andr=C3=A9=20Bombeck?= Date: Wed, 21 Oct 2026 23:20:28 +0200 Subject: [PATCH] fix(medications): widen the medication-effect intervals to a week of lag The automatic Newey-West rule gives a lag of 3 on 30 to 270 daily values, which leaves most of a physiological series' carry-over out of the intervals. The lag is now at least 7. Reference values are recomputed with statsmodels at the new default, and the lag-3 values stay pinned with an explicit lag. A fit whose residual autocorrelation exceeds 0.5 is flagged, and the card, the clinician view and the PDF say the ranges may still be too narrow. --- docs/api/openapi.yaml | 6 ++ messages/de.json | 6 +- messages/en.json | 6 +- messages/es.json | 6 +- messages/fr.json | 6 +- messages/it.json | 6 +- messages/ko.json | 6 +- messages/nl.json | 6 +- messages/pl.json | 6 +- .../clinician/activity-effect-sections.tsx | 5 ++ .../change-at-start-summary.test.tsx | 17 +++++- .../efficacy/change-at-start-summary.tsx | 8 +++ .../activity-effect-sections.test.ts | 31 +++++++++- .../activity-effect-sections.ts | 3 + .../__tests__/segmented-effect.test.ts | 19 +++++++ .../medications/efficacy/segmented-effect.ts | 10 ++++ .../__tests__/segmented-regression.test.ts | 57 +++++++++++++++++-- src/lib/stats/segmented-regression.ts | 23 ++++++-- src/lib/validations/medication-efficacy.ts | 12 +++- 19 files changed, 208 insertions(+), 31 deletions(-) diff --git a/docs/api/openapi.yaml b/docs/api/openapi.yaml index 9c77d2a05..9effbc1d3 100644 --- a/docs/api/openapi.yaml +++ b/docs/api/openapi.yaml @@ -54898,12 +54898,18 @@ components: type: integer minimum: -9007199254740991 maximum: 9007199254740991 + description: "Newey–West lag applied: the automatic rule, at least 7 (a week of daily carry-over)." + autocorrelationCaveat: + type: boolean + description: "True when the lag-1 residual autocorrelation exceeds 0.5: the carry-over is strong enough that the + intervals may still be too narrow, and the read-out says so." required: - n - df - rSquared - residualAutocorrelation - hacLag + - autocorrelationCaveat additionalProperties: false - type: "null" confounders: diff --git a/messages/de.json b/messages/de.json index 47417f3ab..f3e2a8458 100644 --- a/messages/de.json +++ b/messages/de.json @@ -431,7 +431,8 @@ "notTracked": "Einnahme nicht erfasst, daher keine Tage wegen ausgelassener Dosen entfernt", "endDoseChange": "Zeitraum endet am {date}, vor der ersten Dosisänderung", "days": "{before} Tage vor, {after} Tage nach dem Start, {excluded} Tage ausgelassen", - "method": "Segmentierte Regression über die Tageswerte von 90 Tagen vor bis höchstens 180 Tagen nach dem Start, Intervalle nach Newey-West. Tage nach dem Start, an denen in den 14 Tagen davor weniger als 80 % der Dosen eingenommen wurden, sind ausgelassen. Beschreibt einen zeitlichen Zusammenhang, keine nachgewiesene Wirkung des Medikaments." + "method": "Segmentierte Regression über die Tageswerte von 90 Tagen vor bis höchstens 180 Tagen nach dem Start, Intervalle nach Newey-West. Tage nach dem Start, an denen in den 14 Tagen davor weniger als 80 % der Dosen eingenommen wurden, sind ausgelassen. Beschreibt einen zeitlichen Zusammenhang, keine nachgewiesene Wirkung des Medikaments.", + "autocorrelation": "starke Abhängigkeit vom Vortag, die Bereiche können daher noch zu eng sein" } }, "errorBoundary": { @@ -2412,7 +2413,8 @@ "needStart": "Lege ein Startdatum fest, um die Veränderung zum Start zu sehen.", "needBefore": "Für die Veränderung zum Start braucht es Werte an mindestens {days} Tagen, verteilt über {span} Tage vor dem Start.", "needAfter": "Für die Veränderung zum Start braucht es Werte an mindestens {days} Tagen, verteilt über {span} Tage seit dem Start und an Tagen mit eingenommenen Dosen.", - "fitFailed": "Diese Werte ließen sich nicht auswerten, deshalb wird keine Veränderung zum Start gezeigt." + "fitFailed": "Diese Werte ließen sich nicht auswerten, deshalb wird keine Veränderung zum Start gezeigt.", + "autocorrelation": "Die Tageswerte hängen stark vom Vortag ab, deshalb können selbst diese Bereiche zu eng sein." } }, "import": { diff --git a/messages/en.json b/messages/en.json index 8f8fd6fa6..292300819 100644 --- a/messages/en.json +++ b/messages/en.json @@ -431,7 +431,8 @@ "notTracked": "intake not tracked, so no days left out for missed doses", "endDoseChange": "window ends {date}, before the first dose change", "days": "{before} days before, {after} days after the start, {excluded} days left out", - "method": "Segmented regression on daily values from 90 days before to at most 180 days after the start, with Newey-West intervals. Days after the start on which fewer than 80 % of the doses of the preceding 14 days were taken are left out. Describes an association over time, not a proven effect of the medication." + "method": "Segmented regression on daily values from 90 days before to at most 180 days after the start, with Newey-West intervals. Days after the start on which fewer than 80 % of the doses of the preceding 14 days were taken are left out. Describes an association over time, not a proven effect of the medication.", + "autocorrelation": "strong day-to-day carry-over, so the ranges may still be too narrow" } }, "errorBoundary": { @@ -2412,7 +2413,8 @@ "needStart": "Set a start date to see the change at the start.", "needBefore": "The change at the start needs readings on at least {days} days spread over {span} days before the start.", "needAfter": "The change at the start needs readings on at least {days} days spread over {span} days since the start, on days the doses were taken.", - "fitFailed": "These readings could not be fitted, so no change at the start is shown." + "fitFailed": "These readings could not be fitted, so no change at the start is shown.", + "autocorrelation": "The daily values carry over strongly from one day to the next, so even these ranges may be too narrow." } }, "import": { diff --git a/messages/es.json b/messages/es.json index 02f80d63e..2e222ca60 100644 --- a/messages/es.json +++ b/messages/es.json @@ -431,7 +431,8 @@ "notTracked": "toma no registrada, así que no se omitió ningún día por dosis olvidadas", "endDoseChange": "el periodo termina el {date}, antes del primer cambio de dosis", "days": "{before} días antes, {after} días después del inicio, {excluded} días omitidos", - "method": "Regresión segmentada sobre los valores diarios desde 90 días antes hasta como máximo 180 días después del inicio, con intervalos de Newey-West. Se omiten los días posteriores al inicio en los que se tomó menos del 80 % de las dosis de los 14 días anteriores. Describe una asociación en el tiempo, no un efecto demostrado del medicamento." + "method": "Regresión segmentada sobre los valores diarios desde 90 días antes hasta como máximo 180 días después del inicio, con intervalos de Newey-West. Se omiten los días posteriores al inicio en los que se tomó menos del 80 % de las dosis de los 14 días anteriores. Describe una asociación en el tiempo, no un efecto demostrado del medicamento.", + "autocorrelation": "fuerte dependencia del día anterior, así que los rangos pueden seguir siendo demasiado estrechos" } }, "errorBoundary": { @@ -2412,7 +2413,8 @@ "needStart": "Indica una fecha de inicio para ver el cambio al inicio.", "needBefore": "El cambio al inicio necesita valores de al menos {days} días repartidos en {span} días antes del inicio.", "needAfter": "El cambio al inicio necesita valores de al menos {days} días repartidos en {span} días desde el inicio, en días con las dosis tomadas.", - "fitFailed": "No se pudieron ajustar estos valores, así que no se muestra el cambio al inicio." + "fitFailed": "No se pudieron ajustar estos valores, así que no se muestra el cambio al inicio.", + "autocorrelation": "Los valores diarios dependen mucho del día anterior, así que incluso estos rangos pueden ser demasiado estrechos." } }, "import": { diff --git a/messages/fr.json b/messages/fr.json index 70036c66d..ce2b8f4ff 100644 --- a/messages/fr.json +++ b/messages/fr.json @@ -431,7 +431,8 @@ "notTracked": "prises non suivies, aucun jour écarté pour des doses manquées", "endDoseChange": "la période s'arrête le {date}, avant le premier changement de dose", "days": "{before} jours avant, {after} jours après le démarrage, {excluded} jours écartés", - "method": "Régression segmentée sur les valeurs quotidiennes, de 90 jours avant à 180 jours au plus après le démarrage, intervalles de Newey-West. Les jours après le démarrage où moins de 80 % des doses des 14 jours précédents ont été prises sont écartés. Décrit une association dans le temps, pas un effet démontré du médicament." + "method": "Régression segmentée sur les valeurs quotidiennes, de 90 jours avant à 180 jours au plus après le démarrage, intervalles de Newey-West. Les jours après le démarrage où moins de 80 % des doses des 14 jours précédents ont été prises sont écartés. Décrit une association dans le temps, pas un effet démontré du médicament.", + "autocorrelation": "forte dépendance à la veille, les plages peuvent donc rester trop étroites" } }, "errorBoundary": { @@ -2412,7 +2413,8 @@ "needStart": "Indique une date de début pour voir le changement au démarrage.", "needBefore": "Le changement au démarrage demande des valeurs sur au moins {days} jours répartis sur {span} jours avant le démarrage.", "needAfter": "Le changement au démarrage demande des valeurs sur au moins {days} jours répartis sur {span} jours depuis le démarrage, les jours où les doses ont été prises.", - "fitFailed": "Ces valeurs n'ont pas pu être ajustées, aucun changement au démarrage n'est donc affiché." + "fitFailed": "Ces valeurs n'ont pas pu être ajustées, aucun changement au démarrage n'est donc affiché.", + "autocorrelation": "Les valeurs quotidiennes dépendent fortement de la veille, donc même ces plages peuvent être trop étroites." } }, "import": { diff --git a/messages/it.json b/messages/it.json index 362937e74..f7a72c126 100644 --- a/messages/it.json +++ b/messages/it.json @@ -431,7 +431,8 @@ "notTracked": "assunzione non registrata, quindi nessun giorno escluso per dosi saltate", "endDoseChange": "il periodo termina il {date}, prima della prima variazione di dose", "days": "{before} giorni prima, {after} giorni dopo l'inizio, {excluded} giorni esclusi", - "method": "Regressione segmentata sui valori giornalieri da 90 giorni prima a non oltre 180 giorni dopo l'inizio, con intervalli di Newey-West. Sono esclusi i giorni dopo l'inizio in cui è stato assunto meno dell'80 % delle dosi dei 14 giorni precedenti. Descrive un'associazione nel tempo, non un effetto dimostrato del farmaco." + "method": "Regressione segmentata sui valori giornalieri da 90 giorni prima a non oltre 180 giorni dopo l'inizio, con intervalli di Newey-West. Sono esclusi i giorni dopo l'inizio in cui è stato assunto meno dell'80 % delle dosi dei 14 giorni precedenti. Descrive un'associazione nel tempo, non un effetto dimostrato del farmaco.", + "autocorrelation": "forte dipendenza dal giorno precedente, quindi gli intervalli possono restare troppo stretti" } }, "errorBoundary": { @@ -2412,7 +2413,8 @@ "needStart": "Imposta una data di inizio per vedere la variazione all'inizio.", "needBefore": "La variazione all'inizio richiede valori in almeno {days} giorni distribuiti su {span} giorni prima dell'inizio.", "needAfter": "La variazione all'inizio richiede valori in almeno {days} giorni distribuiti su {span} giorni dall'inizio, nei giorni con le dosi assunte.", - "fitFailed": "Non è stato possibile elaborare questi valori, quindi non viene mostrata alcuna variazione all'inizio." + "fitFailed": "Non è stato possibile elaborare questi valori, quindi non viene mostrata alcuna variazione all'inizio.", + "autocorrelation": "I valori giornalieri dipendono molto dal giorno precedente, quindi anche questi intervalli possono essere troppo stretti." } }, "import": { diff --git a/messages/ko.json b/messages/ko.json index 4ce92c2d9..dce35a661 100644 --- a/messages/ko.json +++ b/messages/ko.json @@ -431,7 +431,8 @@ "notTracked": "복용 미기록으로 놓친 용량에 따른 제외 없음", "endDoseChange": "기간은 첫 용량 변경 전인 {date}에 끝남", "days": "시작 전 {before}일, 시작 후 {after}일, 제외 {excluded}일", - "method": "시작 90일 전부터 최대 180일 후까지의 일별 값에 대한 분절 회귀, Newey-West 구간. 직전 14일 복용률이 80% 미만인 시작 후 날짜는 제외했습니다. 시간에 따른 연관성을 보여 줄 뿐, 약의 효과를 입증하지 않습니다." + "method": "시작 90일 전부터 최대 180일 후까지의 일별 값에 대한 분절 회귀, Newey-West 구간. 직전 14일 복용률이 80% 미만인 시작 후 날짜는 제외했습니다. 시간에 따른 연관성을 보여 줄 뿐, 약의 효과를 입증하지 않습니다.", + "autocorrelation": "전날 값에 크게 좌우되므로 범위가 여전히 너무 좁을 수 있음" } }, "errorBoundary": { @@ -2412,7 +2413,8 @@ "needStart": "시작 시점의 변화를 보려면 시작일을 설정하세요.", "needBefore": "시작 시점의 변화를 보려면 시작 전 {span}일에 걸쳐 최소 {days}일의 측정값이 필요합니다.", "needAfter": "시작 시점의 변화를 보려면 시작 후 {span}일에 걸쳐, 용량을 복용한 날 기준으로 최소 {days}일의 측정값이 필요합니다.", - "fitFailed": "이 측정값으로는 계산할 수 없어 시작 시점의 변화를 표시하지 않습니다." + "fitFailed": "이 측정값으로는 계산할 수 없어 시작 시점의 변화를 표시하지 않습니다.", + "autocorrelation": "하루 값이 전날 값에 크게 좌우되므로 이 범위도 너무 좁을 수 있습니다." } }, "import": { diff --git a/messages/nl.json b/messages/nl.json index dbe122e64..1af10b087 100644 --- a/messages/nl.json +++ b/messages/nl.json @@ -431,7 +431,8 @@ "notTracked": "inname niet bijgehouden, dus geen dagen weggelaten om gemiste doses", "endDoseChange": "periode eindigt op {date}, vóór de eerste dosiswijziging", "days": "{before} dagen voor, {after} dagen na de start, {excluded} dagen weggelaten", - "method": "Gesegmenteerde regressie op de dagwaarden van 90 dagen vóór tot hoogstens 180 dagen na de start, met Newey-West-intervallen. Dagen na de start waarop minder dan 80% van de doses van de 14 dagen ervoor is ingenomen, zijn weggelaten. Beschrijft een samenhang in de tijd, geen bewezen effect van het medicijn." + "method": "Gesegmenteerde regressie op de dagwaarden van 90 dagen vóór tot hoogstens 180 dagen na de start, met Newey-West-intervallen. Dagen na de start waarop minder dan 80% van de doses van de 14 dagen ervoor is ingenomen, zijn weggelaten. Beschrijft een samenhang in de tijd, geen bewezen effect van het medicijn.", + "autocorrelation": "sterke afhankelijkheid van de vorige dag, dus de bereiken kunnen nog te smal zijn" } }, "errorBoundary": { @@ -2412,7 +2413,8 @@ "needStart": "Stel een startdatum in om de verandering bij de start te zien.", "needBefore": "Voor de verandering bij de start zijn waarden op minstens {days} dagen nodig, verspreid over {span} dagen vóór de start.", "needAfter": "Voor de verandering bij de start zijn waarden op minstens {days} dagen nodig, verspreid over {span} dagen sinds de start en op dagen waarop de doses zijn ingenomen.", - "fitFailed": "Deze waarden konden niet worden verwerkt, dus er wordt geen verandering bij de start getoond." + "fitFailed": "Deze waarden konden niet worden verwerkt, dus er wordt geen verandering bij de start getoond.", + "autocorrelation": "De dagwaarden hangen sterk af van de vorige dag, dus zelfs deze bereiken kunnen te smal zijn." } }, "import": { diff --git a/messages/pl.json b/messages/pl.json index 3f5960574..ccbf60069 100644 --- a/messages/pl.json +++ b/messages/pl.json @@ -431,7 +431,8 @@ "notTracked": "przyjmowanie nie jest rejestrowane, więc nie pominięto dni z powodu opuszczonych dawek", "endDoseChange": "okres kończy się {date}, przed pierwszą zmianą dawki", "days": "{before} dni przed, {after} dni po rozpoczęciu, pominięto {excluded} dni", - "method": "Regresja segmentowa na wartościach dziennych od 90 dni przed do najwyżej 180 dni po rozpoczęciu, przedziały Neweya-Westa. Pominięto dni po rozpoczęciu, w których w ciągu poprzednich 14 dni przyjęto mniej niż 80% dawek. Opisuje związek w czasie, a nie udowodnione działanie leku." + "method": "Regresja segmentowa na wartościach dziennych od 90 dni przed do najwyżej 180 dni po rozpoczęciu, przedziały Neweya-Westa. Pominięto dni po rozpoczęciu, w których w ciągu poprzednich 14 dni przyjęto mniej niż 80% dawek. Opisuje związek w czasie, a nie udowodnione działanie leku.", + "autocorrelation": "silna zależność od poprzedniego dnia, więc zakresy mogą nadal być zbyt wąskie" } }, "errorBoundary": { @@ -2412,7 +2413,8 @@ "needStart": "Ustaw datę rozpoczęcia, aby zobaczyć zmianę na początku.", "needBefore": "Do zmiany na początku potrzeba pomiarów z co najmniej {days} dni rozłożonych na {span} dni przed rozpoczęciem.", "needAfter": "Do zmiany na początku potrzeba pomiarów z co najmniej {days} dni rozłożonych na {span} dni od rozpoczęcia, w dni z przyjętymi dawkami.", - "fitFailed": "Nie udało się dopasować tych pomiarów, więc zmiana na początku nie jest pokazana." + "fitFailed": "Nie udało się dopasować tych pomiarów, więc zmiana na początku nie jest pokazana.", + "autocorrelation": "Wartości dzienne silnie zależą od poprzedniego dnia, więc nawet te zakresy mogą być zbyt wąskie." } }, "import": { diff --git a/src/components/clinician/activity-effect-sections.tsx b/src/components/clinician/activity-effect-sections.tsx index 3f62c9878..a920eb3c7 100644 --- a/src/components/clinician/activity-effect-sections.tsx +++ b/src/components/clinician/activity-effect-sections.tsx @@ -167,6 +167,11 @@ export function MedicationEffectSection({ e.days.excludedLowAdherence + e.days.excludedNotTaking, })}

+ {e.fit?.autocorrelationCaveat ? ( +

+ {t("medications.efficacy.segmented.autocorrelation")} +

+ ) : null} {e.confounders.map((c) => (

{ expect(html).toContain("describes the starting dose"); expect(html).toContain("Example B started on"); expect(html).not.toMatch(/work(s|ing)|effective|successful/i); + expect(html).not.toContain("medication-wirkung-autocorrelation"); + }); + + it("warns when the daily values carry over strongly", () => { + const html = render({ + ...present, + fit: { + ...present.fit!, + residualAutocorrelation: 0.7, + autocorrelationCaveat: true, + }, + }); + expect(html).toContain("medication-wirkung-autocorrelation"); + expect(html).toContain("even these ranges may be too narrow"); }); it("names the missing floor instead of a number", () => { diff --git a/src/components/medications/detail/efficacy/change-at-start-summary.tsx b/src/components/medications/detail/efficacy/change-at-start-summary.tsx index 8f05249c0..af2c132fb 100644 --- a/src/components/medications/detail/efficacy/change-at-start-summary.tsx +++ b/src/components/medications/detail/efficacy/change-at-start-summary.tsx @@ -125,6 +125,14 @@ export function ChangeAtStartSummary({ {t("medications.efficacy.segmented.notTracked")}

) : null} + {segmented.fit?.autocorrelationCaveat ? ( +

+ {t("medications.efficacy.segmented.autocorrelation")} +

+ ) : null} {segmented.window.endReason === "dose_change" ? (

{t("medications.efficacy.segmented.endDoseChange")} diff --git a/src/lib/doctor-report-pdf/__tests__/activity-effect-sections.test.ts b/src/lib/doctor-report-pdf/__tests__/activity-effect-sections.test.ts index a2c0ac544..94daa7421 100644 --- a/src/lib/doctor-report-pdf/__tests__/activity-effect-sections.test.ts +++ b/src/lib/doctor-report-pdf/__tests__/activity-effect-sections.test.ts @@ -81,7 +81,8 @@ const EFFECT: EfficacySegmented = { df: 141, rSquared: 0.55, residualAutocorrelation: 0.31, - hacLag: 4, + hacLag: 7, + autocorrelationCaveat: false, }, confounders: [ { medicationName: "Other drug", event: "started", date: "2026-02-14" }, @@ -218,6 +219,34 @@ describe("doctor report v1.44 sections", () => { expect(text).toContain("5 days after the start left out"); expect(text).toContain("Other drug started on"); expect(text).toContain("not a proven effect"); + expect(text).not.toContain("carry-over"); + }); + + it("notes strong day-to-day carry-over in the effect notes", async () => { + const text = ( + await textOf( + render( + data({ + medicationEffects: [ + { + medicationName: "Example ACE inhibitor", + startsOn: "2026-03-02", + targets: [ + { + type: "BLOOD_PRESSURE_SYS", + effect: { + ...EFFECT, + fit: { ...EFFECT.fit!, autocorrelationCaveat: true }, + }, + }, + ], + }, + ], + }), + ), + ) + ).replace(/\s+/g, " "); + expect(text).toContain("strong day-to-day carry-over"); }); it("names the missing floor instead of printing numbers", async () => { diff --git a/src/lib/doctor-report-pdf/activity-effect-sections.ts b/src/lib/doctor-report-pdf/activity-effect-sections.ts index cbf7c0d36..819918bbc 100644 --- a/src/lib/doctor-report-pdf/activity-effect-sections.ts +++ b/src/lib/doctor-report-pdf/activity-effect-sections.ts @@ -232,6 +232,9 @@ export function buildMedicationEffectSection( if (e.adherence.status === "not_tracked") { extras.push(t("doctorReport.medicationEffect.notTracked")); } + if (e.fit?.autocorrelationCaveat) { + extras.push(t("doctorReport.medicationEffect.autocorrelation")); + } if (e.window.endReason === "dose_change") { extras.push( t("doctorReport.medicationEffect.endDoseChange", { diff --git a/src/lib/medications/efficacy/__tests__/segmented-effect.test.ts b/src/lib/medications/efficacy/__tests__/segmented-effect.test.ts index 6f099d6bc..ea4191657 100644 --- a/src/lib/medications/efficacy/__tests__/segmented-effect.test.ts +++ b/src/lib/medications/efficacy/__tests__/segmented-effect.test.ts @@ -173,6 +173,25 @@ describe("buildSegmentedEffect", () => { expect(r.present).toBe(true); expect(r.levelChange!.includesZero).toBe(true); }); + + it("flags strong residual autocorrelation for a caveat, not a clean fit", () => { + // A slow wave around the line: neighbouring residuals share a sign. + const wavy = series(-60, 59).map((p, i) => ({ + day: p.day, + value: p.value + 4 * Math.sin(i / 6), + })); + const r = buildSegmentedEffect(input({ dayValues: wavy })); + expect(r.fit!.residualAutocorrelation).toBeGreaterThan(0.5); + expect(r.fit!.autocorrelationCaveat).toBe(true); + expect(r.fit!.hacLag).toBe(7); + // Alternating noise: negative carry-over, no caveat. + const alternating = series(-60, 59).map((p, i) => ({ + day: p.day, + value: p.value + (i % 2 === 0 ? 2 : -2), + })); + const clean = buildSegmentedEffect(input({ dayValues: alternating })); + expect(clean.fit!.autocorrelationCaveat).toBe(false); + }); }); describe("notTakingPredicate", () => { diff --git a/src/lib/medications/efficacy/segmented-effect.ts b/src/lib/medications/efficacy/segmented-effect.ts index 028f4a55b..f102a8217 100644 --- a/src/lib/medications/efficacy/segmented-effect.ts +++ b/src/lib/medications/efficacy/segmented-effect.ts @@ -48,6 +48,12 @@ export const ADHERENCE_LOOKBACK_DAYS = 14; export const ADHERENCE_FLOOR = 0.8; /** Confidence level of every interval. */ export const SEGMENTED_LEVEL = 0.95; +/** + * Lag-1 residual autocorrelation above which the read-out carries a caveat: + * the carry-over is then strong enough that even the widened intervals may be + * too narrow. + */ +export const AUTOCORRELATION_CAVEAT = 0.5; /** An estimate in display terms: the unit per day for a level, per week for a slope. */ export interface EfficacyEstimate { @@ -106,6 +112,8 @@ export interface EfficacySegmented { rSquared: number; residualAutocorrelation: number; hacLag: number; + /** True above {@link AUTOCORRELATION_CAVEAT}: renderers add the caveat line. */ + autocorrelationCaveat: boolean; } | null; /** Other medication changes inside the window, oldest first. */ confounders: EfficacyConfounder[]; @@ -311,6 +319,8 @@ export function buildSegmentedEffect( rSquared: fit.rSquared, residualAutocorrelation: fit.residualAutocorrelation, hacLag: fit.hacLag, + autocorrelationCaveat: + fit.residualAutocorrelation > AUTOCORRELATION_CAVEAT, }, confounders, }; diff --git a/src/lib/stats/__tests__/segmented-regression.test.ts b/src/lib/stats/__tests__/segmented-regression.test.ts index 57001444b..2c057825b 100644 --- a/src/lib/stats/__tests__/segmented-regression.test.ts +++ b/src/lib/stats/__tests__/segmented-regression.test.ts @@ -2,16 +2,22 @@ * Exact statistics for the segmented (interrupted time series) regression. * * The reference values for the noisy fixtures were computed independently - * with statsmodels 0.14 (`OLS(...).fit(cov_type="HAC", cov_kwds={"maxlags": - * L, "use_correction": True}, use_t=True)`, L = floor(4·(n/100)^(2/9))) and - * scipy's `stats.t`. The fixture generator below is the same deterministic + * with statsmodels 0.14.4 (`OLS(...).fit(cov_type="HAC", cov_kwds={"maxlags": + * L, "use_correction": True}, use_t=True)`) and scipy's `stats.t`, at the + * default lag L = max(7, floor(4·(n/100)^(2/9))) = 7 and, to pin the HAC meat + * at a second lag, at the automatic L = 3 passed explicitly. To recompute: + * build the fixture below in Python (Python's floored `%`), fit the design + * [1, t, D, t·D] with the call above at each L, and read `params`, `bse` and + * `conf_int(0.05)`; the coefficients and R² do not depend on L. The fixture generator below is the same deterministic * formula the reference script used, so a drift in the estimator, the HAC * meat, the small-sample factor or the t quantile fails here. */ import { describe, expect, it } from "vitest"; import { + defaultHacLag, fitSegmentedRegression, + MIN_HAC_LAG, neweyWestLag, studentTCdf, studentTQuantile, @@ -62,10 +68,27 @@ describe("fitSegmentedRegression", () => { expect(fit!.rSquared).toBeCloseTo(1, 12); }); - it("matches statsmodels HAC on a contiguous 90-day series", () => { + it("matches statsmodels HAC on a contiguous 90-day series at the default lag", () => { const fit = fitSegmentedRegression(fixture(40, 50))!; expect(fit.n).toBe(90); expect(fit.df).toBe(86); + expect(fit.hacLag).toBe(7); + expect(fit.levelChange.estimate).toBeCloseTo(-7.336815685253282, 9); + expect(fit.intercept.standardError).toBeCloseTo(0.3926378134611552, 9); + expect(fit.preSlope.standardError).toBeCloseTo(0.015434960504949637, 10); + expect(fit.levelChange.standardError).toBeCloseTo(0.5544100003043471, 9); + expect(fit.slopeChange.standardError).toBeCloseTo(0.02119259135971675, 10); + expect(fit.levelChange.ciLow).toBeCloseTo(-8.43894628913928, 8); + expect(fit.levelChange.ciHigh).toBeCloseTo(-6.234685081367285, 8); + expect(fit.slopeChange.ciLow).toBeCloseTo(-0.07897452738576954, 9); + expect(fit.slopeChange.ciHigh).toBeCloseTo(0.005284427179883318, 9); + expect(fit.rSquared).toBeCloseTo(0.8683372493940129, 10); + }); + + it("matches statsmodels HAC on a contiguous 90-day series at lag 3", () => { + const fit = fitSegmentedRegression(fixture(40, 50), { hacLag: 3 })!; + expect(fit.n).toBe(90); + expect(fit.df).toBe(86); expect(fit.hacLag).toBe(3); expect(fit.intercept.estimate).toBeCloseTo(139.46950523140714, 9); expect(fit.preSlope.estimate).toBeCloseTo(-0.09725785860840647, 10); @@ -82,9 +105,24 @@ describe("fitSegmentedRegression", () => { expect(fit.rSquared).toBeCloseTo(0.8683372493940129, 10); }); - it("matches statsmodels HAC on a series with gaps", () => { + it("matches statsmodels HAC on a series with gaps at the default lag", () => { const fit = fitSegmentedRegression(fixture(30, 30, { gaps: true }))!; expect(fit.n).toBe(48); + expect(fit.hacLag).toBe(7); + expect(fit.levelChange.estimate).toBeCloseTo(-6.123258149435133, 9); + expect(fit.levelChange.standardError).toBeCloseTo(0.7919225818543464, 9); + expect(fit.levelChange.ciLow).toBeCloseTo(-7.71927324237417, 8); + expect(fit.levelChange.ciHigh).toBeCloseTo(-4.527243056496095, 8); + expect(fit.slopeChange.estimate).toBeCloseTo(-0.11857109294626766, 10); + expect(fit.slopeChange.ciLow).toBeCloseTo(-0.22143394263874167, 9); + expect(fit.slopeChange.ciHigh).toBeCloseTo(-0.01570824325379365, 9); + }); + + it("matches statsmodels HAC on a series with gaps at lag 3", () => { + const fit = fitSegmentedRegression(fixture(30, 30, { gaps: true }), { + hacLag: 3, + })!; + expect(fit.n).toBe(48); expect(fit.hacLag).toBe(3); expect(fit.levelChange.estimate).toBeCloseTo(-6.123258149435133, 9); expect(fit.levelChange.standardError).toBeCloseTo(0.9608527406825199, 9); @@ -106,9 +144,16 @@ describe("fitSegmentedRegression", () => { expect(fitSegmentedRegression(fixture(2, 2, { noise: false }))).toBeNull(); }); - it("uses the automatic Newey–West lag", () => { + it("uses the automatic Newey–West lag, never below a week", () => { expect(neweyWestLag(90)).toBe(3); expect(neweyWestLag(30)).toBe(3); expect(neweyWestLag(10)).toBe(2); + expect(MIN_HAC_LAG).toBe(7); + expect(defaultHacLag(90)).toBe(7); + expect(defaultHacLag(270)).toBe(7); + expect(defaultHacLag(100_000)).toBe(neweyWestLag(100_000)); + // A short series clamps to n − 1. + expect(fitSegmentedRegression(fixture(4, 4))!.hacLag).toBe(7); + expect(fitSegmentedRegression(fixture(3, 3))!.hacLag).toBe(5); }); }); diff --git a/src/lib/stats/segmented-regression.ts b/src/lib/stats/segmented-regression.ts index b0451c5cd..e5e81eb23 100644 --- a/src/lib/stats/segmented-regression.ts +++ b/src/lib/stats/segmented-regression.ts @@ -13,11 +13,14 @@ * Daily health series are autocorrelated: today's blood pressure says a lot * about tomorrow's. Plain OLS standard errors ignore that and come out too * narrow, so the covariance here is the Newey–West estimator (Bartlett - * kernel) with the usual automatic lag `floor(4·(n/100)^(2/9))` and the - * small-sample factor n/(n−k). Intervals use Student's t with n−k degrees of - * freedom. This matches `statsmodels` OLS `cov_type="HAC"` with - * `use_correction=True` and `use_t=True`; the unit tests pin it to values - * computed there. + * kernel) and the small-sample factor n/(n−k). The lag is the usual automatic + * `floor(4·(n/100)^(2/9))`, but never below {@link MIN_HAC_LAG}: on 30 to 270 + * daily values the automatic rule gives 3, while a physiological series + * carries over for about a week (weight, resting heart rate, blood pressure + * under a new medication), and three lags leave most of that out. Intervals + * use Student's t with n−k degrees of freedom. This matches `statsmodels` OLS + * `cov_type="HAC"` with `use_correction=True` and `use_t=True` at the same + * `maxlags`; the unit tests pin it to values computed there. * * Observations are taken in the order given, which for a daily series with * gaps treats a gap as adjacent days. That is a known approximation of the @@ -61,11 +64,19 @@ export interface SegmentedFit { const K = 4; +/** Smallest Newey–West lag applied: a week of daily carry-over. */ +export const MIN_HAC_LAG = 7; + /** Automatic Newey–West lag: floor(4·(n/100)^(2/9)). */ export function neweyWestLag(n: number): number { return Math.floor(4 * Math.pow(n / 100, 2 / 9)); } +/** The lag applied by default: the automatic rule, at least {@link MIN_HAC_LAG}. */ +export function defaultHacLag(n: number): number { + return Math.max(MIN_HAC_LAG, neweyWestLag(n)); +} + /** Invert a small symmetric positive-definite matrix (Gauss–Jordan). */ function invert(matrix: number[][]): number[][] | null { const size = matrix.length; @@ -131,7 +142,7 @@ export function fitSegmentedRegression( ); // Meat: Σ e_t² x_t x_tᵀ + Σ_l w_l Σ_t e_t e_{t−l} (x_t x_{t−l}ᵀ + x_{t−l} x_tᵀ) - const lag = Math.max(0, Math.min(options.hacLag ?? neweyWestLag(n), n - 1)); + const lag = Math.max(0, Math.min(options.hacLag ?? defaultHacLag(n), n - 1)); const meat = Array.from({ length: K }, () => new Array(K).fill(0)); for (let l = 0; l <= lag; l++) { const w = l === 0 ? 1 : 1 - l / (lag + 1); diff --git a/src/lib/validations/medication-efficacy.ts b/src/lib/validations/medication-efficacy.ts index 743e3b42a..ba6b075a1 100644 --- a/src/lib/validations/medication-efficacy.ts +++ b/src/lib/validations/medication-efficacy.ts @@ -137,7 +137,17 @@ const segmentedSchema = z df: z.number().int(), rSquared: z.number(), residualAutocorrelation: z.number(), - hacLag: z.number().int(), + hacLag: z + .number() + .int() + .describe( + "Newey–West lag applied: the automatic rule, at least 7 (a week of daily carry-over).", + ), + autocorrelationCaveat: z + .boolean() + .describe( + "True when the lag-1 residual autocorrelation exceeds 0.5: the carry-over is strong enough that the intervals may still be too narrow, and the read-out says so.", + ), }) .nullable(), confounders: z.array(