如何编写自定义函数以从 `effects::Effect()` 中提取预测

Posted 2023-03-12

【中文标题】如何编写自定义函数以从 `effects::Effect()` 中提取预测

【英文标题】：How to write a custom function for extracting predictions from `effects::Effect()`

【发布时间】：2021-01-25 09:00:41

【问题描述】：

我想编写一个函数，它接收数据并运行多项回归（使用nnet::multinom），然后提取焦点预测（使用Effects::effect）。虽然我可以使用常规代码完成它，但自定义函数失败了。

示例

背景

我进行了一项研究，以找出人们最喜欢哪种颜色：红色、绿色或蓝色。我对 200 个人进行抽样，并要求他们选择他们最喜欢的一种颜色。因为我怀疑某些变量可能会混淆结果，所以我也对它们进行了测量：(1) 性别、(2) 色盲和(3) 年龄。

方法

我将使用nnet::multinom 运行多项回归，然后从该模型中提取一个焦点预测（使用Effects::effect），这将解释性别的特定值、色盲和年龄。

数据

library(tidyverse)

set.seed(2020)

df <-
  data.frame(person_id = 1:200,
             chosen_color = sample(c("red", "green", "blue"), size = 200, replace = TRUE),
             age = sample(18:80, size = 200, replace = TRUE),
             is_colorblind = sample(c(0, 1), prob = c(0.2, 0.8), size = 200, replace = TRUE),
             is_female = sample(c(0, 1), prob = c(0.3, 0.7), size = 200, replace = TRUE)
           )

as_tibble(df)

## # A tibble: 200 x 5
##    person_id chosen_color   age is_colorblind is_female
##        <int> <chr>        <int>         <dbl>     <dbl>
##  1         1 blue            57             1         0
##  2         2 blue            51             1         0
##  3         3 blue            38             1         1
##  4         4 red             30             1         1
##  5         5 green           78             1         1
##  6         6 red             72             1         0
##  7         7 green           63             1         1
##  8         8 green           69             0         0
##  9         9 red             57             1         0
## 10        10 blue            20             0         1
## # ... with 190 more rows

每种颜色的流行度比例是多少？

(A) 简单但可能不准确的方法

只要找到chosen color中出现频率最高的颜色：

df %>%
  group_by(chosen_color) %>%
  summarise(n = n()) %>%
  mutate(freq = n / sum(n))

## # A tibble: 3 x 3
##   chosen_color     n  freq
##   <chr>        <int> <dbl>
## 1 blue            76  0.38
## 2 green           60  0.3 
## 3 red             64  0.32

由于我想找到对整个人群普遍的见解，因此我对所获得的表格的准确性几乎没有信心。这是因为我的样本不具有代表性。在我的样本中，20% 的人是色盲，70% 是女性。如果我有理由相信性别和色盲可能会影响颜色流行度，那么这个样本就有问题。

(B) 样本（不）代表性的会计和更正

使用回归我可以：（1）对颜色偏好和人口统计变量之间的关系进行建模，以及（2）根据人口中出现的人口统计值（但不一定在我的样本中）预测“校正的”平均响应。由于我感兴趣的变量是名义变量，因此我使用多项回归（使用 `nnet::multinom`）。

1.拟合模型

library(nnet)

fit <-
  nnet::multinom(chosen_color ~ age + is_colorblind + is_female,
                 data = df)

2.使用恰好在总体级别中的“校正”值定义一个向量，以用于预测步骤。

年龄

-- 我知道人口的平均年龄是 45 岁。

sex

-- 我知道性别大约是 50%，因此是 0.5。

色盲

-- 我知道平均有 2% 的人口是色盲（比如说）。因此为 0.02。

one_average_person <- 
  c(age = 45,
    is_female = 0.5,
    is_colorblind = 0.02
  )

3.给定one_average_person 中的值，使用预测函数获取每种颜色的焦点预测。

我发现只有effects::Effect 可以很好地与nnet::multinom 生成的模型配合使用。尽管如此，由于我找不到一种直接的方法来获得我指定的值的焦点预测，所以我最终找到了一种解决方法。在下面的代码中，age 是“焦点”预测​​器，但我还使用 given.values 参数指定了其他变量。此外，我不能只要求age = 45，因为Effect 不能采用单个值，所以我要求对age = 45 和age = 90 进行预测。然后我删除了90 的预测，因为我不需要它。

library(effects)

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


