blog/binary-logistic-regression-in-r/

[giscus[bot]]

blog/binary-logistic-regression-in-r/

Learn when and how to use a (univariate and multivariate) binary logistic regression in R. Learn also how to interpret, visualize and report results

https://statsandr.com/blog/binary-logistic-regression-in-r/

[aarontimo]

Pure gold - thank you Antoine!

[AntoineSoetewey]

You're welcome Aaron!

[venkatpgi]

Superb Antoine...much needed post for many researchers.
I have come across few situations and trying to find a way out.

1. Variables with unusually weird estimates. When you exponentiate them, the adjusted OR and their C.I are extremely high or low. I read in one post that such variables may be causing a near complete separation or classification of the binary dependent variable and there were some advises as to drop them or perform a Firth's bias reduced logistic regression or some form of penalized regression like LASSO. How do you handle them?
2. Variables to be entered in the model - variable selection. How do you select in your practice and which method do you advise? - filter based (using a 'p' value) or by some random method? I do check for VIF for numerical variables to reduce variables. But could not find a similar approach for categorical variables.

Thanks always.
Venkat

[AntoineSoetewey]

Thanks @venkatpgi!

1. It happens quite often in practice that the width of a CI is huge. I have heard about several methods such as those you mentioned, but very often I simply drop these problematic variables.
2. The method I recommend and use in most cases is the mixed stepwise selection (using AIC). However, I have seen many times and I have already used (in particular in health/medical fields) the following method: first, perform a univariate binary logistic regression separately for each independent variable and second, include significant variables (or close to significance, for instance p-value < 0.2) from the first step in one final model to see whether or not these variables remained significant when using a multivariate binary logistic regression.

Hope this helps.

[tomhiatt]

This is the kind of helpful post that I always wish I had while attempting some past project. Thank you!

Quick mistake: "This means that, the odds of developing a heart disease are (3.574 - 1) × 100 = 257.4% higher for males than for females."

(odds should be replaced with probability)

On that note, using risk instead of odds is easier for non-gamblers to interpret and also easier to communicate findings. Do you have any suggestions on straightforward ways to convert ORs from models to RRs? For example I would love to use tbl_regression() displaying RRs instead of ORs. Is there a way to do that?

[AntoineSoetewey]

Dear @tomhiatt,

Quick mistake: "This means that, the odds of developing a heart disease are (3.574 - 1) × 100 = 257.4% higher for males than for females."
(odds should be replaced with probability)

Good catch, it's corrected, thanks!

On that note, using risk instead of odds is easier for non-gamblers to interpret and also easier to communicate findings. Do you have any suggestions on straightforward ways to convert ORs from models to RRs?

The relationship between OR and RR is rather complex. To convert an OR to a RR, you can use the formula:

$RR = \frac{OR}{(1 – p + (p * OR))},$

where $p$ is the risk in the control group, that is, the prevalence of the outcome in the reference group. This risk in the control group is in our case the chance that the outcome is 1 if the predictor is 0, which can be computed as follows:

$p = \frac{\exp(\hat{\beta}_0)}{1 + \exp(\hat{\beta}_0)}.$

See this post and this thread for more information.

For example I would love to use tbl_regression() displaying RRs instead of ORs. Is there a way to do that?

Based on the package's documentation, I don't think there is an easy way to report RR instead of OR, at least at the moment.

I found this post which stipulates that you can obtain RR with exp(coef(model)) from a glm fit with family = quasipoisson.

As a side note, these two papers may be of interest to you as well:

1. Calculating risk and prevalence ratios and differences in R: developing intuition with a hands-on tutorial and code
2. Common pitfalls in statistical analysis: Odds versus risk

Hope this helps.

[imenedj]

I'm a PhD student currently analyzing my data, and your post is incredibly valuable. I appreciate the effort you put into explaining concepts, which has been particularly helpful for someone like me who is new to data analysis. Could you please provide an explanation for multinomial logistic regression? I need to apply it to my data. Thanks in advance!

[AntoineSoetewey]

Hello @imenedj,

Glad you find it useful!

Writing a tutorial about multinomial logistic regression in R is on my to do list but do not expect it soon because I am quite busy with my PhD thesis at the moment.

In the meantime, here is a great tutorial if you want to learn more: https://bookdown.org/drki_musa/dataanalysis/multinomial-logistic-regression.html

Hope this helps.

Regards,
Antoine

[Eghelakpo]

God bless you for this fantastic piece, Antoine.

[AntoineSoetewey]

Thank you @Eghelakpo !

[nop9898]

This article is incredibly useful (and so is its sibling on linear regression).

Thank you so much for taking the time to put it together.

