I am trying to understand the meaning of the coefficients estimates of the output of flexsurv's flexsurvreg function. For example, let us assume I want to perform the survival analysis and fit of a Weibull model with respect to a covariate. I will call flexsurvreg in order to obtain the Weibull parameters and covariate coefficient of the best fit flexsurvreg can obtain through its standard method.

fit <- flexsurvreg(Surv(time, censored)~covariate, data=struct, dist="weibull")

fit$res.t then returns the estimates of coefficients in such a fashion for example:

              est
shape         3
scale         4
covariate    -0.3

From there on, I want to try to reconstruct the analytic expression of the hazard function. From my understanding, it should be built in this fashion (for a Weibull model):

$h(t) = \frac{\text{shape}}{\text{scale}} \left(\frac{t}{\text{scale}}\right)^{\text{shape}-1}\cdot \exp\left(\text{coefficient} \cdot (\text{covariate}-\mu)\right)$

with $\text{coefficient}$ being the covariate coefficient returned in the flexsurvreg output; $\text{covariate}$ being the covariate value for which I intend to construct the function and $\mu$ the mean value of the covariate over the fitting sample.

However, the plot of this function does not match the one that plot(fit, type="hazard", newdata=list(covariate=cov)) returns. Why is that? I would guess my analytic interpretation of the coefficient estimates is wrong. What is then the mathematical meaning of these estimates?

Update: I have been through the package sources. It is still very unclear to me how the plot is done because I am not at all familiar with R's syntax nor the functions that are called.

However, I found out that I was wrong using fit$res.t, whereas fit$res actually is supposed to give out the parameters along the natural scale, as opposed to a logarithmic scale.

I also noticed that the computation of the analytical curve is likely done in the summary.flexsurvreg function.

Finally, the fit object contains functions that may be related to the computation of the hazard function in fit$dfns. How they take the covariate value into account still is fuzzy to me.