1
$\begingroup$

I've got following data representing survival function.

# A tibble: 53 x 2
   month survival
   <int>    <dbl>
 1     0    1.00 
 2     1    1.00 
 3     2    1.00 
 4     3    1.00 
 5     4    1.00 
 6     5    1.00 
 7     6    0.999
 8     7    0.998
 9     8    0.997
10     9    0.993
11    10    0.984
12    11    0.976
13    12    0.973
14    13    0.971
15    14    0.969
16    15    0.969
17    16    0.969
18    17    0.969
19    18    0.968
20    19    0.968
21    20    0.968
22    21    0.968
23    22    0.968
24    23    0.968
25    24    0.967
26    25    0.966
27    26    0.966
28    27    0.962
29    28    0.957
30    29    0.952
31    30    0.948
32    31    0.944
33    32    0.942
34    33    0.941
35    34    0.941
36    35    0.941
37    36    0.941
38    37    0.940
39    38    0.939
40    39    0.938
41    40    0.938
42    41    0.938
43    42    0.935
44    43    0.934
45    44    0.930
46    45    0.920
47    46    0.910
48    47    0.895
49    48    0.884
50    49    0.881
51    50    0.879
52    51    0.878
53    52    0.878

I would like to fit a distribution to a survival curve. To do so, firstly I plot the survival with respect to month. Then I use fitdist function in order to fit few distributions.

library('fitdistrplus')
library('flexsurv') 
data <- tibble(month = 0:52, survival = c(1, 1, 1, 1, 1, 1, 0.999, 0.998, 
0.997, 0.993, 0.984, 0.976, 0.973, 0.971, 0.969, 0.969, 0.969, 0.969, 0.968, 
0.968, 0.968, 0.968, 0.968, 0.968, 
0.967, 0.966, 0.966, 0.962, 0.957, 0.952, 0.948, 0.944, 
0.942, 0.941, 0.941, 0.941, 0.941, 0.940, 0.939, 0.938, 
0.938, 0.938, 0.935, 0.934, 0.930, 0.920, 0.910, 0.895, 
0.884, 0.881, 0.879, 0.878, 0.878))

data %>% ggplot(aes(month, survival)) + geom_line() 

fit_weibull <- fitdist(data[['survival']], 'weibull')
fit_llogis <- fitdist(data[['survival']], "llogis")
fit_log <- fitdist(data[['survival']], "logis")

fit_weibull$aic
fit_llogis$aic
fit_log$aic

According to AIC I should go for Weibull distribution with a shape = 34.6167936 and scale = 0.9695298. But I've got a problem with understanding how exactly should I use this distribution to calculate my estimated survival. I was confident that because S(t) = 1 - F(t) I should just calculate 1 -pweibull(data[['month']], fit_weibull$estimate[['shape']], fit_weibull$estimate[['scale']]), but it results in following vector:

 [1] 1.00000000 0.05399642 0.00000000 0.00000000 0.00000000 0.00000000 
 0.00000000 0.00000000
 [9] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
 0.00000000 0.00000000
 [17] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
 0.00000000 0.00000000
 [25] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
 0.00000000 0.00000000
 [33] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
 0.00000000 0.00000000
 [41] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 
 0.00000000 0.00000000
 [49] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000

So my understanding seems to be awfully wrong. How should I use fit_weibull to estimate a survival and plot the estimated curve then?

$\endgroup$
4
  • $\begingroup$ I noticed that a scatterplot of the data appears to show considerable seasonal variation. $\endgroup$ Jan 7, 2018 at 0:15
  • $\begingroup$ Indeed, but since these are survival data I would assume this is maybe due to censored information? $\endgroup$
    – jakes
    Jan 7, 2018 at 9:20
  • 2
    $\begingroup$ @jakes your analysis is wrong. 1) your data are the probability of survival at each time. However, this is not the kind of data fitdist takes as input - you need to pass it a vector with the time to event for each of the failure event, not the probabilities. This is the reason why your fit doesn't make any sense. 2) If there is censoring and you don't know when it happens, then you cannot perform a serious survival analysis. $\endgroup$
    – DeltaIV
    Jan 7, 2018 at 18:21
  • 2
    $\begingroup$ I'm voting to close this question as off-topic because it is cross posted on SO, where it already has an accepted answer: stackoverflow.com/questions/48137772/… $\endgroup$ Jan 10, 2018 at 19:43

1 Answer 1

5
$\begingroup$

The setting of your survival analysis is flawed: you have only survival probabilities at each time, but you don't have the time to event for each event, from which these probabilities were created. In other words, the raw data must have been something like this:

