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

# Load the bundled lung dataset as a Polars DataFrame
lung = gw.load_dataset("lung", backend="polars")

# Build a right-censored response; status == 2 marks a death
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

# Display the response summary
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 tick. 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.

# Load data and build the response
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

# Fit a Cox model with age, sex, and ECOG score as covariates
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.

# Print the coefficient table and global test statistics
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.

# Export the coefficient table as a tidy Polars DataFrame
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).

import pandas as pd

# Compare Efron and Breslow tie-handling methods on the same data
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

Covariate scaling

The Cox model’s Newton-Raphson solver computes second derivatives of the log partial likelihood. When covariates differ wildly in scale (for example, a binary indicator with values 0/1 alongside an income variable measured in tens of thousands) the Hessian becomes ill-conditioned and the solver may converge slowly, converge to a poor solution, or produce unreliable standard errors.

Greenwood detects this automatically. If the ratio of the largest to smallest column standard deviation exceeds 100, CoxPH.fit() emits a warning:

UserWarning: Covariates appear to be on very different scales
(max/min standard deviation ratio: 1500). Consider standardizing
covariates before fitting ...

The fix is to standardize the offending columns before fitting. In machine learning this operation is called z-score standardization (or feature standardization): subtract the mean (centering) and divide by the standard deviation (scaling) so that each column has mean 0 and standard deviation 1:

import numpy as np

lung = gw.load_dataset("lung", backend="pandas")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

cols = lung[["age", "sex"]].copy()

# Standardize to zero mean, unit variance
cols_std = (cols - cols.mean()) / cols.std()

gw.CoxPH().fit(y, cols_std).to_frame(format="polars")[["term", "estimate", "conf_low", "conf_high"]]

PolarsRows2Columns4
	term str	estimate f64	conf_low f64	conf_high f64
0	age	0.154660078279	-0.00936337155815	0.318683528116
1	sex	-0.251410124615	-0.41219098161	-0.0906292676199

Standardized coefficients are on a different scale than the original: each coefficient represents the log hazard ratio for a one-standard-deviation increase in the covariate. This makes the coefficients harder to interpret directly but often makes the model fit more stably, especially when variables span very different numeric ranges.

TipPenalized Cox automatically standardizes

CoxNet (the elastic-net Cox model) always standardizes covariates internally before applying the penalty and returns coefficients on the original scale, so you do not need to standardize by hand when using penalized regression.

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.

# Load data and build the response
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

# Fit a Cox model with age, sex, and ECOG score as covariates
cox = gw.CoxPH().fit(y, lung[["age", "sex", "ph.ecog"]])

# Export hazard ratios (exponentiated coefficients) with confidence intervals
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% 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.

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

# Load data and build the response
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

# Fit a Cox model with age, sex, and ECOG score as covariates
cox = gw.CoxPH().fit(y, lung[["age", "sex", "ph.ecog"]])

# Retrieve the one-row model summary: n, events, log-likelihood, AIC
gw.glance(cox, format="polars")

PolarsRows1Columns10
	n i64	nevent i64	loglik f64	aic f64	lr_statistic f64	df i64	lr_p_value f64	frailty_theta null	frailty_lrt_statistic null	frailty_lrt_p_value null
0	227	164	-729.230121375	1464.46024275	30.5006687732	3	1.0828176992e-06	None	None	None

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.

# Compare the three classical global tests (they should be similar)
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.

# Load data and build the response
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

# Fit with Breslow tie handling and show just the term and estimate columns
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.

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

WarningImmortal time bias

A common mistake in observational studies is to assign treatment status based on what happens after study entry. For example: classifying patients as “treated” if they ever received a drug, and then analysing survival from enrolment. The time between enrolment and when they actually started the drug is immortal time. Those patients were guaranteed to be alive during that interval (because they had to survive long enough to receive the drug), yet the time is credited to the treated group. This artificially inflates apparent survival for the treated group.

The fix is to use the counting-process form. Represent each patient as untreated from enrolment to the moment they start treatment, and treated from that moment onward. This way the “treated” rows only enter the risk set after treatment actually begins, and the immortal time is correctly attributed to the untreated period.

# Wrong: classify as treated from day 0 even though treatment starts at day 30
df["treated"] = 1   # all time counted as exposed

