I want to use a Royston-Parmar model (i.e. a flexible parametric proportional‐hazards or proportional‐odds models for censored survival data, ref) for prediction. Rather than a probability at a certain point in time, I am interested in the median survival. In R I am able to produce such predictions as follows:

# Load required packages
require(flexsurv) ; require(survival) ; require(rms); require(Hmisc)

# Use reproducible example dataset 'lung' (from survival package)    
l <- lung[complete.cases(lung[,c("time","status","sex","age","ph.karno")]), ]

# Fit Royston-Parmar model (with 2 knots and odds scale)
mod <- flexsurvspline(Surv(time,status)~age+sex+ph.karno, data=l, k=2, scale="odds")

#flexsurvspline(formula = Surv(time, status) ~ age + sex + ph.karno, 
#    data = l, k = 2, scale = "odds")

#          data mean  est      L95%     U95%     se       exp(est)  L95%     U95%   
#gamma0         NA    -3.2381  -6.8428   0.3666   1.8392       NA        NA       NA
#gamma1         NA     0.9674   0.3791   1.5558   0.3002       NA        NA       NA
#gamma2         NA     0.1506  -0.2267   0.5280   0.1925       NA        NA       NA
#gamma3         NA    -0.3657  -0.9417   0.2102   0.2939       NA        NA       NA
#age       62.4493     0.0143  -0.0132   0.0418   0.0140   1.0144    0.9869   1.0427
#sex        1.3965    -0.9665  -1.4741  -0.4589   0.2590   0.3804    0.2290   0.6320
#ph.karno  81.9383    -0.0405  -0.0613  -0.0196   0.0106   0.9604    0.9406   0.9806

#N = 227,  Events: 164,  Censored: 63
#Total time at risk: 69488
#Log-likelihood = -1132.418, df = 7
#AIC = 2278.836

#         33.33333% 66.66667%           
#1.609438  5.123929  5.825991  6.783325

# Get predicted survival at time=200 for individual cases using summary.flexsurv function
l$pred.30 <- summary(mod, newdata=l, type="survival", t=30, tidy=TRUE, B=0)$est

# Get predicted median survival for individual cases using summary.flexsurv function
l$pred.median <- summary(mod, newdata=l, type="median", tidy=TRUE, B=0)$est

I would like to provide the model formula to calculate this median survival for individual cases in my paper, but I cannot seem to figure out how to formulate this. For the prediction of survival at a certain time point the formula would be:

$log O (t ; x) = γ_0 + γ_1log(t) + γ_2ν_1(log(t)) + γ_3ν_2(log(t)) + β_1χ_1 + β_2χ_2 + β_3χ_3$

where $log O (t ; x)$ is the log (i.e. natural logarithm) cumulative odds survival function over time $t$ (in months) and $x$ represents the case covariates for the predictors $χ_1,χ_2,χ_3$. And where:

$ν_1 (z) = (z – k_j)_+^3 - λ_j (z-k_{min})_+^3 -(1-λ_j)(z – k_{max})_+^3$

$λ_j = (k_{max} - k_j)/(k_{max} - k_{min} )$

$γ_0$ = -3.2381 ; $γ_1$ = 0.9674 ; $γ_2$ = 0.1506 ; $γ_3$ = -0.3657

$k_{min}$ = 1.609438 ; $k_1$ = 5.123929 ; $k_2$ = 5.825991 ; $k_{max}$ = 6.783325

$β_1$ (age) = 0.0143 ; $β_2$ (sex) = -0.9665 ; $β_3$ (ph.karno) = -0.0405 ;

$k_{min}$ and $k_{max}$ represent the boundary knots of a natural cubic spline function. The other internal knots, $k_1$ and $k_2$, were placed at the 33% and 67% quantiles of the log uncensored survival times.

The absolute survival probability at $t$ months can then be given as $(1+ exp(log O (t ; x)))-1$

However, it is not clear to me how to provide a formula for calculating the median survival based on a Royston-Parmar model.

I know that in R the median survival for a Weibull model (another parametric survival model) can be formulated as follows:

# Fit Weibull model using psm formula (rms library)
mod.wb <- psm(Surv(time,status)~age+sex+ph.karno, data=l, x=T, y=T, dist="weibull")

# Calculate linear predictor for mod.wb
l$lp.mod.wb <- predict(mod.wb,type="lp",newdata=l)

# Calculate median survival for patients in l dataset 'manually'
l$pred.median.wb.M <- exp(l$lp.mod.wb + mod.wb$scale * log(-log(1-0.5)))

# Check to see whether predictions are correct by comparing them to those obtained 
# using the Quantile function (rms library)
Q.wb <- Quantile(mod.wb)
l$pred.median.wb.Q <- Q.wb(lp=l$lp.mod.wb)

# Test to see whether function and manual is the same
all.equal(l$pred.median.wb.Q, l$pred.median.wb.M) # TRUE

Given that this is possible for this parametric survival model, I would suspect that it is also possible to formulate this for a Royston-Parmar model (albeit probably more complex due to the restricted cubic splines). Does anyone know how to do this?

Any help would be greatly appreciated!


1 Answer 1


The way to get times for particular probabilities rather than probability for particular time you need the inverse of the survival function which is the quantile function. The flexsurv package also includes a qsurvspline function gives a quantile functional. So computing the survival time for a 50% survival should be fairly straightforward.


However that doesn't seem to take a model argument and apparently only handles the simplest of cases. Looking at the help page for flexsurv we see that there is a summary function that can derive quantiles,and it does accept your model argument:

summary(mod, type="quantile", quantiles=0.5)
  quantile      est      lcl      ucl
1      0.5 319.1645 280.6714 364.0389
  • $\begingroup$ Thanks. I do know how to get times for particular probabilities using the summary function (in this case the median using summary(mod, newdata=l, type="median", tidy=TRUE, B=0)$est as given in my example), but I would like to know how to write out the model formula for the median survival (that I can then provide in a paper so that others can use the model/externally validate it). So I am looking for a fomula for the Royston-Parmar model akin to the formula of the of the Weibull model (median survival = exp(l$lp.mod.wb + mod.wb$scale * log(-log(1-0.5))) as in the above example. $\endgroup$
    – Rob
    Commented Nov 25, 2020 at 21:57
  • $\begingroup$ You would have a separate predicted median survival time for each combination of covariates. If you wanted an analytic expression you would need to solve the inverse of the cubic spline survival formula. I’m guessing that ‘summary ( ... , type =“quantile”)’ computes the quantiles by constructing a relatively dense set of survival estimates and the swapping x (time) and y (survival). I’m pretty sure I’ve demonstrated that strategy over on SO for a different survival model. $\endgroup$
    – DWin
    Commented Nov 28, 2020 at 23:19
  • $\begingroup$ Thanks again. I contacted the author of the package flexsurv and to his knowledge there's no analytical formula for the median survival of a Royston-Parmar model. In the package it's computed by solving the equation CDF(x | parameters) = 0.5, for x, using numerical root-finding methods, since the CDF is known analytically. $\endgroup$
    – Rob
    Commented Nov 29, 2020 at 20:22
  • 1
    $\begingroup$ Pretty much as I said. Here's an example of the method:stackoverflow.com/questions/18813751/… done using predict but in this case you would be using version of summary provided for the package. If you wanted to see the code used type: class(mod) #flexsurvreg then getAnywhere(summary.flexsurvreg) $\endgroup$
    – DWin
    Commented Nov 29, 2020 at 22:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.