在 glmer 中用预测概率和预测带绘制问题答案

【问题标题】：graphing issues in glmer with predicted probabilities and prediction bands在 glmer 中用预测概率和预测带绘制问题
【发布时间】：2020-11-13 20:06:49
【问题描述】：

我有一个在 glmer 中运行的模型，如下所示：

multi.sanctions.bust.full.ag <- glmer(allbuster ~ lageutradeshare100 + lagtradeopenP  + colonial  
                                    + lagtradesharePT + lnlaggdpp + lnlaggdpt  + duration + lndist + nobust + nobustsq + nobustcb + (1 | partnercode) + (1 | caseid),
                                      data=sanctions.data.new.scaled, family=binomial(link="logit"),
                              nAGQ=1,control=glmerControl(optimizer="nlminbwrap",optCtrl=list(maxfun=2e5)))

数据可以访问here。

我已经想通了，使用 R 中的 predictcommand 如何获得预测概率，使用下面的代码。据我所知，预测的概率是正确的：

tmpdat_intraeu <- multi.sanctions.bust.full.ag@frame[, c("caseid", "lndist", "lnlaggdpt", "duration",
"partnercode", "lnlaggdpp", "lagtradeopenP", 
"lageutradeshare100", "lagtradesharePT", "nobust", "nobustsq", 
"nobustcb", "colonial")]

jvalues_intraeu <- with(multi.sanctions.bust.full.ag, seq(from = 
min(multi.sanctions.bust.full.ag@frame[["lageutradeshare100"]]), 
to = max(multi.sanctions.bust.full.ag@frame[["lageutradeshare100"]]), 
length.out = 100))

pp_intraeu <- lapply(jvalues_intraeu, function(j) {
  tmpdat_intraeu$lageutradeshare100 <- j
  predict(multi.sanctions.bust.full.ag, newdata = tmpdat_intraeu, type = "response", re.form = NA)
})

# I don't think that the lines below this point are working for me; this is where I think the problem is:

plotdat_intraeu <- t(sapply(pp_intraeu, function(x) {
    c(M = mean(x), Med = median(x), quantile(x, c(0.25, 0.75), na.rm = TRUE), (mean(x)-(2*sd(x))),
      (mean(x)+(2*sd(x))))
}))

plotdat_intraeu <- as.data.frame(cbind(plotdat_intraeu, jvalues_intraeu))
colnames(plotdat_intraeu) <- c("PredictedProbabilityMean", "PredProbMedian", "quartile1", "quartile3", "lowersd", "uppersd", "lageutradeshare100")
head(plotdat_intraeu)
tail(plotdat_intraeu)

sb_intraeu <- ggplot() + geom_line(data=plotdat_intraeu, aes(x = lageutradeshare100, y = PredictedProbabilityMean), size = 2, color="blue") + 
  geom_ribbon(data=plotdat_intraeu, aes(x = lageutradeshare100, ymin = lowersd, ymax = uppersd),
              fill = "grey50", alpha=.5) +
  ylim(c(-.5, 1)) + 
  geom_hline(yintercept=0) +
    geom_rug(data=subset(multi.sanctions.bust.full.ag@frame,allbuster==0), aes(x=lageutradeshare100), color="black", size=1.0, sides="b", alpha= 3/4, length = unit(0.05, "npc")) +
    geom_rug(data=subset(multi.sanctions.bust.full.ag@frame,allbuster==1), aes(x=lageutradeshare100), color="red", size=1.0, sides="b", alpha = 1) +  
 theme(panel.grid.major = element_line(colour = "gray", linetype = "dotted"), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.title.y = element_text(size=12, face="bold"), axis.title.x = element_text(size=12, face="bold")) +
  xlab("Intra-EU Trade Share") + 
  ylab("Predicted Probability of Sanctions Busting") 

sb_intraeu

我的问题是结果图给了我这样的东西：在将我的论文发送给我的委员会审查时，其中一位教员告诉我，置信区间计算不正确而且太宽了。我同意评估，并且我已经看到在这些类型的模型中很难置信区间，但我不知道如何“修复”图表。我已经看到了使用 predictInterval 和 `bootMER``` 进行引导的建议，但我无法弄清楚如何让它们工作。

任何帮助将不胜感激。我的论文大部分都写好了，但在我更好地可视化我的密钥 IV 的效果之前，我无法提交。

【问题讨论】：

标签： r logistic-regression lme4 mixed-models

【解决方案1】：

我最近一直在使用可能有帮助的 merTools 包。可以在here找到有关如何使用该软件包的教程。

这个包做得很好的一件事是，它允许您指定要在它报告的置信区间中考虑哪些类型的随机性。我会在这里写更多关于它的内容，但是教程比我能更好地解释它，所以我建议只阅读教程的“不确定性”部分。

【讨论】：