library(KMsurv)
data(tongue)
head(tongue)
# type time delta
# 1    1     1
# 1    3     1
# 1    3     1
# 1    4     1
# 1   10     1
# 1   13     1
  • a vector of times to failure or on-study time, in this case in weeks (column time
  • a status vector (here, column delta), indicating whether the recorded time is a time to failure (1) or is an on-study time (0, i.e., right-censored data: right-censored data are the most common, but we can also have left-censored data, interval data, truncated data, etc.)
  • (optionally) other covariates (in this case, the tumor DNA profile type

We start from these raw data, and we fit our distribution (Weibull, exponential, etc.) to it. fitdist is written to handle this user case, and it just doesn't work in your case where you start from the survival probabilities.

We can try another approach, but keep in mind that without the original time to failure data and the left/right/interval censoring tags, you won't obtain reliable results. We start from the observation that the survival function of the Weibull distribution is

$$p = P(T > t) = \begin{cases} 0 & t < 0 \\ -e^{(t/\lambda)^k} & t \geq 0 \end{cases} $$

We want to estimate the scale parameter $\lambda$ and the shape parameter $k$. Take logs of LHS and RHS so that we can use standard Least Squares estimation:

$$-\log{p} = \left(\frac{t}{\lambda}\right)^k $$

Still not linear in the parameters, thus let's take logs again:

$$\log({-\log{p}}) = k\left(\log{t}-\log{\lambda}\right) \implies \frac{1}{k}\log({-\log{p}}) = \log{t}-\log{\lambda} \implies \log{t} = \log{\lambda} + \frac{1}{k}\log({-\log{p}})$$

which is a linear model in the parameters $\log{\lambda}$ and $\frac{1}{k}$. Let's estimate it: for each statistical analysis, the first step is to plot your data and interpret the plot.

library(ggplot2)
library(dplyr)
library(fitdistrplus)
library(flexsurv)

data %>% ggplot(aes(month, survival)) + geom_line()

enter image description here

You see a staircase behavior which of course cannot be matched by a Weibull (or an exponential, logistic, etc.) distribution. This brings us back to the point that you will only get reliable results if you recover the original data together with the censoring/truncation information.

Another problem, is that for the Weibull model, $S(t)=1 \implies t = 0$. You cannot get a 100% probability of survival for all months 0 to month 5. In practice, this means that for those data points, either $\log{t}$ or $\log({-\log{p}})$ are undefined. I solved this by subtracting a very small probability of survival until month 5, but this is just an hack and you should really get those raw data.

mod_data <- data[-1, ]
mod_data <- mutate(mod_data, 
                   p = ifelse(month <= 5, 1 - 0.001/6*month, survival), 
                   log_month = log(month), 
                   loglog_p = log(-log(p)))

mod_data %>% ggplot(aes(month, survival)) +
  geom_point() +
  geom_point(aes(y = p), color = "red")

enter image description here

We can now fit our model:

weibull_model <- lm(log_month ~ loglog_p, data = mod_data) 
summary(weibull_model)
# 
# Call:
#   lm(formula = log_month ~ loglog_p, data = mod_data)
# 
# Residuals:
#   Min       1Q   Median       3Q      Max 
# -0.58000 -0.09654  0.05872  0.15372  0.46430 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  4.90148    0.08388   58.43   <2e-16 ***
#   loglog_p     0.51743    0.02091   24.75   <2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.2474 on 50 degrees of freedom
# Multiple R-squared:  0.9245,  Adjusted R-squared:  0.923 
# F-statistic: 612.6 on 1 and 50 DF,  p-value: < 2.2e-16

Let's make a Weibull plot, which is just the plot of our Weibull fit and the original data on the transformed scale:

log_lambda <- weibull_model$coefficients[1]
k <-  1/weibull_model$coefficients[2]

mod_data <- mutate(mod_data, loglog_weibull = k * (log_month - log_lambda))

ggplot(mod_data, aes(x = log_month, y = loglog_p)) +
         geom_point() +
         geom_line(aes(y = loglog_weibull), color = "red") +
         labs(x = expression(log(month)),
              y = expression(log(-log(p))))

enter image description here

We already see the limitations of our model, most likely due to the lack of censoring information. Let's also plot the fit on the original scale, using the original data: we can do this because on the original scale we don't have to worry about $\log{t}$ or $\log({-\log{p}})$ being undefined.

lambda <- exp(log_lambda)
mod_data <- mutate(mod_data, weibull = exp(-(month/lambda)^k))

ggplot(mod_data, aes(x = month, y = survival)) +
  geom_point() +
  geom_line(aes(y = weibull), color = "red") +
  labs(x = "month",
       y = "survival")

enter image description here

Here the issues with this statistical model are even more evident. The conclusion is clear: go get those raw data!

$\endgroup$

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