# Correct: two rows per patient who switches
# row 1: (0, 30], event=0, treated=0     - untreated until drug starts
# row 2: (30, t], event=event, treated=1 - treated from day 30 onward

The example below shows the correct layout.

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.

# Build the time-varying covariate dataset: one row per (subject, interval)
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],
    }
)

# Preview the interval layout
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.

# Build the counting-process response from the interval endpoints
y_tv = gw.Surv.counting(
    start=intervals["start"], stop=intervals["stop"], event=intervals["event"]
)

# Fit the Cox model on the time-varying covariate layout
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.

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

ImportantSubject 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

PandasRows6Columns6
	subject i64	start i64	stop i64	event i64	start_relative i64	stop_relative i64
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

PolarsRows6Columns6
	subject i64	start i64	stop i64	event i64	start_relative i64	stop_relative i64
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:

# Load data and build the response
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

# Fit with the Lin-Wei sandwich variance for robust standard errors
cox_robust = gw.CoxPH().fit(y, lung[["age", "sex"]], robust=True)

# Display the coefficient table with robust standard errors
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
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:

# Simulate a cluster variable: assign subjects to families of 10
lung_with_cluster = lung.with_columns(family=((pl.int_range(len(lung)) % 10) + 1))

# Fit with clustered sandwich standard errors
cox_cluster = gw.CoxPH().fit(
    y,
    lung_with_cluster[["age", "sex"]],
    robust=True,
    cluster=lung_with_cluster["family"],
)

# Display the coefficient table; SEs reflect within-cluster correlation
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
0	age	0.0170453318454	0.00950511031412	1.79328080181	0.0729280355774	-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.

Shared frailty for clustered hazards

Sometimes cluster correlation is not just a variance issue. If subjects in the same cluster (hospital, family, center, repeated-event unit) share unmeasured risk factors, a shared-frailty model adds a cluster-level random effect directly to the hazard rather than only correcting standard errors.

Greenwood supports two shared-frailty distributions via the frailty= argument: "gamma" and "lognormal". Both require right-censored data with Breslow ties and a frailty_cluster= label.

Gamma frailty

The gamma frailty model places a multiplicative cluster-level random effect on the hazard: h(t \mid x, z) = z \cdot h_0(t) \exp(\beta^\top x) with z \sim \text{Gamma}(1/\theta, 1/\theta). The variance parameter \theta is estimated via an EM algorithm. The frailty_effect_ values are on the multiplicative scale and are always positive.

# Load data and build the response
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

# Shared gamma frailty by institution (inst)
cox_frailty = gw.CoxPH(ties="breslow").fit(
    y,
    lung[["age", "sex"]],
    frailty="gamma",
    frailty_cluster=lung["inst"],
)

# Printed summary includes frailty theta and the variance test
cox_frailty

CoxPH (Cox proportional hazards model, ties='breslow')

       coef  exp(coef)  se(coef)       z         p
age   0.017      1.017  0.009231   1.842   0.06554
sex  -0.511     0.5999    0.1677  -3.047  0.002308

n = 227, events = 164
Likelihood ratio test = 13.94 on 2 df, p = 0.0009402
Shared frailty: gamma (theta = 1.378e-08)
Frailty variance test (theta = 0): LR = 0.0004377, p = 0.4917

To get the frailty-variance inference directly, call frailty_test():

cox_frailty.frailty_test()

{'theta': 1.3779723598335543e-08,
 'lr_statistic': 0.00043770764023065567,
 'df': 1.0,
 'p_value': 0.49165415243657906}

You can also retrieve these fields in glance() output:

gw.glance(cox_frailty, format="polars")[
    ["frailty_theta", "frailty_lrt_statistic", "frailty_lrt_p_value"]
]

PolarsRows1Columns3
	frailty_theta f64	frailty_lrt_statistic f64	frailty_lrt_p_value f64
0	1.37797235983e-08	0.000437707640231	0.491654152437

Interpretation:

• frailty_theta is the estimated cluster-level heterogeneity variance.
• larger values indicate stronger between-cluster hazard heterogeneity.
• the LR test evaluates the null hypothesis theta = 0 (no frailty variance).

Log-normal frailty

The log-normal frailty model adds a normally distributed random effect on the log-hazard scale: h(t \mid x, u) = h_0(t) \exp(\beta^\top x + u) with u \sim \mathcal{N}(0, \sigma^2). Joint optimization over (\beta, u) uses penalized partial likelihood; \sigma^2 is updated via a REML step. The frailty_effect_ values are on the log-hazard scale and can be negative.

