Binary Logistic Regression

May 23, 2023 · 15 min read · Logistic regression Regression ·

Overview

Binary logistic regression is used to model the relationship between a binary dependent variable (e.g. therapy successful? Yes/No) and a set of explanatory variables. As a result of the regression analysis, we identify explanatory variables that influence the binary dependent outcome variable. The binary logistic regression equation can also be used to predict the outcome of the binary dependent variable as a function of values of the explanatory variables.

Binary logistic regression – Principles and an illustrative data analysis using Gnu R

Objective and scope of application

The binary logistic regression is one of the class of general linear models. It is used to model the relationship between a binary dependent variable (e.g. therapy successful? Yes/No) and a set of explanatory variables. As a result of the regression analysis, we identify explanatory variables that influence the binary dependent outcome variable. The binary logistic regression equation can also be used to predict the outcome of the binary dependent variable as a function of values of the explanatory variables.

Motivating example

Microaneurysms are a early visible manifestation of diabetic retinopathy. They are characterised by small red dots on the retina.

Healthy retina	Retina with microaneurysms

In this blog post we use logistic regression to identify predictors for the occurrence of retinal microaneurysms. The predictors we include in the analysis are age, gender, systolic blood pressure, oral glucose tolerance test (OGTT), total cholesterol and body mass index (BMI) of the patients. We only include data from people over 18 years of age in the analysis. The dataset we use for the example analysis in this blog post was pre-compiled from data of the survey “National Health and Nutrition Examination Survey (NHANES), 2007-2008” (United States Department of Health and Human Services. Centers for Disease Control and Prevention. National Center for Health Statistics, 2012).

Analysis script in Gnu R

Used Gnu R toolboxes

For data import, analysis and visualization of the results we use already existing toolboxes for the Gnu R system. The toolboxes in use are described further down in the post. They are integrated by the following code fragment:

 1library(haven)         # import SPSS files
 2library(dplyr)         # functions for dataframe manipulation 
 3library(MASS)          # stepAIC for model selection
 4library(gtsummary)     # publication-ready analytical and summary tables
 5library(ggplot2)       # high quality plots
 6library(ggstatsplot)   # high quality statistics plots
 7library(pander)        # rendering R objects into Pandoc markdown
 8library(gt)            # output of tibbles as table
 9library(gtExtras)      # extension of gt
10library(car)           # vif for multicollinarity test

Import of the data

We use the function read_sav from the software package haven to import the data file.

1# Import a dataset that is in SAV format
2dataIn <- read_sav("retinalMicroaneurysm.sav")

The following variables are included in the dataset we use in the analysis:

BPXSY1: Systolic blood pressure (first reading) mm Hg
RIAGENDR: Gender of the sample person
RIDAGEYR: Age in years of the sample person at time of screening
LBXGLT: Two Hour Glucose(OGTT) (mg/dL)
OPDUMA: Retinal microaneurysms, worse eye
LBDTCSI: Total cholesterol (mmol/L)
BMXBMI: Body Mass Index (kg/m**2)

We use the function as_factor to convert the variables RIAGENDR and OPDUMA into a factors, preserving the variables and value label attributes.

1# Change some variables to factors
2dataAnalysis <- dataIn %>%
3  mutate(RIAGENDR = as_factor(RIAGENDR)) %>%
4  mutate(OPDUMA = as_factor(OPDUMA))

Short data summary

In the table below we give a summary of the data that we include in the further steps of the analysis.

1# Print summary table (table 1)
2dataAnalysis %>%
3  tbl_summary()

Characteristic	N = 1,046¹
Gender
Male	522 (50%)
Female	524 (50%)
Age in years	58 (48, 68)
Systolic: Blood pres (1st rdg) mm Hg	124 (114, 136)
Two Hour Glucose(OGTT) (mg/dL)	122 (99, 158)
Retinal microaneurysms, worse eye	55 (5.3%)
Total Cholesterol( mmol/L)	5.22 (4.63, 5.95)
Body Mass Index (kg/m**2)	27.8 (24.8, 31.4)
¹ n (%); Median (IQR)

