Piecewise exponential models

A piecewise exponential model (PEM) splits the time axis into intervals and assumes a constant hazard within each one. This sits between the Cox model, which leaves the baseline hazard entirely unspecified, and a fully parametric AFT, which assumes a specific distributional shape. The PEM estimates the baseline hazard directly as a step function, giving you interpretable rates in each time window while still allowing covariates to shift the hazard proportionally.

Under the hood the PEM is equivalent to a Poisson GLM on an interval-expanded dataset, with log(exposure) as an offset. Greenwood fits it via Newton-Raphson (iteratively reweighted least squares), so no additional dependencies are needed.

We use the lung dataset throughout:

import greenwood as gw

# Load the lung dataset
lung = gw.load_dataset("lung", backend="polars")

Each model names its response and covariates in a formula, Surv(time, status) ~ age + sex, where a status of 2 marks a death.

Fitting with manual break points

Pass breaks= to set the times where the hazard is allowed to change. Three intervals are a common starting point: early, middle, and late follow-up.

# Fit a PEM with breaks at 180 and 365 days
pem = gw.PiecewiseExponential(breaks=[180, 365]).fit("Surv(time, status) ~ age + sex", data=lung)
pem

PiecewiseExponential (3 intervals)

        coef  se(coef)       z         p
age  0.01579  0.009175   1.721   0.08525
sex  -0.5077    0.1672  -3.037  0.002388

Intervals:
  (0, 180]
  (180, 365]
  (365, inf]

n = 228, events = 165
Log-likelihood = -426.6
AIC = 863.2

The coefficient table looks like a Cox model. The age and sex effects are log hazard ratios, interpreted the same way. The difference is that the baseline hazard is now explicitly estimated rather than left nonparametric.

# Coefficient table with log hazard ratios
pem.to_frame(format="polars")

	term str	estimate f64	std_error f64	statistic f64	p_value f64	conf_low f64	conf_high f64
PolarsRows2Columns7
0	age	0.0157901291643	0.00917500504479	1.72099405802	0.0852519021378	-0.00219255028142	0.0337728086101
1	sex	-0.507742560024	0.167172032413	-3.03724584009	0.00238750650133	-0.835393722775	-0.180091397273

Inspecting the baseline hazard

The baseline_hazard() method returns the estimated constant hazard rate in each interval.

# Constant hazard rate within each interval
pem.baseline_hazard(format="polars")

	start f64	stop f64	hazard f64	log_hazard f64
PolarsRows3Columns4
0	0	180	0.00129910865297	-6.64607690124
1	180	365	0.00221541176629	-6.11231699376
2	365	Inf	0.00239171750055	-6.03574355165

Higher hazard values mean more events per unit time in that interval, adjusting for covariates.

Automatic knot selection

When you do not supply breaks=, Greenwood tries 0 through max_knots internal knots and picks the set that minimizes the AIC (or BIC, with knot_strategy="bic"). Knot candidates are placed at quantiles of the event-time distribution.

# Let AIC choose the best number of knots automatically
pem_auto = gw.PiecewiseExponential(knot_strategy="aic", max_knots=8).fit(
    "Surv(time, status) ~ age + sex", data=lung
)
pem_auto

PiecewiseExponential (1 intervals)

        coef  se(coef)       z         p
age  0.01562  0.009106   1.715    0.0863
sex  -0.4809    0.1671  -2.878  0.003999

Intervals:
  (0, inf]

n = 228, events = 165
Log-likelihood = -281.9
AIC = 569.7

If automatic selection returns zero breaks, the data do not support a more complex baseline than a single exponential. This is a useful model-selection result in itself.

Predicting survival curves

The predict() method returns survival probabilities (or cumulative hazards) at specified times.

# Predicted survival probabilities for the first three subjects
pred = pem.predict(
    newdata=lung.head(3),
    type="survival",
    times=[100, 200, 365, 500, 730],
    format="polars",
)
pred

	time f64	subject_1 f64	subject_2 f64	subject_3 f64
PolarsRows5Columns4
0	100	0.777607528515	0.795490067227	0.827535046911
1	200	0.583594083846	0.612706408535	0.666768710419
2	365	0.287562906672	0.321853848354	0.391419745741
3	500	0.153896140388	0.182261639537	0.24451649032
4	730	0.0530473167585	0.069173978742	0.109693925366

The piecewise constant baseline means survival curves are piecewise exponential: straight segments on a log scale, with slope changes at the break points.

Cumulative hazard is available with type="cumhaz":

# Predicted cumulative hazard for the same subjects
cumhaz = pem.predict(
    newdata=lung.head(3),
    type="cumhaz",
    times=[100, 200, 365, 500, 730],
    format="polars",
)
cumhaz

	time f64	subject_1 f64	subject_2 f64	subject_3 f64
PolarsRows5Columns4
0	100	0.25153334415	0.228796917484	0.189303819865
1	200	0.538549599729	0.489869399806	0.405312054192
2	365	1.24631363686	1.13365772358	0.937974776287
3	500	1.87147731715	1.70231204436	1.40847252727
4	730	2.93657099468	2.67113051682	2.2100612882

Comparing with Cox

PEM covariate estimates should be close to Cox when the interval structure is flexible enough. Large differences may signal that the proportional hazards assumption holds differently once the baseline is parameterized.

# Cox model coefficients for reference
cox = gw.CoxPH().fit("Surv(time, status) ~ age + sex", data=lung)
cox.to_frame(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.00922327347697	1.8480783301	0.0645910121296	-0.00103195198901	0.0351226156798
1	sex	-0.513218517108	0.167457962356	-3.06476031291	0.00217844504671	-0.84143009225	-0.185006941967

# PEM coefficients for side-by-side comparison
pem.to_frame(format="polars")

	term str	estimate f64	std_error f64	statistic f64	p_value f64	conf_low f64	conf_high f64
PolarsRows2Columns7
0	age	0.0157901291643	0.00917500504479	1.72099405802	0.0852519021378	-0.00219255028142	0.0337728086101
1	sex	-0.507742560024	0.167172032413	-3.03724584009	0.00238750650133	-0.835393722775	-0.180091397273

Tidy summaries

The PEM integrates with the tidy/glance framework:

# Tidy coefficient table
gw.tidy(pem, format="polars")

	term str	estimate f64	std_error f64	statistic f64	p_value f64	conf_low f64	conf_high f64
PolarsRows2Columns7
0	age	0.0157901291643	0.00917500504479	1.72099405802	0.0852519021378	-0.00219255028142	0.0337728086101
1	sex	-0.507742560024	0.167172032413	-3.03724584009	0.00238750650133	-0.835393722775	-0.180091397273

# One-row model summary with AIC, BIC, and fit statistics
gw.glance(pem, format="polars")

	n_intervals i64	n i64	nevent i64	loglik f64	aic f64	bic f64	df i64	lr_statistic f64
PolarsRows1Columns8
0	3	228	165	-426.579975304	863.159950608	880.306678752	5	13.4598097006

When to use a PEM

The piecewise exponential model is a good choice when:

You want an explicit baseline hazard estimate (for reporting rates, not just relative effects).
The hazard has clear period structure (early mortality, treatment plateau, late relapse).
You need to extrapolate beyond observed follow-up but do not want to commit to a parametric distribution shape.
You are bridging between a nonparametric analysis and a fully parametric one and want to check sensitivity to baseline hazard assumptions.