The V Lab
AudiencePublic health, epidemiology, and health sciences learners
Study timeApproximately 120–180 minutes
PrerequisitesProbability, confidence intervals, and basic R

About the tutorial data Every record is simulated with fixed random seeds and contains no real personal health information. The simulation mechanism is solely for teaching; associations in the code must not be interpreted as real-world causal or clinical effects.

How to use this tutorial

This tutorial follows one binary outcome: whether a respiratory infection occurred. Read the sections on probability, odds, and the logit in order, then run the model code. Give particular attention to translating model output back to the probability scale. Every example uses only R’s built-in stats and graphics functionality plus knitr, so the document can be knitted in a clean R session.

Learning objectives

After completing this tutorial, you should be able to:

  1. determine whether logistic regression is appropriate for a research question;
  2. distinguish probability, odds, log-odds, and the odds ratio (OR);
  3. fit and interpret crude and adjusted models with glm(..., family = binomial());
  4. transform model coefficients into ORs, confidence intervals, and predicted probabilities;
  5. handle categorical predictors, nonlinear relationships, and interactions;
  6. distinguish calibration, discrimination, and threshold-based classification and evaluate them on a test set;
  7. investigate residuals, leverage, influential observations, sparse data, and separation;
  8. report findings transparently without overstating causality.

1 Why Logistic Regression?

1.1 An ordinary linear model is unsuitable for a binary outcome

Logistic regression applies when each observation has one of two mutually exclusive outcomes, such as:

  • infection / no infection;
  • death / survival;
  • screen positive / screen negative;
  • service uptake / no uptake.

Putting a 0/1 outcome directly into an ordinary linear regression can produce fitted values below 0 or above 1, and the outcome variance changes with its mean. Logistic regression uses a link function that maps any real-valued linear predictor to a probability between 0 and 1.

Define the event first Before modeling, state which value represents the event. In this tutorial, infection_num = 1 means that infection occurred. Factor levels are explicitly ordered as “No” and “Yes,” but the core models use the numeric 0/1 outcome to leave no ambiguity about the event direction.

1.2 Probability, odds, and the logit

If the event probability is pp, then:

odds=p1−p,logit(p)=log⁡(p1−p). \text{odds}=\frac{p}{1-p}, \qquad \text{logit}(p)=\log\left(\frac{p}{1-p}\right).

Probability is “events divided by all people,” whereas odds are “events divided by non-events.” Their numerical values are generally different.

probability_examples <- c(0.05, 0.20, 0.50, 0.80)
probability_table <- data.frame(
  Probability = probability_examples,
  Odds = probability_examples / (1 - probability_examples),
  Log_odds = qlogis(probability_examples)
)
knitr::kable(
  probability_table,
  digits = 3,
  col.names = c("Probability", "Odds", "Log-odds"),
  caption = "Correspondence among probability, odds, and log-odds"
)
Correspondence among probability, odds, and log-odds
Probability Odds Log-odds
0.05 0.053 -2.94
0.20 0.250 -1.39
0.50 1.000 0.00
0.80 4.000 1.39

When p=0.20p=0.20, the odds are 0.20/0.80=0.250.20/0.80=0.25, or approximately 1:4. It would be incorrect to call 0.25 a 25% event probability. The inverse transformation is:

p=odds1+odds=exp⁡(η)1+exp⁡(η). p=\frac{\text{odds}}{1+\text{odds}} =\frac{\exp(\eta)}{1+\exp(\eta)}.

1.3 Model form and likelihood intuition

A logistic regression with kk predictors is:

log⁡(pi1−pi)=β0+β1Xi1+⋯+βkXik. \log\left(\frac{p_i}{1-p_i}\right) =\beta_0+\beta_1X_{i1}+\cdots+\beta_kX_{ik}.

The model specifies a functional relationship between the predictors and log-odds, not a straight-line relationship between predictors and probability. Coefficients are estimated by maximum likelihood: among candidate parameter values, the algorithm finds those under which the observed sequence of 0s and 1s is most likely.

2 Understand the Data and Frame the Question

2.1 Research question and variable roles

The example asks:

In these simulated data, is vaccination status associated with respiratory infection within one year? How does the association look after accounting for age, BMI, smoking, and residential area? How well does the model predict probabilities for new observations?

Infection is the outcome, vaccination status is the primary explanatory variable, and the remaining variables may be used to reduce confounding or improve prediction. Adjustment decisions should follow the research question, temporal ordering, and subject-matter knowledge—not univariable p-value screening.

Before any outcome exploration, 30% of the complete simulated dataset was set aside by stratified sampling as a test set. All descriptions, exploration, functional-form decisions, and model fitting below use only train_data. The first summary of test_data appears in the final evaluation section.

str(train_data)
## 'data.frame':    839 obs. of  10 variables:
##  $ participant_id: chr  "P0001" "P0002" "P0003" "P0004" ...
##  $ infection_num : int  0 0 0 0 0 0 0 0 1 0 ...
##  $ infection     : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 2 1 ...
##  $ age           : num  26 32 61 33 40 32 59 54 41 41 ...
##  $ age10         : num  -2.4 -1.8 1.1 -1.7 -1 -1.8 0.9 0.4 -0.9 -0.9 ...
##  $ bmi           : num  31.5 31.1 28 30.8 30.3 31.9 23.8 27.3 28.8 23 ...
##  $ bmi5          : num  1.3 1.22 0.6 1.16 1.06 1.38 -0.24 0.46 0.76 -0.4 ...
##  $ smoking       : Factor w/ 2 levels "No","Yes": 2 1 1 1 1 1 2 1 1 1 ...
##  $ vaccinated    : Factor w/ 2 levels "No","Yes": 2 1 2 1 1 1 2 2 1 2 ...
##  $ area          : Factor w/ 2 levels "Urban","Rural": 2 1 1 1 1 2 1 2 1 1 ...
data_summary <- data.frame(
  Sample_size = nrow(train_data),
  Infections = sum(train_data$infection_num),
  Infection_proportion = mean(train_data$infection_num),
  Mean_age = mean(train_data$age),
  Mean_BMI = mean(train_data$bmi)
)
knitr::kable(
  data_summary,
  digits = 3,
  col.names = c(
    "Sample size", "Infections", "Infection proportion",
    "Mean age", "Mean BMI"
  ),
  caption = "Overview of the development/training data"
)
Overview of the development/training data
Sample size Infections Infection proportion Mean age Mean BMI
839 169 0.201 48.6 26.4

