# Predicting survival/event probability with multi-level Weibull model and time-dependent covariates

I am still a beginner when it comes to survival analysis. I have fitted a parametric (Weibull) survival regression model with time-dependent covariates using the R package flexsurv via:

fitParam <- flexsurvreg(Surv(tStart, tStop, event) ~ shape(a) + b + cluster(subject),
data = data,
dist = "weibull")


I included a random effect to account for repeated trials by subjects (I assume the cluster() in the formula to take care of that).

When plotting the fitted model object via plot(fitParam), it seems to match decently well with the data:

Now, I would like to predict the event probability (1 - survival probability) and uncertainty for new data (a = 10, b = -1.0, subject = 10), sampled, for instance, at time time <- 0.5. I tried to do it via predict:

time <- 0.5 # Example time
pred <- predict(fitParam,
newdata = data.frame(a = 10, b = -1.0, subject = 10),
type = "survival",
times = c(time),
conf.int = T)


This gives the output:

# A tibble: 1 x 4
.time .pred .pred_lower .pred_upper
<dbl> <dbl>       <dbl>       <dbl>
1   0.5 0.900       0.656       0.985


My questions:

1. Is this a valid way to model and predict event probability, given time-varying covariates?
2. How could I predict for an unknown subject, i.e., if I cannot supply a subject to predict?

Edit-1:

I actually tried fitting a Weibull model only using the shape(a) part, without b in the formula. The resulting model seems to visually fit the data better (see below), but has a higher AIC than the model shown above.

I am trying to fit a model that can predict the probability of an event at any time, given covariates that vary with time. I have repetitions from individuals. My previous question previous question may describe this in more details.

• Where fitParam comes from isn't clear, but there doesn't seem to be a very good match of predicted (red) and actual (black) survival. The predicted curve is close to the observed upper confidence interval (CI) for much of the time span; the observed curve is along the predicted lower CI for most times. And although allowing for a covariate effect on shape is allowed, that's getting beyond standard accelerated failure time modeling. Please say more about what you are trying to accomplish with your model; there might be a better approach that you are missing as "a beginner."
– EdM
Aug 18, 2021 at 15:02
• Thanks. I edited the question with more details. fitParam is a flexsurvreg object, as shown in the first code block. Aug 20, 2021 at 21:50

Trying to make predictions from parametric models with time-varying covariates is not the easiest way for a beginner to start with survival analysis. A few thoughts.

First, think about why you are modeling.

A) If your main interest is in prediction, there might be no need to force your data into a particular standard parametric functional form like the Weibull. In particular, a Weibull model assumes proportional hazards (PH), so a more general PH model, a semi-parametric Cox model or a parametric PH model with a flexible baseline survival, might work better. Alternatively, you could try a non-PH accelerated failure time (AFT) model; only the Weibull family (with its limiting exponential case) meets both PH and AFT. See these course notes for a succinct introduction to parametric survival modeling.

B) If your main interest is interpretable inference, try to use a model in which the coefficients have reasonable interpretations. In your Weibull model, using the parameterization favored by the flexsurv package, covariates modeled as part of the scale parameter (the default, like for b in your first model) have a nice interpretation in terms of speeding up or slowing down the passage of time. You can model covariates as part of the shape parameter (like shape(a) in your first model), but interpretation is not so straightforward. With an AFT model, the shape parameter adjusts the width of the associated error distribution of log-survival times around what you would get from just the linear predictor;* a fixed shape is easy to understand in that context. I have a hard time thinking about what allowing covariates to affect the shape would mean in an AFT context.

Second, if only one event of one type per individual is possible then you don't need to take their id values into account, as explained in the main survival package time-dependence vignette on page 4. If you do have multiple events possible per subject, then you apparently shouldn't use flexsurvreg(). The main flexsurv vignette says (page 3): "The individual survival times are also [assumed] independent, so that flexsurv does not currently support shared frailty, clustered or random effects models."

