# Parametric survival models
The Cox model deliberately avoids specifying the shape of the baseline hazard. Parametric models take the opposite approach: they assume the survival times follow a particular distribution, such as Weibull or log-normal. In exchange for that assumption you gain a fully specified model that can extrapolate beyond the observed follow-up, produce smooth survival and hazard curves, and sometimes fit more efficiently. Greenwood provides these as accelerated failure time models, which describe how covariates stretch or compress the time scale. This page shows how to fit them and how to read their coefficients.
We work from the `lung` outcomes. The response `y` is a [Surv](../reference/Surv.md#greenwood.Surv) object that pairs each follow-up time with an event indicator, and it is what every model on this page is fit against.
``` python
import greenwood as gw
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
y
```
Surv(type=right, n=228, events=165)
The printed response shows 165 events among 228 subjects, with a `+` marking each censored time.
# The accelerated failure time model
An accelerated failure time model, or AFT model, works on the logarithm of survival time. It says that covariates act by multiplying the time scale: a covariate might make time pass twice as fast, halving survival, or twice as slowly, doubling it. This is a different and often more intuitive framing than the Cox model's multiplication of hazards. The model is fit by maximum likelihood, and it includes an intercept because it estimates the actual location of the survival distribution, not just relative effects.
You fit one by naming a distribution. The Weibull is the default and the most common choice.
``` python
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
aft = gw.AFT("weibull").fit(y, lung[["age", "sex"]])
aft
```
AFT (accelerated failure time model, dist='weibull')
coef se(coef) z p
(Intercept) 6.275 0.4814 13.036 7.65e-39
age -0.01226 0.006957 -1.762 0.0781
sex 0.3821 0.1275 2.997 0.002723
Scale = 0.7541
n = 228, events = 165
Log-likelihood = -1147
Printing the fitted model gives a summary in the style of R's `survreg`: the coefficient table, the scale parameter, the sample size, and the log-likelihood. For the coefficients as data, pass the model to `gw.tidy`, which returns a tidy frame.
``` python
gw.tidy(aft, format="polars")
```
PolarsRows3Columns7
|
|
term
str
|
estimate
f64
|
std_error
f64
|
statistic
f64
|
p_value
f64
|
conf_low
f64
|
conf_high
f64
|
| 0 |
(Intercept) |
6.27487873815 |
0.481355210911 |
13.0358591658 |
7.65003888482e-39 |
5.331439861 |
7.21831761531 |
| 1 |
age |
-0.0122574223003 |
0.00695735246721 |
-1.76179406722 |
0.0781040974976 |
-0.0258935825638 |
0.00137873796317 |
| 2 |
sex |
0.3820847039 |
0.127473144911 |
2.99737410706 |
0.00272316307345 |
0.132241930879 |
0.631927476921 |
The coefficients are on the log-time scale. A positive coefficient lengthens survival time, and a negative coefficient shortens it. This is the opposite direction from a Cox hazard ratio, where a positive coefficient means higher risk and shorter survival, so take care when comparing the two.
# Choosing a distribution
Greenwood supports four parametric distributions, each giving a different hazard shape.
- exponential: Constant hazard (h(t) = λ). Simplest; often too restrictive.
- Weibull: monotonic hazard (increasing or decreasing). Default; most flexible while remaining parsimonious.
- log-normal: hazard rises then falls; many peak early. Good for disease incidence.
- log-logistic: similar to log-normal but with heavier tails.
``` python
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
for dist in ("exponential", "weibull", "lognormal", "loglogistic"):
model = gw.AFT(dist).fit(y, lung[["age", "sex"]])
print(f"{dist:20s} loglik={model.loglik_:10.3f} scale={model.scale_:.4f}")
```
exponential loglik= -1156.099 scale=1.0000
weibull loglik= -1147.054 scale=0.7541
lognormal loglik= -1158.750 scale=1.0527
loglogistic loglik= -1152.897 scale=0.5656
To compare and choose among distributions, use the AIC (available from `glance`), where a lower value indicates a better fit:
``` python
import pandas as pd
results = []
for dist in ("exponential", "weibull", "lognormal", "loglogistic"):
model = gw.AFT(dist).fit(y, lung[["age", "sex"]])
glance_result = gw.glance(model, format="pandas")
glance_result["distribution"] = dist
results.append(glance_result)
comparison = pd.concat(results)[["distribution", "loglik", "aic"]]
print(comparison.sort_values("aic"))
```
distribution loglik aic
0 weibull -1147.054431 2302.108863
0 loglogistic -1152.897225 2313.794451
0 exponential -1156.099037 2318.198074
0 lognormal -1158.750143 2325.500285
The scale parameter describes the spread of the distribution. For the exponential, it is fixed at 1 (constant hazard). For Weibull, values \> 1 indicate increasing hazard, \< 1 indicate decreasing. Other distributions use scale differently; see the model's docstring for details.
# Model comparison and selection
Each call to `gw.glance` returns a one-row DataFrame of model-level summaries for a single fit. We can compare all four distributions systematically:
``` python
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
comparisons = []
for dist in ("exponential", "weibull", "lognormal", "loglogistic"):
model = gw.AFT(dist).fit(y, lung[["age", "sex"]])
result = gw.glance(model, format="pandas")
result["distribution"] = dist
comparisons.append(result)
comparison = pd.concat(comparisons)[["distribution", "loglik", "aic"]]
print(comparison.sort_values("aic"))
```
distribution loglik aic
0 weibull -1147.054431 2302.108863
0 loglogistic -1152.897225 2313.794451
0 exponential -1156.099037 2318.198074
0 lognormal -1158.750143 2325.500285
Compare the models by reading the AIC column and preferring the lowest value. A difference of 10 in AIC is substantial; differences of 2-3 are subtle and don't warrant changing models.
**When to use each distribution:**
- exponential: only if you're confident hazard is constant (rare).
- Weibull: safe default that fits many data types well.
- log-normal/log-logistic: When you expect a hazard peak early in follow-up.
- Gompertz: for aging-related or mortality data where hazard accelerates exponentially.
> **Note: Parametric or Cox?**
>
> Choose a parametric model when you need to extrapolate, want smooth hazard or quantile predictions, or have a distributional form suggested by theory. Choose the Cox model when you want to avoid distributional assumptions and care mainly about relative effects. Both are valid; they answer slightly different questions.
> **Tip: Interpreting the acceleration factor**
>
> Exponentiating an AFT coefficient gives a time-acceleration factor. A factor of 1.5 for a covariate means subjects with a one-unit-higher value survive 1.5 times as long on average, holding other covariates fixed.
# Predicting survival and quantiles
Because a parametric model specifies the whole distribution, it can predict smoothly at any time and extrapolate beyond the observed follow-up. We predict for two new subjects, a 50-year-old and a 70-year-old, both with `sex` coded 1.
``` python
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
aft = gw.AFT("weibull").fit(y, lung[["age", "sex"]])
newdata = pd.DataFrame({"age": [50, 70], "sex": [1, 1]})
aft.predict(newdata, type="quantile", p=[0.25, 0.5, 0.75], format="polars")
```
PolarsRows3Columns3
|
|
p
f64
|
subject_1
f64
|
subject_2
f64
|
| 0 |
0.25 |
164.781523474 |
128.956101802 |
| 1 |
0.5 |
319.808645462 |
250.278522567 |
| 2 |
0.75 |
539.364324762 |
422.100241011 |
With `type="quantile"` the model returns survival-time quantiles: each column is a subject and each row a failure probability. The middle row (`p = 0.5`) is the predicted median survival time, and the older subject's is shorter, as expected. These quantiles match R's `survreg`.
For the survival curve itself, use `type="survival"` with the times you want.
``` python
aft.predict(newdata, type="survival", times=[180, 365, 730], format="polars")
```
PolarsRows3Columns3
|
|
time
f64
|
subject_1
f64
|
subject_2
f64
|
| 0 |
180 |
0.723657400537 |
0.639099212255 |
| 1 |
365 |
0.437821037604 |
0.31877897223 |
| 2 |
730 |
0.126066878662 |
0.0568943508148 |
Each column is a subject and each row a requested time, giving the estimated probability of surviving past that time.
# Flexible parametric (Royston-Parmar) models
The AFT distributions impose a fixed hazard shape. When none of them fits well but you still want the smooth, extrapolatable curves of a parametric model, a Royston-Parmar model is a good middle ground. It models the log cumulative hazard as a restricted cubic spline in log time, so the baseline shape is estimated from the data rather than assumed. The flexibility is set by `df`, the number of spline degrees of freedom.
``` python
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
rp = gw.RoystonParmar(df=3).fit(y, lung[["age", "sex"]])
rp.to_frame(format="polars")
```
PolarsRows6Columns7
|
|
term
str
|
estimate
f64
|
std_error
f64
|
statistic
f64
|
p_value
f64
|
conf_low
f64
|
conf_high
f64
|
| 0 |
gamma0 |
-7.22845418511 |
1.32694246946 |
-5.44745107755 |
5.10967630237e-08 |
-9.8292136348 |
-4.62769473541 |
| 1 |
gamma1 |
1.02751158554 |
0.297232335012 |
3.45693070539 |
0.000546365323608 |
0.444946913874 |
1.6100762572 |
| 2 |
gamma2 |
-0.0964201003026 |
0.130227015742 |
-0.740400137045 |
0.459057235695 |
-0.351660360972 |
0.158820160366 |
| 3 |
gamma3 |
0.117232405393 |
0.185653198447 |
0.631459120415 |
0.527740370406 |
-0.246641177177 |
0.481105987964 |
| 4 |
age |
0.0161468803859 |
0.00919444985347 |
1.75615514176 |
0.0790619218346 |
-0.00187391018459 |
0.0341676709563 |
| 5 |
sex |
-0.510127303404 |
0.167176594587 |
-3.05142777112 |
0.00227755838203 |
-0.837787407852 |
-0.182467198956 |
The `gamma` terms are the spline coefficients for the baseline log cumulative hazard, and the named terms are the covariate effects on that scale.
## Choosing the degrees of freedom
The `df` parameter controls spline flexibility. `df=1` is exactly a Weibull model (no spline); higher values let the hazard adapt to the data shape. More flexibility fits the observed data better but risks overfitting and unstable extrapolation.
``` python
# Compare different df values using loglik
results = []
for df in (1, 2, 3, 4, 5):
rp_df = gw.RoystonParmar(df=df).fit(y, lung[["age", "sex"]])
results.append({"df": df, "loglik": rp_df.loglik_})
comparison = pd.DataFrame(results)
print(comparison)
```
df loglik
0 1 -1147.054431
1 2 -1146.810064
2 3 -1146.570537
3 4 -1146.438891
4 5 -1146.513762
Notice that loglik generally improves up to `df=3`, after which gains diminish, suggesting that additional flexibility is not justified by the data. Two or three degrees of freedom are typical defaults. Choose based on the trade-off between fit and complexity:
- **df=1**: equivalent to Weibull; simplest.
- **df=2-3**: common choices; balance parsimony and fit.
- **df=4+**: use only if AIC clearly favors it and validation shows good extrapolation.
Prediction works as for the AFT model, with survival, hazard, or cumulative hazard curves.
``` python
newdata = pd.DataFrame({"age": [50, 70], "sex": [1, 1]})
rp.predict(newdata, type="survival", times=[180, 365, 730], format="polars")
```
PolarsRows3Columns3
|
|
time
f64
|
subject_1
f64
|
subject_2
f64
|
| 0 |
180 |
0.732463257148 |
0.650496589053 |
| 1 |
365 |
0.43369897704 |
0.315422138026 |
| 2 |
730 |
0.123266741824 |
0.0555002980129 |
# Next steps
You can now fit, compare, and predict from parametric survival models.
- [Cox regression](cox-regression.md) is the semiparametric alternative, useful when you prefer not to assume a distribution.
- [Prediction performance](prediction-performance.md) evaluates any fitted model, parametric or not.
- [Competing risks](competing-risks.md) handles the case of more than one type of event.