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 the Greenwood library
import greenwood as gw

# Load the bundled lung dataset as a Polars DataFrame
lung = gw.load_dataset("lung", backend="polars")

# Build a right-censored response (status 2 = dead in this R-originating dataset)
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

# Display the response summary
y

Surv(type=right, n=228, events=165)

The printed response shows a right-censored response with 165 events among 228 subjects.

Fitting a distribution to survival data

Before adding covariates, it is often useful to ask: what parametric distribution fits this data best? The Parametric class fits a standalone survival distribution by maximum likelihood, giving you parameter estimates, standard errors, and model-selection statistics.

# Fit a Weibull distribution by maximum likelihood (no covariates)
fit = gw.Parametric("weibull").fit(y)

# Display parameter estimates, standard errors, and model-selection statistics
fit

Parametric (weibull distribution)

       estimate  std_error
shape     1.317    0.08221
scale     417.8       24.7

n = 228, events = 165
Log-likelihood = -1154
AIC = 2312, BIC = 2319

The summary reports the natural parameters of the chosen family (here, Weibull shape and scale), the log-likelihood, and AIC/BIC. The estimated shape parameter tells you about the hazard trend: shape > 1 means the hazard increases over time, shape < 1 means it decreases, and shape = 1 is a constant hazard (exponential).

Comparing distributions

Rather than fitting each family one at a time, compare_distributions fits all four and returns a model-selection table sorted by AIC (where lower is better):

# Fit all four families at once and rank by AIC (lower is better)
gw.compare_distributions(y, format="polars")

	dist str	n_params i64	loglik f64	aic f64	bic f64
PolarsRows4Columns5
0	weibull	2	-1153.85118809	2311.70237618	2318.56106744
1	loglogistic	2	-1160.93062351	2325.86124702	2332.71993828
2	exponential	1	-1162.33817579	2326.67635157	2330.1056972
3	lognormal	2	-1169.26905531	2342.53811061	2349.39680187

This immediately tells you which distributional assumption is most supported by the data. A difference of 2 or less in AIC is negligible whereas differences above 10 are decisive.

Predictions from a fitted distribution

A fitted Parametric object gives you the full distributional toolkit without needing to specify covariates:

# Refit the Weibull for prediction
fit = gw.Parametric("weibull").fit(y)

# Query the survival function at 6 months, 1 year, and 2 years
fit.survival([180, 365, 730])

array([0.71893491, 0.43295354, 0.12425187])

# Instantaneous hazard rate at the same time points (risk of event per unit time)
fit.hazard([180, 365, 730])

array([0.00241409, 0.00302016, 0.00376191])

# Time by which 25%, 50%, and 75% of subjects have experienced the event
fit.quantile([0.25, 0.5, 0.75])

array([162.19120953, 316.26369238, 535.36451702])

# Single-number summaries of the fitted distribution
print(f"Mean:   {fit.mean():.1f} days")
print(f"Median: {fit.median():.1f} days")

Mean:   384.9 days
Median: 316.3 days

These predictions describe the population-level survival distribution. Once you are satisfied with a distributional family, you can move to an AFT model to understand how covariates shift that distribution.

From exploration to regression

The Parametric class is for exploration: which distribution shape fits my data? Once you’ve chosen, use AFT with the same distribution to model covariate effects. For example, if compare_distributions() picks Weibull, follow up with gw.AFT("weibull").fit(y, covariates).

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.

# Fit a Weibull AFT model with age and sex as covariates
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

# The model estimates how each covariate stretches or compresses survival time
aft = gw.AFT("weibull").fit(y, lung[["age", "sex"]])

# Display the coefficient table, scale, and log-likelihood
aft

AFT (accelerated failure time model, dist='weibull')

                 coef  se(coef)       z          p
(Intercept)     6.275    0.4814  13.036  7.671e-39
age          -0.01226  0.006957  -1.762    0.07811
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 DataFrame.

# Extract the coefficient table as a tidy DataFrame (one row per term)
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.27487656662	0.481362779739	13.035649682	7.67107986691e-39	5.33142285484	7.21833027841
1	age	-0.0122574075016	0.00695741143991	-1.76177700679	0.0781069810601	-0.0258936833495	0.00137886834625
2	sex	0.382085592498	0.127473120177	2.99738165949	0.00272309560252	0.132242867954	0.631928317043

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 and often too restrictive.
Weibull: monotonic hazard (increasing or decreasing). The default and 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.

# Fit each distribution to the same data and print the log-likelihood and scale
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

# Collect glance summaries from each distribution into one table
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)

# Combine and sort by AIC (lower is better)
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 so see the model’s help page 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:

# Fit each distribution and collect glance results for ranking
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)

# Rank distributions by AIC: the lowest value wins
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 and 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 as 1.

# Fit a Weibull AFT and define two new subjects for prediction
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"]])

# Two new subjects: a 50-year-old and a 70-year-old, both sex = 1
newdata = pd.DataFrame({"age": [50, 70], "sex": [1, 1]})

