Solved – Fit a Weibull distribution to…right-censored data

censoringrsurvivalweibull distribution

Premise: I know next to nothing about survival analysis, I just started. I have a vector of failure times for some machines. Most of the machines (3362 vs 2694) are still running today, so I know they haven't failed, and the failure time is at some unknown time in the future. If I'm not mistaken, we call this right-censored data. Correct? Now, I would like to start simple and fit a Weibull distribution to these data. No regression, no covariates – just fit a distribution to these data, and see what it looks like. Is this possible? What are the methods used to fit a Weibull distribution to right-censored data, and are there some of these methods available in R?

EDIT adding a data sample to show what I would like to do

# sample data set
foo <- structure(list(time = c(416L, 307L, 552L, 740L, 1126L, 442L, 
765L, 375L, 275L, 623L, 455L, 695L, 724L, 1653L, 529L, 615L, 
15L, 477L, 758L, 353L, 851L, 25L, 524L, 283L, 954L, 183L, 558L, 
851L, 377L, 208L), event = c(1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 
0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1)), .Names = c("time", 
"event"), row.names = c(1170L, 5856L, 1237L, 1448L, 1161L, 4759L, 
2766L, 3321L, 4253L, 5645L, 3875L, 3531L, 238L, 3204L, 714L, 
2161L, 3275L, 5879L, 3758L, 2406L, 1548L, 4234L, 87L, 2548L, 
3958L, 5872L, 3326L, 1549L, 4823L, 1002L), class = "data.frame")

# fit a Weibull distribution
weibull_fit <- survreg(Surv(time, event) ~ 1, data = foo)
summary(weibull_fit)

# Call:
# survreg(formula = Surv(time, event) ~ 1, data = foo)
#             Value Std. Error     z        p
# (Intercept)  6.881      0.132 52.14 0.000000
# Log(scale)  -0.758      0.197 -3.86 0.000115
#
# Scale= 0.469 
# 
# Weibull distribution
# Loglik(model)= -101.2   Loglik(intercept only)= -101.2
# Number of Newton-Raphson Iterations: 8 
# n= 30

R fits the distribution, which is what I wanted. But I hoped there would be some function to plot the survival function, with confidence bounds…it looks like there's none.

EDIT2 for those interested in a simple intro to Survival Analysis, this page is short but sweet:
http://data.princeton.edu/wws509/notes/c7s1.html

Best Answer

Fitting is implemented in R's classic survival package (this function), or flexsurv, which has more flexibility, but a different parametrization for Weibull (flexsurvreg function here). Both also provide a few other distributions to try out.

Although I am not entirely sure about the method, it looks like both packages use likelihood maximization (from inspection of source functions, e.g. body(flexsurv::flexsurvreg)).

For more detailed R instructions and plotting examples, see this question on SO. The idea there is to take various survival fractions $S(t)$, use predict.survreg as $S^{-1}(t)$ to get the corresponding times, and plot those on your own. In the no-covariate case that you have the call can be just predict(weibull_fit, type="quantile", p=x)[1], with x - some quantile.

Note that these are just the predictions given estimated parameters - if you want to incorporate the uncertainty about the estimates as well, you will need to switch to log scale, add se.fit=T to the predict call and use that to calculate the confidence intervals for the predicted times. Something like:

pr = predict(weibull_fit, type="uquantile", p=0.5, se.fit=T)
lims = c(pr$fit[1] + 1.96*pr$se.fit[1], pr$fit[1] - 1.96*pr$se.fit[1])
lims = exp(lims)

should give 95 % CI for predicted median survival time.

Related Solutions

Solved – Right-censored survival fit with JAGS

I was asked to re-post this answer here from my comment at http://doingbayesiandataanalysis.blogspot.com/2012/01/complete-example-of-right-censoring-in.html The specifics of this answer relate to the model in that comment, but the concepts apply to the topic here.

The core of the JAGS model for censored data is this:

isCensored[i] ~ dinterval( y[i] , censorLimitVec[i] )
y[i] ~ dnorm( mu , tau )

The key to understanding what JAGS is doing is that JAGS automatically imputes a random value for any variable that is not specified as a constant in the data. Thus, when y[i] is NA (i.e., a missing value, not a constant), then JAGS imputes a random value for it.

But what value should it generate?

The second line of the model, above, says that y[i] should be randomly generated from a normal distribution with mean mu and precision tau.

But the first line of the model, above, puts another constraint on the randomly generated value of y[i]. That line says that whatever value of y[i] is randomly generated, it must fall on the side of censorLimitVec[i] dictated by the value of isCensored[i].

To understand this part, let's unpack the dinterval() distribution. Suppose that censorLimitVec has 3 values in it, not just 1:

censorLimitVec = c(10,20,30)

Then randomly generated values from dinterval(y,c(10,20,30)) will be either 0, 1, 2, or 3 depending on whether $y<10$, $10 < y < 20$, $20<y<30$, or $30<y$. So, if $y=15$, dinterval(y,c(10,20,30)) has output of $1$ with 100% probability. The trick is this: We instead specify the output of dinterval, and impute a random value of y that could produce it. Thus, if we say

1 ~ dinterval(y,c(10,20,30))

then y is imputed as a random value between 10 and 20.

Putting the two model statements together,

1 ~ dinterval( y , censorLimit )

y ~ dnorm( mu , tau )

means that y comes from a normal density and y must fall above the censorLimit.

Hope that helps!!

R – Relationship between Gumbel and Weibull Distributions and Survival Analysis

The confusion comes from competing definitions of "Gumbel distribution" and competing parameterizations of the Weibull distribution.

(1) It might be best to avoid the term "Gumbel distribution" because it has different interpretations.

One is a maximum extreme value distribution, the definition used in Wikipedia. "This article uses the Gumbel distribution to model the distribution of the maximum value." (Emphasis in original.)

Another is a minimum extreme value distribution, the definition provided by Wolfram. "In this work, the term 'Gumbel distribution' is used to refer to the distribution corresponding to a minimum extreme value distribution." (Emphasis added.) That is used by Mathematica for its GumbelDistribution, which calls the Wikipedia maximum extreme value version the ExtremeValueDistribution.

It's the minimum extreme value version that provides the "standard result" for the association between Weibull and Gumbel distributions. As you used the maximum extreme value version, you got the result that you found.

(2) Continuing from point (1), to make this work you have to alter (a) the relationship between $\alpha$ and $\beta$ to get a mean of 0, and (b) the CDF to match the minimum extreme value Gumbel.

(a) The mean of the minimum extreme value version is $\alpha - \gamma \beta$, where $\gamma$ is Euler's gamma, with $\alpha$ and $\beta$ as represented in the question. That's different from $\alpha + \gamma \beta$ for the maximum extreme value version, as used in the question.

(b) The $q$th quantile (inverse CDF) of the minimum extreme value version is:

$$\alpha +\beta \log (-\log (1-q)).$$

The inverse CDF used in the question's code is for the maximum extreme value version.

I haven't yet done those replacements in the code, but I suspect that (absent other problems) all with then be OK.

(3) The question of "exactly what distribution specification R is using when fitting a Weibull distribution" is not well specified.

R packages can differ in parameterizations, and the same function might use different parameterizations depending on the arguments in the function call. This page provides some examples. Notably, as the manual page for the survreg() function in the survival package explains:

There are multiple ways to parameterize a Weibull distribution. The survreg function embeds it in a general location-scale family, which is a different parameterization than the rweibull function, and often leads to confusion.

survreg's scale  =    1/(rweibull shape)
survreg's intercept = log(rweibull scale)

I don't see any way around these types of confusions, except to be extremely careful in reading specific definitions and manual pages.

Best Answer

Related Solutions

Solved – Right-censored survival fit with JAGS

R – Relationship between Gumbel and Weibull Distributions and Survival Analysis

Related Question