如何从多项回归模型中获得多个预测而不是焦点预测（即，按因子变量拆分）

Question

我想运行多项式回归，以获取每个选择关闭问题的平均频率，除以一个因素（性别：男性/女性）。

背景

我想比较 4 种奶酪来衡量每种奶酪的受欢迎程度，有 4 种可能性：切达干酪、马苏里拉干酪、豪达干酪和布里干酪。我出去向 200 个人要了他们喜欢的奶酪。每个人只从 4 种类型中选择一种。我最终还收集了一些人口统计信息，包括性别、年龄和体重。

完成数据收集后，我想查看每种奶酪类型的受欢迎程度（加起来为 100%）。由于我想控制gender、age 和weight，我认为这里适合使用多项式回归。

但我也很想看看男性和女性之间的结果有何不同，我想在我的模型中加入 gender 作为一个因素。如何根据我的（多项式）模型生成双重预测，分别获得女性和男性的预测值，以便我可以在两个性别水平之间进行比较？

数据

library(truncnorm)
library(tidyverse)

set.seed(999)

cheese_df <-
  tibble(
    age = round(rtruncnorm(
      n = 200,
      a = 20,
      b = 80,
      mean = 25,
      sd = 25.09
    )),
    cheese_response = as_factor(sample(
      c("cheddar", "mozzarella", "gouda", "brie"),
      size = 200,
      replace = TRUE
    )),
    gender = sample(c(0, 1), size = 200, replace = TRUE),
    weight = rtruncnorm(
      n = 200,
      a = 40,
      b = 120,
      mean = 70,
      sd = 25.09
    )
  )


> cheese_df

## # A tibble: 200 x 4
##      age cheese_response gender weight
##    <dbl> <fct>            <dbl>  <dbl>
##  1    45 cheddar              0   62.2
##  2    32 cheddar              0   45.0
##  3    58 cheddar              1   87.6
##  4    28 brie                 0   68.8
##  5    49 gouda                0   88.2
##  6    29 brie                 1   74.5
##  7    49 cheddar              0   74.0
##  8    27 gouda                1   90.3
##  9    28 brie                 0   56.5
## 10    48 mozzarella           0   72.9
## # ... with 190 more rows

如果我只想运行多项回归并控制年龄、性别和体重**而不按性别划分**，我可以这样做：

library(nnet)
library(effects)


fit <- nnet::multinom(cheese_response ~ age + gender + weight, data = cheese_df)

average_person_for_control <-
  c(
    age = 50,
    gender = 0.5,
    weight = 75
  )

prediction <-
  effects::Effect("age",
                  fit,
                  given.values = average_person_for_control,
                  xlevels = list(age =
                                   c(45, 90)))


proportions_for_plot <-
  data.frame(prediction$prob, prediction$lower.prob, prediction$upper.prob) %>% 
  slice(1) %>%
  pivot_longer(., cols = everything(), 
               names_to = c(".value", "response"), 
               names_pattern = "(.*)\\.(.*$)") %>%
  rename("lower_ci" = "L.prob",
         "upper_ci" = "U.prob",
         "estimate" = "prob")


ggplot(proportions_for_plot, aes(x = reorder(response, -estimate), y = estimate)) +
  geom_bar(stat = "identity", width = 0.7, fill = "darkgreen") +
  geom_errorbar(aes(ymin = lower_ci, ymax = upper_ci),
                width = 0.2) +
  geom_text(aes(label = paste0(100*round(estimate,2), "%")),
            vjust = 1.6, 
            color = "white", size = 3) +
  xlab("cheese type") +
  ylab("proportion of people choosing this type")

但是，我有兴趣生成相同的条形图，只是它会拆分男性和女性的条形

这就是我想要得到的情节

（忽略此演示中的值）

一种方法是按性别对数据进行子集化，在每个子集上运行相同的模型，生成两个条形图并将它们合并。但是，我想在模型中加入gender 作为一个因素，然后才输出拆分条形图。这部分得到了解决，因为 gender 已经是模型的一部分： fit <- nnet::multinom(cheese_response ~ age + gender + weight, data = cheese_df).

仍然，至于将预测按性别分开，以便在条形图中并排比较它们，我遇到了麻烦。这是因为effects::Effect() 只接受一个向量到它的given.values 参数中。否则，我会执行以下操作来提供预测（就像我使用 predict 时会做的那样）：

control_by_gender <-
  expand.grid(
    age = 50,
    weight = 75,
    gender = c(0, 1)
  )

> control_by_gender

##   age weight gender
## 1  50     75      0
## 2  50     75      1

如上所示，在处理多项式模型对象时，我知道如何获得这样的多重预测——而不是焦点预测吗？我的最终目标是按性别划分的条形图，就像上面的演示一样。我一直在使用Effects::effect 来生成预测，但我愿意接受任何可以做多重预测技巧的替代品。

Answer 1

为什么不直接将级别应用到effects::Effect 调用中

prediction <- do.call(rbind,lapply(0:1, function(x) {
    eff <- effects::Effect("age",
                  fit,
                  given.values =c(age = 50,
                        weight = 75,
                        gender = x),
                  xlevels = list(age =c(45, 90)))
    data.frame(level=x, eff$prob, eff$lower.prob, eff$upper.prob) %>% slice(1)
    }))


proportions_for_plot <-
  prediction %>% 
  pivot_longer(., cols = !level, 
               names_to = c(".value", "response"), 
               names_pattern = "(.*)\\.(.*$)") %>%
  rename("lower_ci" = "L.prob",
         "upper_ci" = "U.prob",
         "estimate" = "prob")


ggplot(proportions_for_plot, aes(x = as.factor(response), y = estimate, fill=factor(level))) +
  geom_bar(stat = "identity", width = 0.7,position="dodge") +
  geom_errorbar(aes(ymin = lower_ci, ymax = upper_ci), position=position_dodge(.9),
                width = 0.2) +
  geom_text(aes(label = paste0(100*round(estimate,2), "%")),
            vjust = 1.6, 
            color = "white", size = 3, position=position_dodge(.9)) +
  xlab("cheese type") +
  ylab("proportion of people choosing this type")

Answer 2

这个答案使用与@Abdessabour Mtk 相同的直觉，只是用purrr::map 代替并进行了一些重构：

make_eff_df <- function(gender, fit) {
  Effect("age", fit, xlevels = list(age = c(45, 90)),
         given.values = c(age = 50, weight = 75, gender = gender)) %>%
    as_tibble() %>%
    mutate(gender = gender) %>%
    select(gender, matches("[a-z\\.]?prob")) %>%
    slice(1)
}


map_dfr(0:1, make_eff_df, fit) %>% 
  pivot_longer(-gender, names_to = c(".value", "response"), 
               names_pattern = "(.+)\\.(.+$)") %>%
  rename(lower_ci = "L.prob", upper_ci = "U.prob", estimate = "prob") %>%
  mutate(across(1:2, as.factor)) %>%
  ggplot(aes(x = reorder(response, -estimate), y = estimate, fill = gender)) +
  geom_bar(stat = "identity", width = 0.7, position = position_dodge(.9)) +
  geom_errorbar(aes(ymin = lower_ci, ymax = upper_ci), 
                position = position_dodge(.9),
                width = 0.2) +
  geom_text(aes(label = scales::percent(estimate, accuracy = 1)),
            vjust = 1.6, color = "white", size = 3, position=position_dodge(.9)) +
  labs(x = "cheese type",
       y = "proportion of people choosing this type")