GAM Exercises in R: 16 mgcv Practice Problems with Solutions

Explanation: s(Temp) requests a thin-plate regression spline whose smoothness is chosen automatically by minimizing GCV (generalized cross-validation). The estimated 3.69 effective degrees of freedom (edf) is well above 1, which means the Ozone-Temp relationship is meaningfully nonlinear (a straight line would give edf near 1). Switching to lm(Ozone ~ Temp) would force linearity and miss the strong curvature visible at high temperatures.

Explanation: GAMs are additive: each s() term is a separate one-dimensional function of its predictor, and their effects sum. The GCV score dropped from 575.7 to 410.2 versus the single-smooth model, which is strong evidence that Wind adds real explanatory power. The wind smooth uses 2.69 edf, signalling mild but real curvature (Ozone falls sharply at low wind speeds and levels off).

Explanation: Anything outside s() is treated parametrically, just like in lm(). You read the Solar.R coefficient (0.0546 ppb per W/m^2) the same way as a linear regression slope. The smooth-term table is separate because smooths do not have a single slope. Always start parametric if theory supports linearity: each smooth eats 3-5 edf, so over-smoothing inflates variance.

Explanation: k is the maximum basis dimension (the upper bound on edf, roughly k - 1). Setting it small forces a stiffer fit. The penalty still shrinks edf below k - 1 whenever the data permit, so k is a ceiling not a target. A common workflow is to start with k = 10, run gam.check(), and increase k only if its k-index near 1 with a small p-value flags under-smoothing.

Explanation: REML treats smoothing parameters as variance components and integrates out the fixed effects, which makes the marginal likelihood less prone to multiple local minima than GCV. The REML score on this fit is 519.2 (smaller is better but REML scores from different models are only comparable when the fixed-effect structure is identical). For most published GAM work since 2011, method = "REML" is the default recommendation.

Explanation: bs = "cr" places knots at quantiles of the covariate and uses cubic regression splines with a wiggliness penalty on the second derivative. The fit is nearly identical to thin-plate (REML score 519.5 vs 519.2) but scales much better when n runs into the hundreds of thousands. Other useful bases: "ps" for P-splines, "tp" (default thin-plate), "cc" for cyclic smooths, "re" for random effects.

Explanation: The factor is treated parametrically with treatment contrasts: 4-cylinder is the reference, 6-cylinder cars average 3.64 mpg lower at the same weight, 8-cylinder average 6.08 mpg lower. The smooth captures the nonlinear weight effect that is shared across cylinder groups (parallel curves). If you suspect the weight curve has a different shape per cylinder bucket, you need a factor-by smooth (next exercise), not just an additive factor.

Explanation: Each s(wt):cyl_fX is a smooth that operates only on rows where cyl_f == X. With only 32 rows in mtcars, the penalty shrinks every per-group smooth to an effective straight line (edf = 1). For factor-by smooths in real datasets, you also need the parametric main effect cyl_f because by = smooths are mean-zero within each level and cannot capture vertical offsets between groups.

Explanation: bs = "re" makes mgcv treat the factor as a Gaussian random effect, equivalent to a (1 | Chick) term in lme4. The penalty here is the inverse variance of the random effect, so REML estimation here is true REML estimation of variance components. This trick turns gam() into a competent mixed-effects modeller, with the bonus that smooth terms and random effects coexist naturally.

Explanation: te() builds marginal smooths in wt and hp and combines them with an outer-product penalty so the result is a single 2D surface that handles different units cleanly. Prefer te() over s(wt, hp) when the two predictors are measured in incomparable units (kg vs hp). Use s(x1, x2) (an isotropic thin-plate smooth) only when both predictors share a metric, like spatial x-y coordinates.

Explanation: ti(wt, hp) is a tensor-product interaction-only smooth that excludes the main effects (which live in the s() terms). The ti() p-value of 0.34 says the pure 2D interaction is not significant once we account for the marginal smooths. This is the GAM equivalent of testing for an interaction in linear regression by comparing models with and without the cross-term.

Explanation: family = binomial switches the link to logit and the deviance to binomial, identical to glm(), but each predictor can be smooth. The s(wt) edf of 2.34 says the log-odds of manual transmission is a curved function of weight (manuals dominate at low weight, automatics at high weight, with a sharp transition near 3000 lb). Linear logistic regression would have masked that S-shape on the logit scale.

Explanation: Poisson GAMs use a log link, so coefficients are interpreted on the log-rate scale: wool B has exp(-0.206) = 0.81, an 19% lower break rate than wool A. We cap k = 3 because tension only takes three levels (L, M, H) and a higher basis would be unidentified. For real count data check for overdispersion: if the residual deviance is much larger than residual df, switch to family = quasipoisson or family = nb().

Explanation: The Gamma family with log link models the mean of a positive-skew distribution multiplicatively, like Poisson regression but with continuous outcomes. Compared to a Gaussian GAM on the raw Ozone, this respects the lower bound at zero and the long upper tail. A common alternative is family = tw() (Tweedie) which adds a point mass at zero for outcomes that mix zeros with positive values.

Explanation: Point estimates match the Poisson model but the standard errors are larger (0.077 vs 0.052 for woolB) because the negative binomial honestly accounts for over-dispersion. The estimated theta of 13.83 quantifies how much variance exceeds the mean. Rule of thumb: if the Poisson residual deviance is more than twice the residual df, switch to nb() or quasipoisson before trusting the standard errors.

Explanation: k-index near 1 with a p-value comfortably above 0.05 means residual pattern is consistent with k being large enough. A k-index well below 1 with a tiny p-value is the signal to raise k. The k' column is k - 1 (the maximum possible edf). Always check k.check() before trusting a GAM: an under-specified basis silently biases your smooth toward linearity.

Explanation: The standard error swells at the edges (Temp = 60 and Temp = 95) where the data thin out: this is the GAM honestly admitting it has less information out there. For approximate 95% pointwise intervals use fit +/- 1.96 * se_fit. For non-Gaussian families, also pass type = "response" to get predictions on the natural scale (probabilities, rates, expected counts) rather than the linear predictor.

Explanation: AIC drops by 26 from linear to single-smooth (clear evidence the Temp effect is nonlinear) and another 34 with the Wind smooth (clear evidence Wind is a useful second predictor). For nested GAMs estimated by REML you can also use anova(m1, m2, test = "Chisq") for a likelihood-ratio test. For very different model classes (e.g. lm vs gam) stick to AIC since the degrees-of-freedom accounting is uniform there.