# Log-normal shared frailty by institution; effects are on the log-hazard scale
cox_ln = gw.CoxPH(ties="breslow").fit(
    y,
    lung[["age", "sex"]],
    frailty="lognormal",
    frailty_cluster=lung["inst"],
)

# Summary reports sigma2 (variance of the normal random effect) instead of theta
cox_ln

CoxPH (Cox proportional hazards model, ties='breslow')

        coef  exp(coef)  se(coef)       z        p
age  0.01726      1.017   0.00928   1.860  0.06282
sex  -0.5118     0.5994    0.1682  -3.043  0.00234

n = 227, events = 164
Likelihood ratio test = 15.46 on 2 df, p = 0.0004397
Shared frailty: lognormal (sigma2 = 0.01426)
Frailty variance test (sigma2 = 0): LR = 0.05126, p = 0.4104

cox_ln.frailty_test()

{'theta': 0.014255476590025335,
 'lr_statistic': 0.05125652908691336,
 'df': 1.0,
 'p_value': 0.410445534017966}

Choose gamma when the frailty is best understood as a multiplicative hazard multiplier. Choose log-normal when an additive model on the log-hazard scale is more natural or when you need random effects that can shrink below zero on the log scale.

NoteFrailty vs clustered robust SEs

• cluster= / robust=True keeps a standard Cox mean structure and adjusts uncertainty.
• frailty="gamma" or frailty="lognormal" changes the model itself by adding cluster-level random effects.

Use frailty when latent cluster heterogeneity is scientifically part of the estimand, not only a nuisance for inference.

ImportantCurrent frailty scope

Both frailty distributions currently require:

• right-censored Surv.right(...) responses,
• ties="breslow",
• no strata=, and
• no simultaneous robust/cluster sandwich variance options.

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.

# Fit using a Wilkinson formula string; categorical columns are dummy-coded automatically
gw.CoxPH().fit(y, "age + sex", data=lung).to_frame(format="polars")[["term", "estimate"]]

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

# The * operator expands to main effects plus their interaction term
gw.CoxPH().fit(y, "age * sex", data=lung).to_frame(format="polars")[["term", "estimate"]]

PolarsRows3Columns2
	term str	estimate f64
0	age	0.0300312237831
1	sex	0.0953223158055
2	age:sex	-0.00968855903869

Controlling the reference category

By default, formulaic picks the first level (alphabetically) as the reference. To choose a specific baseline, pass contr.treatment(base=...) inside C():

# "female" is the reference and the model reports the HR for "male" vs. "female"
gw.CoxPH().fit(y, "C(sex, contr.treatment(base='female')) + age", data=df)

Important

Do not encode unordered categories as integers. Passing a column of 0/1/2/… values to the model treats the integers as a continuous covariate, implying that category 2 has exactly twice the log-hazard effect of category 1. Use C(col) (or cast the column to a string/category dtype before fitting) so that each level gets its own indicator term.

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.

# Load data and build the response
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

# Select the covariates to include in the penalized model
cols = ["age", "sex", "ph.ecog", "ph.karno", "wt.loss"]

# Fit a lasso Cox model (l1_ratio=1.0 drives weak coefficients to zero)
lasso = gw.CoxNet(penalizer=0.05, l1_ratio=1.0).fit(y, lung[cols])

# Print the surviving (non-zero) coefficients
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.

Choosing the penalty by cross-validation

Rather than guessing a penalizer value, you can use cv_coxnet() to search a full log-spaced grid automatically. It runs stratified k-fold cross-validation across all candidate values and returns the one with the best held-out performance.

# Search a 50-value log-spaced lambda path with 5-fold CV (lasso, concordance metric)
cv_result = gw.cv_coxnet(y, lung[cols], l1_ratio=1.0, seed=23)
cv_result

CoxNetCV (concordance, ↑ higher is better, l1_ratio=1.0, 5-fold)

  best penalizer : 0.10767  (mean concordance: 0.6356)
  1-SE  penalizer : 0.16435  (mean concordance: 0.6205)

  50 penalizers tested, range [0.000218, 0.218]