Analysis by means of binary logistic regression

For the analysis by means of logistic regression we use a generalised linear model, which is provided by the package stats as command glm(). The form of the glm() command is glm(formula, family=familytype(link=linkfunction), data=dataTibble). The meaning of the three parameters formula, family and data is:

formula is a symbolic description of the model to be fitted to the data within the analysis. The variable on the left-hand side of a tilde (~) is called the outcome variable (or dependent variable), while the variables on the right-hand side of a tilde are called the predictors (or independent variables) and are joined by plus signs +. Using binary logistic analysis, the outcome variable must be dichotomous (binary).
family is a description of the distribution of the outcome variables and the link function to be used in the model. In our analysis problem we are dealing with a binary (also dichothomous) outcome variable. We use logistic regression (logit).
data is a tibble containing the data of the variables in the model.

The aim of the analysis is to identify a compact model that represents the collected data as well as possible. We use the Akaike information criterion (AIC) as a metric for the quality of model fit (Akaike, 1974; Burnham & Anderson, 2002; Venables & Ripley, 2002). We use a combination of forward and backward selection approach. In the following, we set up a model that includes all existing predictors, we fit this model to the data:

1# define the model and fit it to the data
2glmlogit <- glm(formula = OPDUMA ~ RIAGENDR + RIDAGEYR + BPXSY1 + LBXGLT +
3                  LBDTCSI + BMXBMI,
4                family = binomial(link="logit"),
5                data = dataAnalysis)

In the next step, we examine multicollinarity in the predictor variables. Multicollinearity is present when two or more of the predictors in a regression model are correlated. It can affect our analysis and thus limit the conclusions we can draw. We use the variance inflation factor (VIF) to measure the strength of the correlation between the predictor variables in the regression model (Fox & Weisberg, 2011; Hair Jr. et al., 2019). A heuristic for the evaluation of VIF values:

A value of 1 indicates there is no correlation between a given and all other predictor variables.
Values between 1 and 5 indicate a moderate correlation between a given and all other predictor variables, but this does not yet require attention.
Values greater than 5 indicate a strong correlation between a given predictor variable and all other predictor variables, the estimated coefficients and p-values are most likely unreliable.

To check whether there is multicollinarity in the predictors, we use the vif function of the car toolbox.

1# multicollinarity test using VIF
2glmlogitvif <- vif(glmlogit)

We reformat the result of the multicollinarity test into tibble format to make it easier to use for subsequent outputs in a table and a visualisation.

 1# reformating the result of the multicollinarity test
 2dfglmlogitvif <- tibble(predictors = names(glmlogitvif),
 3                        vif = glmlogitvif) %>%
 4  mutate(predictors = case_when(predictors == "OPDUMA" ~ "Retinal microaneurysms, worse eye",
 5                                predictors == "RIAGENDR" ~ "Gender",
 6                                predictors == "RIDAGEYR" ~ "Age in years",
 7                                predictors == "BPXSY1" ~ "Systolic: Blood pres (1st rdg) mm Hg",
 8                                predictors == "LBXGLT" ~ "Two Hour Glucose(OGTT) (mg/dL)",
 9                                predictors == "LBDTCSI" ~ "Total Cholesterol( mmol/L)",
10                                predictors == "BMXBMI" ~ "Body Mass Index (kg/m**2)",
11                                .default = predictors))

1# render the result of multicollinarity test as table
2dfglmlogitvif %>%
3  gt() %>%
4  cols_label(
5    predictors = "Predictors",
6    vif = "VIF"
7  )

Predictors	VIF
Gender	1.026555
Age in years	1.223686
Systolic: Blood pres (1st rdg) mm Hg	1.155713
Two Hour Glucose(OGTT) (mg/dL)	1.157670
Total Cholesterol( mmol/L)	1.062960
Body Mass Index (kg/m**2)	1.045360

