Equation for flexible parametric survival model

[mattwarkentin]

Hi @chjackson,

I was wondering if there is a way to write out a closed-form expression for making a "prediction" with the result of a fit by flexsurvspline() for a single point in time.

For example,

library(flexsurv)
#> Loading required package: survival

model <- flexsurvspline(Surv(time, status==2) ~ age, data = cancer, k = 3)

tidy(model)
#> # A tibble: 6 x 5
#>   term   estimate std.error statistic p.value
#>   <chr>     <dbl>     <dbl>     <dbl>   <dbl>
#> 1 gamma0  -7.86     1.30        NA    NA     
#> 2 gamma1   0.963    0.311       NA    NA     
#> 3 gamma2  -0.239    0.311       NA    NA     
#> 4 gamma3   0.421    0.688       NA    NA     
#> 5 gamma4  -0.237    0.524       NA    NA     
#> 6 age      0.0179   0.00913      1.97  0.0247

Can the above info be used to predict survival at a fixed time point (e.g. 200 days) based on a single equation that could be included in a manuscript? I have played around with this to no luck so far.

[chjackson]

The vignette, section 5.1 gives the formula for the survival and cumulative hazard functions for this model. Unless I misunderstand what you want?

[mattwarkentin]

Thanks for the response. I guess what I am hoping to do is present an equation (or set of equations) that describe how to make survival predictions at any covariate values without access to the R model object, in the same way that one would present a linear model equation in a paper: Y^hat = 0.5 + (2 * X1) + ...

But for a flexsurvspline model at a fixed time point:

S(t=200) = exp(-H)

log(H) = s(x, gamma) + BZ

Getting the betas (B) for the predictors is simple. But how does one obtain the value for s(x, gamma) for a single time point to plug into the equation? I'm hoping this question makes some sense.

[chjackson]

Equation (4) at the bottom of page 15 in the vignette defines s() as a linear function of basis terms and the gamma parameters, and the first equation on page 16 defines these basis terms as polynomials. In flexsurv there's a function basis to compute these terms. You'll probablyalso need the knot locations, that are returned in the $knots component of the model object - these come from the quantiles of the uncensored survival times.

[mattwarkentin]

Thanks for pointing me in the right direction, got it sorted out!