Logistic Regression in R: Complete Tutorial with Code (glm)

Logistic regression is one of the most important statistical techniques in clinical psychology and health research. Unlike step-by-step multiple regression, which predicts continuous values, logistic regression allows us to predict the probability of a binary event occurring: whether a patient will drop out of treatment, whether they meet diagnostic criteria, or whether they will respond to an intervention. In this complete tutorial, we will learn how to implement it in R step by step, from data preparation to reporting in APA 7 format.

What is logistic regression?

Binary logistic regression is a generalized regression model (GLM) used when the dependent variable is dichotomous (0/1, yes/no, present/absent). While linear regression directly models the expected value of Y, logistic regression models the logarithm of the odds (log-odds or logit) that Y = 1.

The fundamental difference from linear regression is that in logistic regression, the coefficients (B) represent the change in log-odds for each unit increase in the predictor, not a direct change in the dependent variable. To obtain a more intuitive interpretation, we exponentiate the coefficients to get the Odds Ratios (OR), a key measure in odds ratio and clinical research.

In clinical psychology, logistic regression is widely used to predict diagnoses (e.g., presence or absence of a disorder), treatment response (clinically significant improvement yes/no), therapeutic dropout, and risk factors associated with clinical conditions.

When to use it?

Use binary logistic regression when:

• Your dependent variable is dichotomous (two mutually exclusive categories).
• You want to estimate the probability of belonging to a category based on one or more predictors.
• Your predictors can be continuous, categorical, or a mix of both.
• You seek to identify which variables are significant risk or protective factors.

Typical scenarios in psychology: predicting whether a patient will drop out of therapy, whether an individual will develop PTSD after a traumatic experience, whether a pharmacological treatment will produce a clinical response, or whether a student will be identified with learning difficulties.

Practical example: the dataset

We will work with a realistic clinical example. Suppose a research team collected data from N = 280 adult patients with generalized anxiety disorder (GAD) who started a 12-session group cognitive-behavioral therapy program. The goal is to identify which factors predict treatment dropout (defined as attending fewer than 8 of the 12 sessions).

First, we will create the simulated data in R so you can follow the complete tutorial:

# Create reproducible simulated data
set.seed(2026)
n <- 280

age        <- round(rnorm(n, mean = 38, sd = 11))
bai        <- round(rnorm(n, mean = 28, sd = 9))
depression <- rbinom(n, 1, prob = 0.40)
wai        <- round(rnorm(n, mean = 55, sd = 12))

# Generate dependent variable with known relationships
logit_p <- -1.5 + 0.02 * age + 0.07 * bai + 0.80 * depression - 0.05 * wai
prob    <- plogis(logit_p)
dropout <- rbinom(n, 1, prob)

# Create data frame
datos <- data.frame(dropout, age, bai, depression, wai)

# Ensure realistic ranges
datos$age <- pmax(18, pmin(65, datos$age))
datos$bai <- pmax(0, pmin(63, datos$bai))
datos$wai <- pmax(12, pmin(84, datos$wai))
datos$depression <- factor(datos$depression, levels = c(0, 1),
                           labels = c("No", "Si"))

Step 1: Data exploration and preparation

Before fitting any model, we need to understand the structure of our data, check the distribution of the dependent variable, and look for missing data.

# Data frame structure
str(datos)

# Descriptive summary
summary(datos)

# Distribution of the dependent variable
table(datos$dropout)
prop.table(table(datos$dropout))

# Check for missing data
colSums(is.na(datos))

# Descriptives by group (dropout vs. no dropout)
aggregate(cbind(age, bai, wai) ~ dropout, data = datos, FUN = mean)
aggregate(cbind(age, bai, wai) ~ dropout, data = datos, FUN = sd)

It is crucial to check the class balance of the dependent variable. If one category has fewer than 10-15% of cases, we may need special techniques (oversampling, undersampling, or weighting). In our example, we expect approximately 30-35% dropout, which is reasonable.

Step 2: Fit the model

In R, logistic regression is fitted using the glm() function specifying family = binomial. By default, R uses the logit link function.

# Fit the logistic regression model
modelo <- glm(dropout ~ age + bai + depression + wai,
              data = datos,
              family = binomial(link = "logit"))