2.2 Examine counts and absolute risks first

vaccination_table <- with(
  train_data,
  table(Vaccination = vaccinated, Infection = infection)
)
vaccination_table
##            Infection
## Vaccination  No Yes
##         No  310  98
##         Yes 360  71
infection_by_vaccination <- aggregate(
  infection_num ~ vaccinated,
  data = train_data,
  FUN = function(x) c(n = length(x), events = sum(x), risk = mean(x))
)
infection_by_vaccination <- data.frame(
  Vaccination = infection_by_vaccination$vaccinated,
  Sample_size = infection_by_vaccination$infection_num[, "n"],
  Infections = infection_by_vaccination$infection_num[, "events"],
  Infection_proportion = infection_by_vaccination$infection_num[, "risk"]
)
knitr::kable(
  infection_by_vaccination,
  digits = 3,
  col.names = c(
    "Vaccination", "Sample size", "Infections", "Infection proportion"
  ),
  caption = "Infection counts and proportions by vaccination status"
)
Infection counts and proportions by vaccination status
Vaccination Sample size Infections Infection proportion
No 408 98 0.240
Yes 431 71 0.165 
observed_risk <- tapply(
  train_data$infection_num,
  train_data$vaccinated,
  mean
)
group_n <- table(train_data$vaccinated)
observed_se <- sqrt(observed_risk * (1 - observed_risk) / group_n)

plot(
  seq_along(observed_risk), observed_risk,
  pch = c(16, 17), cex = 1.3,
  col = c(palette_lr["vermillion"], palette_lr["blue"]),
  xaxt = "n", ylim = c(0, max(observed_risk + 2 * observed_se) * 1.10),
  xlab = "Vaccination status", ylab = "Observed infection proportion",
  main = "Begin with absolute risk"
)
axis(1, at = seq_along(observed_risk), labels = names(observed_risk))
arrows(
  seq_along(observed_risk), observed_risk - 1.96 * observed_se,
  seq_along(observed_risk), observed_risk + 1.96 * observed_se,
  angle = 90, code = 3, length = 0.06,
  col = c(palette_lr["vermillion"], palette_lr["blue"])
)
A dot plot compares observed infection proportions in unvaccinated and vaccinated groups, with vertical error bars.

Observed infection proportions by vaccination group; error bars are approximate 95% confidence intervals.

This descriptive comparison is important, but it does not account for compositional differences between groups and cannot automatically be interpreted as an effect caused by vaccination.

3 Fit the First Logistic Regression

3.1 Crude model

crude_model <- glm(
  infection_num ~ vaccinated,
  data = train_data,
  family = binomial(link = "logit")
)
summary(crude_model)
## 
## Call:
## glm(formula = infection_num ~ vaccinated, family = binomial(link = "logit"), 
##     data = train_data)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -1.152      0.116   -9.94   <2e-16 ***
## vaccinatedYes   -0.472      0.174   -2.71   0.0067 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 842.99  on 838  degrees of freedom
## Residual deviance: 835.56  on 837  degrees of freedom
## AIC: 839.6
## 
## Number of Fisher Scoring iterations: 4

By default, glm() reports coefficients on the log-odds scale. Exponentiating the vaccination coefficient gives the crude OR for vaccinated versus unvaccinated participants.

crude_beta <- coef(crude_model)["vaccinatedYes"]
crude_se <- summary(crude_model)$coefficients[
  "vaccinatedYes", "Std. Error"
]
crude_or <- exp(crude_beta)
crude_ci <- exp(crude_beta + qnorm(c(0.025, 0.975)) * crude_se)

crude_result <- data.frame(
  Contrast = "Vaccinated vs unvaccinated",
  OR = crude_or,
  CI_lower = crude_ci[1],
  CI_upper = crude_ci[2]
)
knitr::kable(
  crude_result,
  digits = 3,
  col.names = c("Contrast", "OR", "95% CI lower", "95% CI upper"),
  caption = "Crude odds ratio for vaccination status"
)
Crude odds ratio for vaccination status
Contrast OR 95% CI lower 95% CI upper
vaccinatedYes Vaccinated vs unvaccinated 0.624 0.444 0.877

In the training data, the infection odds ratio for vaccinated versus unvaccinated participants is 0.62 (Wald 95% CI: 0.44 to 0.88). This is a ratio of odds, not a ratio of probabilities or a percentage reduction in risk.

3.2 Adjusted model and meaningful units

Age is expressed per 10 years and centered at 50 years; BMI is expressed per 5 kg/m² and centered at BMI 25. This makes the intercept scientifically plausible and makes continuous-predictor ORs easier to communicate.

adjusted_model <- glm(
  infection_num ~ age10 + bmi5 + smoking + vaccinated + area,
  data = train_data,
  family = binomial()
)

model_coefficient_table <- function(model, labels = NULL) {
  coefficient_matrix <- summary(model)$coefficients
  output <- data.frame(
    Term = rownames(coefficient_matrix),
    Coefficient = coefficient_matrix[, "Estimate"],
    Standard_error = coefficient_matrix[, "Std. Error"],
    OR = exp(coefficient_matrix[, "Estimate"]),
    OR_lower = exp(
      coefficient_matrix[, "Estimate"] -
        qnorm(0.975) * coefficient_matrix[, "Std. Error"]
    ),
    OR_upper = exp(
      coefficient_matrix[, "Estimate"] +
        qnorm(0.975) * coefficient_matrix[, "Std. Error"]
    ),
    P_value = coefficient_matrix[, "Pr(>|z|)"],
    row.names = NULL
  )
  if (!is.null(labels)) {
    replacement <- unname(labels[output$Term])
    output$Term[!is.na(replacement)] <- replacement[!is.na(replacement)]
  }
  output
}

