Logistic regression may give a headache initially. While the structure and idea is the same as “normal” regression, the interpretation of the b’s (ie., the regression coefficients) can be more challenging.

This post provides a convenience function for converting the output of the glm function to a probability. Or more generally, to convert logits (that’s what spit out by glm) to a probabilty.

Note1: The objective of this post is to explain the mechanics of logits. There are more convenient tools out there. (Thanks to Jack’s comment who made me adding this note.)

Note2: I have corrected an error pointed out by Jana’s comment, below (you can always check older versions on the Github repo).

Example time

So, let’s look at an example. First load some data (package need be installed!):

data(titanic_train, package = "titanic")
d <- titanic_train  # less typing

Compute a simple glm:

glm1 <- glm(Survived ~ Pclass, data = d, family = "binomial")

The coeffients are the interesting thing:

coef(glm1)

## (Intercept)      Pclass 
##   1.4467895  -0.8501067

These coefficients are in a form called ‘logits’.

Takeaway

If coefficient (logit) is positive, the effect of this predictor (on survival rate) is positive and vice versa.

Here Pclass coefficient is negative indicating that the higher Pclass the lower is the probability of survival.

Conversion rule

To convert a logit (glm output) to probability, follow these 3 steps:

• Take glm output coefficient (logit)
• compute e-function on the logit using exp() “de-logarithimize” (you’ll get odds then)
• convert odds to probability using this formula prob = odds / (1 + odds). For example, say odds = 2/1, then probability is 2 / (1+2)= 2 / 3 (~.67)

R function to rule ‘em all (ahem, to convert logits to probability)

This function converts logits to probability.

logit2prob <- function(logit){
  odds <- exp(logit)
  prob <- odds / (1 + odds)
  return(prob)
}

For convenience, you can source the function like this:

source("https://sebastiansauer.github.io/Rcode/logit2prob.R")

For our glm:

logit2prob(coef(glm1))

## (Intercept)      Pclass 
##   0.8095038   0.2994105

How to convert logits to probability

How to interpret:

• The survival probability is 0.8095038 if Pclass were zero (intercept).
• However, you cannot just add the probability of, say Pclass == 1 to survival probability of PClass == 0 to get the survival chance of 1st class passengers.

Instead, consider that the logistic regression can be interpreted as a normal regression as long as you use logits. So the general regression formula applies as always:

y = intercept + b*x

That is, in our example

logits_survival = intercept + b_survival*Pclass

where b_survival is given in logits (it’s just the b-coefficient of Pclass).

So, it’ simple to calculate by hand, eg., the survival logits for a 2nd class passenger:

(intercept <- coef(glm1)[1])

## (Intercept) 
##     1.44679

(b_survival <- coef(glm1)[2])

##     Pclass 
## -0.8501067

(logits_survival <- intercept + 2 * b_survival)

## (Intercept) 
##   -0.253424

Thus, the logits of survival are -0.25 Now we can convert to probability:

logit2prob(logits_survival)

## (Intercept) 
##   0.4369809

Lumping logits to probability

Remember that \(e^1 \approx 2.71\). That is, if your logit is 1, your odds will be approx. 2.7 to 1, so the the probability is 2.7 / 3.7, or about 3/4, 75%.

Similarly important, \(e^0 = 1\). Hence, your odds will be 1:1, ie., 50%.

Hence, whenever your logit is negative, the associated probability is below 50% and v.v. (positive logit <–> probability above 50%).

Predict as convenience function

However, more convenient would be to use the predict function instance of glm; this post is aimed at explaining the idea. In practice, rather use:

predict(glm1, data.frame(Pclass = 1), type = "response")

##        1 
## 0.644897

In the 1st class, survival chance is ~65%, and for 2nd class about 44%.

predict(glm1, data.frame(Pclass = 2), type = "response")

##         1 
## 0.4369809

Conversion table

Here’s a look up table for the conversion:

logit_seq <- seq(-10, 10, by = 2)

prob_seq <- round(logit2prob(logit_seq), 3)

df <- data.frame(Logit = logit_seq,
                 Probability = prob_seq)

library(htmlTable)
htmlTable(df)

A handy function is datatable, does not work in this environment however it appears.

library(DT)
datatable(df)

The package mfx provides a convenient functions to get odds out of a logistic regression (Thanks for Henry Cann’s comment for pointing that out!).

Conversion plot

More convenient for an overview is a plot like this.

library(ggplot2)


logit_seq <- seq(-10, 10, by = .1)

prob_seq <- logit2prob(logit_seq)

rm(df)

df <- data.frame(Logit = logit_seq,
                 Probability = prob_seq)

ggplot(df) +
  aes(x = logit_seq, y = prob_seq) +
  geom_point(size = 2, alpha = .3) +
  labs(x = "logit", y = "probability of success")

Takeway

The relationship between logit and probability is not linear, but of s-curve type.

The coefficients in logit form can be be treated as in normal regression in terms of computing the y-value.

Transform the logit of your y-value to probability to get a sense of the probability of the modeled event.

Happy glming!

sessionInfo

sessionInfo()

## R version 3.3.2 (2016-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: macOS Sierra 10.12.4
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] ggplot2_2.2.1 DT_0.2        htmlTable_1.9
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.10        knitr_1.16          magrittr_1.5       
##  [4] munsell_0.4.3       colorspace_1.3-2    R6_2.2.1           
##  [7] rlang_0.1           plyr_1.8.4          stringr_1.2.0      
## [10] highr_0.6           dplyr_0.5.0         tools_3.3.2        
## [13] grid_3.3.2          webshot_0.4.0       gtable_0.2.0       
## [16] checkmate_1.8.2     DBI_0.6-1           htmltools_0.3.6    
## [19] yaml_2.1.14         lazyeval_0.2.0.9000 assertthat_0.2.0   
## [22] rprojroot_1.2       digest_0.6.12       tibble_1.3.0.9002  
## [25] htmlwidgets_0.8     rsconnect_0.8       evaluate_0.10      
## [28] rmarkdown_1.5       labeling_0.3        stringi_1.1.5      
## [31] scales_0.4.1        backports_1.0.5     jsonlite_1.4