Charls空气污染高分方法学史诗级R代码复现，学了就能发n篇公共库文章

环境领域交叉临床流行病学数据，往往联合就可以发高分SCI环境类杂志，IF分数经常>10分。可以说拿着你手中的多中心横断面、纵向队列数据就可以脑洞大开，链接任意空气污染前端数据，分析后端结局的关联关系。 

PM2.5、PM1、O3等等环境数据，是一个宏观万物互联指标可以万能连接任意疾病和结局，样本量足够大的数据集就可以干！本文是国庆期间深度写作的文章，细节超多，内容极其详细，干货满满的，喜欢就多多转发点个赞哦。 

本文内容超长，想跟着我学硬核的技术可以关注即将闪亮登场的超高性价比B款课程哦。爆款！AI+RCT+真实世界+因果机器学习+交叉学科+纵向数据，150小时

一、PM2.5系列环境数据

曾几何时，PM数据被环境专业的人士视若珍宝不肯轻易分享，毕竟这是他们吃饭的家伙。我也曾经在78年前问过很多从业者，都讳莫如深。如今随着时代发展，大家可以十分方便的拿到开源环境数据，这对本就是举步维艰的环境专业又加剧了恶性竞争。 开源数据中，CHAP: https://weijing-rs.github.io/product.html，就是十分优秀的开源集。CHAP 数据集包含 7 种主要空气污染物（即 PM1、PM2.5、PM10、O3、NO2、SO2 和 CO）， PM2.5 化学成分（即 SO42-、NO3-、NH4+、Cl-、BC 和 OM），以及环境多环芳烃（PAHs），包括 7 种致癌 PAHs（即 BaA、Chr、BbF、BkF、BaP、DahA、IcdP）。该 CHAP 数据集为公开且免费向所有用户开放！除了PM1、化学成分、PAHS等涉及作者要拿来当通讯作者筹码的高端数据集，其他作者都很大方的公开了。也足够大家发点小分杂志了。

二、范文

Pubmed搜索charls+PM2.5可以发现许多文章。我这里举几个例子，其中利用我们之前提供的数据集复现一下全文。

 三、R代码 

针对范文5，OR且仅用charls2015数据较为短平快。我复现全文如下。 去chap上下载数据，直接下载太大一百多G，只需下M、Y数据即可，不过只能手动一个个下。当然B款课程都有赠送数据集和代码，让你快速复现charls系列多篇文章。 

3.1 数据整理 

3.2 SCI基线表与纳入排除 

我尝试复现了全文纳入排除，不过还是有2000人左右的出入。流程是没太大问题的，按照纳入排除考虑的也比较周祥。步骤难点在城市获取转为百度 → 先 BD-09 → GCJ-02 → WGS84坐标系，与https://zenodo.org/communities/chap/records?q=&l=list&p=1&s=10&sort=newest 提供的chap数据集.ns进行完美匹配。

3.3 数据整理R代码 

难点在于拼音城市与国际接轨，离线maps::world.cities城市不全，自动版百度地图开放平台AK对中文识别号对拼音和英文识别弱。综上，我尝试了离线和AK在线溯源2种方法组合，以便后续国际公共数据库联合识别经纬度用。离线版完全代码见文末分享。

# #charls PM系列OR范文复现,实战医学统计
#--- 智荟charls_20112020list全集.RData
setwd("D:/聂个人文件/副业项目/广州智荟数据有限公司/A款2024课件/8.公共数据库/CHARLS")
pacman::p_load(haven,dplyr,stringr,purrr,lcmm,LCTMtools)
load("D:/聂个人文件/副业项目/广州智荟数据有限公司/A款2024课件/8.公共数据库/CHARLS/智荟charls_20112020list全集.RData")
##修复中文，修复不了
##charls_v2020[["charls_2011"]][["psu"]]$city   <- iconv(charls_v2020[["charls_2011"]][["psu"]]$city, from = "GBK", to = "UTF-8")

#-----1. 源数据体检16406人-------
#-----选取横断面2015,做体检有16406人
bio2015 <- charls_v2020[["charls_2015"]][["biomarker"]]

