Royston-Parmar flexible parametric survival model (hazard or odds scale).
RoystonParmar(
df=3,
*,
scale="hazard",
conf_level=0.95,
)
The Royston-Parmar model offers a middle ground between rigid parametric models (like Weibull) and fully non-parametric methods (like Kaplan-Meier). It models a transformed baseline survival function as a smooth spline on the log-time scale, combined with covariate effects that are proportional on that transformed scale. This allows flexible baseline shapes while keeping interpretable, constant-over-time covariate effects. It’s a key advantage over fully parametric AFT models.
Two scales are available via scale=:
"hazard" (default): proportional hazards, on the log cumulative hazard scale. exp(coef) is a hazard ratio."odds": proportional odds, on the log odds of failure scale. exp(coef) is an odds ratio. Useful when the proportional-hazards assumption fails but the odds ratio still looks roughly constant over time, which shows up as survival curves that converge rather than stay parallel.
The model uses restricted cubic splines with a fixed number of degrees of freedom (controlled by knots placed at quantiles of event times). A low df value (e.g., df=1) reproduces a 2-parameter AFT distribution exactly (Weibull under "hazard", log-logistic under "odds"); higher df values (e.g., df=3 or 4) provide greater flexibility. Call fit() with a right-censored Surv response and a design matrix. The model reports spline and covariate coefficients, fitted knot locations, log-likelihood, and supports predictions of survival at specified times and covariate values.
The implementation uses maximum likelihood estimation with constraints that ensure the transformed baseline is monotone increasing in time (required for both scales: cumulative hazard and the odds of failure both increase over time for any valid model). Results can be exported to tidy DataFrames or accessed as coefficient arrays.
Parameters
df: int = 3-
Spline degrees of freedom: the number of spline terms beyond the intercept, equal to one more than the number of internal knots. df=1 is a 2-parameter AFT model (Weibull or log-logistic, depending on scale); df=3 (two internal knots) is a common flexible default.
scale: str = "hazard"-
"hazard" (default, proportional hazards) or "odds" (proportional odds).
conf_level: float = 0.95-
Confidence level for coefficient intervals (default 0.95).
Returns
Fitted estimator-
Call
fit() to produce a fitted estimator with cached results (coef_, std_error_, z_, p_value_, conf_low_, conf_high_, knots_, loglik_, aic_, bic_), accessible as arrays or exported to DataFrames.
Details
Call fit(surv, covariates) with a right-censored Surv response and a covariate design (a dataframe, a 2-D array, or a formula string with data). Results are exposed as arrays (coef_, std_error_, …), the fitted knots_, and tidy frames via to_frame() (optionally format=).
Examples
Build a Surv response from the bundled lung dataset and fit a flexible model with three spline degrees of freedom and age and sex as covariates. Printing the fitted object reports the spline and covariate coefficients and the log-likelihood.
import greenwood as gw
# Load data and build a right-censored response
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
# Fit a flexible model with three spline terms
rp = gw.RoystonParmar(df=3).fit(y, lung[["age", "sex"]])
rp
RoystonParmar (flexible parametric survival, df=3, scale='hazard', exp(coef) = hazard ratio)
coef se(coef) z p
gamma0 -7.228 1.327 -5.447 5.11e-08
gamma1 1.028 0.2972 3.457 0.0005464
gamma2 -0.09642 0.1302 -0.740 0.4591
gamma3 0.1172 0.1857 0.631 0.5277
age 0.01615 0.009194 1.756 0.07906
sex -0.5101 0.1672 -3.051 0.002278
n = 228, events = 165
Log-likelihood = -1147
Fit on the proportional-odds scale instead, giving odds-ratio coefficients:
# Fit a flexible proportional-odds 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.988e-06
gamma1 1.012 0.3039 3.330 0.0008675
gamma2 0.198 0.1851 1.069 0.285
gamma3 -0.4249 0.2821 -1.506 0.1319
age 0.02429 0.01376 1.765 0.07753
sex -0.892 0.2569 -3.472 0.0005158
n = 228, events = 165
Log-likelihood = -1146
Methods
|
Name
|
Description
|
|
fit()
|
Fit the Royston-Parmar flexible parametric model to survival data.
|
|
predict()
|
Predict survival probability, hazard, or cumulative hazard from the fitted model.
|
|
predict_expectation()
|
Predict the expected (restricted mean) survival time for each subject.
|
|
predict_median()
|
Predict the median survival time for each subject.
|
|
predict_quantile()
|
Predict survival-time quantiles for each subject.
|
|
to_frame()
|
Return the coefficient table as a DataFrame.
|
fit()
Fit the Royston-Parmar flexible parametric model to survival data.
fit(
surv,
covariates=None,
*,
data=None,
)
Fits a flexible parametric survival model to a right-censored response and optional covariates. The model uses restricted cubic splines on the log-time scale to flexibly estimate the baseline cumulative hazard, combined with proportional-hazards covariate effects. This combines the interpretability of proportional-hazards regression with the flexibility of non-parametric methods.
The spline flexibility is controlled by df (degrees of freedom): df=1 recovers a Weibull model; higher df values provide more flexibility to capture non-standard baseline hazard shapes. An intercept is added automatically. Covariates are optional; if omitted, the fit is a flexible univariate survival model (baseline hazard only).
Parameters
surv: Surv-
A right-censored Surv response. Built with Surv.right().
covariates: Any = None-
Optional. A dataframe (pandas or polars), a 2-D array, or a formula string (e.g., "age + sex") evaluated against the data argument. If None (default), fits a univariate model with no covariates.
data: Any = None-
A dataframe to evaluate the formula string (ignored if
covariates is a dataframe, array, or None).
Returns
RoystonParmar-
The fitted estimator object itself (for method chaining) with cached coefficient arrays (
coef_, std_error_, z_, p_value_), fitted knot locations (knots_), and model fit statistics.
Details
The Royston-Parmar model parameterizes the log cumulative hazard as a restricted cubic spline in log-time, with proportional-hazards covariate effects added linearly. Knots are placed at quantiles of observed event times. Maximum likelihood estimation is used; constraints ensure that the log cumulative hazard is monotone increasing (required for a valid hazard function).
The model is useful when baseline hazard shape is unknown but important, yet you want interpretable proportional-hazards effects of covariates.
Examples
Fit a flexible Royston-Parmar model with three degrees of freedom (two internal knots) on the bundled lung dataset with age and sex as covariates:
import greenwood as gw
# Load data and build a right-censored response
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
# Fit the model with three spline degrees of freedom
rp = gw.RoystonParmar(df=3).fit(y, lung[["age", "sex"]])
rp
RoystonParmar (flexible parametric survival, df=3, scale='hazard', exp(coef) = hazard ratio)
coef se(coef) z p
gamma0 -7.228 1.327 -5.447 5.11e-08
gamma1 1.028 0.2972 3.457 0.0005464
gamma2 -0.09642 0.1302 -0.740 0.4591
gamma3 0.1172 0.1857 0.631 0.5277
age 0.01615 0.009194 1.756 0.07906
sex -0.5101 0.1672 -3.051 0.002278
n = 228, events = 165
Log-likelihood = -1147
Fit a more flexible model with five degrees of freedom:
# Increase spline flexibility to five degrees of freedom
rp_flexible = gw.RoystonParmar(df=5).fit(y, lung[["age", "sex"]])
rp_flexible
RoystonParmar (flexible parametric survival, df=5, scale='hazard', exp(coef) = hazard ratio)
coef se(coef) z p
gamma0 -7.224 1.395 -5.178 2.239e-07
gamma1 1.025 0.3488 2.937 0.003312
gamma2 0.04022 0.4412 0.091 0.9274
gamma3 -0.3697 1.258 -0.294 0.7689
gamma4 0.6169 1.604 0.385 0.7005
gamma5 -0.3493 1.001 -0.349 0.7272
age 0.01616 0.009193 1.758 0.07867
sex -0.5075 0.1673 -3.034 0.002417
n = 228, events = 165
Log-likelihood = -1147
Fit a univariate flexible model without covariates:
# Fit a baseline-only model with no covariates
rp_univariate = gw.RoystonParmar(df=3).fit(y)
rp_univariate
RoystonParmar (flexible parametric survival, df=3, scale='hazard', exp(coef) = hazard ratio)
coef se(coef) z p
gamma0 -6.895 1.16 -5.946 2.753e-09
gamma1 1.029 0.2972 3.463 0.0005346
gamma2 -0.08608 0.1295 -0.665 0.5063
gamma3 0.1022 0.1847 0.553 0.5801
n = 228, events = 165
Log-likelihood = -1153
predict()
Predict survival probability, hazard, or cumulative hazard from the fitted model.
predict(
newdata=None,
*,
type="survival",
times=None,
format=None,
)
Generates predictions from a fitted Royston-Parmar flexible parametric model. Pass newdata=None to predict for a baseline subject (all covariates set to 0, or training data mean if covariates are centered).
The Royston-Parmar model flexibly estimates the baseline cumulative hazard via splines, then multiplies by \exp(\eta) for each subject’s covariate-adjusted log-hazard \eta. This produces smooth, covariate-adjusted survival and hazard curves.
Three prediction types are available:
Survival (type="survival"): Survival probabilities S(t \mid x) at specified times. Useful for survival curves and prognosis.
Hazard (type="hazard"): Instantaneous hazard h(t \mid x) at specified times. Shows the rate of events at each time.
Cumulative hazard (type="cumhaz"): Cumulative hazard H(t \mid x) at specified times. Useful for risk quantification and comparisons.
Parameters
newdata: Any = None-
Covariate values for prediction. A DataFrame (Pandas or Polars), 2-D array, or None (the default). If None, uses baseline (all covariates 0 or the training data mean). Must have the same columns/features as the training data.
type: str = "survival"-
Prediction type (default "survival"):
"survival": Survival probabilities S(t \mid x) = \exp(-H(t \mid x)). Returns a frame with time column and one column per subject."hazard": Instantaneous hazard h(t \mid x) = dH(t \mid x)/dt. Returns a frame with time column and one column per subject."cumhaz": Cumulative hazard H(t \mid x). Returns a frame with time column and one column per subject.
times: Any = None-
Query times at which to evaluate curves. An array-like of floats. Required unless a default grid is used. If None, may raise an error or use a default grid.
format: str | None = None-
Output format for the returned frame:
None (default), "pandas", "polars", or "pyarrow". When None, a backend is auto-detected (Polars, then Pandas, then PyArrow).
Returns
DataFrame-
A DataFrame with columns
time (query times) and subject_1, subject_2, etc. containing predictions for each subject (one row per query time). All three prediction types return the same DataFrame structure with different values.
Raises
ValueError-
If
type= is not one of "survival", "hazard", or "cumhaz".
Details
The Royston-Parmar model represents log cumulative hazard as a smooth spline function in log-time, with proportional-hazards covariate effects: H(t \mid x) = \exp(\eta(t, x)), where \eta(t, x) = \text{spline}(\log t) + x^\top \beta. The spline basis and knot locations are fitted to the training data; predictions use these fixed basis functions.
Hazard is computed numerically as the derivative of cumulative hazard, so predictions may be slightly noisy if times are coarsely spaced. For smooth hazard predictions, use a fine query grid.
Predictions assume the model is well-specified and fit the training data adequately.
Examples
Read predicted survival probabilities off the fitted curves at chosen times. Here are the estimates at 180 and 365 days for the first two subjects; pass format= to choose the backend (here, Polars):
import greenwood as gw
# Load data and build a right-censored response
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"]])
# Predict survival probabilities at 180 and 365 days for two subjects
rp.predict(
lung[["age", "sex"]][:2], type="survival", times=[180, 365], format="polars"
)
shape: (2, 3)| time | subject_1 | subject_2 |
|---|
| f64 | f64 | f64 |
| 180.0 | 0.632099 | 0.659447 |
| 365.0 | 0.292051 | 0.327202 |
Predict the instantaneous hazard (force of mortality) at those same times:
# Predict instantaneous hazard at the same time points
rp.predict(
lung[["age", "sex"]][:2], type="hazard", times=[180, 365], format="polars"
)
shape: (2, 3)| time | subject_1 | subject_2 |
|---|
| f64 | f64 | f64 |
| 180.0 | 0.003529 | 0.003203 |
| 365.0 | 0.004607 | 0.004182 |
Predict cumulative hazard (total risk accumulated by time t):
# Predict cumulative hazard at the same time points
rp.predict(
lung[["age", "sex"]][:2], type="cumhaz", times=[180, 365], format="polars"
)
shape: (2, 3)| time | subject_1 | subject_2 |
|---|
| f64 | f64 | f64 |
| 180.0 | 0.45871 | 0.416354 |
| 365.0 | 1.230827 | 1.117177 |
Predict for a baseline subject (covariates all zero):
# Predict for a baseline subject with all covariates set to zero
rp.predict(type="survival", times=[180, 365], format="polars")
shape: (2, 2)| time | subject_1 |
|---|
| f64 | f64 |
| 180.0 | 0.793509 |
| 365.0 | 0.537618 |
predict_expectation()
Predict the expected (restricted mean) survival time for each subject.
predict_expectation(
newdata=None,
*,
tau,
ci=False,
format=None,
)
The restricted mean survival time (RMST) up to horizon \tau is \mathrm{RMST}_i(\tau) = \int_0^{\tau} S(t \mid x_i)\, dt, where S(t \mid x_i) = \exp(-\exp(s(\log t) + x_i^\top \beta)) is the smooth spline-based survival function. The integral is computed by numerical quadrature.
When ci=True, confidence intervals are computed via the delta method on the log cumulative hazard, integrated over [0, \tau].
Parameters
newdata: Any = None-
Covariate values for prediction. A DataFrame (Pandas or Polars), 2-D array, or None (the default). If None, uses baseline (all covariates 0).
tau: float-
The restriction time (time horizon). Must be positive.
ci: bool = False-
If True, include confidence intervals. Default is False.
format: str | None = None-
Output format:
None (auto-detect), "pandas", "polars", or "pyarrow".
Returns
DataFrame-
A DataFrame with a
tau column and columns subject_1, subject_2, etc. containing the RMST for each subject. With ci=True, additional subject_N_lower and subject_N_upper columns are included.
Examples
import greenwood as gw
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"]])
# Expected survival time up to one year for three subjects
rp.predict_expectation(lung[["age", "sex"]][:3], tau=365, format="polars")
shape: (1, 4)| tau | subject_1 | subject_2 | subject_3 |
|---|
| f64 | f64 | f64 | f64 |
| 365.0 | 231.690441 | 240.410706 | 256.879938 |
With confidence intervals:
rp.predict_expectation(lung[["age", "sex"]][:3], tau=365, ci=True, format="polars")
shape: (1, 10)| tau | subject_1 | subject_1_lower | subject_1_upper | subject_2 | subject_2_lower | subject_2_upper | subject_3 | subject_3_lower | subject_3_upper |
|---|
| f64 | f64 | f64 | f64 | f64 | f64 | f64 | f64 | f64 | f64 |
| 365.0 | 231.690441 | 228.687155 | 234.693728 | 240.410706 | 238.016844 | 242.804567 | 256.879938 | 254.425728 | 259.334147 |
predict_quantile()
Predict survival-time quantiles for each subject.
predict_quantile(
newdata=None,
*,
p=0.5,
ci=False,
format=None,
)
The quantile at failure probability p is the time t at which S(t \mid x) = 1 - p, found by solving \exp(s(\log t) + x^\top \beta) = -\log(1 - p) via root-finding on the smooth spline-based log cumulative hazard.
When ci=True, confidence intervals are computed via the delta method on the log cumulative hazard. At the quantile, the variance of the log-quantile is \text{Var}(\eta) / (\partial \eta / \partial \log t)^2, where the denominator is the spline derivative evaluated at the quantile.
Parameters
newdata: Any = None-
Covariate values for prediction. A DataFrame (Pandas or Polars), 2-D array, or None (the default). If None, uses baseline (all covariates 0).
p: Any = 0.5-
Failure probability or probabilities at which to compute quantiles. Can be a scalar (e.g., 0.5 for median) or array-like (e.g., [0.25, 0.5, 0.75] for quartiles). Must be in (0, 1). Default is 0.5.
ci: bool = False-
If True, include confidence intervals. Default is False.
format: str | None = None-
Output format:
None (auto-detect), "pandas", "polars", or "pyarrow".
Returns
DataFrame-
A DataFrame with a
p column listing the failure probabilities, and columns subject_1, subject_2, etc. containing the corresponding quantile for each subject. With ci=True, additional subject_N_lower and subject_N_upper columns are included.
Examples
import greenwood as gw
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"]])
# Predicted survival-time quartiles for three subjects
rp.predict_quantile(lung[["age", "sex"]][:3], p=[0.25, 0.5, 0.75], format="polars")
shape: (3, 4)| p | subject_1 | subject_2 | subject_3 |
|---|
| f64 | f64 | f64 | f64 |
| 0.25 | 127.493059 | 137.147511 | 158.35802 |
| 0.5 | 241.733051 | 258.958787 | 297.320421 |
| 0.75 | 398.376214 | 428.064439 | 494.974733 |
With confidence intervals at the median:
rp.predict_quantile(lung[["age", "sex"]][:3], p=0.5, ci=True, format="polars")
shape: (1, 10)| p | subject_1 | subject_1_lower | subject_1_upper | subject_2 | subject_2_lower | subject_2_upper | subject_3 | subject_3_lower | subject_3_upper |
|---|
| f64 | f64 | f64 | f64 | f64 | f64 | f64 | f64 | f64 | f64 |
| 0.5 | 241.733051 | 197.459146 | 295.933964 | 258.958787 | 220.380221 | 304.290709 | 297.320421 | 248.181362 | 356.188845 |
to_frame()
Return the coefficient table as a DataFrame.
to_frame(
*,
format=None,
)
Exports one row per spline or covariate term with coefficient estimates, standard errors, Wald statistics, p-values, and confidence limits.
Parameters
format: str | None = None-
Output format:
None (default), "pandas", "polars", or "pyarrow". When None, a backend is auto-detected (Polars, then Pandas, then PyArrow).
Returns
pandas.DataFrame, polars.DataFrame, or pyarrow.Table-
A tidy table with columns
term, estimate, std_error, statistic, p_value, conf_low, and conf_high.
Raises
ImportError-
If the requested (or, when auto-detecting, any) DataFrame library is not installed.
Examples
Fit a Royston-Parmar model and export the coefficient table as a Polars frame:
import greenwood as gw
# Load data and build a right-censored response
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"]])
# Export the coefficient table as a Polars DataFrame
rp.to_frame(format="polars")
shape: (6, 7)| term | estimate | std_error | statistic | p_value | conf_low | conf_high |
|---|
| str | f64 | f64 | f64 | f64 | f64 | f64 |
| "gamma0" | -7.228454 | 1.326942 | -5.447451 | 5.1097e-8 | -9.829214 | -4.627695 |
| "gamma1" | 1.027512 | 0.297232 | 3.456931 | 0.000546 | 0.444947 | 1.610076 |
| "gamma2" | -0.09642 | 0.130227 | -0.7404 | 0.459057 | -0.35166 | 0.15882 |
| "gamma3" | 0.117232 | 0.185653 | 0.631459 | 0.52774 | -0.246641 | 0.481106 |
| "age" | 0.016147 | 0.009194 | 1.756155 | 0.079062 | -0.001874 | 0.034168 |
| "sex" | -0.510127 | 0.167177 | -3.051428 | 0.002278 | -0.837787 | -0.182467 |
Request a different backend with format=:
# Export the same table as a Pandas DataFrame
rp.to_frame(format="pandas")
|
term |
estimate |
std_error |
statistic |
p_value |
conf_low |
conf_high |
| 0 |
gamma0 |
-7.228454 |
1.326942 |
-5.447451 |
5.109676e-08 |
-9.829214 |
-4.627695 |
| 1 |
gamma1 |
1.027512 |
0.297232 |
3.456931 |
5.463653e-04 |
0.444947 |
1.610076 |
| 2 |
gamma2 |
-0.096420 |
0.130227 |
-0.740400 |
4.590572e-01 |
-0.351660 |
0.158820 |
| 3 |
gamma3 |
0.117232 |
0.185653 |
0.631459 |
5.277404e-01 |
-0.246641 |
0.481106 |
| 4 |
age |
0.016147 |
0.009194 |
1.756155 |
7.906192e-02 |
-0.001874 |
0.034168 |
| 5 |
sex |
-0.510127 |
0.167177 |
-3.051428 |
2.277558e-03 |
-0.837787 |
-0.182467 |