【发布时间】：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)
}
有什么想法可以让这个功能发挥作用吗？
【问题讨论】：
标签： r function regression prediction multinomial