#--读取基线数据补齐
data.popu <- charls_v2020[["bdmph"]][["bdmph_population"]]
# ?? 删除 gender 为 NA 的行, 并 id  wave 联合重复选取唯一行
data.popu <- data.popu[!is.na(data.popu$gender) , ]

##21092人 2015基线
data.2015 <- data.popu[data.popu$wave == 2015, ] 

#-----2. 纳入排除-------
#---2015且做biomarker体检 16392人
temp.df <- data.2015 %>%
  semi_join(bio2015, by = "id")

#age>=45, 
temp.df2 <- temp.df %>%
  filter(age >= 45)  #15456人


#血液检测2011+2015
# 21097人2015年抽血，空腹血糖11363人
#---方法1
blood <- charls_v2020[["charls_derived"]][["blood"]]
#11363人
blood2015 <- blood %>%  filter(wave == 2015) %>%
  filter(bl_fasting == 1)  
names(blood2015)

#---方法2，
blood2 <- charls_v2020[["charls_2015"]][["blood"]]
#11364人差不多
blood20152 <- blood2 %>%  filter(bl_fasting == 1)  

#---#仅保留空腹检查者 10705人
temp.df3 <- temp.df2 %>%
  semi_join(blood2015, by = "id") 

#---#进行测量之前 30 分钟内，受试者是否吸烟、运动、进食和饮酒，9198人
pm_ids_clean <- charls_v2020[["charls_derived"]][["physical_measures"]] %>%
  filter(wave == 2015, bpact30 == 0) %>%
  distinct(id)                        # 只要 id
temp.df4 <- temp.df3 %>%
  semi_join(pm_ids_clean, by = "id")

#---#腰围不missing，9147人
pm_ids_clean <- charls_v2020[["charls_derived"]][["physical_measures"]] %>%
  filter(wave == 2015, !is.na(mwaist)) %>%
  distinct(id)                        # 只要 id
temp.df5 <- temp.df4 %>%
  semi_join(pm_ids_clean, by = "id")

#SBP DBP血压都不缺失，9112人
pm_ids_clean <- charls_v2020[["charls_derived"]][["physical_measures"]] %>%
  filter( wave == 2015,!is.na(systo),!is.na(diasto) ) %>%
  distinct(id)
# 只要 id
temp.df6 <- temp.df5 %>%
  semi_join(pm_ids_clean, by = "id")

#--甘油三酯（TG）、高密度脂蛋白（HDL）、空腹血糖（FBG）或 C 反应蛋白（CRP）数据不缺失
blood2015_clean <- blood2015 %>%
  filter( 
    !is.na(bl_tg),          # TG 不缺失
    !is.na(bl_hdl),         # HDL 不缺失
    !is.na(bl_glu),         # 空腹血糖 不缺失
    !is.na(bl_crp)          # CRP 不缺失
  ) %>%
  distinct(id)
# 只要 id, 9104人
temp.df7 <- temp.df6 %>%
  semi_join(blood2015_clean, by = "id")


#---?? .ns文件读取 CHAP_PM1_Y1K_2013_V3.nc 文件是 HDF5 格式的 NetCDF4 文件

#------3. psu 经纬度----
#----1?? 离线版→ psu经纬度----
#离线城市经纬度匹配（不需要上网）
# 思路：用 maps::world.cities 提供的中国（含台湾）城市经纬度做“词典”， 
#将 psu$city_eng 进行规范化处理后，与词典做模糊匹配，得到每个城市的行政中心经纬度。
pacman::p_load(dplyr, stringr, fuzzyjoin, maps, tibble)
#psu抽样省市
psu <- charls_v2020[["charls_2011"]][["psu"]]
#world.cities 数据集来自 maps 包，包含全球城市名及经纬度（long/lat），

