如何解决如何编写自定义函数以从`effects :: Effect`中提取预测 A简单但可能不准确的方法B核算和更正示例不具有代表性
我想编写一个函数，该函数接受数据并运行多项式回归（使用nnet::multinom），然后提取焦点预测（使用Effects::effect）。虽然我可以使用常规代码完成此操作，但是自定义功能失败。
示例
背景
我进行了一项研究，以找出人们最喜欢的颜色类型：红色，绿色或蓝色。我对 200 人进行抽样，并要求他们选择最喜欢的一种颜色。因为我怀疑某些变量可能会使结果混淆，所以我也对它们进行了测量：（1）性别，（2）色盲和（3）年龄。
方法
我将使用nnet::multinom进行多项式回归，然后从该模型（使用Effects::effect提取 focus 预测），该预测将解释特定的性别，色盲和年龄的值。
数据
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,is_colorblind = sample(c(0,1),prob = c(0.2,0.8),is_female = sample(c(0,prob = c(0.3,0.7),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岁。
- 性别-我知道性别之间大约有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构建了函数。实际上，这远非理想，因为我的数据并不总是有年龄的。但是无论如何，我的功能无法正常工作。造成这种困难的原因是，“ 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,xlevels = list(age_var_for_Effect = c(ave_age,90)
)
)
return(prediction)
}
关于使此功能正常工作的任何想法？
解决方法
这是函数的一个版本，如果您将所有变量名都提供为字符串，那么该函数将起作用：
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,is_colorblind = sample(c(0,1),prob = c(0.2,0.8),is_female = sample(c(0,prob = c(0.3,0.7),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中则有效
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。
