Logistic regression

A binary outcome takes two values: infected or not, died or survived, bled or did not. You cannot fit an ordinary straight line to a 0/1 variable, because a line runs past 0 and past 1 and predicts risks that cannot exist. Logistic regression fixes this by modelling the outcome on the log-odds scale, where the numbers are free to run from minus infinity to plus infinity. This lesson builds the model from a single two-group comparison up to a model that adjusts for a confounder, and reads every coefficient along the way. It assumes you can already read a 2x2 table, an odds ratio, and a P-value.

From probability to log odds: the logit

Write the risk of the outcome as , the probability a person has the disease. The odds are . The logit is the natural log of those odds. That is the quantity logistic regression models.

The transform buys you range. A probability is trapped between 0 and 1. The odds open that up to between 0 and infinity. Taking the log removes the last wall, so the log odds run from minus infinity to plus infinity. Now a linear model can sit on the right-hand side and never predict an impossible risk.

The right-hand side is the linear predictor. Each is an exposure or a confounder, and each is its coefficient on the log-odds scale. You read it like multiple linear regression, with one change: the quantity being predicted is the log odds, not the outcome itself.

Two groups: a coefficient is a log odds ratio

Take a cohort of 500 adults in Selangor followed for active tuberculosis, split by current smoking. Code the outcome TB as 1 and no TB as 0, and code smoking as 1 for smokers (the exposed group) and 0 for non-smokers (the baseline group). Always code the outcome event as 1, or the model reports the odds of the wrong thing.

The odds ratio straight from the table is 0.3333 / 0.1111 = 3.0. Smokers have three times the odds of TB. Now read the same numbers as a model. The model says the log odds in each group is built from two parameters: a baseline and an exposure effect.

For non-smokers, , so the log odds is . That is , the log odds in the baseline group. For smokers, , so the log odds is . The difference is , and a difference of two log odds is a log odds ratio. So . On the odds scale the model multiplies.

So returns the odds ratio, and returns the baseline odds. The model has not invented anything. It has repackaged the 2x2 table so that the same machinery extends to several variables at once.

The standard error of the log odds ratio comes from the four cell counts, exactly as it did before you met regression.

Here that is . The Wald test divides the coefficient by its standard error: , which gives P below 0.001. A 95 percent interval for the log odds ratio is , that is 0.604 to 1.593. Exponentiate the ends to get the odds ratio interval, 1.83 to 4.92. Notice it is built on the log scale and then antilogged, never the other way round.

# grouped counts: one row per exposure group
tb    <- c(50, 30)
notb  <- c(150, 270)
smoke <- c(1, 0)            # 1 = smoker, 0 = non-smoker

m <- glm(cbind(tb, notb) ~ smoke, family = binomial)
# inspect coef(m) and exp(coef(m))

m <- glm(cbind(tb, notb) ~ smoke, family = binomial)
coef(m)            # (Intercept) -2.197 , smoke 1.099
exp(coef(m))       # baseline odds 0.111 , OR 3.0
exp(confint.default(m)["smoke", ])   # 1.83 to 4.92

import numpy as np
import statsmodels.api as sm
tb    = np.array([50, 30])
notb  = np.array([150, 270])
smoke = np.array([1, 0])

X = sm.add_constant(smoke)             # adds the intercept column
y = np.column_stack([tb, notb])        # successes, failures
# fit a Binomial GLM and exp the params

X = sm.add_constant(smoke)
y = np.column_stack([tb, notb])
m = sm.GLM(y, X, family=sm.families.Binomial()).fit()
print(m.params)            # const -2.197 , x1 1.099
print(np.exp(m.params))    # baseline odds 0.111 , OR 3.0

Predicted probabilities: from log odds back to risk

Coefficients are easy to compute but hard to feel. Convert the linear predictor back to a probability and the model speaks in plain risk. Invert the logit and you get the logistic function.

For a non-smoker the linear predictor is , so . For a smoker it is , so . Check those against the table: 30/300 = 0.10 and 50/200 = 0.25. The predicted risks land on the observed group rates, which is the sanity check you should run every time.

b0 <- -2.197   # intercept (log odds, non-smokers)
b1 <-  1.099   # smoke coefficient

# plogis() is the logistic function 1/(1+exp(-z))
plogis(b0)        # risk for non-smokers
plogis(b0 + b1)   # risk for smokers

plogis(b0)        # 0.10  matches 30/300
plogis(b0 + b1)   # 0.25  matches 50/200

