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

2023年3月21日 1,319次阅读

我一直在尝试使用R(生存包和幸存函数)复制SAS生命测试的结果 – 特别是计算中位生存时间的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( – 日志(X))

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

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

以下是使用生存包中的肾脏数据的示例：

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代码中“随机”输入参数后,我设法解决了部分问题.

因此,幸存者使用上面给出的公式的对数变换来计算中值时间置信区间,这就是为什么R和SAS的间隔之间存在分歧(默认情况下使用对数 – 对数变换).

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

    `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`

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

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

点赞