mirror of
https://github.com/MBombeck/HealthLog.git
synced 2026-11-09 21:10:59 +00:00
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.
This commit is contained in:
1 parent
99dfd68aad
commit
cfd7c18d48
19 files changed
+208
-31
No files matched your search
@@ -167,6 +167,11 @@ export function MedicationEffectSection({
|
||||
e.days.excludedLowAdherence + e.days.excludedNotTaking,
|
||||
})}
|
||||
</p>
|
||||
{e.fit?.autocorrelationCaveat ? (
|
||||
<p className="text-muted-foreground text-xs">
|
||||
{t("medications.efficacy.segmented.autocorrelation")}
|
||||
</p>
|
||||
) : null}
|
||||
{e.confounders.map((c) => (
|
||||
<p
|
||||
key={`${c.medicationName}:${c.event}:${c.date}`}
|
||||
|
||||
+16
-1
@@ -37,7 +37,8 @@ const present: EfficacySegmented = {
|
||||
df: 146,
|
||||
rSquared: 0.6,
|
||||
residualAutocorrelation: 0.3,
|
||||
hacLag: 4,
|
||||
hacLag: 7,
|
||||
autocorrelationCaveat: false,
|
||||
},
|
||||
confounders: [
|
||||
{ medicationName: "Example B", event: "started", date: "2026-02-10" },
|
||||
@@ -70,6 +71,20 @@ describe("ChangeAtStartSummary", () => {
|
||||
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", () => {
|
||||
|
||||
@@ -125,6 +125,14 @@ export function ChangeAtStartSummary({
|
||||
{t("medications.efficacy.segmented.notTracked")}
|
||||
</p>
|
||||
) : null}
|
||||
{segmented.fit?.autocorrelationCaveat ? (
|
||||
<p
|
||||
className="text-muted-foreground text-xs"
|
||||
data-slot="medication-wirkung-autocorrelation"
|
||||
>
|
||||
{t("medications.efficacy.segmented.autocorrelation")}
|
||||
</p>
|
||||
) : null}
|
||||
{segmented.window.endReason === "dose_change" ? (
|
||||
<p className="text-muted-foreground text-xs">
|
||||
{t("medications.efficacy.segmented.endDoseChange")}
|
||||
|
||||
@@ -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 () => {
|
||||
|
||||
@@ -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", {
|
||||
|
||||
@@ -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", () => {
|
||||
|
||||
@@ -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,
|
||||
};
|
||||
|
||||
@@ -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);
|
||||
});
|
||||
});
|
||||
@@ -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<number>(K).fill(0));
|
||||
for (let l = 0; l <= lag; l++) {
|
||||
const w = l === 0 ? 1 : 1 - l / (lag + 1);
|
||||
|
||||
@@ -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(
|
||||
|
||||
Reference in new issue
Block a user