Cox proportional hazards model.
Usage
The Cox proportional hazards model is the most widely used regression method for survival data. It models the hazard (instantaneous risk of an event) as a multiplicative function of covariates: \(h(t \mid x) = h_0(t) \exp(\beta^\top x)\). The model is semi-parametric: the baseline hazard \(h_0(t)\) is left unspecified (estimated non-parametrically), while covariate effects are estimated parametrically through the log-hazard-ratio coefficients \(\beta\).
To use this model, call fit() with a right-censored or counting-process Surv response and a design matrix of covariates (2-D array or DataFrame). The model automatically handles stratification (via by= in fit), tied event times (via configurable tie-handling methods), and can compute predictions, baseline hazards, and diagnostic residuals. Results include coefficient estimates with confidence intervals, hazard ratios, standard errors, and global significance tests.
The implementation uses maximum partial likelihood to estimate coefficients. Variance estimates use the observed information matrix (Hessian). The model assumes proportional hazards: the ratio of hazards between two subjects remains constant over time. This can be checked using the cox_zph() method for formal tests or diagnostic plots.
ties: str = "efron"Tie-handling method: "efron" (default, as in R) or "breslow".
conf_level: float = 0.95 Fitted estimatorcoef_, hazard_ratio_, std_error_, z_, p_value_, conf_low_, conf_high_, concordance_, lr_stat_, df_), accessible as arrays or exported to DataFrames.
Call fit(surv, covariates) with a Surv response and a design (a 2-D array or a dataframe of covariates). Rows with missing values are dropped (complete-case, as in R’s default na.omit). Results are exposed as arrays (coef_, std_error_, hazard_ratio_, …) and as tidy frames via to_frame() (optionally format=) and greenwood.tidy.
Build a Surv response from the bundled lung dataset and fit the model on age and sex. Printing the fitted object reports the coefficient table and global tests in the style of R’s summary.coxph.
import greenwood as gw
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])
coxCoxPH (Cox proportional hazards model, ties='efron')
coef exp(coef) se(coef) z p
age 0.01705 1.017 0.009223 1.848 0.06459
sex -0.5132 0.5986 0.1675 -3.065 0.002178
n = 228, events = 165
Likelihood ratio test = 14.12 on 2 df, p = 0.0008574Hazard ratios (and their confidence limits) come from tidy with exponentiate=True; pass format= to choose the backend (here, Polars):
| Name | Description |
|---|---|
| baseline_hazard() | Return the uncentered baseline cumulative hazard and survival as a frame. |
| concordance() | Harrell’s concordance index (C-statistic) of the fitted risk scores. |
| cox_zph() | Test the proportional-hazards assumption (Grambsch-Therneau). |
| fit() | Fit the model to a Surv response and a covariate design. |
| predict() | Predict from the fitted model. |
| residuals() | Return diagnostic residuals from the fitted Cox model. |
| to_frame() | Return a tidy coefficient table as a DataFrame (one row per term). |
Return the uncentered baseline cumulative hazard and survival as a frame.
Usage
The baseline hazard represents the hazard rate for a reference subject with all covariates at their mean values. It is useful for understanding the underlying time-to-event distribution estimated by the model, and can be combined with individual covariate values to compute predicted survival probabilities for specific subjects.
In Cox proportional hazards models, the hazard for an individual is modeled as: \(h(t \mid x) = h_0(t) \exp(x^\top \beta)\), where \(h_0(t)\) is the baseline hazard. This method returns the estimated cumulative baseline hazard \(H_0(t)\) at each observed event time, evaluated using the Breslow estimator (non-parametric).
format: str | None = NoneOutput format: None (default), "pandas", "polars", or "pyarrow".
None (default): Auto-detects and tries Polars first, falls back to Pandas, then Pyarrow. Raises an error if no DataFrame library is installed."pandas": returns pandas.DataFrame."polars": returns polars.DataFrame."pyarrow": returns pyarrow.Table. pandas.DataFrame, polars.DataFrame, or pyarrow.TableA DataFrame with one row per event time containing:
time: Event times at which the baseline hazard is evaluated.cumhaz: Cumulative baseline hazard \(H_0(t)\) at each time.survival: Baseline survival probability \(S_0(t) = \exp(-H_0(t))\).strata (if stratified): Stratum label, one baseline hazard per stratum.The baseline hazard is evaluated only at the event times in the training data. The cumulative hazard is non-decreasing by construction. For stratified models, each stratum has its own baseline hazard while coefficients are shared across strata, allowing different baseline risks for different groups.
The baseline survival \(S_0(t)\) is computed from the cumulative hazard using the relationship \(S_0(t) = \exp(-H_0(t))\), consistent with the exponential survival model.
The baseline cumulative hazard (and the implied baseline survival) is reported at every event time. Pass format= to choose the backend (here, Polars):
import greenwood as gw
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])
cox.baseline_hazard(format="polars")| time | cumhaz | survival |
|---|---|---|
| f64 | f64 | f64 |
| 5.0 | 0.002955 | 0.997049 |
| 11.0 | 0.011906 | 0.988164 |
| 12.0 | 0.014928 | 0.985183 |
| 13.0 | 0.021029 | 0.97919 |
| 15.0 | 0.024108 | 0.97618 |
| … | … | … |
| 840.0 | 1.857679 | 0.156034 |
| 883.0 | 2.013058 | 0.13358 |
| 965.0 | 2.013058 | 0.13358 |
| 1010.0 | 2.013058 | 0.13358 |
| 1022.0 | 2.013058 | 0.13358 |
The returned DataFrame shows the estimated hazard and survival trajectory for the reference population (covariates at their means). For stratified models, a separate baseline is provided for each stratum:
cox_stratified = gw.CoxPH().fit(y, lung[["age", "ph.ecog"]], strata=lung["sex"])
cox_stratified.baseline_hazard(format="polars")| time | cumhaz | survival | strata |
|---|---|---|---|
| f64 | f64 | f64 | i64 |
| 11.0 | 0.006805 | 0.993218 | 1 |
| 12.0 | 0.009112 | 0.990929 | 1 |
| 13.0 | 0.01381 | 0.986285 | 1 |
| 15.0 | 0.016196 | 0.983935 | 1 |
| 26.0 | 0.018593 | 0.981579 | 1 |
| … | … | … | … |
| 735.0 | 0.673303 | 0.510021 | 2 |
| 740.0 | 0.673303 | 0.510021 | 2 |
| 765.0 | 0.800708 | 0.449011 | 2 |
| 821.0 | 0.800708 | 0.449011 | 2 |
| 965.0 | 0.800708 | 0.449011 | 2 |
The baseline hazard can be combined with individual predictions to compute personalized survival curves (see predict(type="survival")).
Harrell’s concordance index (C-statistic) of the fitted risk scores.
Usage
The concordance index measures how well the model’s predicted risk scores order subjects by their survival times. It ranges from 0 to 1, where 0.5 indicates predictions are no better than random (coin flip), and 1.0 indicates perfect discrimination (the model always assigns higher risk to subjects who die first).
Pairs of subjects are compared: a subject who experiences an event at time t is considered to have “failed before” another subject still under observation at t (including one censored exactly at t). If the model assigns higher risk to the subject who failed first, the pair is concordant. Ties in predicted risk are treated as half-concordant.
For stratified models, only within-stratum pairs are compared.
floatThe concordance index, a value between 0 and 1. Typical interpretation:
The concordance index is equivalent to the Area Under the Receiver Operating Characteristic curve (AUC) for binary classification problems. It is computed as the fraction of concordant pairs out of all comparable pairs.
Comparable pairs are those where:
Tied event times within the same outcome (both events or both censored at the same time) are excluded from comparison.
Harrell’s C is returned as a single number between 0 and 1. A value of 0.5 means the model is not better than random guessing; 1.0 means perfect discrimination:
Test the proportional-hazards assumption (Grambsch-Therneau).
Usage
The Cox model assumes that the hazard ratio between any two subjects is constant over time (proportional hazards). If this assumption is violated (for example, if a treatment effect diminishes over time) the Cox estimates may be biased. This test checks for violations by regressing scaled Schoenfeld residuals on time.
Large test statistics or small p-values (typically p < 0.05) suggest the proportional-hazards assumption is violated for that covariate. When violated, consider stratified analysis (separate baseline hazards per stratum), time-dependent covariates, or time-varying coefficients.
transform: str = "identity"Transformation to apply to time when computing the test. Options are:
"identity" (default): Use time as-is. Regression on raw time."log": Use log(time). Regression on log-transformed time.cox.zph() (though R defaults to Kaplan-Meier transform; "km" and "rank" are planned).
ZPHResultper_term dict) and a global test (global_test dict) across all covariates. Each includes chi-squared statistic, degrees of freedom, and p-value. Access results via .to_frame() or dictionary keys.
The test uses scaled Schoenfeld residuals, which under the null hypothesis (proportional hazards) have a known asymptotic distribution. The test statistic is approximately chi-squared with 1 df for each term, and chi-squared with degrees of freedom equal to the number of terms for the global test.
Schoenfeld residuals are weighted by the variance-covariance matrix of the risk set at each event time. The regression accounts for the constraint that Schoenfeld residuals sum to zero.
The test returns a ZPHResult summarizing per-term and global p-values:
import greenwood as gw
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])
zph = cox.cox_zph()
zphZPHResult(transform='identity', age: p=0.7065, sex: p=0.0992, GLOBAL p=0.2425)The full statistics are available as a tidy frame, one row per term plus a GLOBAL row; pass format= to choose the backend (here, Polars):
Fit the model to a Surv response and a covariate design.
Usage
covariates is a dataframe or 2-D array, or a right-hand-side formula string (for example "age + sex + C(ph.ecog)") evaluated against data. strata gives per-stratum baseline hazards with shared coefficients. robust=True (or providing cluster ids) reports the Lin-Wei sandwich variance; cluster sums the score residuals within groups before forming the sandwich.
surv: SurvA Surv object representing the response (censoring type must be right-censored or counting-process).
covariates: AnyCovariate design, either a 2-D array, dataframe, or formula string.
data: Any = NoneDataFrame to evaluate formula strings against (required if covariates= is a formula string).
strata: Any = NoneOptional stratification variable, giving each stratum its own baseline hazard while sharing coefficients. Can be a 1-D array or series.
robust: bool = FalseIf True, report Lin-Wei sandwich variance (robust standard errors). Default is False.
cluster: Any = NoneOptional cluster labels for grouped robust variance estimation. Sums score residuals within groups before forming the sandwich.
max_iter: int = 30Maximum number of iterations for the Newton-Raphson solver. Default is 30.
tol: float = 1e-091e-9.
CoxPHcoef_, std_error_, hazard_ratio_, z_, p_value_, and other model diagnostics.
Passing strata= gives each stratum its own baseline hazard while sharing the coefficients. Here we fit age and ph.ecog stratified by sex:
import greenwood as gw
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
gw.CoxPH().fit(y, lung[["age", "ph.ecog"]], strata=lung["sex"]).to_frame(
format="polars"
)| term | estimate | std_error | statistic | p_value | conf_low | conf_high |
|---|---|---|---|---|---|---|
| str | f64 | f64 | f64 | f64 | f64 | f64 |
| "age" | 0.010566 | 0.009241 | 1.143364 | 0.252887 | -0.007547 | 0.028679 |
| "ph.ecog" | 0.462424 | 0.114761 | 4.029453 | 0.000056 | 0.237497 | 0.687352 |
The covariates argument also accepts a right-hand-side formula string (for example "age + sex + C(ph.ecog)"), and robust=True reports the Lin-Wei sandwich variance.
Predict from the fitted model.
Usage
type is one of "lp" (centered linear predictor), "risk" (exp(lp)), or "survival". For "survival", returns a frame of survival probabilities at times (defaulting to the event times), one column per row of newdata. Survival prediction for stratified models is not yet supported.
conditional_after (a scalar or one value per subject) predicts survival conditional on having already survived to that time: the returned value at time \(t\) is \(P(T > t \mid T > c) = S(t) / S(c)\), and is 1 for \(t \le c\).
With ci=True (survival only), the frame also carries _lower and _upper columns per subject: a pointwise confidence band from the cumulative-hazard standard error (the log transform used by R’s survfit), at the model’s conf_level.
newdata: Any = NoneCovariate design for prediction. If None, predictions are made on the fitted data. Can be a 2-D array or dataframe.
type: str = "lp"Type of prediction: "lp" (centered linear predictor, default), "risk" (exp of linear predictor), or "survival" (survival probability).
times: Any = NoneTime points at which to compute survival probabilities (for type="survival"). Defaults to the event times from the fitted model.
conditional_after: Any = NoneOptional scalar or per-subject time for conditional survival prediction. Computes \(P(T > t \mid T > c)\) where \(c\) is the conditional_after time.
ci: bool = FalseIf True (survival only), include confidence intervals (_lower and _upper columns). Default is False.
format: str | None = NoneOutput format (for type="survival" only): None (default), "pandas", "polars", or "pyarrow".
None (default): Auto-detects and tries Polars first, falls back to Pandas, then Pyarrow. Raises an error if no DataFrame library is installed."pandas": returns pandas.DataFrame."polars": returns polars.DataFrame."pyarrow": returns pyarrow.Table. ndarray or DataFrametype="lp" or "risk", returns a 1-D array with one prediction per row. For type="survival", returns a DataFrame with rows for each time point and columns for each subject (named subject_1, subject_2, etc.), optionally with _lower and _upper columns for confidence intervals.
The default type="lp" returns the centered linear predictor as a NumPy array, one value per fitted subject:
import greenwood as gw
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])
cox.predict(type="lp")[:5]array([0.3995047 , 0.2972327 , 0.09268872, 0.10973405, 0.16087005])With type="survival" and newdata, the result is a frame of survival probabilities at the requested times, one column per new subject. Pass format= to choose the backend (here, Polars):
| time | subject_1 | subject_2 | subject_3 |
|---|---|---|---|
| f64 | f64 | f64 | f64 |
| 180.0 | 0.62234 | 0.651705 | 0.705421 |
| 365.0 | 0.268757 | 0.305376 | 0.380304 |
Passing ci=True adds pointwise confidence bands, and conditional_after= gives survival conditional on having already survived to a landmark time.
Return diagnostic residuals from the fitted Cox model.
Usage
Residuals measure the difference between observed events and model predictions, helping diagnose model fit and identify outliers or influential observations. Martingale residuals are individual-level; Schoenfeld residuals are event-level and useful for checking the proportional-hazards assumption. Both types can be visualized against time or other variables to detect systematic deviations.
type: str = "martingale"Type of residuals to return: "martingale" (default) or "schoenfeld".
"martingale": One residual per observation. Ranges from \(-\infty\) to 1. Positive values suggest the model underestimated risk; negative values suggest overestimation. Useful for overall fit assessment."schoenfeld": One row per event with one column per covariate. Useful for checking the proportional-hazards assumption: plot against time to look for trends. Scaled Schoenfeld residuals are used in the cox_zph() test.format: str | None = NoneOutput format (for type="schoenfeld" only): None (default), "pandas", "polars", or "pyarrow".
None (default): Auto-detects and tries Polars first, falls back to Pandas, then Pyarrow. Raises an error if no DataFrame library is installed."pandas": returns pandas.DataFrame."polars": returns polars.DataFrame."pyarrow": returns pyarrow.Table.type="martingale".
ndarray or DataFrametype="martingale": a 1-D array with one residual per observation. For type="schoenfeld": a DataFrame with one row per event and one column per covariate, ordered by stratum and then event time.
Martingale residuals are computed as: \(M_i = \text{event}_i - H_0(t_i) \exp(X_i \beta)\), where \(H_0\) is the baseline cumulative hazard and \(X_i \beta\) is the linear predictor.
Schoenfeld residuals are computed at each event time as \(X_i - \bar{X}\), where \(X_i\) is the covariate vector of the subject with the event and \(\bar{X}\) is the weighted mean covariate vector for the risk set.
Martingale residuals are returned as one value per observation to assess overall model fit. Large negative residuals may indicate overpredicted risk:
import greenwood as gw
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])
cox.residuals("martingale")[:5]array([ 0.00438999, -0.50576203, -3.12981924, 0.53275397, -2.35065745])Schoenfeld residuals are useful for checking the proportional-hazards assumption by plotting against time or other variables:
Return a tidy coefficient table as a DataFrame (one row per term).
Usage
The table contains coefficient estimates, standard errors, test statistics, p-values, and confidence limits. If exponentiate=True, returns hazard ratios instead of log-hazards.
format: str | None = NoneOutput format: None (default), "pandas", "polars", or "pyarrow". When None, a backend is auto-detected (Polars, then Pandas, then PyArrow).
exponentiate: bool = False pandas.DataFrame, polars.DataFrame, or pyarrow.Table ImportErrorFit a Cox model and export the coefficient table as a Polars frame:
import greenwood as gw
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])
cox.to_frame(format="polars")| term | estimate | std_error | statistic | p_value | conf_low | conf_high |
|---|---|---|---|---|---|---|
| str | f64 | f64 | f64 | f64 | f64 | f64 |
| "age" | 0.017045 | 0.009223 | 1.848078 | 0.064591 | -0.001032 | 0.035123 |
| "sex" | -0.513219 | 0.167458 | -3.06476 | 0.002178 | -0.84143 | -0.185007 |
With exponentiate=True, estimates become hazard ratios: