基于Bootstrap的遗传力估计值组间差异检验方法咨询

我尝试比较不同处理组的遗传力估计值，该值由以下函数计算：

data <- 
  data.frame(treatment=c("control", "ivermectin", "pla2guicide", "parasitism", "control", "ivermectin", "pla2guicide", "parasitism"),
             temperature=c("29", "29", "29", "29", "33", "33", "33", "33"), 
             pronotum=c(4.12, 4.21, 4.33, 4.26, 4.34, 4.21, 4.52, 4.35))
data

her_pro <- function(data, indices) {
  mod <- lmer(pronotum ~ 1 + (1 | terrarium), data = data[indices, ])
  df <- as.data.frame(VarCorr(mod))
  family <- df[1, 4]
  residual <- df[2, 4]
  heritability <- family/(family + residual)
  return(heritability)
}

（注：原数据中terrarium列未给出，实际运行时需补充该分组变量）

我通过该函数从混合模型中得到遗传力指数，随后进行Bootstrap抽样：

BC29_p <- boot(data=data, statistic=her_pro, R=100)

现在我想了解不同处理组的遗传力值是否存在统计学差异，请问可以采用什么检验方法？我曾尝试使用finalfit包的boot_predict函数，但不确定操作是否正确。

可行的检验方法 #

1. 基于Bootstrap的组间差异检验 #

既然已经用到Bootstrap，最直接的方式是针对每个处理组独立抽样，通过差值分布判断差异：

• 先按treatment（或结合temperature）拆分数据，对每组单独做Bootstrap抽样，得到各组遗传力的Bootstrap样本分布
• 计算组间遗传力差值的Bootstrap分布，通过分布中极端值的比例得到p值
• 示例代码：

library(boot)
library(lme4)

# 按处理组分拆数据
split_data <- split(data, data$treatment)

# 对每个组执行Bootstrap抽样（建议增大R值到1000+提升精度）
boot_results <- lapply(split_data, function(x) {
  boot(data = x, statistic = her_pro, R = 1000)
})

# 以control和ivermectin组为例，计算差值的Bootstrap分布
control_her_samples <- boot_results$control$t[,1]
ivermectin_her_samples <- boot_results$ivermectin$t[,1]
diff_samples <- ivermectin_her_samples - control_her_samples

# 计算观测到的组间遗传力差值
obs_diff <- boot_results$ivermectin$t0 - boot_results$control$t0

# 双侧检验p值：差值绝对值大于观测差值的样本比例
p_value <- mean(abs(diff_samples) > abs(obs_diff))

2. 带方差结构的混合模型似然比检验 #

不用单独计算遗传力，直接在模型中检验处理组对遗传力的影响：

• 构建允许随机效应方差随处理组变化的混合模型，用似然比检验比较原模型与该模型的拟合差异，若显著则说明遗传力存在组间差异
• 示例代码：

library(lme4)

# 基础模型：所有处理组共享同一随机效应方差
base_mod <- lmer(pronotum ~ treatment + (1 | terrarium), data = data)

# 允许随机效应方差随处理组变化的模型
var_mod <- lmer(pronotum ~ treatment + (1 | terrarium), 
                data = data, weights = varIdent(form = ~1 | treatment))

# 似然比检验
anova(base_mod, var_mod)

若检验结果的p值小于显著性水平（如0.05），则支持不同处理组遗传力存在差异的结论。

3. 关于finalfit的boot_predict #

boot_predict主要用于生成模型预测值的Bootstrap置信区间，并不适合直接比较遗传力这类派生统计量的组间差异，不推荐用这个函数来完成你的需求。

内容的提问来源于stack exchange，提问作者Andrea Esquivel Román