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
ySurv(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 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.
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.
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.
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-06The 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.
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.
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 as this 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 |
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.
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.
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.
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.
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.
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.
Greenwood’s support for both Efron and Breslow tie-handling methods gives you flexibility that is useful in several situations:
ties="breslow" to get results that match exactly.Most analyses will stick with the default Efron method (which matches R), but this flexibility means you’re never locked into one approach.
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.
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 onwardThe 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
intervalsPandasRows8Columns5 | |||||
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.
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.
When a covariate changes several times for the same subject (as in this section), that subject’s own intervals should be measured from their own zero, not from a shared calendar clock. Mixing absolute dates across subjects is a distinct problem from left truncation (late entry), where a single subject legitimately enters the risk set at a delayed start time, for example when age is used as the time scale and follow-up begins at each subject’s enrollment age rather than at birth. Left truncation is fully supported (see Survival data) and does not require every subject’s first interval to start at 0.
Correct (subject-relative): Each subject’s earliest interval starts at 0
Incorrect (calendar time): Subjects enter the study at different 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"])
dfPandasRows6Columns6 | ||||||
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")
)
dfPolarsRows6Columns6 | ||||||
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 warns when it detects widely-spaced entry times in counting-process data, since that pattern often indicates calendar time was used by mistake. This heuristic can also fire on a genuine left-truncation design where subjects legitimately enter at very different times (for example, staggered enrollment ages). If you have confirmed your entry times are intentional late entries rather than an artifact of calendar time, the warning can be safely ignored.
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.
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.
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.
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_frailtyCoxPH (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.4917To get the frailty-variance inference directly, call 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.theta = 0 (no frailty variance).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 and \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_lnCoxPH (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{'theta': 0.014255476590025323,
'lr_statistic': 0.05125652908554912,
'df': 1.0,
'p_value': 0.4104455340191376}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.
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.
Both frailty distributions currently require:
Surv.right(...) responses,ties="breslow",strata=, andInstead 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.
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)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.
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
lassoCoxNet (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 = 4The 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.
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_resultCoxNetCV (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_lassoCoxNet (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 = 2To 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_lassoCoxNet (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 = 2The 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_brierCoxNetCV (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]Cross-validation is the gold standard for tuning the penalty strength, but when refitting is expensive or you want a quick comparison across several candidate models, AIC and BIC offer a closed-form alternative. Both are based on the partial log-likelihood penalized by model complexity.
For a penalized model, the effective number of parameters is not simply the count of covariates. effective_df() returns the effective degrees of freedom: for pure lasso it is the number of non-zero coefficients, and for ridge or elastic-net it is the trace of the hat matrix, which accounts for the continuous shrinkage applied to each coefficient.
# Fit ridge (l1_ratio=0) and lasso (l1_ratio=1) models at the same penalty
ridge = gw.CoxNet(penalizer=0.01, l1_ratio=0.0).fit(y, lung[cols])
lasso = gw.CoxNet(penalizer=0.01, l1_ratio=1.0).fit(y, lung[cols])
# Ridge shrinks continuously: effective df < number of covariates
ridge.effective_df()4.88503624140628aic() and bic() combine the log-likelihood with the effective degrees of freedom. AIC uses a penalty of 2 per effective parameter. BIC uses log(n_events), which penalizes complexity more heavily in larger samples. Lower values indicate a better trade-off between fit and parsimony.
(1326.816320482443, 1327.4381404213532)(np.float64(1341.5558418357482), np.float64(1342.524539605428))You can also sweep a range of penalty values and pick the one that minimizes AIC or BIC, as a faster alternative to cross-validation when a single train/test split is acceptable:
import numpy as np
# Sweep 20 penalty values and collect AIC for each
penalizers = np.logspace(-3, 0, 20)
aic_values = []
for lam in penalizers:
model = gw.CoxNet(penalizer=lam, l1_ratio=1.0).fit(y, lung[cols])
aic_values.append(model.aic())
# The penalty that minimizes AIC
best_idx = int(np.argmin(aic_values))
print(f"Best penalizer by AIC: {penalizers[best_idx]:.4f} (AIC = {aic_values[best_idx]:.1f})")Best penalizer by AIC: 0.0010 (AIC = 1327.0)AIC tends to select models that predict well on new data. BIC tends to select the true model (if it is among the candidates) and favors simpler models as the sample size grows. Cross-validation makes no distributional assumptions and directly estimates out-of-sample performance. When computational cost is not a concern, cross-validation via cv_coxnet() is generally preferred.
The l1_ratio parameter controls the blend between lasso (1.0) and ridge (0.0). Intermediate values give elastic-net regularization, which combines variable selection with grouped shrinkage. This is useful when covariates are correlated: pure lasso tends to pick one from each correlated group and drop the rest, while elastic-net keeps the group together.
# Elastic-net with equal lasso and ridge contributions
enet = gw.CoxNet(penalizer=0.05, l1_ratio=0.5).fit(y, lung[cols])
enetCoxNet (elastic-net Cox, penalizer=0.05, l1_ratio=0.5)
coef
age 0.00961
sex -0.479
ph.ecog 0.4314
ph.karno 0.0007816
wt.loss -0.004202
n = 213, events = 151, nonzero coefficients = 5Cross-validation works the same way. To search over both penalizer and l1_ratio, run cv_coxnet() at several mixing values and compare:
# Compare CV scores at three mixing levels
for ratio in [0.2, 0.5, 1.0]:
cv = gw.cv_coxnet(y, lung[cols], l1_ratio=ratio, seed=23)
print(f"l1_ratio={ratio:.1f} best_penalizer={cv.best_penalizer_:.4f} "
f"mean_score={cv.best_score_:.4f}")l1_ratio=0.2 best_penalizer=0.2660 mean_score=0.6380
l1_ratio=0.5 best_penalizer=0.1411 mean_score=0.6349
l1_ratio=1.0 best_penalizer=0.1077 mean_score=0.6356Once you have a tuned CoxNet, you can generate predictions for new subjects. The predict() method supports three prediction types, matching the interface of CoxPH:
"lp" (default): the centered linear predictor X\beta, a unitless risk score. Higher values mean higher hazard. A score of zero corresponds to a subject with average covariate values."risk": the relative hazard \exp(X\beta). A value of 2.0 means twice the baseline hazard."survival": subject-specific survival probabilities S(t \mid x) at specified time points, computed from the Breslow baseline hazard and each subject’s relative risk.Start by refitting the final model at the CV-selected penalty:
# Refit the final model from CV
final = gw.CoxNet(penalizer=cv_result.best_penalizer_, l1_ratio=1.0).fit(y, lung[cols])
# Linear predictor (risk score) for the first three subjects
final.predict(lung[cols][:3], type="lp")array([ nan, -0.14066634, -0.14066634])The linear predictor is centered so that the average training subject scores near zero. Positive values indicate above-average risk, negative values below-average. Exponentiating gives the relative risk, which is often more interpretable:
array([ nan, 0.86877914, 0.86877914])For survival predictions, pass type="survival" together with the time points of interest. The result is a DataFrame with one column per subject and one row per time point. Each cell is the estimated probability that the subject survives beyond that time:
# Predicted survival curves at 6-month intervals
surv = final.predict(lung[cols][:3], type="survival", times=[180, 365, 540], format="polars")
survPolarsRows3Columns4 | ||||
time f64 |
subject_1 f64 |
subject_2 f64 |
subject_3 f64 |
|
|---|---|---|---|---|
| 0 | 180 | NaN | 0.779641328137 | 0.779641328137 |
| 1 | 365 | NaN | 0.477563524857 | 0.477563524857 |
| 2 | 540 | NaN | 0.32013996854 | 0.32013996854 |
When newdata is omitted, predictions default to the training data. This is useful for computing training-set risk scores for downstream analyses like concordance or calibration:
# Training-set risk scores (one per subject in the original data)
risk = final.predict(type="risk")
risk[:5]array([0.86877914, 0.86877914, 1.11535798, 0.86877914, 1.11535798])These risk scores rank subjects by predicted hazard. You can pass them directly to concordance_index() to measure how well the penalized model discriminates high-risk from low-risk subjects, or feed them into calibration() to check whether predicted and observed risks agree. See Prediction performance for details on both metrics.
Tables of survival probabilities are useful for exact numbers, but a plot makes it much easier to compare how different subjects’ predicted survival diverges over time. plot_predicted_survival() draws one step-function curve per subject, using the same predict() machinery under the hood.
# Predicted survival curves for four subjects from the penalized model
gw.plot_predicted_survival(final, lung[cols][:4])Each curve traces the estimated probability of surviving past each time point. Subjects whose covariates put them at higher risk drop faster. You can add custom labels with the labels= parameter to make the legend more informative.
Plotting how coefficients evolve across a range of penalty values reveals which covariates enter the model first and which are the most robust to regularization. The strongest predictors have non-zero coefficients even at high penalty values.
import numpy as np
# Fit the lasso across a range of penalties and collect coefficients
penalizers = np.logspace(-4, 0, 40)
paths = {name: [] for name in cols}
for lam in penalizers:
m = gw.CoxNet(penalizer=lam, l1_ratio=1.0).fit(y, lung[cols])
for j, name in enumerate(cols):
paths[name].append(m.coef_[j])The loop fits 40 models across a log-spaced grid of penalties, storing each covariate’s coefficient at every step. Reading from right to left (strong to weak penalty), you can see which covariates enter the model first as regularization relaxes. Plotting the result makes this easier to read:
import altair as alt
import pandas as pd
# Reshape into long form for Altair
rows = []
for j, name in enumerate(cols):
for i, lam in enumerate(penalizers):
rows.append({"penalizer": lam, "coefficient": paths[name][i], "term": name})
path_df = pd.DataFrame(rows)
# Plot the coefficient path
(
alt.Chart(path_df)
.mark_line()
.encode(
x=alt.X("penalizer:Q", scale=alt.Scale(type="log"), title="Penalizer (log scale)"),
y=alt.Y("coefficient:Q", title="Coefficient"),
color=alt.Color("term:N", title="Term"),
)
.properties(width=500, height=280, title="Lasso coefficient path")
)Covariates that reach zero early (moving right to left along the x-axis) contribute less to the model. This path is a useful complement to the AIC/BIC sweep for understanding which variables matter.
The tidy() and glance() functions work with CoxNet, following the same pattern as other Greenwood estimators. This makes it easy to extract results programmatically for reporting or for comparing several models in a table.
tidy() returns one row per covariate with the penalized coefficient estimate and its hazard ratio:
# Tidy: one row per covariate with the penalized estimate and hazard ratio
gw.tidy(final, format="polars")PolarsRows5Columns3 | |||
term str |
estimate f64 |
hazard_ratio f64 |
|
|---|---|---|---|
| 0 | age | 0 | 1 |
| 1 | sex | -0.229727659963 | 0.794750015282 |
| 2 | ph.ecog | 0.249841749326 | 1.28382223488 |
| 3 | ph.karno | -0 | 1 |
| 4 | wt.loss | -0 | 1 |
The estimate column is the log hazard ratio on the original covariate scale. Covariates driven to exactly zero by the lasso still appear in the table, making it clear which variables were excluded. Pass exponentiate=True to express the estimates as hazard ratios directly, which is easier to interpret:
PolarsRows5Columns3 | |||
term str |
estimate f64 |
hazard_ratio f64 |
|
|---|---|---|---|
| 0 | age | 1 | 1 |
| 1 | sex | 0.794750015282 | 0.794750015282 |
| 2 | ph.ecog | 1.28382223488 | 1.28382223488 |
| 3 | ph.karno | 1 | 1 |
| 4 | wt.loss | 1 | 1 |
glance() returns a single-row summary of model-level fit statistics. For CoxNet this includes the log-likelihood, AIC, BIC, effective degrees of freedom, the penalty settings used, and the number of non-zero coefficients:
PolarsRows1Columns9 | |||||||||
n i64 |
nevent i64 |
loglik f64 |
aic f64 |
bic f64 |
effective_df f64 |
penalizer f64 |
l1_ratio f64 |
n_nonzero i64 |
|
|---|---|---|---|---|---|---|---|---|---|
| 0 | 213 | 151 | -665.540909058 | 1335.08181812 | 1341.11637779 | 2 | 0.107671018824 | 1 | 2 |
The effective_df column is particularly useful for comparing models across different penalty values. It reflects the true complexity of the fit: for a pure lasso it equals the number of surviving coefficients, while for ridge or elastic-net it accounts for the continuous shrinkage applied to each term.
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
baselinePolarsRows185Columns7 | |||||||
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 evaluatedcumhaz: cumulative baseline hazard H_0(t)cumhaz_lower, cumhaz_upper: confidence interval bounds for cumulative hazardsurvival: baseline survival probability S_0(t) = \exp(-H_0(t))survival_lower, survival_upper: confidence interval bounds for survivalThe 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_stratPolarsRows205Columns8 | ||||||||
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 |
The Cox model assumes that covariates multiply the hazard by a constant factor. The Aalen additive model takes a different approach: it assumes that each covariate adds to the hazard, and that the size of its contribution can change freely over time. The model is:
h(t \mid x) = \beta_0(t) + \beta_1(t) x_1 + \cdots + \beta_p(t) x_p
where each \beta_j(t) is an unspecified function of time. Because the effects are nonparametric and time-varying, the Aalen model is a natural alternative when the proportional hazards assumption is violated or when you want to see how a covariate’s influence evolves.
AalenAdditive is fit by ordinary least squares at each event time. It requires a Surv response and a covariate frame, just like CoxPH. An intercept (\beta_0(t)) is added automatically.
# Load data and build the response
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
# Fit the Aalen additive hazards model
aalen = gw.AalenAdditive().fit(y, lung[["age", "sex"]])
# Print the summary (slope test for each covariate)
aalenAalenAdditive (additive hazards, test='aalen')
slope coef se(coef) z p
Intercept 0.00483 0.005448 0.004703 1.158 0.247
age 8.986e-05 0.0001263 6.903e-05 1.829 0.0674
sex -0.003151 -0.003938 0.001222 -3.222 0.00127
n = 228, events = 165, event times used = 164The printed summary reports a test of whether each covariate’s cumulative effect is significantly different from zero over the follow-up period. The slope column estimates the average rate of change in the cumulative coefficient, and the z and p columns give the corresponding test statistic and p-value.
The model estimates cumulative regression coefficients B_j(t) = \int_0^t \beta_j(s)\,ds at each event time. These are the central output of the Aalen model and show how each covariate’s total effect accumulates over time.
# Cumulative coefficients as a tidy DataFrame (one row per event time)
cumcoef = aalen.cumulative_coefficients(format="polars")
cumcoefPolarsRows164Columns4 | ||||
time f64 |
Intercept f64 |
age f64 |
sex f64 |
|
|---|---|---|---|---|
| 0 | 5 | -0.0250882296543 | 0.000213054848595 | 0.0115932120342 |
| 1 | 11 | -0.0486387455663 | 0.000792563552041 | 0.00568368753452 |
| 2 | 11 | -0.0971721291678 | 0.00175263337549 | 0.000652068767858 |
| 3 | 11 | -0.0957397773905 | 0.00195158095893 | -0.00613536146487 |
| 4 | 12 | -0.120493372499 | 0.00255870331352 | -0.0122860281248 |
| 5 | 13 | -0.152892229773 | 0.0032866444423 | -0.0183010725864 |
| 6 | 13 | -0.145178427114 | 0.00339757279142 | -0.0255431581426 |
| 7 | 15 | -0.152378822785 | 0.00374259235333 | -0.0324718718822 |
| 161 | 765 | -0.00445862876907 | 0.0422876038128 | -0.138155221158 |
| 162 | 791 | 0.406407245594 | 0.0400238799871 | -0.270017134005 |
| 163 | 814 | 0.331079734677 | 0.0476658013845 | -0.480715823961 |
A coefficient that increases steeply means the covariate is actively contributing to hazard during that period. A flat segment means the covariate has no effect at that time. A downward bend (which can happen because the model is additive, not multiplicative) would indicate a protective effect in that interval.
You can plot these curves to visualize how each covariate’s influence changes over time:
import altair as alt
import pandas as pd
# Reshape for plotting (exclude the Intercept for clarity)
cumcoef_pd = cumcoef.to_pandas()
long = cumcoef_pd.melt(id_vars="time", value_vars=["age", "sex"],
var_name="term", value_name="cumulative_coefficient")
alt.Chart(long).mark_line().encode(
x=alt.X("time:Q", title="Time"),
y=alt.Y("cumulative_coefficient:Q", title="Cumulative coefficient"),
color=alt.Color("term:N", title="Term"),
).properties(width=500, height=280, title="Aalen cumulative coefficients")A roughly linear curve suggests a constant effect over time (which a Cox model would capture well). Curvature indicates time-varying effects that the proportional hazards assumption would miss.
predict() generates subject-specific survival probabilities or cumulative hazard estimates at specified time points. It works the same way as for CoxPH and CoxNet:
# Predicted survival curves at 6-month intervals for the first three subjects
aalen.predict(lung[["age", "sex"]][:3], type="survival", times=[180, 365, 540], format="polars")PolarsRows3Columns4 | ||||
time f64 |
subject_1 f64 |
subject_2 f64 |
subject_3 f64 |
|
|---|---|---|---|---|
| 0 | 180 | 0.602004628716 | 0.626021061755 | 0.676966492893 |
| 1 | 365 | 0.300138141161 | 0.320895955823 | 0.366817759964 |
| 2 | 540 | 0.153165162703 | 0.173539991766 | 0.222781295121 |
The survival probabilities are computed from the cumulative hazard, which is built by summing each subject’s covariate-weighted coefficient increments. Pass type="cumhaz" for the cumulative hazard directly:
# Cumulative hazard at the same time points
aalen.predict(lung[["age", "sex"]][:3], type="cumhaz", times=[180, 365, 540], format="polars")PolarsRows3Columns4 | ||||
time f64 |
subject_1 f64 |
subject_2 f64 |
subject_3 f64 |
|
|---|---|---|---|---|
| 0 | 180 | 0.507490144806 | 0.468371263472 | 0.390133500805 |
| 1 | 365 | 1.20351243977 | 1.13663833355 | 1.00289012111 |
| 2 | 540 | 1.87623844504 | 1.75134720607 | 1.50156472813 |
The tidy() and glance() functions work with AalenAdditive. tidy() returns the summary test table with one row per term (intercept plus covariates):
PolarsRows3Columns6 | ||||||
term str |
slope f64 |
coef f64 |
se f64 |
z f64 |
p f64 |
|
|---|---|---|---|---|---|---|
| 0 | Intercept | 0.00482986097649 | 0.00544765328404 | 0.00470324919766 | 1.1582744301 | 0.24675206268 |
| 1 | age | 8.98581665956e-05 | 0.000126254209714 | 6.90309885029e-05 | 1.82894975796 | 0.0674071351672 |
| 2 | sex | -0.00315102875767 | -0.00393789312818 | 0.00122229196633 | -3.22172871674 | 0.00127419718027 |
The slope column is the estimated average rate of change in the cumulative coefficient, and coef is the final cumulative value at the last event time. Together they describe both the trend and total magnitude of each covariate’s effect.
glance() returns a single-row summary of the overall fit:
The Aalen additive model is most useful when:
The trade-off is that the Aalen model has higher variance than Cox because it estimates a separate effect at each event time. With small samples or many covariates, estimates can be noisy. The nmin parameter controls the minimum risk-set size at which estimation stops, which can stabilize the estimates.
You can now fit a Cox model, read hazard ratios, and assess overall fit.