# Predict survival-time quantiles at 25%, 50%, and 75% failure probabilities
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.781414313	128.956054541
1	0.5	319.808460567	250.278451947
2	0.75	539.364048781	422.100149962

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.

# Survival curve at specific times: probability of being alive past each time point
aft.predict(newdata, type="survival", times=[180, 365, 730], format="polars")

	time f64	subject_1 f64	subject_2 f64
PolarsRows3Columns3
0	180	0.723657198394	0.639099089283
1	365	0.437820768409	0.318778859044
2	730	0.126066714828	0.0568943192275

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.

# Fit a Royston-Parmar flexible-parametric model (3 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"]])

# Display the spline and covariate coefficients
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.22846438615	1.32609527414	-5.45093895371	5.0104568809e-08	-9.82756336353	-4.62936540876
1	gamma1	1.02751318673	0.296797465595	3.46200121578	0.000536174733159	0.445800843461	1.60922553
2	gamma2	-0.0964198530341	0.129713726754	-0.743328061314	0.457283086384	-0.350654085773	0.157814379705
3	gamma3	0.117232173258	0.18493276195	0.633917819763	0.526134455489	-0.245229379726	0.479693726242
4	age	0.0161469261877	0.0091925791922	1.75651749635	0.0790000866033	-0.00187019795402	0.0341640503295
5	sex	-0.510126834073	0.167162261846	-3.05168659744	0.00227559559524	-0.837758846865	-0.182494821281

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.

# Sweep across spline degrees of freedom to find the best-fitting complexity
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_})

