Statistical modelling: lm, glm and the formula interface
Fit linear and logistic models with the formula syntax, read the summary critically, predict on new data, and check that the model assumptions hold.
The formula interface
# y depends on x, with an intercept by default
lm(mpg ~ wt, data = mtcars)
# more terms
lm(mpg ~ wt + hp + factor(cyl), data = mtcars)
lm(mpg ~ wt * hp, data = mtcars) # main effects AND interaction
lm(mpg ~ wt + I(wt^2), data = mtcars) # protect arithmetic with I()
lm(mpg ~ ., data = mtcars) # everything else in the frame
lm(mpg ~ wt - 1, data = mtcars) # no intercept| Symbol | Meaning |
|---|---|
~ | Separates response from predictors |
+ | Add a term |
: | Interaction of two terms only |
* | Both main effects plus their interaction |
. | All other columns in the data |
- | Remove a term, including the intercept |
I() | Evaluate arithmetic before it enters the model |
factor(x) | Treat a numeric column as categorical |
Fitting and reading the output
fit <- lm(mpg ~ wt + hp, data = mtcars)
summary(fit) # coefficients, standard errors, t values, R squared
confint(fit) # confidence intervals
anova(fit) # sequential sums of squares
AIC(fit)
# logistic regression for a binary outcome
mtcars$fast <- as.integer(mtcars$qsec < median(mtcars$qsec))
logit <- glm(fast ~ wt + hp, data = mtcars, family = binomial())
exp(coef(logit)) # odds ratios
exp(confint.default(logit))
# prediction, including a transformed response
newdata <- data.frame(wt = c(2.5, 3.5), hp = c(110, 150))
cbind(newdata, predict(fit, newdata, interval = "confidence"))
predict(logit, newdata, type = "response")- A coefficient is the expected change in the response per one-unit change in that predictor, holding the others fixed.
- The reference level of a factor is the baseline, and its coefficient does not appear. Use
relevel()orfct_relevel()to choose it deliberately. - For a
glm,type = "response"gives probabilities; the default gives values on the link scale, which are log odds.
Diagnostics
par(mfrow = c(2, 2))
plot(fit) # residuals, scale-location, leverage, influence
par(mfrow = c(1, 1))
car::vif(fit) # variance inflation: above about 5 is worth a look
shapiro.test(residuals(fit)) # normality of residuals, for small samples
cor(model.matrix(fit)[, -1]) # collinearity between predictors⚠️
A logistic model that does not converge reports coefficients with enormous standard errors rather than an error. Check
fit$converged and the residual deviance; complete separation, where one predictor perfectly splits the outcome, is the usual cause.FAQ
How do I predict when the new data has a factor level the model never saw?
It raises an error about levels. Convert the factor with the same levels as the training data, and decide explicitly what to do with unseen categories: map them to the reference level, or to an explicit
other level added during training.Should I use lm or glm for a proportion?
Use
glm with a binomial family and weights equal to the number of trials, so the variance is modelled correctly. Fitting lm to a proportion ignores that a proportion of 1 from two trials is far less precise than one from 1000 trials.Related
Tidyverse workflow: dplyr and tidyr Debugging R: warnings, factors, NA and performance traps
Last refreshed 2026-09-18.