如何实现数据框多次子抽样并重复执行PGLS模型的循环？

需求概述 #

我手里有1000个物种数据点，需要完成以下操作：

• 重复999次从总数据中随机抽取100个数据点子样本
• 每次对子样本分别运行两个PGLS模型（分别针对分光光度计、人工评分两种性二态性测量方式，分析其与海拔、纬度质心的交互关系）
• 统计两类模型中目标效应（如海拔×纬度质心的交互项）表现出显著性的次数

原模型说明 #

针对两类性二态性评分的PGLS模型如下：

分光光度计评分模型 #

PGLS_VO_Score <- gls(Colour_discriminability_Absolute ~ Altitude_Reported*Centroid.Abs, 
                     correlation = corPagel(1, phy = AvianTreeEdge, form = ~Species), 
                     data = VO_HumanScores_Merged, method = "ML")

人工评分模型 #

PGLS_Human_Score <- gls(Human_Score ~ Altitude_Reported*Centroid.Abs, 
                        correlation = corPagel(1, phy = AvianTreeEdge, form = ~Species), 
                        data = VO_HumanScores_Merged, method = "ML")

数据框VO_HumanScores_Merged包含物种名、人工评分、分光光度计评分、海拔、纬度及对应转换值（如对数转换），用于满足PGLS的模型假设。

R代码实现 #

1. 加载依赖包 #

library(nlme)   # 提供gls函数
library(ape)    # 支持corPagel系统发育相关性结构

2. 定义抽样与模型拟合函数 #

run_subsample_analysis <- function(data, n_sample = 100, phy_tree = AvianTreeEdge) {
  # 随机抽取子样本（默认无放回抽样，若需按物种抽样可修改逻辑）
  subsample <- data[sample(nrow(data), n_sample, replace = FALSE), ]
  
  # 过滤系统发育树中不存在的物种，避免模型报错
  subsample <- subsample[subsample$Species %in% phy_tree$tip.label, ]
  # 若有效样本量过少，返回NA跳过本次
  if(nrow(subsample) < 80) return(list(vo_p = NA, human_p = NA))
  
  # 拟合分光光度计评分模型，捕获报错情况
  vo_p <- tryCatch({
    model_vo <- gls(Colour_discriminability_Absolute ~ Altitude_Reported*Centroid.Abs,
                    correlation = corPagel(1, phy = phy_tree, form = ~Species),
                    data = subsample, method = "ML")
    summary(model_vo)$tTable["Altitude_Reported:Centroid.Abs", "p-value"]
  }, error = function(e) NA)
  
  # 拟合人工评分模型，捕获报错情况
  human_p <- tryCatch({
    model_human <- gls(Human_Score ~ Altitude_Reported*Centroid.Abs,
                       correlation = corPagel(1, phy = phy_tree, form = ~Species),
                       data = subsample, method = "ML")
    summary(model_human)$tTable["Altitude_Reported:Centroid.Abs", "p-value"]
  }, error = function(e) NA)
  
  return(list(vo_p = vo_p, human_p = human_p))
}

3. 重复执行999次抽样分析 #

# 设置随机种子，保证结果可重复
set.seed(123)

# 批量运行999次抽样与模型拟合
results_list <- replicate(999, run_subsample_analysis(data = VO_HumanScores_Merged), simplify = FALSE)

# 将结果整理为数据框
results_df <- do.call(rbind, lapply(results_list, function(x) data.frame(vo_p = x$vo_p, human_p = x$human_p)))

4. 统计显著性次数 #

# 设定显著性阈值（常用0.05）
alpha <- 0.05

# 统计两类模型中交互项显著的次数
vo_significant <- sum(results_df$vo_p < alpha, na.rm = TRUE)
human_significant <- sum(results_df$human_p < alpha, na.rm = TRUE)

# 输出结果
cat("分光光度计评分模型交互项显著次数：", vo_significant, "/999\n")
cat("人工评分模型交互项显著次数：", human_significant, "/999\n")

注意事项 #

• 如果数据存在同一物种的重复记录，需修改抽样逻辑，优先抽取不同物种，避免伪重复
• 若模型频繁报错，可检查子样本与系统发育树的物种匹配度，或调整最小有效样本量阈值
• 可根据研究需求，将统计目标替换为主效应（如Altitude_Reported或Centroid.Abs）的p值

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