Fix for statsmodels 0.15.0 and add end of follow-up outcomes

.github/workflows/python-app.yml

@@ -23,7 +23,7 @@ jobs:
     - uses: actions/checkout@v7
 
     - name: Install uv
-      uses: astral-sh/setup-uv@v8.2.0
+      uses: astral-sh/setup-uv@v10.0.1
 
     - name: Set up Python ${{ matrix.python-version }}
       run: uv python install ${{ matrix.python-version }}

CITATION.cff

@@ -18,7 +18,7 @@ authors:
     affiliation: "CAUSALab, Department of Epidemiology, Harvard T.H. Chan School of Public Health; Department of Biostatistics, Harvard T.H. Chan School of Public Health"
     orcid: https://orcid.org/0000-0003-1619-8456
 title: "pySEQTarget: Sequential Target Trial Emulation"
-version: 0.14.0
+version: 0.14.1
 url: https://pyseqtarget.readthedocs.io/
 repository-code: https://github.com/CausalInference/pySEQTarget
 repository-artifact: https://pypi.org/project/pySEQTarget/

docs/index.rst

@@ -17,4 +17,5 @@
 
    vignettes/getting_started
    vignettes/more_advanced_models
+   vignettes/end_of_followup
    vignettes/exploring_results

pySEQTarget/SEQopts.py

@@ -26,6 +26,26 @@ class SEQopts:
         order) fits a separate denominator model, with its own covariates, in each
         arm; this is only supported for post-expansion weights
         (``weight_preexpansion=False``).
+    :param end_of_fup: Boolean to estimate an end-of-follow-up outcome — one measured
+        at a single follow-up time rather than as a time-to-event — instead of fitting
+        a survival outcome model, default ``False``. The estimate is the weighted
+        average of the outcome within each baseline treatment arm, weighted by the
+        period-trial-specific weight at the time the outcome is taken. Incompatible
+        with ``km_curves`` and ``hazard_estimate``. Results are assembled with
+        ``SEQuential.end_of_followup()``.
+    :param end_of_fup_time: The follow-up time ``k`` (in follow-up periods since trial
+        enrollment) at which the end-of-follow-up outcome is evaluated. Required when
+        ``end_of_fup=True``
+    :param end_of_fup_type: Type of end-of-follow-up outcome, either ``'binary'`` (the
+        default, giving the weighted proportion in each arm) or ``'continuous'``
+        (giving the weighted mean)
+    :param end_of_fup_window: Half-width of the window used when a trial-period has no
+        outcome measurement at exactly ``end_of_fup_time``, default ``0`` (no window).
+        Those trial-periods fall back to the measurement nearest to ``k`` within
+        ``[k - window, k + window]`` (ties — measurements equally far either side of
+        ``k`` — are broken toward the later measurement, so that at least ``k`` of
+        follow-up has elapsed); any with no measurement anywhere in the window are
+        censored, i.e. excluded from the average
     :param excused: Boolean to allow excused conditions when method is censoring
     :param excused_colnames: Column names (at the same length of treatment_level) specifying excused conditions, default ``[]``
     :param expand_only: If True, ``SEQuential.expand()`` returns the expanded dataset and skips weighting,
@@ -91,6 +111,10 @@ class SEQopts:
     covariates: Optional[str] = None
     cox_package: Literal["lifelines", "scikit-survival"] = "lifelines"
     denominator: Optional[Union[str, List[str]]] = None
+    end_of_fup: bool = False
+    end_of_fup_time: Optional[float] = None
+    end_of_fup_type: Literal["binary", "continuous"] = "binary"
+    end_of_fup_window: float = 0.0
     excused: bool = False
     excused_colnames: List[str] = field(default_factory=lambda: [])
     expand_only: bool = False
