# 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:
``` python
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.
``` python
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.
``` python
pem.to_frame(format="polars")
```
PolarsRows2Columns7
|
|
term
str
|
estimate
f64
|
std_error
f64
|
statistic
f64
|
p_value
f64
|
conf_low
f64
|
conf_high
f64
|
| 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.
``` python
pem.baseline_hazard(format="polars")
```
PolarsRows3Columns4
|
|
start
f64
|
stop
f64
|
hazard
f64
|
log_hazard
f64
|
| 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.
``` python
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.
``` python
pred = pem.predict(
newdata=lung[["age", "sex"]].head(3),
type="survival",
times=[100, 200, 365, 500, 730],
format="polars",
)
pred
```
PolarsRows5Columns4
|
|
time
f64
|
subject_1
f64
|
subject_2
f64
|
subject_3
f64
|
| 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"`:
``` python
cumhaz = pem.predict(
newdata=lung[["age", "sex"]].head(3),
type="cumhaz",
times=[100, 200, 365, 500, 730],
format="polars",
)
cumhaz
```
PolarsRows5Columns4
|
|
time
f64
|
subject_1
f64
|
subject_2
f64
|
subject_3
f64
|
| 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.
``` python
cox = gw.CoxPH().fit(y, lung[["age", "sex"]])
cox.to_frame(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.00922327347697 |
1.8480783301 |
0.0645910121296 |
-0.00103195198901 |
0.0351226156798 |
| 1 |
sex |
-0.513218517108 |
0.167457962356 |
-3.06476031291 |
0.00217844504671 |
-0.84143009225 |
-0.185006941967 |
``` python
pem.to_frame(format="polars")
```
PolarsRows2Columns7
|
|
term
str
|
estimate
f64
|
std_error
f64
|
statistic
f64
|
p_value
f64
|
conf_low
f64
|
conf_high
f64
|
| 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:
``` python
gw.tidy(pem, format="polars")
```
PolarsRows2Columns7
|
|
term
str
|
estimate
f64
|
std_error
f64
|
statistic
f64
|
p_value
f64
|
conf_low
f64
|
conf_high
f64
|
| 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 |
``` python
gw.glance(pem, format="polars")
```
PolarsRows1Columns8
|
|
n_intervals
i64
|
n
i64
|
nevent
i64
|
loglik
f64
|
aic
f64
|
bic
f64
|
df
i64
|
lr_statistic
f64
|
| 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.