Royston-Parmar flexible parametric survival model (proportional hazards scale).
RoystonParmar(
df=3,
*,
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 the log baseline cumulative hazard as a smooth spline function on the log-time scale, combined with proportional-hazards covariate effects. This allows flexible baseline shapes while maintaining interpretable proportional-hazards covariate coefficients. It’s a key advantage over fully parametric AFT models.
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) approaches a Weibull fit; 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 monotone increasing log cumulative hazard (valid hazard functions). The flexible baseline makes this model useful when baseline hazard shape is unknown but important, yet you want interpretable proportional-hazards effects of covariates. 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 Weibull model; df=3 (two internal knots) is a common flexible default.
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
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
rp = gw.RoystonParmar(df=3).fit(y, lung[["age", "sex"]])
rp
RoystonParmar (flexible parametric survival, df=3)
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
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.
|
|
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
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
rp = gw.RoystonParmar(df=3).fit(y, lung[["age", "sex"]])
rp
RoystonParmar (flexible parametric survival, df=3)
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:
rp_flexible = gw.RoystonParmar(df=5).fit(y, lung[["age", "sex"]])
rp_flexible
RoystonParmar (flexible parametric survival, df=5)
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:
rp_univariate = gw.RoystonParmar(df=3).fit(y)
rp_univariate
RoystonParmar (flexible parametric survival, df=3)
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
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
rp = gw.RoystonParmar(df=3).fit(y, lung[["age", "sex"]])
rp.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:
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):
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):
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 |
to_frame()
Return the coefficient table as a DataFrame.
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
lung = gw.load_dataset("lung", backend="polars")
y = gw.Surv.right(lung["time"], event=(lung["status"] == 2))
rp = gw.RoystonParmar(df=3).fit(y, lung[["age", "sex"]])
rp.to_frame(format="polars")
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=:
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 |