Survfit()对区间删失数据生成的难以置信的宽置信区间

svujldwt  于 2023-04-18  发布在  其他
关注(0)|答案(2)|浏览(177)

我有一个数据是由间歇性访谈产生的,其中一个人被问到他们是否正在经历某种症状。每个人最后一次知道 * 没有 * 有这种特定的症状,表示为tstart。如果适用,观察到个人出现症状的时间是tstop。使用R中的survival包,使用Surv函数创建生存对象,指定这是区间删失数据。我想要生存函数的非参数最大似然估计。这可以使用survfit函数完成,这似乎将调用传递给了一个内部函数survfitTurnbull。由此产生的置信区间宽得令人难以置信。我无法弄清楚为什么会这样。

# A random sample of the data using dput()
structure(list(tstart = c(0.01, 38, 0.01, 0.01, 23, 26, 0.01, 
19, 0.01, 0.01, 22, 6, 0.01, 14, 16, 0.01, 0.01, 0.01, 0.01, 
21, 15, 0.01, 0.01, 13, 10, 0.01, 0.01, 19, 0.01, 0.01, 0.01, 
0.01, 22, 17, 27, 14, 16, 0.01, 20, 27, 10, 0.01, 0.01, 16, 20, 
7, 6, 15, 0.01, 0.01), tstop = c(4.01, NA, 5.01, 8.01, NA, NA, 
5.01, NA, 3.01, 16.01, NA, 6.01, 8.01, NA, NA, 7.01, 16.01, 1.01, 
10.01, NA, NA, 5.01, 8.01, NA, NA, 2.01, 3.01, NA, 7.01, 5.01, 
2.01, 9.01, NA, NA, NA, NA, NA, 10.01, NA, NA, NA, 5.01, 10.01, 
NA, NA, NA, 7.01, NA, 14.01, 4.01)), row.names = c(NA, -50L), class = "data.frame")

survObj <- with(temp_df, Surv(time = tstart, time2 = tstop, type = "interval2"))
survFit <- survfit(SurvObj ~ 1))
summary(survFit)

置信区间并没有随着时间的推移而变窄。使用整个数据集(包含大约10倍的事件数量),它并没有变窄。我无法弄清楚哪里出了问题。

nbysray5

nbysray51#

值得注意的是,这看起来并不像是软件中的错误,而是使用像非参数最大似然估计量这样灵活的东西的潜在限制(NPMLE,也称为Turnbull估计器,如果你给予它区间删失数据,那么survfit是拟合的),用于估计生存曲线。这个答案的TLDR版本是,我建议您使用参数模型,如Weibull,使用survival::survregicenReg::ic_paricenReg::ic_bayes。我是icenReg的作者。
关于NPMLE的一个有点技术性但非常相关的注意是,它仅将正概率质量分配给Turnbull间隔,Turnbull间隔是被定义为具有间隔的左侧是某个观察间隔的左侧并且Turnbull间隔的右侧是观察间隔中的 * 任何 * 的下一个最接近的右侧的间隔。我已经画出了你的观察间隔和相应的特恩布尔间隔。

请注意,在最后两个Turnbull间隔之间有一个巨大的差距!这导致了一个非常“跳跃”的NPMLE,这也导致了跳跃之间的相当多的误差。
在花了很长时间思考这个问题之后,我的总结是,这是只有少量信息数据和太多灵活性的结果。在大多数生存分析案例中,假设平滑的生存曲线是合理的,例如参数分布。只要分布不是过于严格(阅读:单参数指数分布),这种温和的平滑假设允许您从数据中获得更多的信息,而不会引入**太多的偏差。
为了说明,我附上了一个Weibull拟合+置信区间的图,旁边是拟合的NPMLE。

仅供参考,您看到的NPMLE框 * 不是 * 置信区间,而是NPMLE仅在分配给每个Turnbull区间的概率内是唯一的,但该概率在区间内的分布 * 不会影响对数似然。因此,任何通过该框的生存曲线都会使对数似然最大化。

x6492ojm

x6492ojm2#

我还发现默认的对于区间删失数据,Survfit的置信区间(NPLME)令人惊讶地宽,即使使用具有〉50000个观察值的数据集。(例如,观测次数多得多的地层有时比观测次数为十分之一的地层具有更宽的CI)。我不是数学家,但使用method = "rothman"method = "loglog"km.ci::km.ci产生更窄的置信限,以更符合我预期的方式响应样本量的变化。Sachs等人的这篇文章。(2022)似乎也很相关:https://www.nature.com/articles/s41416-022-01920-5

相关问题