I have used flexsurvreg to estimate the parameters of a Weibull distribution and got the following output. I would like to reconstruct the survival function to estimate the survival rate at a given time t.

flexsurvreg(formula = Surv(time, Censor) ~ 1, data = CC, dist = "weibull")

shape 1.02
scale 2630.00

I have read the flexsurvreg document but it is still not very clear to me. Could anyone please point me in the right direction?

My plan is then to use Excel weibull.dist function to estimate the survival rate at a given point in time.


Weibull distribution has probability density function

$$ f(x;\lambda,k) = \begin{cases} \frac{k}{\lambda}\left(\frac{x}{\lambda}\right)^{k-1}e^{-(x/\lambda)^{k}} & x\geq0 ,\\ 0 & x<0 \end{cases} $$

where $\lambda>0$ is scale parameter and $k>0$ is shape parameter. Different values of parameters are presented on the plot below.

Parameters of Weibull distribution

Basically, as the names suggests, shape parameter controls it's shape and scale parameter makes it wider or narrower (notice the $x/\lambda$ parts of probability density function). To gain more intuition you can plot yourself the distribution with different parameter values and check what happens when you change them.

The description of shape parameter is nicely summarized on Wikipedia:

If the quantity $X$ is a "time-to-failure", the Weibull distribution gives a distribution for which the failure rate is proportional to a power of time. The shape parameter, $k$, is that power plus one, and so this parameter can be interpreted directly as follows:

  • A value of $k < 1$ indicates that the failure rate decreases over time. This happens if there is significant "infant mortality", or defective items failing early and the failure rate decreasing over time as the defective items are weeded out of the population. In the context of the diffusion of innovations, this means negative word of mouth: the hazard function is a monotonically decreasing function of the proportion of adopters;
  • A value of $k = 1$ indicates that the failure rate is constant over time. This might suggest random external events are causing mortality, or failure. The Weibull distribution reduces to an exponential distribution;
  • A value of $k > 1$ indicates that the failure rate increases with time. This happens if there is an "aging" process, or parts that are more likely to fail as time goes on. In the context of the diffusion of innovations, this means positive word of mouth: the hazard function is a monotonically increasing function of the proportion of adopters. The function is first concave, then convex with an inflexion point at $(e^{1/k} - 1)/e^{1/k}, k > 1$.
| cite | improve this answer | |
  • $\begingroup$ +1. And since the OP's k is roughly 1, this says that they might be able to simplify their model to an exponential. $\endgroup$ – Wayne Feb 2 '17 at 20:59
  • $\begingroup$ @Tim thanks for your reply. My apologies for not being clear about what my confusion came from. I knew the basic of the Weibull distribution but just was not sure if we can just take the estimated scale and shape as they were without needing to do log transformation. I know in SAS the shape is actually 1/scale (or something like that). But looks like with R, I do not need to do any further transformation before using the parameters. $\endgroup$ – C.Woo Feb 2 '17 at 21:52
  • $\begingroup$ @Wayne Great idea. it turns out that exponential distribution works just fine. $\endgroup$ – C.Woo Feb 2 '17 at 21:55
  • 1
    $\begingroup$ @C.Woo: Depending on the distribution you're using, you may have something like rate = 1/scale and need to specify shape= and scale= rather than specifying by-position so you don't accidentally specify the scale as the rate. Note also that if you use coef on your flexsurvreg result, you may see results that need to be exp'd or something. (I've been working with a Burr distribution.) $\endgroup$ – Wayne Feb 2 '17 at 22:22

dist = "weibullPH" rather than dist = "weibull" provides scale and shape parameters that can be used as λ and γ, respectively in the parameterisation below to calculate the survival function. (k AKA γ).

Using larynx data:

larynx$late_stage <-  factor(larynx$stage == 3 | larynx$stage == 4) #create a binary variable for illustration

swei  <- flexsurvreg(s ~ late_stage,dist='weibullPH', data=larynx)
plot(swei, type = "survival", xlim = c(0,10)) #plot the output from the WeibullPH model

#Now extract the lambda and gamma from swei$res

lambda <- swei$res[2]
gamma <- swei$res[1]
t<- larynx$time

You can then plot the calculated survival function for group 1 using $ S(t) = exp(-\lambda t^{\gamma}) $ and for group 2 by introducing the effect of the coefficient through raising the previous equation to the power of the linear predictor $ S(t) = exp(-\lambda t^{\gamma})^{exp(\hat{\beta}x_{i}) } $

group1 <- exp(-lambda*(t^gamma))
group2 <- group1^(exp(swei$coefficients[3]))

plot(swei, type = "survival", xlim = c(0,10))
plot(t, group1 , ylim=c(0,1), xlim = c(0,10))
plot(t, group2, ylim=c(0,1), xlim = c(0,10))

Weibull plot

You can see that your manual survival plot overlies that of the function perfectly.

You could achieve the same thing using using the scale and shape values provided with dist="Weibull" using the alternative parameterisation $ S(t) = exp(-(\frac{t}{\lambda})^{\gamma }) $

| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.