如何编写自定义函数以从 `effects::Effect()` 中提取预测
How to write a custom function for extracting predictions from `effects::Effect()`
我想编写一个函数来接收数据和 运行 多项式回归（使用 nnet::multinom），然后提取焦点预测（使用 Effects::effect）。虽然我能够使用常规代码完成它，但自定义函数失败了。
示例
背景
我进行了一项研究，以确定人们最喜欢哪种颜色：红色、绿色或蓝色。我对 200 个人进行了抽样，并要求他们选择他们最喜欢的一种颜色。因为我怀疑某些变量可能会混淆结果，所以我也测量了它们：(1) sex、(2) color blindness，以及 (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
由于我想找到对所有人都通用的见解，所以我不太相信我得到的 table 的准确性。这是因为我的样本不具有代表性。在我的样本中，20% 的人是色盲，70% 是女性。如果我有理由相信性别和色盲可能会影响颜色流行度，那么这个样本就有问题。
(B) 样本（非）代表性的核算和更正
使用回归，我可以：(1) 为颜色偏好和人口统计变量之间的关系建模，以及 (2) 根据人口中出现的人口统计值（但不一定在我的样本中）预测 "corrected" 平均响应。由于我感兴趣的变量是名义变量，因此我使用了多项式回归（使用 `nnet::multinom`）。
1.拟合模型
library(nnet)
fit <-
nnet::multinom(chosen_color ~ age + is_colorblind + is_female,
data = df)
2。定义一个带有“校正”值的向量，因为它们恰好处于人口水平，用于预测步骤。
- age -- 我知道人口的平均年龄是 45.
- sex -- 我知道 sex 大约是 50% split，因此是 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
table中的值反映了每种颜色的受欢迎程度，同时考虑了人口级别的情况。
写一个函数来简化上面的回归+预测过程
虽然我不得不用 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)
}
关于使此功能起作用的任何想法？
如果您以字符串形式提供所有变量名称，这是您的函数的一个版本：
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
我想编写一个函数来接收数据和 运行 多项式回归（使用 nnet::multinom），然后提取焦点预测（使用 Effects::effect）。虽然我能够使用常规代码完成它，但自定义函数失败了。
示例
背景
我进行了一项研究，以确定人们最喜欢哪种颜色：红色、绿色或蓝色。我对 200 个人进行了抽样，并要求他们选择他们最喜欢的一种颜色。因为我怀疑某些变量可能会混淆结果，所以我也测量了它们：(1) sex、(2) color blindness，以及 (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
由于我想找到对所有人都通用的见解，所以我不太相信我得到的 table 的准确性。这是因为我的样本不具有代表性。在我的样本中，20% 的人是色盲，70% 是女性。如果我有理由相信性别和色盲可能会影响颜色流行度，那么这个样本就有问题。
(B) 样本（非）代表性的核算和更正
使用回归，我可以：(1) 为颜色偏好和人口统计变量之间的关系建模，以及 (2) 根据人口中出现的人口统计值（但不一定在我的样本中）预测 "corrected" 平均响应。由于我感兴趣的变量是名义变量，因此我使用了多项式回归（使用 `nnet::multinom`）。1.拟合模型
library(nnet)
fit <-
nnet::multinom(chosen_color ~ age + is_colorblind + is_female,
data = df)
2。定义一个带有“校正”值的向量，因为它们恰好处于人口水平，用于预测步骤。
- age -- 我知道人口的平均年龄是 45.
- sex -- 我知道 sex 大约是 50% split，因此是 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
table中的值反映了每种颜色的受欢迎程度，同时考虑了人口级别的情况。
写一个函数来简化上面的回归+预测过程
虽然我不得不用 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)
}
关于使此功能起作用的任何想法？
如果您以字符串形式提供所有变量名称，这是您的函数的一个版本：
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