data("world.cities", package = "maps")
gaz <- world.cities %>%
  filter(country.etc %in% c("China")) %>%   # 视分析口径而定，也可只用 China
  transmute(
    city_gaz_raw  = name,                                                # 原始城市名（英文）
    city_gaz_norm = str_to_lower(str_replace_all(name, "[^a-zA-Z0-9]", "")),  # 规范化：小写并去掉非字母数字
    longitude     = long,                                                # 经度
    latitude      = lat,                                                 # 纬度
    pop           = pop                                                  # 城市人口（用于平手时优先人口更多者）
  ) %>%
  distinct(city_gaz_norm, .keep_all = TRUE)                              # 去重：同名城市取一条（通常保留人口更大/更常见者）
# 对 psu$city_eng 做规范化
# 目标：把各种英文写法（大小写、空格、连字符、备注后缀）统一到“可比”的形式，
# 便于后续模糊匹配提高命中率。
stopifnot("city_eng" %in% names(psu))  # 保护性检查：必须有 city_eng 列
psu_norm <- psu %>%
  mutate(
    city_eng_raw  = city_eng,                                                     # 保留原始英文城市名
    city_eng_norm = str_to_lower(str_replace_all(city_eng, "[^a-zA-Z0-9]", "")), # 规范化：小写+去掉非字母数字
    # 去掉常见后缀或民族形容词（prefecture/region/league/autonomous 等），
    # 例如“chuxiong yi autonomous prefecture”去掉后缀，便于匹配到“chuxiong”
    city_eng_norm = str_replace_all(
      city_eng_norm,
      c("autonomousprefecture" = "",
        "mongolian"            = "",
        "uygur"                = "",
        "kazakh"               = "",
        "prefecture"           = "",
        "region"               = "",
        "league"               = "",
        "airport"              = "",
        "city"                 = "")
    ) %>%
      str_squish()                                                                  # 再次压缩多余空格（保险）
  )
……

#------------4. 测试.ns读取-----
#------3?? .ns翻转栅格PM1 2013-------
#r_rot（翻转/旋转之后的 PM1 栅格，已经设置了 extent 和 CRS = WGS84）
# 加载必要包,时间较久！
# options(install.packages.compile.from.source = "never")  # 禁止从源码编译
# install.packages("ncdf4")                                # 直接装ncdf4，还是后续错误
pacman::p_load(ncdf4,raster,terra,rhdf5) 
# 1. 先定义文件路径 
nc_file <- "D:/聂个人文件/副业项目/广州智荟数据有限公司/B款环境与真实世界纵向数据/环境交叉/CHAP/pm1/CHAP_PM1_Y1K_2013_V3.nc"
#--rhdf5 把.nc 里的经纬度与变量读出来交给terra 设置 extent/CRS
# 1) 看看文件里都有什么节点（找 lon/lat/变量 名称）
h5ls(nc_file)
# 2) 读经纬度与污染变量
lon <- h5read(nc_file, "lon")        # 或 "longitude"、"/lon"
lat <- h5read(nc_file, "lat")        # 或 "latitude"、"/lat"
pm  <- h5read(nc_file, "PM1")        # 比如 "PM1"；若是 3D (time×lat×lon)，再选时间层
# 3) 构造 SpatRaster
#   假设 pm 为 2D 矩阵，维度为 [lat, lon]（先确认 dim(pm) == c(length(lat), length(lon))）
r <- rast(pm)
# 4) 设定空间范围和坐标系（WGS84）
ext(r) <- c(min(lon), max(lon), min(lat), max(lat))
crs(r) <- "EPSG:4326"

# 若 lat 是从北到南（递减），需要翻转一下：
if (length(lat) > 1 && lat[2] < lat[1]) r <- flip(r, direction = "vertical")

# 5) 可视化/裁剪/输出,反了
plot(r)
# 假设原数据覆盖 70E-140E, 15N-55N
m <- as.matrix(r, wide=TRUE)
m_rot <- t(apply(m, 2, rev))   # 顺时针旋转 90°
r_rot <- rast(m_rot)
ext(r_rot) <- c(70, 140, 15, 55)  # 手动补上经纬度范围
crs(r_rot) <- "EPSG:4326"
plot(r_rot)


