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.
generalized gamma: when none of the above is clearly right and you want a data-driven test of which simpler family is adequate (see below).

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.

To compare predicted curves visually, plot_predicted_survival() draws one step-function curve per subject. This is especially useful for parametric models, where the smooth survival function can be evaluated at many time points.

# Visualize predicted survival for the two new subjects
gw.plot_predicted_survival(aft, newdata)

The older subject’s curve drops faster, reflecting the shorter predicted survival. You can pass labels= to name the curves, or switch to type="cumulative_hazard" for the cumulative-hazard scale.

Median survival time

The predict_median() method is a convenient shortcut for the predicted median survival time. For AFT models this is computed in closed form rather than by inverting the survival curve.

# Predicted median survival time for each subject
aft.predict_median(newdata, format="polars")

	p f64	subject_1 f64	subject_2 f64
PolarsRows1Columns3
0	0.5	319.808460567	250.278451947

This is equivalent to aft.predict(newdata, type="quantile", p=0.5) but returns a single-row DataFrame and supports confidence intervals via the delta method on the log scale.

# Median survival with confidence intervals
aft.predict_median(newdata, ci=True, format="polars")

	p f64	subject_1 f64	subject_1_lower f64	subject_1_upper f64	subject_2 f64	subject_2_lower f64	subject_2_upper f64
PolarsRows1Columns7
0	0.5	319.808460567	253.670571603	403.19005395	250.278451947	211.881100683	295.634217997

Survival quantiles

predict_median() is really a shortcut for predict_quantile(p=0.5). Call predict_quantile() directly with a list of probabilities to get several survival-time quantiles at once, still in closed form and still with confidence intervals available.