wrangled_prediction_data <-
  data.frame(prediction$prob, prediction$lower.prob, prediction$upper.prob) %>% 
  slice(1) %>%  ## <----- here I remove the unnecessary prediction for age = 90
  pivot_longer(., cols = everything(), 
               names_to = c(".value", "response"), 
               names_pattern = "(.*)\\.(.*$)") %>%
  rename("lower_ci" = "L.prob",
         "upper_ci" = "U.prob",
         "estimate" = "prob")


> wrangled_prediction_data

## # A tibble: 3 x 4
##   response estimate lower_ci upper_ci
##   <chr>       <dbl>    <dbl>    <dbl>
## 1 blue        0.474    0.328    0.625
## 2 green       0.290    0.172    0.445
## 3 red         0.236    0.129    0.391

表中的值反映了每种颜色的流行度，考虑到人口水平的情况。

编写一个函数来简化上面的回归+预测过程

虽然我不得不用Effect 做一些体操来获得我需要的东西（如果你看到比我笨拙的代码更好的方法，请提供反馈），我想编写一个函数来使这项工作更简洁。

我不成功的功能

如您所见，我仅限于使用age 作为预测器，所以我最终围绕age 构建了函数。实际上，这远非理想，因为我的数据中并不总是有年龄。但无论如何，我的功能都不起作用。造成这种困难的原因是“年龄”在focal.predictors 参数中作为字符串输入，但在xlevels 中作为变量输入（在列表中）。我尝试使用双花括号 (of tidy evaluation)，但仍然不成功。

require(dplyr)
require(nnet)
require(effects)

analyze_multiple_choice_w_age <-
  function(data,
           vars_demog,
           vars_dv,
           age_var_for_Effect,
           ave_age,
           one_ave_person_vec) 
    fit <-
      data %>%
      nnet::multinom(
        data = .,
        formula = as.formula(
        paste(
          vars_dv,
          paste(names(select( data , vars_demog )), collapse = " + "),
          sep = " ~ "
        )) 
        )
    
    prediction <-
      effects::Effect(
        focal.predictors = age_var_for_Effect,
        mod = fit,
        given.values = one_average_person,
        xlevels = list(age_var_for_Effect = c(ave_age, 90)
        )
      )
    
    return(prediction)

  

有什么想法可以让这个功能发挥作用吗？

【参考方案1】：

如果您将所有变量名都提供为字符串，则以下是您的函数的一个版本：

set.seed(2020)

df <-
  data.frame(person_id = 1:200,
             chosen_color = sample(c("red", "green", "blue"), size = 200, replace = TRUE),
             age = sample(18:80, size = 200, replace = TRUE),
             is_colorblind = sample(c(0, 1), prob = c(0.2, 0.8), size = 200, replace = TRUE),
             is_female = sample(c(0, 1), prob = c(0.3, 0.7), size = 200, replace = TRUE)
  )

require(dplyr)
require(nnet)
require(effects)
library(rlang)

analyze_multiple_choice_w_age <-
  function(data,
           vars_demog,
           vars_dv,
           age_var_for_Effect,
           ave_age,
           one_ave_person_vec) 
    fit <-
      data %>%
      nnet::multinom(
        data = .,
        formula = as.formula(
          paste(
            vars_dv,
            paste(vars_demog, collapse = " + "),
            sep = " ~ "
          )) 
      )
    
    prediction <-
      effects::Effect(
        focal.predictors = age_var_for_Effect,
        mod = fit,
        given.values = one_ave_person_vec,
        xlevels = list2(!!age_var_for_Effect := c(ave_age, 90)
        )
      )
    
    return(prediction)
    
  

test <- analyze_multiple_choice_w_age(
  data = df,
  vars_demog = c("age", "is_colorblind", "is_female"),
  vars_dv = "chosen_color",
  age_var_for_Effect = "age",
  ave_age = 45,
  one_ave_person_vec = c(age = 45,
                         is_female = 0.5,
                         is_colorblind = 0.02
  )
)


test

age effect (probability) for blue
age
       45        90 
0.3030466 0.2604459 

age effect (probability) for green
age
       45        90 
0.3992617 0.5270109 

age effect (probability) for red
age
       45        90 
0.2976917 0.2125432 

我改变了什么：

as.formula 可以直接处理字符串，所以我简化了这个 来自rlang，我使用!! 来强制评估age_var_for_Effect 将其用作列表中的变量名。您可以使用rlang 中的:= 将（强制）名称指定为列表的变量名称，但这在普通list 中不起作用，但在rlang::list2 中不起作用