[AntoineSoetewey]

Thanks @nop9898

I'm glad you found these two posts useful!

[raffaem]

I don't think the normal distribution of the dependent variable is a prerequisite for the linear model

[AntoineSoetewey]

Dear @raffaem

Indeed it's the residuals which should follow approximately a normal distribution. This is explained in more details in the post about linear regression: https://statsandr.com/blog/multiple-linear-regression-made-simple/#conditions-of-application

I have updated the text accordingly. Thanks for your input!

Regards,
Antoine

[hokiestudent]

This post was extremely helpful for a graduate class project, especially for someone relatively new to R and regressions. Thank you for the detailed explanations!

[AntoineSoetewey]

You're welcome @hokiestudent !

[Urtidi]

Dear Antoine, thank you very much for the best article on binary logistic regression in R I have been able to find!

I would like to ask two questions - I apologise in advance if the answer should be obvious:

1. For multivariable binary logistic regression, how can you find out the overall significance of your model?

2. For multivariable binary logistic regression, is it a problem if (some of) the independent variables are known to logically impact each other - or would that not be an issue, as long as the VIF (variance inflation factor) stays below 5 or 10?

Thank you so much for the brilliant article :)

[AntoineSoetewey]

If the main effects show no multicollinearity, and if the interactions are theoretically justified and statistically significant, you can keep the interaction model. Just interpret it with care, as coefficients are conditional on each other (in the presence of interactions, each effect depends on the level of the other variable it’s interacting with).

However, if the interaction terms are not significant, it’s usually better to remove them. This simplifies the model, avoids multicollinearity, and makes interpretation more straightforward, especially if the simpler model performs similarly.

So no need to abandon the full model entirely, but simplify if it adds no value.

Regards,
Antoine

[Urtidi]

Hi Antoine,

Thank you very much for replying on a Sunday morning! What you wrote has really helped me with my work today.

I applied AIC and mixed selection, and ended up with three instead of six interactions. Only two of the interactions are significant at p<0.05; I tried looking for a way to specify a minimum significance level in the AIC/mixed selection step process, but could not find anything.

The simplified model has less multicollinearity than the full one, and one independent variable which had borderline significance in the full model has a p-value >0.05 in the simplified model.

One thing I don't really understand: How can the second order interaction of two independent variables be significant when neither of the two independent variables on their own have p-values <0.05?

Have a good start to the week!

Thank you,

Urtidi

[AntoineSoetewey]

How can an interaction be significant when neither main effect is?

This happens because an interaction term captures a relationship that only emerges when both variables are considered together. Even if each variable alone doesn’t have a significant effect on the outcome, their combination might.

A classic example:

• Suppose neither wind nor temp alone affects the outcome,
• But the effect of wind becomes important only at certain temperatures. That joint effect can show up as a significant interaction.

This is common and not a statistical error. It just means the relationship is more complex than additive effects.

I applied AIC and mixed selection, and ended up with three instead of six interactions. Only two of the interactions are significant at p<0.05; I tried looking for a way to specify a minimum significance level in the AIC/mixed selection step process, but could not find anything.

Also, AIC-based selection doesn’t directly filter based on p-values. It selects models based on overall fit and complexity, not statistical significance. If you decide to simplify further for interpretability, just be transparent about the criteria you used, and always remove terms one at a time (not all at once) because a variable that appears non-significant may become significant once another is removed.

[Urtidi]

Hi Antoine,

Thank you very much for your explanations. They have been a massive help.

Regarding the interactions, e.g. "wind:temp" stands for the effect which wind and temperature have on the dependent variable together, with both wind and temp weighted equally (so "wind:temp" means the exact same thing as "temp:wind")? And comparing this to the effect which they each have individually, you can deduce, for example, that higher temperatures increase (or even create) the effect of wind on the dependent variable?

Have a nice evening!

Urtidi

[AntoineSoetewey]

Regarding the interactions, e.g. "wind:temp" stands for the effect which wind and temperature have on the dependent variable together, with both wind and temp weighted equally (so "wind:temp" means the exact same thing as "temp:wind")?

You're right: an interaction captures the combined effect of the two variables on the outcome, and yes wind:temp is the same as temp:wind, the order doesn’t matter.

And comparing this to the effect which they each have individually, you can deduce, for example, that higher temperatures increase (or even create) the effect of wind on the dependent variable?

Your interpretation is also correct: if the interaction is significant, it suggests that the effect of one variable (e.g. wind) on the dependent variable depends on the level of the other (e.g. temperature). So, for example, higher temperatures might amplify or reduce the effect of wind on the outcome, or even create an effect where there wasn’t one before.