# Predicted survival-time quantiles at the 25th, 50th, and 75th percentiles
aft.predict_quantile(newdata, 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

# Survival quantiles with confidence intervals
aft.predict_quantile(newdata, p=[0.25, 0.5, 0.75], ci=True, format="polars")

	p 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	0.25	164.781414313	130.703845308	207.74380768	128.956054541	109.171806696	152.325627891
1	0.5	319.808460567	253.670571603	403.19005395	250.278451947	211.881100683	295.634217997
2	0.75	539.364048781	427.821034859	679.988951953	422.100149962	357.342167003	498.593653508

Expected survival time

predict_expectation() computes the restricted mean survival time (RMST): the area under the predicted survival curve up to a time horizon tau. It answers a different question than the quantile methods above. Instead of “when does survival reach some probability”, it asks “how much event-free time is expected, on average, up to tau”. Unlike predict_median() and predict_quantile(), tau is required and must be specified explicitly, since RMST is only defined relative to a chosen horizon.

# Expected survival time (RMST) up to 365 days for each subject
aft.predict_expectation(newdata, tau=365, format="polars")

	tau f64	subject_1 f64	subject_2 f64
PolarsRows1Columns3
0	365	263.588858998	236.115103815

# Expected survival time with confidence intervals
aft.predict_expectation(newdata, tau=365, ci=True, format="polars")

	tau f64	subject_1 f64	subject_1_lower f64	subject_1_upper f64	subject_2 f64	subject_2_lower f64	subject_2_upper f64
PolarsRows1Columns7
0	365	263.588858998	237.722196949	285.681738357	236.115103815	215.461039459	255.214223481

Generalized gamma regression

The four standard AFT distributions each impose a fixed hazard shape. Picking among them by AIC works, but it treats the candidates as unrelated. The generalized gamma distribution offers a more principled approach: it is a three-parameter family that nests Weibull (Q = 1), log-normal (Q approaches 0), and exponential (Q = 1, sigma = 1) as special cases. Fitting it gives a formal likelihood ratio test of whether a simpler sub-model is adequate.

# Fit a generalized gamma AFT model
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

gg = gw.AFT("gengamma").fit(y, lung[["age", "sex"]])
gg

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

                 coef  se(coef)       z          p
(Intercept)     6.275    0.4809  13.049  6.446e-39
age          -0.01226  0.006947  -1.764    0.07765
sex            0.3821    0.1269   3.012   0.002596

Scale = 0.7541
Q = 1
n = 228, events = 165
Log-likelihood = -1147

The extra shape parameter Q controls the distributional family. When Q is near 1 the data favor a Weibull shape. When it is near 0 the data favor a log-normal shape. Values far from both suggest that neither simple family is adequate.

Testing nested sub-models

The test_distributions() method fits each nested sub-model (Weibull, log-normal, exponential, log-logistic) and compares it against the generalized gamma via a likelihood ratio test with one degree of freedom. A large p-value means the simpler model cannot be rejected.

# Compare the generalized gamma against its nested sub-models
gg.test_distributions(format="polars")

	dist str	loglik f64	aic f64	lr_statistic f64	df i64	p_value f64
PolarsRows4Columns6
0	weibull	-1147.05443143	2302.10886287	0	1	1
1	lognormal	-1158.75014259	2325.50028519	23.3914223228	1	1.321665652e-06
2	exponential	-1156.09903714	2318.19807428	18.0892114111	1	2.10792418187e-05
3	loglogistic	-1152.89722533	2313.79445065	11.6855877835	1	0.000629860851085

A non-significant p-value for a sub-model (say, Weibull) means the extra Q parameter does not improve the fit enough to justify the added complexity. You can then use the simpler model with confidence. The AIC column gives a penalty-adjusted comparison that also accounts for sample size.

Predictions

Predictions from the generalized gamma work the same way as from any other AFT model. All prediction types (survival, quantiles, mean, RMST) are supported.

# Predict survival probabilities at specific times
newdata = pd.DataFrame({"age": [50, 70], "sex": [1, 1]})
gg.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

# Predicted median survival time
gg.predict(newdata, type="quantile", p=0.5, format="polars")

	p f64	subject_1 f64	subject_2 f64
PolarsRows1Columns3
0	0.5	319.808460567	250.278451947

When to use the generalized gamma

Use "gengamma" when you want a data-driven answer to “which parametric family fits best?” rather than comparing AICs across separate fits. If test_distributions() shows that Weibull is adequate (large p-value for Q = 1), use the simpler Weibull model for final inference. If no sub-model is adequate, the generalized gamma itself may be the right choice.

Three-parameter (threshold) distributions

The Weibull, exponential, log-normal, and log-logistic distributions all put positive density arbitrarily close to time zero. Some failure processes have a genuine minimum lifetime instead: a component with a guaranteed burn-in period, or a disease with an incubation time before death is even possible. A three-parameter model captures this with a location (threshold) parameter \gamma: T = \gamma + T', where T' follows the ordinary 2-parameter distribution. Pass threshold=True to add it.

# Fit a Weibull AFT model with a threshold parameter
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

aft_threshold = gw.AFT("weibull", method="mps", threshold=True).fit(y, lung[["age", "sex"]])
aft_threshold

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

                 coef   se(coef)         z  p
(Intercept)     6.344    0.06134   103.422  0
age          -0.01354   2.68e-05  -505.186  0
sex              0.39  0.0003594  1085.191  0

Scale = 0.7872
Threshold = 2.311
n = 228, events = 165
Log-likelihood = -1147

The fitted threshold_ is the estimated minimum survival time before which no death can occur, here just over two days.

Why `method="mps"`

Ordinary maximum likelihood can behave badly with a threshold parameter: as \gamma approaches the smallest observed follow-up time, the likelihood can grow without bound for some distributions, pushing the estimate toward a meaningless boundary value. method="mps" (maximum product of spacings) replaces each observation’s density contribution with the gap between consecutive order statistics on the CDF scale, which stays bounded where the ordinary likelihood does not. Combining threshold=True with the default method="mle" raises a warning for exactly this reason:

import warnings

with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always")
    gw.AFT("weibull", threshold=True).fit(y, lung[["age", "sex"]])
print(caught[0].message)

threshold=True with method='mle' can drive the threshold toward the smallest observed time, where the likelihood is unbounded for some distributions. method='mps' is recommended instead.

When to use a threshold

Use threshold=True only when there is a substantive reason to believe survival time has a genuine floor. Reliability engineering with a burn-in guarantee is the classic case. Always pass method="mps" alongside it: method="mle" stays the default for backward compatibility with ordinary 2-parameter fits, not because it is the right choice once a threshold is added.

Predictions with a threshold

Every prediction type accounts for the threshold automatically: threshold_ is added back onto every time-scale prediction, and survival is exactly 1 for any query time at or before it.

# Survival at times before, at, and after the fitted threshold
query_times = [aft_threshold.threshold_ * 0.5, aft_threshold.threshold_, 100, 365]
aft_threshold.predict(
    lung[["age", "sex"]][:2], type="survival", times=query_times, format="polars"
)

	time f64	subject_1 f64	subject_2 f64
PolarsRows4Columns3
0	1.15528153244	1	1
1	2.31056306488	1	1
2	100	0.792986660706	0.811226908917
3	365	0.293005226994	0.330480255205

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

Royston-Parmar models also support predict_median(), which finds the time at which the smooth survival curve crosses 0.5 via root-finding on the spline. Confidence intervals use the delta method on the log cumulative hazard.

# Median survival time from the Royston-Parmar model, with confidence intervals
rp.predict_median(newdata, ci=True, format="polars")

	p f64	subject_1 f64	subject_1_lower f64	subject_1_upper f64	subject_2 f64	subject_2_lower f64	subject_2_upper f64
PolarsRows1Columns7
0	0.5	318.737847932	252.375332743	402.550497311	253.082870007	213.059801327	300.62423175

Survival quantiles

As with the AFT models, predict_quantile() generalizes predict_median() to arbitrary failure probabilities, root-finding on the spline for each requested p.

# Predicted survival-time quantiles at the 25th, 50th, and 75th percentiles
rp.predict_quantile(newdata, p=[0.25, 0.5, 0.75], format="polars")

	p f64	subject_1 f64	subject_2 f64
PolarsRows3Columns3
0	0.25	169.979666666	133.861669027
1	0.5	318.737847932	253.082870007
2	0.75	532.603653513	417.905940907

# Survival quantiles with confidence intervals
rp.predict_quantile(newdata, p=[0.25, 0.5, 0.75], ci=True, format="polars")

	p 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	0.25	169.979666666	131.548067949	219.639007478	133.861669027	106.244819946	168.657130238
1	0.5	318.737847932	252.375332743	402.550497311	253.082870007	213.059801327	300.62423175
2	0.75	532.603653513	420.202701188	675.070985821	417.905940907	352.302775853	495.725232429

Expected survival time

predict_expectation() works the same way for the Royston-Parmar model, integrating the smooth spline-based survival curve up to tau to give the restricted mean survival time. As with the AFT models, tau is required and must be specified explicitly.

# Expected survival time (RMST) up to 365 days for each subject
rp.predict_expectation(newdata, tau=365, format="polars")

	tau f64	subject_1 f64	subject_2 f64
PolarsRows1Columns3
0	365	264.580086123	237.537308382

# Expected survival time with confidence intervals
rp.predict_expectation(newdata, tau=365, ci=True, format="polars")

	tau f64	subject_1 f64	subject_1_lower f64	subject_1_upper f64	subject_2 f64	subject_2_lower f64	subject_2_upper f64
PolarsRows1Columns7
0	365	264.580086123	261.623116437	267.537055809	237.537308382	234.98053358	240.094083183

Selecting degrees of freedom with AIC and BIC

The loglik sweep above shows how fit improves with spline complexity, but it does not penalize for the additional parameters. AIC and BIC do: they balance log-likelihood against model size, so you can pick the df that fits best without overfitting. The glance() function computes both automatically for a fitted RoystonParmar.

# Compare AIC and BIC across spline degrees of freedom
rows = []
for df in (1, 2, 3, 4, 5):
    m = gw.RoystonParmar(df=df).fit(y, lung[["age", "sex"]])
    g = gw.glance(m, format="polars")
    rows.append({"df": df, "aic": g["aic"][0], "bic": g["bic"][0]})

pd.DataFrame(rows)

	df i64	aic f64	bic f64
PandasRows5Columns3
0	1	2302.10886287	2314.53264476
1	2	2303.62012843	2319.1498558
2	3	2305.14107377	2323.77674661
3	4	2306.87778236	2328.61940068
4	5	2309.02752368	2333.87508747

Lower values are better. AIC and BIC often agree, but BIC penalizes complexity more heavily and may favor a simpler model. When they disagree, AIC is oriented toward prediction accuracy while BIC favors the most parsimonious model that explains the data.

Summarizing with tidy and glance

tidy() returns the coefficient table from to_frame(), giving one row per parameter (spline gammas and covariate effects) with the estimate, standard error, z-statistic, and p-value:

# Tidy: one row per parameter (spline + covariates)
gw.tidy(rp, 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 rows are the restricted cubic spline coefficients that define the baseline log cumulative hazard shape. The named rows (age, sex) are the covariate effects on the log cumulative hazard scale. A positive coefficient means higher cumulative hazard, and therefore shorter survival.

glance() returns a single-row summary of the overall fit:

# Glance: model-level summary
gw.glance(rp, format="polars")

	scale str	n i64	nevent i64	loglik f64	aic f64	bic f64	df i64	n_knots i64
PolarsRows1Columns8
0	hazard	228	165	-1146.57053688	2305.14107377	2323.77674661	6	3

The df column is the total number of estimated parameters (spline coefficients plus covariate effects), and n_knots is the number of interior knots used by the spline. These summaries are useful for building comparison tables across models with different df values or different covariate sets.

Proportional odds

Everything above assumes proportional hazards: the ratio of two subjects’ hazards stays constant over time. That assumption can fail even when the data are otherwise well-behaved, for example when survival curves for two groups converge over time instead of staying a fixed multiple apart. A proportional-odds model relaxes the assumption by putting the spline on the log odds of failure scale instead of the log cumulative hazard, so it is the odds ratio between groups that stays constant. Pass scale="odds":

# Fit a proportional-odds Royston-Parmar model
rp_odds = gw.RoystonParmar(df=3, scale="odds").fit(y, lung[["age", "sex"]])
rp_odds

RoystonParmar (flexible parametric survival, df=3, scale='odds', exp(coef) = odds ratio)

           coef  se(coef)       z          p
gamma0   -7.252     1.525  -4.755  1.985e-06
gamma1    1.012     0.304   3.330  0.0008687
gamma2    0.198    0.1852   1.069     0.2852
gamma3  -0.4249    0.2822  -1.506     0.1321
age     0.02429   0.01376   1.765    0.07751
sex      -0.892    0.2569  -3.472  0.0005161

n = 228, events = 165
Log-likelihood = -1146

The coefficients are on a different scale now: exp(coef) is an odds ratio, not a hazard ratio. A coefficient of -0.89 for sex, for instance, means the odds of having died by any given time are exp(-0.89) ≈ 0.41 times as high for that group, holding age fixed, rather than the hazard being scaled by that factor. Everything else works exactly as before: predict(), predict_quantile(), predict_expectation(), tidy(), and glance() all behave the same way, just evaluated on the odds scale internally.

# Predicted survival at 6 months, 1 year, and 2 years from the proportional-odds model
rp_odds.predict(newdata, type="survival", times=[180, 365, 730], format="polars")

	time f64	subject_1 f64	subject_2 f64
PolarsRows3Columns3
0	180	0.727063061609	0.62102287528
1	365	0.418814555531	0.307140481762
2	730	0.113734981618	0.0731672389127

Choosing hazard vs. odds

Just as scale="hazard" with df=1 is exactly a Weibull model, scale="odds" with df=1 is exactly a log-logistic model, so the two scales are directly comparable at every df via AIC or BIC. There is no substitute for checking the proportional-hazards assumption first (see Cox model diagnostics); reach for scale="odds" when that check fails but the data still look well-behaved on the odds scale, not as a default choice.

Rank-based regression (Buckley-James)

Every model so far still commits to a distributional family (or, for RoystonParmar, a smooth spline shape). BuckleyJames fits the same log-linear form, \log(T) = X\beta + \varepsilon, but makes no assumption about \varepsilon at all. It alternates two steps until the coefficients stabilize: impute each censored observation’s log-time using the Kaplan-Meier curve of the current residuals, then refit by ordinary least squares.

# Fit a Buckley-James rank-based AFT model
bj = gw.BuckleyJames().fit(y, lung[["age", "sex"]])
bj

BuckleyJames (rank-based accelerated failure time regression)

                 coef  se(coef)   z   p
(Intercept)     6.379        NA  NA  NA
age          -0.02306        NA  NA  NA
sex            0.4912        NA  NA  NA

n = 228, events = 165
iterations = 13, converged = True

The iterations and converged lines report the iterative fit’s status. This is a genuinely different estimation strategy from everything else on this page, worth sanity-checking against: on data with no censoring at all, every residual is already exact, so there is nothing to impute and the fit converges in a single iteration to plain ordinary least squares of log(time) on the covariates.

Standard errors via bootstrap

Buckley-James has no closed-form variance formula, so by default std_error_, p_value, and the confidence limits are nan. Pass n_boot= to estimate them by case-resampling bootstrap instead (refitting the whole iterative algorithm on each resample, so larger values of n_boot cost proportionally more time):

# Bootstrap standard errors (200 resamples)
bj_boot = gw.BuckleyJames(n_boot=200, seed=0).fit(y, lung[["age", "sex"]])
bj_boot.to_frame(format="polars")

	term str	estimate f64	std_error f64	statistic f64	p_value f64	conf_low f64	conf_high f64
PolarsRows3Columns7
0	(Intercept)	6.37909324299	0.562108221274	11.3485143991	7.5433445446e-30	5.27738137388	7.4808051121
1	age	-0.023064121313	0.0090905534195	-2.53715260762	0.0111758224591	-0.0408812786147	-0.00524696401122
2	sex	0.491199287137	0.150234189599	3.26955727221	0.00107715929125	0.196745686277	0.785652887997

Predicting survival

predict() supports the linear predictor (type="lp") and survival probabilities (type="survival"). Since there is no fitted distribution to evaluate, survival predictions reuse the empirical Kaplan-Meier curve of the fitted residuals, shifted by each subject’s linear predictor, rather than a smooth closed-form curve.

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

	time f64	subject_1 f64	subject_2 f64
PolarsRows3Columns3
0	180	0.745255252889	0.594394467194
1	365	0.498909301396	0.263334089298
2	730	0.156168564262	0.0265384488288

When to reach for Buckley-James

Buckley-James predates the maximum-likelihood AFT and is largely superseded by it in modern practice: AFT is more efficient when its distributional assumption is reasonable, has closed-form standard errors, and converges directly rather than iteratively. Reach for BuckleyJames only when you specifically want to avoid committing to any distributional family and can tolerate the computational cost of bootstrap inference, not as a default alternative to AFT or RoystonParmar.

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 behavior, 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 rmst_test() and 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.

Residuals

After fitting a parametric model, residuals tell you whether the distributional assumption holds and whether any observations are unusual or overly influential. The AFT residuals() method supports seven types, each answering a different diagnostic question.

Response and Cox-Snell residuals

Response residuals are the simplest: the difference between the observed time and the predicted time on the original scale. Cox-Snell residuals are the estimated cumulative hazard at each observation’s event or censoring time. Under a correctly specified model, the Cox-Snell residuals for uncensored observations follow an Exp(1) distribution, making them the standard tool for checking overall model fit.

# Fit a Weibull AFT for residual diagnostics
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"]])

# Response residuals: time - predicted time (one per subject)
aft.residuals("response")[:5]

array([  -8.16364231,  116.86052814,  618.27998426, -176.94785081,
        510.02264437])

# Cox-Snell residuals: cumulative hazard at observed time (always non-negative)
aft.residuals("cox_snell")[:5]

array([0.96568596, 1.48239478, 3.51167345, 0.44462163, 3.13585876])

Martingale and deviance residuals

Martingale residuals are the event indicator minus the Cox-Snell residual: they measure how much “more” or “less” each subject experienced the event than the model expected. Deviance residuals are a likelihood-ratio-based transformation that is more symmetrically distributed, making them better suited for identifying outliers.

# Martingale residuals: event - Cox-Snell (range: (-inf, 1])
aft.residuals("martingale")[:5]

array([ 0.03431404, -0.48239478, -3.51167345,  0.55537837, -2.13585876])

# Deviance residuals: symmetric, better for spotting outliers
aft.residuals("deviance")[:5]

array([-0.03471458,  0.42127403,  2.65015979, -0.71435742,  1.40922369])

Score, dfbeta, and dfbetas residuals

Score residuals are the per-observation derivatives of the log-likelihood with respect to the coefficients. They have one row per observation and one column per term (including the intercept). Dfbeta residuals approximate the change in each coefficient if observation i were deleted. Dfbetas standardizes dfbeta by the coefficient standard error so that values are comparable across covariates on different scales.

# Score residuals: derivative of each observation's log-likelihood contribution
aft.residuals("score", format="polars")

	(Intercept) f64	age f64	sex f64
PolarsRows228Columns3
0	-0.0455062483222	-3.36746237584	-0.0455062483222
1	0.639737389142	43.5021424617	0.639737389142
2	4.65707525874	260.79621449	4.65707525874
3	-0.736526025753	-41.9819834679	-0.736526025753
4	2.83251138002	169.950682801	2.83251138002
5	6.33854022688	469.051976789	6.33854022688
6	-0.614144062977	-41.7617962825	-1.22828812595
7	-0.411225275382	-29.1969945521	-0.822450550763
225	0.189838671143	14.2379003357	0.379677342286
226	0.53188538103	35.104435148	0.53188538103
227	0.287822663117	16.6937144608	0.575645326234

# Dfbeta: approximate leave-one-out influence on each coefficient
aft.residuals("dfbeta", format="polars")

	(Intercept) f64	age f64	sex f64
PolarsRows228Columns3
0	0.00101257876635	-2.31564602155e-05	0.00022239890343
1	-0.00229865735491	0.000139737753364	-0.00326689452397
2	0.15705271881	-0.00168789647631	-0.0258254409606
3	-0.0225478235422	0.000231292332548	0.00405741412322
4	0.0602888767913	-0.000478169626026	-0.0152931625417
5	-0.141041537812	0.00322545056998	-0.0309778205803
6	0.0168494651195	-0.000156604734966	-0.00684327112463
7	0.0151186299984	-0.000164577891734	-0.00462730373423
225	-0.00934076513693	0.000112732991416	0.0021639228683
226	0.00139689354712	6.46872192033e-05	-0.00275503427957
227	0.00105385461693	-6.59283707542e-05	0.00310189636853

# Dfbetas: standardized dfbeta, comparable across covariates
aft.residuals("dfbetas", format="polars")

	(Intercept) f64	age f64	sex f64
PolarsRows228Columns3
0	0.0021035668086	-0.00332831548278	0.00174467294062
1	-0.00477531178493	0.0200847333194	-0.0256281051207
2	0.326266851988	-0.242604090744	-0.202595189674
3	-0.0468416431251	0.0332440210767	0.0318295662456
4	0.125246236994	-0.068728093797	-0.119971665559
5	-0.293004660411	0.463599227648	-0.243014531512
6	0.0350036725496	-0.0225090518677	-0.0536840324855
7	0.0314079746809	-0.0236550465867	-0.0363002312
225	-0.019404834628	0.0162032952039	0.0169755228811
226	0.00290195587593	0.00929759864885	-0.0216126684256
227	0.00218931471498	-0.00947599136885	0.0243337290576

Use martingale and deviance residuals to check whether the distributional assumption is reasonable and to spot outliers. Use dfbeta and dfbetas to identify observations that disproportionately influence the coefficient estimates. Large values in either direction warrant closer investigation.

Residuals across distributions

All seven residual types work for every AFT distribution (Weibull, exponential, log-normal, log-logistic, and generalized gamma). The underlying math adapts to each distribution automatically via the standardized error distribution, so you can compare residual patterns across distributional choices to see which assumption fits best.

Weibull diagnostic plots

AIC and likelihood ratio tests give you a number, but a picture often tells you more. The Weibull plot is the standard visual diagnostic for checking whether a parametric distribution fits the data. It transforms the Kaplan-Meier survival estimates onto axes where the assumed distribution appears as a straight line. Points that fall on a line confirm the assumption. Systematic curvature suggests the distribution is wrong.

The plot_weibull() function draws this plot. By default it uses the Weibull (complementary log-log) transform: \log(t) on the x-axis and \log(-\log(S(t))) on the y-axis.

# Load data and build the response
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))

