从 Lexis 模型生成粗发生率 table（按因子变量分层）

Question

我正在使用 R 中的 'Epi' 包对研究中的后续数据进行建模。我对声明 Lexis 模型或运行 Poisson 和（结合生存包）Cox 回归没有任何问题。

作为初始数据审查的一部分，我想找到一种简单的方法来从 R 中的词汇模型中的数据中得出 table 未经调整的 incidence/event 比率（预拟合任何 poisson/cox款）。

我找到了一种编码方法，它允许我执行此操作并将变量分层作为探索性数据分析的一部分：

    #Generic Syntax Example
total <-cbind(tapply(lexis_model$lex.Xst,lexis_model$stratifying_var,sum),tapply(lexis_model$lex.dur,lexis_model$stratifying_var,sum))
    #Add up the number of events within the stratifying variable
    #Add up the amount of follow-up time within the stratifying the variable
    rates <- tapply(lexis_model$lex.Xst,lexis_model$stratifying_var,sum)/tapply(lexis_model$lex.dur,lexis_model$stratifying_var,sum)*10^3
    #Given rates per 1,000 person years
    ratetable <- (cbind(totals,rates))

#Specific Example based on the dataset
totals <-cbind(tapply(lexis_model$lex.Xst,lexis_model$grade,sum),tapply(lexis_model$lex.dur,lexis_model$grade,sum))
rates <- tapply(lexis_model$lex.Xst,lexis_model$grade,sum)/tapply(lexis_model$lex.dur,lexis_model$grade,sum)*10^3
ratetable <- (cbind(totals,rates))
ratetable

                 rates
1 90 20338.234 4.4251630
2 64  7265.065 8.8092811
#Shows number of events, years follow-up, number of events per 1000 years follow-up, stratified by the stratifying variable

请注意，这是粗略的 unadjusted/absolute 比率 - 不是泊松模型的输出。虽然我很欣赏上面的代码确实产生了所需的输出（并且非常简单），但我想看看人们是否知道可以获取词汇数据集并输出它的命令。我查看了 Epi 和 epitools 包中的可用命令 - 可能遗漏了一些东西但看不到明显的方法来执行此操作。

因为这是一件很常见的事情，所以我想知道是否有人知道 package/function 可以通过简单地指定词汇数据集和分层变量（或者实际上是单个函数）来做到这一点一次执行上述步骤）。

理想情况下，输出如下所示（取自 STATA，我正试图摆脱它以支持 R！）：

这里是实际数据的前二十行左右的副本（数据已经使用 Epi 包放入词汇模型中，因此所有相关的词汇变量都在那里）： https://www.dropbox.com/s/yjyz1kzusysz941/rate_table_data.xlsx?dl=0

Answer 1

我会简单地使用 tidyverse R 包来做到这一点：

library(tidyverse)
lexis_model %>%
  group_by(grade) %>%
  summarise(sum_Xst = sum(lex.Xst), sum_dur = sum(lex.dur)) %>%
  mutate(rate = sum_Xst/sum_dur*10^3) -> rateable
rateable

#  A tibble: 2 x 4
#   grade sum_Xst   sum_dur     rate
#   <dbl>   <int>     <dbl>    <dbl>
# 1     1       2 375.24709 5.329821
# 2     2       0  92.44079 0.000000

您可以自己将其包装到一个函数中：

rateFunc <- function(data, strat_var)
{
  lexis_model %>%
    group_by_(strat_var) %>%
    summarise(sum_Xst = sum(lex.Xst), sum_dur = sum(lex.dur)) %>%
    mutate(rate = sum_Xst/sum_dur*10^3)
}

然后你会调用：

rateFunc(lexis_model, "grade")

这很有用，因为结合使用 tidyverse summarise 和 mutate 可以很容易地向 table 添加更多摘要统计信息。

编辑：在澄清问题后，可以使用 popEpi 包使用 rate 命令来完成：

popEpi::rate(lexis_model, obs = lex.Xst, pyrs = lex.dur, print = grade)

# Crude rates and 95% confidence intervals:
#    grade lex.Xst  lex.dur       rate     SE.rate     rate.lo   rate.hi
# 1:     1       2 375.2472 0.00532982 0.003768752 0.001332942 0.0213115
# 2:     2       0  92.4408 0.00000000 0.000000000 0.000000000       NaN

从 Lexis 模型生成粗发生率 table（按因子变量分层）

Generate a crude incidence rate table (stratified by a factor variable) from a Lexis Model

r

rate

survival-analysis