[venkatpgi]

Dear Antoine, I do have a request. If you can do a blog on model calibration and validation after the model output - it will be helpful. Regards Venkat

…

On Thu, May 22, 2025 at 11:59 PM Urtidi ***@***.***> wrote: Dear Antoine, thank you very much for the best article on binary logistic regression in R I have been able to find! I would like to ask two questions - I apologise in advance if the answer should be obvious: 1. For multivariable binary logistic regression, how can you find out the overall significance of your model? 2. For multivariable binary logistic regression, is it a problem if (some of) the independent variables are known to logically impact each other - or would that not be an issue, as long as the VIF (variance inflation factor) stays below 5 or 10? Thank you so much for the brilliant article :) — Reply to this email directly, view it on GitHub <#111 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AQO6R25XI4QFSNHLJSGKLBL27YJRDAVCNFSM6AAAAAB2UQKC5KVHI2DSMVQWIX3LMV43URDJONRXK43TNFXW4Q3PNVWWK3TUHMYTGMRTHAYDEMA> . You are receiving this because you were mentioned.Message ID: ***@***.*** com>

-- Dr.S.Venkataseshan Professor Division of Neonatology Department of Pediatrics PGIMER, Chandigarh India-160012 phone: +91-172-4735794 (Residence) +91-172-2755340 (Office) +91-9478001129 (Personal Mobile) +91-7087008487 (Official Mobile) Mail: ***@***.***

[AntoineSoetewey]

Dear @venkatpgi,

It's now on my to do list, but in the meantime here are some thoughts so you don't have to wait too long.

Model calibration refers to how well the predicted probabilities from your model agree with the observed outcomes. Poor calibration means your model’s predictions are systematically too high or too low. Calibration can be assessed thanks to a calibration plot (also called a reliability diagram) and/or a Hosmer-Lemeshow test for goodness-of-fit.

Model validation checks how well your model performs on new or unseen data. It’s about generalizability. There are usually 2 main types:

1. Internal validation: assesses performance on the training dataset (often via resampling, e.g., cross-validation or bootstrapping)
2. External validation: assesses performance on a completely independent dataset

Common metrics are:

• Discrimination: How well the model distinguishes between events and non-events (usually measured with AUC (Area Under ROC Curve))
• Calibration: (again!) Assesses agreement between predicted and observed probabilities

Here’s a reproducible example in R using simulated data, covering:

1. Fitting a binary logistic regression model
2. Evaluating calibration (calibration plot + Hosmer-Lemeshow test)
3. Evaluating validation (ROC curve + AUC)
4. Performing cross-validation with the {caret} package

# --- Load required packages ---
library(rms)
library(ResourceSelection)
library(pROC)
library(caret)

# --- Simulate a binary outcome dataset ---
set.seed(123)
n <- 500
x1 <- rnorm(n)
x2 <- rbinom(n, 1, 0.5)
x3 <- rnorm(n)

# Generate outcome
lin_pred <- -1 + 0.8 * x1 - 1.2 * x2 + 0.5 * x3
prob <- 1 / (1 + exp(-lin_pred))
y <- rbinom(n, 1, prob)

# Create dataset
sim_data <- data.frame(y = y, x1 = x1, x2 = x2, x3 = x3)

# --- Step 1: Fit the logistic regression model ---
model <- glm(y ~ x1 + x2 + x3, data = sim_data, family = "binomial")

# --- Step 2: Calibration ---

# 2.1 Calibration plot
dd <- datadist(sim_data)
options(datadist = "dd")
val.prob(fitted(model), sim_data$y, m = 10, cex = 0.7)

This will divide data into 10 groups (m = 10) based on predicted probabilities, and plot observed proportions vs predicted probabilities. A well-calibrated model will show points close to the diagonal line.

# 2.2 Hosmer-Lemeshow test
hoslem.test(sim_data$y, fitted(model), g = 10)

g = 10: number of groups used to group predicted probabilities. Null hypothesis: the model fits the data well. A non-significant p-value (e.g., > 0.05) = good model fit.

# --- Step 3: Validation (same dataset) ---

# ROC and AUC
roc_obj <- roc(sim_data$y, fitted(model))
plot(roc_obj, col = "blue", main = "ROC Curve")
auc(roc_obj)

ROC curve shows the trade-off between sensitivity and specificity. AUC closer to 1 indicates better discrimination.

# --- Step 4: Cross-validation with caret ---

# Convert outcome to factor
sim_data$y_factor <- factor(ifelse(sim_data$y == 1, "Yes", "No"))

# Set up 10-fold CV
train_control <- trainControl(
  method = "cv",
  number = 10,
  classProbs = TRUE,
  summaryFunction = twoClassSummary
)

