如何将5种鸟类的存在/缺失纳入线性混合效应模型？

一、简化建模的核心思路 #

你目前把每种鸟类拆成单独的0/1变量纳入模型，不仅冗余，还不利于整体分析鸟类种类的效应。正确的做法是先将这些鸟类存在/缺失变量合并为单一的分类因子变量，再直接将该因子作为固定效应纳入模型，无需手动拆分哑变量。

二、具体实现步骤 #

1. 数据预处理：合并鸟类变量为因子 #

假设你的数据集里Purple、Green、Crested、Blue、Quit、Bull是0/1变量（1代表该鸟类访问过该植株），可以用以下代码生成一个名为bird_species的因子列，标记每株花被哪种鸟类访问（或未被访问）：

# 生成鸟类种类因子
data$bird_species <- apply(data[, c("Purple", "Green", "Crested", "Blue", "Quit", "Bull")], 
                          1, 
                          function(x) {
                            birds <- names(x)[x == 1]
                            if(length(birds) == 0) return("None")
                            else return(birds) # 若每株仅被单一种鸟类访问，直接返回对应种类
                          })
# 转换为因子类型
data$bird_species <- as.factor(data$bird_species)

如果存在单株被多种鸟类访问的情况，可将return(birds)替换为return(paste(birds, collapse = ","))，将访问的鸟类种类合并为一个类别。

2. 拟合优化后的线性混合模型 #

现在可以直接用合并后的因子变量建模，代码更简洁，完全匹配你最初的需求：

library(nlme)
# 拟合模型：PD.adjusted ~ VR + bird_species，Plant.ID为随机截距
nlm1 <- lme(PD.adjusted ~ VR + bird_species, 
            data = data,   
            random = ~ 1 | Plant.ID)
# 查看模型结果
m1_total <- summary(nlm1)
print(m1_total)

模型中bird_species作为分类因子，R会自动以其中一个水平为参考组（默认是字母排序最靠前的类别，比如"Blue"或"None"），输出其他水平相对于参考组的效应值，无需手动设置多个变量。

三、模型结果解读建议 #

1. 固定效应部分：

   • VR的系数：表示在控制鸟类种类的情况下，访问率每增加1单位，花粉沉积量(PD.adjusted)的平均变化量。
   • bird_species的各水平系数：比如bird_speciesPurple的系数，代表Purple鸟类访问的植株，相比参考组（比如未被任何鸟类访问的"None"），PD.adjusted的平均差异（已控制VR和植株随机效应）。
   • 每个系数的p值：用于判断该效应是否具有统计学显著性。

2. 随机效应部分：

   • Plant.ID的方差：表示不同植株之间PD.adjusted的差异程度，方差越大说明植株个体间的变异越明显。

3. 若需比较所有鸟类水平的差异，可使用事后多重比较：

library(emmeans)
# 比较不同鸟类种类的PD.adjusted均值
emmeans(nlm1, pairwise ~ bird_species)

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