Extracting transition-specific estimates from a parametric multi-state survival model

Question

I fit a multi-state survival model following Weibull distribution using function flexsurvreg. This data is about the depression trajectory of people with mild depression. For a follow-up economic evaluaiton, I need to extract two types of the data: 1) transition-specific probabilities to be used in a Markov model, 2) transition-specific shape and scale parameters to be used in a Discrete Event Simulation model to sample time-to-event using rweibull function. In another words, I need to extract parameters that would allow me subsequently sample time to each specific transition (i.e. mild-moderate, mild-severe, mild-remission) from the distribution.

I manage to get the probabilities using function pmatrix.fs (task 1). However, I struggle to extract transition-specific parameters (task 2). Right now, I obtain only one coefficient for "trans" variable regardless to where the person transited from the starting state. One solution would be to stratify the data according to the end state, however, this would not account for a competing risk that is present in the data.

Some simplified code that I use:

Survival analysis

crwei <- flexsurvreg(Surv(time, status) ~ trans + shape (trans), data = data.mild,
                          dist = "weibull")

Extracting transition probabilities at 26 weeks

tmat <- trans.comprisk(3, names = c("mild", "remission", "moderate", "severe")
probs_wei      <- pmatrix.fs(crwei,  t = 26, trans = tmat)

Below is the matrix with transition probabilities (column 1 represents staying in mild depression, column 2 tranisiting mild-remission, 3 mild-moderate, 4 mild-severe).

Unfortunately, I do not know how to extract from the crwei object transition-specific parameters for the Weibull distribution as the coefficients do not seems to be transition-specific (see output picture). My guess is that it shall be possible, because otherwise I would not be able to obtain transition-specific probabilities and I also found that it is possible to extract also transition-specific hazard using function msfit.flexsurvreg (see this article by Jackson). An advise on how to extract transition-specific coefficients is appreciated.

Answer 1

With this type of data, you probably need trans to be coded as a factor to provide different coefficients for each transition (trans). I suspect that you have coded it as numeric, so that flexsurvreg() interpreted it as a continuous predictor and thus only provided a single coefficient value for each of trans and shape(trans). Try using factor(trans), or give the levels of trans interpretable names if you wish to proceed with a parametric survival model.

You also should seriously consider Frank Harrell's suggestion to handle this instead with an ordinal longitudinal model.