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

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

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.

pem = gw.PiecewiseExponential(breaks=[180, 365]).fit(y, lung[["age", "sex"]])
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.

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.

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.

pem_auto = gw.PiecewiseExponential(knot_strategy="aic", max_knots=8).fit(
    y, lung[["age", "sex"]]
)
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.

pred = pem.predict(
    newdata=lung[["age", "sex"]].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":

cumhaz = pem.predict(
    newdata=lung[["age", "sex"]].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 = gw.CoxPH().fit(y, lung[["age", "sex"]])
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.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:

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

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.