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 object that pairs each follow-up time with an event indicator, and it is what every model on this page is fit against.

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.

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.

gw.tidy(aft, format="polars")

	term str	estimate f64	std_error f64	statistic f64	p_value f64	conf_low f64	conf_high f64
PolarsRows3Columns7
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.

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:

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:

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.

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.

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.

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")

	p f64	subject_1 f64	subject_2 f64
PolarsRows3Columns3
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.

aft.predict(newdata, type="survival", times=[180, 365, 730], format="polars")

	time f64	subject_1 f64	subject_2 f64
PolarsRows3Columns3
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.

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")

	term str	estimate f64	std_error f64	statistic f64	p_value f64	conf_low f64	conf_high f64
PolarsRows6Columns7
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.

# 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.

newdata = pd.DataFrame({"age": [50, 70], "sex": [1, 1]})

rp.predict(newdata, type="survival", times=[180, 365, 730], format="polars")

	time f64	subject_1 f64	subject_2 f64
PolarsRows3Columns3
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 is the semiparametric alternative, useful when you prefer not to assume a distribution.
Prediction performance evaluates any fitted model, parametric or not.
Competing risks handles the case of more than one type of event.