# Fit an unstratified KM and draw the Weibull diagnostic plot
km = gw.KaplanMeier().fit(y)
gw.plot_weibull(km)

If the points fall approximately on a straight line, the Weibull assumption is reasonable. The slope of the line is the reciprocal of the Weibull scale parameter, and the intercept reflects the location.

Checking the proportional hazards assumption

The Weibull plot is also used to check the proportional hazards assumption between groups. If two groups have proportional hazards, their transformed curves will be parallel (shifted vertically but with the same slope). Fit a stratified Kaplan-Meier and pass it to plot_weibull():

# Stratified KM by sex: parallel lines suggest proportional hazards
km_sex = gw.KaplanMeier().fit(y, by=lung["sex"])
gw.plot_weibull(km_sex, title="Weibull plot by sex")

Lines that are roughly parallel indicate that the PH assumption holds between the groups. If the lines cross or diverge, the assumption may be violated and a time-varying or stratified model may be more appropriate.

Overlaying a fitted parametric model

When you have both a Kaplan-Meier estimate and a fitted AFT model, you can overlay the theoretical parametric line on the empirical points to visualize how well the model fits. Pass the fitted AFT as the aft= argument:

# Fit a Weibull AFT and overlay the theoretical line on the KM scatter
km_overall = gw.KaplanMeier().fit(y)
aft = gw.AFT("weibull").fit(y, lung[["age", "sex"]])