adjusted_labels <- c(
  "(Intercept)" = "Intercept: age 50, BMI 25, non-smoker, unvaccinated, urban",
  age10 = "Age (per 10-year increase)",
  bmi5 = "BMI (per 5 kg/m² increase)",
  smokingYes = "Smoking: Yes vs No",
  vaccinatedYes = "Vaccination: Yes vs No",
  areaRural = "Area: Rural vs Urban"
)

adjusted_results <- model_coefficient_table(adjusted_model, adjusted_labels)
knitr::kable(
  adjusted_results,
  digits = 3,
  col.names = c(
    "Term", "Coefficient", "SE", "OR",
    "OR 95% lower", "OR 95% upper", "p-value"
  ),
  caption = "Multivariable logistic regression for infection (Wald intervals)"
)
Multivariable logistic regression for infection (Wald intervals)
Term Coefficient SE OR OR 95% lower OR 95% upper p-value
Intercept: age 50, BMI 25, non-smoker, unvaccinated, urban -1.600 0.158 0.202 0.148 0.275 0.000
Age (per 10-year increase) 0.424 0.065 1.527 1.344 1.735 0.000
BMI (per 5 kg/m² increase) 0.405 0.105 1.499 1.219 1.843 0.000
Smoking: Yes vs No 0.748 0.208 2.112 1.406 3.173 0.000
Vaccination: Yes vs No -0.652 0.188 0.521 0.361 0.753 0.001
Area: Rural vs Urban 0.468 0.190 1.596 1.100 2.315 0.014

The vaccination interpretation must include the condition “holding the other variables in the model constant.” Interpretations of age and BMI must state their units.

Among observations with the same age, BMI, smoking status, and area, the estimated infection odds ratio for vaccinated versus unvaccinated participants is 0.52. The adjusted OR for each 10-year increase in age is 1.53. These are conditional associations; this simulated observational analysis cannot itself establish causality.

3.2.1 How should the intercept be interpreted?

The intercept is the log-odds when all numeric predictors equal zero and all categorical predictors are at their reference levels. Because age and BMI were centered, it describes a 50-year-old, BMI-25, non-smoking, unvaccinated urban resident. Exponentiating the intercept gives this reference profile’s baseline odds, not an OR comparing two groups.

3.3 Why an OR is not a risk ratio

If the unexposed-group risk is p0p_0 and the OR is θ\theta, the corresponding exposed-group risk is:

p1=θp01−p0+θp0. p_1=\frac{\theta p_0}{1-p_0+\theta p_0}. 

baseline_risks <- c(0.05, 0.20, 0.50)
example_or <- 2
risk_from_or <- example_or * baseline_risks /
  (1 - baseline_risks + example_or * baseline_risks)

or_risk_table <- data.frame(
  Baseline_risk = baseline_risks,
  OR = example_or,
  Corresponding_risk = risk_from_or,
  Risk_ratio = risk_from_or / baseline_risks,
  Risk_difference = risk_from_or - baseline_risks
)
knitr::kable(
  or_risk_table,
  digits = 3,
  col.names = c(
    "Baseline risk", "OR", "Corresponding risk",
    "Risk ratio", "Risk difference"
  ),
  caption = "One OR implies different risk ratios and differences at different baseline risks"
)
One OR implies different risk ratios and differences at different baseline risks
Baseline risk OR Corresponding risk Risk ratio Risk difference
0.05 2 0.095 1.91 0.045
0.20 2 0.333 1.67 0.133
0.50 2 0.667 1.33 0.167

Only when the outcome is rare can an OR and risk ratio be numerically close. Even then, use the correct measure name when reporting results.

4 Return from Coefficients to Predicted Probabilities

4.1 Predict for specific profiles

profile_data <- expand.grid(
  age10 = c(0, 2),
  bmi5 = 0,
  smoking = factor("No", levels = levels(train_data$smoking)),
  vaccinated = factor(
    c("No", "Yes"),
    levels = levels(train_data$vaccinated)
  ),
  area = factor("Urban", levels = levels(train_data$area))
)
profile_data$Age <- profile_data$age10 * 10 + 50
profile_data$BMI <- profile_data$bmi5 * 5 + 25

link_prediction <- predict(
  adjusted_model,
  newdata = profile_data,
  type = "link",
  se.fit = TRUE
)
profile_data$Predicted_probability <- plogis(link_prediction$fit)
profile_data$Probability_lower <- plogis(
  link_prediction$fit - qnorm(0.975) * link_prediction$se.fit
)
profile_data$Probability_upper <- plogis(
  link_prediction$fit + qnorm(0.975) * link_prediction$se.fit
)

knitr::kable(
  profile_data[, c(
    "Age", "BMI", "smoking", "vaccinated", "area",
    "Predicted_probability", "Probability_lower", "Probability_upper"
  )],
  digits = 3,
  col.names = c(
    "Age", "BMI", "Smoking", "Vaccination", "Area",
    "Predicted probability", "Mean probability 95% lower",
    "Mean probability 95% upper"
  ),
  caption = "Model-predicted probabilities for specific covariate profiles"
)
Model-predicted probabilities for specific covariate profiles
Age BMI Smoking Vaccination Area Predicted probability Mean probability 95% lower Mean probability 95% upper
50 25 No No Urban 0.168 0.129 0.216
70 25 No No Urban 0.320 0.242 0.410
50 25 No Yes Urban 0.095 0.070 0.129
70 25 No Yes Urban 0.197 0.144 0.263

Construct the confidence interval on the link scale first, then transform it back with plogis(). These intervals reflect uncertainty in the mean probability estimate for each profile. They are not probability intervals for an individual’s eventual binary outcome.

4.2 Standardized predictions and an average risk difference

ORs are often difficult to interpret. We can copy each person in the training data twice, set everyone to unvaccinated in one copy and vaccinated in the other, and average the two sets of model-predicted probabilities. This procedure is called model standardization or predictive marginalization.

