Beyond OLS: GLMs, logistic regression, and the modelling ecosystem

So far every regression in this course has ended at lm() — ordinary least squares on a continuous outcome. Real applied work reaches further: outcomes that are 0/1 (did the household adopt M-PESA?), counts (how many loan defaults this month?), endogenous regressors that need instruments, panels that need fixed effects, and series that need forecasting. This module is the map of that wider territory, starting with the single most important step past OLS: the generalised linear model.

Why OLS is not enough for a 0/1 outcome

Suppose the outcome is binary: adopted <- 1 if a household uses M-PESA, 0 otherwise. Fit lm(adopted ~ income) and you get a linear probability model. It runs, and the coefficients are sometimes a fine first look — but the fitted line can predict probabilities below 0 or above 1, and the errors are heteroskedastic by construction. For a proper probability model we need to bend the linear predictor onto the 0-1 interval. That is exactly what a generalised linear model does.

The generalised linear model

A GLM keeps the familiar linear predictor beta0 + beta1*x but passes it through a link function that maps it to the natural scale of the outcome. lm() is the special case with a Gaussian family and an identity link. Change the family and you change the model: binomial for 0/1 data, poisson for counts, Gamma for positive skewed data. The interface is one function — glm() — with the same formula syntax you already know.

Logistic regression — the logit link

For a binary outcome the logistic model writes the probability that y = 1 as a logistic (S-shaped) function of the linear predictor. This guarantees a value strictly between 0 and 1:

Rearranged, the model is linear not in the probability but in the log-odds — the logit. The link function is the log of the odds, and that is the scale the coefficients live on:

So beta1 is the change in log-odds for a one-unit change in x. Log-odds are hard to feel intuitively, which is why we almost always exponentiate them into odds ratios — see below.

Fitting a logistic model with glm()

# Toy example: does M-PESA adoption rise with income?
set.seed(1)
hh <- data.frame(
    income_k = round(runif(400, 5, 120)),          # monthly income, KES thousands
    urban    = rbinom(400, 1, 0.5)
)
# true model: higher income and urban households more likely to adopt
p <- plogis(-2 + 0.03 * hh$income_k + 0.8 * hh$urban)
hh$adopted <- rbinom(400, 1, p)

model <- glm(adopted ~ income_k + urban,
             data = hh,
             family = binomial("logit"))
summary(model)

The coefficients in summary(model) are on the log-odds scale. A positive coefficient means higher x raises the probability that y = 1; the size, though, is not a probability and not directly comparable to an lm() slope.

Coefficients as log-odds; exp() to odds ratios

Because logit(p) = beta0 + beta1*x is additive in log-odds, exponentiating turns a coefficient into a multiplicative effect on the odds — an odds ratio. If beta1 = 0.03, then exp(0.03) is about 1.03: each extra KES 1,000 of monthly income multiplies the odds of adoption by about 1.03, a 3% rise in the odds. An odds ratio of 1 means no effect; above 1 raises the odds, below 1 lowers them.

# Odds ratios by hand
exp(coef(model))

# Tidy odds-ratio table with confidence intervals
library(broom)
tidy(model, exponentiate = TRUE, conf.int = TRUE)
#   term        estimate  std.error  statistic  p.value  conf.low  conf.high
#   (Intercept)   ...
#   income_k      ~1.03    ...        ...        ...      ...       ...
#   urban         ~2.2     ...        ...        ...      ...       ...

Fitted probabilities with predict(type = "response")

predict() on a glm defaults to type = "link", returning log-odds. To get probabilities on the 0-1 scale — the thing you usually want — pass type = "response", which applies the inverse link for you.

# Probability of adoption for a new household
newhh <- data.frame(income_k = c(20, 60, 100), urban = c(0, 1, 1))

predict(model, newdata = newhh, type = "link")      # log-odds
predict(model, newdata = newhh, type = "response")  # probabilities, 0-1

# Same by hand: the inverse logit is plogis()
plogis(predict(model, newdata = newhh))

# Fitted probabilities and residuals as a tidy data frame
library(broom)
augment(model, type.predict = "response")

Poisson regression for counts (briefly)

When the outcome is a non-negative count — number of defaults, number of mobile-money transactions in a day — the Poisson family with its log link is the natural GLM. The same glm() call, a different family. Coefficients are on the log scale, so exp() turns them into multiplicative rate ratios, exactly as odds ratios work for logistic.

counts <- data.frame(
    defaults = rpois(200, 3),
    branch_age = runif(200, 1, 30)
)
pois_model <- glm(defaults ~ branch_age,
                  data = counts,
                  family = poisson("log"))
exp(coef(pois_model))   # rate ratios: multiplicative effect on the expected count

The wider ecosystem — where to go next

GLMs cover a lot, but three whole families of applied methods live in dedicated packages. You do not need to master them today; you need to know they exist and what to type when the question arrives.

• Instrumental variables (endogenous regressors): AER::ivreg() or the standalone ivreg::ivreg()
• Panel / fixed effects: fixest::feols() (fast, high-dimensional FE, cluster-robust SEs) or the classic plm::plm()
• Time series and forecasting: the fable package (tidyverse-native) and the older forecast package, with Arima() / ARIMA() for ARIMA models

Instrumental variables — ivreg()

When a regressor is correlated with the error term (simultaneity, omitted variables, measurement error), OLS is biased and you instrument. ivreg uses a formula with a vertical bar separating the structural equation from the instruments: y ~ x | z, read as 'regress y on x, using z as the instrument for the endogenous part'.

library(ivreg)   # or: library(AER); ivreg is exported there too

# Structural eq: consumption ~ mpesa_use (endogenous)
# Instrument: distance to nearest agent (affects use, not consumption directly)
iv <- ivreg(consumption ~ mpesa_use + income |     # structural equation
                          agent_distance + income,  # instruments (income is exogenous, in both)
            data = survey)
summary(iv, diagnostics = TRUE)   # includes weak-instrument and Wu-Hausman tests

Fixed effects — fixest::feols()

For panel data — the same units observed over time — fixed effects absorb every time-invariant difference between units. fixest::feols() puts the fixed-effect factors after a vertical bar and is fast enough for millions of rows and multiple high-dimensional FE dimensions, with cluster-robust standard errors built in.

library(fixest)

# Two-way fixed effects: household and year absorbed
fe <- feols(consumption ~ mpesa_use + income | hh_id + year,
            data = panel,
            cluster = ~hh_id)   # cluster-robust SEs at the household level
summary(fe)

# broom works on feols too
library(broom)
tidy(fe, conf.int = TRUE)

Time series — fable and Arima()

For forecasting, the modern tidyverse-native tool is fable, built on tsibble (a time-aware tibble); the older, widely cited toolkit is Hyndman's forecast package with Arima() and auto.arima(). Both fit ARIMA, ETS, and related models and produce forecasts with prediction intervals.

library(forecast)

# mpesa monthly volume as a time series, 12 periods per year
vol <- ts(mpesa$volume_bn, frequency = 12, start = c(2018, 1))
fit <- auto.arima(vol)          # picks the ARIMA(p,d,q) orders for you
forecast(fit, h = 12)           # 12-month-ahead forecast with intervals

# Or fit a specific ARIMA(1,1,1)
Arima(vol, order = c(1, 1, 1))