statspai.survival

survival

Survival and duration analysis models.

CausalSurvivalForestResult

Bases: GRFFamilyForest

A fitted causal survival forest.

Attributes:

Name	Type	Description
ate	float	Doubly-robust average effect (RMST difference or survival- probability difference at horizon).
ate_rmst	float	Same as ate when target="RMST" (kept for compatibility); NaN otherwise.
se, ci, pvalue	(float, tuple, float)
cate	ndarray	Out-of-bag tau(X_i).
cate_variance	ndarray or None
horizon	float
target	str
n_obs, n_trees	int
detail	dict

Examples:

>>> import numpy as np, statspai as sp
>>> rng = np.random.default_rng(0)
>>> X = rng.normal(size=(600, 2)); W = rng.binomial(1, 0.5, 600)
>>> T = rng.exponential(1 / np.exp(0.5 * X[:, 0] - 0.5 * W))
>>> C = rng.exponential(3, 600)
>>> fit = sp.causal_survival_forest(time=np.minimum(T, C),
...                                 event=(T <= C).astype(int), treat=W,
...                                 covariates=X, horizon=1.0,
...                                 n_estimators=200)
>>> fit.horizon
1.0

predict

predict(newdata: Any = None, estimate_variance: bool = False) -> DataFrame

tau(x) at newdata (OOB for the training rows when None).

get_scores

get_scores() -> ndarray

Doubly-robust scores tau(X_i) + (A_i - tau B_i) / (e_i (1 - e_i)).

CumIncResult dataclass

Bases: ResultProtocolMixin

Cumulative incidence functions for competing risks.

Attributes:

Name	Type	Description
cif_table	DataFrame	Long table with columns group ("all" when ungrouped), cause, time, cif, se, ci_lower, ci_upper.
causes	list	The competing-cause labels (excluding the 0 censoring code).
gray_test	dict or None	cause -> {"statistic", "df", "p_value"} from Gray's K-sample test, or None when group was not supplied.
alpha	float	Significance level used for the confidence bands.

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 200
>>> t1 = rng.exponential(scale=1.0, size=n)
>>> t2 = rng.exponential(scale=1.5, size=n)
>>> cens = rng.exponential(scale=2.0, size=n)
>>> time = np.minimum(np.minimum(t1, t2), cens)
>>> status = np.where((t1 <= t2) & (t1 <= cens), 1,
...                   np.where((t2 < t1) & (t2 <= cens), 2, 0))
>>> df = pd.DataFrame({"time": time, "status": status})
>>> ci = sp.cuminc(df, duration="time", event="status")
>>> isinstance(ci, sp.CumIncResult)
True
>>> ci.causes
[1, 2]
>>> at1 = ci.cif_at(1.0, cause=1)  # CIF for cause 1 at time t=1
>>> bool((at1["cif"] >= 0).all())
True

cif_at

cif_at(time: float, cause: Optional[int] = None) -> DataFrame

Cumulative incidence (last value at or before time).

plot

plot(cause: Optional[int] = None, ax: Any = None, **kwargs: Any) -> Any

Step-plot the cumulative incidence function(s).

FineGrayResult dataclass

Bases: ResultProtocolMixin

Fine-Gray proportional subdistribution hazards model result.

Attributes:

Name	Type	Description
params	ndarray	Estimated coefficients (log subdistribution hazard ratios).
bse	ndarray	Standard errors (Fine-Gray sandwich by default; see vce).
covariates	list of str	Covariate names aligned with params.
cause	int	The cause of interest whose subdistribution was modelled.
n_obs, n_events	int	Sample size and number of cause-of-interest events.
vcov	ndarray	Covariance matrix behind bse.
vce	str	"robust" (Fine-Gray sandwich), "robust+n/(n-1)" (Stata stcrreg scaling) or "model" (inverse information).

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 300
>>> x = rng.normal(size=n)
>>> t1 = rng.exponential(scale=np.exp(-0.5 * x))
>>> t2 = rng.exponential(scale=1.5)
>>> cens = rng.exponential(scale=2.0)
>>> time = np.minimum(np.minimum(t1, t2), cens)
>>> status = np.where((t1 <= t2) & (t1 <= cens), 1,
...                   np.where((t2 < t1) & (t2 <= cens), 2, 0))
>>> df = pd.DataFrame({"time": time, "status": status, "x": x})
>>> res = sp.finegray(df, duration="time", event="status", x=["x"])
>>> isinstance(res, sp.FineGrayResult)
True
>>> res.covariates
['x']
>>> res.shr.shape   # one subdistribution hazard ratio per covariate
(1,)