standardized_no <- train_data
standardized_yes <- train_data
standardized_no$vaccinated <- factor(
  "No", levels = levels(train_data$vaccinated)
)
standardized_yes$vaccinated <- factor(
  "Yes", levels = levels(train_data$vaccinated)
)

standardized_risk_no <- mean(predict(
  adjusted_model, newdata = standardized_no, type = "response"
))
standardized_risk_yes <- mean(predict(
  adjusted_model, newdata = standardized_yes, type = "response"
))

standardized_results <- data.frame(
  Scenario = c(
    "Everyone set to unvaccinated",
    "Everyone set to vaccinated",
    "Vaccinated minus unvaccinated"
  ),
  Estimate = c(
    standardized_risk_no,
    standardized_risk_yes,
    standardized_risk_yes - standardized_risk_no
  )
)
knitr::kable(
  standardized_results,
  digits = 3,
  caption = "Standardized probabilities and average risk difference from the adjusted model"
)
Standardized probabilities and average risk difference from the adjusted model
Scenario Estimate
Everyone set to unvaccinated 0.250
Everyone set to vaccinated 0.158
Vaccinated minus unvaccinated -0.093

Standardization does not automatically create a causal effect Setting a predictor to two levels is a calculation. Interpreting the resulting contrast causally requires additional assumptions, including consistency, exchangeability, positivity, correct model form, and reliable measurement.

5 Categorical Predictors, Nonlinearity, and Interaction

5.1 Categorical predictors and reference levels

R creates one indicator variable for a two-level factor. Its coefficient compares the non-reference level with the reference level. Inspect levels explicitly:

lapply(
  train_data[c("smoking", "vaccinated", "area", "infection")],
  levels
)
## $smoking
## [1] "No"  "Yes"
## 
## $vaccinated
## [1] "No"  "Yes"
## 
## $area
## [1] "Urban" "Rural"
## 
## $infection
## [1] "No"  "Yes"
model.matrix(
  ~ smoking + vaccinated + area,
  data = train_data[1:6, ]
)
##    (Intercept) smokingYes vaccinatedYes areaRural
## 1            1          1             1         1
## 2            1          0             0         0
## 3            1          0             1         0
## 4            1          0             0         0
## 5            1          0             0         0
## 10           1          0             0         1
## attr(,"assign")
## [1] 0 1 2 3
## attr(,"contrasts")
## attr(,"contrasts")$smoking
## [1] "contr.treatment"
## 
## attr(,"contrasts")$vaccinated
## [1] "contr.treatment"
## 
## attr(,"contrasts")$area
## [1] "contr.treatment"

Use relevel() to change a reference group. Releveling changes the coefficient representation but not the fitted probability for any observation.

5.2 Linearity of continuous predictors on the logit scale

Logistic regression does not require age itself to follow a normal distribution. The basic model does, however, specify a straight-line relationship between age and the logit. Directly comparing unadjusted age-group infection rates with a conditional prediction for one fixed covariate profile would mix different quantities. Instead, the next code uses a component-plus-residual plot to examine the age term on the same adjusted, conditional logit scale as the fitted model.

age_component <-
  unname(coef(adjusted_model)["age10"]) * train_data$age10
age_partial_residual <-
  age_component + residuals(adjusted_model, type = "working")

age_group <- cut(
  train_data$age,
  breaks = quantile(train_data$age, probs = seq(0, 1, 0.1)),
  include.lowest = TRUE,
  ordered_result = TRUE
)
age_partial_data <- data.frame(
  age = train_data$age,
  partial_residual = age_partial_residual,
  age_group = age_group
)
age_partial_summary <- aggregate(
  cbind(age, partial_residual) ~ age_group,
  data = age_partial_data,
  FUN = mean
)

age_sequence <- seq(
  min(train_data$age),
  max(train_data$age),
  length.out = 200
)
linear_age_component <-
  unname(coef(adjusted_model)["age10"]) * ((age_sequence - 50) / 10)
age_smooth <- lowess(
  train_data$age,
  age_partial_residual,
  f = 2 / 3
)
plot(
  train_data$age, age_partial_residual,
  pch = 16, cex = 0.45,
  col = rgb(107/255, 114/255, 128/255, 0.28),
  xlab = "Age (years)", ylab = "Age component + working residual",
  main = "Check age form on the adjusted conditional logit scale"
)
points(
  age_partial_summary$age,
  age_partial_summary$partial_residual,
  pch = 16, cex = 1.15, col = palette_lr["blue"]
)
lines(
  age_smooth$x, age_smooth$y,
  lwd = 3, lty = 2, col = palette_lr["teal"]
)
lines(
  age_sequence, linear_age_component,
  lwd = 3, col = palette_lr["vermillion"]
)
legend(
  "topleft",
  legend = c(
    "Age-decile means", "LOWESS smooth", "Model-specified linear age term"
  ),
  pch = c(16, NA, NA), lty = c(NA, 2, 1), lwd = c(NA, 3, 3),
  col = c(
    palette_lr["blue"],
    palette_lr["teal"],
    palette_lr["vermillion"]
  ),
  bty = "n"
)
A scatter-and-smooth plot compares age partial residuals with the model-specified linear age term on the adjusted conditional logit scale.

Component-plus-residual plot for age in the adjusted model. Systematic departure of grouped means and the smooth from the model line suggests that the functional form should be reconsidered.

Working residuals can become large when fitted probabilities approach 0 or 1, so this plot is a diagnostic clue rather than a formal verdict. If the smooth persistently departs from the model line, use subject-matter knowledge to compare prespecified quadratic, piecewise-linear, or spline terms, and use resampling to assess whether predictive performance truly improves.

If scientific knowledge or the diagnostic suggests curvature, a squared term, piecewise-linear term, or restricted cubic spline can be prespecified. The example compares a linear age term with a quadratic age term. It demonstrates functional form; it should not become an open-ended search for the smallest p-value.

quadratic_model <- glm(
  infection_num ~ age10 + I(age10^2) + bmi5 + smoking + vaccinated + area,
  data = train_data,
  family = binomial()
)

