R中基于dredge的二项式GLM模型比较与结果保存方法咨询

问题描述 #

正在对含多参数多种汇总方式的数据集进行探索，目标为测试数据预测筛选优先级变量。计划通过dredge工具遍历所有变量组合及二阶交互项，缩小探索范围。变量分为diet、海表温度（SST）、NAO三类，要求每次建模每类仅纳入一个变量（因SST不同汇总指标高度相关，需对比方差、均值等指标效用）。目前已构建约70个全局二项式GLM模型，并完成dredge操作。

需求：

1. 提取每个全局模型中最优及AICc差值≤2的所有模型公式与AIC值，生成指定格式表格；
2. 可访问这些最优模型的汇总信息用于后续探索；
3. 咨询当前模型构建与dredge执行方法是否最优；
4. 如何保存所需输出结果。

数据集结构（前5行）：

structure(list(DP = c(1, 1, 0, 0, 0), HI = structure(c(2L, 2L, 
2L, 2L, 2L), levels = c("MR", "MSI", "SINWR"), class = "factor"), 
    NIP = c(1, 1, 1, 0, 0), PP = structure(c(2L, 1L, 1L, 1L, 
    1L), levels = c("0", "1"), class = "factor"), MPFHMOAP = c(0.750452915231012, 
    0.453965698407607, 0.441925608384321, 0.722610207052164, 
    0.441925608384321), MHVHMOAP = c(-0.0524813927346641, 1.01282281910627, 
    0.757232590524714, 0.99109247578557, 0.757232590524714), 
    MLVHMOAP = c(0.73056230124052, 1.59501115848116, 0.598198676596868, 
    0.449493863725604, 0.598198676596868), MPFHMORAP = c(0.523489867811386, 
    0.192245391224712, 0.178793838876015, 0.492383153005024, 
    0.178793838876015), MHVHMORAP = c(-0.329148414641852, 0.788691843441836, 
    0.520497046213404, 0.765889856632941, 0.520497046213404), 
    MLVHMORAP = c(0.688793519232913, 1.70845779144675, 0.532663451086749, 
    0.357258065885503, 0.532663451086749), MPFAAP = c(0.43070772291489, 
    1.13584913195813, -0.426234961686277, -0.203429863689975, 
    -0.426234961686277), MHVAAP = c(-0.0579325416076462, 0.849585831032181, 
    0.121001819024069, -0.201997642218822, 0.121001819024069), 
    MLVAAP = c(-1.35485853425042, 0.484252451276881, 0.768258290849834, 
    1.48519986147912, 0.768258290849834), MPFVAP = c(-0.244597407127982, 
    1.27217782014326, 0.247421666481899, 0.774014226262795, 0.247421666481899
    ), MHVVAP = c(-0.37754544487536, 0.758649910763721, 0.681804652562789, 
    0.887267156622061, 0.681804652562789), MLVVAP = c(-0.572980090345317, 
    1.62888586457038, -0.0115436727418632, 0.0160827313504334, 
    -0.0115436727418632), SSTANAP = c(1.07683379079566, -1.38798265446896, 
    -0.559237858951006, 0.0697567434994638, -0.559237858951006
    ), SSTVNAP = c(2.09070187356431, 0.144862160369841, -0.853503337708263, 
    -1.68867970657618, -0.853503337708263), SSTABAAP = c(0.173411254232863, 
    -1.11203442428642, -0.473195106714148, 0.827784658554605, 
    -0.473195106714148), SSTVBAAP = c(-0.629121122035084, 0.392528189012985, 
    0.379784122762904, 0.490232696930265, 0.379784122762904), 
    NAOA = c(0.897128780719261, -0.204661067596483, -1.061094353052, 
    -0.824996528412911, -1.061094353052), NAOV = c(-0.028733326081941, 
    0.855031609762443, -0.028481970866968, 0.766860571745026, 
    -0.028481970866968)), row.names = c(NA, -5L), class = c("tbl_df", 
"tbl", "data.frame"))

已运行代码片段：

# 构建diet/sst/nao所有组合的全局模型
model5_diet_HI_MPFHMOAP_SSTANAP_NAOA_DRsatr_autoglm=glm(DP~(HI+MPFHMOAP+SSTANAP+NAOA)^2, data = Model5_Diet_EnvData_SummarizedByRegionHY_RecruitmentSummaries_SubAbbrev, family = binomial(link = "logit"))

model5_diet_HI_MPFHMOAP_SSTVNAP_NAOA_DRsatr_autoglm=glm(DP~(HI+MPFHMOAP+SSTVNAP+NAOA)^2, data = Model5_Diet_EnvData_SummarizedByRegionHY_RecruitmentSummaries_SubAbbrev, family = binomial(link = "logit"))

model5_diet_HI_MPFHMOAP_SSTABAAP_NAOA_DRsatr_autoglm=glm(DP~(HI+MPFHMOAP+SSTABAAP+NAOA)^2, data = Model5_Diet_EnvData_SummarizedByRegionHY_RecruitmentSummaries_SubAbbrev, family = binomial(link = "logit"))


# 对全局模型执行dredge操作
options(na.action = "na.fail") # dredge运行必需设置