# Fit cross-validated logistic regression
cv_model <- train(
  y_factor ~ x1 + x2 + x3,
  data = sim_data,
  method = "glm",
  family = "binomial",
  metric = "ROC",
  trControl = train_control
)

# Print results
print(cv_model)

This performs 10-fold CV with ROC as the metric, computes AUC in each fold, and reports the mean cross-validated AUC, which is a better estimate of how your model would perform on unseen data. Caret expects the response variable to be a factor with levels like Yes/No for ROC/AUC computation, hence the transformation.

[venkatpgi]

This is great, thank you so much!

[venkatpgi]

Hi Antoine, I exactly followed the steps you have shown. When I try to execute the val.prob function, I get the following error message and a warning message. Error in attributes(.Data) <- c(attributes(.Data), attrib) : 'names' attribute [2] must be the same length as the vector [1] In addition: Warning message: In Ops.factor(p, y) : ‘+’ not meaningful for factors

[AntoineSoetewey]

Do you have a reproducible example?

[venkatpgi]

Interesting discussion. In this relation, I have another query and we all face the problem. Feature selection - selection of variables to be entered in a model. Recently, in one of my analysis, I used the LASSO method for penalised regression for variable selection. But I later read in a article/blog by Frank Harell that LASSO is not apt for variable selection due to high noise. Many papers recommend against using stepwise selection methods.

[AntoineSoetewey]

Dear @venkatpgi,

While lasso is often used for variable selection, especially in high-dimensional data, I agree with experts like Frank Harrell who say we should not rely on it purely for selecting "important" predictors. This is because lasso, like stepwise methods, can be unstable: small changes in the data can lead to different sets of selected variables, especially when predictors are correlated or there’s a lot of noise.

That being said, even if lasso is less reliable for inference or deciding which variables are truly "important", it's still useful for prediction, especially when the goal is parsimony and when you’re dealing with many predictors.

Also, if prediction is your main goal (actually you're doing variable selection for better prediction), then partial least squares (PLS) regression might also be an option.

If you're interested in reading more about ridge, lasso, PCR or PLS regressions, I highly recommend to have a look at sections 6.2 to 6.4 of the well-known book "An Introduction to Statistical Learning" by Gareth James and al. You have probably heard of it (or even read it) but I still mention it because it gave me a good overview of these methods, and it even includes practical examples in R.

Finally, I’m not an expert in this area, but here are a couple of alternatives I’ve heard about:

• Elastic Net: a combination of lasso and ridge regression that is often considered more stable when predictors are correlated.
• Bootstrap stability selection: instead of relying on a single run of lasso, this involves resampling the data multiple times and keeping variables that are consistently selected. It seems to provide a more robust set of predictors.

Of course, no method is perfect, but these seem to be more reliable in many cases.

Regards,
Antoine

[joshuamarie]

Hi again, I am from minutes ago from your other post. From this code you made:

# recode sex
dat$sex <- factor(dat$sex,
  levels = c(FALSE, TRUE),
  labels = c("female", "male")
)

# recode chest_pain
dat$chest_pain <- factor(dat$chest_pain,
  levels = 1:4,
  labels = c("typical angina", "atypical angina", "non-anginal pain", "asymptomatic")
)

Just a suggestion: Since you loaded the entire {dplyr} package, you can take advantage of using mutate() function.
Here's how:

dat |> 
    mutate(
        sex = factor(
            sex, 
            levels = c(FALSE, TRUE),
            labels = c("female", "male")
        ), 
        chest_pain = factor(
            dat,
            levels = 1:4,
            labels = c("typical angina", "atypical angina", "non-anginal pain", "asymptomatic")
        )
    )

This is way more readable and less boilerplate (Shameless promotion: I literally have a live course tutorial about R's metaprogramming, which is the core of {tidyverse} API)

[AntoineSoetewey]

You're right @joshuamarie!

Feel free to propose a merge request with these changes and I'll be happy to merge it.

[AlbertoF53]

A well-thought-out and comprehensive treatment of logistic regression in R. I was impressed by the use of gt() and gtsummary(). Particularly valuable is the illustration of alternative methods for communicating the results of logistic regression in an accessible manner. I would love to see a follow-up in terms of multinomial logistic regression.

[AntoineSoetewey]

Dear @AlbertoF53,

Thank you for the feedback!

A post on multinomial logistic regression is a great suggestion and definitely something I may cover in the future.

Regards,
Antoine

[AlbertoF53]

If you ever consider writing a book on logistic and multinomial regression, I will buy it. Your approach is quite pedagogical.