gw.plot_weibull(km_overall, aft=aft, title="Weibull fit: KM vs. parametric")

The dashed line is the fitted Weibull model. When it tracks the empirical points closely, the parametric assumption is well supported. Gaps between the line and the points indicate regions where the model fits poorly.

AFT overlay requires an unstratified KM

The parametric overlay draws a single theoretical line from the AFT model’s fitted parameters. This is only meaningful for an unstratified (single-population) Kaplan-Meier. If you pass a stratified KM together with an AFT, plot_weibull() raises an error.

Other distributions

The dist= parameter switches the y-axis transform to match other AFT distributions. This lets you check log-normal or log-logistic assumptions the same way:

dist="lognormal": plots \Phi^{-1}(1 - S(t)) vs \log(t), where \Phi^{-1} is the standard normal quantile function (probit). Linear points indicate a log-normal distribution.
dist="loglogistic": plots \log((1 - S(t)) / S(t)) vs \log(t) (the logit of the failure probability). Linear points indicate a log-logistic distribution.

# Check the log-normal assumption
gw.plot_weibull(km_overall, dist="lognormal", title="Log-normal diagnostic")

# Check the log-logistic assumption
gw.plot_weibull(km_overall, dist="loglogistic", title="Log-logistic diagnostic")

Comparing the three plots side by side tells you which distributional family linearizes the data best. This visual approach complements the AIC comparison from compare_distributions() and the likelihood ratio tests from test_distributions().

