将分布拟合到生存曲线

Question

我有以下代表生存函数的数据。

# 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

我想将分布拟合到生存曲线。为此，我首先绘制相对于月份的生存率。然后我使用fitdist函数来适应一些分布。

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

根据 AIC，我应该选择带有 和 的威布尔shape = 34.6167936分布scale = 0.9695298。但我在理解应该如何使用这个分布来计算我的估计生存期时遇到了问题。我很有信心，因为S(t) = 1 - F(t)我应该只计算1 -pweibull(data[['month']], fit_weibull$estimate[['shape']], fit_weibull$estimate[['scale']])，但它会产生以下向量：

 [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

所以我的理解似乎是非常错误的。那么我应该如何fit_weibull估计生存并绘制估计曲线呢？

Answer 1

您需要在这里处理非标准版本的生存分析。通常生存分析数据是根据离散事件（个体死亡的时间）来记录的 - 这就是flexsurv（您加载的但据我所知没有使用）所期望的。

不幸的是fitdistrplus::fitdist，也不适用于您的数据 - 这将期望生存时间的分布。此外，即使您确实有独立生存时间的数据，您的数据也会受到审查（在该时间段结束时只有 12% 的个体死亡/失败）；我不知道是否fitdist允许审查。

您可能无法对曲线之间的差异做出非常有力的统计结论，因为您不知道（或者至少您没有说过）该生存曲线实际上代表了多少独立试验 - 例如是初始试验由 10、100 或 10^6 个人组成的队列...？

但是，您可以按如下方式拟合曲线：

dat <- data.frame(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))

通过非线性最小二乘拟合（不是一个很好的统计模型，但足够了）。另外：需要良好的起始值。

n1 <- nls(survival~pweibull(month,exp(logshape),exp(logscale),
                      lower.tail=FALSE),
    start=list(logshape=0,logscale=log(20)),data=dat)
n2 <- nls(pmin(survival,0.999)~plogis(month,location,exp(logscale),
                            lower.tail=FALSE),
          start=list(location=40,logscale=log(20)),data=dat)

情节结果：

par(bty="l",las=1)
plot(survival~month,data=dat,type="l")
lines(dat$month,predict(n1),col="red")
lines(dat$month,predict(n2),col="blue")