每组的 R 风险时间
R time at risk for each group
我一直在用 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 行数据并且没有对所需输出的明确描述（至少对于非 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 函数提供粗略的事件率（在技术意义上使用 'crude'，无意贬低）是微不足道的：
by(dat, dat$GROUP, function(d) sum(d$INJURY)/sum(d$FOLLOWUP) )
#----------------
dat$GROUP: 1
[1] 0.0003636364
------------------------------------------------------
dat$GROUP: 2
[1] 0.0008032129
对应的time at risk函数（或者都打印到控制台）将是一个非常简单的修改。 'Epi' 或 'epiR' 包或其他专门用于教授基本流行病学的包之一可能会为此设计功能。 'survival' 和 'coin' 作者可能认为没有必要编写和记录如此简单的函数。
当我需要在因子协变量层内汇总实际事件与预期事件的比率时，我需要构建一个函数来正确创建事件分层表（以支持置信度估计），总和 "expecteds"（根据年龄、性别和观察持续时间计算），除以实际 A/E 比率。我 assemble 将它们放入列表对象并将比率四舍五入到小数点后两位。当我完成它时，我发现这些最有用的方法是对我使用 'survival' 和 'rms' 回归方法得到的结果进行敏感性检查。它们还有助于向更熟悉表格方法而不是回归的非统计受众解释结果。我现在将它作为我的 Startup .profile 的一部分。
我一直在用 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 行数据并且没有对所需输出的明确描述（至少对于非 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 函数提供粗略的事件率（在技术意义上使用 'crude'，无意贬低）是微不足道的：
by(dat, dat$GROUP, function(d) sum(d$INJURY)/sum(d$FOLLOWUP) )
#----------------
dat$GROUP: 1
[1] 0.0003636364
------------------------------------------------------
dat$GROUP: 2
[1] 0.0008032129
对应的time at risk函数（或者都打印到控制台）将是一个非常简单的修改。 'Epi' 或 'epiR' 包或其他专门用于教授基本流行病学的包之一可能会为此设计功能。 'survival' 和 'coin' 作者可能认为没有必要编写和记录如此简单的函数。
当我需要在因子协变量层内汇总实际事件与预期事件的比率时，我需要构建一个函数来正确创建事件分层表（以支持置信度估计），总和 "expecteds"（根据年龄、性别和观察持续时间计算），除以实际 A/E 比率。我 assemble 将它们放入列表对象并将比率四舍五入到小数点后两位。当我完成它时，我发现这些最有用的方法是对我使用 'survival' 和 'rms' 回归方法得到的结果进行敏感性检查。它们还有助于向更熟悉表格方法而不是回归的非统计受众解释结果。我现在将它作为我的 Startup .profile 的一部分。