Skills › Scientific Agent Skills › scikit-survival › competing_risk_cif.py
competing_risk_cif.py
317 lines · skills/scikit-survival/scripts/competing_risk_cif.py at 5f2bc03
1#!/usr/bin/env python3
2"""Estimate nonparametric competing-risk cumulative incidence from local data."""
4from __future__ import annotations
6import argparse
7from typing import Any
9from _common import (
10 DEFAULT_SEED,
11 CliError,
12 atomic_save_npz,
13 bounded_int,
14 emit_json,
15 parse_floats,
16 probability,
17 read_csv,
18)
21VARIANCE_CHOICES = ("Aalen", "Dinse", "Dinse_Approx")
22MAX_CIF_ROWS = 5_000
23MAX_CAUSES = 32
26def synthetic_competing_risks(
27 *, rows: int = 300, seed: int = DEFAULT_SEED
28) -> tuple[Any, Any]:
29 """Create deterministic, non-clinical outcomes with three competing causes."""
31 try:
32 import numpy as np
33 import pandas as pd
34 except ImportError as exc:
35 raise CliError("NumPy and pandas are required") from exc
36 rng = np.random.default_rng(seed)
37 latent = rng.normal(size=rows)
38 cause_times = np.column_stack(
39 [
40 rng.exponential(scale=7.0 * np.exp(-0.25 * latent)),
41 rng.exponential(scale=10.0 * np.exp(0.15 * latent)),
42 rng.exponential(scale=14.0, size=rows),
43 ]
44 )
45 censor_time = rng.exponential(scale=16.0, size=rows)
46 cause = np.argmin(cause_times, axis=1) + 1
47 event_time = cause_times[np.arange(rows), cause - 1]
48 event = np.where(event_time <= censor_time, cause, 0)
49 time = np.minimum(event_time, censor_time) + 0.05
50 for cause_code in (1, 2, 3):
51 event[cause_code - 1] = cause_code
52 time[cause_code - 1] = float(cause_code)
53 frame = pd.DataFrame({"status": event, "time": time})
54 return frame, {"kind": "synthetic", "network_used": False, "seed": seed}
57def normalize_competing_event(values: Any):
58 """Return non-negative contiguous integer cause codes."""
60 try:
61 import numpy as np
62 import pandas as pd
63 except ImportError as exc:
64 raise CliError("NumPy and pandas are required") from exc
65 series = pd.Series(values)
66 if series.isna().any():
67 raise CliError("competing-risk event codes contain missing values")
68 try:
69 numeric = pd.to_numeric(series, errors="raise").to_numpy(dtype=float)
70 except (TypeError, ValueError) as exc:
71 raise CliError("event codes must be integers") from exc
72 if not np.isfinite(numeric).all() or (numeric < 0).any():
73 raise CliError("event codes must be finite non-negative integers")
74 if not np.equal(numeric, np.floor(numeric)).all():
75 raise CliError("event codes must be integers")
76 if (numeric > MAX_CAUSES).any():
77 raise CliError(f"cause codes must not exceed {MAX_CAUSES}")
78 event = numeric.astype(int)
79 causes = sorted(set(event.tolist()) - {0})
80 if len(causes) < 2:
81 raise CliError("at least two competing causes are required")
82 if len(causes) > MAX_CAUSES:
83 raise CliError(f"at most {MAX_CAUSES} competing causes are allowed")
84 expected = list(range(1, max(causes) + 1))
85 if causes != expected:
86 raise CliError(
87 "positive cause codes must be contiguous and every cause must occur"
88 )
89 return event
92def normalize_times(values: Any):
93 """Return a finite, strictly positive time vector."""
95 from _common import normalize_positive_times
97 return normalize_positive_times(values, label="competing-risk time")
100def horizon_values(
101 unique_times: Any,
102 cumulative_incidence: Any,
103 horizons: list[float],
104) -> list[dict[str, Any]]:
105 """Evaluate right-continuous step estimates at requested horizons."""
107 import numpy as np
109 records: list[dict[str, Any]] = []
110 for horizon in horizons:
111 if horizon <= 0 or not np.isfinite(horizon):
112 raise CliError("horizons must be finite and strictly positive")
113 if horizon > unique_times[-1]:
114 raise CliError("report horizons must not exceed observed follow-up")
115 index = int(np.searchsorted(unique_times, horizon, side="right") - 1)
116 values = (
117 np.zeros(cumulative_incidence.shape[0], dtype=float)
118 if index < 0
119 else cumulative_incidence[:, index]
120 )
121 records.append(
122 {
123 "cause_specific": {
124 str(cause): float(values[cause])
125 for cause in range(1, cumulative_incidence.shape[0])
126 },
127 "horizon": float(horizon),
128 "total_risk": float(values[0]),
129 }
130 )
131 return records
134def estimate_cif(
135 event: Any,
136 time: Any,
137 *,
138 horizons: list[float] | None = None,
139 confidence: bool = False,
140 confidence_level: float = 0.95,
141 variance_type: str = "Aalen",
142 time_min: float | None = None,
143) -> tuple[dict[str, Any], dict[str, Any]]:
144 """Estimate CIF curves and return a bounded summary plus numeric arrays."""
146 try:
147 import numpy as np
148 from sksurv.nonparametric import (
149 cumulative_incidence_competing_risks,
150 kaplan_meier_estimator,
151 )
152 except ImportError as exc:
153 raise CliError("the pinned scikit-survival stack is required") from exc
154 checked_event = normalize_competing_event(event)
155 checked_time = normalize_times(time)
156 if len(checked_event) != len(checked_time):
157 raise CliError("event and time arrays must have equal length")
158 if len(checked_event) > MAX_CIF_ROWS:
159 raise CliError(f"CIF estimation is limited to {MAX_CIF_ROWS} rows")
160 if confidence and variance_type == "Dinse" and len(checked_event) > 2_000:
161 raise CliError("Dinse variance is limited to 2,000 rows")
162 if time_min is not None and (
163 not np.isfinite(time_min) or time_min < 0 or time_min >= checked_time.max()
164 ):
165 raise CliError("time_min must be finite, non-negative, and below max time")
166 if time_min is not None and time_min > checked_time.min():
167 raise CliError(
168 "scikit-survival 0.28.0 has a conditional CIF shape defect; "
169 "time_min that removes observed times is unsupported"
170 )
171 if not np.isfinite(confidence_level) or not 0 < confidence_level < 1:
172 raise CliError("confidence_level must be strictly between zero and one")
174 result = cumulative_incidence_competing_risks(
175 checked_event,
176 checked_time,
177 time_min=time_min,
178 conf_level=confidence_level,
179 conf_type="log-log" if confidence else None,
180 var_type=variance_type,
181 )
182 if confidence:
183 unique_times, cumulative_incidence, intervals = result
184 # 0.28.0 reverses the total-CIF endpoints and ignores conf_level for
185 # that row. Reconstruct it from the public all-event KM estimator.
186 km_times, _, km_intervals = kaplan_meier_estimator(
187 checked_event > 0,
188 checked_time,
189 conf_type="log-log",
190 conf_level=confidence_level,
191 )
192 if not np.array_equal(km_times, unique_times):
193 raise CliError("total-risk confidence grid differs from the CIF grid")
194 intervals[0] = 1.0 - km_intervals[::-1]
195 if not np.isfinite(intervals).all() or (intervals[:, 0] > intervals[:, 1]).any():
196 raise CliError("CIF confidence intervals are non-finite or reversed")
197 else:
198 unique_times, cumulative_incidence = result
199 intervals = None
200 if not np.isfinite(cumulative_incidence).all() or (
201 (cumulative_incidence < -1e-10) | (cumulative_incidence > 1 + 1e-10)
202 ).any():
203 raise CliError("CIF estimates must be finite probabilities")
204 if (np.diff(cumulative_incidence, axis=1) < -1e-10).any():
205 raise CliError("estimated cumulative incidence unexpectedly decreased")
206 if not np.allclose(
207 cumulative_incidence[0],
208 cumulative_incidence[1:].sum(axis=0),
209 rtol=1e-8,
210 atol=1e-10,
211 ):
212 raise CliError("cause-specific CIFs do not sum to total risk")
214 if horizons is None:
215 horizons = [
216 float(value) for value in np.quantile(checked_time, [0.25, 0.5, 0.75])
217 ]
218 if len(horizons) > 32:
219 raise CliError("at most 32 reporting horizons are allowed")
220 summary = {
221 "cause_counts": {
222 str(code): int((checked_event == code).sum())
223 for code in range(0, int(checked_event.max()) + 1)
224 },
225 "confidence": {
226 "enabled": confidence,
227 "level": confidence_level if confidence else None,
228 "type": "log-log" if confidence else None,
229 "variance_type": variance_type if confidence else None,
230 },
231 "event_coding": "0=censored; positive contiguous integers=causes",
232 "horizons": horizon_values(
233 unique_times, cumulative_incidence, sorted(horizons)
234 ),
235 "n_causes": int(checked_event.max()),
236 "rows": int(len(checked_event)),
237 "schema_version": "1.0",
238 "time_min": time_min,
239 "warning": (
240 "Cumulative incidence is cause-specific absolute event probability, "
241 "not a cause-specific hazard or proof of treatment benefit."
242 ),
243 }
244 arrays = {
245 "cumulative_incidence": cumulative_incidence,
246 "time": unique_times,
247 }
248 if intervals is not None:
249 arrays["confidence_interval"] = intervals
250 return summary, arrays
253def build_parser() -> argparse.ArgumentParser:
254 parser = argparse.ArgumentParser(
255 description=(
256 "Estimate nonparametric cumulative incidence for competing causes. "
257 "No input uses deterministic synthetic data; no network is used."
258 )
259 )
260 parser.add_argument("--input", help="Optional bounded local .csv")
261 parser.add_argument("--event-column", default="status")
262 parser.add_argument("--time-column", default="time")
263 parser.add_argument("--horizons", help="Optional comma-separated report horizons")
264 parser.add_argument("--time-min", type=float)
265 parser.add_argument("--confidence", action="store_true")
266 parser.add_argument("--confidence-level", type=probability, default=0.95)
267 parser.add_argument("--variance-type", choices=VARIANCE_CHOICES, default="Aalen")
268 parser.add_argument("--seed", type=bounded_int(0, 2**32 - 1), default=DEFAULT_SEED)
269 parser.add_argument("--synthetic-rows", type=bounded_int(80, 20_000), default=300)
270 parser.add_argument(
271 "--curve-output", help="Optional local .npz containing full CIF arrays"
272 )
273 parser.add_argument("--output", help="Optional local .json summary")
274 parser.add_argument("--force", action="store_true")
275 return parser
278def main(argv: list[str] | None = None) -> int:
279 parser = build_parser()
280 args = parser.parse_args(argv)
281 try:
282 if args.input:
283 frame, name = read_csv(args.input)
284 source = {"kind": "local_csv", "name": name}
285 else:
286 frame, source = synthetic_competing_risks(
287 rows=args.synthetic_rows, seed=args.seed
288 )
289 missing = [
290 name
291 for name in (args.event_column, args.time_column)
292 if name not in frame.columns
293 ]
294 if missing:
295 raise CliError(f"missing columns: {', '.join(missing)}")
296 horizons = parse_floats(args.horizons) or None
297 summary, arrays = estimate_cif(
298 frame[args.event_column],
299 frame[args.time_column],
300 horizons=horizons,
301 confidence=args.confidence,
302 confidence_level=args.confidence_level,
303 variance_type=args.variance_type,
304 time_min=args.time_min,
305 )
306 summary["source"] = source
307 if args.curve_output:
308 atomic_save_npz(arrays, args.curve_output, force=args.force)
309 summary["curve_output"] = args.curve_output
310 emit_json(summary, output=args.output, force=args.force)
311 except CliError as exc:
312 parser.error(str(exc))
313 return 0
316if __name__ == "__main__":
317 raise SystemExit(main())