Я хочу написать функцию, которая принимает данные и запускает мультиномиальную регрессию (используя 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) Accounting and correcting for sample (un)representativeness

Используя регрессию, я могу: (1) моделировать взаимосвязь между цветовыми предпочтениями и демографическими переменными и (2) прогнозировать «скорректированный» средний ответ на основе демографических значений, которые встречаются в популяции (но не обязательно в моей выборке). Поскольку моя переменная, представляющая интерес, является номинальной, я использую полиномиальную регрессию (с `nnet :: multinom`).

1. Fit the model

 library(nnet)

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

2. Define a vector with the «corrected» values as they happen to be in the population level, to use in the prediction step.

• age — I know that the average age in the population is 45.
• sex — I know that sex is roughly at 50% split, thus 0.5.
• color blindness — I know that on average, 2% of the population is color blinded (say). Hence 0.02.

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

3. Use a prediction function to get a focal prediction for each color, given the values in one_average_person .

I’ve found only effects::Effect to play well with a model generated from nnet::multinom . Still, since I couldn’t find a straightforward way to get a focal prediction for the values I specify, I ended up with a workaround. In the following code, age is the «focal» predictor, but I also specify the other variables using the given.values argument. Furthermore, I can’t just ask for age = 45 because Effect cannot take a single value, so I ask for a prediction for both age = 45 and age = 90 . Then I remove the prediction for 90 because I don’t need it.

 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)

  }
  

Есть идеи, как заставить эту функцию работать?

Вот версия вашей функции, которая работает, если вы предоставляете все имена переменных в виде строк:

 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 использовать this в качестве имени переменной в списке. Вы можете использовать := from rlang для назначения (принудительного) имени в качестве имени переменной списка, однако это работает не в обычном list , а в rlang::list2