# View the full summary
summary(modelo)

The formula dropout ~ age + bai + depression + wai indicates that we want to predict dropout from the four predictors in an additive manner. The model summary shows:

• Estimate: B coefficient on the log-odds scale. A positive value indicates that as the predictor increases, the probability of dropout = 1 increases.
• Std. Error: Standard error of the coefficient.
• z value: Wald statistic (Estimate / Std. Error). Tests whether the coefficient is significantly different from zero.
• Pr(>|z|): p-value associated with the Wald test.

At the end of the summary, R displays the null deviance (model without predictors, intercept only) and the residual deviance (model with predictors). The difference between them indicates how much variation the predictors explain. The AIC (Akaike Information Criterion) is also reported, which is useful for comparing models: a lower AIC indicates better relative fit.

Step 3: Interpret the Odds Ratios

Log-odds coefficients are difficult to interpret directly. By exponentiating them, we obtain the Odds Ratios (OR), which are much more clinically intuitive. If you need to quickly verify your calculations, you can use our odds ratio calculator.

# Odds Ratios
OR <- exp(coef(modelo))
round(OR, 3)

# 95% confidence intervals for the OR
IC <- exp(confint(modelo))
round(IC, 3)

# Combined table
tabla_or <- data.frame(
  B       = round(coef(modelo), 3),
  OR      = round(OR, 3),
  IC_inf  = round(IC[, 1], 3),
  IC_sup  = round(IC[, 2], 3),
  p       = round(summary(modelo)$coefficients[, 4], 4)
)
print(tabla_or)

The expected results table would look approximately like this:

How to interpret each OR:

• age (OR = 1.02): Each additional year of age multiplies the odds of dropout by 1.02 (2% more), but the effect is not statistically significant (p = .102). Age is not a relevant predictor in this model.
• bai (OR = 1.07): Each additional point on the BAI multiplies the odds of dropout by 1.07. That is, higher baseline anxiety is significantly associated with a greater probability of treatment dropout (p < .001).
• depressionSi (OR = 2.23): Patients with comorbid depression have 2.23 times the odds of dropping out of treatment compared to those without depression (p = .003). This is a clinically important risk factor.
• wai (OR = 0.95): Each additional point in therapeutic alliance multiplies the odds of dropout by 0.95, meaning it reduces them by 5%. Higher therapeutic alliance is a significant protective factor (p = .001).

Step 4: Evaluate model fit

Unlike linear regression, there is no direct R-squared in logistic regression. We use pseudo-R-squared and other fit measures.

# Nagelkerke pseudo-R squared
install.packages("DescTools")  # only the first time
library(DescTools)
PseudoR2(modelo, which = "Nagelkerke")

# Hosmer-Lemeshow test
install.packages("ResourceSelection")
library(ResourceSelection)
hoslem.test(datos$dropout, fitted(modelo), g = 10)

# Confusion matrix (cutoff = 0.5)
predicciones <- ifelse(fitted(modelo) > 0.5, 1, 0)
confusion    <- table(Observado = datos$dropout, Predicho = predicciones)
print(confusion)

# Classification metrics
accuracy    <- sum(diag(confusion)) / sum(confusion)
sensitivity <- confusion[2, 2] / sum(confusion[2, ])
specificity <- confusion[1, 1] / sum(confusion[1, ])

cat("Accuracy:", round(accuracy, 3), "\n")
cat("Sensitivity:", round(sensitivity, 3), "\n")
cat("Specificity:", round(specificity, 3), "\n")

# Model AIC
AIC(modelo)

Nagelkerke's pseudo-R-squared ranges between 0 and 1 and is interpreted analogously to R-squared in linear regression, although it is not exactly equivalent. Values between .20 and .40 are considered acceptable in social sciences. The Hosmer-Lemeshow test evaluates whether predicted probabilities match observed frequencies: a non-significant result (p > .05) indicates good fit. Sensitivity is the proportion of true dropouts correctly identified, and specificity is the proportion of non-dropouts correctly identified.

Step 5: ROC curve and AUC

The ROC (Receiver Operating Characteristic) curve is the standard tool for evaluating the discriminative ability of a logistic model. It plots sensitivity against 1 - specificity for all possible cutoff points.