anova(adjusted_model, quadratic_model, test = "Chisq")
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
833 748 NA NA NA
832 748 1 0.004 0.952 
AIC(adjusted_model, quadratic_model)
df AIC
adjusted_model 6 760
quadratic_model 7 762

5.3 Interaction: an association that depends on context

If the vaccination association may differ by residential area, include a product term between vaccination and area:

interaction_model <- glm(
  infection_num ~ age10 + bmi5 + smoking + vaccinated * area,
  data = train_data,
  family = binomial()
)

interaction_coefficients <- coef(interaction_model)
interaction_vcov <- vcov(interaction_model)

# In urban residents, the vaccination OR uses only vaccinatedYes.
urban_contrast <- c(
  "vaccinatedYes" = 1,
  "vaccinatedYes:areaRural" = 0
)

# In rural residents, the vaccination log(OR) is the sum of the main
# vaccination coefficient and the interaction coefficient.
rural_contrast <- c(
  "vaccinatedYes" = 1,
  "vaccinatedYes:areaRural" = 1
)

contrast_or <- function(model, contrast) {
  all_contrast <- setNames(rep(0, length(coef(model))), names(coef(model)))
  all_contrast[names(contrast)] <- contrast
  estimate <- sum(all_contrast * coef(model))
  standard_error <- sqrt(
    as.numeric(t(all_contrast) %*% vcov(model) %*% all_contrast)
  )
  c(
    OR = exp(estimate),
    lower = exp(estimate - qnorm(0.975) * standard_error),
    upper = exp(estimate + qnorm(0.975) * standard_error)
  )
}

area_specific_or <- rbind(
  Urban = contrast_or(interaction_model, urban_contrast),
  Rural = contrast_or(interaction_model, rural_contrast)
)
knitr::kable(
  data.frame(
    Area = rownames(area_specific_or),
    area_specific_or,
    row.names = NULL
  ),
  digits = 3,
  col.names = c("Area", "Vaccination OR", "95% lower", "95% upper"),
  caption = "Area-specific vaccination odds ratios from the interaction model"
)
Area-specific vaccination odds ratios from the interaction model
Area Vaccination OR 95% lower 95% upper
Urban 0.406 0.254 0.65
Rural 0.780 0.433 1.40

With an interaction in the model, vaccinatedYes describes the vaccination contrast only in the reference area, Urban. The rural contrast must also include the interaction term. Do not interpret a so-called “main effect” in isolation.

6 Model Fit, Comparison, and Diagnostics

6.1 Likelihood, deviance, and pseudo-R2R^2

Lower residual deviance and AIC can help compare candidate models fit to the same outcome and observations. A likelihood-ratio test is available for nested models. None of these replaces external validation or scientific plausibility.

null_model <- glm(
  infection_num ~ 1,
  data = train_data,
  family = binomial()
)

mcfadden_r2 <- 1 - as.numeric(logLik(adjusted_model) / logLik(null_model))
fit_statistics <- data.frame(
  Metric = c(
    "Null-model deviance",
    "Adjusted-model deviance",
    "AIC",
    "McFadden pseudo-R²"
  ),
  Value = c(
    deviance(null_model),
    deviance(adjusted_model),
    AIC(adjusted_model),
    mcfadden_r2
  )
)
knitr::kable(
  fit_statistics,
  digits = 3,
  caption = "Fit statistics in the training data"
)
Fit statistics in the training data
Metric Value
Null-model deviance 842.992
Adjusted-model deviance 748.125
AIC 760.125
McFadden pseudo-R² 0.113 
anova(null_model, adjusted_model, test = "Chisq")
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
838 843 NA NA NA
833 748 5 94.9 0

McFadden’s pseudo-R2R^2 is not the proportion of variance explained from linear regression, and values from different definitions are not directly interchangeable. State the definition when reporting it.

6.2 Logistic regression does not require normal, homoscedastic residuals

The logistic response follows a Bernoulli distribution, whose variance is naturally p(1−p)p(1-p). Do not import linear regression requirements of residual normality and constant variance. Instead, investigate functional form, unusual residuals, leverage, influence, correlated observations, overdispersion, separation, and model specification.

diagnostic_data <- data.frame(
  fitted_probability = fitted(adjusted_model),
  deviance_residual = residuals(adjusted_model, type = "deviance"),
  pearson_residual = residuals(adjusted_model, type = "pearson"),
  leverage = hatvalues(adjusted_model),
  cooks_distance = cooks.distance(adjusted_model)
)

diagnostic_summary <- data.frame(
  Criterion = c(
    "|Deviance residual| > 2",
    "Leverage > 2p/n",
    "Cook's distance > 4/n"
  ),
  Count = c(
    sum(abs(diagnostic_data$deviance_residual) > 2),
    sum(
      diagnostic_data$leverage >
        2 * length(coef(adjusted_model)) / nrow(train_data)
    ),
    sum(diagnostic_data$cooks_distance > 4 / nrow(train_data))
  )
)
knitr::kable(
  diagnostic_summary,
  caption = "Heuristic thresholds for locating observations that require review"
)
Heuristic thresholds for locating observations that require review
Criterion Count
|Deviance residual| > 2 25
Leverage > 2p/n 72
Cook’s distance > 4/n 59 
old_par <- par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))

plot(
  diagnostic_data$fitted_probability,
  diagnostic_data$deviance_residual,
  pch = 16, cex = 0.55,
  col = rgb(0, 114/255, 178/255, 0.40),
  xlab = "Fitted probability", ylab = "Deviance residual",
  main = "Residuals"
)
abline(h = c(-2, 0, 2), lty = c(2, 1, 2), col = palette_lr["grey"])

plot(
  diagnostic_data$leverage,
  type = "h", col = palette_lr["teal"],
  xlab = "Training observation index", ylab = "Leverage",
  main = "Leverage"
)
abline(
  h = 2 * length(coef(adjusted_model)) / nrow(train_data),
  lty = 2, col = palette_lr["vermillion"]
)

