Phase 2: incidence() (daily + weekly MMWR)

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.

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)

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",
 ]

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:

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}")

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:

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)
 

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"}