能否使用MuMIn::predict.averaging()结合glmmTMB模型与栅格栈进行预测？

核心问题梳理 #

使用glmmTMB构建含随机效应的混合模型，经MuMIn::dredge模型选择、model.avg模型平均后，调用terra::predict预测栅格栈时出现以下问题：

• 参数full、re.form未被识别
• 预测仍受随机效应水平约束
• 出现「无效因子水平」「对比矩阵行数错误」警告

问题根源与解决办法 #

1. 参数full被忽略 #

model.avg默认仅保留入选模型中出现的变量（子集模型集），此时predict时的full=TRUE参数无效。需在模型平均阶段就指定full=TRUE，确保保留全模型的所有变量：

avgm <- model.avg(confset.95p, full = TRUE)

2. 随机效应与因子水平不匹配 #

模型包含(1|Species)随机效应，即使设置re.form=~0，glmmTMB的预测函数仍要求输入数据的因子水平与原模型匹配。若要做边际预测（忽略随机效应），需：

• 在预测时传递re.form=~0到底层的glmmTMB::predict，同时设置allow.new.levels=TRUE（允许新的/缺失的因子水平）
• 确保预测数据的Species因子水平与原模型一致，或在处理栅格时将数值转换为对应因子

3. 栅格数据的因子类型问题 #

原数据中Species是因子，而生成的栅格spec是数值型，直接用于预测会导致因子转换错误。需将栅格转换为因子类型，并匹配原模型的因子水平：

# 原模型的Species因子水平
spec_levels <- levels(data$Species)
# 将数值栅格转为因子，指定水平
spec <- as.factor(spec, levels = spec_levels)

修改后的完整代码 #

library(glmmTMB)
library(MuMIn)
library(tidyverse)
library(terra)
data(iris)

# 数据预处理：创建变量、中心化/标准化
data <- iris %>%
  mutate(Species = as.factor(Species),
         Sepal.Length.2 = Sepal.Length^2,
         Sepal.Width.2 = Sepal.Width^2,
         Petal.Width.2 = Petal.Width^2,
         Sepal.Length.Log = log(Sepal.Length),
         Sepal.Width.Log = log(Sepal.Width),
         Petal.Width.Log = log(Petal.Width)) %>%
  mutate(across(c(-Species, -Petal.Length), ~c(scale(.x))))

# 构建全模型
mod <- glmmTMB(Petal.Length ~ Sepal.Length + Sepal.Length.2 + Sepal.Length.Log +
                 Sepal.Width + Sepal.Width.2 + Sepal.Width.Log +
                 Petal.Width + Petal.Width.2 + Petal.Width.Log +
                 (1|Species), data = data, na.action = 'na.fail')

# 模型选择
mmodels <- dredge(mod, m.lim = c(2,3), fixed = ~(1|Species))
confset.95p <- get.models(mmodels, subset = cumsum(weight) <= .95)

# 模型平均：指定full=TRUE保留全变量
avgm <- model.avg(confset.95p, full = TRUE)

# 数据框测试预测：正确传递参数到glmmTMB
test <- avgm$x %>%
  as.data.frame() %>%
  select(-1) %>%
  mutate(Species = NA)
# 使用参数传递，确保re.form生效
predict(avgm, test, type = 'link', se.fit = TRUE, 
        re.form = ~0, allow.new.levels = TRUE, full = TRUE)

# 生成栅格数据：修正Species因子类型
sw <- sw2 <- swl <- sl <- sl2 <- sll <- pw <- rast(ncol=10, nrow=10)
set.seed(4)
values(sw) <- runif(ncell(sw), min = data$Sepal.Width, max = data$Sepal.Width)
values(sw2) <- runif(ncell(sw2), min = data$Sepal.Width.2, max = data$Sepal.Width.2)
values(swl) <- runif(ncell(swl), min = data$Sepal.Width.Log, max = data$Sepal.Width.Log)
values(sl) <- runif(ncell(sl), min = data$Sepal.Length, max = data$Sepal.Length)
values(sl2) <- runif(ncell(sl2), min = data$Sepal.Length.2, max = data$Sepal.Length.2)
values(sll) <- runif(ncell(sll), min = data$Sepal.Length.Log, max = data$Sepal.Length.Log)
values(pw) <- runif(ncell(pw), min = data$Petal.Width, max = data$Petal.Width)

# 将Species栅格转为因子，匹配原模型水平
spec <- (pw*0) + 1
spec_levels <- levels(data$Species)
spec <- as.factor(spec, levels = spec_levels)

# 组合栅格栈并命名
x <- c(sw, sw2, swl, sl, sl2, sll, pw, spec)
names(x) <- c("Sepal.Width", "Sepal.Width.2", "Sepal.Width.Log",
              "Sepal.Length", "Sepal.Length.2", "Sepal.Length.Log",
              "Petal.Width", "Species")

# 栅格预测：自定义函数传递参数到predict.averaging，处理随机效应
out <- terra::predict(x, avgm, type = 'link', se.fit = TRUE, 
                      re.form = ~0, allow.new.levels = TRUE, full = TRUE,
                      fun = function(model, ...) {
                        pred <- MuMIn:::predict.averaging(model, ...)
                        # 返回拟合值或标准误，这里返回标准误
                        pred$se.fit
                      })

额外说明 #

• re.form=~0用于忽略随机效应，实现边际预测；若要包含随机效应，需确保栅格中Species的水平与原模型完全匹配
• allow.new.levels=TRUE允许预测数据中出现原模型未有的因子水平（包括NA），避免因子水平报错
• 模型平均时指定full=TRUE，确保预测时能处理全变量集，避免full参数被忽略

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