plot(
  diagnostic_data$cooks_distance,
  type = "h", col = palette_lr["orange"],
  xlab = "Training observation index", ylab = "Cook's distance",
  main = "Influence"
)
abline(
  h = 4 / nrow(train_data),
  lty = 2, col = palette_lr["vermillion"]
)
Three panels show fitted probability versus deviance residual, observation index versus leverage, and observation index versus Cook's distance.

Deviance-residual, leverage, and Cook’s-distance diagnostics for the logistic regression. Heuristic thresholds are screening aids, not automatic deletion rules. 

par(old_par)

A large residual means that an outcome is discordant with its fitted probability; high leverage indicates an uncommon covariate pattern; and a large Cook’s distance means that deleting an observation might noticeably change the fit. Review data quality, study processes, and sensitivity analyses rather than mechanically deleting observations that cross a cutoff.

6.3 Sparse cells, complete separation, and convergence

If everyone at one predictor level experiences the event—or if no one does—the maximum-likelihood estimate may tend toward positive or negative infinity. This is complete separation. Common signals include:

  • a zero cell in a cross-tabulation;
  • extremely large absolute coefficients and standard errors;
  • a glm() warning about fitted probabilities of 0 or 1, or failure to converge;
  • highly unstable results across software or model variations.
with(train_data, table(infection, smoking, area))
## , , area = Urban
## 
##          smoking
## infection  No Yes
##       No  407  76
##       Yes  69  30
## 
## , , area = Rural
## 
##          smoking
## infection  No Yes
##       No  151  36
##       Yes  45  25
c(
  converged = adjusted_model$converged,
  iterations = adjusted_model$iter,
  events = sum(train_data$infection_num),
  non_events = sum(train_data$infection_num == 0),
  parameters = length(coef(adjusted_model))
)
##  converged iterations     events non_events parameters 
##          1          5        169        670          6

Possible responses include removing unnecessary parameters, combining categories when scientifically defensible, collecting more informative observations, or using penalized or bias-reduced methods. A fixed “10 events per parameter” rule is an oversimplification: event and non-event counts, effect sizes, predictor distributions, missingness, and validation strategy all affect reliability.

6.4 Overdispersion and correlated observations

For independent, individual-level binary data, the Bernoulli variance is determined by the mean. With grouped counts, repeated measures, household or community clustering, a simple logistic regression’s independence and variance specification may be wrong. Depending on the design, consider quasi-binomial methods, generalized estimating equations, mixed-effects models, or design-based uncertainty estimates.

7 Predictive Performance: Calibration and Discrimination

7.1 Evaluate on data that were not used for fitting

The training set was used to estimate coefficients; the test set is used only for final evaluation. The 70/30 split was stratified by event status, so both sets contain events and non-events.

split_summary <- data.frame(
  Dataset = c("Training", "Test"),
  Sample_size = c(nrow(train_data), nrow(test_data)),
  Events = c(sum(train_data$infection_num), sum(test_data$infection_num)),
  Event_proportion = c(
    mean(train_data$infection_num),
    mean(test_data$infection_num)
  )
)
knitr::kable(
  split_summary,
  digits = 3,
  col.names = c("Dataset", "Sample size", "Events", "Event proportion"),
  caption = "Composition of the training and test sets"
)
Composition of the training and test sets
Dataset Sample size Events Event proportion
Training 839 169 0.201
Test 361 73 0.202 
test_probability <- predict(
  adjusted_model,
  newdata = test_data,
  type = "response"
)

7.2 Calibration: do predicted probabilities agree with observed proportions?

The Brier score is the mean squared error between predicted probabilities and 0/1 outcomes; lower is better. It reflects both calibration and discrimination and depends on the outcome’s baseline frequency.

test_brier <- mean((test_data$infection_num - test_probability)^2)
test_log_loss <- -mean(
  test_data$infection_num * log(clamp_probability(test_probability)) +
    (1 - test_data$infection_num) *
      log(1 - clamp_probability(test_probability))
)

calibration_group <- cut(
  test_probability,
  breaks = unique(quantile(
    test_probability,
    probs = seq(0, 1, 0.1)
  )),
  include.lowest = TRUE,
  ordered_result = TRUE
)
calibration_data <- aggregate(
  cbind(predicted = test_probability, observed = test_data$infection_num) ~
    calibration_group,
  FUN = mean
)

calibration_metrics <- data.frame(
  Metric = c("Test-set Brier score", "Test-set log loss"),
  Value = c(test_brier, test_log_loss)
)
knitr::kable(
  calibration_metrics,
  digits = 3,
  caption = "Probability-error metrics in the test set"
)
Probability-error metrics in the test set
Metric Value
Test-set Brier score 0.156
Test-set log loss 0.486 
plot(
  calibration_data$predicted,
  calibration_data$observed,
  pch = 16, cex = 1.15, col = palette_lr["blue"],
  xlim = c(0, 1), ylim = c(0, 1), asp = 1,
  xlab = "Mean predicted probability within group",
  ylab = "Observed event proportion within group",
  main = "Test-set calibration"
)
abline(0, 1, lty = 2, lwd = 2, col = palette_lr["grey"])
A calibration scatterplot compares mean predicted probability with observed event proportion in each group and includes an ideal diagonal line.

Test-set calibration plot grouped by deciles of predicted probability. Ideal calibration lies near the 45-degree line.

A grouped calibration plot depends on the grouping rule and sample size. Do not treat the p-value from a single calibration test as proof that a model “passes.” Examine overall calibration, the probability range relevant to decisions, and important subgroups.

7.3 Discrimination: do cases receive higher scores?

An ROC curve displays the trade-off between sensitivity and 1−specificity over all possible thresholds. The AUC is the probability that the model assigns a higher score to a randomly selected event than to a randomly selected non-event, with ties counting as one-half.

auc_rank <- function(y, probability) {
  n_event <- sum(y == 1)
  n_nonevent <- sum(y == 0)
  rank_sum <- sum(rank(probability, ties.method = "average")[y == 1])
  (rank_sum - n_event * (n_event + 1) / 2) /
    (n_event * n_nonevent)
}

