使用R survival package计算中位生存时间的proc lifetest 95％CI

Posted 2021-04-20

我一直在尝试使用R（proc lifetest包和survival函数）在SAS中复制survifit的结果 - 特别是计算中位存活时间的95％置信区间。

我知道SAS正在使用以下公式来计算中位数的置信区间：

*abs(g(S(t))-g(1-0.5)/g'(S(t))σ(S(t)))<=1.96*

g'（x）是g（x）的一阶导数，σ（S（t））是生存曲线的标准误差，SAS中g的默认变换是g(x)=log(-log(x))

所以绝对内部的公式变为：

(log(-log(S(t)))-log(-log(0.5)))*S(t)*log(S(t))/σ(S(t))

以下是使用kidney包中的survival数据的示例：

fit1 = survfit(Surv(kidney$time,kidney$status)~kidney$sex, data=kidney)
print(fit1)
BCinds<-abs((log(-log(fit1$surv))-log(-log(0.5)))*fit1$surv*log(fit1$surv)/fit1$std.err)<=1.96

当我运行我从print(fit1)获得的代码时：

                n events median 0.95LCL 0.95UCL
kidney$sex=1 20     18     22      12      63
kidney$sex=2 56     40    130      66     190

然而，当我通过BCinds计算它时，我得到一个非常不同且更宽的CI（9,154），性别= 1，而性别= 2，CI是（39,511）。

sex=1 95%CI: (9, 154)  sex=2 95%CI: (39, 511)

SAS代码还为同一数据集的中位生存时间生成不同的置信区间：

    ods graphics on;
proc lifetest data=work.test
    plots=survival(nocensor cb=hw cl strata=panel);
    strata sex/group=sex;
    time time*status(0);
    run;
ods graphics off;

结果如下：

 sex=1: median=22 and 95%CI: (12, 30)
 sex=2: median=130 and 95%CI: (58,185)

为什么我会得到如此不同的结果的任何想法？你也可以建议我如何自动化方法的最后一步？目前我是在视觉上做的，但我想把它放在一个循环中，所以我需要自动完成。

谢谢！

答案

更新

因此，在R代码中“随机”输入参数后，我设法解决了部分问题。

所以survfit使用上面给出的公式的对数变换来计算中值时间置信区间，这就是为什么R和SAS的区间之间存在分歧（默认情况下使用对数 - 对数变换）。

因此，通过在R代码中添加一个参数，我们可以强制R以与SAS相同的方式计算置信区间。因此，对于上面我给出的kidney数据的示例，我们有：

    `survfit(Surv(kidney$time,kidney$status)~kidney$sex, conf.type="log-log"
    + )
    Call: survfit(formula = Surv(kidney$time, kidney$status) ~ kidney$sex, 
        conf.type = "log-log")

              n events median 0.95LCL 0.95UCL
kidney$sex=1 20     18     22      12      30
kidney$sex=2 56     40    130      58     185`

我们可以从survfit获得的其他置信区间类型是：“log”, “log-log”, “plain”, “none”

我仍然没有弄清楚我用来获得置信区间的代码有什么问题，所以如果有人知道它有什么问题我会很感激任何反馈。

另一答案

我想这是因为fit1$std.err中的BCinds部分。在这里你应该适应S(t)的标准误差 - 但是fit1$std.err（根据survfit.object的R文档）给出了累积危险或-log（生存）的标准误差。尝试使用summary(fit1)$std.err代替。