Я хочу написать функцию, которая принимает данные и запускает полиномиальную регрессию (используя 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) Учет и корректировка выборки (не) репрезентативности
Using regression I can: (1) model the relationship between color preference and demographics variables, and (2) predict a "corrected" average response based on the demographic values that occur in the population (but not necessarily in my sample). Since my variable of interest is nominal, I'm using a multinomial regression (with `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. На самом деле это далеко не идеально, потому что не всегда в моих данных будет возраст. Но моя функция не работает, несмотря на это. Причина этой трудности в том, что возраст вводится как строка в focal.predictors аргументе, но как переменная в xlevels (в списке). Я пробовал использовать двойные фигурные скобки (оценки приборки), но все равно безуспешно.
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)
}
Есть идеи, как заставить эту функцию работать?
