问题描述
我想在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()