# Display the improvement (or lack thereof) as df increases
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 (balances 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.

# Two new subjects for prediction from the Royston-Parmar model
newdata = pd.DataFrame({"age": [50, 70], "sex": [1, 1]})

# Predicted survival at 6 months, 1 year, and 2 years
rp.predict(newdata, type="survival", times=[180, 365, 730], format="polars")

	time f64	subject_1 f64	subject_2 f64
PolarsRows3Columns3
0	180	0.732463442055	0.650496559624
1	365	0.433699254697	0.315422083546
2	730	0.123266928691	0.0555002672218

Survival predictions with confidence intervals

When predicting survival probabilities from an AFT model, you can also compute confidence intervals that quantify uncertainty in the predictions. This is especially valuable when making personalized survival predictions for clinical use.

Confidence intervals are computed via the delta-method, which propagates coefficient uncertainty through the survival function. Wider intervals indicate greater model uncertainty, often due to covariate values that are far from the training data mean.

Pass ci=True to get confidence bounds alongside the survival probability:

# Refit the Weibull AFT for confidence interval demonstration
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"]])

# Two subjects with different covariate profiles
newdata = pd.DataFrame({"age": [50, 70], "sex": [1, 2]})

# Predict survival with delta-method confidence intervals (log-log transform)
predictions = aft.predict(
    newdata, type="survival", times=[180, 365, 730],
    ci=True, conf_type="log-log", format="polars"
)

# Display the point estimates and their confidence bounds
predictions

	time f64	subject_1 f64	subject_1_lower f64	subject_1_upper f64	subject_2 f64	subject_2_lower f64	subject_2_upper f64
PolarsRows3Columns7
0	180	0.723657198394	0.644182108512	0.788300132548	0.763589325374	0.694594995229	0.81903096516
1	365	0.437820768409	0.325293389151	0.544736721602	0.502186717452	0.394310422343	0.600619587507
2	730	0.126066714828	0.0598546494462	0.21803422453	0.177812008035	0.0969662403084	0.278528071515

The returned frame includes one column per subject with point estimates and confidence bounds:

subject_1: survival probability for subject 1
subject_1_lower, subject_1_upper: lower and upper confidence interval bounds for subject 1
(same pattern for subject 2, etc.)

The conf_type= parameter controls the confidence interval transform, just as with baseline hazard confidence intervals:

"log-log" (default): Uses a log-log transformation that ensures bounds respect the survival probability constraint (between 0 and 1). Recommended.
"plain": Wald confidence intervals without transformation. Simpler but may produce invalid bounds (survival < 0 or > 1).

Confidence intervals are valuable for:

clinical decision-making: quantifying uncertainty in individual prognoses
study planning: assessing how precisely your model can predict for new subjects
model diagnostics: detecting when predictions are made far from the training data (wider intervals indicate extrapolation)

The width of the confidence interval grows as covariate values move further from the training mean, reflecting increased model uncertainty in extrapolated regions.

Mean survival and expected remaining lifetime

Survival curves and quantiles describe what happens at specific time points, but sometimes you want a single summary number for each subject: how long are they expected to survive overall? How much longer should a censored subject expect to live given they have already reached their last follow-up time? These questions are answered by mean survival predictions and conditional expectations.

Unconditional mean survival

Using type="mean" returns the expected survival time E[T] for each subject under the fitted distribution. Each distribution has a closed-form formula, so the calculation is exact and fast:

Weibull/Exponential: E[T] = e^{\mu}\,\Gamma(1+\sigma)
Log-normal: E[T] = e^{\mu + \sigma^2/2}
Log-logistic (\sigma < 1): E[T] = e^{\mu}\,\pi\sigma / \sin(\pi\sigma)

where \mu = X\beta is the log-time location for each subject and \sigma is the scale parameter.

# Fit a log-normal AFT model (different distributional choice for illustration)
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
aft = gw.AFT("lognormal").fit(y, lung[["age", "sex"]])

# Two subjects for prediction
newdata = pd.DataFrame({"age": [50, 70], "sex": [1, 1]})

# Unconditional expected survival time E[T] for each subject
aft.predict(newdata, type="mean")

array([551.96643754, 345.96965546])

The returned array has one value per subject. The younger subject has a substantially higher expected survival time, as expected.

Log-logistic mean only exists when \sigma < 1

For log-logistic, the mean is only finite when \sigma < 1. When \sigma \ge 1 the distribution has a heavy enough right tail that the mean diverges. In that case type="mean" returns np.inf. Use type="rmst" with a finite restriction time as an alternative summary.

Conditional mean survival given landmark time

In follow-up studies, subjects are often censored at different times. For a censored subject who has already survived to time t_0, the relevant question is not “how long will they survive in total?” but “how much longer can they expect to live?” This is the conditional mean E[T \mid T > t_0], computed by passing conditional_after:

# Conditional mean given that each subject has already survived one year
aft.predict(newdata, type="mean", conditional_after=365)

array([1013.98560697,  837.9523237 ])

Compare with the unconditional mean to see how much being alive at one year updates the prognosis. The conditioning time can also be a per-subject array. This is useful when each subject has a different follow-up time:

import numpy as np

# Each subject has a different last-known follow-up time (their landmark)
landmark_times = np.array([200.0, 400.0])

# Conditional mean given each subject's own landmark
aft.predict(newdata, type="mean", conditional_after=landmark_times)

array([768.57588911, 889.45750323])

Expected remaining lifetime

type="mean_remaining" computes the expected remaining lifetime E[T - t_0 \mid T > t_0] directly. This is how much longer a subject is expected to survive past their landmark time:

# Expected additional lifetime beyond the landmark: E[T - t0 | T > t0]
remaining = aft.predict(newdata, type="mean_remaining", conditional_after=landmark_times)
remaining

array([568.57588911, 489.45750323])

The relationship between the two quantities is exact:

E[T \mid T > t_0] = t_0 + E[T - t_0 \mid T > t_0]

so you can recover either from the other. The remaining-lifetime form is often more interpretable in a clinical setting: “given that you are alive today, you can expect approximately X more months.”

The memoryless exponential

Under an exponential (constant-hazard) model, E[T - t_0 \mid T > t_0] = E[T] for all t_0. The expected remaining life never changes, regardless of how long the subject has already survived (the defining property of the exponential distribution).

Restricted mean survival time

The restricted mean survival time (RMST) at horizon \tau is the area under the survival curve from 0 to \tau:

\text{RMST}(\tau) = E[\min(T, \tau)] = \int_0^\tau S(t)\,dt

It answers the question “on average, how many days of event-free follow-up does a subject accumulate within the first \tau days?” Unlike the mean, it is always finite regardless of the distribution’s tail behaviour, and it avoids the need to extrapolate beyond the observation window.

Use type="rmst" with the tau argument:

# RMST within one year: expected event-free days accumulated in [0, 365]
aft.predict(newdata, type="rmst", tau=365)

array([261.92109853, 212.64413244])

The result is one RMST value per subject, reflecting how covariate differences shift the entire survival curve within the follow-up window.

Because RMST is bounded by \tau, it is straightforward to compare predictions across subjects or models on an absolute scale (days, or whatever the time unit is). Comparing RMST values at several horizons also reveals how quickly a covariate effect unfolds:

# RMST at several horizons: shows how the covariate gap widens over time
horizons = [180, 365, 730]

for tau in horizons:
    rmst = aft.predict(newdata, type="rmst", tau=tau)
    print(f"tau={tau:4d}: subject_1={rmst[0]:.1f}  subject_2={rmst[1]:.1f}")

tau= 180: subject_1=157.7  subject_2=140.2
tau= 365: subject_1=261.9  subject_2=212.6
tau= 730: subject_1=375.6  subject_2=277.2

Subject-specific vs. group RMST

The type="rmst" prediction from AFT.predict() gives a subject-specific RMST for covariate-adjusted predictions from a parametric model. Greenwood also provides gw.rmst_test() and gw.rmst_diff() for nonparametric group comparisons of RMST from Kaplan-Meier curves. The two approaches are complementary: parametric RMST leverages distributional assumptions for individual-level predictions; nonparametric RMST makes no assumptions and tests group-level differences.

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.