Weibull plot + AIC = confidence in your distribution choice

Use compare_distributions() or test_distributions() for a numeric answer and plot_weibull() for the visual confirmation. If AIC picks Weibull and the Weibull plot is linear, you can proceed with confidence. If the numbers and the picture disagree, investigate further.

Mixture cure models

Standard survival models assume that every subject will eventually experience the event if followed long enough. In many settings this is not realistic: some patients are cured and will never relapse, some parolees will never reoffend, and some machines will never fail in their operating lifetime. A mixture cure model handles this by splitting the population into two groups:

Susceptible subjects who will eventually experience the event. Their time-to-event follows a standard survival distribution (the latency submodel).
Cured subjects who will never experience the event. They contribute only through the censoring mechanism.

The model estimates both the incidence (what fraction is cured, as a function of covariates) and the latency (how long susceptible subjects survive, also as a function of covariates). Greenwood’s MixtureCure uses a logistic regression for incidence and a Cox PH model for latency, fitted jointly via the EM algorithm.

Fitting a cure model

The fit() method takes a Surv response and two covariate frames: one for the cure (incidence) submodel and one for the latency submodel. The covariates can overlap or differ.

# Load the E1684 melanoma interferon trial
e1684 = gw.load_dataset("e1684", backend="polars")
y_mel = gw.Surv.right(e1684["FAILTIME"], event=(e1684["FAILCENS"] == 1))

