如何循环抽样数据集并获取compar.gee的多次p值结果

需求目标 #

重复999次随机抽取100种鸟类样本，每次运行ape包中的compar.gee模型（公式：Dichromatism ~ Temp*Precip），提取结果中Pr(T > |t|)列的p值，最终输出每行对应一次抽样的Temp、Precip、Temp*Precip三个变量的p值。

遇到的问题 #

1. 误用sample_N函数导致Error in d:nrow(100) : argument of length 0错误，该函数并非R基础函数。
2. replicate代码块格式错误，触发unexpected '{'语法报错。
3. 示例进化树的edge维度与节点数量不匹配，引发维度兼容错误。

修正后的完整代码 #

1. 数据与进化树预处理（确保名称匹配） #

# 加载依赖包
library(ape)

# 示例数据（实际使用你的完整BirdData数据集）
BirdData <- structure(
    list(
      Species = c(
        "Acanthagenys_rufogularis",
        "Acanthiza_apicalis",
        "Acanthiza_chrysorrhoa",
        "Acanthiza_lineata",
        "Acanthiza_nana",
        "Acanthiza_uropygialis"
      ),
      Temp = c(27.1, 27.63, 28.8, 32.16, 29.95, 28.24),
      Precip = c(630, 475, 765, 900, 780, 590),
      Dichromatism = c(0.45, 0.3, 0.55, 0.69, 0.24, 0.58),
      row.names = c(
        "Acanthagenys_rufogularis",
        "Acanthiza_apicalis",
        "Acanthiza_chrysorrhoa",
        "Acanthiza_lineata",
        "Acanthiza_nana",
        "Acanthiza_uropygialis"
      )
    )) |> as.data.frame()

# 生成符合要求的示例进化树（实际保留你的共识树构建流程即可）
# 确保tip标签与BirdData的Species完全匹配
TrimAvianTree <- rtree(nrow(BirdData))
TrimAvianTree$tip.label <- BirdData$Species
TrimAvianTree <- compute.brlen(TrimAvianTree, method = "Grafen", power = 1)

# 检查名称匹配，确保无遗漏或不匹配项
name.check(TrimAvianTree, BirdData)

2. 重复抽样与模型运行核心代码 #

# 定义单次抽样、建模、提取p值的函数
run_single_gee <- function() {
  # 随机抽取100个样本（替换错误的sample_N为基础sample函数）
  sampled_data <- BirdData[sample(nrow(BirdData), size = 100, replace = FALSE), ]
  # 修剪进化树，匹配当前抽样的物种
  trimmed_tree <- keep.tip(TrimAvianTree, sampled_data$Species)
  # 运行compar.gee模型
  gee_result <- compar.gee(Dichromatism ~ Temp*Precip, data = sampled_data, phy = trimmed_tree)
  # 提取目标p值，排除截距项，重命名列名
  p_values <- coef(summary(gee_result))[-1, "Pr(T > |t|)"]
  names(p_values) <- c("Temp", "Precip", "Temp*Precip")
  return(p_values)
}

# 设置随机种子保证结果可重复
set.seed(123)
# 重复999次运行
all_p_values <- replicate(n = 999, run_single_gee())

# 转换为目标格式的数据框
result_df <- as.data.frame(t(all_p_values))

# 查看结果示例
head(result_df)

关键错误修正说明 #

• 抽样函数修正：用R基础的sample函数替代不存在的sample_N，正确实现行抽样。
• 代码块格式修正：将多步逻辑封装为单独函数，避免replicate内代码块的语法混乱；也可在replicate的大括号内用分号串联所有步骤。
• 进化树匹配修正：每次抽样后必须用keep.tip修剪进化树，确保树的tip标签与抽样数据的物种完全一致，避免模型运行时的名称不兼容问题。
• 示例进化树修正：原示例中edge矩阵维度与tip数量不匹配，改用rtree生成符合结构要求的示例树，实际使用时保留你的共识树构建流程即可。

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