Generalized Additive Models
Last modified: 2026-10-25
Example 2 (The age at enrollment curve on the logit and risk scales, for each behavior pattern)
Example 2 (The age at enrollment curve on the logit and risk scales, for each behavior pattern) Plot the age at enrollment curve of Example 1, and turn it into a predicted probability of CHD. We draw one curve for each behavior pattern, type A and type B, from the one fitted model. The additive model is linear in its terms only on the logit scale, so there the s(age) curve is a partial effect. It adds to the others. The plotted log odds differ from s(age) only by a constant, the intercept plus the other terms at each curve’s covariate pattern. On the probability scale the effect of age depends on the other predictors, so each curve is drawn at a covariate pattern. That pattern fixes the other predictors. The coefficient of Example 1 for behavior pattern B shifts the whole curve, and both curves hold body mass index and systolic blood pressure fixed at their medians in the WCGS data. Both scales use predict() with type = "link", and the probability scale applies the inverse logit, plogis(), to the link-scale ends of the interval, so that the interval is built on the scale where the fit is approximately normal and stays inside (0, 1).
Show code
library(ggplot2)
library(patchwork)
# Two covariate patterns, one per behavior pattern:
# age varies,
# and the other predictors sit at typical values.
grid <- expand.grid(
age = seq(min(wcgs$age), max(wcgs$age), length.out = 100),
dibpat = levels(wcgs$dibpat)
)
grid$bmi <- median(wcgs$bmi)
grid$sbp <- median(wcgs$sbp)
grid$dibpat <- factor(grid$dibpat, levels = levels(wcgs$dibpat))
# Fit and standard error on the logit (link) scale,
# with the intercept and the other terms held at each pattern's values.
pred <- predict(m_wcgs, grid, type = "link", se.fit = TRUE)
link_scale <- data.frame(
age = grid$age, dibpat = grid$dibpat, fit = pred$fit,
lower = pred$fit - 2 * pred$se.fit, upper = pred$fit + 2 * pred$se.fit
)
# The inverse logit maps the fit and both interval ends to probabilities.
prob_scale <- within(link_scale, {
fit <- plogis(fit)
lower <- plogis(lower)
upper <- plogis(upper)
})
# Both panels share the bands, lines and rug.
# Only the y scale differs.
curve_panel <- function(df) {
ggplot(df, aes(age, fit, colour = dibpat, fill = dibpat)) +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2, colour = NA) +
geom_line(linewidth = 0.9) +
geom_rug(aes(x = age), data = wcgs, inherit.aes = FALSE, alpha = 0.05,
sides = "b") +
scale_colour_manual(values = c(A = "#D55E00", B = "#0072B2")) +
scale_fill_manual(values = c(A = "#D55E00", B = "#0072B2")) +
labs(x = "age at enrollment (years)", colour = "behavior pattern",
fill = "behavior pattern")
}
# The left panel is drawn on the log-odds scale,
# with its ticks put at the log odds of round probabilities,
# so that they read as risks.
risk_ticks <- c(0.03, 0.05, 0.1, 0.2)
p_link <- curve_panel(link_scale) +
scale_y_continuous(breaks = qlogis(risk_ticks),
labels = scales::label_percent()(risk_ticks)) +
labs(y = "probability of CHD (log-odds scale)")
p_prob <- curve_panel(prob_scale) +
scale_y_continuous(labels = scales::label_percent()) +
labs(y = "probability of CHD")
p_link + p_prob + plot_layout(guides = "collect")
# Predicted risks at two ages, for each behavior pattern.
pred_pts <- expand.grid(age = c(40, 55), dibpat = levels(wcgs$dibpat))
pred_pts$bmi <- median(wcgs$bmi)
pred_pts$sbp <- median(wcgs$sbp)
pred_pts$dibpat <- factor(pred_pts$dibpat, levels = levels(wcgs$dibpat))
pred_two <- predict(m_wcgs, pred_pts, type = "link", se.fit = TRUE)
risk_two <- plogis(pred_two$fit)
# The interval is built on the link scale, then mapped through the inverse logit.
risk_lower <- plogis(pred_two$fit - qnorm(0.975) * pred_two$se.fit)
risk_upper <- plogis(pred_two$fit + qnorm(0.975) * pred_two$se.fit)
At the median body mass index and systolic blood pressure, the model gives a CHD probability of 6.6% (95% interval [5.2, 8.4]%) at age 40 for behavior pattern A and 14.7% ([11.9, 18.0]%) at age 55. For behavior pattern B the probabilities are 3.2% ([2.4, 4.3]%) at age 40 and 7.5% ([5.6, 9.8]%) at age 55. The type B curve lies below the type A curve at every age by the same amount on the logit scale, the type B coefficient, yet the gap on the probability scale grows with age, because the inverse logit steepens as the probability grows toward one half. The same logit-scale gap is a larger probability gap at a higher baseline risk.
Show code
library(ggplot2)
model_full <- backfit(bases_ad, penalties_ad, y_ad, lam_ad)
grids_ad <- lapply(1:3, function(j) {
seq(min(X_ad[, j]), max(X_ad[, j]), length.out = 200)
})
# The centered component functions of a fitted model, stacked over the three grids.
component_curves <- function(model) {
unlist(lapply(1:3, function(j) {
drop(natural_basis(grids_ad[[j]], knots_ad[[j]]) %*% model$coefs[[j]]) -
model$centers[j]
}))
}
# A bootstrap band: refit the model on cars drawn with replacement, and take the
# middle 95% of the refitted curves at each grid point.
set.seed(1)
boot_curves <- replicate(200, {
rows <- sample(nrow(X_ad), replace = TRUE)
component_curves(backfit(
lapply(bases_ad, function(B) B[rows, ]), penalties_ad, y_ad[rows], lam_ad
))
})
band_ad <- t(apply(boot_curves, 1, quantile, c(0.025, 0.975)))
# The model works on standardized predictors; the plot shows each one in its
# original units (horsepower, pounds, model year).
ad_center <- attr(X_ad, "scaled:center")
ad_scale <- attr(X_ad, "scaled:scale")
term_labels <- c(horsepower = "horsepower", weight = "weight (lb)",
year = "model year")
to_original <- function(z, j) z * ad_scale[j] + ad_center[j]
effects_ad <- data.frame(
x = unlist(lapply(1:3, function(j) to_original(grids_ad[[j]], j))),
effect = component_curves(model_full),
lower = band_ad[, 1],
upper = band_ad[, 2],
term = factor(rep(predictors_ad, each = 200), levels = predictors_ad,
labels = term_labels[predictors_ad])
)
rug_ad <- data.frame(
x = unlist(lapply(1:3, function(j) to_original(X_ad[, j], j))),
term = factor(rep(predictors_ad, each = nrow(X_ad)), levels = predictors_ad,
labels = term_labels[predictors_ad])
)
ggplot(effects_ad, aes(x, effect)) +
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.2) +
geom_line(linewidth = 0.9) +
geom_rug(aes(x = x), data = rug_ad, inherit.aes = FALSE, alpha = 0.2,
sides = "b") +
facet_wrap(~term, scales = "free_x") +
# Whole-number ticks, so that model years read 70, 75, 80.
scale_x_continuous(
breaks = function(lims) {
ticks <- scales::extended_breaks()(lims)
ticks[ticks == round(ticks)]
},
labels = scales::label_comma()
) +
labs(x = NULL, y = "change in predicted mpg")
Each curve shows how the prediction changes with one predictor while the others are held at their averages, and is read like a coefficient in a linear model, except that its slope varies. The additive form has one limit: it cannot represent an effect of one predictor that depends on another. A two-variable smooth term, or a product feature, restores an interaction, but the model then loses the simple one-curve-per-predictor picture.
In practice, software chooses the smoothing parameters automatically and handles non-Gaussian responses: the mgcv package in R [2], and the pygam and statsmodels packages in Python.