如何使用基于nnet :: multinom模型的{ggeffects}获取预测概率图的置信区间？

如何解决如何使用基于nnet :: multinom模型的{ggeffects}获取预测概率图的置信区间？

我想在R中绘制拟合nnet::multinom()函数的多项式模型的预测概率。我有对数刻度上的数字预测变量。

即使{ggeffects}应该与multinom()兼容，该图也不像线性模型那样显示置信区间。

我不熟悉R和这个社区，所以我很抱歉，如果这个问题很基础或缺少必要的内容。这是一个小例子：

library(tidyverse)
library(nnet)
library(effects)
library(ggeffects)


df <- data.frame(response = c("1 Better","1 Better","2 Medium","3 Worse","3 Worse"),count = c(1000,2000,4000,6000,10000,3000,5000,11000))

mod1 <- multinom(response ~ log(count),data = df)
summary(mod1)

effects::effect(mod1,term="log(count)",se=TRUE,confidence.level=.95) %>% plot() # Produces CIs.

ggeffects::ggpredict(mod1,terms = "count") %>% plot() + theme_bw() # No confidence intervals.

如果其他人正在寻找{ggeffects}的替代品，我在寻找解决方案时尝试了以下方法：

使用effects::effect()：有效，包括置信区间，但是外观不是那么可定制。

结合{ggeffects}和{effects}：参见此post on R Studio Community，其中将效果包的置信区间与ggeffects组合在一起以创建图。我得到了错误

Error in FUN(X[[i]],...) : object 'L' not found

但这对那个人有用。

使用{MNLpred}软件包及其 mnl_pred_ova() ：对我来说不起作用，因为我的预测变量在对数范围内。我收到以下错误：

Error in eval(parse(text = paste0("data$",xvari))) : attempt to apply non-function

使用mnlAveEffPlot()中的{damisc}函数：可行，但是绘图不像我想要的那样可定制。

解决方法

您可以使用ggeffects::ggemmeans()进行此操作。

library(tidyverse)
library(ggthemes)
library(nnet)
library(ggeffects) # package version used: v0.16.0

df <- data.frame(response = c("1 Better","1 Better","2 Medium","3 Worse","3 Worse"),count = c(1000,2000,4000,6000,10000,3000,5000,11000))

mod1 <- multinom(response ~ log(count),data = df)

ggemmeans(mod1,terms = "count") %>% plot() + ggthemes::theme_tufte()

有关如何使用{ggeffects}的更多信息，您可能还想看看package documentation，尤其是ggemmeans()和ggpredict()等之间的区别（例如here）。

{ggeffects}软件包借鉴了{effects}创建的输出，但是，我相信这就是您要寻找的内容，它使使用标准ggplot命令自定义绘图变得更加容易。

MNLpred 包无法处理回归函数内的 log()，但在您预先计算对数尺度时可以使用。

# Packages
library(tidyverse)
library(nnet)
library(MASS)
library(MNLpred)
library(scales)
library(ggeffects)
library(ggthemes)

df <- data.frame(response = c("1 Better",data = df)
summary(mod1)

# Log-scaled
df$count_log <- log(df$count)

# Regression
mod2 <- multinom(response ~ count_log,data = df,Hess = TRUE)

# The models are identical:
coef(mod1) == coef(mod2)

在此步骤之后，您可以使用 mnl_pred_ova 或 mnl_fd2_ova 函数进行预测概率或一阶差分/预测边际效应。

# 10 steps for predictions
steps <- (max(df$count_log) - min(df$count_log))/9

pred1 <- mnl_pred_ova(mod2,by = steps,x = "count_log")

x_breaks <- seq(from = min(df$count_log),to = max(df$count_log),length.out = 5)

x_labels <- seq(from = min(df$count),to = max(df$count),length.out = 5)

pred1$plotdata %>%
  ggplot(aes(x = count_log,y = mean,ymin = lower,ymax = upper)) +
  facet_wrap(. ~ response) +
  geom_line() +
  geom_ribbon(alpha = 0.2) +
  scale_y_continuous(labels = percent_format()) +
  scale_x_continuous(breaks = x_breaks,labels = x_labels) +
  theme_bw()

或预测的边际效应：

pred_fd <- mnl_fd2_ova(model = mod2,x = "count_log",value1 = min(df$count_log),value2 = max(df$count_log),data = df)

pred_fd$plotdata_fd %>%
  ggplot(aes(x = categories,ymax = upper)) +
  geom_pointrange() +
  scale_y_continuous(labels = percent_format()) +
  labs(title = "Predicted effect of Count on responses",x = "Categories",y = "Predicted marginal effect") +
  theme_bw()