Signs of logistic regression coefficients

[pfc5098]

It seems that, using tidy(logistic_fit), 1/[1 + exp(-Xb)] generates (1-p), not p.

[topepo]

Can you please provide a minimal reprex (reproducible example)? The goal of a reprex is to make it as easy as possible for me to recreate your problem so that I can fix it: please help me help you!

If you've never heard of a reprex before, start by reading "What is a reprex", and follow the advice further down that page.

[pfc5098]

# An example based on the titanic data set

library(titanic)
library(tidyverse)
titanic_df <- titanic_train %>% mutate(Survived = as.factor(Survived))
titanic_df %>% glimpse()
#> Rows: 891
#> Columns: 12
#> $ PassengerId <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17…
#> $ Survived    <fct> 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, …
#> $ Pclass      <int> 3, 1, 3, 1, 3, 3, 1, 3, 3, 2, 3, 1, 3, 3, 3, 2, 3, 2, 3, …
#> $ Name        <chr> "Braund, Mr. Owen Harris", "Cumings, Mrs. John Bradley (F…
#> $ Sex         <chr> "male", "female", "female", "female", "male", "male", "ma…
#> $ Age         <dbl> 22, 38, 26, 35, 35, NA, 54, 2, 27, 14, 4, 58, 20, 39, 14,…
#> $ SibSp       <int> 1, 1, 0, 1, 0, 0, 0, 3, 0, 1, 1, 0, 0, 1, 0, 0, 4, 0, 1, …
#> $ Parch       <int> 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 1, 0, 0, 5, 0, 0, 1, 0, 0, …
#> $ Ticket      <chr> "A/5 21171", "PC 17599", "STON/O2. 3101282", "113803", "3…
#> $ Fare        <dbl> 7.2500, 71.2833, 7.9250, 53.1000, 8.0500, 8.4583, 51.8625…
#> $ Cabin       <chr> "", "C85", "", "C123", "", "", "E46", "", "", "", "G6", "…
#> $ Embarked    <chr> "S", "C", "S", "S", "S", "Q", "S", "S", "S", "C", "S", "S…

library(tidymodels)
set.seed(123)
titanic_split <- initial_split(titanic_df, 0.75, strata = Survived)
titanic_training <- titanic_split %>% training()
titanic_testing <- titanic_split %>% testing()

logistic_model <- logistic_reg() %>% 
  set_engine("glm") %>% 
  set_mode("classification")
titanic_recipe <- recipe(
  Survived ~ Fare + Sex + SibSp + Parch,
  data = titanic_training
) %>% 
  step_corr(all_numeric(), threshold = 0.8) %>% 
  step_normalize(all_numeric()) %>% 
  step_dummy(all_nominal(), -all_outcomes())

titanic_recipe_prep <- titanic_recipe %>% prep(titanic_training)
titanic_training_prep <- titanic_recipe_prep %>% bake(titanic_training)
titanic_testing_prep <- titanic_recipe_prep %>% bake(titanic_testing)

logistic_fit <- logistic_model %>% 
  fit(Survived ~ ., titanic_training_prep)


# The following says males are more likely to survive
levels(titanic_df$Survived)
#> [1] "0" "1"
tidy(logistic_fit)
#> # A tibble: 5 x 5
#>   term        estimate std.error statistic  p.value
#>   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
#> 1 (Intercept)    1.18      0.163      7.29 3.13e-13
#> 2 Fare           0.691     0.148      4.66 3.15e- 6
#> 3 SibSp         -0.390     0.125     -3.11 1.89e- 3
#> 4 Parch         -0.170     0.106     -1.61 1.08e- 1
#> 5 Sex_male      -2.72      0.215    -12.6  1.12e-36

# But simple summary statistics says otherwise
titanic_training %>% 
  count(Sex, Survived) %>% 
  group_by(Sex) %>% 
  mutate(pct = n / sum(n)) %>% 
  ungroup()
#> # A tibble: 4 x 4
#>   Sex    Survived     n   pct
#>   <chr>  <fct>    <int> <dbl>
#> 1 female 0           60 0.252
#> 2 female 1          178 0.748
#> 3 male   0          352 0.817
#> 4 male   1           79 0.183


# A mannual check shows that
#  if we use the coefficients from tidy(logistic_fit)
#  to construct a logistic function
#  it outputs the predicted probabilties for
#  the 2nd level of factor "Survived"
b <- tidy(logistic_fit) %>% pull(estimate)
names(b) <- tidy(logistic_fit) %>% pull(term)
dim(b) <- c(5, 1)
X <- bind_cols(intercept = 1, titanic_training_prep) %>% 
  select(-Survived) %>%
  head() %>% 
  as.matrix() 