# Fit a mixture cure model:
#   - Incidence (cure) submodel: treatment and age predict who is cured
#   - Latency submodel: treatment and age predict when susceptible subjects fail
mc = gw.MixtureCure().fit(
    y_mel,
    latency=e1684[["AGE", "TRT"]],
    cure=e1684[["AGE", "TRT"]],
    nboot=200,
)

mc

MixtureCure (logistic incidence + Cox PH latency)

Cure probability model:
                coef  se(coef)       z        p
(Intercept)    1.354    0.3189   4.244  2.2e-05
AGE          0.02083   0.01489   1.399    0.162
TRT          -0.6049     0.347  -1.743   0.0813

Failure time distribution model:
         coef  se(coef)        z      p
AGE  -0.00752  0.006843   -1.099  0.272
TRT    -0.142    0.1665  -0.8528  0.394

n = 284, events = 196, EM iterations = 38

The printed summary shows two coefficient tables. The cure probability model (logistic) estimates the log-odds of being susceptible. A positive coefficient increases the probability of being susceptible (i.e., decreases the cure probability). The failure time distribution model (Cox PH) estimates log hazard ratios among susceptible subjects, with the same interpretation as an ordinary Cox model.

Standard errors are computed by bootstrap resampling (controlled by the nboot parameter). More resamples give more stable SEs at the cost of longer fitting time.