roc_coordinates <- function(y, probability) {
  thresholds <- c(
    Inf,
    sort(unique(probability), decreasing = TRUE),
    -Inf
  )
  true_positive_rate <- vapply(thresholds, function(threshold) {
    predicted <- probability >= threshold
    sum(predicted & y == 1) / sum(y == 1)
  }, numeric(1))
  false_positive_rate <- vapply(thresholds, function(threshold) {
    predicted <- probability >= threshold
    sum(predicted & y == 0) / sum(y == 0)
  }, numeric(1))
  data.frame(
    threshold = thresholds,
    false_positive_rate = false_positive_rate,
    true_positive_rate = true_positive_rate
  )
}

test_auc <- auc_rank(test_data$infection_num, test_probability)
test_roc <- roc_coordinates(test_data$infection_num, test_probability)
plot(
  test_roc$false_positive_rate,
  test_roc$true_positive_rate,
  type = "l", lwd = 3, col = palette_lr["blue"],
  xlim = c(0, 1), ylim = c(0, 1), asp = 1,
  xlab = "False-positive rate (1 - specificity)",
  ylab = "True-positive rate (sensitivity)",
  main = paste0("Test-set ROC: AUC = ", round(test_auc, 3))
)
abline(0, 1, lty = 2, col = palette_lr["grey"])
An ROC curve displays the trade-off between false-positive and true-positive rates, with a diagonal line representing random ranking.

ROC curve for the adjusted model in the held-out test set.

The AUC is 0.648. AUC measures ranking, not probability accuracy. A model with a high AUC can still systematically overestimate risk.

7.4 Thresholds, confusion matrices, and decision consequences

Turning probabilities into 0/1 classifications requires a threshold. A value of 0.5 is not inherently optimal. The threshold should reflect the consequences of missed events and false alarms, resource constraints, acceptable workload, and equity considerations.

classification_metrics <- function(y, probability, threshold) {
  predicted <- as.integer(probability >= threshold)
  tp <- sum(predicted == 1 & y == 1)
  fp <- sum(predicted == 1 & y == 0)
  tn <- sum(predicted == 0 & y == 0)
  fn <- sum(predicted == 0 & y == 1)
  data.frame(
    Threshold = threshold,
    True_positive = tp,
    False_positive = fp,
    True_negative = tn,
    False_negative = fn,
    Accuracy = (tp + tn) / (tp + fp + tn + fn),
    Sensitivity = tp / (tp + fn),
    Specificity = tn / (tn + fp),
    Positive_predictive_value = tp / (tp + fp)
  )
}

threshold_results <- do.call(
  rbind,
  lapply(c(0.20, 0.30, 0.50), function(threshold) {
    classification_metrics(
      test_data$infection_num,
      test_probability,
      threshold
    )
  })
)
knitr::kable(
  threshold_results,
  digits = 3,
  col.names = c(
    "Threshold", "True positive", "False positive",
    "True negative", "False negative", "Accuracy",
    "Sensitivity", "Specificity", "Positive predictive value"
  ),
  caption = "Consequences of alternative classification thresholds in the test set"
)
Consequences of alternative classification thresholds in the test set
Threshold True positive False positive True negative False negative Accuracy Sensitivity Specificity Positive predictive value
0.2 42 107 181 31 0.618 0.575 0.628 0.282
0.3 21 43 245 52 0.737 0.288 0.851 0.328
0.5 6 6 282 67 0.798 0.082 0.979 0.500

Lowering the threshold generally increases sensitivity while reducing specificity. Accuracy can also be dominated by the majority class. With a rare outcome, predicting non-event for everyone may produce deceptively high accuracy.

Model metrics are not permission to deploy Before real use, a model still requires external validation, monitoring for data drift, subgroup performance and calibration checks, an assessment of operational feasibility, and shared decisions about the consequences of false positives and false negatives. Excluding protected characteristics from a model does not guarantee fair results.

8 Missing Data, Confounding, and Reproducibility

8.1 Complete-case analysis changes the represented population

By default, glm() excludes rows with missing values in any model variable. Complete-case analysis is defensible only when the missingness mechanism and the population represented by complete records can be justified.

set.seed(20260812)
logistic_data_missing <- train_data
missing_rows <- sample(seq_len(nrow(logistic_data_missing)), 36)
logistic_data_missing$bmi5[missing_rows] <- NA

missing_summary <- data.frame(
  Total_rows = nrow(logistic_data_missing),
  Missing_BMI = sum(is.na(logistic_data_missing$bmi5)),
  Complete_cases = sum(complete.cases(
    logistic_data_missing[c(
      "infection_num", "age10", "bmi5",
      "smoking", "vaccinated", "area"
    )]
  ))
)
knitr::kable(
  missing_summary,
  col.names = c("Total rows", "Missing BMI", "Complete cases"),
  caption = "Sample composition after introducing missing values"
)
Sample composition after introducing missing values
Total rows Missing BMI Complete cases
839 36 803

Report missing counts for every variable, the number of observations actually used by each model, and differences between included and excluded participants. Depending on the missingness mechanism, consider multiple imputation or sensitivity analysis. Do not fill missing predictors with the outcome mean.

8.2 Variable selection should not rely only on p-values

For causal interpretation, choose confounders using causal structure, temporal order, and domain knowledge; avoid adjusting for mediators or colliders. For prediction, use prespecified candidate predictors, regularization, and resampling validation to control overfitting. Stepwise regression ignores model-selection uncertainty and can yield overly optimistic p-values and performance estimates.

8.3 A minimal reproducibility record

