Formula for predicting median survival using a Royston-Parmar model

Question

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")
mod

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

#Estimates: 
#          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

mod$knots
#         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!

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.

?qsurvspline

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