Predicting cure probability

predict_cure_prob() returns the estimated probability that each subject is susceptible (uncured), based on their incidence covariates. Subtract from 1 to get the cure probability:

# Susceptibility probability for three subjects
suscept = mc.predict_cure_prob(e1684[["AGE", "TRT"]][:3])

# Cure probability = 1 - susceptibility
import numpy as np

for i, p in enumerate(suscept):
    print(f"Subject {i + 1}: susceptible = {p:.3f}, cured = {1 - p:.3f}")

Subject 1: susceptible = 0.627, cured = 0.373
Subject 2: susceptible = 0.655, cured = 0.345
Subject 3: susceptible = 0.863, cured = 0.137

Subjects with a low susceptibility probability are more likely to be in the cured fraction. These estimates depend only on the incidence covariates, not on the latency model.

Predicting population survival

predict_survival() returns the population survival curve, which blends the cured and susceptible groups. For a subject with incidence covariates z and latency covariates x:

S_{\text{pop}}(t \mid x, z) = (1 - \pi(z)) + \pi(z) \cdot S(t \mid x)

where \pi(z) is the susceptibility probability and S(t \mid x) is the latency survival. At t = 0 this equals 1, and as t \to \infty it plateaus at the cure fraction (1 - \pi(z)) rather than falling to zero.

