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.

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 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.

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 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.

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 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.

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(), which returns a tidy frame.

cox.to_frame(format="polars")

	term str	estimate f64	std_error f64	statistic f64	p_value f64	conf_low f64	conf_high f64
PolarsRows3Columns7
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).

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"],
})

	Efron f64	Breslow f64
PandasRows2Columns2
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.

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")

	term str	estimate f64	std_error f64	statistic f64	p_value f64	conf_low f64	conf_high f64
PolarsRows3Columns7
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.

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.

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")

	n i64	nevent i64	loglik f64	aic f64	lr_statistic f64	df i64	lr_p_value f64
PolarsRows1Columns7
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.

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.

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"]]

	term str	estimate f64
PolarsRows2Columns2
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.

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.

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

	subject i64	start i64	stop i64	event i64	treated i64
PandasRows8Columns5
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.

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"]]

	term str	estimate f64	p_value f64
PolarsRows1Columns3
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.

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.

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:

# 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

	subject i64	start i64	stop i64	event i64	start_relative i64	stop_relative i64
PandasRows6Columns6
0	1	0	10	0	0	10
1	1	10	25	1	10	25
2	2	365	375	0	0	10
3	2	375	390	1	10	25
4	3	730	740	0	0	10
5	3	740	755	1	10	25

Polars:

import polars as pl

# Example data in calendar time (subjects enroll at different dates)
df = pl.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 (polars)
df = df.with_columns([
    pl.col('start').min().over('subject').alias('entry_date')
]).with_columns([
    (pl.col('start') - pl.col('entry_date')).alias('start_relative'),
    (pl.col('stop') - pl.col('entry_date')).alias('stop_relative')
]).drop('entry_date')

df

	subject i64	start i64	stop i64	event i64	start_relative i64	stop_relative i64
PolarsRows6Columns6
0	1	0	10	0	0	10
1	1	10	25	1	10	25
2	2	365	375	0	0	10
3	2	375	390	1	10	25
4	3	730	740	0	0	10
5	3	740	755	1	10	25

Now you can use start_relative and stop_relative in Surv.counting().

Note: Greenwood will warn if it detects counting-process data where subjects don’t start at 0, as this usually indicates an error in data preparation.

This matters because risk sets are formed by subjects at risk at each event time. Using calendar time instead of subject-relative time creates an imbalance where early times have many subjects at risk and late times have few, distorting the baseline hazard and model coefficients.

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:

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")

	term str	estimate f64	std_error f64	statistic f64	p_value f64	conf_low f64	conf_high f64
PolarsRows2Columns7
0	age	0.0170453318454	0.00948922108421	1.79628356154	0.0724494308764	-0.00155319972097	0.0356438634118
1	sex	-0.513218517108	0.159919372277	-3.20923293909	0.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:

# 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")

	term str	estimate f64	std_error f64	statistic f64	p_value f64	conf_low f64	conf_high f64
PolarsRows2Columns7
0	age	0.0170453318454	0.00950511031412	1.79328080181	0.0729280355773	-0.00158434203935	0.0356750057302
1	sex	-0.513218517108	0.158290649924	-3.24225415308	0.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, installed with the formula extra.

gw.CoxPH().fit(y, "age + sex", data=lung).to_frame(format="polars")[["term", "estimate"]]

	term str	estimate f64
PolarsRows2Columns2
0	age	0.0170453318454
1	sex	-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.

gw.CoxPH().fit(y, "age * sex", data=lung).to_frame(format="polars")[["term", "estimate"]]

	term str	estimate f64
PolarsRows3Columns2
0	age	0.0300312237831
1	sex	0.0953223158055
2	age: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.

Penalized regression

When you have many covariates, or they are collinear, an unpenalized fit can overfit or become unstable. 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.

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 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 accepts a CoxNet directly, so you can compare candidate values by held-out concordance.

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

Standardization

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.

Next steps

You can now fit a Cox model, read hazard ratios, and assess overall fit.

Cox model diagnostics checks the proportional hazards assumption, computes residuals, predicts survival curves, and adds robust variance and stratification.
Parametric survival models offers an alternative modeling approach when you are willing to specify the shape of the survival distribution.
Prediction performance evaluates how well a fitted model discriminates and calibrates.