diff --git a/PR_PHASE2_DEBATE.md b/PR_PHASE2_DEBATE.md new file mode 100644 index 0000000..22c3862 --- /dev/null +++ b/PR_PHASE2_DEBATE.md @@ -0,0 +1,33 @@ +Adds basic incidence computation from a line list (daily and weekly CDC/MMWR). + +Key points +- Add `epydem.incidence(df, date_col=..., freq=...)`. +- Supports: + - `freq="D"` -> daily counts by calendar date + - `freq="W-MMWR"` -> weekly counts by CDC/MMWR epiweek (adds epi_year + epi_week) +- Optional stratification via `by=[...]`. + +Why this implementation +- Incidence is the most common first step for line-list analysis; we want a minimal, composable primitive. +- Using `epi_year` + `epi_week` avoids ambiguity at year boundaries. +- Returning a tidy DataFrame (group columns + `cases`) makes plotting and further transforms straightforward in pandas. +- We intentionally do not implement missing-week filling, rolling averages, or cumulative sums yet to keep the first API small. + +Multi-role debate (differences, not consensus) + +Role A — pragmatic developer +- 👍 Likes: minimal API, returns tidy DF, works with `by` strata. +- ⚠️ Concern: performance if we iterate row-by-row for epiweek on large datasets; may need vectorization/caching later. + +Role B — architecture +- 👍 Likes: keeps time semantics centralized in `epydem.time` and reuses `epiweek()`. +- ⚠️ Concern: frequency naming: `W-MMWR` is clear but may expand to ISO/WHO; may want a more general `freq="W"` + `system=`. + +Role C — developer user (DX) +- 👍 Likes: one-liner to get epicurve-like table. +- ⚠️ Concern: expects convenience features soon (fill missing weeks, ordering, start/end bounds, cumulative) and clear examples in docs. + +Points of divergence to revisit later +1) Performance: pure-Python loop vs vectorized epiweek mapping. +2) Output schema: keep tidy long DF vs allow pivot/wide convenience. +3) Feature growth: add `fill_missing`, `cumulative`, `rolling`, and `freq` generalization. diff --git a/README.md b/README.md index e279414..59597c4 100644 --- a/README.md +++ b/README.md @@ -30,10 +30,14 @@ print(year, week) week_only = epydem.epiweek_number("2024-01-01") print(week_only) + +# Incidence (line list -> counts) +# df = pandas.DataFrame({"onset_date": [...], "sex": [...]}) +# weekly = epydem.incidence(df, date_col="onset_date", freq="W-MMWR", by=["sex"]) ``` ## Roadmap (high level) - CDC/MMWR epiweek (week starts Sunday; week 1 contains Jan 4) ✅ (in progress) -- Incidence / epicurves from line lists (pandas) +- Incidence / epicurves from line lists (pandas) ✅ (in progress) - Summary statistics (descriptive epi) diff --git a/epydem/__init__.py b/epydem/__init__.py index e5d148f..aa7026f 100644 --- a/epydem/__init__.py +++ b/epydem/__init__.py @@ -1,6 +1,7 @@ """epydem public API.""" from .epiweek import calculate +from .incidence import incidence from .time import epiweek, epiweek_number, mmwr_week, mmwr_week1_start, parse_ymd __all__ = [ @@ -10,4 +11,5 @@ "mmwr_week", "mmwr_week1_start", "parse_ymd", + "incidence", ] diff --git a/epydem/epiweek.py b/epydem/epiweek.py index 98139bb..4df3c32 100644 --- a/epydem/epiweek.py +++ b/epydem/epiweek.py @@ -10,7 +10,7 @@ from __future__ import annotations -from .time import epiweek, epiweek_number, mmwr_week, mmwr_week1_start, parse_ymd +from .time import epiweek_number def calculate(date_str: str) -> int: diff --git a/epydem/incidence.py b/epydem/incidence.py new file mode 100644 index 0000000..da37631 --- /dev/null +++ b/epydem/incidence.py @@ -0,0 +1,173 @@ +from __future__ import annotations + +from collections.abc import Sequence +from dataclasses import dataclass +from datetime import timedelta +from typing import Literal + +import pandas as pd + +from .time import epiweek, mmwr_week1_start + +EpiFreq = Literal["D", "W-MMWR"] +OutputFormat = Literal["long", "wide"] + + +@dataclass(frozen=True) +class IncidenceSpec: + date_col: str = "date" + freq: EpiFreq = "W-MMWR" + by: tuple[str, ...] = () + + +def incidence( + df: pd.DataFrame, + *, + date_col: str, + freq: EpiFreq = "W-MMWR", + by: Sequence[str] | None = None, + count_col: str = "cases", + output: OutputFormat = "wide", + fill_missing: bool = True, +) -> pd.DataFrame: + """Compute incidence counts from a line list. + + Args: + df: Line list dataframe. + date_col: Column containing dates. Supports values acceptable to `epydem.time.epiweek`. + freq: + - "D": daily counts by calendar date + - "W-MMWR": weekly counts by CDC/MMWR epiweek (returns epi_year + epi_week) + by: Optional stratification columns. + count_col: Name of count column in the output. + output: + - "wide" (default): pivot table style (DX-friendly) + - "long": tidy long-form table + fill_missing: If True, fill missing dates/weeks with 0 counts. + + Returns: + DataFrame in the requested output format. + + Notes: + - Performance: for weekly counts we compute epiweek for *unique* dates and map back, + avoiding a pure Python loop per row. + """ + + if by is None: + by_cols: list[str] = [] + else: + by_cols = list(by) + + if date_col not in df.columns: + raise KeyError(f"date_col not found: {date_col}") + + work = df.copy() + + if freq == "D": + work["date"] = pd.to_datetime(work[date_col]).dt.date + group_cols = by_cols + ["date"] + long = work.groupby(group_cols, dropna=False).size().rename(count_col).reset_index() + + if fill_missing: + all_dates = pd.date_range(long["date"].min(), long["date"].max(), freq="D").date + + if not by_cols: + long = ( + long.set_index("date") + .reindex(all_dates, fill_value=0) + .rename_axis("date") + .reset_index() + ) + else: + # Fill missing dates per stratum combination. + strata = long[by_cols].drop_duplicates() + strata["_k"] = 1 + full_dates = pd.DataFrame({"date": list(all_dates)}) + full_dates["_k"] = 1 + full = strata.merge(full_dates, on="_k", how="outer").drop(columns=["_k"]) + + long = full.merge(long, on=by_cols + ["date"], how="left") + long[count_col] = long[count_col].fillna(0).astype(int) + + if output == "long": + return long + + # wide + if by_cols: + wide = long.pivot_table( + index="date", + columns=by_cols, + values=count_col, + aggfunc="sum", + fill_value=0, + ) + else: + wide = long.set_index("date")[[count_col]] + + return wide.sort_index() + + if freq == "W-MMWR": + # Compute (epi_year, epi_week) for unique dates, then map back for performance. + uniq = pd.unique(work[date_col]) + mapping = {v: epiweek(v) for v in uniq} + epi_pairs = work[date_col].map(mapping) + + work["epi_year"] = epi_pairs.map(lambda t: t[0]) + work["epi_week"] = epi_pairs.map(lambda t: t[1]) + + group_cols = by_cols + ["epi_year", "epi_week"] + long = work.groupby(group_cols, dropna=False).size().rename(count_col).reset_index() + + if fill_missing: + # Fill missing epiweeks between min and max observed. + long = long.sort_values(["epi_year", "epi_week"]).reset_index(drop=True) + start_y, start_w = int(long.iloc[0]["epi_year"]), int(long.iloc[0]["epi_week"]) + end_y, end_w = int(long.iloc[-1]["epi_year"]), int(long.iloc[-1]["epi_week"]) + + # Build the full sequence of (y,w) by stepping Sundays. + start_date = mmwr_week1_start(start_y) + timedelta(days=(start_w - 1) * 7) + end_date = mmwr_week1_start(end_y) + timedelta(days=(end_w - 1) * 7) + + full_pairs = [] + d = start_date + while d <= end_date: + full_pairs.append(epiweek(d)) + d += timedelta(days=7) + + idx = pd.MultiIndex.from_tuples(full_pairs, names=["epi_year", "epi_week"]) + + if not by_cols: + long = ( + long.set_index(["epi_year", "epi_week"]) + .reindex(idx, fill_value=0) + .reset_index() + ) + else: + # Fill missing epiweeks per stratum combination. + strata = long[by_cols].drop_duplicates() + strata["_k"] = 1 + full_time = idx.to_frame(index=False) + full_time["_k"] = 1 + full = strata.merge(full_time, on="_k", how="outer").drop(columns=["_k"]) + + long = full.merge(long, on=by_cols + ["epi_year", "epi_week"], how="left") + long[count_col] = long[count_col].fillna(0).astype(int) + + if output == "long": + return long + + # wide + if by_cols: + wide = long.pivot_table( + index=["epi_year", "epi_week"], + columns=by_cols, + values=count_col, + aggfunc="sum", + fill_value=0, + ) + else: + wide = long.set_index(["epi_year", "epi_week"])[[count_col]] + + return wide.sort_index() + + raise ValueError(f"Unknown freq: {freq}") diff --git a/epydem/time.py b/epydem/time.py index 7251461..3aa8904 100644 --- a/epydem/time.py +++ b/epydem/time.py @@ -2,9 +2,9 @@ import re from datetime import date, datetime, timedelta -from typing import Literal, Tuple, Union +from typing import Literal -DateLike = Union[str, date, datetime] +DateLike = str | date | datetime EpiWeekSystem = Literal["mmwr"] _DATE_PATTERN = re.compile(r"^\d{4}-\d{2}-\d{2}$") @@ -47,7 +47,7 @@ def mmwr_week1_start(year: int) -> date: return jan4 - timedelta(days=days_since_sunday) -def mmwr_week(value: DateLike) -> Tuple[int, int]: +def mmwr_week(value: DateLike) -> tuple[int, int]: """Compute CDC/MMWR epidemiological week for a date. Returns: @@ -56,24 +56,32 @@ def mmwr_week(value: DateLike) -> Tuple[int, int]: Notes: - Weeks start Sunday. - Week 1 is the week containing Jan 4. - - Dates in early January can belong to the previous MMWR year. + - The MMWR year is *not* always the calendar year of the date. + Example: 2023-12-31 is the start of 2024 week 1. + + Implementation rule: + - Find the unique `mmwr_year` such that: + week1_start(mmwr_year) <= d < week1_start(mmwr_year + 1) """ d = parse_ymd(value) - start_this_year = mmwr_week1_start(d.year) - if d < start_this_year: - mmwr_year = d.year - 1 - start = mmwr_week1_start(mmwr_year) - else: - mmwr_year = d.year - start = start_this_year + year = d.year + start = mmwr_week1_start(year) + start_next = mmwr_week1_start(year + 1) + + if d < start: + year -= 1 + start = mmwr_week1_start(year) + elif d >= start_next: + year += 1 + start = start_next week = ((d - start).days // 7) + 1 - return mmwr_year, week + return year, week -def epiweek(value: DateLike, system: EpiWeekSystem = "mmwr") -> Tuple[int, int]: +def epiweek(value: DateLike, system: EpiWeekSystem = "mmwr") -> tuple[int, int]: """Compute epidemiological week. Currently supported systems: diff --git a/tests/test_epiweek.py b/tests/test_epiweek.py index 2ab4980..425ccc9 100644 --- a/tests/test_epiweek.py +++ b/tests/test_epiweek.py @@ -40,11 +40,12 @@ def test_known_boundaries_tuple_output(self): # New year transitions. assert epydem.epiweek("2023-01-01") == (2023, 1) - assert epydem.epiweek("2023-12-31") == (2023, 53) + # 2023-12-31 is the start of 2024 week 1 under CDC/MMWR. + assert epydem.epiweek("2023-12-31") == (2024, 1) assert epydem.epiweek("2024-01-01") == (2024, 1) assert epydem.epiweek("2024-01-07") == (2024, 2) - assert epydem.epiweek("2024-12-31") == (2024, 53) + assert epydem.epiweek("2024-12-31") == (2025, 1) assert epydem.epiweek("2026-01-01") == (2025, 53) diff --git a/tests/test_incidence.py b/tests/test_incidence.py new file mode 100644 index 0000000..fb30e56 --- /dev/null +++ b/tests/test_incidence.py @@ -0,0 +1,65 @@ +from __future__ import annotations + +import pandas as pd + +import epydem + + +def test_incidence_daily_basic(): + df = pd.DataFrame( + { + "onset_date": [ + "2024-01-01", + "2024-01-01", + "2024-01-02", + ], + "sex": ["M", "F", "M"], + } + ) + + out = epydem.incidence(df, date_col="onset_date", freq="D") + # default output is wide + assert list(out.columns) == ["cases"] + + assert out.loc[pd.to_datetime("2024-01-01").date(), "cases"] == 2 + assert out.loc[pd.to_datetime("2024-01-02").date(), "cases"] == 1 + + +def test_incidence_weekly_mmwr_with_strata_fills_missing(): + # Create a gap: only week 1 has cases for stratum B. + df = pd.DataFrame( + { + "onset_date": [ + "2024-01-01", # (2024,1) + "2024-01-01", + "2024-01-08", # (2024,2) + ], + "province": ["A", "B", "A"], + } + ) + + out = epydem.incidence( + df, + date_col="onset_date", + freq="W-MMWR", + by=["province"], + fill_missing=True, + ) + + # Both strata should have rows for both weeks. + assert out.loc[(2024, 1), "A"] == 1 + assert out.loc[(2024, 2), "A"] == 1 + + assert out.loc[(2024, 1), "B"] == 1 + assert out.loc[(2024, 2), "B"] == 0 + + # Long format remains available + long = epydem.incidence( + df, + date_col="onset_date", + freq="W-MMWR", + by=["province"], + output="long", + fill_missing=True, + ) + assert set(long.columns) == {"province", "epi_year", "epi_week", "cases"}