import numpy as np
b0, b1 = -2.197, 1.099

def sigmoid(z):
    return 1 / (1 + np.exp(-z))
# evaluate at x = 0 and x = 1

def sigmoid(z):
    return 1 / (1 + np.exp(-z))
print(sigmoid(b0))        # 0.10
print(sigmoid(b0 + b1))   # 0.25

For a continuous exposure the same formula traces an S-shaped curve. Plot it once and the shape stays with you: risk rises slowly, then fast through the middle, then flattens as it approaches 1.

b0 <- -4; b1 <- 1.2          # log odds at x = 0, slope per unit
x  <- seq(0, 6, by = 0.1)    # cigarettes per day

p <- plogis(b0 + b1 * x)
# plot p against x

p <- plogis(b0 + b1 * x)
plot(x, p, type = "l", lwd = 2,
     xlab = "cigarettes per day", ylab = "P(TB)",
     main = "Logistic dose-response")

import numpy as np
import matplotlib.pyplot as plt
b0, b1 = -4, 1.2
x = np.linspace(0, 6, 100)

p = 1 / (1 + np.exp(-(b0 + b1 * x)))
# plot p vs x

p = 1 / (1 + np.exp(-(b0 + b1 * x)))
plt.plot(x, p, lw=2)
plt.xlabel("cigarettes per day"); plt.ylabel("P(TB)")
plt.title("Logistic dose-response"); plt.show()

More than two groups: indicator variables

When an exposure has more than two levels, you compare each non-baseline level with one reference group. The model does this with indicator variables: variables that take only 0 or 1, one for each non-baseline level. Most packages create them automatically when you declare the variable as a factor. Each indicator carries the log odds ratio for its level against the baseline.

Stay with the same cohort and add household crowding at three levels: low, medium, high. Low crowding is the baseline because it is the natural reference.

Odds: low 20/400 = 0.05, medium 35/350 = 0.10, high 60/300 = 0.20. The odds ratios are 0.10/0.05 = 2.0 for medium and 0.20/0.05 = 4.0 for high. The model produces two coefficients, and , plus the baseline. The original three-level variable is never entered on its own, only its indicators.

The likelihood ratio test for a variable

The honest test for a categorical variable as a block is the likelihood ratio test. Fit the model with crowding, fit the model without it, and compare how well each explains the data. The measure of fit is the log-likelihood, written . The test statistic is twice the gap, which for grouped data equals the deviance form below.

Here is each observed cell count and is the count expected under the no-crowding model, which spreads the overall TB rate 115/1165 = 0.0987 across the groups. Summing over the six cells gives . You compare it with a chi-squared distribution on 2 degrees of freedom, one per indicator, which gives P below 0.001. So crowding belongs in the model, tested as a single variable.

tb    <- c(20, 35, 60)
notb  <- c(400, 350, 300)
crowd <- factor(c("low","med","high"),
                levels = c("low","med","high"))

m1 <- glm(cbind(tb, notb) ~ crowd, family = binomial)
m0 <- glm(cbind(tb, notb) ~ 1,     family = binomial)
# compare m0 and m1

m1 <- glm(cbind(tb, notb) ~ crowd, family = binomial)
exp(coef(m1))             # baseline 0.05 , OR med 2.0 , OR high 4.0
m0 <- glm(cbind(tb, notb) ~ 1, family = binomial)
anova(m0, m1, test = "LRT")   # deviance drop 31.0 on 2 df, P < 0.001

import numpy as np, pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
df = pd.DataFrame({"tb":[20,35,60], "notb":[400,350,300],
    "crowd": pd.Categorical(["low","med","high"],
                            categories=["low","med","high"])})
df["n"] = df.tb + df.notb
df["prop"] = df.tb / df.n
m0 = smf.glm("prop ~ 1", data=df, freq_weights=df.n, family=sm.families.Binomial()).fit()
m1 = smf.glm("prop ~ crowd", data=df, freq_weights=df.n, family=sm.families.Binomial()).fit()

m1 = smf.glm("prop ~ crowd", data=df, freq_weights=df.n,
             family=sm.families.Binomial()).fit()
m0 = smf.glm("prop ~ 1", data=df, freq_weights=df.n,
             family=sm.families.Binomial()).fit()
# LR statistic = 2*(llf_big - llf_small)