#-----4?? 链接 psu_geo_pm1 2013-------
pacman::p_load(terra, sf, dplyr)
# 1) 把 psu_geo 转换为 SpatVector 点图层
pts <- terra::vect(
  psu_geo,
  geom = c("wgs_lon", "wgs_lat"),   # 列名要和 psu_geo 里的一致
  crs = "EPSG:4326",
  keepgeom = TRUE
)
# 2) 从栅格提取对应点的 PM1 值
pm1_2013 <- terra::extract(r_rot, pts, ID=FALSE)  
# 返回一个数据框，按点顺序
# 3) 合并回 psu_geo
psu_geo_pm1 <- psu_geo %>%
  bind_cols(pm1_2013 = pm1_2013[,1])  

#------------6.  最终纳排→ temp.df8-----
# temp.df7 合并psu_geo_all
pacman::p_load(dplyr, stringr)

# 1) 构造 temp.df7 id匹配前7位数
temp.df7_keyed <- temp.df7 %>%
  mutate(
    id_chr = as.character(id),
    comm7  = str_sub(id_chr, 1, 7)   # 取前7位
  )
# 2) 统一 psu_geo_all 的键类型
psu_geo_all_keyed <- psu_geo_all %>%
  mutate(communityid = as.character(communityid))
# 3) 如 communityid 有重复，先去重（保留第一次；可改为其它规则）
dup_check <- psu_geo_all_keyed %>%
  count(communityid, name = "n") %>%
  filter(n > 1)
psu_geo_all_dedup <- psu_geo_all_keyed %>%
  distinct(communityid, .keep_all = TRUE)
# 4) 左连接：把 psu_geo_all 的全部列并到 temp.df7
temp.df8 <- temp.df7_keyed %>%
  left_join(psu_geo_all_dedup, by = c("comm7" = "communityid")) %>%
  dplyr::select(-id_chr, -comm7)   # 清理临时列
# # 可选：看看匹配率
# match_rate <- mean(!is.na(temp.df8$province_eng))  # 用其中一列做判断
# message(sprintf("匹配率：%.1f%%", 100 * match_rate))




#-----------7. ?? paper数据 -----
##-----7.1 x covs-------
names(temp.df8) 
#新增东中西地区region，省份划分
temp.df8 <- temp.df8 %>%
  mutate(region = case_when(
    # 东部
    province_eng %in% c("Beijing","Tianjin","Hebei","Shanghai",
                        "Jiangsu","Zhejiang","Fujian","Shandong",
                        "Guangdong","Hainan") ~ "East",
    # 中部
    province_eng %in% c("Shanxi","Anhui","Jiangxi",
                        "Henan","Hubei","Hunan") ~ "Midland",
    # 西部
    province_eng %in% c("Inner Mongolia","Guangxi","Chongqing","Sichuan",
                        "Guizhou","Yunnan","Tibet","Shaanxi",
                        "Gansu","Qinghai","Ningxia","Xinjiang") ~ "West",
    # 东北（根据你需要决定并到哪一类）
    province_eng %in% c("Liaoning","Jilin","Heilongjiang") ~ "Midland", 
    #province_eng %in% c("Liaoning","Jilin","Heilongjiang") ~ "Northeast"
    TRUE ~ NA_character_
  ))
paper.1 <-  temp.df8 %>%
  dplyr::select(
    id,
    age,
    gender,
    residence = residence_tag,
    education = edu_1_tag,   # 或改成 edu_2_tag
    marital = marry_1_tag,   # 或组合 marry_1/2/3
    household_exp = income_household,
    smoking = smoke_1_tag,   # 或结合 smoke_2_tag
    drinking = drink_1_tag,  # 或结合 drink_2/3
    physical_activity = physical_activity_tag,
   ## cooking_fuel = Cooking_fuel_type, # 待确认 
    region,
    waist,
    tg = bl_tg,
    hdl = bl_hdl,
    fbg = bl_glu,
    sbp,
    dbp,
    crp = bl_crp, #核心x变量
    anti_htn = itvtn_hypertension_drug_tag,
    anti_dm = itvtn_insulin_or_diabetes_drug_tag 
   ## mets  #你的Y变量 
  )