The summary shows two selected penalizers. The best penalizer maximises the mean concordance across folds. The 1-SE penalizer is the largest value (most regularized, fewest non-zero coefficients) whose mean concordance falls within one standard error of the best (following the 1-SE rule). It often produces a sparser, more interpretable model at no meaningful loss of predictive accuracy.

Inspect the full path as a tidy DataFrame to see how the score and sparsity evolve:

# All 50 tested values: penalizer, mean CV score, its SD, and mean non-zero count
cv_result.to_frame(format="polars")

PolarsRows50Columns4
	penalizer f64	mean_score f64	std_score f64	n_nonzero f64
0	0.217881959052	0.533798259659	0.0324174616548	0.6
1	0.189232959572	0.592455616024	0.0447627069837	1.2
2	0.164350977677	0.620501014326	0.0503977774928	1.8
3	0.1427406934	0.629349328932	0.041545668588	2
4	0.123971915715	0.629349328932	0.041545668588	2
5	0.107671018824	0.635581263659	0.0389162401061	2.2
6	0.093513504472	0.629311167196	0.0410072880256	2.4
7	0.0812175422333	0.634739395274	0.048139137214	2.8
47	0.000288848589473	0.61803357802	0.036493617691	5
48	0.000250868285249	0.61803357802	0.036493617691	5
49	0.000217881959052	0.61803357802	0.036493617691	5

Once you have selected a penalizer, fit the final model on all of the data:

# Fit a final lasso model at the best CV penalizer
final_lasso = gw.CoxNet(penalizer=cv_result.best_penalizer_, l1_ratio=1.0).fit(y, lung[cols])
final_lasso

CoxNet (elastic-net Cox, penalizer=0.10767101882396805, l1_ratio=1.0)

             coef
age             0
sex       -0.2297
ph.ecog    0.2498
ph.karno       -0
wt.loss        -0

n = 213, events = 151, nonzero coefficients = 2

To prefer parsimony, use the 1-SE penalizer instead:

# Sparser model: same l1_ratio, stricter regularisation
sparse_lasso = gw.CoxNet(penalizer=cv_result.penalizer_1se_, l1_ratio=1.0).fit(y, lung[cols])
sparse_lasso

CoxNet (elastic-net Cox, penalizer=0.16435097767696763, l1_ratio=1.0)

              coef
age              0
sex       -0.05629
ph.ecog     0.1203
ph.karno        -0
wt.loss         -0

n = 213, events = 151, nonzero coefficients = 2

TipThe 1-SE rule

The best_penalizer_ minimises the CV error. penalizer_1se_ is the largest penalizer whose mean score is within one standard error of the minimum (the 1-SE rule). It gives a sparser, more interpretable model whose performance is statistically indistinguishable from the best in terms of fold-to-fold variability. When in doubt, prefer penalizer_1se_ for final reporting.

The cv_coxnet() function also supports the Brier score metric for calibration-focused tuning. Pass metric="brier" together with evaluation time points:

# Tune for calibration (integrated Brier score) instead of discrimination
cv_brier = gw.cv_coxnet(
    y,
    lung[cols],
    l1_ratio=1.0,
    metric="brier",
    times=[180, 365, 540],
    seed=23,
)
cv_brier

CoxNetCV (brier, ↓ lower is better, l1_ratio=1.0, 5-fold)

  best penalizer : 0.053208  (mean brier: 0.2066)
  1-SE  penalizer : 0.18923  (mean brier: 0.2193)

  50 penalizers tested, range [0.000218, 0.218]

NoteStandardization

CoxNet standardizes covariates to unit variance before applying the penalty 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.

Baseline hazard with confidence intervals

The baseline hazard is the estimated hazard rate for a reference subject with all covariates at their mean values. It describes the underlying risk trajectory estimated by the model and is useful for predictions and visualization. While the point estimate is always available, confidence intervals quantify the uncertainty in the baseline hazard estimate.

You can retrieve the baseline hazard with confidence intervals using the ci=True parameter:

# Load data and build the response
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

# Fit a Cox model with age, sex, and ECOG score as covariates
cox = gw.CoxPH().fit(y, lung[["age", "sex", "ph.ecog"]])

# Retrieve baseline cumulative hazard and survival with 95% CI (log-log transform)
baseline = cox.baseline_hazard(ci=True, conf_type="log-log", format="polars")

