10. Remedial measures and transformations

A biologist lines up the brain and body weights of 62 land mammals, from a 3-gram mouse to a five-ton African elephant, and plots one against the other. The picture is almost useless. Sixty of the points are crushed into the bottom-left corner, the elephant and a few whales float off alone in the top right, and any straight line you draw is dragged around by the giants. Figure 1 is that plot. The spread of the points also grows with body size, so the tidy constant-variance assumption from Chapter 2 is nowhere in sight.

Nothing is wrong with the data. What is wrong is the scale. Take the logarithm of both axes and the same 62 points snap into a clean, straight, evenly scattered band, shown in Figure 2. The relationship was always simple; it was simple in percentage terms, not in raw grams and kilograms, and the log scale is where percentages live. This is the central move of the chapter. When the diagnostics of Chapter 9 fire, you rarely need to abandon regression. More often you need to fit it to a transformed version of the problem, or weight the observations, or replace the loss function, so that the model’s assumptions become true.

Chapter 9 taught you to detect trouble: curved residual plots, funnel shapes, heavy tails, one island or one country bending the whole fit. This chapter is the toolbox of fixes. You will read log coefficients as percentages, pick a transformation with the Box-Cox likelihood, weight observations that carry unequal noise, and resist a stubborn outlier with robust regression. At the end you will meet a dataset that no transformation can save, and learn to recognize when the honest fix is a different model entirely.

10.1 A menu of remedies¶

Before any single technique, it helps to see the whole menu and the logic that picks among its items. Chapter 9 ended by naming problems, from curved residual plots and funnel-shaped scatter to the high-leverage cases of 9.3 Influence: which points actually change the fit. Each named problem has a matching family of fixes.

If the residual plot curves, the mean function is wrong: a straight line is being asked to trace a bend. The fix is to bend the model back straight, either by transforming the predictor $X$ (leaving the error alone) or by transforming the response $Y$. If the residual plot fans out, the variance is not constant: the model assumes one $\sigma^2$ but the data carry more noise in some places than others. The fix is either a variance-stabilizing transformation of $Y$ or weighted least squares, which tells the fit to trust the low-noise points more. If the residuals have heavy tails or a lone outlier, one or two points are dominating the squared-error loss. The fix is a loss that grows more slowly than the square, which is robust regression. And if the response is a count, a proportion, or a strictly positive amount with structural skew, the normal-errors model may be the wrong frame altogether, and the fix is a different model, which is where Chapters 13 and 14 go.

The whole logic fits on one page. Figure 3 lays it out as a decision chart: start at the residual plot, read the pattern, and follow the arrow to the fix. Keep it beside you for the rest of the chapter, which walks the four columns from left to right.

Transformations of the predictor and the response do different jobs, and it is worth keeping them apart. Transforming $X$ changes the shape of the mean curve without touching the errors, so it is the tool for a curved relationship whose scatter is already even. Transforming $Y$ changes both the shape and the spread at once, so it is the tool when a curve and a funnel appear together, exactly as in the mammal data. A rough guide, drawn in Figure 4, is the “ladder of powers”: if a plot bends one way, step down the ladder from $Y$ toward $\sqrt{Y}$, $\log Y$, $1/Y$; if it bends the other way, step up toward $Y^2$. The Box-Cox method in 10.3 Box-Cox: letting the data choose the power turns this guessing into an estimate.

The rest of the chapter takes these remedies one at a time. Throughout, keep the Chapter 2 model in view: $Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i$ with errors that should average to zero, share one variance $\sigma^2$, and be uncorrelated. Every remedy here is an effort to make one of those three conditions true after a diagnostic showed it was false.

10.2 The log transformation and reading its coefficients¶

The log transformation earns its own section because it is the most common remedy and the one most often misread. We can now fit the mammal line; the harder question is what its slope of 0.75 actually means, because a coefficient on a logged variable is not a change in grams per kilogram. It is a percentage.

Intuition¶

Logarithms turn multiplication into addition, so they turn “grows by a percentage” into “grows by a constant.” A quantity that doubles every time some driver doubles looks curved on ordinary axes and straight on log axes. Brain weight relates to body weight that way: across species, a body ten times heavier tends to carry a brain a fixed multiple larger, not a fixed number of grams larger. That is why the raw plot curved and the log-log plot is straight. The price of the straight line is that the slope now speaks in percents, and you have to translate.

Formula¶

Two log models cover almost every case. The log-log model logs both sides,

and the log-linear model logs only the response,

Here $\log$ is the natural logarithm (base $e$) throughout this book. The slope means something different in each model, and each meaning carries a name (Definition 10.1, Definition 10.2).

In words: on the log-log scale the slope answers “if $X$ rises 1%, about how many percent does $Y$ rise,” and on the log-linear scale it answers “if $X$ rises by one of its own units, about how many percent does $Y$ rise.” Both readings have an exact form and an approximate form, and the next derivation produces both.

