# Cox regression The Kaplan-Meier estimator and the log-rank test describe and compare whole groups. To measure the effect of a continuous predictor, or to adjust for several variables at once, you need a regression model. The Cox proportional hazards model is the most widely used tool for this. It models the hazard, the instantaneous risk of the event among those still at risk, and expresses how covariates multiply that risk without requiring you to specify the shape of the baseline hazard. This page covers fitting the model, reading hazard ratios, and the tests it reports. ``` python import greenwood as gw lung = gw.load_dataset("lung", backend="polars") y = gw.Surv.right(lung["time"], event=(lung["status"] == 2)) y ``` Surv(type=right, n=228, events=165) The response `y` is a [Surv](../reference/Surv.md#greenwood.Surv) object that pairs each subject's follow-up time with an event indicator, marking censored times with a `+`. It records 165 events among 228 subjects, and it is the outcome the Cox model below regresses on the covariates. # The proportional hazards idea The Cox model assumes that each covariate multiplies the hazard by a constant factor that does not change over time. If being in a treatment group halves the hazard at one month, the model assumes it halves the hazard at every month. That multiplicative factor is the hazard ratio, and estimating it is the point of the model. The word "proportional" refers to this constant-factor assumption, which you should check after fitting; see [Cox model diagnostics](cox-diagnostics.md). The great convenience of the Cox model is that it estimates these hazard ratios without assuming any particular form for how the baseline risk changes over time. This is why it is called semiparametric, and why it is the default choice for most analyses. # Fitting a model You fit the model with a [Surv](../reference/Surv.md#greenwood.Surv) response and a set of covariates. Covariates can be a pandas or Polars data frame; numeric columns are used directly and non-numeric columns are turned into indicator variables automatically. Rows with missing covariate values are dropped, as in a standard complete-case analysis. ``` python 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", "ph.ecog"]]) ``` The call returns a fitted [CoxPH](../reference/CoxPH.md#greenwood.CoxPH) estimator, held in `cox`. Printing it gives a compact summary modeled on R's `coxph` output: the coefficient table, the sample size and event count, and the overall likelihood ratio test. ``` python cox ``` CoxPH (Cox proportional hazards model, ties='efron') coef exp(coef) se(coef) z p age 0.01107 1.011 0.009267 1.194 0.2324 sex -0.5526 0.5754 0.1677 -3.294 0.0009861 ph.ecog 0.4637 1.59 0.1136 4.083 4.447e-05 n = 227, events = 164 Likelihood ratio test = 30.5 on 3 df, p = 1.083e-06 The printed summary is meant for reading. When you want the numbers as data, to filter, join, or plot them, ask for the coefficient table with [to_frame()](../reference/AFT.md#greenwood.AFT.to_frame), which returns a tidy frame. ``` python cox.to_frame(format="polars") ```
PolarsRows3Columns7
term str
estimate f64
std_error f64
statistic f64
p_value f64
conf_low f64
conf_high f64
0 age 0.0110667645601 0.0092674110137 1.19415924725 0.232415681 -0.00709702725671 0.0292305563769
1 sex -0.552612395704 0.167739053787 -3.29447664826 0.000986051372138 -0.881374899927 -0.22384989148
2 ph.ecog 0.46372847537 0.113577266162 4.0829339448 4.44706665186e-05 0.24112112423 0.68633582651
The coefficient table reports, for each covariate, the estimated log hazard ratio (`estimate`), its standard error, a Wald z-statistic and p-value, and a confidence interval. The coefficients are on the log scale, which is convenient for the arithmetic of the model but not for interpretation. ## Handling ties in event times When multiple subjects experience events at the same time, the Cox model needs a method to handle these "ties". Greenwood supports two tie-handling methods, controlled by the `ties` parameter in the constructor: - **"efron"** (default): Efron's method. Recommended; matches R's survival package default. - **"breslow"**: Breslow's method. Computationally simpler but less accurate when ties are common. The choice usually has minimal impact on results unless ties are very common. When in doubt, stick with the default (Efron). ``` python import pandas as pd # Compare methods cox_efron = gw.CoxPH(ties="efron").fit(y, lung[["age", "sex"]]) cox_breslow = gw.CoxPH(ties="breslow").fit(y, lung[["age", "sex"]]) # Coefficients are nearly identical unless ties are very common pd.DataFrame({ "Efron": cox_efron.to_frame(format="polars")["estimate"], "Breslow": cox_breslow.to_frame(format="polars")["estimate"], }) ```
PandasRows2Columns2
Efron f64
Breslow f64
0 0.0170453318454 0.0170128891984
1 -0.513218517108 -0.512564791519
# Hazard ratios and their interpretation To interpret the model you exponentiate the coefficients, which turns log hazard ratios into hazard ratios. Greenwood does this for you through the tidy layer. ``` python 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", "ph.ecog"]]) gw.tidy(cox, exponentiate=True, format="polars") ```
PolarsRows3Columns7
term str
estimate f64
std_error f64
statistic f64
p_value f64
conf_low f64
conf_high f64
0 age 1.01112822772 0.0092674110137 1.19415924725 0.232415681 0.99292809717 1.02966196224
1 sex 0.57544455619 0.167739053787 -3.29447664826 0.000986051372138 0.414213018549 0.799435127387
2 ph.ecog 1.58999119005 0.113577266162 4.0829339448 4.44706665186e-05 1.27267517777 1.98642358129
A hazard ratio above 1 means the covariate increases the hazard, and a value below 1 means it decreases it. For a continuous covariate such as `age`, the hazard ratio is the multiplicative change in hazard per one-unit increase; a hazard ratio of 1.02 for age means roughly a 2 percent higher hazard for each additional year. For the `sex` indicator, coded 1 and 2 in this dataset, a hazard ratio below 1 means the higher-coded group has lower risk. > **Note: Hazard ratio, not risk ratio** > > A hazard ratio compares instantaneous rates among those still at risk, not the probability of the event over the whole study. The two are related but not identical, and the difference matters when events are common. Report hazard ratios as such, and consider an absolute measure like a survival difference or a restricted mean difference alongside them. # Model-fit statistics The `glance` view gives one-row summary statistics for the whole model, including the log-likelihood, the AIC, and the likelihood-ratio test of the null hypothesis that all coefficients are zero. ``` python 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", "ph.ecog"]]) gw.glance(cox, format="polars") ```
PolarsRows1Columns7
n i64
nevent i64
loglik f64
aic f64
lr_statistic f64
df i64
lr_p_value f64
0 227 164 -729.230121375 1464.46024275 30.5006687732 3 1.0828176992e-06
The model reports three classical global tests, all of which assess whether the covariates jointly improve the fit: the likelihood-ratio test, the Wald test, and the score test. They usually agree closely, and the likelihood-ratio test is generally preferred. ``` python print("likelihood ratio:", round(cox.lr_stat_, 3)) print("Wald:", round(cox.wald_stat_, 3)) print("score:", round(cox.score_stat_, 3)) ``` likelihood ratio: 30.501 Wald: 29.929 score: 30.5 # Handling tied event times When two or more events occur at exactly the same recorded time, the partial likelihood must account for the tie. Greenwood defaults to the Efron approximation, which is accurate and is also the default in R. The Breslow approximation is available and is faster but slightly less accurate when ties are common. ``` python lung = gw.load_dataset("lung", backend="polars") y = gw.Surv.right(lung["time"], event=(lung["status"] == 2)) gw.CoxPH(ties="breslow").fit(y, lung[["age", "sex"]]).to_frame(format="polars")[["term", "estimate"]] ```
PolarsRows2Columns2
term str
estimate f64
0 age 0.0170128891984
1 sex -0.512564791519
Unless you have a specific reason to match another tool's Breslow output, the Efron default is the right choice. > **Tip: Flexibility in tie-handling methods** > > Greenwood's support for both Efron and Breslow tie-handling methods gives you flexibility that is useful in several situations: > > - matching textbook examples: if you're learning from a textbook that uses SAS (which defaults to Breslow), you can use `ties="breslow"` to get results that match exactly. > - comparing across tools: if you're validating results across R, SAS, and Python, you can switch the tie method in Greenwood to match any tool's output with a single parameter change, rather than rewriting your analysis in a different language. > - sensitivity analysis: you can easily compare whether your results change substantially between the two methods, which is a good diagnostic when ties are very common. > > Most analyses will stick with the default Efron method (which matches R), but this flexibility means you're never locked into one approach. # Time-varying covariates Some covariates change during follow-up: a treatment that starts partway through, a lab value that is remeasured, or a status that switches once. The Cox model handles these through the counting-process form, where each subject contributes one row per interval over which their covariates are constant. Each row records the interval `(start, stop]`, whether the event happened at its end, and the covariate values that held during it. We build a small illustrative dataset where a treatment switches on for some subjects. The response is a `gw.Surv.counting` object rather than `gw.Surv.right`. ``` python intervals = pd.DataFrame( { "subject": [1, 1, 2, 3, 3, 4, 5, 5], "start": [0, 4, 0, 0, 5, 0, 0, 3], "stop": [4, 10, 7, 5, 14, 9, 3, 12], "event": [0, 1, 1, 0, 0, 1, 0, 1], "treated": [0, 1, 0, 0, 1, 0, 0, 1], } ) intervals ```
PandasRows8Columns5
subject i64
start i64
stop i64
event i64
treated i64
0 1 0 4 0 0
1 1 4 10 1 1
2 2 0 7 1 0
3 3 0 5 0 0
4 3 5 14 0 1
5 4 0 9 1 0
6 5 0 3 0 0
7 5 3 12 1 1
Subject 1, for example, is untreated over `(0, 4]` and treated over `(4, 10]`, with the event at day 10. We pass the interval endpoints to `gw.Surv.counting` and fit as usual. ``` python y_tv = gw.Surv.counting(start=intervals["start"], stop=intervals["stop"], event=intervals["event"]) gw.CoxPH().fit(y_tv, intervals[["treated"]]).to_frame(format="polars")[["term", "estimate", "p_value"]] ```
PolarsRows1Columns3
term str
estimate f64
p_value f64
0 treated -22.3171831456 0.999461642315
The risk set at each event time correctly includes only the intervals that span it, using each subject's covariate values as of that moment. This is exactly R's start-stop `coxph`, and Greenwood matches it to tolerance. > **Note: One row per interval, not per subject** > > The only change from an ordinary Cox fit is the data layout: a subject with a covariate that changes `k` times contributes `k + 1` rows, and the event indicator is 1 only on the interval where the event occurred. Left truncation and delayed entry use the same counting-process response. > **Important: Subject timelines must start at `0`** > > Each subject needs their own timeline starting at 0 (subject-relative time), not calendar time. > > **Correct** (subject-relative): Each subject's earliest interval starts at 0 > > - Subject 1: (0, 4\], (4, 10\] (personal follow-up time) > - Subject 2: (0, 7\] (personal follow-up time) > > **Incorrect** (calendar time): Subjects enter the study at different dates > > - Subject 1: (2024-01-01, 2024-01-04\], (2024-01-04, 2024-01-10\] > - Subject 2: (2024-06-15, 2024-06-22\] (different calendar dates) > > If your data uses calendar time (e.g., from a wide dataset where subjects enroll on different dates), subtract each subject's entry date from their start/stop times. > > **Pandas**: > >
> > ``` python > # Example data in calendar time (subjects enroll at different dates) > df = pd.DataFrame({ > 'subject': [1, 1, 2, 2, 3, 3], > 'start': [0, 10, 365, 375, 730, 740], # Different enrollment dates > 'stop': [10, 25, 375, 390, 740, 755], > 'event': [0, 1, 0, 1, 0, 1] > }) > > # Convert calendar time to subject-relative time (pandas) > df['entry_date'] = df.groupby('subject')['start'].transform('min') > df['start_relative'] = df['start'] - df['entry_date'] > df['stop_relative'] = df['stop'] - df['entry_date'] > > # Drop helper column > df = df.drop(columns=['entry_date']) > df > ``` > >
> >
> >
> > PandasRows6Columns6 > >
> > # Robust and clustered standard errors > > By default, the Cox model uses the model-based (Fisher information) standard errors. When the model may be misspecified, or when observations are not independent (e.g., repeated measures, family data), use the sandwich (robust) variance estimator by setting `robust=True`: > >
> > ``` python > lung = gw.load_dataset("lung", backend="polars") > > y = gw.Surv.right(lung["time"], event=(lung["status"] == 2)) > > cox_robust = gw.CoxPH().fit(y, lung[["age", "sex"]], robust=True) > > gw.tidy(cox_robust, format="polars") > ``` > >
> >
> > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > >
> PolarsRows2Columns7 >
>
> term >
>
> str >
>
>
> estimate >
>
> f64 >
>
>
> std_error >
>
> f64 >
>
>
> statistic >
>
> f64 >
>
>
> p_value >
>
> f64 >
>
>
> conf_low >
>
> f64 >
>
>
> conf_high >
>
> f64 >
>
0age0.01704533184540.009489221084211.796283561540.0724494308764-0.001553199720970.0356438634118
1sex-0.5132185171080.159919372277-3.209232939090.00133089620495-0.826654727201-0.199782307016
> >
> >
> >
> > The standard errors adjust upward (more conservative) if there is overdispersion, or downward if there is underdispersion. For clustered data (e.g., family members, multiple measurements per subject), pass the cluster variable to ensure standard errors account for within-cluster dependence: > >
> > ``` python > # Simulated example: suppose subjects cluster by family > lung_with_cluster = lung.with_columns( > family=((pl.int_range(len(lung)) % 10) + 1) > ) > > cox_cluster = gw.CoxPH().fit(y, lung_with_cluster[["age", "sex"]], > robust=True, cluster=lung_with_cluster["family"]) > > gw.tidy(cox_cluster, format="polars") > ``` > >
> >
> > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > >
> PolarsRows2Columns7 >
>
> term >
>
> str >
>
>
> estimate >
>
> f64 >
>
>
> std_error >
>
> f64 >
>
>
> statistic >
>
> f64 >
>
>
> p_value >
>
> f64 >
>
>
> conf_low >
>
> f64 >
>
>
> conf_high >
>
> f64 >
>
0age0.01704533184540.009505110314121.793280801810.0729280355773-0.001584342039350.0356750057302
1sex-0.5132185171080.158290649924-3.242254153080.00118588214896-0.823462490048-0.202974544168
> >
> >
> >
> > Clustered sandwich standard errors are wider than model-based ones, reflecting the loss of information from within-cluster correlation. > > Instead of selecting columns yourself, you can describe the model with a formula, passing the right-hand side as a string and the data frame as `data`. This is convenient for categorical variables, interactions, and transformations, and it mirrors the notation used in R. The formula support uses [formulaic](https://matthewwardrop.github.io/formulaic/), installed with the `formula` extra. > >
> > ``` python > gw.CoxPH().fit(y, "age + sex", data=lung).to_frame(format="polars")[["term", "estimate"]] > ``` > >
> >
> > > > > > > > > > > > > > > > > > > > > > > > >
> PolarsRows2Columns2 >
>
> term >
>
> str >
>
>
> estimate >
>
> f64 >
>
0age0.0170453318454
1sex-0.513218517108
> >
> >
> >
> > Categorical columns are expanded into indicator terms automatically, and interactions are written with `*` (which expands to the main effects plus their product) or `:` (the product alone). A term like `C(ph.ecog)` forces a column to be treated as categorical. > >
> > ``` python > gw.CoxPH().fit(y, "age * sex", data=lung).to_frame(format="polars")[["term", "estimate"]] > ``` > >
> >
> > > > > > > > > > > > > > > > > > > > > > > > > > > > > >
> PolarsRows3Columns2 >
>
> term >
>
> str >
>
>
> estimate >
>
> f64 >
>
0age0.0300312237831
1sex0.0953223158055
2age:sex-0.00968855903869
> >
> >
> >
> > Rows with a missing value in any formula term are dropped, the same complete-case rule used when you pass columns directly. The formula interface is also available on [AFT](../reference/AFT.md#greenwood.AFT). > >
> > # Penalized regression > > When you have many covariates, or they are collinear, an unpenalized fit can overfit or become unstable. [CoxNet](../reference/CoxNet.md#greenwood.CoxNet) fits an elastic-net penalized Cox model: it shrinks the coefficients and, for the lasso, sets some of them to exactly zero, which selects variables. The `penalizer` argument sets the overall strength, and `l1_ratio` mixes the lasso (`1.0`) and ridge (`0.0`) penalties. > >
> > ``` python > lung = gw.load_dataset("lung", backend="polars") > > y = gw.Surv.right(lung["time"], event=(lung["status"] == 2)) > > cols = ["age", "sex", "ph.ecog", "ph.karno", "wt.loss"] > > lasso = gw.CoxNet(penalizer=0.05, l1_ratio=1.0).fit(y, lung[cols]) > lasso > ``` > >
> > CoxNet (elastic-net Cox, penalizer=0.05, l1_ratio=1.0) > > coef > age 0.006174 > sex -0.4093 > ph.ecog 0.3673 > ph.karno 0 > wt.loss -0.0006976 > > n = 213, events = 151, nonzero coefficients = 4 > >
> >
> > The printed summary reports how many coefficients survived: the lasso has driven the weakest ones to exactly zero, leaving a smaller model. A pure ridge penalty (`l1_ratio=0.0`) instead shrinks every coefficient toward zero without removing any. Because penalized estimates are biased on purpose, [CoxNet](../reference/CoxNet.md#greenwood.CoxNet) reports coefficients for prediction and selection but not p-values, and with `penalizer=0` it reduces to the ordinary Cox fit. > > The penalty strength is usually chosen by cross-validation. [cross_validate](../reference/cross_validate.md#greenwood.cross_validate) accepts a [CoxNet](../reference/CoxNet.md#greenwood.CoxNet) directly, so you can compare candidate values by held-out concordance. > >
> > ``` python > for lam in [0.02, 0.05, 0.1]: > cv = gw.cross_validate(gw.CoxNet(penalizer=lam, l1_ratio=1.0), y, lung[cols], k=5, seed=1) > print(f"penalizer={lam}: CV concordance = {cv['mean']:.4f}") > ``` > >
> > penalizer=0.02: CV concordance = 0.6175 > penalizer=0.05: CV concordance = 0.6216 > penalizer=0.1: CV concordance = 0.6096 > >
> >
> > > **Note: Standardization** > > > > [CoxNet](../reference/CoxNet.md#greenwood.CoxNet) standardizes covariates to unit variance before applying the penalty, as glmnet does, so that the penalty treats variables on comparable scales. Coefficients are returned on the original scale, so you read and use them exactly as you would from [CoxPH](../reference/CoxPH.md#greenwood.CoxPH). > >
> > # Next steps > > You can now fit a Cox model, read hazard ratios, and assess overall fit. > > - [Cox model diagnostics](cox-diagnostics.md) checks the proportional hazards assumption, computes residuals, predicts survival curves, and adds robust variance and stratification. > - [Parametric survival models](parametric-models.md) offers an alternative modeling approach when you are willing to specify the shape of the survival distribution. > - [Prediction performance](prediction-performance.md) evaluates how well a fitted model discriminates and calibrates. > >
> >
> >
> >
> >
> >
> >