# Population survival at 1, 2, 5, and 10 years for three subjects
times = [1.0, 2.0, 5.0, 10.0]
surv = mc.predict_survival(times, e1684[["AGE", "TRT"]][:3], e1684[["AGE", "TRT"]][:3])

pd.DataFrame(surv, columns=[f"t={t}" for t in times],
             index=["Subject 1", "Subject 2", "Subject 3"])

	t=1.0 f64	t=2.0 f64	t=5.0 f64	t=10.0 f64
PandasRows3Columns4
0	0.614460552634	0.503387627309	0.412570482772	0.373131382767
1	0.607737616311	0.490588293513	0.391317180961	0.344834454577
2	0.505919640588	0.350161225668	0.210825084399	0.137475488224

The survival curves flatten rather than approaching zero, reflecting the cured fraction. A subject with a high cure probability will have a higher survival plateau.

Summarizing with tidy and glance

tidy() returns a stacked coefficient table with both submodels. The submodel column distinguishes the incidence (“cure”) and latency rows:

# Tidy: stacked cure + latency coefficients
gw.tidy(mc, format="polars")

	submodel str	term str	coef f64	se f64	z f64	p f64
PolarsRows5Columns6
0	cure	(Intercept)	1.35354721221	0.318936283753	4.24394238335	2.19626764797e-05
1	cure	AGE	0.020825999797	0.0148915337583	1.39851274792	0.161959147254
2	cure	TRT	-0.604906265204	0.347002786785	-1.74323172101	0.0812931466248
3	latency	AGE	-0.0075197731212	0.00684347033122	-1.09882453744	0.271844607182
4	latency	TRT	-0.141986488987	0.166490219362	-0.852821802569	0.393758131342

The cure rows include an intercept (the baseline log-odds of susceptibility), while the latency rows do not (following Cox PH convention). Positive cure coefficients increase susceptibility. Negative latency coefficients decrease the hazard among susceptible subjects.

glance() returns a single-row model summary:

# Glance: model-level summary
gw.glance(mc, format="polars")

	n i64	nevent i64	n_iter i64
PolarsRows1Columns3
0	284	196	38

The n_iter column shows how many EM iterations were needed for convergence. If n_iter equals the maximum (emmax, default 50), the algorithm may not have converged and you should increase emmax.

When to consider a cure model

A cure model is appropriate when the Kaplan-Meier curve for your data levels off well above zero, suggesting that a fraction of subjects will never experience the event. If the KM curve continues to decline throughout follow-up, a standard Cox or AFT model is more appropriate.

Next steps

You can now fit, compare, predict from, and diagnose parametric survival models.

Cox model diagnostics covers residuals and influence diagnostics for the semiparametric alternative, with the same conceptual framework.
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.

Fitting a distribution to survival data

Comparing distributions

Predictions from a fitted distribution

The accelerated failure time model

Choosing a distribution

Model comparison and selection

Predicting survival and quantiles

Median survival time

Survival quantiles

Expected survival time

Generalized gamma regression

Testing nested sub-models

Predictions

Three-parameter (threshold) distributions

Why method="mps"

Predictions with a threshold

Flexible parametric (Royston-Parmar) models

Choosing the degrees of freedom

Survival quantiles

Expected survival time

Selecting degrees of freedom with AIC and BIC

Summarizing with tidy and glance

Proportional odds

Rank-based regression (Buckley-James)

Standard errors via bootstrap

Predicting survival

Survival predictions with confidence intervals

Mean survival and expected remaining lifetime

Unconditional mean survival

Conditional mean survival given landmark time

Expected remaining lifetime

Restricted mean survival time

Residuals

Response and Cox-Snell residuals

Martingale and deviance residuals

Score, dfbeta, and dfbetas residuals

Weibull diagnostic plots

Checking the proportional hazards assumption

Overlaying a fitted parametric model

Other distributions

Mixture cure models

Fitting a cure model

Predicting cure probability

Predicting population survival

Summarizing with tidy and glance

Next steps

Why `method="mps"`