Derivation (exact and approximate percent readings)¶

Proof. Undo the log on the response. The log-log model says $\log Y = \beta_0 + \beta_1 \log X$ for the mean, so exponentiating,

In words: the mean of $Y$ is a constant times a power of $X$. Now multiply $X$ by a factor $c$ (for a 1% increase, $c = 1.01$; for a doubling, $c = 2$). The new mean divided by the old is

So multiplying $X$ by $c$ multiplies $Y$ by $c^{\beta_1}$, exactly. The exact percent change in $Y$ for a 1% rise in $X$ is $100\,(1.01^{\beta_1} - 1)$. For the approximate reading, use $c^{\beta_1} = e^{\beta_1 \log c} \approx 1 + \beta_1 \log c$ when $\log c$ is small; with $c = 1.01$, $\log c \approx 0.01$, giving a percent change of about $100 \cdot \beta_1 \cdot 0.01 = \beta_1$ percent. That is the elasticity headline: on the log-log scale, $\beta_1$ is roughly the percent change in $Y$ per percent change in $X$. $\blacksquare$

Proof (log-linear semi-elasticity). In the log-linear model the mean satisfies $\log Y = \beta_0 + \beta_1 X$, so $Y = e^{\beta_0} e^{\beta_1 X}$. Raise $X$ by one unit:

So each one-unit step in $X$ multiplies $Y$ by $e^{\beta_1}$, exactly, giving an exact percent change of $100\,(e^{\beta_1} - 1)$. Since $e^{\beta_1} \approx 1 + \beta_1$ for small $\beta_1$, the approximate percent change is $100\,\beta_1$ percent per unit. $\blacksquare$

The gap between exact and approximate widens as the coefficient grows, drawn in Figure 5. For a coefficient near $\pm 0.1$ the two agree to a rounding error, which is why the “percent per unit” shortcut is safe for small effects and misleading for large ones. This same $e^{\beta}$ machinery returns twice later: as the odds ratio in logistic regression (13.3 Reading coefficients as odds ratios) and as the percent growth per month in the logged airline series (15.5 Forecasting honestly: test on the future).

R¶

The mammal fit is a log-log model, so its slope is an elasticity.

mammals <- read.csv("data/mammals.csv")
dim(mammals)
head(mammals, 3)

[1] 62  3
          species  body brain
1      Arctic fox 3.385  44.5
2      Owl monkey 0.480  15.5
3 Mountain beaver 1.350   8.1

The fastest way to feel why an elasticity is a percentage is to set one yourself and read the translation as it changes.

What to notice. The exact reading $100(1.10^{b_1} - 1)$ and the shortcut $10\,b_1$ never part company by more than a tenth of a point anywhere on the slider, which is exactly what the derivation promises for a small coefficient. Try this. Push $b_1$ down from 1 until SSE stops falling, land near 0.75, and read off the percent that 10.2 The log transformation and reading its coefficients asks you to report.

10.3 Box-Cox: letting the data choose the power¶

Guessing a transformation from the ladder of powers works, but it feels arbitrary, and two analysts can disagree. Box and Cox proposed a principled fix: put the candidate powers into one family indexed by a number $\lambda$, then let maximum likelihood pick the $\lambda$ that makes the normal model fit best. We can now motivate the method with a dataset that curves and fans at once: the classic cars data, where a car’s stopping distance is measured against its speed.

Intuition¶

The powers $\sqrt{Y}$, $\log Y$, $Y$, and $Y^2$ are all special cases of raising $Y$ to a power. If you write the transformation as one smooth function of a power parameter $\lambda$, then choosing a transformation becomes choosing a number, and choosing a number is something likelihood does well. Box-Cox writes down the likelihood of the data as a function of $\lambda$, with the regression coefficients and error variance profiled out, and reports the $\lambda$ that maximizes it, plus a confidence interval so you know how sharply the data pin it down.

Formula¶

• $\lambda$ is the power parameter to be estimated; $\lambda = 1$ is essentially no transformation, $\lambda = 0.5$ is the square root, $\lambda = 0$ is the log, $\lambda = -1$ is the reciprocal.

• The “-1” and division by $\lambda$ make the family continuous at $\lambda = 0$, where it becomes $\log Y$, so the ladder has no gap.

In words: one dial $\lambda$ slides smoothly through all the rungs of the ladder of powers, and the log sits exactly at the $\lambda = 0$ notch.

Derivation (Box-Cox as a profile likelihood)¶

Proof. Assume that for the correct power the transformed response follows the normal linear model,