#引入固态燃烧
fuel <- charls_v2020[["charls_derived"]][["financial_housing"]]
# 1) 取 fuel 中 wave==2015，留下 id 和 cooking_fuel_type，并对 id 去重
fuel_2015 <- fuel %>%
  filter(wave == 2015) %>%
  select(id, cooking_fuel_type) %>%
  distinct(id, .keep_all = TRUE)
# 2) 与 paper 按 id 左连接，只新增 cooking_fuel_type
paper.2 <- paper.1 %>%
  left_join(fuel_2015, by = "id")


##-----7.2 y→ mets------
library(dplyr)
library(stringr)

df <- temp.df8 %>%
  mutate(
    # 0) 自动判断单位（粗略但实用）
    tg_unit  = ifelse(median(bl_tg,  na.rm=TRUE)  > 10, "mgdl", "mmol"),
    hdl_unit = ifelse(median(bl_hdl, na.rm=TRUE)  > 10, "mgdl", "mmol"),
    glu_unit = ifelse(median(bl_glu, na.rm=TRUE)  > 20, "mgdl", "mmol"),
    
    tg_thr   = ifelse(tg_unit  == "mgdl", 150, 1.7),
    hdl_m_thr= ifelse(hdl_unit == "mgdl",  40, 1.04),
    hdl_f_thr= ifelse(hdl_unit == "mgdl",  50, 1.29),
    glu_thr  = ifelse(glu_unit == "mgdl", 100, 5.6),
    
    # 1) 性别到 1/2
    gender_num = case_when(
      is.numeric(gender) ~ as.numeric(gender),
      is.character(gender) & str_to_lower(gender) %in% c("male","m","男") ~ 1,
      is.character(gender) & str_to_lower(gender) %in% c("female","f","女") ~ 2,
      TRUE ~ NA_real_
    ),
    
    # 2) 五个成分
    high_wc = case_when(
      !is.na(gender_num) & !is.na(waist) &
        ((gender_num==1 & waist>=85) | (gender_num==2 & waist>=80)) ~ 1L,
      !is.na(gender_num) & !is.na(waist) ~ 0L,
      TRUE ~ NA_integer_
    ),
    
    # 药物字段更稳健：允许 1/0 或 "Yes/No"
    anti_htn_flag = case_when(
      is.na(itvtn_hypertension_drug_tag) ~ NA,
      itvtn_hypertension_drug_tag %in% c(1,"1","Yes","yes","Y","y","是","有") ~ 1L,
      itvtn_hypertension_drug_tag %in% c(0,"0","No","no","N","n","否","无") ~ 0L,
      TRUE ~ NA_integer_
    ),
    anti_dm_flag = case_when(
      is.na(itvtn_insulin_or_diabetes_drug_tag) ~ NA,
      itvtn_insulin_or_diabetes_drug_tag %in% c(1,"1","Yes","yes","Y","y","是","有") ~ 1L,
      itvtn_insulin_or_diabetes_drug_tag %in% c(0,"0","No","no","N","n","否","无") ~ 0L,
      TRUE ~ NA_integer_
    ),
    
    high_bp = case_when(
      (!is.na(sbp) & sbp>=130) | (!is.na(dbp) & dbp>=85) | anti_htn_flag==1L ~ 1L,
      (!is.na(sbp) | !is.na(dbp) | !is.na(anti_htn_flag)) ~ 0L,
      TRUE ~ NA_integer_
    ),
    
    low_hdl = case_when(
      !is.na(gender_num) & !is.na(bl_hdl) &
        ((gender_num==1 & bl_hdl < hdl_m_thr) | (gender_num==2 & bl_hdl < hdl_f_thr)) ~ 1L,
      !is.na(gender_num) & !is.na(bl_hdl) ~ 0L,
      TRUE ~ NA_integer_
    ),
    
    high_tg = case_when(
      !is.na(bl_tg) & bl_tg >= tg_thr ~ 1L,
      !is.na(bl_tg) ~ 0L,
      TRUE ~ NA_integer_
    ),
    
    high_glu = case_when(
      (!is.na(bl_glu) & bl_glu >= glu_thr) | anti_dm_flag==1L ~ 1L,
      (!is.na(bl_glu) | !is.na(anti_dm_flag)) ~ 0L,
      TRUE ~ NA_integer_
    )
  ) %>%
  mutate(
    # 缺失不计分（若要“有缺失则NA”，用你注释里那段替换）
    MetS_risk_num = rowSums(across(c(high_wc,high_bp,low_hdl,high_tg,high_glu)), na.rm = TRUE),
    mets = case_when(
      !is.na(MetS_risk_num) & MetS_risk_num >= 3 ~ 1L,
      !is.na(MetS_risk_num) ~ 0L,
      TRUE ~ NA_integer_
    )
  )