Third, things that work in one software package might not in another. You might have avoided the confusion with differing Weibull parameterizations, but you got caught in your use of a cluster() term in flexsurvreg(). Although it's not supposed to deal with clustered data, in version 2.0 it seems to accept them in a strange way. In other R survival software, cluster() terms affect the coefficient covariance matrix, leaving coefficient point estimates unaltered. When I added a cluster() term to a flexsurvreg() model, I got a coefficient estimate for cluster(id) and altered values for the other coefficients! Not sure how to interpret that. You should omit your cluster() term if you use flexsurv.(That also solves one of your problems.)

Fourth, a single plot or a simple test value like AIC isn't sufficient to validate a survival model. For a parametric model, you need to verify that you chose the correct general form; for a PH model, you need to document how well the PH assumption holds. You need to check that you have included an appropriate set of predictors, possibly including interactions among them. For continuous predictors you need to document that you used an appropriate functional form. For the main survival package there are several types of "residuals" available to help with such checks when event times are censored, but (as you probably know) the standard survreg() function doesn't accept the counting-process data format used for time-varying covariates.

Harrell's course notes and book show ways to validate both PH and AFT models. To apply those techniques to parametric models from flexsurv or the eha package, which can accept counting-process data, you will have to do some coding. The only "residual" reported directly at this time by either package seems to be the "response" residual in version 2.0 of flexsurv, the difference between predicted* and observed event times. That type of residual can be misleading if used improperly. More useful for many AFT models are the scaled residuals,* if you also keep information about censoring versus events at the observation times. That information is in the models, but for now you'll have to pull it out yourself.

Fifth, be very very careful about trying to make predictions involving time-varying covariates. Although your model was built from data including time-varying covariates, you only specified a single constant set of values for a and b that are assumed to hold from time = 0 to time = 0.5 in your example. That's OK. But if you start specifying sets of time-varying covariates for predictions you can get into serious logical problems. That's allowed by coxph() models in R, but the author of the Python lifelines package made a deliberate choice not to allow such predictions. It's easy for predictions based on time-dependent covariates to suffer from survivorship bias.

*AFT models can be written as

$$\log T = X' \beta + \sigma W$$

where $$X' \beta$$ is the linear predictor based on covariates $$X$$ and coefficients $$\beta$$, $$W$$ is a specific probability distribution associated with the type of AFT model, and $$\sigma$$ is a shape parameter. An estimated scaled residual is $$(\log T - X' \hat\beta)/\hat\sigma$$, which should have the same (censored) distribution as $$W$$.

If $$W$$ is a symmetric distribution, then with $$W=0$$ the linear predictor estimates mean $$\log T$$. If $$W$$ isn't symmetric--like with the minimum extreme value distribution assumed by a Weibull model--then you have to be careful just what is being "predicted" by a software package. It isn't necessarily the mean (log) survival time. Or it might be. Read the manual carefully.

• Thank you very much for the elaborated answer and sorry for the late follow-up. Many thanks also for the helpful references. I am actually interested in both, A) prediction and B) interpretable inference. But say I am more interested in A) prediction, is there some metric that I could evaluate/compare different models (AFT, PH) on? Good catch for the flexsurv package, it indeed does not seem to support random effects and using the cluster() term seems wrong as you say. Do you know any other package/method where I could include random effects (repeated trials)? Oct 11, 2021 at 13:04
• @Tester01 can a single individual experience more than 1 event in your study? Is there any other form of correlation other than within a single individual (e.g., multiple individuals treated at a single hospital)? Just what do you mean by "repeated trials"?
– EdM
Oct 11, 2021 at 13:17
• Yes, one single individual experienced multiple events in the study. The data are from an experiment with a repeated measures design, that is why I referred to "repeated trials." I am not aware of any other (obvious) form of correlation. Oct 11, 2021 at 14:22
• @Tester01 coxph() fits Cox models to such data, correcting for repeated events in individuals with an id (random effect-like) or a cluster (sandwich variance) argument to the function. aftreg() in the eha package allows for both time-varying coefficients and id terms in parametric models. See the Fourth section of the answer for evaluating/comparing models. There are many tools for evaluating coxph() models; for aftreg() you might have to do some coding. Comparing model types is tricky. Think hard about your comparison criteria.
– EdM
Oct 11, 2021 at 14:48
• Thanks for the suggestions, I will try coxph() and aftreg() with the id parameter. Oct 14, 2021 at 8:40