where $\mathbf{x}_i'$ is the $i$-th row of the design matrix $\mathbf{X}$ from 7.1 The model and least squares in matrix form. In words: once $Y$ is raised to the right power, the result obeys the ordinary normal regression model. The parameters are $\lambda$, $\boldsymbol{\beta}$, and $\sigma^2$. The subtle point is that the likelihood must be written for the original data $Y_i$, not for the transformed $Y_i^{(\lambda)}$, because we are comparing different values of $\lambda$ and each puts $Y^{(\lambda)}$ on a different scale. Changing variables from $Y_i^{(\lambda)}$ to $Y_i$ brings in the Jacobian, the derivative factor that any change of variables introduces, here $\mathrm{d}Y_i^{(\lambda)}/\mathrm{d}Y_i = Y_i^{\lambda - 1}$. The log-likelihood of the observed $Y_i$ is therefore

In words: this is the usual normal log-likelihood for the transformed response, plus one extra term that corrects for stretching the $Y$ scale (the sum of the log-Jacobians). Now profile out $\boldsymbol{\beta}$ and $\sigma^2$: for any fixed $\lambda$, the values that maximize $\ell$ are exactly the least-squares fit of $Y^{(\lambda)}$ on $\mathbf{X}$, giving residual sum of squares $\mathrm{RSS}(\lambda)$ and $\hat\sigma^2(\lambda) = \mathrm{RSS}(\lambda)/n$ (the same argument as the normal maximum-likelihood estimate of the variance in Chapter 2). Substituting these back and dropping constants leaves the profile log-likelihood, a function of $\lambda$ alone:

In words: for each $\lambda$, fit the transformed response by ordinary least squares, read off its residual sum of squares, penalize by the Jacobian term, and the best $\lambda$ is the one that maximizes the result. The maximizer $\hat\lambda$ is the Box-Cox estimate. A $100(1-\alpha)$% confidence interval is the set of $\lambda$ with $\ell_p(\lambda) \ge \ell_p(\hat\lambda) - \tfrac{1}{2}\chi^2_{1,\,1-\alpha}$, the usual likelihood-ratio interval. $\blacksquare$

R¶

R’s MASS::boxcox computes $\ell_p(\lambda)$ over a grid and draws it. The plot in Figure 8 is the profile log-likelihood for the cars model, with its peak and 95% interval marked.

cars <- read.csv("data/cars.csv")
fit_lin <- lm(dist ~ speed, data = cars)
coef(fit_lin)

(Intercept)       speed
 -17.579095    3.932409

The straight-line fit is not obviously wrong, but its residuals fan out as the fitted distance grows, the funnel shown in Figure 7. That widening scatter is the signal that the response needs a transformation, and Box-Cox will say which one.

The profile likelihood is easier to trust once you have watched it fight the residual plot for a while, so take the dial yourself.

What to notice. The profile log-likelihood and the shape of the residual plot agree: $\ell_p$ climbs from -135.63 at $\lambda = 1$ to its peak -126.73 at $\lambda = 0.43$ exactly as the funnel flattens, and $\mathrm{RSS}(\lambda)$ alone would have told you nothing, since it falls all the way to 9.56 at the log. Try this. Find the peak by hand, then keep going down past $\lambda = 0.2$ and watch a new low outlier open up, the picture of over-transforming that 10.3 Box-Cox: letting the data choose the power warns about.

10.4 Weighted least squares¶

Transformations attack a curved mean or a variance that grows with the mean. Sometimes, though, the variance is unequal for a reason that has nothing to do with the mean: some observations are simply measured more precisely than others. Averaging ten readings gives a more reliable point than a single reading; a well-calibrated instrument beats a rough one. When you know the relative precision of each point, you should not throw that knowledge away by treating every point equally. Weighted least squares is how you use it.

Intuition¶

Ordinary least squares gives every observation the same vote. If one point carries ten times the variance of another, that is like letting a noisy witness testify as loudly as a careful one. Weighted least squares turns down the volume on the noisy points and turns it up on the precise ones, in exact inverse proportion to their variances. Figure 11 shows the idea: low-noise points are drawn large because they count more.

Formula¶

Keep the linear model but let the error variances differ:

Here $w_i$ is known up to the common constant $\sigma^2$.

In words: penalize each squared residual by its weight, so precise points contribute more to the total the fit is trying to shrink.

Derivation (WLS from Gauss-Markov)¶

Proof. The trick is to transform the unequal-variance model into an equal-variance one and then quote the Gauss-Markov result we already proved. Multiply the $i$-th observation through by $\sqrt{w_i}$. Collecting these into matrices with $\mathbf{W}^{1/2} = \operatorname{diag}(\sqrt{w_1},\dots,\sqrt{w_n})$, define

The transformed model is $\mathbf{Y}^{\ast} = \mathbf{X}^{\ast}\boldsymbol{\beta} + \boldsymbol{\varepsilon}^{\ast}$, and its errors are now homoscedastic:

The starred model satisfies every ordinary least squares assumption, so by the matrix Gauss-Markov theorem (7.6 The Gauss-Markov theorem) the best linear unbiased estimator is ordinary least squares on the starred data:

