0

I am having an issue with tidymodels that I can't seem to figure out. Not sure if this is the intended behavior or an issue, but either way I would appreciate any help!

I am building a logistic regression prediction model with a two-level factor as the outcome, and per tidymodels convention have set the "positive class" as the first level.

The base R stats::glm() assumes exactly the opposite: that the "positive class" is the second level, and the "reference" is the first level.

With that in mind, I anticipated that fitting a model with a tidymodels workflow vs. stats::glm() would result in estimated coefficients with similar magnitude and opposite directions. However, it seems that in reality, tidymodels is behaving as stats::glm() and treating the second level as the positive class.

library(tidymodels)

#build model to predict "manual" (am == 1)
#Positive class is first level of factors per tidymodels convention
df <- 
  mtcars %>% 
  as_tibble() %>% 
  mutate(am = factor(am, levels = c("1", "0")))

#tidymodels
recipe <- recipe(df) %>% 
  update_role("am", new_role = "outcome") %>% 
  update_role("mpg", new_role = "predictor")

glm_model <- 
  logistic_reg() %>% 
  set_engine("glm") %>% 
  set_mode("classification")

glm_wf <- 
  workflow() %>% 
  add_recipe(recipe) %>% 
  add_model(glm_model)

glm_fit <-
  glm_wf %>%
  fit(df)

glm_fit %>%
  extract_fit_parsnip() %>%
  tidy(exponentiate = T)

# A tibble: 2 × 5
  term        estimate std.error statistic p.value
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1 (Intercept)  738.        2.35       2.81 0.00498
2 mpg            0.736     0.115     -2.67 0.00751

#base R
glm(am ~ 
      mpg,
    family = "binomial",
    data = df) %>% 
  tidy(exponentiate = T)

# A tibble: 2 × 5
  term        estimate std.error statistic p.value
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1 (Intercept)  738.        2.35       2.81 0.00498
2 mpg            0.736     0.115     -2.67 0.00751

#base R (treats second level as positive class) and tidymodels (treats first level as positive class) have the same output!

Any ideas? This is causing a lot of havoc when I try to report ORs and then use yardstick for performance assessment (yardstick assumes positive class is first). Thanks so much for the help, loving tidymodels overall.

out_of_step
  • 1
  • 1
  • 2
    In base R, the first factor level is assumed to be the event. Using the yardstick package, you have the [flexibility to select which event level](https://stackoverflow.com/questions/63211597/how-do-i-set-which-level-is-the-event-in-my-outcome-variable-using-tidymodels) you wish to evaluate. – Desmond May 05 '22 at 03:11
  • 1
    I understand that I can manually choose to treat the second level as the positive class in yardstick, but why would parsnip and yardstick treat factor levels differently? It’s not consistent - would require a lot of manual adjustment of metric specifications when really I just want to use the default and be consistent within tidymodels. – out_of_step May 05 '22 at 12:07
  • How is it not consistent? From the post I linked, both base R and tidymodels treats the first level as the positive class. By default, yardstick functions [generally defaults to first level too](https://yardstick.tidymodels.org/reference/roc_curve.html#:~:text=The%20default%20uses%20an%20internal%20helper%20that%20generally%20defaults%20to%20%22first%22) – Desmond May 06 '22 at 02:58
  • In the documentation to `glm` it states: *For binomial and quasibinomial families the response can also be specified as a factor (when the first level denotes failure and all others success)*. `glm` will model the second level as the positive class or event (not first). `yardstick` evaluates the first level as the positive class. In order to get coefficients, predicted probabilities, and performance metrics for the positive class that make sense, you have to modify the default of `yardstick`, or transform the coefficients from `glm`. – out_of_step May 06 '22 at 12:31
  • Where do you see that it says that parsnip sets the first class as the "positive class", I fail to find that reference. – EmilHvitfeldt May 06 '22 at 16:51
  • 1
    `parsnip` does not state this explicitly, but maybe it should? I think I am confused because it seems `parsnip` and `yardstick` seem assume different levels to be the positive class, yet both exist within the `tidymodels` meta-package. – out_of_step May 06 '22 at 17:05
  • 1
    Refer to the link in my first answer where Julia Silge states "You don't need to set which level of your outcome variable is the "event" until it is time to evaluate your model.". I also found some [related reading](https://www.tmwr.org/performance.html#:~:text=There%20is%20some%20heterogeneity%20in%20R%20functions) that may be helpful for you. – Desmond May 12 '22 at 01:43

0 Answers0