# 各成分阳性率
sapply(df[,c("high_wc","high_bp","low_hdl","high_tg","high_glu")], function(x) mean(x==1, na.rm=TRUE))
# 总体 MetS 率
mean(df$mets==1, na.rm=TRUE)
# 组成数分布
table(df$MetS_risk_num, useNA="ifany")


#----?? 最终论文数据-----
paper <- paper.2 %>%
  left_join(df %>% select(id, mets), by = "id")

#----?? 最终论文数据-----
paper <- paper.2 %>%
  left_join(df %>% select(id, mets), by = "id")

3.4 统计R代码 

最终分析数据与过程如下。

OK，完美复刻且更优美，你可以用代码替换任意结局变量。

#-----------8. ?? 智荟 table1 -----
str(paper)
convars <- c("age", "waist", "sbp", "dbp", "tg", "hdl", "fbg", "crp")

#--physical_activity >40% 删除此变量
catvars <- c("gender", "residence", "education", "marital",
             "smoking", "drinking",  ##"physical_activity"移出
             "region", "anti_htn", "anti_dm", "cooking_fuel_type")
paper[convars] <- lapply(paper[convars], as.numeric)
paper[catvars] <- lapply(paper[catvars], factor)
#分组yvar
yvar = "mets"
##time="ftime"
##xvar="sleeptime"
# 确保 time 变量不包含在 convars 中!
##convars <- setdiff(convars, "time")
#除了y 和 time变量= vars,x在var里面！
##vars <- c(catvars,convars)
# covs是排除主要因素 xvar
##covs <- setdiff(vars, xvar)
# 查看缺失比例 
colMeans(is.na(paper[, c(convars, catvars, yvar)]))
library(VIM)
aggr(paper[, c(convars, catvars, yvar)], numbers = TRUE, prop = TRUE)
…
#执行 5 重插补 + 合并（连续取均值、因子取众数）
MI_data <- fix_na(factorCols = factorCols, data = missdata, n = 5)
paper_MI <- paper
paper_MI[, vars_use] <- MI_data[, vars_use]

## 简单检查，缺失0
#sapply(MI_data[, convars], function(x) mean(is.na(x)))  # 连续变量缺失应接近 0
#sapply(MI_data[, factorCols], function(x) mean(is.na(x)))  # 因子缺失应接近 0
#table(paper_MI$mets, useNA = "ifany")

四、发散思维

“

1.任意队列或公共数据库其实都可以结合chap进行PM系列探索，上游是万能的，下游结局任意y就可无限重复刷文章。

“

2.OR HR RR都可以，PM系列的选择也可以多样，可以传统aOR、aHR、也可广义倾向评分更前沿技术探索PM系列因果推断(B款特色，近期热门发文)。

3.复现charls，原文数据不一致，也是公共数据库矛盾之一。严谨溯源耐心缺一不可。 

4.案例过程成熟溯源，举一反三，实战公共数据库系列，且可复制推广超多idea复刻。B款课程提供公共数据库详细介绍与前沿方法，让你脑洞大开。

五、小结

“

Charls等等公共库+PM2.5，环境交叉发高分神器，跟着代码一路next狂奔吧。

福利，本文文献与chap部分代码获取方法，发送公众号聊天：1009