This is the WLS estimator, and because it is OLS on a model that truly has constant variance, it inherits every good property: it is unbiased and it is BLUE, with $\operatorname{Var}\{\mathbf{b}_W\} = \sigma^2 (\mathbf{X}'\mathbf{W}\mathbf{X})^{-1}$. Plain OLS on the original data is still unbiased here, but it is no longer best: it wastes the precision information in the weights. $\blacksquare$

The derivation also tells you how to compute WLS with nothing but an OLS routine: form $\sqrt{w_i}\,Y_i$ and $\sqrt{w_i}\,\mathbf{x}_i$ and regress. We use that as a check below.

R¶

The strongx data are a physics experiment: at ten beam energies, a scattering cross-section was measured, and each measurement comes with its own known standard deviation sd. Known standard deviations are the textbook case for WLS, with weights $w_i = 1/\text{sd}_i^2$.

strongx <- read.csv("data/strongx.csv")
strongx[, c("energy", "crossx", "sd")]

   energy crossx sd
1   0.345    367 17
2   0.287    311  9
3   0.251    295  9
4   0.225    268  7
5   0.207    253  7
6   0.186    239  6
7   0.161    220  6
8   0.132    213  6
9   0.084    193  5
10  0.060    192  5

When the variances are not known, you estimate the weights, usually by modeling how the spread grows (for example, regressing the absolute residuals on a predictor and using the fitted values to build weights). That two-step method, iterated, is common and effective, and it recycles the variance-modeling instinct that Chapter 14 reuses to handle overdispersion in count models. The strongx case is cleaner because physics handed us the variances directly.

Weights are easiest to believe when you can turn them up gradually and watch the line answer, so slide between ordinary and weighted least squares yourself.

What to notice. At $p = 0$ the weighted line sits exactly on the ordinary one, which is the derivation’s statement that equal weights give back plain least squares. Try this. Raise $p$ to 1 and watch the slope fall from 619.7 to 530.8 while $\mathrm{SE}(b_0)$ drops from 10.08 to 8.08: the numbers of Example 10.3, arriving as a rotation you can see. Keep going to $p = 2$ and ask yourself whether you would still trust a fit resting on so few points (10.4 Weighted least squares).

10.5 When a transformation is not enough: the Galapagos counts¶

Every remedy so far assumed the underlying model was sound and only its scale or its weighting was off. Sometimes that assumption is wrong, and a transformation only hides the problem. The Galapagos species data from Chapter 9 are a case in point, and following them here shows you how to tell a fixable scale problem from a model that needs replacing.

Recall the setup from 9.3 Influence: which points actually change the fit: for 30 Galapagos islands, the number of plant species is regressed on geographic predictors (area, elevation, and distances), and one island, Isabela, has an enormous leverage value because it dwarfs the others in area. The raw fit has a residual funnel, because islands with more species also scatter more widely, the signature of count data. A natural remedy is a variance-stabilizing transformation. For counts, the square root is the classic choice, since a Poisson count with mean $\mu$ has variance $\mu$, and $\sqrt{\text{count}}$ has roughly constant variance.

gala <- read.csv("data/gala.csv")
fit_g_lin  <- lm(Species ~ Area + Elevation + Nearest + Scruz + Adjacent,
                 data = gala)
fit_g_sqrt <- lm(sqrt(Species) ~ Area + Elevation + Nearest + Scruz + Adjacent,
                 data = gala)
round(c(R2_linear = summary(fit_g_lin)$r.squared,
        R2_sqrt   = summary(fit_g_sqrt)$r.squared), 4)

R2_linear   R2_sqrt
   0.7658    0.7827

The square root helps a little, and Figure 14 shows the funnel easing. But look closer and the cracks show. Isabela still has a leverage value of 0.97 on the transformed fit, essentially unchanged, because transforming $Y$ does nothing about an $X$ that is far from the others.

h <- hatvalues(fit_g_sqrt)
names(h) <- gala$island
round(sort(h, decreasing = TRUE)[1:3], 3)

   Isabela Fernandina     Darwin
     0.969      0.950      0.466

Figure 15 shows why the gap opens, using the squaring curve that undoes a square-root fit. Take two predictions equal distances above and below the average. Squaring bends the curve upward, so the high one climbs more than the low one drops. Their average (on the straight chord) lands above the square of the average (on the curve), and that vertical gap is the bias. Any curved back-transform, including the log, does the same thing.

Step back and count the warning signs. The response is a count. Several islands have very few species (six islands have fewer than ten). The variance grows with the mean by construction. Back-transformed predictions are biased, and a plain linear fit to the counts predicts a negative number of species for the smallest islands, which is impossible. Every one of these says the same thing: the normal, constant-variance linear model is the wrong frame, and no power of $Y$ will make it right. The honest fix is not to transform the data until they fit the model, but to change the model to one built for counts. That model is Poisson regression, and 14.1 Why counts break the linear model returns to these exact islands and fits them properly, recalling both Isabela’s high-leverage status from Chapter 9 and this failed transformation. For now, the lesson is diagnostic: learn to tell a scale problem, which a transformation fixes, from a distributional problem, which it cannot.

10.6 A first look at robust regression¶

The last remedy addresses a different failure: not a curved mean or unequal variance, but a handful of points with outsized residuals that drag the fit toward themselves. Least squares is exquisitely sensitive to such points because it squares residuals, so a single large one can dominate the entire sum. Robust regression replaces the square with a loss that grows more gently, so no one point can take over.

Figure 16 shows the payoff on a simple cloud with one point knocked far off. The least-squares line tilts up to chase the stray point; the Huber line shrugs it off and stays with the crowd. That is the entire promise of robust regression in one picture, and the rest of the section explains how the capped loss delivers it.

Intuition¶

Picture the sum of squared residuals as a tug of war in which each point pulls with a force proportional to its residual. A point twice as far off pulls four times as hard, because of the square. Huber’s idea is to cap that escalation: inside a normal range, keep the squared loss and its efficiency, but once a residual is large enough to look like an outlier, let its pull grow only linearly, not quadratically. Figure 17 draws the resulting weight each point receives.

Formula¶

M-estimation, where the M signals a maximum-likelihood-style objective, replaces the least-squares objective with

where $\rho_c$ is the Huber loss.

• $u = (Y_i - \mathbf{x}_i'\boldsymbol{\beta})/s$ is the residual standardized by a resistant scale estimate $s$.

• $c$ is a tuning constant; $c = 1.345$ gives about 95% of the efficiency of OLS when the errors really are normal, while bounding the influence of outliers.

In words: penalize ordinary residuals by the familiar square, but penalize far-out residuals only linearly, so an outlier is expensive rather than catastrophic.

Derivation (IRLS and the honest limitation)¶

Proof. Differentiate the objective and set it to zero. With $\psi_c = \rho_c'$, the estimating equations are $\sum_i \psi_c(u_i)\,\mathbf{x}_i = \mathbf{0}$. Write $\psi_c(u) = w(u)\,u$ with the weight

Then the estimating equations become $\sum_i w(u_i)\,(Y_i - \mathbf{x}_i'\boldsymbol{\beta})\,\mathbf{x}_i = \mathbf{0}$, which are exactly the weighted least squares normal equations from 10.4 Weighted least squares with weights $w(u_i)$. The catch is that the weights depend on the residuals, which depend on the fit. So we iterate: start from the OLS fit, compute residuals and hence weights, solve the WLS problem, recompute residuals and weights, and repeat until the coefficients stop moving. This is iteratively reweighted least squares (IRLS). A residual inside the band ($|u| \le c$) keeps weight 1; a wild residual gets weight $c/|u|$, shrinking as it grows. $\blacksquare$

The derivation exposes the method’s honest boundary. The weight $w(u)$ depends only on the residual, so Huber regression protects against outliers in the response direction, points with a large vertical miss. It does nothing about a point whose predictor values are extreme, a high-leverage point, unless that point also happens to have a large residual. A high-leverage point often pulls the line onto itself and so keeps a small residual, sailing through with full weight. Guarding against that needs a bounded-influence or MM-estimator, which we do not develop here. The savings data make the limitation concrete.

R¶

Fit the Chapter 8 savings model both ways, ordinary and Huber, on the 50 countries.

10.7 Chapter summary¶

You can now respond to a failed diagnostic instead of just naming it. You read a menu of remedies and match each to the pattern that calls for it: transform the predictor for a curved mean with even scatter, transform the response for a curve and a funnel together, weight for known unequal precision, use a slower-growing loss for vertical outliers, and change the model when the response is a count or a proportion. You interpret log coefficients as percentages on both the exact and the approximate scale, choose a power with the Box-Cox profile likelihood and its confidence interval for $\lambda$, derive weighted least squares from the Gauss-Markov argument and compute it from known weights, and read Huber regression as iteratively reweighted least squares while stating its honest limitation. And you can recognize, as with the Galapagos counts, when no transformation will do and a different model is the right answer.

Key results at a glance

Key terms. Remedial measure, log-log model, log-linear model, elasticity, semi-elasticity, Box-Cox family, profile likelihood, Jacobian term, weighted least squares, weight, whitening, back-transformation bias, robust regression, M-estimation, Huber loss, iteratively reweighted least squares.

You should now be able to

• Choose a remedy (transform $X$, transform $Y$, weight, or change the model) from the diagnostic pattern you see.

• Interpret a log-transformed coefficient on both the exact and the approximate percent scale, and derive both readings.

• Derive the Box-Cox profile log-likelihood, including the Jacobian term, and use it to pick a response transformation with a confidence interval for $\lambda$.

• Derive weighted least squares from the Gauss-Markov argument and compute it from known weights.

• Explain and demonstrate Huber robust regression, and state what it protects against and what it does not.

• Decide when a transformation is a patch and a different model is the honest fix.

Where this fits. In the workflow spine of The modeling workflow, this chapter lives at the seam between CHECK and FIT. Chapter 9 (CHECK) told you an assumption was broken; the remedies here send you back to FIT with a better-specified model, a transformed scale, or a weighting scheme, after which you CHECK again and only then USE the result. That loop, diagnose then remedy then re-diagnose, is the working core of applied regression, and it rarely ends after one pass. The threads continue: the log-interpretation machinery of 10.2 The log transformation and reading its coefficients returns as odds ratios in 13.3 Reading coefficients as odds ratios and as monthly growth in 15.5 Forecasting honestly: test on the future, the variance-modeling instinct of 10.4 Weighted least squares returns as overdispersion in 14.6 One family: the generalized linear model, and the Galapagos islands that defeated every transformation here get their honest resolution in 14.1 Why counts break the linear model.

10.8 Frequently asked questions¶

Q1. Should I transform $X$, transform $Y$, or both? Look at the residual plot. A curved mean with even scatter is a predictor problem: transform $X$ and leave the errors alone. A curve together with a funnel is a response problem: transform $Y$, which fixes both at once. When only the variance fans out but the mean is straight, prefer weighting or a variance-stabilizing transform of $Y$. The mammal data needed both variables logged because the relationship was multiplicative in both.

Q2. When is the “percent per unit” shortcut safe? When the coefficient is small in magnitude, roughly under 0.1. Then the exact factor $e^{\beta}$ and the approximation $1 + \beta$ agree to a fraction of a percent (Figure 5). For a coefficient like 0.5 or larger, report the exact percent $100(e^{\beta}-1)$, since the shortcut can be off by several points.

Q3. Box-Cox gave $\hat\lambda = 0.43$. Why did we use $\lambda = 0.5$ instead? Because 0.5 is inside the 95% confidence interval $[0.23, 0.66]$, so the data do not distinguish it from the maximum, and the square root is far easier to interpret and to justify than a raw power of 0.43. Round to a nearby interpretable value whenever the interval allows it.

Q4. My response has zeros or negatives. Can I still use Box-Cox or a log? Not directly: both need $Y > 0$. Common fixes are a small shift ($\log(Y + a)$) or, better, a model built for the data type, since zeros in a count or a proportion usually signal that a generalized linear model (Chapters 13 and 14) is the right tool, not a shifted transform.

Q5. If OLS is unbiased even under unequal variances, why bother with WLS? OLS stays unbiased but stops being efficient: it no longer has the smallest variance, and its reported standard errors are wrong because they assume one common $\sigma^2$. WLS restores both the efficiency and the honest standard errors by using the known precision of each point, as the strongx intercept’s tighter standard error showed.

Q6. Does robust regression replace diagnostics? No. It guards against vertical outliers automatically, but it is blind to high-leverage points (the Libya lesson), and its standard errors are approximate. Use it as one tool alongside the Chapter 9 diagnostics, and always investigate why a point is unusual before deciding to downweight it. A downweighted point is sometimes the most interesting observation in the data.

Q7. How do I get a prediction on the original scale from a transformed model? Undo the transform, but know that naively back-transforming the fitted value is biased, as squaring Isabela’s prediction overshot by 46 species. The mean of a nonlinear function is not that function of the mean. For the log model there is a standard smearing correction; for careful prediction intervals, transform the interval endpoints rather than the point estimate.

10.9 Practice problems¶

1. (A) For each diagnostic pattern, name the remedy this chapter recommends: (i) a curved residual plot with even scatter; (ii) a curve and a funnel together; (iii) a single point with a huge residual; (iv) a count response with many small values.

2. (A) Explain the difference between transforming a predictor and transforming the response in terms of what each does to the mean function and to the error variance.

3. (A) A log-log model has slope $\beta_1 = 1.3$. State in words what a 1% increase in $X$ does to $Y$, and whether $Y$ grows more or less than proportionally.

4. (A) Distinguish an elasticity from a semi-elasticity, and say which log model produces each.

5. (A) In the Box-Cox family, what transformation does each of $\lambda = 1, 0.5, 0, -1$ correspond to? Why is the family written with the “-1” and the division by $\lambda$?

6. (A) Explain in one or two sentences why the Box-Cox profile log-likelihood includes the Jacobian term $(\lambda-1)\sum \log Y_i$, and what goes wrong without it.

7. (A) State the weighted least squares criterion in words, and explain why a point with small variance should get a large weight.

8. (A) Why does ordinary least squares remain unbiased under unequal error variances but stop being the best linear unbiased estimator?

9. (A) Describe the difference between a vertical outlier and a high-leverage point, and say which one Huber regression protects against.

10. (A) Give two features of a response variable that should push you toward a different model rather than a transformation, and name the chapter that supplies that model.

11. (B) Starting from the log-log mean relation $Y = e^{\beta_0} X^{\beta_1}$, derive that multiplying $X$ by $c$ multiplies $Y$ by $c^{\beta_1}$ exactly (Theorem 10.3), and derive the approximate “$\beta_1$ percent per percent” reading.

12. (B) For the log-linear model, derive the exact factor $e^{\beta_1}$ per unit of $X$ and its approximation $1 + \beta_1$, and find the coefficient value at which the approximation understates the true percent change by exactly one percentage point.

13. (B) Derive the Box-Cox profile log-likelihood $\ell_p(\lambda)$ (Theorem 10.6) from the full log-likelihood, showing how $\boldsymbol{\beta}$ and $\sigma^2$ are profiled out and where the Jacobian term comes from.

14. (B) Show that the Box-Cox family is continuous at $\lambda = 0$, that is, $\lim_{\lambda \to 0}(Y^\lambda - 1)/\lambda = \log Y$.

15. (B) Derive the weighted least squares estimator $\mathbf{b}_W = (\mathbf{X}'\mathbf{W}\mathbf{X})^{-1}\mathbf{X}'\mathbf{W}\mathbf{Y}$ (Theorem 10.8) by transforming the model with $\mathbf{W}^{1/2}$ and applying the Gauss-Markov theorem. State the variance of $\mathbf{b}_W$.

16. (B) Show that WLS minimizes $Q_w(\boldsymbol{\beta}) = \sum w_i (Y_i - \mathbf{x}_i'\boldsymbol{\beta})^2$ by differentiating and obtaining the weighted normal equations $\mathbf{X}'\mathbf{W}\mathbf{X}\,\mathbf{b}_W = \mathbf{X}'\mathbf{W}\mathbf{Y}$.

17. (B) For simple linear regression through weighted least squares, derive the closed form $b_{W,1} = \sum w_i (X_i - \bar{X}_w)(Y_i - \bar{Y}_w) / \sum w_i (X_i - \bar{X}_w)^2$, where $\bar{X}_w$ and $\bar{Y}_w$ are the weighted means. Define the weighted means.

18. (B) Starting from the Huber loss $\rho_c$ (Definition 10.10), compute $\psi_c = \rho_c'$ and show that the estimating equations can be written as weighted normal equations with weight $w(u) = \min(1, c/|u|)$ (Theorem 10.11).

19. (B) Explain, using the weight $w(u) = \min(1, c/|u|)$, why a high-leverage point with a small residual is not downweighted by Huber regression. Contrast with a low-leverage point that has a large residual.

20. (B) Suppose row $i$ of a dataset is the average of $n_i$ independent readings each of variance $\tau^2$. Show that the correct WLS weight is $w_i = n_i$, and identify the common constant $\sigma^2$ in the model $\operatorname{Var}\{\varepsilon_i\} = \sigma^2/w_i$.

21. (C) Fit the mammal log-log model in R or Python, report $b_0$, $b_1$, and $R^2$, and translate the slope into the exact percent change in brain weight for a 25% increase in body weight.

22. (C) Using cars.csv, reproduce the Box-Cox estimate $\hat\lambda$ and its 95% interval, then fit both dist ~ speed and sqrt(dist) ~ speed and compare their residual-versus-fitted plots. Which model better satisfies constant variance?

23. (C) Using strongx.csv, fit OLS and WLS with weights $1/\text{sd}^2$, report both slopes and both intercept standard errors, and confirm the WLS coefficients by regressing the $\sqrt{w}$-scaled variables with no intercept.

24. (C) On strongx.csv, add a quadratic term (crossx ~ energy + I(energy^2)) to the WLS fit and compare it to the linear WLS fit. Does the curvature term improve the fit, and what does that say about the linear model?

25. (C) Using gala.csv, fit Species and sqrt(Species) on the five geographic predictors, compare $R^2$ and the residual plots, and report Isabela’s leverage value in both fits. Explain why the transformation does not reduce Isabela’s leverage value.

26. (C) Predict the species count for Fernandina from the square-root gala model, back-transform by squaring, and compare to its actual count. Comment on the direction of the back-transformation bias.

27. (C) Using savings.csv, fit OLS and Huber regression of sr on the four predictors, list the four countries with the smallest Huber weights, and cross-tabulate each against its leverage value. Identify a high-leverage country that Huber does not downweight.

28. (C) Simulate heteroscedastic data (seed 4210): 40 points with $\operatorname{Var}\{\varepsilon_i\}$ proportional to $X_i^2$, fit OLS and WLS with weights $1/X_i^2$, and compare the two slope estimates and their standard errors across 500 replications. Which is more precise, and does either show bias?

10.10 Exam practice¶

These five questions are written in the style of the course exams: each asks you to explain, evaluate, or interpret in full sentences, not just to produce a number. Where a question shows software output, the numbers were produced on the course machine (R 4.6.0) reading the same CSV files from data/ you used all term; Python with statsmodels gives the same values. Read the output, then answer in complete sentences with units. Try each question before opening its model answer.

EP 10.1 (evaluate a claim). A classmate fits a log-linear model of Galapagos plant species on island elevation and reads the slope as a percent effect.

gala <- read.csv("data/gala.csv")
fit <- lm(log(Species) ~ Elevation, data = gala)
round(coef(fit), 6)

(Intercept)   Elevation
   2.591399    0.002490

The classmate writes: “The coefficient 0.00249 means about $0.25\%$ more species per meter of elevation, so an island 100 meters higher has about $25\%$ more species.” Explain what the coefficient means on the percent scale, then evaluate both readings: the per-meter one and the 100-meter one. Which is sound, which is not, and what is the correct figure for a 100-meter difference?

EP 10.2 (interpret output in context). Box-Cox run on the raw mammal regression of brain weight on body weight returns the output below.

mammals <- read.csv("data/mammals.csv")
library(MASS)
bc <- boxcox(lm(brain ~ body, data = mammals),
             lambda = seq(-0.5, 0.5, by = 0.01), plotit = FALSE)
lambda_hat <- bc$x[which.max(bc$y)]
ci <- range(bc$x[bc$y > max(bc$y) - 0.5 * qchisq(0.95, 1)])
round(c(lambda_hat = lambda_hat, ci_low = ci[1], ci_high = ci[2]), 3)

lambda_hat     ci_low    ci_high
      0.07      -0.03       0.18

Interpret $\hat\lambda$ and its interval in context. Which transformation of the response does the output recommend, and why? The chapter’s opening logged both brain and body weight; does this Box-Cox result, which transforms only the response, justify that whole log-log move? Explain what Box-Cox does and does not tell you here.

EP 10.3 (what would change if). A student fits weighted least squares to the strongx physics data with weights $w_i = 1/\text{sd}_i^2$, then multiplies every weight by 10 and refits, expecting the line to “trust the precise points ten times more” and shift.

strongx <- read.csv("data/strongx.csv")
w <- 1 / strongx$sd^2
fit1 <- lm(crossx ~ energy, data = strongx, weights = w)
fit2 <- lm(crossx ~ energy, data = strongx, weights = 10 * w)
rbind(w = coef(fit1), tenw = coef(fit2))
c(se_energy_w = summary(fit1)$coef["energy", "Std. Error"],
  se_energy_10w = summary(fit2)$coef["energy", "Std. Error"])

     (Intercept)   energy
w       148.4732 530.8354
tenw    148.4732 530.8354
  se_energy_w se_energy_10w
        47.55         47.55

The slope, the intercept, and the slope’s standard error are identical to the last digit shown. Explain why multiplying every weight by the same constant changes nothing you would report, referring to the WLS model $\operatorname{Var}\{\varepsilon_i\} = \sigma^2/w_i$. What single quantity does change, and why is it not a result?

EP 10.4 (interpret output in context). The output below shows the fitted values from the plain linear regression of Galapagos plant species on the five geographic predictors, for the five islands with the smallest fitted values.

gala <- read.csv("data/gala.csv")
fit_lin <- lm(Species ~ Area + Elevation + Nearest + Scruz + Adjacent, data = gala)
pred <- predict(fit_lin)
head(data.frame(island = gala$island, actual = gala$Species,
                fitted = round(pred, 1))[order(pred), ], 5)

      island actual fitted
   Coamano      2  -36.4
    Darwin     10   -9.0
 Bartolome     31   -7.3
  Gardner1     58   -4.0
  Genovesa     40   -0.5

Interpret this output in the context of choosing a remedy. What is impossible about these predictions, why does it happen, and what does it tell you about whether any transformation of the response could rescue the linear normal model? Name the honest fix and the chapter that supplies it.

EP 10.5 (evaluate a claim). An analyst refits the savings model with Huber robust regression and tells a colleague: “The robust regression fit automatically downweights every unusual country, so I no longer need the Chapter 9 diagnostics.” The output lists the four highest-leverage countries with their Huber weights; for contrast, Zambia has a leverage value of 0.064 and Huber weight 0.472.

savings <- read.csv("data/savings.csv")
fit_ols <- lm(sr ~ pop15 + pop75 + dpi + ddpi, data = savings)
fit_hub <- MASS::rlm(sr ~ pop15 + pop75 + dpi + ddpi, data = savings)
h <- hatvalues(fit_ols); wt <- fit_hub$w
data.frame(country = savings$country, leverage_value = round(h, 3),
           huber_weight = round(wt, 3))[order(-h)[1:4], ]

       country leverage_value huber_weight
         Libya          0.531        1.000
 United States          0.334        1.000
         Japan          0.223        0.879
       Ireland          0.212        1.000

Evaluate the analyst’s claim using these numbers. Which country is the clearest counterexample, what general limitation of Huber regression does it expose, and what should the analyst still do?

Chapter game¶