如何用marginaleffects实现仅截距多项模型的两两比较？

当使用nnet::multinom()拟合仅含截距的多项logit模型时，marginaleffects包的avg_comparisons()函数会因模型无有效预测变量报错，但emmeans可以通过emm_pairs()轻松完成分类概率的两两比较。以下是两种用marginaleffects实现相同效果的方法：

示例数据与模型 #

library(tidyverse)
library(nnet)
library(marginaleffects)
library(emmeans)
library(magrittr)

# 构造数据
d1 <- tibble(
  dv = letters[1:3] %>% 
    extract(1:3 %>% rep(c(600, 300, 100))) %>% 
    factor(letters[1:3])
)

# 拟合仅含截距的多项logit模型
mn1 <- multinom(dv ~ 1, data = d1)

用avg_predictions()可以正常获取各分类的预测概率：

mn1 %>% avg_predictions()

输出：

Group Estimate Std. Error    z Pr(>|z|)     S  2.5 % 97.5 %
 a      0.6    0.01549 38.7   <0.001   Inf 0.5696  0.630
 b      0.3    0.01449 20.7   <0.001 313.9 0.2716  0.328
 c      0.1    0.00949 10.5   <0.001  83.9 0.0814  0.119

但直接调用avg_comparisons()会报错：

mn1 %>% avg_comparisons()

错误信息：

Error: There is no valid predictor variable. Please change the `variables` 
argument or supply an alternative data frame to the `newdata` argument.

方法一：构造虚拟变量触发两两比较 #

通过给数据集添加一个取值完全相同的虚拟变量，重新拟合模型后，利用这个变量触发avg_comparisons()对因变量分类进行两两比较：

# 添加取值统一的虚拟变量
d1$dummy <- factor(rep("x", nrow(d1)))

# 重新拟合模型（虚拟变量不影响模型结果）
mn1_dummy <- multinom(dv ~ dummy, data = d1)

# 指定对dv的分类做两两比较
avg_comparisons(
  mn1_dummy,
  variables = list(dv = "pairwise"),
  by = NULL
)

输出结果与emmeans的pairs()结果一致：

contrast estimate std.error     z Pr(>|z|)  2.5 % 97.5 %
1      a-b      0.3    0.0285 10.54   <0.001 0.2441 0.3559
2      a-c      0.5    0.0212 23.57   <0.001 0.4584 0.5416
3      b-c      0.2    0.0190 10.54   <0.001 0.1628 0.2372

方法二：手动计算两两对比与统计量 #

先提取预测概率和模型的方差协方差矩阵，再手动构造对比矩阵计算差值、标准误和p值：

# 获取各分类的预测概率
preds <- predictions(mn1, type = "response")

# 提取模型的方差协方差矩阵
vcov_mat <- vcov(mn1)

# 构造两两对比矩阵
contrast_mat <- matrix(
  c(1, -1, 0,
    1, 0, -1,
    0, 1, -1),
  nrow = 3, byrow = TRUE,
  dimnames = list(c("a-b", "a-c", "b-c"), c("a", "b", "c"))
)

# 计算差值、标准误、z值和p值
comparisons <- tibble(
  contrast = rownames(contrast_mat),
  estimate = as.vector(contrast_mat %*% preds$estimate),
  std.error = sqrt(diag(contrast_mat %*% vcov_mat %*% t(contrast_mat)))
) %>%
  mutate(
    z = estimate / std.error,
    p.value = 2 * pnorm(-abs(z))
  )

print(comparisons)

输出：

# A tibble: 3 × 5
  contrast estimate std.error     z p.value
  <chr>       <dbl>     <dbl> <dbl>   <dbl>
1 a-b           0.3    0.0285  10.5  2.78e-26
2 a-c           0.5    0.0212  23.6  1.29e-122
3 b-c           0.2    0.0190  10.5  2.78e-26

内容的提问来源于stack exchange，提问作者tomw