[使用下面的示例数据，我试图使用R中的ageadjust.indirect和ageadjust.direct函数，按年，地区和性别等子组来计算直接和间接的标准化比率，但是运气不高。感谢您的帮助。
#study Data
data <- data.frame (
year = rep(c(rep("2010",20), rep("2011",20)),1),
subregion = rep(c(rep("A",10), rep("B",10)),2),
gender=rep(c(rep("M",5), rep("F",5)),4),
agegr = rep(c('00-14','15-34','35-54','55-74','75+'),8),
pop = round((runif(40, min = 10000, max = 99999)), digits = 0),
count = round((runif(40, min = 100, max = 999)), digits = 0)
)
data$rate = data$count / data$pop
# standard data
stdata <- data.frame(age=rep(c('00-14','15-34','35-54','55-74','75+'),2),
sex=rep(c('M','F'),c(5,5)),
pop=c(308543,401996,409450,199486,610631,
293991,388762,418814,227170,104944)
)
#implement direct age standardization using 'ageadjust.direct'
library(epitools)
dsr<- ageadjust.direct(count = data$count, pop = data$pop,
, stdpop = stdata$pop)
round(100000*dsr, 2) ##rate per 100,000 per year using standard pop
如何按年份分地区和性别获取dsr？
#implement indirect age standardization using 'ageadjust.indirect'
dataB<-subset(data,subregion == 'B')
isr <- ageadjust.indirect(count = data$count, pop = data$pop,
stdcount = dataB$count, stdpop = dataB$pop)
round(isr$sir, 2) ##standarized incidence ratio
round(100000*isr$rate, 2) ##rate per 100,000 per year
如何按年份，分区和性别获得isr？
我正在尝试使用“ epitools”软件包按年份，分区和性别来计算dsr和isr。使用以下Laura的帖子：Calculate age standardised rates by sub-group with confidence intervals in R
I used Arkun code which worked for dsr, but I could not run it for isr. Please look at below the errors.
# how to get dsr by year subregion, and gender?
library(epitools)
library(tidyverse)
test1a<-data %>%
group_by(year,subregion,gender) %>%
summarise(age_adjust = list(ageadjust.direct(count = count,
pop = pop, stdpop = stdata$pop))) %>%
mutate(age_adjust = map(age_adjust, as.data.frame.list)) %>%
unnest(c(age_adjust))
> test1a
# A tibble: 8 x 7
# Groups: year, subregion [4]
year subregion gender crude.rate adj.rate lci uci
<fct> <fct> <fct> <dbl> <dbl> <dbl> <dbl>
1 2010 A F 0.0132 0.0205 0.0199 0.0211
2 2010 A M 0.0110 0.0147 0.0142 0.0153
3 2010 B F 0.00662 0.00780 0.00751 0.00809
4 2010 B M 0.0207 0.0196 0.0190 0.0201
5 2011 A F 0.0149 0.0206 0.0201 0.0212
6 2011 A M 0.00901 0.00911 0.00884 0.00938
7 2011 B F 0.0108 0.0119 0.0115 0.0123
8 2011 B M 0.0160 0.0283 0.0275 0.0290
> round(100000*test1a$crude.rate, 1) ##rate per 100,000 per year
[1] 1323.6 1099.0 661.6 2068.4 1490.9 900.5 1077.7 1604.1
# how to get isr by year, subregion, and sex?
library(epitools)
library(tidyverse)
test1b<-data %>%
group_by(year,subregion,gender) %>%
summarise(age_adjustind = list(ageadjust.indirect(count = data$count, pop = data$pop,
stdcount = dataB$count, stdpop = dataB$pop))) %>%
mutate(age_adjustind = map(age_adjustind, as.data.frame.list)) %>%
unnest(c(age_adjustind))
Error in (function (..., row.names = NULL, check.rows = FALSE, check.names = TRUE, :
arguments imply differing number of rows: 5, 4
I tried using Parfait code for dsr and isr, but I could not run it for isr as well. Please look at the errors below.
# how to get dsr by year subregion, and gender?
age_adjust_test_list <- by(data, data[,c("year", "subregion", "gender")], function(sub) {
tmp <- ageadjust.direct(count = sub$count, pop = sub$pop,
rate = sub$rate, stdpop = stdata$pop)
data.frame(year = max(as.character(sub$year)),
subregion = max(as.character(sub$subregion)),
gender = max(as.character(sub$gender)),
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)
> final_df
year subregion gender crude_rate adj_rate lower_CI upper_CI
1 2010 A F 0.013236260 0.020528363 0.019932496 0.02113975
2 2011 A F 0.014908717 0.020642121 0.020074390 0.02122418
3 2010 B F 0.006615836 0.007795789 0.007513140 0.00808890
4 2011 B F 0.010776963 0.011882004 0.011472889 0.01230821
5 2010 A M 0.010989512 0.014720892 0.014195311 0.01526214
6 2011 A M 0.009005453 0.009105862 0.008839246 0.00938093
7 2010 B M 0.020683981 0.019563740 0.019025361 0.02011709
8 2011 B M 0.016041262 0.028282079 0.027540780 0.02904046
# how to get isr by year, subregion, and sex?
age_adjustin_test_list <- by(data, data[,c("year", "subregion", "gender")], function(sub) {
tmp <- ageadjust.indirect(count = sub$count, pop = sub$pop,
stdcount = dataB$count, stdpop = dataB$pop)
data.frame(year = max(as.character(sub$year)),
subregion = max(as.character(sub$subregion)),
gender = max(as.character(sub$gender)),
observed = tmp[[1]],
exp = tmp[[2]],
sir = tmp[[3]],
lci = tmp[[4]],
uci = tmp[[5]])
})
Error in tmp[[3]] : subscript out of bounds
final_df <- do.call(rbind, age_adjustin_test_list)
#Individual region isr calculation.
dataA<-subset(data,subregion == 'A'& year=='2010')
dataB<-subset(data,subregion == 'B'& year=='2010')
dataB <- rename(dataB, count.b='count', pop.b='pop')
# region A
isr1 <- ageadjust.indirect(count = dataA$count, pop = dataA$pop,
stdcount = dataB$count.b, stdpop = dataB$pop.b)
round(isr1$sir, 2) ##standarized incidence ratio
observed exp sir lci uci
5852.00 9659.27 0.61 0.59 0.62
# region B
isr2 <- ageadjust.indirect(count = dataB$count.b, pop = dataB$pop.b,
stdcount = dataB$count.b, stdpop = dataB$pop.b)
round(isr2$sir, 2) ##standarized incidence ratio
round(isr2$sir, 2) ##standarized incidence ratio
observed exp sir lci uci
4757.00 4757.00 1.00 0.97 1.03
## isr by region: it does not calculate expected count an isr for region A. Please look at the output below.
data2<-subset(data,year=='2010')
library(epitools)
library(dplyr)
library(purrr)
left_join(data2, dataB %>%
select(-subregion)) %>%
group_by(subregion) %>%
summarise(age_adjustind = list(ageadjust.indirect(count = count,
pop = pop, stdcount = count.b, stdpop = pop.b))) %>%
mutate(age_adjustind = map(age_adjustind, ~
map_dfr(.x, as.data.frame.list))) %>%
unnest(c(age_adjustind))
Joining, by = c("year", "gender", "agegr", "rate")
# A tibble: 4 x 8
subregion observed exp sir lci uci crude.rate adj.rate
<fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 A 5852 NA NA NA NA NA NA
2 A NA NA NA NA NA 0.00945 NA
3 B 4757 4757 1 0.972 1.03 NA NA
4 B NA NA NA 0.0111 0.0117 0.0114 0.0114
group_split，然后在每个list元素中应用该函数library(epitools)
library(dplyr)
library(purrr)
out <- data %>%
group_split(subregion) %>%
map_dfr(~
.x %>%
group_by(year, gender) %>%
summarise(age_adjustind =
list(ageadjust.indirect(count = count, pop = pop,
stdcount = .x$count, stdpop =.x$pop))))