# Install and load pROC
install.packages("pROC")
library(pROC)

# Calculate the ROC curve
roc_obj <- roc(datos$dropout, fitted(modelo))

# Plot the ROC curve
plot(roc_obj,
     col = "#2f6b57",
     lwd = 2,
     main = "ROC Curve - Dropout Model",
     print.auc = TRUE,
     auc.polygon = TRUE,
     auc.polygon.col = "#2f6b5722")

# AUC value with confidence interval
auc(roc_obj)
ci.auc(roc_obj)

# Optimal cutoff point (maximizes sensitivity + specificity)
corte_optimo <- coords(roc_obj, "best", ret = c("threshold", "sensitivity", "specificity"))
print(corte_optimo)

The optimal cutoff point is not always 0.5. The coords() function with the "best" criterion uses Youden's index (maximizes sensitivity + specificity - 1) to find the threshold that best separates both groups. In clinical contexts, it may be preferable to prioritize sensitivity (detecting all possible dropouts) over specificity.

Step 6: Check assumptions

Although logistic regression has fewer assumptions than linear regression, there are several we must verify:

1. Linearity of the logit

For continuous variables, the relationship between the predictor and the log-odds must be linear. This is tested using the Box-Tidwell test, which includes the interaction between each continuous predictor and its natural logarithm:

# Box-Tidwell test for linearity of the logit
# Create the interaction terms with the log
datos$age_log <- datos$age * log(datos$age)
datos$bai_log <- datos$bai * log(datos$bai + 1)  # +1 to avoid log(0)
datos$wai_log <- datos$wai * log(datos$wai)

modelo_bt <- glm(dropout ~ age + bai + depression + wai +
                           age_log + bai_log + wai_log,
                 data = datos, family = binomial)
summary(modelo_bt)

# If the *_log terms are NOT significant,
# the linearity assumption is met

2. Absence of multicollinearity

# VIF (Variance Inflation Factor)
library(car)
vif(modelo)

# VIF > 5 indicates problematic multicollinearity
# VIF > 10 indicates severe multicollinearity

3. Independence of observations

Each observation must be independent of the others. This is a design requirement, not something easily verified statistically. In our example, each patient appears only once, so this is met. If we had repeated measures or clustered data (patients within therapists), we would need a multilevel model.

4. Adequate sample size

As we already mentioned, the EPV (events per variable) rule requires at least 10 cases of the less frequent category per predictor. We verify:

# Check EPV
n_eventos   <- min(table(datos$dropout))
n_predictores <- 4
epv <- n_eventos / n_predictores
cat("Events per variable:", epv, "\n")
cat("Criterion met:", epv >= 10, "\n")

Step 7: Stepwise model (optional)

In an exploratory approach, we can use stepwise selection based on AIC to identify the best subset of predictors. However, this technique has important limitations and should not replace theory-driven selection.

# Full model
modelo_full <- glm(dropout ~ age + bai + depression + wai,
                   data = datos, family = binomial)

# Backward selection with AIC
modelo_step <- step(modelo_full, direction = "backward", trace = 1)
summary(modelo_step)

# Compare models with likelihood ratio test
modelo_reducido <- glm(dropout ~ bai + depression + wai,
                       data = datos, family = binomial)
anova(modelo_reducido, modelo_full, test = "Chisq")

The likelihood ratio test compares the reduced model (without age) against the full model. If the result is not significant (p > .05), the simpler model is preferable because it explains practically the same with fewer parameters.

How to report the results in APA 7

Before you submit: reporting B instead of OR, giving an odds ratio with no confidence interval, or fitting more predictors than your events per variable can carry are exactly the things a methodological reviewer picks up on the first read. If your manuscript is already written, run it through the AI Paper Reviewer: a free Reviewer 2 style pre-review that tells you what they will challenge about your logistic model while you can still fix it.

Common mistakes

These are the most common mistakes when performing and interpreting a logistic regression:

Complete reproducible code

Below is the complete script that you can copy and paste directly into R or RStudio:

# =============================================================
# LOGISTIC REGRESSION IN R: COMPLETE TUTORIAL
# Prediction of therapeutic dropout in patients with GAD
# =============================================================

