# Flexible parametric survival (Royston-Parmar) model for non-proportional hazards

I am reading more about flexible parametric survival models and their implementation in the flexsurvspline function in the flexsurv package in R (there also seems to be a very similar command in Stata - stpm2)

I have a Cox model with several predictors that violate proportional hazards. I have tried transformations, time interactions and splitting time to generate time-dependent coefficients, without much luck for remedying the proportionality of the hazards.

If I use a model such as:

mod <- flexsurvspline(Surv(time, event) ~ x1 + x2 + x3 + x4, anc = list(gamma1 = ~ x1 + x2), data = dat, scale = "hazard")


am I correct in my understanding that the model will not force the hazards for x1 and x2 to be proportional (thus creating a time-varying HR for those two predictors)? This would also be a reasonable incidental way then to circumvent the problem of non-proportional hazards, right?

You can try the Royston-Parmar spline-modeling approach to accomplish what you want, but it takes some care. As I understand the implementation in the flexsurv package, you first specify the overall number of interior knots for the baseline hazard via the k argument, then specify which of the spline parameters gammaN (where N represents a spline-parameter number, with maximum k+1) are allowed to vary with covariate values.
Your example seems to use the default of 0 interior knots, which is just a Weibull model in which both the shape and scale parameters could vary with x1 and x2. See Section 5.1 and Table 2 of the vignette.
Given the difficulties you have had with the proportional hazards assumptions and these data so far, you thus might need to specify more knots overall and more corresponding spline coefficients gammaN to model as functions of the covariates x1 and x2. You also might want to consider the default placement of knots, depending on your particular data set.