按 R 中的置信区间按子组计算年龄标准化率答案

【问题标题】：Calculate age standardised rates by sub-group with confidence intervals in R按 R 中的置信区间按子组计算年龄标准化率
【发布时间】：2018-05-14 03:24:21
【问题描述】：

我有一个如下所示的数据框：

df <- data.frame (
  time = rep(c("2010", "2011", "2012", "2013", "2014"),4),
  age = rep(c("40-44", "45-49", "50-54", "55-59", "60-64"),4),
 weight = rep(c(0.38, 0.23, 0.19, 0.12, 0.08),4),
 ethnic = rep(c(rep("M",5),rep("NM",5)),2),
 gender = c(rep("M",10), rep("F",10)),
 pop = round((runif(10, min = 10000, max = 99999)), digits = 0),
 count = round((runif(10, min = 100, max = 999)), digits = 0)
)
df$rate = df$count / df$pop

我想计算直接年龄标准化发病率，其中发病率 = 计数/人口），以及这些的置信区间；对于每个子组。因此，对于时间、性别、种族、年龄的每种组合，我都会有一个标准化的费率。有没有办法在 R 中做到这一点？

我尝试使用 R 包 {epitools} 中的函数 ageadjust.direct，如下所示：

age_adjust_test <- ageadjust.direct(count = df$count, pop = df$pop, 
                                rate = df$rate, stdpop = df$weight)

由此得出的输出是总体调整率、置信区间和粗率。有没有办法让每个子组得到这个输出？

【问题讨论】：

包是epitools
谢谢，您指的子组是什么。是按种族、性别和年龄划分的吗？
我想要每个年龄、时间、种族和性别组合的发病率。所以对于性别 = M，种族 = M，年龄 = 40-44，年份 = 2010；性别 = M，种族 = M，年龄 = 40-44，年份 = 2011 等等

标签： r confidence-interval rate

【解决方案1】：

我们可以将summarise 组合成list，然后将unnest 和list 组件组合成单独的列

library(tidyverse)
df %>% 
   group_by(time,age, ethnic, gender) %>% 
   summarise(age_adjust = list(ageadjust.direct(count = count,
            pop = pop, rate = rate, stdpop = weight))) %>%
   mutate(age_adjust = map(age_adjust, as.data.frame.list))  %>% 
   unnest
# A tibble: 20 x 8
# Groups:   time, age, ethnic [10]
#   time  age   ethnic gender crude.rate adj.rate     lci     uci
#   <fct> <fct> <fct>  <fct>       <dbl>    <dbl>   <dbl>   <dbl>
# 1 2010  40-44 M      F         0.00763  0.00763 0.00709 0.00820
# 2 2010  40-44 M      M         0.00763  0.00763 0.00709 0.00820
# 3 2010  40-44 NM     F         0.0281   0.0281  0.0257  0.0306 
# 4 2010  40-44 NM     M         0.0281   0.0281  0.0257  0.0306 
# 5 2011  45-49 M      F         0.0145   0.0145  0.0136  0.0155 
# 6 2011  45-49 M      M         0.0145   0.0145  0.0136  0.0155 
# 7 2011  45-49 NM     F         0.0425   0.0425  0.0399  0.0453 
# 8 2011  45-49 NM     M         0.0425   0.0425  0.0399  0.0453 
# 9 2012  50-54 M      F         0.0116   0.0116  0.0109  0.0124 
#10 2012  50-54 M      M         0.0116   0.0116  0.0109  0.0124 
#11 2012  50-54 NM     F         0.00708  0.00708 0.00607 0.00821
#12 2012  50-54 NM     M         0.00708  0.00708 0.00607 0.00821
#13 2013  55-59 M      F         0.0251   0.0251  0.0232  0.0271 
#14 2013  55-59 M      M         0.0251   0.0251  0.0232  0.0271 
#15 2013  55-59 NM     F         0.00733  0.00733 0.00678 0.00792
#16 2013  55-59 NM     M         0.00733  0.00733 0.00678 0.00792
#17 2014  60-64 M      F         0.0101   0.0101  0.00944 0.0109 
#18 2014  60-64 M      M         0.0101   0.0101  0.00944 0.0109 
#19 2014  60-64 NM     F         0.00916  0.00916 0.00852 0.00984
#20 2014  60-64 NM     M         0.00916  0.00916 0.00852 0.00984

【讨论】：

谢谢。运行此代码时出现此错误：“mutate_impl(.data, dots) 中的错误：列 adjust 的长度必须为 1（组大小），而不是 4”
当我在完整数据集上运行此代码时，它可以工作，但它会创建一个名为“age_adjust”的新变量，该变量的值是所有四个输出值（速率、调整后的速率）的列表、LCI 和 UCI）。有没有办法获得四个新变量，分别对应汇率、调整汇率、LCI 和 UCI？

【解决方案2】：

只需使用by 按一个或多个因素对数据帧进行子集化，然后将子集传递到您的函数中。在这里，by 将使用docs page 上显示的函数值返回数据帧列表。在by 之外，您可以使用do.call(rbind,...) 将所有 dfs 绑定到一个最终数据帧中。

age_adjust_test_list <- by(df, df[,c("time", "gender", "ethnicity", "age")], function(sub) {
                            tmp <- ageadjust.direct(count = sub$count, pop = sub$pop, 
                                                    rate = sub$rate, stdpop = sub$weight)

                            data.frame(time = max(sub$time),
                                       gender = max(sub$gender),
                                       ethnicity = max(sub$ethnicity),
                                       age = max(sub$age),
                                       crude_rate = tmp[[1]],
                                       adj_rate = tmp[[2]],
                                       lower_CI = tmp[[3]],
                                       upper_CI = tmp[[4]])
                           })

final_df <- do.call(rbind, age_adjust_test_list)

NULL 将显示在数据框中未表示的组合。所以考虑根据需要过滤掉：

age_adjust_test_list <- Filter(function(x) !is.null(x), age_adjust_test_list)

【讨论】：

非常感谢！有没有办法将此输出作为表格或数据框？
查看更新构建数据框列表。我正在使用docs 页面，假设从方法返回了一个向量列表。
谢谢。当我运行此代码时，我收到以下错误：“tmp$crude.rate 中的错误：$ 运算符对原子向量无效”我认为我需要一种不同的方式来为临时对象下标。你知道怎么做吗？
我通过更改 crude_rate = tmp[[1]], adj_rate = tmp[[2]] 等来实现此功能 - 通过正确索引临时对象 (tmp)
知道了。因此，它们没有像文档中提到的那样被命名为向量。这完全是猜测，因为我没有安装软件包。相应地进行了编辑。很高兴它成功了。编码愉快！

【解决方案3】：

这里有个方便的data.table方式，一行就够了。

library(data.table)
library(epitools)
# convert df to data.table
setDT(df)
# define subgroups
group_by<-c('time','age', 'ethnic', 'gender')

# ageadjust.direct by subgroups. The trick is to include as.list()
df[, as.list(ageadjust.direct(count = count, pop = pop, rate = rate, stdpop = weight)), by=group_by]

【讨论】：