# --- 0. Install and load packages ---
# install.packages(c("DescTools", "ResourceSelection", "pROC", "car"))
library(DescTools)
library(ResourceSelection)
library(pROC)
library(car)

# --- 1. Create simulated data ---
set.seed(2026)
n <- 280

age        <- round(rnorm(n, mean = 38, sd = 11))
bai        <- round(rnorm(n, mean = 28, sd = 9))
depression <- rbinom(n, 1, prob = 0.40)
wai        <- round(rnorm(n, mean = 55, sd = 12))

logit_p <- -1.5 + 0.02 * age + 0.07 * bai + 0.80 * depression - 0.05 * wai
prob    <- plogis(logit_p)
dropout <- rbinom(n, 1, prob)

datos <- data.frame(dropout, age, bai, depression, wai)
datos$age <- pmax(18, pmin(65, datos$age))
datos$bai <- pmax(0, pmin(63, datos$bai))
datos$wai <- pmax(12, pmin(84, datos$wai))
datos$depression <- factor(datos$depression, levels = c(0, 1),
                           labels = c("No", "Si"))

# --- 2. Data exploration ---
str(datos)
summary(datos)
table(datos$dropout)
prop.table(table(datos$dropout))
colSums(is.na(datos))
aggregate(cbind(age, bai, wai) ~ dropout, data = datos, FUN = mean)

# --- 3. Fit the model ---
modelo <- glm(dropout ~ age + bai + depression + wai,
              data = datos, family = binomial(link = "logit"))
summary(modelo)

# --- 4. Odds Ratios and confidence intervals ---
OR <- exp(coef(modelo))
IC <- exp(confint(modelo))
tabla_or <- data.frame(
  B      = round(coef(modelo), 3),
  OR     = round(OR, 3),
  IC_inf = round(IC[, 1], 3),
  IC_sup = round(IC[, 2], 3),
  p      = round(summary(modelo)$coefficients[, 4], 4)
)
print(tabla_or)

# --- 5. Model fit evaluation ---
PseudoR2(modelo, which = "Nagelkerke")
hoslem.test(datos$dropout, fitted(modelo), g = 10)

predicciones <- ifelse(fitted(modelo) > 0.5, 1, 0)
confusion    <- table(Observado = datos$dropout, Predicho = predicciones)
print(confusion)

accuracy    <- sum(diag(confusion)) / sum(confusion)
sensitivity <- confusion[2, 2] / sum(confusion[2, ])
specificity <- confusion[1, 1] / sum(confusion[1, ])
cat("Accuracy:", round(accuracy, 3), "\n")
cat("Sensitivity:", round(sensitivity, 3), "\n")
cat("Specificity:", round(specificity, 3), "\n")

# --- 6. ROC curve and AUC ---
roc_obj <- roc(datos$dropout, fitted(modelo))
plot(roc_obj, col = "#2f6b57", lwd = 2,
     main = "ROC Curve - Dropout Model",
     print.auc = TRUE, auc.polygon = TRUE,
     auc.polygon.col = "#2f6b5722")
auc(roc_obj)
ci.auc(roc_obj)
corte_optimo <- coords(roc_obj, "best",
                       ret = c("threshold", "sensitivity", "specificity"))
print(corte_optimo)

# --- 7. Check assumptions ---
# Multicollinearity
vif(modelo)

# Linearity of the logit (Box-Tidwell)
datos$age_log <- datos$age * log(datos$age)
datos$bai_log <- datos$bai * log(datos$bai + 1)
datos$wai_log <- datos$wai * log(datos$wai)
modelo_bt <- glm(dropout ~ age + bai + depression + wai +
                           age_log + bai_log + wai_log,
                 data = datos, family = binomial)
summary(modelo_bt)

# EPV
n_eventos <- min(table(datos$dropout))
cat("EPV:", n_eventos / 4, "\n")

# --- 8. Stepwise model (optional) ---
modelo_step <- step(modelo, direction = "backward", trace = 1)
summary(modelo_step)

# Compare full vs. reduced model
modelo_reducido <- glm(dropout ~ bai + depression + wai,
                       data = datos, family = binomial)
anova(modelo_reducido, modelo, test = "Chisq")