r - 为什么 effect() 和 predict() 会产生不同的模型预测?
问题描述
这篇文章的数据可在此处获得,R 脚本和数据可在此处获得(R 脚本也在下面的帖子中)。提前感谢您的帮助。
我在glmmTMB
. 我最好的两个模型如下。
igm_20 <- glmmTMB(igm_pres ~ fRHDV2_arrive_cat + fseason + sage + save_ajust_abun + fseason*fRHDV2_arrive_cat + (1 | fsite), data = edit_pp_dat, family = binomial)
igm_21 <- glmmTMB(igm_pres ~ fRHDV2_arrive_cat + fseason + sage + save_ajust_abun + fseason*fRHDV2_arrive_cat + sage*fRHDV2_arrive_cat + (1 | fsite), data = edit_pp_dat, family = binomial)
我对交互特别感兴趣fseason*fRHDV2_arrive_cat
,因此在构建模型后,我创建了effect()
显示这种交互对我在两个模型中的结果变量的影响的图。
ef_1 <- effect(term = "fRHDV2_arrive_cat*fseason", mod = igm_20)
windows();plot(ef_1, xlab = "Season", ylab = "Predicted probability of IgM antibody presence", main = "", factor.names = FALSE)
ef_2 <- effect(term = "fRHDV2_arrive_cat*fseason", mod = igm_21)
windows();plot(ef_2, xlab = "Season", ylab = "Predicted probability of IgM antibody presence", main = "", factor.names = FALSE)
效果图 1 效果图 2 (很抱歉提供图的链接,我没有足够的声誉来发布实际的图)
从效果图中可以看出,交互的影响fseason*fRHDV2_arrive_cat
在两个模型中非常相似,这并不奇怪。然后我将这两个模型平均如下:
mod_ave_list_1 <- list(igm_20, igm_21)
mod_ave_1 <- model.avg(mod_ave_list_1, rank = AICc)
summary(mod_ave_1)
根据模型的平均结果,我尝试创建与effect()
上述类似的图。但是,由于该effect()
函数不适用于平均模型,并且模型没有实现产生总体平均模型预测的re.form = NA
能力,因此我首先必须在另一个包中重新创建和重新平均我的两个模型,如下所示:predict()
glmmTMB
predict_1 <- glmer(igm_pres ~ fRHDV2_arrive_cat + fseason + sage + save_ajust_abun + fseason*fRHDV2_arrive_cat + (1 | fsite), data = edit_pp_dat, family = binomial)
predict_2 <- glmer(igm_pres ~ fRHDV2_arrive_cat + fseason + sage + save_ajust_abun + fseason*fRHDV2_arrive_cat + sage*fRHDV2_arrive_cat + (1 | fsite), data = edit_pp_dat, family = binomial)
predict_list_1 <- list(predict_1, predict_2)
ave_predict <- model.avg(predict_list_1, rank = AICc)
然后,我创建了一个newdata
框架,我从中制作并绘制了模型预测,作为生成与effect()
上述类似的图的一种手段。在进行模型预测时,我使用了数值预测变量的平均值,因为这是另一篇文章建议在调用effect()
. 我包含re.form = NA
在predict()
函数中,以便我得到人口平均预测,因为我的模型包含随机效应。
a <- as.data.frame(c("Summer", "Autumn", "Winter", "Spring", "Summer", "Autumn", "Winter", "Spring"))
a$fRHDV2_arrive_cat <- c("Pre-RHDV2 arrival", "Pre-RHDV2 arrival", "Pre-RHDV2 arrival", "Pre-RHDV2 arrival", "Post-RHDV2 arrival", "Post-RHDV2 arrival", "Post-RHDV2 arrival", "Post-RHDV2 arrival")
mean(edit_pp_dat$sage, na.rm = TRUE) #4.659477e-17
mean(edit_pp_dat$save_ajust_abun, na.rm = TRUE) #-3.004684e-17
a$sage <- c(4.659477e-17, 4.659477e-17, 4.659477e-17, 4.659477e-17, 4.659477e-17, 4.659477e-17, 4.659477e-17, 4.659477e-17)
a$save_ajust_abun <- c(-3.004684e-17, -3.004684e-17, -3.004684e-17, -3.004684e-17, -3.004684e-17, -3.004684e-17, -3.004684e-17, -3.004684e-17)
a$fsite <- c(NA, NA, NA, NA, NA, NA, NA, NA)
colnames(a) <- c("fseason", "fRHDV2_arrive_cat", "sage", "save_ajust_abun", "fsite")
predict.values <- predict(ave_predict, backtransform = TRUE, newdata = a, se.fit = TRUE, re.form = NA)
a$estimates <- predict.values$fit
a$se <- predict.values$se.fit
a$lci <- a$estimates - 1.96*a$se
a$uci <- a$estimates + 1.96*a$se
a$fseason <- factor(a$fseason, levels = c("Summer", "Autumn", "Winter", "Spring"))
a$fRHDV2_arrive_cat <- factor(a$fRHDV2_arrive_cat, levels = c("Pre-RHDV2 arrival", "Post-RHDV2 arrival"))
ggplot(a, aes(x = fseason, y = estimates, colour = fRHDV2_arrive_cat, group = fRHDV2_arrive_cat)) + geom_line(size = 1) + geom_point(size = 3) + geom_errorbar(aes(ymin = lci, ymax = uci), width = .2) + labs(x = "Season", y = "Predicted probability of IgM seropositivity", colour = "RHDV2 arrival category") + scale_color_manual(labels = c("Pre-arrival", "Post-arrival"), values = c("red", "blue")) + theme(axis.title.x = element_text(face = "bold", size = 16), axis.title.y = element_text(face = "bold", size = 16), axis.text.x = element_text(size = 14), axis.text.y = element_text(size = 14), legend.title = element_text(face = "bold", size = 14), legend.text = element_text(size = 12))
为什么最后一个情节与effect()
上面制作的两个情节如此不同?我期待它们非常相似。例如,在这两个effect()
图中,在 RHDV2 到来后的夏季和冬季,igm 抗体存在的预测概率要低得多,但是在最后一个图中predict()
,使用平均模型,igm 抗体存在的预测概率在RHDV2 到达后的夏季和冬季 RHDV2 到达前和到达后的类似情况。
我注意到这里有一个类似的帖子,但这并没有帮助我解决我的问题。
解决方案
对于那些可能感兴趣的人,我想出了如何解决我的问题。edit_pp_dat$sage
和edit_pp_dat$save_ajust_abun
是标准化变量,因此它们的平均值为 0。因此,a$sage
应该a$save_ajust_abun
如下所示:
a$sage <- c(0, 0, 0, 0, 0, 0, 0, 0)
a$save_ajust_abun <- c(0, 0, 0, 0, 0, 0, 0, 0)
我在计算机上也遇到了困难,因为edit_pp_dat$sage
它是一个矩阵,根据提供给模型的数据是在矩阵还是数据框中edit_pp_dat$save_ajust_abun
,它的操作似乎有所不同。predict()
我不确定为什么mean(edit_pp_dat$sage, na.rm = TRUE)
,mean(edit_pp_dat$save_ajust_abun, na.rm = TRUE)
也不给 0。
推荐阅读
- kedro - 在 DataCatalog 中指定 Kedro 数据版本?
- ionic3 - 离子可选择打开时如何显示加载程序
- django - 如何在 Django 中过滤和显示“多对多”元素
- javascript - 节点 - 在后台不规则地运行脚本
- wireshark - 如何将 .pcapng 文件(wireshark)数据转换为 .bin 文件
- laravel - 如何在 laravel 8 中使用复合主键?
- javascript - Discord JS clearTimeout 未重置 setTimeout
- java - 尝试在 cmd 中运行 java 类时出错 - 找不到该类
- kotlin - 有没有一种简单的方法来优化 Kotlin 中的代码?
- module - 供应商的 Terraform 模块