You need to enable JavaScript to run this app.
最新活动
大模型
产品
解决方案
定价
生态与合作
支持与服务
开发者
了解我们

glmmTMB秩亏模型预测组均值/标准误及DHARMa、performance工具调用的错误问题

glmmTMB秩亏模型预测组均值/标准误及DHARMa、performance工具调用的错误问题

你好!我仔细梳理了你的问题和可复现代码,这个问题的核心是秩亏模型的参数化逻辑与预测时模型矩阵的匹配冲突,再加上部分工具包对glmmTMB秩亏模型的支持限制。下面分点给你具体的解决思路:

一、预测组均值/标准误报错的原因与解决方法

为什么会出现这个错误?

你的模型公式是count~time*lake(包含主效应和交互项),但你删除了time=Alake=ref1的所有观测,导致模型的设计矩阵出现秩亏(列线性相关)。虽然你用rank_check="adjust"让glmmTMB自动处理了秩亏(丢弃了冗余的系数列),但predict函数的校验逻辑很严格:它要求新数据的模型矩阵(由time*lake生成的列)必须和拟合时的模型矩阵列空间完全一致。

哪怕你的pred_data里的组合在原始数据中存在,秩亏后的参数化方式也会让某些组合的预测依赖于被丢弃的冗余系数,最终触发Prediction is not possible for unknown fixed effects的报错。

具体解决方法

方法1:用emmeans包(最推荐)

emmeans是专门处理边际/组均值计算的工具,能自动适配秩亏模型的参数化逻辑,不需要手动构建newdata,直接就能输出你需要的组均值、标准误和置信区间:

library(emmeans)

# 计算time*lake所有组合的预测均值,offset设为log(1)对应你要的min=1场景
emm_results <- emmeans(model, ~time*lake, offset = log(1), type = "response")

# 转成数据框方便后续处理
emm_df <- as.data.frame(emm_results)
head(emm_df)

输出结果里的response就是组均值,SE是标准误,还有现成的置信区间,完全避开了手动构建newdata的坑。

方法2:手动调整模型公式,消除秩亏

既然你明确知道缺失的是time=A+lake=ref1的组合,直接在模型公式里删除这个冗余的交互项,让模型成为满秩模型:

# 去掉缺失的交互项,模型设计矩阵变为满秩
model_fullrank <- glmmTMB(count~time*lake - timeA:lakeref1, 
                          family=nbinom1, 
                          offset=log(min), 
                          data=dat2)

之后再用你原来的predict代码就不会报错了,因为满秩模型的所有系数都有明确的估计值,新数据的组合都能找到对应的系数计算预测值。

方法3:限制newdata为拟合数据中实际存在的独立组合

如果你一定要用predict函数,可以从拟合数据dat2中提取所有不重复的time*lake组合作为新数据,确保它的模型矩阵完全在拟合时的列空间内:

# 提取dat2中所有存在的time和lake组合
pred_data_safe <- dplyr::distinct(dat2, time, lake)
# 设置min=1(对应你原来的需求)
pred_data_safe$min <- 1

# 现在运行预测就不会报错了
pred_data_safe$pred_Freq <- predict(model_fullrank, newdata = pred_data_safe, type = 'response', se.fit=FALSE)
pred_se <- predict(model_fullrank, newdata = pred_data_safe, type = 'response', se.fit=TRUE)$se.fit
pred_data_safe$pred_Freq_lwr <- pred_data_safe$pred_Freq - pred_se
pred_data_safe$pred_Freq_upr <- pred_data_safe$pred_Freq + pred_se

二、DHARMa模拟残差报错的解决方法

DHARMa的simulateResiduals函数内部会调用模型的预测逻辑,所以同样会遇到秩亏模型的参数化冲突。解决思路和上面一致:

  1. 优先用满秩模型:用方法2中调整后的model_fullrank来运行DHARMa:
res <- DHARMa::simulateResiduals(model_fullrank, n = 100)
plot(res) # 正常绘制残差诊断图
  1. 如果保留原秩亏模型:可以手动指定newdata为拟合数据dat2,避免函数自动生成包含冗余组合的新数据:
res <- DHARMa::simulateResiduals(model, n = 100, newdata = dat2)

三、performance包check_model报错的说明

这个错误的原因非常直接:目前performance包还没有实现对glmmTMB模型的check_model支持。你可以用以下替代方案做模型诊断:

  • 用DHARMa的残差诊断图(上面的plot(res))来检查模型拟合质量;
  • 用glmmTMB自带的summary(model)查看系数估计、标准误和显著性;
  • 用performance包的model_performance(model)来获取模型的拟合指标(这个函数是支持glmmTMB的):
library(performance)
model_performance(model)

备注:内容来源于stack exchange,提问作者E. McCallum

火山引擎 最新活动