1# visualise the result of multicollinarity test as table
2dfglmlogitvif %>%
3  ggplot(aes(y=predictors, x=vif)) +
4  geom_col(fill ="gray") +
5  theme_bw() +
6  labs(y = "Predictors", x = "VIF") +
7  geom_vline(xintercept = 1, color = "green", linewidth=1.) +
8  geom_vline(xintercept = 5, color = "red", linewidth=1.)

As a result of the test, no multicollinarity could be detected.

Subsequently we give a summary of the results of the binary logistic modelling using this “large” model.

1# summary of the results of the binary logistic modelling in a table
2glmlogit %>%
3  tbl_regression(estimate_fun = purrr::partial(style_ratio, digits = 3),
4                 pvalue_fun = purrr::partial(style_sigfig, digits = 3)) %>%
5  add_significance_stars(hide_ci = FALSE, hide_p = FALSE,
6                         pattern = "{p.value}{stars}")

Characteristic	log(OR)¹	SE¹	95% CI¹	p-value²
Gender
Male	—	—	—
Female	-0.563	0.294	-1.155, 0.002	0.055
Age in years	-0.006	0.013	-0.032, 0.019	0.623
Systolic: Blood pres (1st rdg) mm Hg	0.004	0.008	-0.012, 0.019	0.617
Two Hour Glucose(OGTT) (mg/dL)	0.007	0.002	0.003, 0.012	0.000***
Total Cholesterol( mmol/L)	-0.081	0.139	-0.358, 0.186	0.560
Body Mass Index (kg/m**2)	-0.035	0.029	-0.092, 0.020	0.228
¹ OR = Odds Ratio, SE = Standard Error, CI = Confidence Interval
² p<0.05; p<0.01; **p<0.001

As can be seen in the column headed p-value, not all predictors contribute significantly to the model.

For model selection and refinement we use the stepAIC function of the MASS Toolbox and approach of stepwise regression, which includes both the forward and the backward selection of the predictors (Venables & Ripley, 2002). In the following, the function stepAIC is applied. The trace of the model selection steps with the respective AIC value and the finally selected model is displayed.

 1# model selection the AIC criterion
 2glmlogitaic <- glmlogit %>% stepAIC(trace = T, direction = 'both')
 3## Start:  AIC=427.03
 4## OPDUMA ~ RIAGENDR + RIDAGEYR + BPXSY1 + LBXGLT + LBDTCSI + BMXBMI
 5## 
 6##            Df Deviance    AIC
 7## - RIDAGEYR  1   413.27 425.27
 8## - BPXSY1    1   413.28 425.28
 9## - LBDTCSI   1   413.38 425.38
