将 Weibull 分布拟合到...右删失数据？

机器算法验证 r 生存 审查 威布尔分布
2022-04-07 05:04:38

前提：我对生存分析几乎一无所知，我才刚刚开始。我有一些机器的故障时间向量。大多数机器（3362 vs 2694）今天仍在运行，所以我知道它们没有发生故障，并且故障时间在未来某个未知时间。如果我没记错的话，我们称之为右删失数据。正确的？现在，我想从简单的开始，为这些数据拟合 Weibull 分布。没有回归，没有协变量——只需对这些数据进行分布拟合，看看它是什么样子。这可能吗？将 Weibull 分布拟合到右删失数据的方法有哪些，R 中是否有这些方法可用？

编辑添加数据样本以显示我想做什么

# 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 适合分布，这正是我想要的。但我希望有一些函数来绘制生存函数，有置信度……看起来没有。

EDIT2对于那些对生存分析的简单介绍感兴趣的人，这个页面很短但很甜蜜： http ://data.princeton.edu/wws509/notes/c7s1.html

1个回答

拟合是在 R 的经典survival包（这个函数）中实现的，或者flexsurv，它具有更大的灵活性，但是 Weibull 的参数化不同（这里flexsurvreg的函数）。两者还提供了一些其他发行版供您试用。

虽然我不完全确定该方法，但看起来这两个包都使用似然最大化（通过检查源函数，例如body(flexsurv::flexsurvreg)）。

有关更详细的 R 说明和绘图示例，请参阅SO 上的这个问题。这里的想法是采用各种生存分数，用作以获得相应的时间，并自行绘制这些时间。在没有协变量的情况下，您可以调用只是，带有- 一些分位数。 S(t)predict.survregS−1(t)predict(weibull_fit, type="quantile", p=x)[1]x

请注意，这些只是给定估计参数的预测 - 如果您还想包含估计的不确定性，则需要切换到对数刻度，添加se.fit=T到predict调用中并使用它来计算预测时间的置信区间。就像是：

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)

应为预测的中位生存时间提供 95 % CI。

其它你可能感兴趣的问题
上一篇回归：量化的重尾自变量中的零 下一篇在哪些假设下，弱平稳性意味着强平稳性