Chapter 6 asked whether individual observations change a fitted model. A different problem arises when many observations follow a pattern that the model misses. We begin with pictures and numbers: when should we change the predictor scale, and when should we change the response scale? After choosing and checking a plausible model form, we introduce variable selection: which candidate terms should remain in the model?
Learning Objectives
By the end of this chapter, students should be able to:
- distinguish a curved relationship from a model nonlinear in its parameters;
- use curvature and changing residual spread to motivate a transformation;
- distinguish transforming a predictor from transforming the response;
- interpret predictions after transforming the response;
- carry out forward selection, backward elimination, and stepwise selection;
- explain ridge and lasso, and choose a penalty using cross-validation;
- separate model development from final prediction assessment.
Start with a fitted straight line and ask two questions: does it follow the center of the responses, and is the spread around that center roughly constant? A systematic curve points to a mean-shape problem. A widening spread points to a variance problem. A transformation may help with one or both, but the choice must match the pattern and the meaning of the variables. NIST’s transformation examples illustrate these two goals.
The following datasets are deliberately constructed to make the differences easy to see. Each has 15 observations, three at each of five predictor values; they illustrate model behavior, rather than provide evidence about a population.
Example 7.1 (A Curved Mean with Similar Spread: Transform the Predictor) Suppose we measure a response at input levels \(x=1,2,4,8,16\). In the top-left panel, the response increases quickly at first and then more slowly. The straight line misses the center of several groups. In the top-right panel, doubling \(x\) becomes a one-unit step on the \(\log_2 x\) axis, and the same observations follow a straight line.
library(ggplot2)
ch7_curve <- data.frame(x = rep(c(1, 2, 4, 8, 16), each = 3))
ch7_curve$y <- 10 + 4 * log2(ch7_curve$x) + rep(c(-0.6, 0, 0.6), 5)
ch7_curve_raw <- lm(y ~ x, data = ch7_curve)
ch7_curve_logx <- lm(y ~ log2(x), data = ch7_curve)
Read the residuals. With the original line, low-\(x\) responses lie below the line, middle responses tend to lie above it, and high-\(x\) responses fall below it again. That shared pattern is not an isolated outlier. After transforming \(x\), every group has residuals \(-0.6,0,0.6\); no curved mean pattern remains in this constructed example.
Put numbers on the change. The group means increase by 4 whenever \(x\) doubles. The sample SD within each group is 0.6, so the motivation here is the mean shape, not increasing spread.
| x | Responses | Mean response | Fitted: x | Fitted: log2(x) |
|---|---|---|---|---|
| 1 | 9.4, 10, 10.6 | 10 | 12.97 | 10 |
| 2 | 13.4, 14, 14.6 | 14 | 13.94 | 14 |
| 4 | 17.4, 18, 18.6 | 18 | 15.87 | 18 |
| 8 | 21.4, 22, 22.6 | 22 | 19.74 | 22 |
| 16 | 25.4, 26, 26.6 | 26 | 27.48 | 26 |
The estimated equations are
\[ \begin{aligned} \text{Original predictor:}\quad &\hat y=12.000+0.968x,\\ \text{Transformed predictor:}\quad &\hat y=10.000+4.000\log_2 x. \end{aligned} \]
At \(x=4\), the original line predicts 15.87 and the transformed model predicts 18.00, matching the observed group mean. The fitted slope 4 means a doubling of \(x\) adds 4 response units, since \(\log_2(2x)-\log_2(x)=1\).
The training RMSE, \(\sqrt{\mathrm{SSE}/15}\), falls from 2.090 to 0.490. This comparison is meaningful because both fits use the same observations and response scale, and both have two coefficients. It describes training fit, not held-out prediction accuracy. Only the predictor changed; the response and its units did not.
Example 7.2 (Increasing Level and Spread: Transform the Response) Now suppose a positive amount \(y\) is measured at times \(x=0,1,2,3,4\). The middle response doubles at each step, and the three observations spread farther apart in the original units. A straight line with constant error variance is missing both features. Taking logs compresses large values: the top-right panel becomes straight, and the bottom-right residuals have the same spread at every time.
ch7_growth <- data.frame(x = rep(0:4, each = 3))
ch7_growth$y <- 10 * 2^ch7_growth$x * rep(c(0.8, 1, 1.25), 5)
ch7_growth_raw <- lm(y ~ x, data = ch7_growth)
ch7_growth_logy <- lm(log(y) ~ x, data = ch7_growth)
Put numbers on the spread. At \(x=0\), the observations are \(8,10,12.5\); at \(x=4\), they are \(128,160,200\). Both triples use the same multipliers, \(0.8,1,1.25\), around their middle value. The absolute spread grows, but the proportional spread stays the same.
| x | Responses | Mean response | SD of y | SD of log(y) |
|---|---|---|---|---|
| 0 | 8, 10, 12.5 | 10.17 | 2.25 | 0.223 |
| 1 | 16, 20, 25 | 20.33 | 4.51 | 0.223 |
| 2 | 32, 40, 50 | 40.67 | 9.02 | 0.223 |
| 3 | 64, 80, 100 | 81.33 | 18.04 | 0.223 |
| 4 | 128, 160, 200 | 162.67 | 36.07 | 0.223 |
Each SD is computed from the three observations at that time, not from all 15 pooled together. The original-scale SD rises from 2.25 to 36.07, while the log-scale SD stays at 0.223. A nonlinear mean alone would follow the upward curve but would not remove this increasing original-scale variation.
R gives the fitted equations
\[ \begin{aligned} \text{Original response:}\quad &\hat y=-10.167+36.600x,\\ \text{Log response:}\quad &\widehat{\log y}=2.303+0.693x. \end{aligned} \]
The original line even predicts a negative amount at \(x=0\). The log fit has a different interpretation: one time step adds approximately \(0.693=\log 2\) on the log scale, so exponentiating the fitted value gives the curve \(10\times2^x\). At \(x=2\), that curve gives 40, whereas the observed arithmetic mean is 40.67. Exponentiating an estimated log mean is not the same as estimating an arithmetic mean; Section 7.3 explains the required back-transformation and its assumptions.
Here log() means the natural logarithm; log2() in the previous example made the doubling interpretation especially direct. Do not compare the raw SSE or RMSE of \(y\) with that of \(\log y\): the residuals have different units. Compare model adequacy on each fitted scale and, for prediction, evaluate the same target in common original units.
When should we keep the original scale? If the fitted mean is plausible and residuals show roughly constant spread without a systematic pattern, there may be no reason to transform. A skewed histogram of \(Y\) alone is not enough: regression assumptions concern errors conditional on the predictors. An isolated unusual point calls for the investigation in Chapter 6, not an automatic transformation of all observations.
Transformations are candidates, not mandatory fixes. Curvature may instead be represented by polynomial terms or splines; nonconstant variance may be modeled with weighted least squares. If both shape and spread need attention, choose a meaningful response scale and then recheck which predictor functions are needed. NIST discusses transformation and weighting. In real data, always refit and recheck: these tidy constructed examples do not establish normality, independence, or performance on new observations.
Definition 7.1 (Linear in the Parameters) A regression model is linear in its parameters when it can be written as \(Y_i=\beta_0+\sum_{j=1}^d\beta_j h_j(\mathbf x_i)+\varepsilon_i\), where the functions \(h_j\) are specified independently of the unknown coefficients. The design row is \((1\;h_1(\mathbf x_i)\;\cdots\;h_d(\mathbf x_i))\); the response need not be a straight-line function of the original predictors.
| Relationship | Model fitted by least squares | What changes? |
|---|---|---|
| Curvature in \(x\) | \(Y=\beta_0+\beta_1x+\beta_2x^2+\varepsilon\) | Add a polynomial column; retain the response scale. |
| Diminishing change, \(x>0\) | \(Y=\beta_0+\beta_1\log x+\varepsilon\) | Transform a predictor. |
| Exponential growth, \(Y>0\) | \(\log Y=\beta_0+\beta_1x+\varepsilon\) | Transform the response; specify errors on the log scale. |
| Power relationship, \(Y,x>0\) | \(\log Y=\beta_0+\beta_1\log x+\varepsilon\) | Transform both variables. |
In R, write a quadratic as lm(y ~ x + I(x^2), data = dat); I() makes the square arithmetic. For a log-log model, use lm(log(y) ~ log(x), data = dat). An interaction from Chapter 5 is another constructed column, \(x_1x_2\).
\(Y=a\exp(bx)\exp(\varepsilon)\) becomes \(\log Y=\log a+bx+\varepsilon\). In contrast, taking logs of \(Y=a\exp(bx)+\varepsilon\) does not produce that same additive-error model. Likewise, a saturation model such as \(Y=ax/(b+x)+\varepsilon\) is nonlinear in \(a,b\); reciprocating it changes the errors and may be inappropriate. Some problems need nonlinear least squares or another response model, not forced linearization.
Exercise 7.1 (Linear in Which Sense?) Which can be fitted by lm() using known predictor columns?
Solution 7.1. Models 1 and 3 are linear in their coefficients. Model 2 has an unknown exponent inside the predictor function; ordinary lm() cannot estimate all three coefficients as written. Model 3 models the conditional mean of \(\log Y\), not directly the conditional mean of \(Y\).
Motivation: measuring tree volume. Measuring a standing tree’s diameter and height is easier than measuring its timber volume. The trees from Chapter 5 suggest a physical starting point: volume is roughly proportional to cross-sectional area times height, so a power relationship such as \(V\approx aD^bH^c\) may be more natural than an additive straight-line model. The exponents are to be estimated, not fixed by the cylinder analogy. Data documentation.
Apply what the two examples showed.
Example 7.3 (From Tree Geometry to a Log-Linear Model) For the 31 trees, consider
\[ V_i=aD_i^{b}H_i^{c}\exp(\varepsilon_i) \quad\Longrightarrow\quad \log V_i=\beta_0+b\log D_i+c\log H_i+\varepsilon_i, \qquad \beta_0=\log a. \]
Here \(D\) is diameter in inches, \(H\) is height in feet, and \(V\) is volume in cubic feet. We use all 31 trees to illustrate specification and diagnostics; this is not a held-out performance claim.
library(ggplot2)
ch7_trees <- transform(datasets::trees, Diameter = Girth)
ch7_raw <- lm(Volume ~ Diameter + Height, data = ch7_trees)
ch7_log <- lm(log(Volume) ~ log(Diameter) + log(Height), data = ch7_trees)
knitr::kable(data.frame(Term = names(coef(ch7_log)),
Estimate = unname(coef(ch7_log))), digits = 3, row.names = FALSE)| Term | Estimate |
|---|---|
| (Intercept) | -6.632 |
| log(Diameter) | 1.983 |
| log(Height) | 1.117 |
The estimated diameter and height exponents are about 1.983 and 1.117, close to the geometric motivation but not constrained to it. Holding height fixed, multiplying diameter by \(r\) changes the fitted log-volume by \(\hat b\log r\). Under the lognormal model below, the fitted median volume is multiplied by \(r^{\hat b}\). Thus a 1% diameter increase corresponds approximately to a \(\hat b\)% change in that median; the exact percentage is \(100(1.01^{\hat b}-1)\).
The log model is a candidate to examine, not a guaranteed cure. Both panels still deserve scrutiny. With only 31 observations, a large search over powers and interactions would be hard to justify.
Definition 7.2 (Back-Transformation Under Lognormal Errors) If \(\log Y=\eta(\mathbf x)+\varepsilon\) with \(\varepsilon\mid\mathbf x\sim N(0,\sigma^2)\), then
\[ \operatorname{Median}(Y\mid\mathbf x)=\exp\{\eta(\mathbf x)\}, \qquad \mathbb E(Y\mid\mathbf x)=\exp\{\eta(\mathbf x)+\sigma^2/2\}. \]
Exponentiating a fitted log mean therefore estimates a conditional median, not the original-scale mean. A plug-in mean prediction uses an estimate of \(\sigma^2\); it is not an exactly unbiased estimator.
ch7_newtree <- data.frame(Diameter = 15, Height = 80)
ch7_eta <- as.numeric(predict(ch7_log, newdata = ch7_newtree))
knitr::kable(data.frame(
`Median volume` = exp(ch7_eta),
`Mean volume, lognormal correction` = exp(ch7_eta + sigma(ch7_log)^2 / 2),
check.names = FALSE), digits = 2,
caption = "Predicted cubic feet for diameter 15 inches and height 80 feet.")| Median volume | Mean volume, lognormal correction |
|---|---|
| 37.82 | 37.94 |
With homoscedastic errors independent of the predictors, a less distribution-specific alternative is the smearing estimate \(\exp(\hat\eta)\,n^{-1}\sum_i\exp(e_i)\), where \(e_i\) are training log residuals. A single correction can fail if the error distribution changes with the predictors. Estimate any correction within each training fold when validating.
Other powers. For positive responses, Box-Cox considers \(g_\kappa(y)=(y^\kappa-1)/\kappa\) for \(\kappa\ne0\) and \(g_0(y)=\log y\). Thus \(\kappa=1\) leaves the response shape unchanged, \(\kappa=1/2\) gives a square-root scale, and \(\kappa=0\) gives logs. MASS::boxcox(ch7_raw) profiles the transformation parameter for the specified mean model. Treat it as guidance, then refit and inspect; choosing the power is itself model selection. Box-Cox documentation.
Do not rank an lm(y ~ ...) fit against an lm(log(y) ~ ...) fit using their raw SSE, \(R^2\), or unadjusted software AIC. Their responses differ. Compare held-out predictions of the same target in the same units, or use likelihoods expressed on a common response scale with the transformation Jacobian. Logs require positive values; adding an arbitrary constant changes the model. Binary or count responses may instead call for a GLM.
Exercise 7.2 (What Does Exponentiating Predict?) Suppose \(\log Y=2+0.4x+\varepsilon\), with \(\varepsilon\sim N(0,0.25)\). At \(x=0\), give the conditional median and mean of \(Y\). Does increasing \(x\) by one add 0.4 units to the original-scale mean?
Solution 7.2. The median is \(e^2\approx7.39\) and the mean is \(e^{2+0.25/2}\approx8.37\). A one-unit increase in \(x\) multiplies both by \(e^{0.4}\approx1.49\) because the log-error variance is constant. It does not add 0.4 response units.
Recheck influence after transforming. A new response scale changes residuals, and new predictor columns can change leverage. Apply the diagnostics from Chapter 6 to the fitted log model as well.
Example 7.4 (Does One Tree Drive the Power Relationship?) Return to the log-volume model. A wrongly recorded diameter could distort a timber prediction rule, but an accurately measured unusual tree may reveal where the rule needs improvement. First identify observations worth inspecting:
ch7_diag <- data.frame(Tree = seq_len(nrow(ch7_trees)),
Studentized = rstudent(ch7_log), Leverage = hatvalues(ch7_log),
Cook = cooks.distance(ch7_log))
ch7_top <- head(order(ch7_diag$Cook, decreasing = TRUE), 3)
| Tree | Studentized | Leverage | Cook |
|---|---|---|---|
| 18 | -2.326 | 0.126 | 0.224 |
| 17 | 1.598 | 0.116 | 0.106 |
| 11 | 1.703 | 0.072 | 0.070 |
Refit to understand sensitivity. Temporarily omit the largest-\(D_i\) tree and compare the same coefficients; keep the original fit as the reference.
ch7_flag <- which.max(ch7_diag$Cook)
ch7_log_without <- update(ch7_log, data = ch7_trees[-ch7_flag, ])
knitr::kable(data.frame(Term = names(coef(ch7_log)),
All_trees = unname(coef(ch7_log)),
Omit_flagged_tree = unname(coef(ch7_log_without))),
digits = 3, row.names = FALSE,
col.names = c("Term", "All trees", "Omit flagged tree"))| Term | All trees | Omit flagged tree |
|---|---|---|
| (Intercept) | -6.632 | -7.185 |
| log(Diameter) | 1.983 | 1.958 |
| log(Height) | 1.117 | 1.261 |
Tree 18 has \(D_i\approx0.224\). Omitting it changes the diameter exponent from 1.983 to 1.958, and the height exponent from 1.117 to 1.261. The estimated height contribution is more sensitive here; influence can matter differently for different scientific questions. This calculation supplies no evidence of a recording error. Report the sensitivity and investigate the measurement before deciding how to handle it.
Why recheck before selection? One influential observation can make a predictor appear useful or change which term a search retains. Resolve data-quality questions and inspect model form first, then assess whether the selection is stable. Recheck the final selected model as well. Ridge and lasso still use squared-error loss here; shrinkage is not an automatic remedy for outliers. For prediction assessment, all outcome-driven diagnostic and model-development decisions belong within training data, not the held-out test set.
Transformation asks how a variable should enter the model; variable selection asks which candidate terms to retain. In Section 7.1, replacing \(x\) by \(\log_2 x\) kept the same predictor but changed its relationship with the response. A selection method choosing among only untransformed straight-line terms could not discover that relationship.
Now suppose several predictors are available. Adding columns cannot increase training SSE, even when they represent noise, but estimating more coefficients can make predictions less stable. Selection balances fit against complexity within a scientifically defensible candidate set. Specify plausible transformations and required terms first; if their choice is data-driven, include it in the validation procedure in Section 7.6.
Definition 7.3 (Variable Selection) Variable selection chooses a subset of candidate predictors or model terms using a specified rule. A term may contain several columns, such as a factor. The candidate set, terms that must remain, comparison criterion, and stopping rule are part of the procedure.
Choose the rule before the search. For Gaussian models fitted to the same observations and response, \(\mathrm{AIC}=-2\hat\ell+2k\) and \(\mathrm{BIC}=-2\hat\ell+k\log n\), where \(k=p+1\) includes \(p\) regression coefficients and the variance parameter. Smaller values are preferred. Ordinary \(R^2\) alone cannot stop a search because it never decreases when columns are added. These criteria are scores, not significance levels. R documentation.
| Method | Start | One step | Stop |
|---|---|---|---|
| Forward selection | Required terms, often just the intercept | Try each admissible addition; take the greatest score improvement. | No addition improves the score. |
| Backward elimination | Full candidate model | Try each admissible removal; take the greatest improvement. | No removal improves the score. |
| Stepwise selection | A specified starting model | Reconsider additions and removals as the fit changes. | Neither move improves the score. |
These are greedy searches; different starting points or directions can yield different models. Best subset selection evaluates every allowed subset and finds the best score within that finite set, but \(d\) optional terms can require \(2^d\) fits. Backward elimination requires an estimable starting model with positive residual degrees of freedom for the usual Gaussian criterion.
Keep treatment assignments, blocking variables, and scientifically required adjustment terms when the question requires them. Retain parent terms with interactions and lower powers with polynomial terms. Search a multi-level factor as a whole term. R’s formula-based step() respects interaction hierarchy, but it will not infer every scientific constraint, such as retaining x whenever I(x^2) is present.
Example 7.5 (Which Sensor Measurements Help Predict an Outcome?) Imagine a manufacturing outcome and eight candidate sensor readings. We simulate data so the useful predictors are known: \(Y=5+3x_1-2x_2+1.5x_3+\varepsilon\), with independent \(N(0,4)\) errors. Predictors \(x_1,x_2\) are correlated; \(x_4,\ldots,x_8\) have no additional mean contribution. All methods receive the same eight candidates. This is a teaching simulation, not an industrial dataset.
set.seed(8561)
ch7_n <- 160
ch7_X <- matrix(rnorm(ch7_n * 8), nrow = ch7_n)
ch7_X[, 2] <- 0.7 * ch7_X[, 1] + sqrt(1 - 0.7^2) * ch7_X[, 2]
colnames(ch7_X) <- paste0("x", 1:8)
ch7_dat <- data.frame(y = 5 + 3 * ch7_X[, 1] - 2 * ch7_X[, 2] +
1.5 * ch7_X[, 3] + rnorm(ch7_n, sd = 2), ch7_X)
ch7_id <- sample(seq_len(ch7_n), 120)
ch7_train <- ch7_dat[ch7_id, ]
ch7_test <- ch7_dat[-ch7_id, ]
ch7_folds <- sample(rep(1:5, length.out = nrow(ch7_train)))We reserve 40 observations now and fit all searches on the other 120. The candidate set contains only main effects, so no polynomial or interaction hierarchy needs to be imposed in this example.
See the first forward step. Start with an intercept-only model. Fit each of the eight one-predictor alternatives on the same 120 training observations; the table shows the starting model and the three best additions. Use AIC() consistently, and choose the lowest score only if it improves on the starting score.
ch7_null <- lm(y ~ 1, data = ch7_train, na.action = na.fail)
ch7_full <- lm(y ~ ., data = ch7_train, na.action = na.fail)
ch7_one_term <- lapply(paste0("x", 1:8), function(term)
lm(reformulate(term, response = "y"), data = ch7_train, na.action = na.fail))
ch7_first_scores <- data.frame(
Model = c("Intercept only", paste("Add", paste0("x", 1:8))),
AIC = c(AIC(ch7_null), vapply(ch7_one_term, AIC, numeric(1))))
ch7_first_scores$Change <- ch7_first_scores$AIC - AIC(ch7_null)
ch7_best_additions <- order(ch7_first_scores$AIC[-1]) + 1
knitr::kable(ch7_first_scores[c(1, head(ch7_best_additions, 3)), ],
digits = 2, row.names = FALSE, align = "c",
col.names = c("Candidate model", "AIC", "Change from starting AIC"))| Candidate model | AIC | Change from starting AIC |
|---|---|---|
| Intercept only | 646.33 | 0.00 |
| Add x1 | 604.57 | -41.75 |
| Add x3 | 623.40 | -22.93 |
| Add x7 | 647.05 | 0.72 |
The first move adds x1. AIC falls from 646.33 to 604.57. The next step must compare additions conditional on that retained term; it cannot simply reuse the remaining one-predictor rankings.
Continue the searches.
ch7_scope <- list(lower = formula(ch7_null), upper = formula(ch7_full))
ch7_forward <- step(ch7_null, scope = ch7_scope,
direction = "forward", trace = 0)
ch7_backward <- step(ch7_full, scope = ch7_scope,
direction = "backward", trace = 0)
ch7_both <- step(ch7_null, scope = ch7_scope,
direction = "both", trace = 0)
ch7_searches <- list(Forward = ch7_forward, Backward = ch7_backward,
Stepwise = ch7_both)
knitr::kable(data.frame(
Search = names(ch7_searches),
Selected = vapply(ch7_searches, function(m)
paste(attr(terms(m), "term.labels"), collapse = ", "), character(1)),
AIC = vapply(ch7_searches, AIC, numeric(1))), digits = 2, row.names = FALSE)| Search | Selected | AIC |
|---|---|---|
| Forward | x1, x2, x3, x4 | 494.56 |
| Backward | x1, x2, x3, x4 | 494.56 |
| Stepwise | x1, x2, x3, x4 | 494.56 |
Here all three searches retain \(x_1,x_2,x_3,x_4\). The first three generate the mean, but \(x_4\) is noise: selection can retain a chance association.
step() uses AIC by default; k = log(nrow(ch7_train)) changes the penalty to the BIC-style rule. Fix the analysis rows across all candidates, especially when missing values occur. The trace uses extractAIC(), whose additive constant differs from AIC() for lm; compare scores computed with the same function. Selection documentation.
Exercise 7.3 (Take One Forward Step) The current AIC is 210. Adding \(x_2\), \(x_3\), or \(x_4\) gives 207, 209, or 212. Which move is made? At the next step, the best available AIC is 207.4. Continue or stop? Would a different answer prove that one variable has no scientific role?
Solution 7.3. Add \(x_2\), producing the largest AIC decrease. Stop at the next step because 207.4 is worse than 207. This is a conditional search result under the chosen candidates and score, not proof that an excluded variable has no scientific role.
Subset selection makes discrete include-or-exclude decisions. Regularization instead penalizes coefficient size: it accepts some bias to reduce sensitivity to the sample. Lasso can also produce exact zeros.
Definition 7.4 (Penalized Least Squares) Let \(\mathbf x_i\) contain \(d\) candidate predictor columns, excluding the intercept, and let \(\boldsymbol\beta=(\beta_1,\ldots,\beta_d)\). On \(n_{\mathrm{tr}}\) training observations, estimate the intercept and slopes by minimizing
\[ \frac{1}{2n_{\mathrm{tr}}}\sum_i(y_i-\beta_0-\mathbf x_i^\top\boldsymbol\beta)^2 +\lambda\left\{\frac{1-\alpha}{2}\sum_j\beta_j^2 +\alpha\sum_j|\beta_j|\right\}. \]
The intercept is not penalized. Ridge uses \(\alpha=0\), lasso uses \(\alpha=1\), and elastic net uses \(0<\alpha<1\). Larger \(\lambda\) puts more weight on shrinkage. At \(\lambda=0\) the objective is OLS; uniqueness still depends on the design rank.
Why can lasso give zeros? Near zero, a squared penalty becomes nearly flat, while an absolute-value penalty keeps a nonzero slope on each side of its corner. A small improvement in fit may therefore fail to move a lasso coefficient away from zero. Ridge usually reduces its magnitude continuously instead.
| Method | Main behavior | Interpretation |
|---|---|---|
| Ridge | Shrinks slopes, usually without exact zeros. | Useful when many correlated columns carry information. |
| Lasso | Can set slopes exactly to zero. | A sparse prediction rule; correlated predictors can compete for selection. |
| Elastic net | Combines both penalties. | Can retain correlated predictors together; both \(\alpha\) and \(\lambda\) need a selection rule. |
Why scaling matters. Converting a predictor from meters to millimeters changes its coefficient by a factor of 1,000. Standardize predictor columns using training data so the penalty does not favor a unit convention. In glmnet, standardize = TRUE performs this within each fit, including CV training folds, while predictions use the original input units. An ordinary lasso selects columns; it does not automatically preserve an entire factor or interaction hierarchy. Package documentation.
Choose the penalty on training data. Divide the 120 training observations into five folds. Fit on four and score predictions on the fifth, repeating for all folds. lambda.min minimizes average CV error; lambda.1se is the largest penalty within one estimated standard error of the minimum. The latter is a simplicity rule, not a hypothesis test. A single cv.glmnet() call tunes \(\lambda\) at a fixed \(\alpha\). CV documentation.
Install glmnet once with install.packages("glmnet") if needed. Continue the same sensor example:
library(glmnet)Loading required package: MatrixLoaded glmnet 5.1ch7_xtrain <- as.matrix(ch7_train[paste0("x", 1:8)])
ch7_xtest <- as.matrix(ch7_test[paste0("x", 1:8)])
ch7_ridge <- cv.glmnet(ch7_xtrain, ch7_train$y,
alpha = 0, foldid = ch7_folds, standardize = TRUE, type.measure = "mse")
ch7_lasso <- cv.glmnet(ch7_xtrain, ch7_train$y,
alpha = 1, foldid = ch7_folds, standardize = TRUE, type.measure = "mse")
From one column to many. Ridge remains computable with more columns than observations; ordinary full-model OLS then cannot have full column rank. Lasso can give sparse fits in that setting, but it does not guarantee recovery of the true variables. Here we focus on a modest candidate set so the classical searches and penalized fits are directly comparable.
Our choices are now fixed: AIC for the three subset searches and the one-SE rule for ridge and lasso. Evaluate all five on the same untouched 40 observations. Add full-model OLS and a training-mean baseline for context. For \(m\) test observations, \(\mathrm{RMSE}=\{m^{-1}\sum_i(y_i-\hat y_i)^2\}^{1/2}\), in response units.
ch7_predictions <- cbind(
Baseline = rep(mean(ch7_train$y), nrow(ch7_test)),
Full_OLS = predict(ch7_full, newdata = ch7_test),
Forward = predict(ch7_forward, newdata = ch7_test),
Backward = predict(ch7_backward, newdata = ch7_test),
Stepwise = predict(ch7_both, newdata = ch7_test),
Ridge = as.numeric(predict(ch7_ridge, newx = ch7_xtest, s = "lambda.1se")),
Lasso = as.numeric(predict(ch7_lasso, newx = ch7_xtest, s = "lambda.1se")))
knitr::kable(data.frame(Method = colnames(ch7_predictions),
Test_RMSE = sqrt(colMeans((ch7_predictions - ch7_test$y)^2))),
digits = 3, row.names = FALSE)| Method | Test_RMSE |
|---|---|
| Baseline | 2.972 |
| Full_OLS | 2.182 |
| Forward | 2.185 |
| Backward | 2.185 |
| Stepwise | 2.185 |
| Ridge | 2.095 |
| Lasso | 2.131 |
ch7_lasso_coef <- as.matrix(coef(ch7_lasso, s = "lambda.1se"))[-1, 1]
knitr::kable(data.frame(Predictor = names(ch7_lasso_coef),
Lasso_coefficient = unname(ch7_lasso_coef)), digits = 3, row.names = FALSE)| Predictor | Lasso_coefficient |
|---|---|
| x1 | 2.562 |
| x2 | -1.445 |
| x3 | 1.267 |
| x4 | -0.078 |
| x5 | 0.000 |
| x6 | 0.000 |
| x7 | 0.000 |
| x8 | 0.000 |
Ridge and lasso have slightly smaller test RMSE than full OLS on this split. Lasso keeps a small coefficient for the noise predictor \(x_4\), while setting \(x_5,\ldots,x_8\) to zero. Equal search results are possible, and no method must win on every split. Report the comparison without revising the methods after looking at test outcomes. If choosing a winning method from this table, its same test error is no longer an independent final assessment of that choice.
If CV is used to evaluate a selection procedure, rerun the entire variable search inside each training fold. Selecting on all observations first and then cross-validating only the selected formula leaks information. Learned transformations, imputation, scaling, and tuning follow the same rule. Final performance claims require an untouched test set or appropriately nested resampling.
Prediction and inference answer different questions. A selected coefficient is not automatically a scientific discovery; ordinary fixed-model \(t\) tests and confidence intervals do not automatically remain valid after selection or lasso followed by OLS refitting. Likewise, a good predictive score does not establish a causal effect.
Exercise 7.4 (Find the Leakage) A student selects predictors using all 160 outcomes, then splits the observations and reports test RMSE. Another student reserves 40 observations first and runs selection only on the other 120. Which analysis gives an honest held-out assessment of its fitted procedure? Does that make its selected-coefficient p-values valid?
Solution 7.4. Only the second keeps test outcomes out of model development, provided the test results are not used to revise the procedure. That separation assesses prediction; it does not repair ordinary post-selection coefficient inference.
| What do you see or need? | What should you try? | What must you check? |
|---|---|---|
| Curved mean, roughly constant spread | A predictor transformation or suitable polynomial terms; keep \(Y\) in its units when appropriate. | Does the residual mean pattern disappear? Interpret the new predictor scale. |
| Positive responses with proportional variation | A log response, possibly together with transformed predictors. | Recheck shape and spread; distinguish the original-scale mean from an exponentiated log prediction. |
| Many plausible candidate terms | Forward, backward, or stepwise selection with a stated score and required terms. | Preserve factors and hierarchy; selected terms and ordinary p-values are not automatic discoveries. |
| Unstable estimates or many correlated columns | Ridge for shrinkage; lasso when sparsity is also useful; consider elastic net. | Scale within training data and tune penalties by CV. |
| A final claim about prediction | Evaluate the entire fixed procedure on untouched observations. | Use the same target and units; learn transformations and selection within training folds. |
A transformation is unnecessary when the existing model is adequate for the purpose. A curved relationship can still be linear in its coefficients. Whatever the chosen form, recheck the Chapter 6 diagnostics.
The next design unit asks which terms must be retained because of randomization, blocking, and the scientific question. That provides a useful boundary for automatic selection. Dose-specific comparisons from Chapter 5 are retained in Contrasts and Estimable Comparisons (Optional).
Exercise 7.5 (Practice: Specify, Select, and Assess)
lambda.min with lambda.1se. Where should the folds, transformation choices, and final test observations enter the procedure?Solution 7.5.
lambda.min minimizes CV error; lambda.1se prefers the largest penalty within one estimated SE of that minimum. Choose transformations and tune models within training data/folds, then evaluate the fixed procedure on an untouched test set.