print(np.exp(m1.params))          # Intercept 0.05 , med 2.0 , high 4.0
lr = 2 * (m1.llf - m0.llf)
df_diff = m1.df_model - m0.df_model
print(round(lr, 1), int(df_diff)) # 31.0  2

Adjusting for a confounder with multiple logistic regression

A confounder is a variable tied to both the exposure and the outcome that distorts the crude association. Suppose smokers in this sample cluster in high-crowding housing, and crowding itself raises TB. Then the crude smoking odds ratio mixes the smoking effect with the crowding effect. Multiple logistic regression separates them by putting both on the linear predictor and estimating each while holding the other fixed.

Within each crowding stratum the smoking odds ratio is exactly 2.0: (20×80)/(80×10) = 2.0 and (40×90)/(60×30) = 2.0. Collapse the strata and the crude odds ratio drops to (60×170)/(140×40) = 1.82, because crowding pulls the two together. The model term for smoking, fitted alongside crowding, recovers the within-stratum value near 2.0. That is the adjusted odds ratio.

tb    <- c(20, 10, 40, 30)
notb  <- c(80, 80, 60, 90)
smoke <- c(1, 0, 1, 0)
crowd <- c(0, 0, 1, 1)        # 0 = low, 1 = high

m <- glm(cbind(tb, notb) ~ smoke + crowd, family = binomial)
# adjusted OR for smoking = exp(smoke coefficient)

m <- glm(cbind(tb, notb) ~ smoke + crowd, family = binomial)
exp(coef(m))["smoke"]        # adjusted OR near 2.0 (crude was 1.82)

# Mantel-Haenszel summary OR over the two strata
arr <- array(c(20,10,80,80, 40,30,60,90), dim = c(2,2,2))
mantelhaen.test(arr, correct = FALSE)   # common OR 2.0

import numpy as np, pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
d = pd.DataFrame({"tb":[20,10,40,30], "notb":[80,80,60,90],
                  "smoke":[1,0,1,0], "crowd":[0,0,1,1]})
d["n"] = d.tb + d.notb
d["prop"] = d.tb / d.n

m = smf.glm("prop ~ smoke + crowd", data=d, freq_weights=d.n,
            family=sm.families.Binomial()).fit()
# exponentiate to read the adjusted OR for smoking

m = smf.glm("prop ~ smoke + crowd", data=d, freq_weights=d.n,
            family=sm.families.Binomial()).fit()
print(np.exp(m.params))      # smoke OR near 2.0 , crowd OR its own effect

The adjusted smoking odds ratio of about 2.0 sits close to the Mantel-Haenszel summary, which is the link back to Chapter 18: logistic regression and Mantel-Haenszel answer the same confounding question, and the model just does it for any number of variables at once. The model assumes the smoking odds ratio is the same in every crowding stratum, that is, no interaction. If that assumption is wrong you fit an interaction term instead.

Common mistakes

• Reading exp(beta) as a risk ratio. It is an odds ratio. When the outcome is common, the odds ratio is further from 1 than the risk ratio, so this overstates the effect. Only when the outcome is rare are the two close.
• Testing a multi-level variable by one indicator. A single indicator P-value compares one level with the baseline, not the variable as a whole. Use the likelihood ratio test on the full set of indicators.
• Getting the outcome coding backwards. If you code the event as 0 and the non-event as 1, every odds ratio flips to its reciprocal. Confirm that the event is coded 1 before you trust a sign.
• Reporting coefficients on the log scale. A coefficient of 1.099 means nothing to a clinician. Exponentiate to 3.0 and report the odds ratio with its 95 percent interval.
• Adjusting for a variable on the causal pathway. If crowding caused smoking to cause TB along a single chain, adjusting for a mediator removes part of the real effect. Adjust for confounders, not mediators.

Tips

• Exponentiate coefficients and confidence limits, then check that lower × upper around the point estimate is symmetric on the ratio scale. It is a quick audit of a published table.
• Choose the baseline level deliberately. A large, clinically natural reference group gives tighter, more interpretable odds ratios.
• For an ordered exposure such as age group, also fit it as one linear term and compare with the indicator model. A single coefficient per step is simpler when the log odds rise roughly in a straight line.
• Use the likelihood ratio test to decide whether a variable stays in the model, and the Wald z for a single coefficient. They usually agree, and the LRT is the more reliable of the two for small samples.
• Keep one dataset and one set of variable names across the whole analysis, so a reviewer never has to relearn your coding.