From 9c5917662beab047a502de718b90f3d74a921c6b Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Fri, 14 Aug 2026 13:18:51 +0100 Subject: [PATCH 1/9] Bump astral-sh/setup-uv to v10.0.1 --- .github/workflows/python-app.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/python-app.yml b/.github/workflows/python-app.yml index 2fd9136..f104b35 100644 --- a/.github/workflows/python-app.yml +++ b/.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 }} From 8307262ef30474d4c72d0598a7aeb0d9a58566e2 Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Fri, 28 Aug 2026 08:34:43 +0100 Subject: [PATCH 2/9] Add end-of-follow-up outcomes Single-time measurements averaged per arm with window fallback, paired bootstrap contrasts, accounting diagnostics, and a new end_of_followup method. --- pySEQTarget/SEQopts.py | 44 +++ pySEQTarget/SEQoutput.py | 49 ++- pySEQTarget/SEQuential.py | 74 ++++- pySEQTarget/analysis/__init__.py | 7 + pySEQTarget/analysis/_endoffup.py | 374 ++++++++++++++++++++++ pySEQTarget/error/_data_checker.py | 28 +- pySEQTarget/error/_param_checker.py | 33 ++ pySEQTarget/expansion/_binder.py | 19 +- pySEQTarget/expansion/_diagnostics.py | 9 +- pySEQTarget/helpers/_output_files.py | 11 + tests/test_end_of_fup.py | 441 ++++++++++++++++++++++++++ 11 files changed, 1070 insertions(+), 19 deletions(-) create mode 100644 pySEQTarget/analysis/_endoffup.py create mode 100644 tests/test_end_of_fup.py diff --git a/pySEQTarget/SEQopts.py b/pySEQTarget/SEQopts.py index cf63e76..e6f542f 100644 --- a/pySEQTarget/SEQopts.py +++ b/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 ( diff --git a/pySEQTarget/SEQoutput.py b/pySEQTarget/SEQoutput.py index bcce6a7..1d8ab5a 100644 --- a/pySEQTarget/SEQoutput.py +++ b/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: diff --git a/pySEQTarget/SEQuential.py b/pySEQTarget/SEQuential.py index c120b93..eafb5d3 100644 --- a/pySEQTarget/SEQuential.py +++ b/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, ) diff --git a/pySEQTarget/analysis/__init__.py b/pySEQTarget/analysis/__init__.py index f39b8c1..7e7097c 100644 --- a/pySEQTarget/analysis/__init__.py +++ b/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", diff --git a/pySEQTarget/analysis/_endoffup.py b/pySEQTarget/analysis/_endoffup.py new file mode 100644 index 0000000..0a3d431 --- /dev/null +++ b/pySEQTarget/analysis/_endoffup.py @@ -0,0 +1,374 @@ +"""End-of-follow-up outcomes. + +An end-of-follow-up outcome is measured once, at follow-up time +``end_of_fup_time`` (``k``), rather than as a time-to-event. Instead of fitting +a survival outcome model, the estimate in each baseline treatment arm is the +weighted average of the outcome read at ``k`` — weighted by the +period-trial-specific weight at the time the measurement was taken. Ported from +SEQTaRget's ``internal_endoffup.R``. +""" + +import numpy as np +import polars as pl +from scipy import stats + +from ._risk_estimates import _ci_label + + +def _eof_frame(self): + """The frame the estimate is read from. + + Under ``method='censoring'`` the artificially censored (treatment switch) + rows are not measurements, so a subject who deviates before ``k`` is + correctly excluded rather than contributing a carried-forward value. + """ + DT = self.DT + if self.method == "censoring" and "switch" in DT.columns: + DT = DT.filter(pl.col("switch") != 1) + return DT + + +def _eof_group_cols(self, DT): + tx_bas = f"{self.treatment_col}{self.indicator_baseline}" + cols = [tx_bas] + if self.subgroup_colname is not None and self.subgroup_colname in DT.columns: + cols.append(self.subgroup_colname) + return cols + + +def _eof_measure(self, DT): + """Select the end-of-follow-up measurement for each (id, trial). + + The measurement at exactly ``k`` when one exists, otherwise — if + ``end_of_fup_window`` is non-zero — the measurement nearest to ``k`` within + ``[k - window, k + window]``, with ties (measurements equally far either + side of ``k``) broken toward the later one, so that at least ``k`` of + follow-up has elapsed. Trial-periods with no measurement anywhere in the + window contribute no row, i.e. they are censored out of the estimate. + """ + k = self.end_of_fup_time + w = self.end_of_fup_window + tx_bas = f"{self.treatment_col}{self.indicator_baseline}" + + cols = [self.id_col, "trial", "followup", tx_bas, self.outcome_col] + if self.subgroup_colname is not None and self.subgroup_colname in DT.columns: + cols.append(self.subgroup_colname) + if "weight" in DT.columns: + cols.append("weight") + + candidates = DT.filter( + pl.col(self.outcome_col).is_not_null() + & (pl.col("followup") >= k - w) + & (pl.col("followup") <= k + w) + ).select(cols) + + # Nearest to k wins (exact k has distance 0); equidistant ties break toward + # the later measurement so at least k of follow-up has elapsed. followup is + # UInt32 in the expanded frame, so cast before differencing — unsigned + # subtraction wraps for followups below k and .abs() cannot recover it. + measured = ( + candidates.with_columns( + (pl.col("followup").cast(pl.Int64) - k).abs().alias("_dist") + ) + .sort( + [self.id_col, "trial", "_dist", "followup"], + descending=[False, False, False, True], + ) + .group_by([self.id_col, "trial"], maintain_order=True) + .first() + .drop("_dist") + .rename({self.outcome_col: "eof_value"}) + ) + + # Unweighted analysis = equal weights + if "weight" not in measured.columns: + measured = measured.with_columns(pl.lit(1.0).alias("weight")) + return measured + + +def _eof_period_flags(self, DT): + """One row per trial-period, flagged with what it has available.""" + k = self.end_of_fup_time + w = self.end_of_fup_window + valid = pl.col(self.outcome_col).is_not_null() + in_window = valid & (pl.col("followup") >= k - w) & (pl.col("followup") <= k + w) + + return DT.group_by( + [self.id_col, "trial"] + _eof_group_cols(self, DT) + ).agg( + [ + (valid & (pl.col("followup") == k)).any().alias("at_k"), + in_window.any().alias("in_window"), + valid.any().alias("measured"), + ] + ) + + +def _eof_estimate(self): + """Weighted end-of-follow-up average within each treatment arm. + + Weight truncation (``weight_min`` / ``weight_max``, including the bounds + ``weight_p99`` resolves to) is applied here as it is for the outcome model, + since this average is the estimator in end_of_fup mode. + + Alongside the estimate this counts the trial-periods censored for want of a + measurement in the window — those measured at some point but not within + ``[k - window, k + window]`` — so that share can be reported next to the + estimate they were dropped from. Trial-periods never measured at all are + counted separately rather than folded in, so the analysed, censored and + never-measured counts partition the eligible total. + """ + DT = _eof_frame(self) + by = _eof_group_cols(self, DT) + + measured = _eof_measure(self, DT).with_columns( + pl.col("weight").clip( + lower_bound=self.weight_min, upper_bound=self.weight_max + ) + ) + + est = measured.group_by(by).agg( + [ + ((pl.col("weight") * pl.col("eof_value")).sum() / pl.col("weight").sum()) + .alias("estimate"), + pl.len().alias("n"), + pl.col(self.id_col).n_unique().alias("n_subjects"), + ] + ) + totals = _eof_period_flags(self, DT).group_by(by).agg( + [ + pl.len().alias("n_eligible"), + (pl.col("measured") & ~pl.col("in_window")).sum().alias("n_censored"), + (~pl.col("measured")).sum().alias("n_nomeasure"), + ] + ) + return est.join(totals, on=by, how="inner").sort(by) + + +def _eof_counts(self, DT, unique): + """Account for every trial-period at the end-of-follow-up time. + + Four mutually exclusive categories against the Eligible total: measured + ``At k``; measured ``In window`` (no measurement at k but one within the + window); ``Excluded (outside window)`` (measured somewhere, but not within + the window); ``Excluded (no measurement)``. Trial-period counts partition + Eligible; subject counts (``unique=True``) need not, since one subject can + fall into different categories for different trials. + """ + by = _eof_group_cols(self, DT) + levels = ["At k", "In window", "Excluded (outside window)", "Excluded (no measurement)"] + + flags = _eof_period_flags(self, DT).with_columns( + pl.when(pl.col("at_k")) + .then(pl.lit(levels[0])) + .when(pl.col("in_window")) + .then(pl.lit(levels[1])) + .when(pl.col("measured")) + .then(pl.lit(levels[2])) + .otherwise(pl.lit(levels[3])) + .alias("_category") + ) + + counter = ( + pl.col(self.id_col).n_unique() if unique else pl.len() + ) + counted = flags.group_by(by + ["_category"]).agg(counter.alias("N")) + wide = counted.pivot(on="_category", index=by, values="N").fill_null(0) + for lv in levels: + if lv not in wide.columns: + wide = wide.with_columns(pl.lit(0).alias(lv)) + totals = flags.group_by(by).agg(counter.alias("Eligible")) + return wide.join(totals, on=by).select(by + ["Eligible"] + levels).sort(by) + + +def _eof_summary(self, DT): + """N, mean and SD of the raw selected measurements per baseline arm — the + unweighted analogue of the outcome count tables, reported for continuous + outcomes where event counts have no meaning.""" + by = _eof_group_cols(self, DT) + measured = _eof_measure(self, DT) + return ( + measured.group_by(by) + .agg( + [ + pl.len().alias("N"), + pl.col("eof_value").mean().alias("Mean"), + pl.col("eof_value").std().alias("SD"), + ] + ) + .sort(by) + .rename({by[0]: "A"}) + ) + + +def _create_endoffup(self): + """Assemble end-of-follow-up estimates and bootstrap confidence intervals. + + Mirrors the risk assembly: ``eof_data`` holds the per-arm estimate and + ``eof_comparison`` the pairwise between-arm contrasts, both with bootstrap + confidence intervals when bootstrapped. Contrasts are paired by bootstrap + iteration, so the interval accounts for the correlation between arms. + """ + tx_bas = f"{self.treatment_col}{self.indicator_baseline}" + sub = self.subgroup_colname + is_binary = self.end_of_fup_type == "binary" + label = "Proportion" if is_binary else "Mean" + ci_label = _ci_label(self.bootstrap_CI) + z = stats.norm.ppf(1 - (1 - self.bootstrap_CI) / 2) + alpha = (1 - self.bootstrap_CI) / 2 + use_se = self.bootstrap_CI_method == "se" + + full = self.outcome_model[0]["eof"] + boots = [ + m["eof"] for m in self.outcome_model[1:] if m["eof"] is not None and m["eof"].height > 0 + ] + has_ci = len(boots) > 1 + sub_in_full = sub is not None and sub in full.columns + key_cols = ([sub] if sub_in_full else []) + [tx_bas] + + def boot_estimates(key): + """Per-iteration estimates for one (subgroup, arm) key, NaN when absent.""" + out = np.full(len(boots), np.nan) + for i, b in enumerate(boots): + f = b + for col, val in zip(key_cols, key): + f = f.filter(pl.col(col) == val) + if f.height == 1: + out[i] = f["estimate"][0] + return out + + keys = [tuple(row) for row in full.select(key_cols).iter_rows()] + + # Per-arm table ========================================================== + data = full.rename( + { + tx_bas: "A", + "estimate": label, + "n_eligible": "Trial-periods (Eligible)", + "n": "Trial-periods (Analysed)", + "n_censored": "Trial-periods (Censored)", + "n_nomeasure": "Trial-periods (No measurement)", + "n_subjects": "Subjects", + } + ).with_columns( + [ + pl.lit(self.end_of_fup_type).alias("Type"), + pl.lit(float(self.end_of_fup_time)).alias("Time"), + ( + 100 + * pl.col("Trial-periods (Censored)") + / pl.col("Trial-periods (Eligible)") + ).alias("% Censored"), + ] + ) + + if has_ci: + se_rows, lci_rows, uci_rows = [], [], [] + for key in keys: + draws = boot_estimates(key) + se = float(np.nanstd(draws, ddof=1)) + point = float( + full.filter( + pl.all_horizontal( + [pl.col(c) == v for c, v in zip(key_cols, key)] + ) + )["estimate"][0] + ) + if use_se: + lci, uci = point - z * se, point + z * se + else: + lci = float(np.nanquantile(draws, alpha)) + uci = float(np.nanquantile(draws, 1 - alpha)) + if is_binary: + # Proportions are bounded, so clamp the binary interval to [0, 1] + lci, uci = max(0.0, lci), min(1.0, uci) + se_rows.append(se) + lci_rows.append(lci) + uci_rows.append(uci) + ci_frame = pl.DataFrame( + { + **{c: [k[i] for k in keys] for i, c in enumerate(key_cols)}, + "SE": se_rows, + f"{ci_label} LCI": lci_rows, + f"{ci_label} UCI": uci_rows, + } + ) + rename_back = {key_cols[-1]: "A"} + data = data.join(ci_frame.rename(rename_back), on=(["A"] if not sub_in_full else [sub, "A"])) + + lead = ["Type", "Time"] + ([sub] if sub_in_full else []) + [ + "A", + label, + "Trial-periods (Eligible)", + "Trial-periods (Analysed)", + "Trial-periods (Censored)", + "Trial-periods (No measurement)", + "% Censored", + "Subjects", + ] + data = data.select(lead + [c for c in data.columns if c not in lead]) + + # Pairwise between-arm contrasts, both directions ======================== + rows = [] + subgroup_vals = ( + sorted(set(full[sub].to_list())) if sub_in_full else [None] + ) + for g in subgroup_vals: + sub_frame = full.filter(pl.col(sub) == g) if g is not None else full + arms = sub_frame[tx_bas].to_list() + point = dict(zip(arms, sub_frame["estimate"].to_list())) + for a_x in arms: + for a_y in arms: + if a_x == a_y: + continue + row = {"Time": float(self.end_of_fup_time)} + if g is not None: + row[sub] = g + row["A_x"] = a_x + row["A_y"] = a_y + diff = point[a_y] - point[a_x] + ratio = point[a_y] / point[a_x] if point[a_x] != 0 else float("inf") + row["Difference"] = diff + if has_ci: + kx = ((g,) if g is not None else ()) + (a_x,) + ky = ((g,) if g is not None else ()) + (a_y,) + d = boot_estimates(ky) - boot_estimates(kx) + d_se = float(np.nanstd(d, ddof=1)) + with np.errstate(divide="ignore", invalid="ignore"): + r = boot_estimates(ky) / boot_estimates(kx) + r_valid = r[np.isfinite(r) & (r > 0)] + r_logse = ( + float(np.std(np.log(r_valid), ddof=1)) + if r_valid.size > 1 + else float("nan") + ) + if use_se: + d_lci, d_uci = diff - z * d_se, diff + z * d_se + if np.isfinite(r_logse) and ratio > 0 and np.isfinite(ratio): + r_lci = float(np.exp(np.log(ratio) - z * r_logse)) + r_uci = float(np.exp(np.log(ratio) + z * r_logse)) + else: + r_lci = r_uci = float("nan") + else: + d_lci = float(np.nanquantile(d, alpha)) + d_uci = float(np.nanquantile(d, 1 - alpha)) + if r_valid.size > 1: + r_lci = float(np.quantile(r_valid, alpha)) + r_uci = float(np.quantile(r_valid, 1 - alpha)) + else: + r_lci = r_uci = float("nan") + row[f"Difference {ci_label} LCI"] = d_lci + row[f"Difference {ci_label} UCI"] = d_uci + row["Difference SE"] = d_se + # Ratios need an outcome bounded away from zero, so they are + # reported for proportions only + if is_binary: + row["Ratio"] = ratio + if has_ci: + row[f"Ratio {ci_label} LCI"] = r_lci + row[f"Ratio {ci_label} UCI"] = r_uci + row["log(Ratio) SE"] = r_logse + rows.append(row) + + comparison = pl.DataFrame(rows) if rows else pl.DataFrame() + return {"eof_data": data, "eof_comparison": comparison} diff --git a/pySEQTarget/error/_data_checker.py b/pySEQTarget/error/_data_checker.py index ec6bc8a..f71da0a 100644 --- a/pySEQTarget/error/_data_checker.py +++ b/pySEQTarget/error/_data_checker.py @@ -4,14 +4,38 @@ def _check_binary(data, col): unique_vals = set(data[col].drop_nulls().unique().to_list()) if not unique_vals.issubset({0, 1}): + # Cap the listing — a continuous column can hold thousands of values + offending = sorted(unique_vals - {0, 1}) + shown = ", ".join(str(v) for v in offending[:10]) + if len(offending) > 10: + shown += f", ... ({len(offending)} distinct values)" raise ValueError( - f"Column '{col}' must be binary (0/1) but contains values: {sorted(unique_vals)}" + f"Column '{col}' must be binary (0/1) but contains values: {shown}" ) def _data_checker(self): _check_binary(self.data, self.eligible_col) - _check_binary(self.data, self.outcome_col) + + # end_of_fup treats a null outcome as "not measured at this time"; in every + # other mode the outcome must be complete. + if not self.end_of_fup and self.data[self.outcome_col].null_count() > 0: + raise ValueError( + f"Column '{self.outcome_col}' contains missing values; missing " + "outcome measurements are only permitted with end_of_fup=True." + ) + # Continuous end-of-follow-up outcomes are averaged, not modelled, so + # non-binary values are expected there. + if not (self.end_of_fup and self.end_of_fup_type == "continuous"): + try: + _check_binary(self.data, self.outcome_col) + except ValueError as e: + if self.end_of_fup: + raise ValueError( + f"{e} For an outcome that is not 0/1, set " + "end_of_fup_type='continuous' in SEQopts." + ) from None + raise if self.cense_eligible_colname is not None: _check_binary(self.data, self.cense_eligible_colname) diff --git a/pySEQTarget/error/_param_checker.py b/pySEQTarget/error/_param_checker.py index 941cd85..0afe9c6 100644 --- a/pySEQTarget/error/_param_checker.py +++ b/pySEQTarget/error/_param_checker.py @@ -41,6 +41,39 @@ def _param_checker(self): if self.km_curves and self.hazard_estimate: raise ValueError("km_curves and hazard cannot both be set to True.") + # End-of-follow-up outcomes replace the survival outcome model with a direct + # weighted average at a single follow-up time, so the survival-based outputs + # have no meaning there and the requested time must lie inside the expansion. + if self.end_of_fup: + if self.km_curves or self.hazard_estimate: + raise ValueError( + "end_of_fup is not compatible with km_curves or hazard_estimate: " + "an end-of-follow-up outcome is evaluated at a single time, so " + "there is no survival curve or hazard to estimate." + ) + if self.method == "dose-response": + raise ValueError( + "end_of_fup is not supported for the dose-response method." + ) + if self.compevent_colname is not None: + raise ValueError( + "end_of_fup is not compatible with compevent_colname: competing " + "events are a survival-outcome concept." + ) + upper = self.end_of_fup_time + self.end_of_fup_window + if upper > self.followup_max: + raise ValueError( + f"end_of_fup_time plus end_of_fup_window ({upper}) exceeds the " + f"maximum follow-up ({self.followup_max}); widen followup_max or " + "narrow the window." + ) + if self.end_of_fup_time - self.end_of_fup_window < self.followup_min: + raise ValueError( + "end_of_fup_time minus end_of_fup_window " + f"({self.end_of_fup_time - self.end_of_fup_window}) is below the " + f"minimum follow-up ({self.followup_min})." + ) + if self.hazard_estimate and self.method == "dose-response": raise ValueError( "Hazard ratio estimation is not supported for method='dose-response': " diff --git a/pySEQTarget/expansion/_binder.py b/pySEQTarget/expansion/_binder.py index 4e10271..3c58e40 100644 --- a/pySEQTarget/expansion/_binder.py +++ b/pySEQTarget/expansion/_binder.py @@ -101,13 +101,16 @@ def _binder(self, kept_cols): # Truncate each (id, trial) at the first outcome event so that subjects who # experience the outcome early are not carried forward with subsequent rows. - DT = DT.filter( - pl.col(self.outcome_col) - .fill_null(0) - .cum_max() - .shift(1, fill_value=0) - .over([self.id_col, "trial"]) - == 0 - ) + # end_of_fup outcomes are measurements, not events — truncating would + # discard the rows the estimate is read from. + if not self.end_of_fup: + DT = DT.filter( + pl.col(self.outcome_col) + .fill_null(0) + .cum_max() + .shift(1, fill_value=0) + .over([self.id_col, "trial"]) + == 0 + ) return DT diff --git a/pySEQTarget/expansion/_diagnostics.py b/pySEQTarget/expansion/_diagnostics.py index 56dcb0c..bd4dcf3 100644 --- a/pySEQTarget/expansion/_diagnostics.py +++ b/pySEQTarget/expansion/_diagnostics.py @@ -2,9 +2,12 @@ def _diagnostics(self): - unique_out = _outcome_diag(self, unique=True) - nonunique_out = _outcome_diag(self, unique=False) - out = {"unique_outcomes": unique_out, "nonunique_outcomes": nonunique_out} + # Outcome tables count outcome == 1 rows — meaningless for a continuous + # end-of-follow-up outcome, so they are suppressed there. + out = {} + if not (self.end_of_fup and self.end_of_fup_type == "continuous"): + out["unique_outcomes"] = _outcome_diag(self, unique=True) + out["nonunique_outcomes"] = _outcome_diag(self, unique=False) unique_fu = _followup_diag(self, unique=True) nonunique_fu = _followup_diag(self, unique=False) diff --git a/pySEQTarget/helpers/_output_files.py b/pySEQTarget/helpers/_output_files.py index dd4cbb7..8c85164 100644 --- a/pySEQTarget/helpers/_output_files.py +++ b/pySEQTarget/helpers/_output_files.py @@ -49,6 +49,17 @@ def _model_section(title, kind): _model_section("Outcome Model", "outcome") + if getattr(self.options, "end_of_fup", False) and self.eof_data is not None: + lines.append("### End-of-Follow-up Outcome") + lines.append("") + lines.append(self.eof_data.to_pandas().to_markdown(index=False)) + lines.append("") + if self.eof_comparison is not None and self.eof_comparison.height > 0: + lines.append("#### Between-arm Comparison") + lines.append("") + lines.append(self.eof_comparison.to_pandas().to_markdown(index=False)) + lines.append("") + if self.options.hazard_estimate and self.hazard is not None: lines.append("### Hazard") lines.append("") diff --git a/tests/test_end_of_fup.py b/tests/test_end_of_fup.py new file mode 100644 index 0000000..7946032 --- /dev/null +++ b/tests/test_end_of_fup.py @@ -0,0 +1,441 @@ +"""End-of-follow-up outcomes: a single measurement taken at end_of_fup_time, +averaged within each baseline arm using the weight at that time, rather than a +time-to-event fitted with a survival outcome model. Ported from SEQTaRget's +test_end_of_fup.R (PR #168).""" + +import numpy as np +import polars as pl +import pytest +from pytest import approx + +from pySEQTarget import SEQopts, SEQuential +from pySEQTarget.data import load_data +from pySEQTarget.analysis._endoffup import _eof_measure + + +def _eof_run(data=None, outcome="outcome", method="ITT", run=True, **opts): + s = SEQuential( + data if data is not None else load_data("SEQdata"), + id_col="ID", + time_col="time", + eligible_col="eligible", + treatment_col="tx_init", + outcome_col=outcome, + time_varying_cols=["N", "L", "P"], + fixed_cols=["sex"], + method=method, + parameters=SEQopts(seed=42, **opts), + ) + if run: + s.expand() + if s.bootstrap_nboot > 0: + s.bootstrap() + s.fit() + s.end_of_followup() + return s + + +def _continuous_data(seed=42): + rng = np.random.default_rng(seed) + d = load_data("SEQdata") + return d.with_columns( + ( + 10 + + 2 * pl.col("tx_init") + + pl.col("N") + + pl.Series(rng.standard_normal(d.height)) + ).alias("cont") + ) + + +def test_unweighted_estimate_is_mean_of_selected_measurements(): + k, w = 12, 3 + s = _eof_run(end_of_fup=True, end_of_fup_time=k, end_of_fup_window=w) + + # Independent re-implementation of the selection rule: the non-missing value + # nearest to k within [k - w, k + w], ties broken toward the later one. + manual = ( + s.DT.filter( + pl.col("outcome").is_not_null() + & (pl.col("followup") >= k - w) + & (pl.col("followup") <= k + w) + ) + .with_columns((pl.col("followup").cast(pl.Int64) - k).abs().alias("d")) + .sort(["ID", "trial", "d", "followup"], descending=[False, False, False, True]) + .group_by(["ID", "trial"], maintain_order=True) + .first() + .group_by("tx_init_bas") + .agg([pl.col("outcome").mean().alias("manual"), pl.len().alias("n")]) + .sort("tx_init_bas") + ) + + got = s.eof_data.sort("A") + assert got["Proportion"].to_list() == approx(manual["manual"].to_list()) + assert got["Trial-periods (Analysed)"].to_list() == manual["n"].to_list() + + +def test_window_only_adds_trial_periods_without_measurement_at_k(): + exact = _eof_run(end_of_fup=True, end_of_fup_time=12) + windowed = _eof_run(end_of_fup=True, end_of_fup_time=12, end_of_fup_window=3) + + n_exact = exact.eof_data.sort("A")["Trial-periods (Analysed)"].to_list() + n_win = windowed.eof_data.sort("A")["Trial-periods (Analysed)"].to_list() + # Widening the window can only ever add contributors, never remove them + assert all(w >= e for w, e in zip(n_win, n_exact)) + assert any(w > e for w, e in zip(n_win, n_exact)) + + # A window of 0 is the same as no window at all + zero = _eof_run(end_of_fup=True, end_of_fup_time=12, end_of_fup_window=0) + assert zero.eof_data.equals(exact.eof_data) + + +def test_window_takes_nearest_measurement_not_earliest(): + # Hand-built trial-periods where 'nearest' and 'earliest' disagree, driven + # through the selection helper directly so the rule is tested in isolation. + s = _eof_run(end_of_fup=True, end_of_fup_time=3, end_of_fup_window=2, run=False) + + DT = pl.DataFrame( + { + "ID": [1, 1, 2, 2, 3, 3], + "trial": [0] * 6, + "followup": [1, 4, 2, 3, 2, 4], + "tx_init_bas": [0, 0, 1, 1, 0, 0], + "outcome": [10.0, 20.0, 30.0, 40.0, 50.0, 60.0], + } + # followup is UInt32 in the real expanded frame; the distance-to-k + # computation must not wrap for followups below k (unsigned subtraction) + ).with_columns(pl.col("followup").cast(pl.UInt32)) + got = _eof_measure(s, DT).sort("ID") + + # ID 1: |1-3|=2 vs |4-3|=1, so the later measurement is nearer — the + # earliest rule would have taken followup 1 + # ID 2: measured at exactly k, which always wins + # ID 3: |2-3|=|4-3|=1, an equidistant tie broken toward the later, so that + # at least k of follow-up has elapsed + assert got["followup"].to_list() == [4, 3, 4] + assert got["eof_value"].to_list() == [20.0, 40.0, 60.0] + + +def test_continuous_outcome_reported_as_mean_not_clamped(): + s = _eof_run( + data=_continuous_data(), + outcome="cont", + end_of_fup=True, + end_of_fup_time=12, + end_of_fup_type="continuous", + end_of_fup_window=3, + bootstrap_nboot=5, + ) + assert "Mean" in s.eof_data.columns + assert "Proportion" not in s.eof_data.columns + assert (s.eof_data["Mean"] > 1).all() + # A continuous mean is unbounded, so its interval must not be clamped to [0, 1] + lci_col = next(c for c in s.eof_data.columns if "LCI" in c) + assert (s.eof_data[lci_col] > 1).all() + + +def test_binary_intervals_clamped_to_unit_range(): + s = _eof_run( + end_of_fup=True, end_of_fup_time=12, end_of_fup_window=3, bootstrap_nboot=5 + ) + lci = next(c for c in s.eof_data.columns if "LCI" in c) + uci = next(c for c in s.eof_data.columns if "UCI" in c) + assert (s.eof_data[lci] >= 0).all() + assert (s.eof_data[uci] <= 1).all() + + +def test_weighted_censoring_end_of_fup_runs(): + s = _eof_run( + method="censoring", + end_of_fup=True, + end_of_fup_time=12, + end_of_fup_window=3, + weighted=True, + weight_preexpansion=True, + ) + assert s.eof_data.height == 2 + assert s.eof_data["Proportion"].is_between(0, 1).all() + # Weighted and unweighted estimates should differ (weights are not all 1) + unw = _eof_run( + method="censoring", end_of_fup=True, end_of_fup_time=12, end_of_fup_window=3 + ) + assert s.eof_data["Proportion"].to_list() != approx( + unw.eof_data["Proportion"].to_list(), rel=1e-12 + ) + + +def test_bootstrap_gives_per_arm_and_paired_ci_with_requested_level(): + s = _eof_run( + end_of_fup=True, + end_of_fup_time=12, + end_of_fup_window=3, + bootstrap_nboot=5, + bootstrap_CI=0.9, + ) + assert "90% LCI" in s.eof_data.columns and "90% UCI" in s.eof_data.columns + assert "SE" in s.eof_data.columns + comp = s.eof_comparison + for col in ("Difference 90% LCI", "Difference 90% UCI", "Difference SE"): + assert col in comp.columns + # Paired contrasts are antisymmetric + comp = comp.sort(["A_x", "A_y"]) + assert comp["Difference"][0] == approx(-comp["Difference"][1]) + + +def test_binary_reports_ratio_continuous_does_not(): + b = _eof_run( + end_of_fup=True, end_of_fup_time=12, end_of_fup_window=3, bootstrap_nboot=5 + ) + assert "Ratio" in b.eof_comparison.columns + assert "log(Ratio) SE" in b.eof_comparison.columns + + c = _eof_run( + data=_continuous_data(), + outcome="cont", + end_of_fup=True, + end_of_fup_time=12, + end_of_fup_type="continuous", + end_of_fup_window=3, + bootstrap_nboot=5, + ) + assert "Difference" in c.eof_comparison.columns + assert "Ratio" not in c.eof_comparison.columns + + +def test_subgroups_produce_one_estimate_set_each(): + s = _eof_run( + end_of_fup=True, + end_of_fup_time=12, + end_of_fup_window=3, + subgroup_colname="sex", + ) + assert "sex" in s.eof_data.columns + # One row per (subgroup, arm) + assert s.eof_data.height == 4 + assert "sex" in s.eof_comparison.columns + assert s.eof_comparison.height == 4 # both directions x two subgroups + + +def test_incompatible_options_raise(): + with pytest.raises(ValueError, match="km_curves or hazard"): + _eof_run(end_of_fup=True, end_of_fup_time=12, km_curves=True, run=False) + with pytest.raises(ValueError, match="km_curves or hazard"): + _eof_run(end_of_fup=True, end_of_fup_time=12, hazard_estimate=True, run=False) + with pytest.raises(ValueError, match="dose-response"): + _eof_run( + end_of_fup=True, end_of_fup_time=12, method="dose-response", run=False + ) + with pytest.raises(ValueError, match="compevent"): + _eof_run( + end_of_fup=True, + end_of_fup_time=12, + compevent_colname="outcome", + run=False, + ) + + +def test_seqopts_validates_end_of_fup_arguments(): + with pytest.raises(ValueError, match="end_of_fup_time"): + SEQopts(end_of_fup=True) + with pytest.raises(ValueError, match="end_of_fup_type"): + SEQopts(end_of_fup=True, end_of_fup_time=12, end_of_fup_type="count") + with pytest.raises(ValueError, match="end_of_fup_window"): + SEQopts(end_of_fup=True, end_of_fup_time=12, end_of_fup_window=-1) + with pytest.raises(ValueError, match="non-negative"): + SEQopts(end_of_fup=True, end_of_fup_time=-3) + + +def test_requested_time_must_lie_within_expanded_followup(): + with pytest.raises(ValueError, match="exceeds the maximum follow-up"): + _eof_run(end_of_fup=True, end_of_fup_time=10_000, run=False) + with pytest.raises(ValueError, match="exceeds the maximum follow-up"): + _eof_run( + end_of_fup=True, + end_of_fup_time=12, + end_of_fup_window=3, + followup_max=13, + run=False, + ) + + +def test_expansion_not_truncated_at_first_event(): + # The survival path cuts each trial at its first outcome row. An + # end-of-follow-up outcome is a status measured repeatedly and read at a + # fixed time, so that truncation must not apply or the measurement at k + # would be discarded for anyone whose status was ever 1 earlier. SEQdata + # cannot show this — every subject's series already ends at their event — + # so this uses a dataset where measurement genuinely continues afterwards. + rng = np.random.default_rng(7) + n_id, n_t = 40, 20 + d = pl.DataFrame( + { + "ID": np.repeat(np.arange(1, n_id + 1), n_t), + "time": np.tile(np.arange(n_t), n_id), + "eligible": np.ones(n_id * n_t, dtype=int), + "tx_init": rng.binomial(1, 0.5, n_id * n_t), + "outcome": rng.binomial(1, 0.3, n_id * n_t), # recurring status + "N": rng.standard_normal(n_id * n_t), + "L": rng.standard_normal(n_id * n_t), + "P": rng.standard_normal(n_id * n_t), + "sex": np.repeat(rng.binomial(1, 0.5, n_id), n_t), + } + ) + + eof = _eof_run(data=d, end_of_fup=True, end_of_fup_time=5, run=False) + eof.expand() + surv = _eof_run(data=d, run=False) + surv.expand() + + assert eof.DT.height > surv.DT.height + # Truncation leaves at most one outcome row per trial; without it a trial + # keeps several + events_per_trial = pl.col("outcome").fill_null(0).sum() + surv_max = ( + surv.DT.group_by(["ID", "trial"]).agg(events_per_trial.alias("n"))["n"].max() + ) + eof_max = ( + eof.DT.group_by(["ID", "trial"]).agg(events_per_trial.alias("n"))["n"].max() + ) + assert surv_max == 1 + assert eof_max > 1 + + # And the estimate is still readable at k for trials whose status was 1 earlier + eof.fit() + eof.end_of_followup() + assert eof.eof_data["Proportion"].is_finite().all() + + +def test_estimates_table_reports_censored_share_and_partition(): + s = _eof_run(end_of_fup=True, end_of_fup_time=12, end_of_fup_window=3) + d = s.eof_data + assert "% Censored" in d.columns + # Eligible = Analysed + Censored + No measurement, per arm + total = ( + d["Trial-periods (Analysed)"] + + d["Trial-periods (Censored)"] + + d["Trial-periods (No measurement)"] + ) + assert total.to_list() == d["Trial-periods (Eligible)"].to_list() + assert d["% Censored"].to_list() == approx( + ( + 100 * d["Trial-periods (Censored)"] / d["Trial-periods (Eligible)"] + ).to_list() + ) + + +def test_counts_table_accounts_for_every_trial_period(): + s = _eof_run(end_of_fup=True, end_of_fup_time=12, end_of_fup_window=3) + counts = s.diagnostics["nonunique_eof"].sort("tx_init_bas") + levels = [ + "At k", + "In window", + "Excluded (outside window)", + "Excluded (no measurement)", + ] + # Trial-period categories partition Eligible + total = sum(counts[lv] for lv in levels) + assert total.to_list() == counts["Eligible"].to_list() + # At k + In window equals the analysed trial-periods in the estimate table + analysed = (counts["At k"] + counts["In window"]).to_list() + assert analysed == s.eof_data.sort("A")["Trial-periods (Analysed)"].to_list() + # The subject-level table exists too + assert s.diagnostics["unique_eof"].height == 2 + + +def test_zero_window_puts_every_contributor_at_k(): + s = _eof_run(end_of_fup=True, end_of_fup_time=12) + counts = s.diagnostics["nonunique_eof"] + assert (counts["In window"] == 0).all() + assert (counts["At k"] > 0).all() + + +def test_counts_tables_absent_unless_end_of_fup(): + s = _eof_run(run=False) + s.expand() + s.fit() + assert "unique_eof" not in s.diagnostics + assert "nonunique_eof" not in s.diagnostics + + +def test_missing_outcomes_permitted_only_in_end_of_fup_mode(): + rng = np.random.default_rng(0) + d = load_data("SEQdata") + mask = pl.Series(rng.random(d.height) < 0.2) + d = d.with_columns( + pl.when(mask).then(None).otherwise(pl.col("outcome")).alias("outcome") + ) + + # Permitted (null = "not measured at this time") in eof mode + s = _eof_run(data=d, end_of_fup=True, end_of_fup_time=12, end_of_fup_window=3) + assert s.eof_data.height == 2 + + # Rejected in ordinary survival mode + with pytest.raises(ValueError, match="missing"): + _eof_run(data=d, run=False) + + +def test_nonbinary_outcome_error_hints_at_end_of_fup_type(): + with pytest.raises(ValueError, match="continuous"): + _eof_run( + data=_continuous_data(), + outcome="cont", + end_of_fup=True, + end_of_fup_time=12, + run=False, + ) + + +def test_continuous_outcome_diagnostics(): + s = _eof_run( + data=_continuous_data(), + outcome="cont", + end_of_fup=True, + end_of_fup_time=12, + end_of_fup_type="continuous", + end_of_fup_window=3, + ) + # Outcome event tables are suppressed for a continuous outcome only... + assert "unique_outcomes" not in s.diagnostics + # ...replaced by the N/Mean/SD summary of the analysed measurements + summary = s.diagnostics["eof_summary"] + assert set(summary.columns) >= {"A", "N", "Mean", "SD"} + assert (summary["Mean"] > 1).all() + # Follow-up tables survive + assert "unique_followup" in s.diagnostics + + # Binary eof keeps the outcome tables + b = _eof_run(end_of_fup=True, end_of_fup_time=12) + assert "unique_outcomes" in b.diagnostics + + +def test_collect_and_retrieve_data_expose_eof_tables(): + s = _eof_run( + end_of_fup=True, end_of_fup_time=12, end_of_fup_window=3, bootstrap_nboot=3 + ) + out = s.collect() + assert out.outcome_models is None + assert out.retrieve_data("eof_data").height == 2 + assert out.retrieve_data("eof_comparison").height == 2 + assert out.retrieve_data("nonunique_eof").height == 2 + with pytest.raises(ValueError, match="not created"): + out.retrieve_data("eof_summary") # binary outcome: no summary table + + +def test_survival_and_hazard_blocked_in_eof_mode(): + s = _eof_run(end_of_fup=True, end_of_fup_time=12) + with pytest.raises(ValueError, match="end_of_followup"): + s.survival() + with pytest.raises(ValueError, match="end_of_followup"): + s.hazard() + + +def test_end_of_followup_requires_eof_mode_and_fit(): + plain = _eof_run(run=False) + plain.expand() + with pytest.raises(ValueError, match="end_of_fup=False"): + plain.end_of_followup() + + s = _eof_run(end_of_fup=True, end_of_fup_time=12, run=False) + s.expand() + with pytest.raises(ValueError, match="fit"): + s.end_of_followup() From b74bd437b4e299d673fd4645365e4f7b06d7a905 Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Fri, 28 Aug 2026 08:51:57 +0100 Subject: [PATCH 3/9] Allow for statsmodels 0.15.0 changes --- pySEQTarget/analysis/_survival_pred.py | 9 +++++---- pySEQTarget/helpers/_fix_categories.py | 26 ++++++++++++++++++++------ pyproject.toml | 3 +++ 3 files changed, 28 insertions(+), 10 deletions(-) diff --git a/pySEQTarget/analysis/_survival_pred.py b/pySEQTarget/analysis/_survival_pred.py index 725f412..81cefa8 100644 --- a/pySEQTarget/analysis/_survival_pred.py +++ b/pySEQTarget/analysis/_survival_pred.py @@ -2,7 +2,8 @@ import polars as pl from patsy import PatsyError, dmatrix -from ..helpers._fix_categories import _fix_categories_for_predict +from ..helpers._fix_categories import (_fix_categories_for_predict, + _model_design_info) from ..helpers._predict_model import _safe_predict from ._outcome_fit import _cast_categories @@ -58,7 +59,7 @@ def _cached_predict(model, X_cached, ref_column_names, data): _safe_predict on mismatch (e.g. a bootstrap resample that dropped a categorical level). """ - dinfo = model.model.data.design_info + dinfo = _model_design_info(model) if list(dinfo.column_names) == ref_column_names: probs = np.asarray(model.predict(X_cached, transform=False)) if not np.any(np.isnan(probs)): @@ -78,14 +79,14 @@ def _get_outcome_predictions(self, TxDT, idx=None): main = self.outcome_model[0] main_dict = main[idx] if idx is not None else main main_outcome = self._offloader.load_model(main_dict["outcome"]) - outcome_dinfo = main_outcome.model.data.design_info + outcome_dinfo = _model_design_info(main_outcome) X_outcome = _build_design_matrix(outcome_dinfo, data) outcome_cols = list(outcome_dinfo.column_names) X_compevent = compevent_cols = None if self.compevent_colname is not None: main_compevent = self._offloader.load_model(main_dict["compevent"]) - compevent_dinfo = main_compevent.model.data.design_info + compevent_dinfo = _model_design_info(main_compevent) X_compevent = _build_design_matrix(compevent_dinfo, data) compevent_cols = list(compevent_dinfo.column_names) diff --git a/pySEQTarget/helpers/_fix_categories.py b/pySEQTarget/helpers/_fix_categories.py index f1a343f..ea4fe2d 100644 --- a/pySEQTarget/helpers/_fix_categories.py +++ b/pySEQTarget/helpers/_fix_categories.py @@ -1,16 +1,30 @@ import pandas as pd +def _model_design_info(model): + """The patsy DesignInfo a fitted model was built with, or None. + + statsmodels 0.15 renamed ``model.data.design_info`` to the engine-neutral + ``model.data.model_spec`` (still a patsy DesignInfo under the default patsy + engine); older statsmodels and the glum/jax wrappers use ``design_info``. + Check both so every backend and version resolves through one place. + """ + inner = getattr(model, "model", None) + data = getattr(inner, "data", None) + if data is None: + return None + design_info = getattr(data, "design_info", None) + if design_info is None: + design_info = getattr(data, "model_spec", None) + return design_info + + def _fix_categories_for_predict(model, newdata): """ Fix categorical column ordering in newdata to match what the model expects. """ - if ( - hasattr(model, "model") - and hasattr(model.model, "data") - and hasattr(model.model.data, "design_info") - ): - design_info = model.model.data.design_info + design_info = _model_design_info(model) + if design_info is not None: for factor, factor_info in design_info.factor_infos.items(): if factor_info.type == "categorical": col_name = factor.name() diff --git a/pyproject.toml b/pyproject.toml index 4a3be99..f8a04a0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -98,4 +98,7 @@ pythonpath = ["."] testpaths = ["tests"] filterwarnings = [ "ignore:FigureCanvasAgg is non-interactive:UserWarning", + # joblib internals still assign array.shape; deprecated in NumPy 2.5, + # awaiting an upstream joblib fix + "ignore:Setting the shape on a NumPy array has been deprecated:DeprecationWarning", ] From c92b32dbd6673f266dc1a91bb41410cd82626af2 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" Date: Fri, 28 Aug 2026 07:52:13 +0000 Subject: [PATCH 4/9] Auto-format code --- pySEQTarget/analysis/_endoffup.py | 82 +++++++++++++++++-------------- tests/test_end_of_fup.py | 10 ++-- 2 files changed, 48 insertions(+), 44 deletions(-) diff --git a/pySEQTarget/analysis/_endoffup.py b/pySEQTarget/analysis/_endoffup.py index 0a3d431..aa6ff47 100644 --- a/pySEQTarget/analysis/_endoffup.py +++ b/pySEQTarget/analysis/_endoffup.py @@ -93,9 +93,7 @@ def _eof_period_flags(self, DT): valid = pl.col(self.outcome_col).is_not_null() in_window = valid & (pl.col("followup") >= k - w) & (pl.col("followup") <= k + w) - return DT.group_by( - [self.id_col, "trial"] + _eof_group_cols(self, DT) - ).agg( + return DT.group_by([self.id_col, "trial"] + _eof_group_cols(self, DT)).agg( [ (valid & (pl.col("followup") == k)).any().alias("at_k"), in_window.any().alias("in_window"), @@ -122,25 +120,28 @@ def _eof_estimate(self): by = _eof_group_cols(self, DT) measured = _eof_measure(self, DT).with_columns( - pl.col("weight").clip( - lower_bound=self.weight_min, upper_bound=self.weight_max - ) + pl.col("weight").clip(lower_bound=self.weight_min, upper_bound=self.weight_max) ) est = measured.group_by(by).agg( [ - ((pl.col("weight") * pl.col("eof_value")).sum() / pl.col("weight").sum()) - .alias("estimate"), + ( + (pl.col("weight") * pl.col("eof_value")).sum() / pl.col("weight").sum() + ).alias("estimate"), pl.len().alias("n"), pl.col(self.id_col).n_unique().alias("n_subjects"), ] ) - totals = _eof_period_flags(self, DT).group_by(by).agg( - [ - pl.len().alias("n_eligible"), - (pl.col("measured") & ~pl.col("in_window")).sum().alias("n_censored"), - (~pl.col("measured")).sum().alias("n_nomeasure"), - ] + totals = ( + _eof_period_flags(self, DT) + .group_by(by) + .agg( + [ + pl.len().alias("n_eligible"), + (pl.col("measured") & ~pl.col("in_window")).sum().alias("n_censored"), + (~pl.col("measured")).sum().alias("n_nomeasure"), + ] + ) ) return est.join(totals, on=by, how="inner").sort(by) @@ -156,7 +157,12 @@ def _eof_counts(self, DT, unique): fall into different categories for different trials. """ by = _eof_group_cols(self, DT) - levels = ["At k", "In window", "Excluded (outside window)", "Excluded (no measurement)"] + levels = [ + "At k", + "In window", + "Excluded (outside window)", + "Excluded (no measurement)", + ] flags = _eof_period_flags(self, DT).with_columns( pl.when(pl.col("at_k")) @@ -169,9 +175,7 @@ def _eof_counts(self, DT, unique): .alias("_category") ) - counter = ( - pl.col(self.id_col).n_unique() if unique else pl.len() - ) + counter = pl.col(self.id_col).n_unique() if unique else pl.len() counted = flags.group_by(by + ["_category"]).agg(counter.alias("N")) wide = counted.pivot(on="_category", index=by, values="N").fill_null(0) for lv in levels: @@ -220,7 +224,9 @@ def _create_endoffup(self): full = self.outcome_model[0]["eof"] boots = [ - m["eof"] for m in self.outcome_model[1:] if m["eof"] is not None and m["eof"].height > 0 + m["eof"] + for m in self.outcome_model[1:] + if m["eof"] is not None and m["eof"].height > 0 ] has_ci = len(boots) > 1 sub_in_full = sub is not None and sub in full.columns @@ -269,9 +275,7 @@ def boot_estimates(key): se = float(np.nanstd(draws, ddof=1)) point = float( full.filter( - pl.all_horizontal( - [pl.col(c) == v for c, v in zip(key_cols, key)] - ) + pl.all_horizontal([pl.col(c) == v for c, v in zip(key_cols, key)]) )["estimate"][0] ) if use_se: @@ -294,25 +298,29 @@ def boot_estimates(key): } ) rename_back = {key_cols[-1]: "A"} - data = data.join(ci_frame.rename(rename_back), on=(["A"] if not sub_in_full else [sub, "A"])) - - lead = ["Type", "Time"] + ([sub] if sub_in_full else []) + [ - "A", - label, - "Trial-periods (Eligible)", - "Trial-periods (Analysed)", - "Trial-periods (Censored)", - "Trial-periods (No measurement)", - "% Censored", - "Subjects", - ] + data = data.join( + ci_frame.rename(rename_back), on=(["A"] if not sub_in_full else [sub, "A"]) + ) + + lead = ( + ["Type", "Time"] + + ([sub] if sub_in_full else []) + + [ + "A", + label, + "Trial-periods (Eligible)", + "Trial-periods (Analysed)", + "Trial-periods (Censored)", + "Trial-periods (No measurement)", + "% Censored", + "Subjects", + ] + ) data = data.select(lead + [c for c in data.columns if c not in lead]) # Pairwise between-arm contrasts, both directions ======================== rows = [] - subgroup_vals = ( - sorted(set(full[sub].to_list())) if sub_in_full else [None] - ) + subgroup_vals = sorted(set(full[sub].to_list())) if sub_in_full else [None] for g in subgroup_vals: sub_frame = full.filter(pl.col(sub) == g) if g is not None else full arms = sub_frame[tx_bas].to_list() diff --git a/tests/test_end_of_fup.py b/tests/test_end_of_fup.py index 7946032..10fd2d2 100644 --- a/tests/test_end_of_fup.py +++ b/tests/test_end_of_fup.py @@ -9,8 +9,8 @@ from pytest import approx from pySEQTarget import SEQopts, SEQuential -from pySEQTarget.data import load_data from pySEQTarget.analysis._endoffup import _eof_measure +from pySEQTarget.data import load_data def _eof_run(data=None, outcome="outcome", method="ITT", run=True, **opts): @@ -222,9 +222,7 @@ def test_incompatible_options_raise(): with pytest.raises(ValueError, match="km_curves or hazard"): _eof_run(end_of_fup=True, end_of_fup_time=12, hazard_estimate=True, run=False) with pytest.raises(ValueError, match="dose-response"): - _eof_run( - end_of_fup=True, end_of_fup_time=12, method="dose-response", run=False - ) + _eof_run(end_of_fup=True, end_of_fup_time=12, method="dose-response", run=False) with pytest.raises(ValueError, match="compevent"): _eof_run( end_of_fup=True, @@ -317,9 +315,7 @@ def test_estimates_table_reports_censored_share_and_partition(): ) assert total.to_list() == d["Trial-periods (Eligible)"].to_list() assert d["% Censored"].to_list() == approx( - ( - 100 * d["Trial-periods (Censored)"] / d["Trial-periods (Eligible)"] - ).to_list() + (100 * d["Trial-periods (Censored)"] / d["Trial-periods (Eligible)"]).to_list() ) From f9725c2bd0fb5b629fc2c6d4ee47188037a0f90e Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Fri, 28 Aug 2026 10:28:18 +0100 Subject: [PATCH 5/9] Bump version --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index f8a04a0..44e87be 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "pySEQTarget" -version = "0.14.0" +version = "0.14.1" description = "Sequential Target Trial Emulation" readme = "README.md" license = {text = "MIT"} From 502d47c55975cba71cf0e656ed1426579b1237e1 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" Date: Fri, 28 Aug 2026 09:28:34 +0000 Subject: [PATCH 6/9] Auto-format code --- CITATION.cff | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CITATION.cff b/CITATION.cff index 4f3ca12..59d3f48 100644 --- a/CITATION.cff +++ b/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/ From 700fc2ffae8723fb7a6e38d0da2b880d6ccc4f88 Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Fri, 28 Aug 2026 12:35:45 +0100 Subject: [PATCH 7/9] Add end-of-follow-up outcomes vignette --- docs/index.rst | 1 + docs/vignettes/end_of_followup.md | 280 ++++++++++++++++++++++++++++++ 2 files changed, 281 insertions(+) create mode 100644 docs/vignettes/end_of_followup.md diff --git a/docs/index.rst b/docs/index.rst index 9c5e394..f7c2ca0 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -17,4 +17,5 @@ vignettes/getting_started vignettes/more_advanced_models + vignettes/end_of_followup vignettes/exploring_results diff --git a/docs/vignettes/end_of_followup.md b/docs/vignettes/end_of_followup.md new file mode 100644 index 0000000..b25b977 --- /dev/null +++ b/docs/vignettes/end_of_followup.md @@ -0,0 +1,280 @@ +# End-of-Follow-up Outcomes + +## What this is for + +The rest of pySEQTarget estimates survival outcomes: a binary event that may occur at any point during follow-up, summarised through risks, survival curves or a hazard ratio. + +An end-of-follow-up outcome is different. It is measured once, at a single follow-up time chosen by the user — a biomarker at 12 months, disease status at two years, a questionnaire score at the end of the trial. There is no time-to-event to model, and the quantity of interest is simply the average outcome in each treatment arm. + +The {py:class}`~pySEQTarget.SEQopts` option `end_of_fup=True` switches to that estimand. For each trial-period the outcome is read at the requested follow-up time and averaged within each baseline treatment arm, weighted by the period-trial-specific weight at the time the measurement was taken. For a binary outcome this is the weighted proportion in each arm; for a continuous outcome, the weighted mean. + +Because there is no outcome model, `end_of_fup` cannot be combined with `km_curves`, `hazard_estimate`, `compevent_colname`, or the dose-response method. + +## A minimal example + +`end_of_fup_time` is the follow-up time `k` at which the outcome is evaluated, counted in follow-up periods since trial enrollment (not calendar time). After {py:meth}`~pySEQTarget.SEQuential.fit`, the estimates are assembled with {py:meth}`~pySEQTarget.SEQuential.end_of_followup` and land in the `eof_data` and `eof_comparison` attributes. + +```python +from pySEQTarget import SEQopts, SEQuential +from pySEQTarget.data import load_data + +options = SEQopts(end_of_fup=True, + # evaluate the outcome 12 follow-up periods after enrollment + end_of_fup_time=12, + # "binary" reports a proportion, "continuous" a mean + end_of_fup_type="binary", + bootstrap_nboot=20, + # fixes the bootstrap resamples, so the intervals below are + # reproducible + seed=1636) + +model = SEQuential(load_data("SEQdata"), + id_col="ID", + time_col="time", + eligible_col="eligible", + treatment_col="tx_init", + outcome_col="outcome", + time_varying_cols=["N", "L", "P"], + fixed_cols=["sex"], + method="ITT", + parameters=options) +model.expand() +model.bootstrap() +model.fit() +model.end_of_followup() + +model.eof_data.select(["A", "Proportion", "SE", "95% LCI", "95% UCI", + "% Censored", "Subjects"]) +``` + +```text +┌─────┬────────────┬──────────┬──────────┬──────────┬────────────┬──────────┐ +│ A ┆ Proportion ┆ SE ┆ 95% LCI ┆ 95% UCI ┆ % Censored ┆ Subjects │ +╞═════╪════════════╪══════════╪══════════╪══════════╪════════════╪══════════╡ +│ 0 ┆ 0.019865 ┆ 0.003014 ┆ 0.013956 ┆ 0.025773 ┆ 19.687712 ┆ 266 │ +│ 1 ┆ 0.027256 ┆ 0.002342 ┆ 0.022666 ┆ 0.031847 ┆ 28.46412 ┆ 236 │ +└─────┴────────────┴──────────┴──────────┴──────────┴────────────┴──────────┘ +``` + +`eof_data` gives the weighted proportion (or mean) in each arm with its bootstrap confidence interval, and an account of how much of the arm it rests on. `Trial-periods (Eligible)` is every trial-period that reached the follow-up time, and the next three partition it: `(Analysed)` contribute to the estimate, `(Censored)` were measured at some point but not within the window, and `(No measurement)` were never measured at all — so the three sum back to the eligible total. `% Censored` gives the censored share of that total. `Subjects` counts the distinct people behind the analysed trial-periods; one subject contributes several trial-periods and can be analysed in some and censored in others. + +```python +model.eof_comparison.select(["A_x", "A_y", "Difference", "Difference 95% LCI", + "Difference 95% UCI", "Difference SE"]) +``` + +```text +┌─────┬─────┬────────────┬────────────────────┬────────────────────┬───────────────┐ +│ A_x ┆ A_y ┆ Difference ┆ Difference 95% LCI ┆ Difference 95% UCI ┆ Difference SE │ +╞═════╪═════╪════════════╪════════════════════╪════════════════════╪═══════════════╡ +│ 0 ┆ 1 ┆ 0.007392 ┆ -0.000471 ┆ 0.015255 ┆ 0.004012 │ +│ 1 ┆ 0 ┆ -0.007392 ┆ -0.015255 ┆ 0.000471 ┆ 0.004012 │ +└─────┴─────┴────────────┴────────────────────┴────────────────────┴───────────────┘ +``` + +`eof_comparison` gives the pairwise between-arm contrast: the difference in proportions here, or the difference in means for a continuous outcome, with its standard error and confidence interval. For a binary outcome the ratio of proportions is reported alongside it, with an interval computed on the log scale and a `log(Ratio) SE` for inverse-variance pooling. Contrasts are paired by bootstrap iteration, so their intervals account for the correlation between arms. Both directions of each arm pair are reported, so the row you want is the one whose `A_x` is your reference arm. + +## Missing measurements and the time window + +Outcomes measured at particular visits are rarely available for everyone at exactly time `k`. Encode "not measured at this time" as a missing value (null) in the outcome column — `end_of_fup` is the one mode that accepts missing outcomes, precisely because missingness is meaningful here. + +`end_of_fup_window` sets the half-width of a window used when a trial-period has no measurement at exactly `k`: + +```python +options = SEQopts(end_of_fup=True, + end_of_fup_time=12, + # accept a measurement anywhere in [9, 15] when there is none at 12 + end_of_fup_window=3, + bootstrap_nboot=20, + seed=1636) + +windowed = SEQuential(load_data("SEQdata"), + id_col="ID", time_col="time", eligible_col="eligible", + treatment_col="tx_init", outcome_col="outcome", + time_varying_cols=["N", "L", "P"], fixed_cols=["sex"], + method="ITT", parameters=options) +windowed.expand() +windowed.bootstrap() +windowed.fit() +windowed.end_of_followup() + +windowed.eof_data.select(["A", "Proportion", "95% LCI", "95% UCI", + "Trial-periods (Eligible)", "Trial-periods (Analysed)", + "Trial-periods (Censored)", "% Censored"]) +``` + +```text +┌─────┬────────────┬──────────┬──────────┬──────────────┬──────────────┬──────────────┬────────────┐ +│ A ┆ Proportion ┆ 95% LCI ┆ 95% UCI ┆ Trial-period ┆ Trial-period ┆ Trial-period ┆ % Censored │ +│ ┆ ┆ ┆ ┆ s (Eligible) ┆ s (Analysed) ┆ s (Censored) ┆ │ +╞═════╪════════════╪══════════╪══════════╪══════════════╪══════════════╪══════════════╪════════════╡ +│ 0 ┆ 0.075704 ┆ 0.058602 ┆ 0.092805 ┆ 2946 ┆ 2523 ┆ 423 ┆ 14.358452 │ +│ 1 ┆ 0.102475 ┆ 0.0864 ┆ 0.11855 ┆ 6257 ┆ 4889 ┆ 1368 ┆ 21.863513 │ +└─────┴────────────┴──────────┴──────────┴──────────────┴──────────────┴──────────────┴────────────┘ +``` + +```python +windowed.eof_comparison.select(["A_x", "A_y", "Difference", "Ratio", + "Ratio 95% LCI", "Ratio 95% UCI", "log(Ratio) SE"]) +``` + +```text +┌─────┬─────┬────────────┬──────────┬───────────────┬───────────────┬───────────────┐ +│ A_x ┆ A_y ┆ Difference ┆ Ratio ┆ Ratio 95% LCI ┆ Ratio 95% UCI ┆ log(Ratio) SE │ +╞═════╪═════╪════════════╪══════════╪═══════════════╪═══════════════╪═══════════════╡ +│ 0 ┆ 1 ┆ 0.026771 ┆ 1.353635 ┆ 1.024306 ┆ 1.788848 ┆ 0.142237 │ +│ 1 ┆ 0 ┆ -0.026771 ┆ 0.738752 ┆ 0.559019 ┆ 0.976271 ┆ 0.142237 │ +└─────┴─────┴────────────┴──────────┴───────────────┴───────────────┴───────────────┘ +``` + +The selection rule, applied to each trial-period independently, is: + +1. If there is a measurement at exactly `k`, use it. +2. Otherwise, use the measurement *nearest* to `k` within `[k - window, k + window]`. Where two measurements are equally far either side of `k`, the *later* one is taken, so that at least `k` of follow-up has elapsed. +3. If there is no measurement anywhere in the window, the trial-period is *censored* — it contributes nothing to the average. + +The weight used is always the weight at the time the chosen measurement was taken, not the weight at `k`. + +A window is not free. Widening it recovers trial-periods that would otherwise be dropped, but the measurements it recovers are taken further from the time you actually care about, and the trial-periods it recovers are not a random subset — a trial-period with no measurement at `k` is often one whose follow-up ended early. Treat the window as a trade-off between precision and how literally the estimate answers "the outcome at time `k`", and check how much of the estimate rests on it using the accounting table below. + +## Checking what contributed + +The diagnostics report where every trial-period went. `nonunique_eof` counts trial-periods and `unique_eof` counts distinct subjects: + +```python +windowed.diagnostics["nonunique_eof"] +``` + +```text +┌─────────────┬──────────┬──────┬───────────┬───────────────────┬───────────────────────────┐ +│ tx_init_bas ┆ Eligible ┆ At k ┆ In window ┆ Excluded (outside ┆ Excluded (no measurement) │ +│ ┆ ┆ ┆ ┆ window) ┆ │ +╞═════════════╪══════════╪══════╪═══════════╪═══════════════════╪═══════════════════════════╡ +│ 0 ┆ 2946 ┆ 2366 ┆ 157 ┆ 423 ┆ 0 │ +│ 1 ┆ 6257 ┆ 4476 ┆ 413 ┆ 1368 ┆ 0 │ +└─────────────┴──────────┴──────┴───────────┴───────────────────┴───────────────────────────┘ +``` + +The four categories are mutually exclusive, so the trial-period counts partition `Eligible`: + +- *At k* — contributed, using a measurement at exactly `k`. +- *In window* — contributed, having fallen back to the window. +- *Excluded (outside window)* — measured at some point, but not within the window. +- *Excluded (no measurement)* — never measured at any follow-up time. Under `method="censoring"` this also picks up trial-periods artificially censored before any measurement was taken. + +`At k` plus `In window` is exactly the number of trial-periods behind the estimate, so the two tables always reconcile. The subject counts in `unique_eof` need *not* sum to `Eligible`, because one subject can fall into different categories for different trials. + +If `In window` is large relative to `At k`, or `Excluded` dominates, the estimate is resting on much less — or much more indirect — data than the arm totals alone suggest. + +These tables are also available from {py:meth}`~pySEQTarget.SEQoutput.retrieve_data` after {py:meth}`~pySEQTarget.SEQuential.collect`, as `"unique_eof"`, `"nonunique_eof"`, `"eof_data"` and `"eof_comparison"`. + +## Continuous outcomes + +Set `end_of_fup_type="continuous"` for an outcome that is not 0/1. The estimate becomes a weighted mean, reported in a `Mean` column rather than `Proportion`, and its confidence interval is not clamped to `[0, 1]`. The between-arm contrast is the difference in means; no ratio is reported, since a continuous outcome need not be bounded away from zero. + +```python +import numpy as np +import polars as pl + +rng = np.random.default_rng(42) +data = load_data("SEQdata") +data = data.with_columns( + (10 + 2 * pl.col("tx_init") + pl.col("N") + + pl.Series(rng.standard_normal(data.height))).alias("biomarker") +) + +continuous = SEQuential(data, + id_col="ID", time_col="time", eligible_col="eligible", + treatment_col="tx_init", outcome_col="biomarker", + time_varying_cols=["N", "L", "P"], fixed_cols=["sex"], + method="ITT", + parameters=SEQopts(end_of_fup=True, + end_of_fup_time=12, + end_of_fup_type="continuous", + end_of_fup_window=3, + bootstrap_nboot=20, + seed=1636)) +continuous.expand() +continuous.bootstrap() +continuous.fit() +continuous.end_of_followup() + +continuous.eof_data.select(["A", "Mean", "SE", "95% LCI", "95% UCI"]) +``` + +```text +┌─────┬───────────┬──────────┬───────────┬───────────┐ +│ A ┆ Mean ┆ SE ┆ 95% LCI ┆ 95% UCI │ +╞═════╪═══════════╪══════════╪═══════════╪═══════════╡ +│ 0 ┆ 21.703792 ┆ 0.105694 ┆ 21.496635 ┆ 21.910949 │ +│ 1 ┆ 21.767593 ┆ 0.091838 ┆ 21.587594 ┆ 21.947591 │ +└─────┴───────────┴──────────┴───────────┴───────────┘ +``` + +```python +continuous.eof_comparison +``` + +```text +┌──────┬─────┬─────┬────────────┬────────────────────┬────────────────────┬───────────────┐ +│ Time ┆ A_x ┆ A_y ┆ Difference ┆ Difference 95% LCI ┆ Difference 95% UCI ┆ Difference SE │ +╞══════╪═════╪═════╪════════════╪════════════════════╪════════════════════╪═══════════════╡ +│ 12.0 ┆ 0 ┆ 1 ┆ 0.063801 ┆ -0.16125 ┆ 0.288851 ┆ 0.114824 │ +│ 12.0 ┆ 1 ┆ 0 ┆ -0.063801 ┆ -0.288851 ┆ 0.16125 ┆ 0.114824 │ +└──────┴─────┴─────┴────────────┴────────────────────┴────────────────────┴───────────────┘ +``` + +Here `Difference` is the difference in means, and there is no `Ratio` column. + +Note that the usual outcome diagnostic tables count `outcome == 1` rows, which has no meaning for a continuous outcome, so they are omitted. In their place the diagnostics report the N, mean and SD of the raw analysed measurements per arm in `eof_summary`; the follow-up and end-of-follow-up tables remain available. + +```python +continuous.diagnostics["eof_summary"] +``` + +```text +┌─────┬──────┬───────────┬──────────┐ +│ A ┆ N ┆ Mean ┆ SD │ +╞═════╪══════╪═══════════╪══════════╡ +│ 0 ┆ 2523 ┆ 21.703792 ┆ 5.067746 │ +│ 1 ┆ 4889 ┆ 21.767593 ┆ 5.180225 │ +└─────┴──────┴───────────┴──────────┘ +``` + +## Per-protocol effects + +`end_of_fup` composes with weighting in the usual way, so a per-protocol end-of-follow-up effect is the censoring method with `weighted=True`: + +```python +perprotocol = SEQuential(load_data("SEQdata"), + id_col="ID", time_col="time", eligible_col="eligible", + treatment_col="tx_init", outcome_col="outcome", + time_varying_cols=["N", "L", "P"], fixed_cols=["sex"], + method="censoring", + parameters=SEQopts(weighted=True, + numerator="sex", + denominator="N + L + P + sex", + end_of_fup=True, + end_of_fup_time=12, + end_of_fup_window=3, + bootstrap_nboot=20, + seed=1636)) +perprotocol.expand() +perprotocol.bootstrap() +perprotocol.fit() +perprotocol.end_of_followup() + +perprotocol.eof_data.select(["A", "Proportion", "95% LCI", "95% UCI", + "% Censored", "Subjects"]) +``` + +```text +┌─────┬────────────┬──────────┬──────────┬────────────┬──────────┐ +│ A ┆ Proportion ┆ 95% LCI ┆ 95% UCI ┆ % Censored ┆ Subjects │ +╞═════╪════════════╪══════════╪══════════╪════════════╪══════════╡ +│ 0 ┆ 0.045069 ┆ 0.011689 ┆ 0.07845 ┆ 82.959946 ┆ 86 │ +│ 1 ┆ 0.088733 ┆ 0.065572 ┆ 0.111894 ┆ 49.592456 ┆ 211 │ +└─────┴────────────┴──────────┴──────────┴────────────┴──────────┘ +``` + +Subjects who deviate from their assigned strategy are artificially censored at the point of deviation, and their outcome is missing from then on. A trial-period that deviates before `k` therefore has no measurement to contribute and is excluded rather than carried forward — it appears under `Excluded (outside window)` or `Excluded (no measurement)` in the accounting table, depending on what it had measured earlier. That is why `% Censored` is so much larger here than under ITT. From f8d9e8c9e87e6062737f9b5021ceb9cba6dd6d48 Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Fri, 28 Aug 2026 12:45:59 +0100 Subject: [PATCH 8/9] Annotate the two expected SingularMatrixWarnings from statsmodels 0.15.0 --- tests/test_categorical_covariates.py | 9 ++++++++- tests/test_followup_options.py | 11 +++++++++++ 2 files changed, 19 insertions(+), 1 deletion(-) diff --git a/tests/test_categorical_covariates.py b/tests/test_categorical_covariates.py index 5aa4019..740f15b 100644 --- a/tests/test_categorical_covariates.py +++ b/tests/test_categorical_covariates.py @@ -1,5 +1,7 @@ import numpy as np import polars as pl +import pytest +from statsmodels.tools.sm_exceptions import SingularMatrixWarning from pySEQTarget import SEQopts, SEQuential from pySEQTarget.data import load_data @@ -63,7 +65,12 @@ def test_rare_level_not_dropped_by_bootstrap_resample(): s = _model(data, bootstrap_nboot=8, bootstrap_sample=0.5, seed=7) s.expand() s.bootstrap() - s.fit() + # The frozen level set makes resamples lacking level "c" carry an all-zero + # design column, so statsmodels (>= 0.15) reports the deliberate rank + # deficiency; pinv resolves it with a zero coefficient. The warning is the + # expected signature of the mechanism under test. + with pytest.warns(SingularMatrixWarning): + s.fit() s.survival() # previously raised ValueError on NaN predictions rd = s.risk_estimates["risk_difference"] diff --git a/tests/test_followup_options.py b/tests/test_followup_options.py index c2005d5..4cd6c65 100644 --- a/tests/test_followup_options.py +++ b/tests/test_followup_options.py @@ -1,3 +1,5 @@ +import pytest + from pySEQTarget import SEQopts, SEQuential from pySEQTarget.data import load_data @@ -38,6 +40,15 @@ def test_followup_class(): assert [round(x, 3) for x in matrix] == [round(x, 3) for x in expected] +# patsy's natural-cubic cr() basis spans the constant function, so cr() plus +# the model Intercept is always collinear by exactly one dimension. Fitted +# values (all the spline path is used for) are unique regardless; only the +# individual spline coefficients are non-identified, resolved deterministically +# by statsmodels' pinv. statsmodels >= 0.15 reports the redundancy as a +# SingularMatrixWarning. +@pytest.mark.filterwarnings( + "ignore::statsmodels.tools.sm_exceptions.SingularMatrixWarning" +) def test_followup_spline(): data = load_data("SEQdata") From c9615925c1fb314d949de07bf481eb1e5f8390a6 Mon Sep 17 00:00:00 2001 From: Tom Palmer Date: Wed, 2 Sep 2026 13:02:20 +0100 Subject: [PATCH 9/9] Shorten comments --- pySEQTarget/analysis/_endoffup.py | 6 ++---- pySEQTarget/error/_param_checker.py | 5 ++--- 2 files changed, 4 insertions(+), 7 deletions(-) diff --git a/pySEQTarget/analysis/_endoffup.py b/pySEQTarget/analysis/_endoffup.py index aa6ff47..2fdf095 100644 --- a/pySEQTarget/analysis/_endoffup.py +++ b/pySEQTarget/analysis/_endoffup.py @@ -62,10 +62,8 @@ def _eof_measure(self, DT): & (pl.col("followup") <= k + w) ).select(cols) - # Nearest to k wins (exact k has distance 0); equidistant ties break toward - # the later measurement so at least k of follow-up has elapsed. followup is - # UInt32 in the expanded frame, so cast before differencing — unsigned - # subtraction wraps for followups below k and .abs() cannot recover it. + # Nearest to k wins; equidistant ties break toward the later measurement. + # Cast before differencing - followup is UInt32 and wraps below k. measured = ( candidates.with_columns( (pl.col("followup").cast(pl.Int64) - k).abs().alias("_dist") diff --git a/pySEQTarget/error/_param_checker.py b/pySEQTarget/error/_param_checker.py index 0afe9c6..4ef1637 100644 --- a/pySEQTarget/error/_param_checker.py +++ b/pySEQTarget/error/_param_checker.py @@ -41,9 +41,8 @@ def _param_checker(self): if self.km_curves and self.hazard_estimate: raise ValueError("km_curves and hazard cannot both be set to True.") - # End-of-follow-up outcomes replace the survival outcome model with a direct - # weighted average at a single follow-up time, so the survival-based outputs - # have no meaning there and the requested time must lie inside the expansion. + # end_of_fup replaces the outcome model with an average at one time - the + # survival outputs have no meaning and k +/- window must fit the expansion. if self.end_of_fup: if self.km_curves or self.hazard_estimate: raise ValueError(