shr property

shr: ndarray

Subdistribution hazard ratios, exp(coef).

FrailtyResult dataclass

Bases: ResultProtocolMixin

lr_theta0 property

lr_theta0: Tuple[float, float]

Likelihood-ratio test of θ = 0: (chibar2, p).

The null sits on the boundary, so the p-value is half a chi2(1) tail (Stata's chibar2(01)).

CoxResult

Bases: EconometricResults

Result from Cox proportional hazards estimation.

Extends EconometricResults with survival-specific methods: .plot(), .ph_test(), .baseline_hazard(), .concordance.

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 120
>>> x = rng.normal(size=n)
>>> time = rng.exponential(scale=np.exp(-0.5 * x))
>>> event = (rng.random(n) < 0.8).astype(int)
>>> df = pd.DataFrame({"time": time, "status": event, "x": x})
>>> res = sp.cox(data=df, duration="time", event="status", x=["x"])
>>> type(res).__name__
'CoxResult'
>>> list(res.params.index)
['x']
>>> bool(isinstance(float(res.concordance), float))
True
>>> list(res.baseline_hazard().columns)
['time', 'baseline_cumhaz', 'baseline_survival']
>>> isinstance(res.summary(), str)
True

concordance property

concordance: float

Harrell's C-statistic (concordance index).

baseline_hazard

baseline_hazard() -> DataFrame

Baseline cumulative hazard (Breslow estimator).

Returns:

Type	Description
DataFrame	Columns time, baseline_cumhaz, baseline_survival.

ph_test

ph_test(transform: Any = 'km', method: str = 'score') -> DataFrame

Test the proportional hazards assumption (Grambsch and Therneau).

For each covariate the model is embedded in one whose coefficient varies with time, beta_j(t) = beta_j + theta_j * g(t), and theta_j = 0 is tested by a score test at the fitted model; the GLOBAL row tests all the theta_j jointly. No larger model is fitted. This is the test of survival::cox.zph in R (version 3.0 and later), and the statistics agree with it under either tie rule and with strata.

Parameters:

Name	Type	Description	Default
transform	(km, rank, identity, log)	The function of time g. 'km' (default, as in cox.zph) is one minus the left-continuous Kaplan-Meier estimate of the pooled sample, which spreads the event times evenly and keeps a few late events from dominating; 'identity' is time itself (the default of Stata estat phtest); 'rank' is the rank of the time among all observations (among the failures under method='approx', as in Stata); 'log' needs positive times. A callable is applied to the array of event times.	'km'
method	(score, approx)	'score' is the score test described above. 'approx' is the approximation of Grambsch and Therneau (1994) that Stata estat phtest and versions of survival before 3.0 report: it replaces the information matrix at each event time by the average over event times, which turns the test into a regression of the scaled Schoenfeld residuals on g(t). With transform='identity' it reproduces estat phtest, and with the other three its km, rank and log options. Unweighted right-censored fits only.	'score'

Returns:

Type	Description
DataFrame	One row per covariate and a final GLOBAL row, with columns variable, rho (correlation of the scaled Schoenfeld residuals with g(t); missing for GLOBAL), chi2, df and p_value.

Notes

Weighted fits (weights=) and counting-process fits (start=) get the score test of cox.zph for such fits: the case weights enter the score and the information as they enter the partial likelihood, the information is the model-based one whatever variance the fit reports, and g ignores the weights ('km' is the unweighted Kaplan-Meier of the (start, stop] data, 'rank' ranks the stop times of all rows). Fits with design= are refused.

The two methods differ most when the covariate distribution of the risk set changes over time, which is when the average information matrix is a poor stand-in for each event time's own.

Before 1.39 this method returned n * rho^2 with rho the Spearman correlation between time and the unscaled residuals. That number is not the Grambsch-Therneau statistic and was several times too small in places.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 400
>>> x = rng.binomial(1, 0.5, n)
>>> t = rng.exponential(1.0, n) ** np.where(x == 1, 2.0, 1.0)
>>> df = pd.DataFrame({"t": t, "d": 1, "x": x})
>>> fit = sp.cox(data=df, duration="t", event="d", x=["x"])
>>> bool(fit.ph_test().loc[0, "p_value"] < 0.01)
True

References

[@grambsch1994proportional]

plot

plot(kind: str = 'survival', ax: Any = None, **kwargs: Any) -> Any

Plot survival-related curves.

Parameters:

Name	Type	Description	Default
kind	str	'survival' — baseline survival curve, 'hazard' — baseline cumulative hazard, 'hr' — hazard ratio forest plot.	'survival'
ax	Axes		None
**kwargs	Any	type= is accepted as an alias of kind (the spelling other result classes use); the rest go to matplotlib.	{}

KMResult

Bases: ResultProtocolMixin

Kaplan-Meier survival analysis result.

Attributes:

Name	Type	Description
survival_table	DataFrame	Life table with time, n_risk, n_event, n_censor, survival, std_err, ci_lower, ci_upper.
median_survival	float or dict	Median survival time (per group if groups present).

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 120
>>> df = pd.DataFrame({
...     "time": rng.exponential(10, n).round(2),
...     "status": (rng.random(n) < 0.7).astype(int),
... })
>>> km = sp.kaplan_meier(data=df, duration="time", event="status",
...                      conf_type="log-log")
>>> type(km).__name__
'KMResult'
>>> "survival" in km.survival_table.columns
True
>>> bool(0.0 <= km.survival_table["survival"].iloc[-1] <= 1.0)
True
>>> isinstance(km.summary(), str)
True

survival_table property

survival_table: DataFrame

Life table (combined if multiple groups).

median_survival property

median_survival: Union[float, Dict[str, float]]

Median survival time (scalar or dict by group).

The first time at which S(t) <= 0.5, as Stata stci reports it. When S(t) equals 0.5 exactly on a step, R survfit reports the midpoint between that time and the next event time instead (deaths at 1, 2, 3, 4: Stata 2, R 2.5).

summary

summary() -> str

Formatted summary of the Kaplan-Meier analysis.

plot

plot(ax: Any = None, **kwargs: Any) -> Any

Plot Kaplan-Meier survival curves with confidence bands.

Parameters:

Name	Type	Description	Default
ax	Axes		None

causal_survival_forest

causal_survival_forest(data: Optional[DataFrame] = None, time: Any = None, event: Any = None, treat: Any = None, covariates: Any = None, horizon: Optional[float] = None, *, target: str = 'RMST', failure_times: Any = None, clusters: Any = None, weights: Any = None, equalize_cluster_weights: bool = False, W_hat: Any = None, n_estimators: int = 2000, min_samples_leaf: int = 5, max_samples: float = 0.5, mtry: Optional[int] = None, honest: bool = True, honesty_fraction: float = 0.5, honesty_prune_leaves: bool = True, split_alpha: float = 0.05, imbalance_penalty: float = 0.0, stabilize_splits: bool = True, ci_group_size: int = 2, max_depth: Optional[int] = None, propensity_bounds: Optional[Sequence[float]] = None, random_state: Optional[int] = 42, n_jobs: int = 1, alpha: float = 0.05) -> CausalSurvivalForestResult

Causal survival forest: heterogeneous effects on censored survival.

Parameters:

Name	Type	Description	Default
data	DataFrame	Input data; when omitted, the other inputs are arrays.	None
time	str or array - like	Observed time min(T, C), non-negative.	None
event	str or array - like	1 = event observed, 0 = censored.	None
treat	str or array - like	Binary treatment.	None
covariates	str, list of str or 2-D array	Covariates X (confounders and effect modifiers).	None
horizon	float	h of the estimand. Defaults to the 80th percentile of observed event times (recorded in detail); grf requires it -- choose it from the study design.	None
target	('RMST', 'survival_probability')		"RMST"
failure_times	array - like	Grid for the survival nuisance forests (default: observed times).	None
clusters	Any	As in :func:statspai.causal_forest.	None
weights	Any	As in :func:statspai.causal_forest.	None
equalize_cluster_weights	Any	As in :func:statspai.causal_forest.	None
W_hat	array - like	Precomputed propensities; default out-of-bag regression forest.	None
n_estimators	int		2000
min_samples_leaf	int		5
max_samples	float		0.5
mtry	float		0.5
honest	float		0.5
honesty_fraction	float		0.5
honesty_prune_leaves	float		0.5
imbalance_penalty	float	Forest options (grf names).	0.0
stabilize_splits	float	Forest options (grf names).	0.0
ci_group_size	float	Forest options (grf names).	0.0
max_depth	float	Forest options (grf names).	0.0
split_alpha	float	grf alpha; each child must also hold split_alpha of the parent's size in failures.	0.05
propensity_bounds	(float, float)	Clip e in the average-effect score denominators (off by default, as in grf).	None
random_state	int		42
n_jobs	int		1
alpha	float	Significance level of the reported interval.	0.05

Returns:

Type	Description
CausalSurvivalForestResult	.ate / .se / .ci / .pvalue, .cate (OOB), .predict(newdata), .best_linear_projection(A).

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 600; x = rng.normal(size=(n, 2)); w = rng.binomial(1, 0.5, n)
>>> T = rng.exponential(1 / np.exp(0.5 * x[:, 0] - 0.5 * w))
>>> C = rng.exponential(3, n)
>>> df = pd.DataFrame({"t": np.minimum(T, C), "d": (T <= C).astype(int),
...                    "w": w, "x1": x[:, 0], "x2": x[:, 1]})
>>> csf = sp.causal_survival_forest(df, time="t", event="d", treat="w",
...                                 covariates=["x1", "x2"], horizon=1.0,
...                                 n_estimators=200)
>>> csf.target
'RMST'

References

[@cui2023estimating], [@athey2019generalized]

cuminc

cuminc(data: DataFrame, duration: str, event: str, group: Optional[str] = None, alpha: float = 0.05, variance: str = 'delta', conf_type: str = 'linear', rho: float = 0.0) -> CumIncResult

Cumulative incidence functions for competing risks (Aalen-Johansen).

Parameters:

Name	Type	Description	Default
data	DataFrame	Input data.	required
duration	str	Column name for the follow-up time.	required
event	str	Column name for the event indicator. 0 = censored; 1, 2, ... = competing causes.	required
group	str	Column name for a grouping variable. When supplied, CIFs are estimated per group and Gray's K-sample test is reported per cause.	None
alpha	float	Significance level for the confidence bands.	0.05
variance	('delta', 'gray')	Pointwise variance of the CIF. "delta" is the Marubini-Valsecchi delta-method estimator (the one Stata's stcompet reports); "gray" is Gray's asymptotic variance, the var component of R cmprsk::cuminc.	"delta"
conf_type	('linear', 'log-log')	"linear" is F +/- z se clipped to [0, 1]; "log-log" is F^exp(+/- z se / (F log F)), stcompet's bounds.	"linear"
rho	float	Power of the (1 - F(t-))^rho weight in Gray's test (cmprsk's rho).	0

Returns:

Type	Description
CumIncResult	With .cif_table, .gray_test, .summary(), .plot().

Notes

The cumulative incidence for a single cause is not 1 - KM applied to that cause; treating competing events as censoring over-states the risk. The Aalen-Johansen estimator weights each cause-specific increment by the overall (all-cause) survival probability, so the CIFs of all causes plus the overall survival sum to one at every time.

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 200
>>> arm = rng.integers(0, 2, n)
>>> t1 = rng.exponential(scale=1.0 + 0.5 * arm, size=n)   # cause 1
>>> t2 = rng.exponential(scale=1.5, size=n)               # cause 2
>>> cens = rng.exponential(scale=2.0, size=n)
>>> time = np.minimum(np.minimum(t1, t2), cens)
>>> status = np.where((t1 <= t2) & (t1 <= cens), 1,
...                   np.where((t2 < t1) & (t2 <= cens), 2, 0))
>>> df = pd.DataFrame({"time": time, "status": status, "arm": arm})
>>> ci = sp.cuminc(df, duration="time", event="status", group="arm")
>>> isinstance(ci, sp.CumIncResult)
True
>>> ci.causes
[1, 2]
>>> ci.summary()
>>> ci.plot(cause=1)

finegray

finegray(data: DataFrame, duration: str, event: str, x: Sequence[str], cause: int = 1, alpha: float = 0.05, max_iter: int = 50, tol: float = 1e-10, vce: str = 'robust', small_sample: bool = False, weights: Optional[str] = None) -> FineGrayResult

Fine & Gray (1999) proportional subdistribution hazards model.

Models the effect of covariates on the cumulative incidence of cause through its subdistribution hazard, so coefficients exponentiate to subdistribution hazard ratios that map monotonically to the CIF (unlike cause-specific Cox coefficients).

Parameters:

Name	Type	Description	Default
data	DataFrame	Input data.	required
duration	str	Follow-up-time column.	required
event	str	Event indicator: 0 = censored, 1, 2, ... = causes.	required
x	sequence of str	Covariate column names.	required
cause	int	Cause of interest (default 1).	1
alpha	float	Significance level for confidence intervals.	0.05
max_iter	(int, float)	Newton-Raphson controls (tol is on the largest coefficient step).	50
tol	(int, float)	Newton-Raphson controls (tol is on the largest coefficient step).	50
vce	('robust', 'model')	"robust" is Fine & Gray's sandwich variance, whose score residuals include the term for estimating the censoring distribution: the var of R cmprsk::crr. "model" is the inverse information of the weighted partial likelihood (crr's invinf), which ignores both the weighting and the estimation of Ĝ and is not a valid variance for this model.	"robust"
small_sample	bool	Multiply the robust variance by n / (n - 1), as Stata's stcrreg does.	False
weights	str	Column of sampling weights (strictly positive), Stata stset [pw=]. Every sum of the fit is weighted: the censoring distribution Ĝ is the weighted Kaplan-Meier, the risk sets and event counts carry the weights, and the robust variance uses the weighted score residuals, including the term for estimating Ĝ. The weights are rescaled to mean one (their scale does not change the estimate or the variance; the log pseudo-likelihood is then Stata's). With small_sample=True coefficients and standard errors equal Stata stcrreg after stset ... [pw=w]. vce='model' is refused with weights (it is not a variance for sampling weights).	None

Returns:

Type	Description
FineGrayResult	With .shr, .tidy(), .summary().

Notes

Subjects who fail from a competing cause are retained in the risk set with time-decaying inverse-probability-of-censoring weights w_i(t) = Ĝ(t-) / Ĝ(T_i-) (Ĝ = KM estimate of the censoring survival, evaluated at left limits). The weighted partial likelihood is maximised by Newton-Raphson with the Breslow tie approximation. Coefficients, both variances and the log pseudo-likelihood match R cmprsk::crr; with small_sample=True the standard errors match Stata stcrreg.

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 300
>>> x = rng.normal(size=n)
>>> t1 = rng.exponential(scale=np.exp(-0.5 * x))   # cause of interest
>>> t2 = rng.exponential(scale=1.5, size=n)        # competing cause
>>> cens = rng.exponential(scale=2.0, size=n)
>>> time = np.minimum(np.minimum(t1, t2), cens)
>>> status = np.where((t1 <= t2) & (t1 <= cens), 1,
...                   np.where((t2 < t1) & (t2 <= cens), 2, 0))
>>> df = pd.DataFrame({"time": time, "status": status, "x": x})
>>> res = sp.finegray(
...     df, duration="time", event="status", x=["x"], cause=1
... )
>>> res.cause
1
>>> res.tidy()["term"].tolist()
['x']
>>> bool(res.shr[0] > 0)  # subdistribution hazard ratio
True

References

fine1999proportional

cox_frailty

cox_frailty(formula: str, data: DataFrame, cluster: str, alpha: float = 0.05, maxiter: int = 50, tol: float = 1e-10, theta: Optional[float] = None, ties: str = 'efron') -> FrailtyResult

Cox proportional hazards with shared gamma frailty.

Parameters:

Name	Type	Description	Default
formula	str	"duration + event ~ x1 + x2" (like R's Surv(time, event) ~ x).	required
data	DataFrame		required
cluster	str	Column identifying clusters (e.g. hospital, site).	required
alpha	float	Significance level (kept for API symmetry; intervals are built by the result's accessors).	0.05
maxiter	(int, float)	Newton-Raphson controls for the penalised fit at each θ.	50
tol	(int, float)	Newton-Raphson controls for the penalised fit at each θ.	50
theta	float	Fix the frailty variance instead of estimating it (R frailty(id, theta=)). 0 fits the ordinary Cox model.	None
ties	('efron', 'breslow')	Tie handling in the partial likelihood (R's default is Efron, Stata's stcox default is Breslow; identical without ties).	"efron"

Returns:

Type	Description
FrailtyResult	beta / se (log hazard ratios; SEs conditional on θ, from the inverse of the full penalised information), theta (frailty variance), frailties (exp of the log-frailties log_frailties), log_likelihood (the integrated log-likelihood maximised over θ) and lr_theta0.

Notes

θ is the maximiser of the integrated likelihood. At a fixed θ the coefficients, standard errors and log-frailties are those of Stata stcox, shared(id) and of R coxph(... + frailty(id, theta=, sparse=FALSE)). R's default sparse=TRUE approximates the frailty block of the information by its diagonal, which moves the standard errors in the fourth digit, and its default method="em" stops the θ search on a 1e-5 relative change in the likelihood, so its reported θ is not the exact maximiser.

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> cluster = np.repeat(np.arange(20), 15)
>>> n = cluster.size
>>> frailty = rng.gamma(2.0, 0.5, 20)[cluster]
>>> x = rng.normal(size=n)
>>> rate = frailty * np.exp(0.7 * x)
>>> T = rng.exponential(1.0 / rate)
>>> C = rng.exponential(2.0, n)
>>> df = pd.DataFrame({"time": np.minimum(T, C),
...                    "event": (T <= C).astype(int),
...                    "x": x, "hospital": cluster})
>>> res = sp.cox_frailty("time + event ~ x", df, cluster="hospital")
>>> res.var_names
['x']
>>> bool(res.theta >= 0)  # gamma frailty variance
True

References

[@therneau2000modeling]

cox

cox(formula: Optional[str] = None, data: DataFrame = None, duration: Optional[str] = None, event: Optional[str] = None, x: Optional[List[str]] = None, ties: str = 'efron', strata: Optional[str] = None, robust: str = 'nonrobust', cluster: Optional[str] = None, hazard_ratio: bool = True, alpha: float = 0.05, weights: Optional[str] = None, start: Optional[str] = None, ssc: Optional[str] = None, weights_as: Optional[str] = None, design: Optional[Any] = None) -> CoxResult

Cox Proportional Hazards model via partial likelihood.

Parameters:

Name	Type	Description	Default
formula	str	Formula of the form 'duration ~ x1 + x2'. If given, duration is inferred from the LHS.	None
data	DataFrame	Input data.	None
duration	str	Column name for follow-up time (overrides formula LHS).	None
event	str	Column name for event indicator (1 = event, 0 = censored).	None
x	list of str	Covariate column names (overrides formula RHS).	None
ties	str	Tie-handling method: 'efron' or 'breslow'.	``'efron'``
strata	str	Column name for stratification variable.	None
robust	str	'hc0' (or 'robust' / True) for the sandwich SE of R coxph(robust = TRUE); with ssc='stata' it carries Stata's N/(N-1). Other values raise ('cluster' only together with cluster=).	``'nonrobust'``
cluster	str	Column name for cluster-robust SE.	None
hazard_ratio	bool	If True, report hazard ratios in the summary alongside coefficients.	True
alpha	float	Significance level for confidence intervals.	0.05
weights	str	Column of non-negative case weights (inverse-probability-of- treatment weights, survey weights, frequency counts). Weighted partial likelihood as in R coxph(weights=): risk sets sum w_k exp(x_k'b) and each Efron term of an event time carries the mean weight of its deaths. The model-based variance reads the weights as frequencies; for sampling / IPT weights pass weights_as='sampling' (or robust='hc0' / cluster=), which is what R reports by default for non-integer weights and Stata stcox [pw=] always does.	None
start	str	Column with the start of each row's interval: counting-process data (start, duration] as in R Surv(start, stop, event) and Stata stset stop, time0(start) / enter(). A row is at risk at t iff start < t <= duration, which covers delayed entry (left truncation) and time-varying covariates split into intervals. With several rows per subject pass cluster=<subject id> for the robust variance by subject (coxph(cluster=id); Stata clusters on the stset id automatically).	None
ssc	(r, stata)	Small-sample factor of the robust and clustered variances. 'r' applies none, as R coxph(robust=TRUE) / coxph(cluster=). 'stata' multiplies the robust variance by N/(N-1) and the clustered one by G/(G-1), as Stata stcox, vce(robust) / vce(cluster) (N counts observations, not weights). Not given: the robust variance has no factor and the clustered one G/(G-1) (the behaviour before this option existed). The model-based variance and the coefficients are unaffected.	'r'
weights_as	(frequency, sampling)	How to read weights=. 'frequency' keeps the model-based variance; 'sampling' (Stata [pw=]) asks for the robust one when neither robust= nor cluster= chose a variance.	'frequency'
design	SurveyDesign	A complex survey design from :func:sp.svydesign declared on data (same rows). The partial likelihood is weighted by the design weights and the covariance is I^-1 M I^-1 with M the survey linearisation covariance of the weighted score residuals over strata, PSUs and the finite population correction -- R survey::svycoxph (Binder 1992). Inference uses t with the design degrees of freedom. With counting-process rows (start=) the subject's rows must lie in one PSU. Not with weights=, cluster=, robust= or ssc=. A replicate-weight design (:func:sp.svrepdesign) is accepted too: the estimator is refitted with every replicate's weights and the replicate variance replaces the linearisation (R svrepdesign / Stata svy brr / svy jackknife).	None

Returns:

Type	Description
CoxResult	Result object extending EconometricResults with .concordance, .baseline_hazard(), .ph_test(), .plot().

Notes

With start= or weights= the concordance is computed over the risk sets of the fit (pairs in which one row fails while the other is at risk, weighted by the product of their weights), as R concordance() does for counting-process and weighted fits; Harrell's pairwise C of right-censored data is the special case without either.

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 150
>>> age = rng.normal(50, 10, n)
>>> treatment = rng.integers(0, 2, n).astype(float)
>>> time = rng.exponential(
...     scale=np.exp(-0.03 * (age - 50) - 0.5 * treatment))
>>> status = (rng.random(n) < 0.8).astype(int)
>>> df = pd.DataFrame({"time": time, "status": status,
...                    "age": age, "treatment": treatment})
>>> res = sp.cox(formula="time ~ age + treatment", data=df,
...              event="status")
>>> bool(list(res.params.index) == ["age", "treatment"])
True
>>> bool(isinstance(res.summary(), str))
True

kaplan_meier

kaplan_meier(data: DataFrame, duration: str, event: str, group: Optional[str] = None, alpha: float = 0.05, conf_type: Optional[str] = 'log-log', weights: Optional[str] = None, start: Optional[str] = None, cluster: Optional[str] = None, se_method: Optional[str] = None, design: Optional[Any] = None) -> KMResult

Kaplan-Meier non-parametric survival function estimator.

Parameters:

Name	Type	Description	Default
data	DataFrame	Input data.	required
duration	str	Column name for duration / follow-up time.	required
event	str	Column name for event indicator (1 = event, 0 = censored).	required
group	str	Column name for group variable (stratification).	None
alpha	float	Significance level for confidence intervals (Greenwood formula).	0.05
conf_type	(log - log, log, plain)	Scale on which the interval is symmetric. "plain" is S(t) +/- z se truncated to [0, 1]; "log" is symmetric for log S(t) (the default of R survival::survfit); "log-log" for log(-log S(t)) (the default of Stata sts list), which always stays inside (0, 1) and has the best coverage in the tails. The estimate and its standard error are the same under all three. .. versionchanged:: 1.40.0 The default is "log-log" (it was "plain"): the plain interval is the one of the three that leaves [0, 1] and has to be truncated, and it is neither reference package's default. Pass conf_type="plain" for the interval of 1.39 and earlier.	"log-log"
weights	str	Column of non-negative case weights (inverse-probability-of- treatment weights give the IPTW-adjusted curve of each group). Deaths and the number at risk are weighted sums, as in R survfit(weights=); n_risk / n_event / n_censor are then weighted and reported as floats.	None
start	str	Column with each row's entry time: counting-process data (start, duration] (delayed entry), as R Surv(start, stop, event). A row is at risk at t iff start < t <= duration.	None
cluster	str	Column whose levels are independent units (subject id with several rows); the robust variance sums the influence within each.	None
se_method	(greenwood, ij, xie_liu)	'greenwood' is Greenwood's formula on the (weighted) counts; 'ij' the infinitesimal jackknife (robust) variance of R survfit(robust=TRUE), which treats the weights as fixed sampling weights; 'xie_liu' Greenwood's formula with the effective number at risk (sum w)^2 / sum w^2 in the denominator, the variance R adjustedCurves reports for IPTW Kaplan-Meier curves (adjustedsurv(method = "iptw_km")). Not given: R survfit's default, 'ij' with cluster= or non-integer weights and 'greenwood' otherwise.	'greenwood'
design	SurveyDesign	A complex survey design from :func:sp.svydesign declared on data (same rows). The curve is weighted by the design weights and std_err is the survey linearisation of the influence of each row on S(t) over strata, PSUs and the finite population correction (R survey::svyrecvar applied to the weighted infinitesimal-jackknife influence that survival's resid(survfit, type = "surv") reports), with intervals on the t distribution with the design degrees of freedom. R svykm(se = TRUE) is not this estimator: its curve differs from the product-limit one. Each group's curve is a domain: its standard errors keep every PSU of the design. Not with weights=, cluster= or se_method=.	None

Returns:

Type	Description
KMResult	Object with .survival_table, .median_survival, .plot(), .summary().

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 150
>>> df = pd.DataFrame({
...     "time": rng.exponential(10, n).round(2),
...     "status": (rng.random(n) < 0.7).astype(int),
... })
>>> km = sp.kaplan_meier(data=df, duration="time", event="status",
...                      conf_type="log-log")
>>> bool("survival" in km.survival_table.columns)
True
>>> bool(isinstance(km.median_survival, float))
True

logrank_test

logrank_test(data: DataFrame, duration: str, event: str, group: str) -> dict

Log-rank test for equality of survival distributions across groups.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
duration	str	Column names.	required
event	str	Column names.	required
group	str	Column names.	required

Returns:

Type	Description
dict	Keys: test_statistic, p_value, df, n_groups, expected_events (per group), observed_events (per group).

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 150
>>> treatment = rng.integers(0, 2, n)
>>> time = rng.exponential(scale=np.exp(-0.5 * treatment))
>>> status = (rng.random(n) < 0.8).astype(int)
>>> df = pd.DataFrame({"time": time, "status": status,
...                    "treatment": treatment})
>>> res = sp.logrank_test(data=df, duration="time", event="status",
...                       group="treatment")
>>> bool("p_value" in res and "test_statistic" in res)
True

References

mantel1959statistical

survreg

survreg(formula: Optional[str] = None, data: DataFrame = None, duration: Optional[str] = None, event: Optional[str] = None, x: Optional[List[str]] = None, dist: str = 'weibull', robust: str = 'nonrobust', cluster: Optional[str] = None, alpha: float = 0.05, frailty: Optional[str] = None) -> EconometricResults

Parametric survival model (AFT parameterization).

The data may also come first: sp.survreg(df, duration=..., event=..., x=[...]) and sp.survreg(df, "t ~ x1 + x2", event=...).

Parameters:

Name	Type	Description	Default
formula	str	Formula 'duration ~ x1 + x2'.	None
data	DataFrame		None
duration	str	Follow-up time column (or formula LHS).	None
event	str	Event indicator column.	None
x	list of str	Covariate columns (or formula RHS).	None
dist	str	Distribution: 'weibull', 'exponential', 'lognormal', 'loglogistic'.	``'weibull'``
robust	str	'nonrobust' (observed information), 'robust' (sandwich, Stata's vce(robust)) or 'cluster'. vce= is an alias.	``'nonrobust'``
cluster	str	Cluster column; implies the cluster sandwich.	None
alpha	float		0.05
frailty	(None, gamma)	Unobserved heterogeneity. 'gamma' multiplies each observation's hazard by a gamma variable with mean one and variance theta, as Stata's streg, frailty(gamma). The survivor function becomes [1 - theta log S(t)]^(-1/theta). log(theta) is reported with the other parameters, and model_info carries theta and the likelihood-ratio test of theta = 0, whose p-value is halved because the null is on the boundary. Ignoring heterogeneity biases a Weibull fit towards negative duration dependence.	None

Returns:

Type

Description

EconometricResults

Fitted parametric survival model. Parameters include covariates and log(sigma) (scale). diagnostics['Log-likelihood'] is the likelihood of the durations, as R survival::survreg reports it; model_info['ll_log_time'] is the likelihood of the log durations, which is what Stata's streg prints (the two differ by the sum of log(t) over the events), and what :func:sp.aft reports for the same fit. model_info ['ll_scale'] names the convention of 'Log-likelihood'.

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 150
>>> age = rng.normal(50, 10, n)
>>> treatment = rng.integers(0, 2, n).astype(float)
>>> time = rng.exponential(
...     scale=np.exp(0.02 * (age - 50) + 0.5 * treatment))
>>> status = (rng.random(n) < 0.8).astype(int)
>>> df = pd.DataFrame({"time": time, "status": status,
...                    "age": age, "treatment": treatment})
>>> res = sp.survreg("time ~ age + treatment", data=df,
...                  event="status", dist="weibull")
>>> bool("log(sigma)" in list(res.params.index))
True