# Preview the baseline hazard frame
baseline

PolarsRows185Columns7
	time f64	cumhaz f64	survival f64	cumhaz_lower f64	cumhaz_upper f64	survival_lower f64	survival_upper f64
0	5	0.00276281507377	0.997240997987	0.00027467523281	0.0277897175284	0.972592864531	0.999725362487
1	11	0.011118182539	0.988943396027	0.00234289270617	0.0527612650144	0.948606450999	0.997659849725
2	12	0.0139379475929	0.986158735888	0.00312286426338	0.0622077575966	0.939687639264	0.996882006805
3	13	0.0196518735056	0.980539965838	0.00474849859508	0.0813301561637	0.921889273769	0.9952627577
4	15	0.0225413114407	0.977710845715	0.0055786533136	0.0910812507795	0.912943530512	0.994436878477
5	26	0.0254410963375	0.974879801265	0.0064168407806	0.100867296694	0.904052995747	0.993603703176
6	30	0.0283685830848	0.97203002694	0.00726692692643	0.110745094094	0.89516690318	0.992759413344
7	31	0.0313239923381	0.969161521296	0.00812825275623	0.120713826873	0.886287554986	0.991904692168
182	965	2.07762788042	0.125226913368	0.602533738076	7.16397661526	0.000773970642048	0.547422851307
183	1010	2.07762788042	0.125226913368	0.602533738076	7.16397661526	0.000773970642048	0.547422851307
184	1022	2.07762788042	0.125226913368	0.602533738076	7.16397661526	0.000773970642048	0.547422851307

The returned DataFrame includes:

• time: event times at which the baseline hazard is evaluated
• cumhaz: cumulative baseline hazard H_0(t)
• cumhaz_lower, cumhaz_upper: confidence interval bounds for cumulative hazard
• survival: baseline survival probability S_0(t) = \exp(-H_0(t))
• survival_lower, survival_upper: confidence interval bounds for survival

The conf_type= parameter controls the confidence interval transform:

• "log-log" (default): uses a log-log transformation that ensures bounds respect the survival probability constraint (staying between 0 and 1). Recommended for most applications.
• "plain": Wald confidence intervals without transformation. Simpler but may produce invalid bounds (negative cumulative hazard or survival > 1).

The baseline hazard can be combined with individual predictions to compute personalized survival curves with uncertainty bands. For stratified models, each stratum has its own baseline hazard and corresponding confidence intervals.

# For stratified models, each stratum gets its own baseline hazard and CI
cox_strat = gw.CoxPH().fit(y, lung[["age", "ph.ecog"]], strata=lung["sex"])
baseline_strat = cox_strat.baseline_hazard(ci=True, format="polars")

# Preview the per-stratum baseline frame
baseline_strat

PolarsRows205Columns8
	time f64	cumhaz f64	survival f64	cumhaz_lower f64	cumhaz_upper f64	survival_lower f64	survival_upper f64	strata i64
0	11	0.00680528755742	0.993217815974	0.00201397599597	0.0229952784104	0.977267098017	0.997988050693	1
1	12	0.00911229164612	0.990929099465	0.00306148138876	0.0271221178573	0.973242383993	0.996943200167	1
2	13	0.0138099489996	0.986284970897	0.00529883925165	0.0359917865617	0.964648216529	0.994715174833	1
3	15	0.0161956325551	0.983934811544	0.00646151326631	0.0405939759076	0.960218922826	0.99355931742	1
4	26	0.0185931740912	0.981578612636	0.0076409686919	0.0452437559588	0.955764480149	0.992388149299	1
5	30	0.0210224764464	0.979196955454	0.00884699710557	0.0499541833987	0.951273007598	0.99119202242	1
6	31	0.023484045062	0.976789559161	0.0100786043272	0.0547199150369	0.946750281339	0.989972014607	1
7	53	0.0284573689275	0.971943728266	0.012590303582	0.0643210738331	0.937703879014	0.987488622708	1
202	765	0.800707896517	0.449010998265	0.249157461855	2.57320463442	0.0762906695591	0.77945722894	2
203	821	0.800707896517	0.449010998265	0.249157461855	2.96496039698	0.0515625112356	0.805544638503	2
204	965	0.800707896517	0.449010998265	0.249157461855	3.47088969895	0.0310893581836	0.83133932731	2

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.