10## - BMXBMI    1   414.55 426.55
11## <none>          413.03 427.03
12## - RIAGENDR  1   416.84 428.84
13## - LBXGLT    1   424.43 436.43
14## 
15## Step:  AIC=425.27
16## OPDUMA ~ RIAGENDR + BPXSY1 + LBXGLT + LBDTCSI + BMXBMI
17## 
18##            Df Deviance    AIC
19## - BPXSY1    1   413.42 423.42
20## - LBDTCSI   1   413.53 423.53
21## - BMXBMI    1   414.63 424.63
22## <none>          413.27 425.27
23## + RIDAGEYR  1   413.03 427.03
24## - RIAGENDR  1   417.16 427.16
25## - LBXGLT    1   424.52 434.52
26## 
27## Step:  AIC=423.42
28## OPDUMA ~ RIAGENDR + LBXGLT + LBDTCSI + BMXBMI
29## 
30##            Df Deviance    AIC
31## - LBDTCSI   1   413.66 421.66
32## - BMXBMI    1   414.77 422.77
33## <none>          413.42 423.42
34## + BPXSY1    1   413.27 425.27
35## + RIDAGEYR  1   413.28 425.28
36## - RIAGENDR  1   417.39 425.39
37## - LBXGLT    1   425.90 433.90
38## 
39## Step:  AIC=421.66
40## OPDUMA ~ RIAGENDR + LBXGLT + BMXBMI
41## 
42##            Df Deviance    AIC
43## - BMXBMI    1   415.00 421.00
44## <none>          413.66 421.66
45## + LBDTCSI   1   413.42 423.42
46## + BPXSY1    1   413.53 423.53
47## + RIDAGEYR  1   413.58 423.58
48## - RIAGENDR  1   418.07 424.07
49## - LBXGLT    1   426.07 432.07
50## 
51## Step:  AIC=421
52## OPDUMA ~ RIAGENDR + LBXGLT
53## 
54##            Df Deviance    AIC
55## <none>          415.00 421.00
56## + BMXBMI    1   413.66 421.66
57## + LBDTCSI   1   414.77 422.77
58## + BPXSY1    1   414.87 422.87
59## + RIDAGEYR  1   414.99 422.99
60## - RIAGENDR  1   419.49 423.49
61## - LBXGLT    1   426.62 430.62

Below we give a summary of the results of the reduced model.

1# summary of the refined model
2glmlogitaic %>%
3  tbl_regression(estimate_fun = purrr::partial(style_ratio, digits = 3),
4                 pvalue_fun = purrr::partial(style_sigfig, digits = 3)) %>%
5  add_significance_stars(hide_ci = FALSE, hide_p = FALSE,
6                         pattern = "{p.value}{stars}")

Characteristic	log(OR)¹	SE¹	95% CI¹	p-value²
Gender
Male	—	—	—
Female	-0.603	0.290	-1.186, -0.045	0.037*
Two Hour Glucose(OGTT) (mg/dL)	0.007	0.002	0.003, 0.011	0.000***
¹ OR = Odds Ratio, SE = Standard Error, CI = Confidence Interval
² p<0.05; p<0.01; **p<0.001

If we look again at the column with the p-values, the two remaining predictors are significant. In the second column, the coefficients of the model are listed in terms of the logarithm of the odds ratio. An interpretation of the result is easier if we look at the odds ratio itself.

1# summary of the refined model using odds ratio
2glmlogitaic %>%
3  tbl_regression(exponentiate = TRUE,
4                 estimate_fun = purrr::partial(style_ratio, digits = 3),
5                 pvalue_fun = purrr::partial(style_sigfig, digits = 3)) %>%
6  add_significance_stars(hide_ci = FALSE, hide_p = FALSE,
7                         pattern = "{p.value}{stars}")

Characteristic	OR¹	SE¹	95% CI¹	p-value²
Gender
Male	—	—	—
Female	0.547	0.290	0.305, 0.956	0.037*
Two Hour Glucose(OGTT) (mg/dL)	1.007	0.002	1.003, 1.011	0.000***
¹ OR = Odds Ratio, SE = Standard Error, CI = Confidence Interval
² p<0.05; p<0.01; **p<0.001

The odds ratios can be determined from the log(OR) using the exponential function. The odds ratio describes the effect on the output variable if only one of the predictors changes and all others remain constant.

Let us first look at the odds ratio of the OGTT value: The OR is 1.007. This means that with an increase in the OGTT value by one unit (1mg/dl), on average the odds of developing a retinal microaneurysm increase by a factor of 1.007. Or in other words: The odds increase by 0.7% per one-unit OGTT.

For the predictor gender, the odds ratio for women compared to men is 0.547. This means that on average, women have a factor of 0.547 lower odds of developing retinal microaneurysms. Or in other words, the mean risk for women compared to men is reduced by 45.3%.

