此帖子的数据可用here和R脚本，而数据可用的here（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)
Effect plot 1Effect plot 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()图表的一种方式。在进行模型预测时，我使用了数字预测变量的平均值，因为这是在调用another post时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))
Model averaged prediction plot
为什么最后一个图与上面生成的两个effect()图如此不同？我期望他们会非常相似。例如，在两个effect()图中，在RHDV2夏季和冬季到来后，igm抗体存在的预测概率要低得多，但是在使用平均模型从predict()中生成的最后一个图中，预测的概率是igDV抗体的存在在RHDV2的夏季到来较高，在冬季对于RHDV2的到来和之后都相似。
我注意到有类似的帖子here，但这并没有帮助我解决我的问题。
对于那些可能感兴趣的人，我想出了解决问题的方法。 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。