年龄性别标化率Age-sex standardized rate

R实现

pacman::p_load(
  rio,                 # import/export data
  xlsx,haven,readr,
  here,                # locate files
  tidyverse,           # data management and visualization
  stringr,             # cleaning characters and strings
  frailtypack,         # needed for dsr, for frailty models
  dsr,                 # standardise rates
  PHEindicatormethods) # alternative for rate standardisation
#出错使用旧版dsr包
packageurl <- "https://cran.r-project.org/src/contrib/Archive/dsr/dsr_0.2.2.tar.gz"
install.packages(packageurl, repos=NULL, type="source")

# conntries*age*sex 组合population每行唯一, 分母
pop_countries <- A_demo %>%  # begin with country A dataset
  bind_rows(B_demo) %>%        # bind rows, because cols are identically named
  pivot_longer(                       # pivot longer
    cols = c(m, f),                   # columns to combine into one
    names_to = "Sex",                 # name for new column containing the category ("m" or "f") 
    values_to = "Population") %>%     # name for new column containing the numeric values pivoted
  mutate(Sex = recode(Sex,            # re-code values for clarity
                      "m" = "Male",
                      "f" = "Female"))
# conntries*age*sex 组合death每行唯一, 分子
deaths_countries <- A_deaths %>%    # begin with country A deaths dataset
  bind_rows(B_deaths) %>%        # bind rows with B dataset, because cols are identically named
  pivot_longer(                  # pivot longer
    cols = c(Male, Female),        # column to transform into one
    names_to = "Sex",              # name for new column containing the category ("m" or "f") 
    values_to = "Deaths") %>%      # name for new column containing the numeric values pivoted
  rename(age_cat5 = AgeCat)      # rename for clarity

#生成
country_data <- pop_countries %>% 
  left_join(deaths_countries, by = c("age_cat5", "Sex"))

#DT::datatable(country_data)

#按国家、年龄类别和性别对其进行排序; 目的是进行AB数据库合并,
country_data <- country_data %>% 
  mutate(
    Country = fct_relevel(Country, "A", "B"),
      Sex = fct_relevel(Sex, "Male", "Female"),
        age_cat5 = fct_relevel(
      age_cat5,
      "0-4", "5-9", "10-14", "15-19",
      "20-24", "25-29",  "30-34", "35-39",
      "40-44", "45-49", "50-54", "55-59",
      "60-64", "65-69", "70-74",
      "75-79", "80-84", "85")) %>%  #其实85就是85+
arrange(Country, age_cat5, Sex)
#===========AB合并 country_data===========
#===========标准人口standard_pop_clean ===========
#DT::datatable(standard_pop_data)
# 移除一些不一致多余的字符, 标准世界人口85+保存为85
library(stringr)
standard_pop_clean <- standard_pop_data %>%
  mutate(
    AgeGroup = str_replace_all(AgeGroup, "years", ""),   # 去除 "year"
    AgeGroup = str_replace_all(AgeGroup, "plus", ""),    # 去除 "plus"
    AgeGroup = str_replace_all(AgeGroup, " ", ""))%>%    # 去除 " " space
    rename(pop = WorldStandardPopulation) %>%   # change col name to "pop", as this is expected by dsr package
  rename(age_cat5 = AgeGroup)
#    AgeGroup = str_replace_all(AgeGroup, "\\+", ""))  %>%   # 去除+号,比较特别\\
#重新赋值pop 取代 worldstandardpopulation

all_data <- left_join(country_data, standard_pop_clean, by=c("age_cat5", "Sex"))
#DT::datatable(all_data)

#===========dsr包进行标化 age*sex，需要country_data  ===========
# ①  直接法标化率 age+sex 
mortality_rate <- dsr::dsr(
  data = country_data,  # AB合并库
  event = Deaths,        #分子 
  fu = Population,    #“follow-up” 设为分母  
  subgroup = Country,   
  age_cat5,             
  Sex,
  refdata = standard_pop_clean, # 标化人口库
  method = "gamma",      # gamma 分布估计 95% CI
  sig = 0.95,            
  mp = 100000,           # 计算10w分之率
  decimals = 2)          # 小数位2)

DT::datatable(mortality_rate)


#===========dsr包还可以计算 RR RD  ===========
# A 相对B 的RR 95%CI，mortality_rr，说明AB标化率有差异,A>B
mortality_rr <- dsr::dsrr(
  data = country_data,  # AB合并库
  event = Deaths,        #分子 
  fu = Population,    #“follow-up” 设为分母  
  subgroup = Country,   
  age_cat5,             
  Sex,
  refdata = standard_pop_clean, # 标化人口库
  refgroup = "B",
  estimate = "ratio",
  sig = 0.95,            
  mp = 100000,           # 计算10w分之率
  decimals = 2)          # 小数位2)

DT::datatable(mortality_rr)

#===========dsr包还可以计算 RR RD  ===========
# A 相对B 的RD 95%CI，mortality_rd，说明AB标化率有差异,A>B
mortality_rd <- dsr::dsrr(
  data = country_data,  # AB合并库
  event = Deaths,        #分子 
  fu = Population,    #“follow-up” 设为分母  
  subgroup = Country,   
  age_cat5,             
  Sex,
  refdata = standard_pop_clean, # 标化人口库
  refgroup = "B",
  estimate = "difference",
  sig = 0.95,            
  mp = 100000,           # 计算10w分之率
  decimals = 2)          # 小数位2)

DT::datatable(mortality_rd)

#------------ PHEindicatormethods包 可以直接法+间接法--------------
# 直接法， 利用合并AB+标化人口数据的合集库 all_data 
library(dplyr)
#有 + 号需要处理掉；dobson法
dat <-  filter(all_data, !grepl("\\+",age_cat5))
mortality_ds_rate_phe <- dat %>%  
  group_by(Country) %>%
  PHEindicatormethods::phe_dsr(
    x = Deaths,                 
    n = Population,             
    stdpop = pop,              
    stdpoptype = "field")        

DT::datatable(mortality_ds_rate_phe)

#间接法，#使用国家 B 作为参考人口计算国家 A 的比率
refpopCountryB <- country_data %>% 
  filter(Country == "B") 
mortality_is_rate_phe_A <- country_data %>%
  filter(Country == "A") %>%
  PHEindicatormethods::phe_isr(
    x = Deaths,                 # column with observed number of events
    n = Population,             # column with non-standard pops for each stratum
    x_ref = refpopCountryB$Deaths,  # reference number of deaths for each stratum
    n_ref = refpopCountryB$Population)  # reference population for each stratum

#knitr::kable(mortality_is_rate_phe_A)
DT::datatable(mortality_is_rate_phe_A)