如何为phylolm包的phyloglm获取加权Bootstrap估计值与置信区间

我用phylolm包的phyloglm()函数构建了针对二元响应变量的系统发育逻辑回归模型，简化后的代码如下：

library(phylolm)
phyloglm(y ~ x1 + x2 + x3 + x4 + X5 + X6 + x1:x2 + x1:x3 + x1:x4 + x1:X5 + x1:x6,
         method="logistic_MPLE",          
         phy=phy,
         data=dat,
         boot=10000,
         btol=30) # needed to increase btol

数据集样本量为1000，每行对应一个物种：

• y是二元响应变量，其中1占比15%、0占比85%；
• x1、x2、x3、x4、x5、x6包含定量和定性预测变量，定量变量已做0中心化缩放。

我想验证数据集更平衡时研究结果是否变化，计划执行加权Bootstrap：按响应变量y的占比反向设置抽样概率——抽中y=1的概率为0.85，抽中y=0的概率为0.15，通过有放回抽样让样本中y=1和y=0的占比平均达到50%。

对于普通LM或GLM，我可以在sample函数中添加prob参数实现加权Bootstrap，但phylolm::phyloglm没有内置加权Bootstrap选项，自行编码时遇到了核心问题：模拟样本会包含同一物种的多行数据，我不知道如何为每次模拟更新系统发育树phy。

具体疑问：

1. 能否在phylolm::phyloglm的Bootstrap中指定权重或抽样概率？
2. 如果不能，如何在包含同一物种多行数据的模拟样本上运行phyloglm？（核心是更新系统发育树的代码）
3. 是否有其他方法获取phyloglm的加权Bootstrap估计值与置信区间？

一、phyloglm内置Bootstrap不支持自定义权重/抽样概率 #

目前phyloglm的boot参数仅支持基于物种的简单有放回抽样（即每个物种被抽中的概率相等），没有内置选项允许指定自定义抽样权重或概率，无法直接通过函数参数实现加权Bootstrap。

二、处理重复物种样本的系统发育树适配方案 #

当模拟样本包含重复物种时，不需要修改系统发育树的拓扑结构——进化树描述的是物种间的进化关系，重复抽样只是增加该物种的样本权重，而非改变进化关系。这里提供两种可行思路：

思路1：借助boot包手动实现加权Bootstrap

不需要修改系统发育树，直接通过boot包自定义抽样逻辑，重复行的物种会自动在模型拟合中被赋予更高的权重，效果等价于加权拟合：

library(boot)
library(phylolm)

# 定义Bootstrap拟合函数
weighted_phyloglm_boot <- function(data, indices, phy, formula) {
  # 根据加权抽样索引获取样本
  boot_dat <- data[indices, ]
  # 直接用原始树拟合，重复物种会自动匹配树的节点
  model <- phyloglm(formula, method = "logistic_MPLE", phy = phy, data = boot_dat, btol = 30)
  # 返回模型系数
  return(coef(model))
}

# 设置抽样概率并归一化
sampling_probs <- ifelse(dat$y == 1, 0.85, 0.15)
sampling_probs <- sampling_probs / sum(sampling_probs)

# 执行加权Bootstrap
boot_result <- boot(
  data = dat, 
  statistic = weighted_phyloglm_boot, 
  R = 10000, 
  phy = phy,
  formula = y ~ x1 + x2 + x3 + x4 + X5 + X6 + x1:x2 + x1:x3 + x1:x4 + x1:X5 + x1:x6,
  prob = sampling_probs
)

# 计算置信区间
boot.ci(boot_result, type = "bca")

关键注意：确保dat中的物种ID与phy$tip.label完全匹配，否则模型无法正确匹配物种和进化树节点。

思路2：调整系统发育树枝长（可选）

如果需要从进化角度严格反映重复抽样的权重，可以给重复抽样的物种对应的枝长乘以重复次数，这种方法适用于对进化信号有特殊要求的场景：

# 假设某次Bootstrap抽样后的数据集为boot_dat，且dat包含species_id列对应phy的tip标签
species_counts <- table(boot_dat$species_id)
# 复制原始树
boot_phy <- phy
# 调整对应物种的枝长
for(sp in names(species_counts)){
  tip_idx <- which(boot_phy$tip.label == sp)
  boot_phy$edge.length[boot_phy$edge[,2] == tip_idx] <- boot_phy$edge.length[boot_phy$edge[,2] == tip_idx] * species_counts[sp]
}
# 用调整后的树拟合模型
model <- phyloglm(formula, method = "logistic_MPLE", phy = boot_phy, data = boot_dat, btol = 30)

三、替代方案：加权系统发育逻辑回归 #

如果核心需求是验证平衡数据集下的结果，也可以直接给y=1的样本分配更高权重，使用支持加权的系统发育模型包（如brms的贝叶斯框架），无需Bootstrap即可得到等效结果：

library(brms)
# 设置权重：y=1的权重为0.85/0.15≈5.67，y=0的权重为1
dat$weight <- ifelse(dat$y == 1, 0.85/0.15, 1)
# 拟合加权系统发育逻辑回归
model <- brm(
  y ~ x1 + x2 + x3 + x4 + X5 + X6 + x1:x2 + x1:x3 + x1:x4 + x1:X5 + x1:x6 + (1|gr(species_id, cov = phy)),
  data = dat,
  family = bernoulli(),
  weights = weight,
  chains = 4,
  iter = 2000
)
# 查看结果与置信区间
summary(model)

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