0

我一直在准备 R 中的生存分析和 cox 回归。但是，我的直线经理是 Stata 用户，并且希望以与 Stata 显示它的方式类似的方式显示输出，例如

# Stata code
. strate
. stsum, by (GROUP)

stsum将为每个组输出一个风险时间和一个发生率，我不知道如何用 R 来实现这一点。

数据大致如下（我无法获取它，因为它处于安全环境中）：

PERS GROUP INJURY FOLLOWUP
111  1     0      2190
222  2     1      45
333  1     1      560
444  2     0      1200

到目前为止，我一直在使用相当混乱的标准代码：

library(survival)
library(coin)
# survival analysis
table(data$INJURY, data$GROUP)
survdiff(Surv(FOLLOWUP, INJURY)~GROUP, data=data)
surv_test(Surv(FOLLOWUP, INJURY)~factor(GROUP), data=data)
surv.all <- survfit(Surv(FOLLOWUP, INJURY)~GROUP, data=data)
print(sur.all, print.rmean=TRUE)
# cox regression
cox.all<- coxph(Surv(FOLLOWUP, INJURY)~GROUP, data=data))
summary(cox.all)
4

1 回答 1

1

目前我们有 4 行数据，但没有明确的描述（至少对于非 Stata 用户）所需的输出：

dat <- read.table(text="PERS GROUP INJURY FOLLOWUP
111  1     0      2190
222  2     1      45
333  1     1      560
444  2     0      1200",header=TRUE)

我不知道硬币或生存包中是否有功能可以为此类数据提供粗略的事件率。使用普通 R 函数提供粗略的事件率（在技术意义上使用“粗略”，无意贬低）是微不足道的：

 by(dat, dat$GROUP, function(d) sum(d$INJURY)/sum(d$FOLLOWUP) )
#----------------
dat$GROUP: 1
[1] 0.0003636364
------------------------------------------------------ 
dat$GROUP: 2
[1] 0.0008032129

风险时间的相应功能（或两者都打印到控制台）将是非常简单的修改。“Epi”或“epiR”包或其他专门用于教授基本流行病学的包之一可能为此设计了功能。“生存”和“硬币”的作者可能没有看到需要编写和记录如此简单的函数。

当我需要聚合因子协变量层内实际事件与预期事件的比率时，我需要构建一个函数来正确创建事件分层表（以支持置信度估计）、“预期”总和（根据年龄计算） 、性别和观察时间），并除以实际的 A/E 比。我将它们组合成一个列表对象，并将比率四舍五入到小数点后 2 位。当我完成它时，我发现这些最有用的方法是对我使用的“生存”和“有效值”回归方法得到的结果进行敏感性检查。它们还有助于向更熟悉表格方法而不是回归的非统计受众解释结果。我现在将它作为我的 Startup 的一部分.profile。

于 2015-02-07T17:47:20.997 回答