土壤有毒元素与表观遗传年龄加速关联分析的技术问询

研究背景 #

本研究旨在探究BKI、DAUS、SRINA、PALLA四个城市中，土壤有毒元素（共10种，采用ICP-MS测定）与表观遗传年龄加速（经实际年龄回归计算得到）的关联。

核心技术问题 #

• 四个城市样本量差异显著（BKI1000例、DAUS250例、SRINA200例、PALLA100例），如何处理该差异？
• 需要获取每个城市中10种有毒元素与表观遗传年龄的关联结果。当前使用的lm模型仅能得到整体关联，无法分城市呈现；尝试带城市交互项的lm模型后，不知如何合并回归结果与置信区间；同时不确定是否需采用GLM模型，若使用高斯族，链接函数的选择及作用也不明确。

现有代码 #

基础线性回归模型 #

library("QuantPsyc")
data <- read.table("clr.clean.file2.txt", header=T, sep=",")

model1 <- lm (epigenetic_age_acceleration ~ 
   As + Se + Fe + Co + Zn + Mn + Hg + Sb + Mo + Pb + 
   Smoking_Status + Sex + Age, data = data)
model1
model1.stat<-lm.beta (model1)

model1.stat

As
0.0256056478741109
Se
0.00499178037586947
Fe
0.00210283404005497
Co
-0.00916637143431217
Zn
0.0639371964557919
Mn
-0.0213600659139311
Hg
0.0328431516176923
Sb
0.000169338014091565
Mo
-0.0200956999960768

带城市交互项的线性回归尝试 #

IEAA_elements <- lm (epigentic_age_acceleration ~ 0 + city + 
  city:(As +Se + Fe + Co + Zn +Mn + Hg + Sb+ Mo +Pb +   
        Smoking_Status + Sex + Age), data = data)

cf <- confint(IEAA_elements)
cf

GLM模型尝试 #

# Fit a GLM
formula <- as.formula(epigentic_age_acceleration ~ 0 + city + 
  city:(As +Se + Fe + Co + Zn +Mn + Hg + Sb+ Mo +Pb +   
        Smoking_Status + Sex + Age))
model <- glm(formula, data = data, family = gaussian(link="identity"))

针对性解决方案 #

1. 样本量差异的处理方式 #

• 加权回归：给样本量小的城市赋予更高权重，抵消样本量差异带来的影响。可以用weights参数在lm或glm中实现，权重可设为各城市样本量的倒数，或基于逆概率加权（IPW）计算。
• 分层分析+元合并：对每个城市单独拟合模型，之后用固定效应/随机效应元分析方法合并结果，既能保留城市特异性，又能整合整体趋势。
• 混合效应模型：将城市设为随机效应，元素设为固定效应，同时加入城市与元素的交互项，这种模型能同时考虑组内变异和组间差异，适合样本量不均的场景。

2. 分城市获取关联结果的方法 #

方法一：拆分数据单独建模（最直观）

直接按城市拆分数据集，每个城市拟合独立线性模型，提取系数和置信区间：

# 按城市拆分数据
city_data <- split(data, data$city)

# 遍历每个城市拟合模型并整理结果
results <- lapply(city_data, function(df) {
  model <- lm(epigenetic_age_acceleration ~ As + Se + Fe + Co + Zn + Mn + Hg + Sb + Mo + Pb + Smoking_Status + Sex + Age, data = df)
  # 提取元素相关的系数与置信区间
  elem_names <- c("As", "Se", "Fe", "Co", "Zn", "Mn", "Hg", "Sb", "Mo", "Pb")
  res_df <- data.frame(
    Element = elem_names,
    Beta = coef(model)[elem_names],
    CI_Lower = confint(model)[elem_names, 1],
    CI_Upper = confint(model)[elem_names, 2],
    City = unique(df$city)
  )
  return(res_df)
})

# 合并所有城市结果
final_results <- do.call(rbind, results)
print(final_results)

方法二：处理交互项模型的结果

针对你已拟合的带城市交互项的lm模型，通过提取系数名拆分各城市结果：

# 提取系数和置信区间
coef_df <- data.frame(
  Coefficient = names(coef(IEAA_elements)),
  Beta = coef(IEAA_elements),
  CI_Lower = confint(IEAA_elements)[,1],
  CI_Upper = confint(IEAA_elements)[,2]
)

# 筛选元素与城市的交互项
elem_interaction <- grepl("city.*:(As|Se|Fe|Co|Zn|Mn|Hg|Sb|Mo|Pb)", coef_df$Coefficient)
elem_results <- coef_df[elem_interaction, ]

# 拆分城市与元素名称
elem_results$City <- gsub("city(.*):.*", "\\1", elem_results$Coefficient)
elem_results$Element <- gsub(".*:(.*)", "\\1", elem_results$Coefficient)

# 整理成清晰格式
elem_results <- elem_results[, c("City", "Element", "Beta", "CI_Lower", "CI_Upper")]
print(elem_results)

3. GLM模型的选择与链接函数 #

当因变量（表观遗传年龄加速）是连续型变量时，glm的高斯族（family=gaussian）和lm本质等价——因为高斯族默认的**恒等链接（identity）**就是直接拟合线性关系，此时glm结果和lm完全一致。只有当因变量不符合正态分布（如计数、二分类）时，才需要选择其他族和链接函数。你的场景用lm足够，无需额外使用GLM。

内容的提问来源于stack exchange，提问作者Yogesh Gupta