1# visualization of the odds ratios of the refined model
2ggcoefstats(glmlogitaic, exponentiate = TRUE) +
3  # further modification with the ggplot2 commands
4  # note the order in which the labels are entered
5  ggplot2::scale_y_discrete(labels = c("(Intercept)", "Gender/female",
6                                       "Two Hour Glucose(OGTT) (mg/dL)")) +  ggplot2::labs(y = "Predictors", x = "Odds ratio")

Gnu R toolboxes that we used

For the statistical analysis of the data and the visualisation of results and intermediate results, we used a number of toolboxes. Below is a list of these toolboxes with a short description and links to more detailed information:

tidyverse — This toolbox includes a number of Gnu R software packages that use tidy data data structures and in which the design philosophy and grammar of tidy data are implemented, e.g. pipes. Information on tidy data: H. Wickham, Tidy data, The Journal of Statistical Software, Vol. 59, 2014, Information on the toolbox: https://www.tidyverse.org/.
haven — The toolbox haven enables Gnu R to read and write various data formats used by other statistical software systems. The following formats are supported: SAS, SPSS and Stata. Information on the toolbox: https://haven.tidyverse.org/.
ggplot2 — In the ggplot2 toolbox a system for the declarative creation of graphics is realised, which is based on ‘The Grammar of Graphics’. Information on the toolbox: https://ggplot2.tidyverse.org/.
gtsummary — The package gtsummary enables the creation of publication-ready analytical and summary tables using the R programming language. Information on the toolbox, website: https://www.danieldsjoberg.com/gtsummary/, cheat sheet: https://github.com/rstudio/cheatsheets/blob/main/gtsummary.pdf
ggstatsplot — The toolbox ggstatsplot can be used to annotate visualisations of data analysis results with the results of statistical tests. Information on the toolbox, website: https://indrajeetpatil.github.io/ggstatsplot/.
MASS — The toolbox MASS provides functions and datasets to support Venables and Ripley, “Modern Applied Statistics with S” (4th edition, 2002) (Venables & Ripley, 2002). Information on the toolbox, website: https://cran.r-project.org/web/packages/MASS/index.html.
car — The toolbox car provides functions and datasets associated with the book Fox and Weisberg “An R Companion to Applied Regression” (3rd Edition, 2011) (Fox & Weisberg, 2011). Information on the toolbox, website: https://cran.r-project.org/web/packages/car/index.html.

Information about the hard- and software configuration of the computer on which this post was authored:

1pander(sessionInfo())

R version 4.4.0 (2024-04-24 ucrt)

Platform: x86_64-w64-mingw32/x64

locale: LC_COLLATE=German_Austria.utf8, LC_CTYPE=German_Austria.utf8, LC_MONETARY=German_Austria.utf8, LC_NUMERIC=C and LC_TIME=German_Austria.utf8

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: car(v.3.1-2), carData(v.3.0-5), gtExtras(v.0.5.0), gt(v.0.10.1), pander(v.0.6.5), ggstatsplot(v.0.12.3), ggplot2(v.3.5.1), gtsummary(v.1.7.2), MASS(v.7.3-60.2), dplyr(v.1.1.4) and haven(v.2.5.4)

loaded via a namespace (and not attached): tidyselect(v.1.2.1), farver(v.2.1.2), statsExpressions(v.1.5.4), fastmap(v.1.2.0), TH.data(v.1.1-2), blogdown(v.1.19), bayestestR(v.0.13.2), broom.helpers(v.1.15.0), labelled(v.2.13.0), digest(v.0.6.35), estimability(v.1.5.1), lifecycle(v.1.0.4), survival(v.3.5-8), magrittr(v.2.0.3), compiler(v.4.4.0), rlang(v.1.1.3), sass(v.0.4.9), tools(v.4.4.0), utf8(v.1.2.4), yaml(v.2.3.8), knitr(v.1.46), labeling(v.0.4.3), xml2(v.1.3.6), multcomp(v.1.4-25), abind(v.1.4-5), withr(v.3.0.0), purrr(v.1.0.2), grid(v.4.4.0), datawizard(v.0.10.0), fansi(v.1.0.6), xtable(v.1.8-4), colorspace(v.2.1-0), paletteer(v.1.6.0), emmeans(v.1.10.2), scales(v.1.3.0), zeallot(v.0.1.0), insight(v.0.19.11), cli(v.3.6.2), mvtnorm(v.1.2-5), rmarkdown(v.2.27), generics(v.0.1.3), performance(v.0.11.0), rstudioapi(v.0.16.0), tzdb(v.0.4.0), parameters(v.0.21.7), commonmark(v.1.9.1), cachem(v.1.1.0), stringr(v.1.5.1), splines(v.4.4.0), effectsize(v.0.8.8), vctrs(v.0.6.5), Matrix(v.1.7-0), sandwich(v.3.1-0), jsonlite(v.1.8.8), bookdown(v.0.39), hms(v.1.1.3), patchwork(v.1.2.0), ggrepel(v.0.9.5), correlation(v.0.8.4), fontawesome(v.0.5.2), tidyr(v.1.3.1), jquerylib(v.0.1.4), glue(v.1.7.0), rematch2(v.2.1.2), codetools(v.0.2-20), stringi(v.1.8.4), gtable(v.0.3.5), prismatic(v.1.1.2), munsell(v.0.5.1), tibble(v.3.2.1), pillar(v.1.9.0), htmltools(v.0.5.8.1), R6(v.2.5.1), evaluate(v.0.23), lattice(v.0.22-6), highr(v.0.10), markdown(v.1.12), readr(v.2.1.5), backports(v.1.4.1), broom(v.1.0.6), bslib(v.0.7.0), Rcpp(v.1.0.12), coda(v.0.19-4.1), xfun(v.0.44), zoo(v.1.8-12), forcats(v.1.0.0) and pkgconfig(v.2.0.3)

Theoretical foundations of the logistic regression

Binary logit regression model

Phenomena such as the presence of a disease (yes/no) can be coded by means of binary (dichotomous) variables. Binary logistic regression is used when the dependence of a binary (dichotomous) outcome variable (dependent variable) on predictors (independent variable) is to be modelled. The variable $p$ indicates the probability that an event (e.g. presence of the disease) will occur. The odds of having the disease are defined as the quotient of probability and counter probability $\frac{p}{1 - p}$. By means of logarithm (natural logarithm), the probability values $(0, 1)$ can be mapped into the range of real numbers $(-\infty, \infty)$. The binary logistic regression model is defined by the equation $$ \DeclareMathOperator{\logit}{logit} \logit(p) = \ln\left(\frac{p}{1-p}\right) = \beta_0 + \beta_1 x_1 + \beta_2 x_2+ \cdots +\beta_k x_k , $$ with

$x_i$ the predictors (independent variables)
$\beta_0$ the estimated intercept term
$\beta_i$ the estimated binary logistic regression coefficients.

If we solve this equation for $p$, we get $$ p = \frac{\exp(\beta_0 + \beta_1 x_1 + \beta_2 x_2+ \cdots +\beta_k x_k)}{1 + \exp(\beta_0 + \beta_1 x_1 + \beta_2 x_2+ \cdots +\beta_k x_k)} $$ (Petrie & Sabin, 2005; Rosner, 2016).

Interpretation of the regression coefficients

Continuous predictors

In a binary logistic regression model, consider a single continuous predictor $x_i$ separately that differs between two individuals by a $\beta_{i\Delta}$ and all other predictors in the model have the same values for both individuals. The odds ratio between the two individuals is estimated by $$\hat{OR} = \exp(\hat{\beta}_{i\Delta}).$$

Dichotomous predictors

Suppose there is a dichotomous predictor $x_i$, coded 1 for property present and 0 property absent. If two individuals differ exclusively by this predictor, then the odds ratio can be determined by $$\hat{OR} = \exp(\hat{\beta}_{i}).$$