analysis_record <- list(
  outcome_definition = "Respiratory infection within one year: 1=Yes, 0=No",
  model_formula = formula(adjusted_model),
  training_n = nobs(adjusted_model),
  event_count = sum(model.response(model.frame(adjusted_model))),
  factor_levels = lapply(
    train_data[c("smoking", "vaccinated", "area")],
    levels
  ),
  seed_for_split = 20260811,
  r_version = R.version.string
)
analysis_record
## $outcome_definition
## [1] "Respiratory infection within one year: 1=Yes, 0=No"
## 
## $model_formula
## infection_num ~ age10 + bmi5 + smoking + vaccinated + area
## <environment: 0x70731f188>
## 
## $training_n
## [1] 839
## 
## $event_count
## [1] 169
## 
## $factor_levels
## $factor_levels$smoking
## [1] "No"  "Yes"
## 
## $factor_levels$vaccinated
## [1] "No"  "Yes"
## 
## $factor_levels$area
## [1] "Urban" "Rural"
## 
## 
## $seed_for_split
## [1] 20260811
## 
## $r_version
## [1] "R version 4.6.1 (2026-06-24)"

9 Reporting Template

A clear report answers at least these questions:

  1. Population and outcome: Who was included? How was the event defined? What was the follow-up window?
  2. Model: Which link function was used? How were continuous predictors scaled or transformed? What were the reference categories?
  3. Adjustment set: Why were these variables selected? How many rows were excluded because of missingness?
  4. Estimates: Report ORs, 95% confidence intervals, and explicit units; preferably also report meaningful predicted probabilities or absolute contrasts.
  5. Diagnostics: Were nonlinearity, sparse cells, separation, influential observations, and correlated observations considered?
  6. Validation: Was performance estimated in the training set, an internal test set, or external data? Report both calibration and discrimination.
  7. Limitations: Which biases, unmeasured confounding, measurement errors, and generalizability concerns remain?

Four-sentence example

In a simulated cohort, we used multivariable logistic regression with a binomial distribution and logit link to examine infection within one year, prespecifying age (per 10 years), BMI (per 5 kg/m²), smoking, vaccination, and area. Holding the other model variables constant, the infection odds ratio for vaccinated versus unvaccinated participants was 0.52. From the same adjusted model, standardized predicted infection probabilities in the training population were 25% with everyone set to unvaccinated and 15.8% with everyone set to vaccinated, although this contrast is not automatically causal. In the held-out test set, the model had an AUC of 0.65 and a Brier score of 0.156; external validation and evaluation of unmeasured confounding, data drift, and subgroup calibration would still be required.

10 Common Errors at a Glance

Common statement or practice Problem Better approach
“An OR of 0.70 means a 30% risk reduction” Confuses odds with risk Say that the odds are 30% lower, and report predicted risks or a risk ratio
Interpret glm() coefficients directly Default coefficients are log-odds Use exp(coef) for ORs or transform to probabilities
Report a continuous-predictor OR without units The size of the contrast is unknown State per 1, 5, or 10 units
Treat 0.5 as the default best threshold Ignores costs and the baseline rate Choose from decision consequences and report several thresholds
Report only accuracy or AUC Ignores calibration and class imbalance Also report probability error, calibration, and threshold metrics
Evaluate performance in the training data Performance is generally optimistic Use resampling, a test set, or external validation
Select significant variables stepwise Estimates and uncertainty are unstable Prespecify from the question or use validated regularization
Delete every influential observation A threshold is not a deletion rule Verify data, explain its source, and perform sensitivity analysis
Require normally distributed residuals Imports a linear-model assumption Check logit functional form, separation, influence, and independence

11 Exercises and Answers

Exercise 1: From probability to odds

A person’s predicted infection probability is 0.25. What are the corresponding odds and log-odds?

Answer: The odds are 0.25/(1−0.25)=1/3≈0.3330.25/(1-0.25)=1/3\approx0.333; the log-odds are log⁡(1/3)≈−1.10\log(1/3)\approx-1.10.
Exercise 2: Interpret a continuous-predictor OR

The OR for age10 is 1.40. How should it be interpreted?

Answer: Among observations with the other model variables held constant, a 10-year increase in age multiplies the event odds by 1.40, or raises the odds by 40%. This does not mean the event probability increases by 40%, nor does it automatically establish a causal effect of age.
Exercise 3: Probability output

Why can predict(adjusted_model, newdata = x) not be treated directly as a probability?

Answer: The default prediction from a binomial glm is on the link scale, here log-odds. Use type = "response" or apply plogis() to a link-scale prediction.
Exercise 4: Choose a threshold

If missing a genuinely high-risk person is far more costly than a false alarm, in which direction should the threshold generally move?

Answer: Usually lower the threshold to increase sensitivity. This creates more false positives and lowers specificity, so the final choice must also consider resources, harms of follow-up interventions, and equity.
Exercise 5: Interpret an interaction

With vaccinated * area in the model, can vaccinatedYes be interpreted as the average vaccination OR across all areas?

Answer: No. It is the vaccination OR in the reference area, Urban. The rural contrast also includes vaccinatedYes:areaRural. For a population-average contrast, calculate standardized predictions in a clearly defined target population.
Exercise 6: Diagnose separation

All 18 members of a rare exposure group are non-events, and the model returns an extremely large negative coefficient and standard error. What is the likely problem?

Answer: Complete or quasi-complete separation is likely. Check cross-tabulations and data quality first, then consider reducing parameters, obtaining more information, or using penalized or bias-reduced methods. Do not treat an enormous OR as stable evidence.

12 Quick Reference

Goal Base R expression
Fit a model glm(y ~ x1 + x2, family = binomial(), data = d)
Inspect log-odds coefficients coef(model)
Calculate ORs exp(coef(model))
Wald intervals for ORs exp(coef(model) + outer(sqrt(diag(vcov(model))), qnorm(c(.025, .975)))) (then arrange dimensions carefully)
Predict probabilities predict(model, newdata = d_new, type = "response")
Link-scale prediction and SE predict(model, newdata = d_new, type = "link", se.fit = TRUE)
Deviance residuals residuals(model, type = "deviance")
Leverage hatvalues(model)
Cook’s distance cooks.distance(model)
Compare nested models anova(model_small, model_large, test = "Chisq")
Information criterion AIC(model)

Final checklist

This tutorial provides a framework for learning logistic regression, not a clinical or public health decision rule. Real analyses should incorporate the study design, sampling structure, subject-matter knowledge, preregistered plans, and appropriate software validation.