将分布拟合到生存曲线

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 的 Weibull 分布。但是我在理解我应该如何使用这个分布来计算我的估计生存率时遇到了问题。我有信心，因为 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% 的人有 died/failed）；我不知道 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")

将分布拟合到生存曲线

fitting a distribution to survival curve

statistics

r

curve-fitting

survival-analysis