get_prob <- function(x) 1 / (1 + exp(-x))
get_prob(X %*% b)
#>           [,1]
#> [1,] 0.1246483
#> [2,] 0.8499279
#> [3,] 0.7514858
#> [4,] 0.8115059
#> [5,] 0.1667158
#> [6,] 0.1675730

logistic_fit %>% 
  predict(titanic_training_prep %>% head(), type = "prob")
#> # A tibble: 6 x 2
#>   .pred_0 .pred_1
#>     <dbl>   <dbl>
#> 1   0.875   0.125
#> 2   0.150   0.850
#> 3   0.249   0.751
#> 4   0.188   0.812
#> 5   0.833   0.167
#> 6   0.832   0.168

^{Created on 2021-05-14 by the reprex package (v1.0.0)}

[pfc5098]

Hi Max, please see the example above. I was confused by the model coefficient of "Sex_male" then realized what went wrong.

[hfrick]

Hi @pfc5098, thanks for providing a reprex! Yes, glm() models the success probability (success here being survival), see its documentation

For binomial and quasibinomial families the response can also be specified as a factor (when the first level denotes failure and all others success)

Therefore the tidy() method shows a negative coefficient, as returned by glm(). So I would say all is working as intended/documented, even if the glm() behaviour may have been a little surprising.

library(titanic)
library(tidymodels)
#> Registered S3 method overwritten by 'tune':
#>   method                   from   
#>   required_pkgs.model_spec parsnip

titanic_df <- titanic_train %>% mutate(Survived = as.factor(Survived))

set.seed(123)
titanic_split <- initial_split(titanic_df, 0.75, strata = Survived)
titanic_training <- titanic_split %>% training()

titanic_recipe <-
  recipe(
    Survived ~ Fare + Sex + SibSp + Parch,
    data = titanic_training
  ) %>%
  step_corr(all_numeric(), threshold = 0.8) %>%
  step_normalize(all_numeric()) %>%
  step_dummy(all_nominal_predictors())

titanic_recipe_prep <- titanic_recipe %>% prep(titanic_training)
titanic_training_prep <- titanic_recipe_prep %>% bake(titanic_training)


glm_mod <- glm(Survived ~ ., data = titanic_training_prep, family = binomial())
summary(glm_mod)
#>
#> Call:
#> glm(formula = Survived ~ ., family = binomial(), data = titanic_training_prep)
#>
#> Deviance Residuals:
#>     Min       1Q   Median       3Q      Max  
#> -2.1178  -0.6231  -0.6088   0.7432   2.1602  
#>
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)   1.1994     0.1654   7.252 4.11e-13 ***
#> Fare          0.6647     0.1513   4.394 1.11e-05 ***
#> SibSp        -0.3174     0.1264  -2.511    0.012 *  
#> Parch        -0.2113     0.1044  -2.024    0.043 *  
#> Sex_male     -2.7133     0.2181 -12.441  < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for binomial family taken to be 1)
#>
#>     Null deviance: 888.31  on 666  degrees of freedom
#> Residual deviance: 640.07  on 662  degrees of freedom
#> AIC: 650.07
#>
#> Number of Fisher Scoring iterations: 5

^{Created on 2021-05-17 by the reprex package (v2.0.0)}

[pfc5098]

Hi @hfrick,

Thanks for the response and it makes great sense.

So, for a 2-level logistic regression model, glm always models the probability of the 2nd level of the response factor, right?

How about other 2-level classification models that parsnip provides APIs for? It would be great if we can have consistent behaviors across all the models. This consistency will make people give more trust to real model parameters, instead of things like Shapley values (which gives good model interpretation no matter which factor level is being modeled but is computationally more expensive) ... or, I don't know, add a feature showing which level is being modeled along with the tidy() output?😂

[juliasilge]

In almost all cases, we don't change the default behavior like that of underlying models, so that you get the same result with the same data whether training a model via tidymodels or directly via the underlying functions.

We have support in yardstick, the package for model metrics and measuring model performance, for choosing a base level when computing metrics. For example, notice the event_level argument of roc_auc() and similar functions.

[pfc5098]

Got it 😃 I was thinking about something like:

tidy.parsnip_glm_logistic <- function(mod) {
glue("Modeling the probability of {levels(mod$response)[2]}")
tidy(mod$glm)
}

But yeah I can always manually check the default modeling level.

By the way, thanks for the amazing work you guys are doing!

[juliasilge]

Thank you so much for the kind words! 🙌

[github-actions]

This issue has been automatically locked. If you believe you have found a related problem, please file a new issue (with a reprex: https://reprex.tidyverse.org) and link to this issue.