@@ -137,6 +161,7 @@ class SEQopts:
 
     def _validate_bools(self):
         bools = [
+            "end_of_fup",
             "excused",
             "expand_only",
             "followup_class",
@@ -186,6 +211,23 @@ def _validate_ranges(self):
             if any(not isinstance(t, (int, float)) or t < 0 for t in times):
                 raise ValueError("risk_times values must be non-negative numbers.")
 
+        if self.end_of_fup:
+            if self.end_of_fup_time is None or not isinstance(
+                self.end_of_fup_time, (int, float)
+            ):
+                raise ValueError(
+                    "end_of_fup_time must be a single non-missing follow-up time "
+                    "when end_of_fup=True."
+                )
+            if self.end_of_fup_time < 0:
+                raise ValueError("end_of_fup_time must be non-negative.")
+            if not isinstance(self.end_of_fup_window, (int, float)) or (
+                self.end_of_fup_window < 0
+            ):
+                raise ValueError(
+                    "end_of_fup_window must be a single non-negative number."
+                )
+
     def _validate_choices(self):
         if self.plot_type not in ["risk", "survival", "incidence"]:
             raise ValueError(
@@ -197,6 +239,8 @@ def _validate_choices(self):
             raise ValueError("glm_package must be 'statsmodels', 'glum', or 'jax'")
         if self.cox_package not in ["lifelines", "scikit-survival"]:
             raise ValueError("cox_package must be 'lifelines' or 'scikit-survival'")
+        if self.end_of_fup_type not in ["binary", "continuous"]:
+            raise ValueError("end_of_fup_type must be 'binary' or 'continuous'")
 
     def _normalize_formulas(self):
         for i in (

pySEQTarget/SEQoutput.py

@@ -39,6 +39,18 @@ class SEQoutput:
     :type risk_ratio: pl.DataFrame or None
     :param risk_difference: Dataframe of risk differences, compared between treatments and subgroups
     :type risk_difference: pl.DataFrame or None
+    :param eof_data: Per-arm end-of-follow-up estimates (``end_of_fup=True``): the
+        weighted proportion (binary) or mean (continuous) read at
+        ``end_of_fup_time``, with bootstrap confidence intervals when
+        bootstrapped, the eligible trial-periods partitioned into analysed,
+        censored (measured, but not within the window) and never measured, the
+        censoring share, and the distinct contributing subjects
+    :type eof_data: pl.DataFrame or None
+    :param eof_comparison: Pairwise between-arm end-of-follow-up contrasts: the
+        difference in proportions/means with its bootstrap SE and confidence
+        interval (paired by iteration), plus — for a binary outcome only — the
+        ratio of proportions with a log-scale interval and ``log(Ratio) SE``
+    :type eof_comparison: pl.DataFrame or None
     :param time: Timings for every step of the process completed thus far
     :type time: dict or None
     :param diagnostic_tables: Diagnostic tables (outcome, follow-up, switch, and
@@ -63,6 +75,8 @@ class SEQoutput:
     km_graph: matplotlib.figure.Figure = None
     risk_ratio: pl.DataFrame = None
     risk_difference: pl.DataFrame = None
+    eof_data: pl.DataFrame = None
+    eof_comparison: pl.DataFrame = None
     time: dict = None
     diagnostic_tables: dict = None
 
@@ -132,6 +146,11 @@ def retrieve_data(
                 "nonunique_compevent",
                 "unique_switches",
                 "nonunique_switches",
+                "eof_data",
+                "eof_comparison",
+                "unique_eof",
+                "nonunique_eof",
+                "eof_summary",
             ]
         ] = None,
     ) -> pl.DataFrame:
@@ -150,12 +169,25 @@ def retrieve_data(
           (expanded rows). The nonunique count is much larger because each subject
           contributes one row per follow-up period; it is the denominator that,
           with ``nonunique_outcomes``, gives the per-arm event rate.
+        - ``unique_eof`` / ``nonunique_eof`` (``end_of_fup=True`` only): account
+          for every trial-period at the end-of-follow-up time across four
+          mutually exclusive categories — measured ``At k``, measured
+          ``In window``, ``Excluded (outside window)`` and ``Excluded (no
+          measurement)`` — against the ``Eligible`` total. The nonunique
+          (trial-period) counts partition ``Eligible``; the unique (subject)
+          counts may overlap, since one subject can fall into different
+          categories for different trials.
+        - ``eof_summary`` (continuous ``end_of_fup`` only): N/Mean/SD of the
+          analysed measurements per arm, standing in for the suppressed outcome
+          count tables.
 
         :param type: Data which you would like to access, ['km_data', 'hazard',
             'risk_ratio', 'risk_difference', 'unique_outcomes',
             'nonunique_outcomes', 'unique_followup', 'nonunique_followup',
             'unique_compevent', 'nonunique_compevent',
-            'unique_switches', 'nonunique_switches']
+            'unique_switches', 'nonunique_switches',
+            'eof_data', 'eof_comparison', 'unique_eof', 'nonunique_eof',
+            'eof_summary']
         :type type: str
         """
         match type:
@@ -166,9 +198,10 @@ def retrieve_data(
             case "risk_difference":
                 data = self.risk_difference
             case "unique_outcomes":
-                data = self.diagnostic_tables["unique_outcomes"]
+                # Absent for continuous end-of-follow-up outcomes
+                data = self.diagnostic_tables.get("unique_outcomes")
             case "nonunique_outcomes":
-                data = self.diagnostic_tables["nonunique_outcomes"]
+                data = self.diagnostic_tables.get("nonunique_outcomes")
             case "unique_followup":
                 data = self.diagnostic_tables["unique_followup"]
             case "nonunique_followup":
@@ -181,6 +214,16 @@ def retrieve_data(
                 data = self.diagnostic_tables.get("unique_switches")
             case "nonunique_switches":
                 data = self.diagnostic_tables.get("nonunique_switches")
+            case "eof_data":
+                data = self.eof_data
+            case "eof_comparison":
+                data = self.eof_comparison
+            case "unique_eof":
+                data = self.diagnostic_tables.get("unique_eof")
+            case "nonunique_eof":
+                data = self.diagnostic_tables.get("nonunique_eof")
+            case "eof_summary":
+                data = self.diagnostic_tables.get("eof_summary")
             case _:
                 data = self.km_data
         if data is None:

pySEQTarget/SEQuential.py

@@ -9,8 +9,9 @@
 import polars as pl
 
 from .analysis import (_calculate_hazard, _calculate_survival, _clamp,
-                       _outcome_fit, _pred_risk, _risk_estimates,
-                       _subgroup_fit)
+                       _create_endoffup, _eof_counts, _eof_estimate,
+                       _eof_frame, _eof_summary, _outcome_fit, _pred_risk,
+                       _risk_estimates, _subgroup_fit)
 from .error import _data_checker, _param_checker
 from .expansion import _binder, _diagnostics, _dynamic, _random_selection
 from .helpers import Offloader, _col_string, _format_time, bootstrap_loop
@@ -322,6 +323,24 @@ def fit(self) -> None:
         is_boot = boot_idx is not None
         start = getattr(self, "_outcome_start_params", None) if is_boot else None
 
+        if self.end_of_fup:
+            # An end-of-follow-up outcome is a weighted average read at a single
+            # follow-up time, not a fitted model — skip the outcome model.
+            est = _eof_estimate(self)
+            if not is_boot:
+                # Counted from the same data as the estimate so the two always
+                # reconcile; only the main fit populates the diagnostics.
+                DT_eof = _eof_frame(self)
+                diag = dict(getattr(self, "diagnostics", None) or {})
+                diag["unique_eof"] = _eof_counts(self, DT_eof, unique=True)
+                diag["nonunique_eof"] = _eof_counts(self, DT_eof, unique=False)
+                if self.end_of_fup_type == "continuous":
+                    # Mean/SD of the analysed measurements stands in for the
+                    # (suppressed) outcome count tables for continuous outcomes
+                    diag["eof_summary"] = _eof_summary(self, DT_eof)
+                self.diagnostics = diag
+            return {"eof": est}
+
         if self.subgroup_colname is not None:
             models_list = _subgroup_fit(self, start_params=start)
             if not is_boot:
@@ -383,6 +402,12 @@ def survival(self, **kwargs) -> None:
             else:
                 raise ValueError(f"Unknown or misplaced argument: {key}")
 
+        if self.end_of_fup:
+            raise ValueError(
+                "Survival curves are not available with end_of_fup=True; use "
+                "the 'end_of_followup' method instead."
+            )
+
         if not hasattr(self, "outcome_model") or not self.outcome_model:
             raise ValueError(
                 "Outcome model not found. Please run the 'fit' method before calculating survival."
@@ -401,12 +426,47 @@ def survival(self, **kwargs) -> None:
         end = time.perf_counter()
         self._survival_time = _format_time(start, end)
 
+    def end_of_followup(self) -> None:
+        """
+        Assembles the end-of-follow-up estimates (``end_of_fup=True``):
+        the per-arm weighted proportion or mean read at ``end_of_fup_time``
+        (``eof_data``) and the pairwise between-arm contrasts
+        (``eof_comparison``), with bootstrap confidence intervals when
+        bootstrapped. Contrasts are paired by bootstrap iteration, so the
+        interval accounts for the correlation between arms.
+        """
+        start = time.perf_counter()
+
+        if not self.end_of_fup:
+            raise ValueError(
+                "End-of-follow-up estimates were not created as a result of "
+                "end_of_fup=False."
+            )
+        if not hasattr(self, "outcome_model") or not self.outcome_model:
+            raise ValueError(
+                "End-of-follow-up estimates not found. Please run the 'fit' "
+                "method before assembling them."
+            )
+
+        eof = _create_endoffup(self)
+        self.eof_data = eof["eof_data"]
+        self.eof_comparison = eof["eof_comparison"]
+
+        end = time.perf_counter()
+        self._eof_time = _format_time(start, end)
+
     def hazard(self) -> None:
         """
         Uses fit outcome models (outcome, competing event) to estimate hazard ratios
         """
         start = time.perf_counter()
 
+        if self.end_of_fup:
+            raise ValueError(
+                "Hazard ratios are not available with end_of_fup=True; use "
+                "the 'end_of_followup' method instead."
+            )
+
         if self.method == "dose-response":
             raise NotImplementedError(
                 "Hazard ratio estimation is not supported for method='dose-response': "
@@ -457,6 +517,8 @@ def collect(self) -> SEQoutput:
             "_model_time",
             "_expansion_time",
             "weight_stats",
+            "eof_data",
+            "eof_comparison",
         ]
         for attr in generated:
             if not hasattr(self, attr):
@@ -479,7 +541,11 @@ def collect(self) -> SEQoutput:
             "collection_time": self._time_collected,
         }
 
-        if self.outcome_model is not None:
+        if self.end_of_fup:
+            # No outcome model is fit for an end-of-follow-up outcome
+            outcome_models = None
+            compevent_models = None
+        elif self.outcome_model is not None:
             outcome_models = [model["outcome"] for model in self.outcome_model]
             if self.compevent_colname is not None:
                 compevent_models = [model["compevent"] for model in self.outcome_model]
@@ -510,6 +576,8 @@ def collect(self) -> SEQoutput:
             km_graph=self.km_graph,
             risk_ratio=risk_ratio,
             risk_difference=risk_difference,
+            eof_data=self.eof_data,
+            eof_comparison=self.eof_comparison,
             time=time,
             diagnostic_tables=self.diagnostics,
         )

pySEQTarget/analysis/__init__.py

@@ -1,3 +1,5 @@
+from ._endoffup import (_create_endoffup, _eof_counts, _eof_estimate,
+                        _eof_frame, _eof_summary)
 from ._hazard import _calculate_hazard
 from ._outcome_fit import _outcome_fit
 from ._risk_estimates import _risk_estimates
@@ -7,6 +9,11 @@
 
 __all__ = [
     "_calculate_hazard",
+    "_create_endoffup",
+    "_eof_counts",
+    "_eof_estimate",
+    "_eof_frame",
+    "_eof_summary",
     "_outcome_fit",
     "_risk_estimates",
     "_subgroup_fit",