model5_diet_HI_MPFHMOAP_SSTANAP_NAOA_DRdredge=dredge(model5_diet_HI_MPFHMOAP_SSTANAP_NAOA_DRsatr_autoglm, beta = "none", evaluate = TRUE, rank = AICc)

model5_diet_HI_MPFHMOAP_SSTVNAP_NAOA_DRdredge=dredge(model5_diet_HI_MPFHMOAP_SSTVNAP_NAOA_DRsatr_autoglm, beta = "none", evaluate = TRUE, rank = AICc)

model5_diet_HI_MPFHMOAP_SSTABAAP_NAOA_DRdredge=dredge(model5_diet_HI_MPFHMOAP_SSTABAAP_NAOA_DRsatr_autoglm, beta = "none", evaluate = TRUE, rank = AICc)

options(na.action = "na.omit") # 恢复默认设置

期望输出格式：

解决方案 #

1. 提取最优模型的公式与AIC值 #

将所有dredge结果存入列表，批量筛选并整理成指定表格：

library(MuMIn)
# 收集所有dredge结果到列表
dredge_results <- list(
  model5_diet_HI_MPFHMOAP_SSTANAP_NAOA_DRdredge,
  model5_diet_HI_MPFHMOAP_SSTVNAP_NAOA_DRdredge,
  model5_diet_HI_MPFHMOAP_SSTABAAP_NAOA_DRdredge
  # 其余70个模型的dredge结果依次加入
)

# 批量筛选AICc差值≤2的模型
optimal_models_list <- lapply(dredge_results, function(x) subset(x, delta <= 2))

# 整理成指定格式表格
all_optimal_models <- do.call(rbind, lapply(optimal_models_list, function(df) {
  data.frame(
    Model = sapply(df$model, function(m) paste0("DP ~ ", as.character(formula(m))[3])),
    AIC = round(df$AICc, 1),
    stringsAsFactors = FALSE
  )
}))

# 输出表格
knitr::kable(all_optimal_models, col.names = c("Model", "AIC"))

2. 访问最优模型的汇总信息 #

dredge结果的model列存储完整模型对象，直接提取即可查看汇总：

# 查看第一个dredge结果中最优模型的汇总
summary(optimal_models_list[[1]]$model[[1]])

# 批量提取所有最优模型的汇总
all_model_summaries <- lapply(optimal_models_list, function(df) lapply(df$model, summary))

# 访问第二个dredge结果中第二个最优模型的汇总
all_model_summaries[[2]][[2]]

3. 模型构建与dredge执行的优化建议 #

批量构建模型

手动编写70个模型代码效率极低，通过变量组合循环自动生成：

# 定义三类变量全集
diet_vars <- c("HI", "NIP", "PP")
sst_vars <- c("SSTANAP", "SSTVNAP", "SSTABAAP", "SSTVBAAP") # 补充所有SST变量
nao_vars <- c("NAOA", "NAOV") # 补充所有NAO变量

# 生成所有合法变量组合
var_combinations <- expand.grid(diet = diet_vars, sst = sst_vars, nao = nao_vars, stringsAsFactors = FALSE)

# 批量构建全局模型
global_models <- list()
for(i in 1:nrow(var_combinations)) {
  vars <- c(var_combinations$diet[i], var_combinations$sst[i], var_combinations$nao[i])
  formula_str <- paste0("DP~(", paste(vars, collapse = "+"), ")^2")
  global_models[[i]] <- glm(formula_str, data = Model5_Diet_EnvData_SummarizedByRegionHY_RecruitmentSummaries_SubAbbrev, family = binomial(link = "logit"))
  names(global_models)[i] <- paste0("model_", paste(vars, collapse = "_"))
}

# 批量执行dredge
options(na.action = "na.fail")
dredge_results <- lapply(global_models, function(m) dredge(m, beta = "none", evaluate = TRUE, rank = AICc))
options(na.action = "na.omit")

控制计算量

• 先通过单变量模型筛选每类中初步有效的变量，减少组合数；
• 在dredge中设置subset参数，强制每类至少保留一个变量：

  dredge(m, beta = "none", evaluate = TRUE, rank = AICc,
         subset = (HI + NIP + PP >=1) & (SSTANAP + SSTVNAP + SSTABAAP + SSTVBAAP >=1) & (NAOA + NAOV >=1))
  

并行加速

使用parallel包启用并行计算，提升dredge运行效率：

library(parallel)
cl <- makeCluster(detectCores() - 1)
clusterEvalQ(cl, library(MuMIn))
dredge_results <- parLapply(cl, global_models, function(m) dredge(m, beta = "none", evaluate = TRUE, rank = AICc))
stopCluster(cl)

4. 保存输出结果 #

保存表格

将最优模型表格保存为CSV或Excel：

# 保存为CSV
write.csv(all_optimal_models, "optimal_models_aic.csv", row.names = FALSE)

# 保存为Excel（需安装openxlsx包）
library(openxlsx)
write.xlsx(all_optimal_models, "optimal_models_aic.xlsx")

保存模型对象

将dredge结果、最优模型等保存为RDS文件，方便后续加载：

# 保存所有dredge结果
save