Choosing a Transformation to Linearize Data

Part 1

Why straighten a curved scatter plot?

Ordinary least-squares regression assumes four things about the model \( y = \beta_0 + \beta_1 x + \varepsilon \): the mean of \(y\) is a straight-line function of \(x\), the errors \(\varepsilon\) have constant variance, they are independent, and for intervals and tests they are roughly normal. A curved scatter plot breaks the first assumption. A plot that fans out breaks the second. A transformation re-expresses one or both variables so that both assumptions hold on the new scale.

Straightening pays off in three ways. The model fits the data, so residuals carry no leftover pattern. The slope becomes a single, interpretable number such as a growth rate or an exponent. And confidence and prediction intervals have the right width everywhere, not too narrow at one end and too wide at the other.

You do not always need to linearize. If you only need predictions and have plenty of data, a nonlinear model or a flexible method may serve you better. If physics hands you a nonlinear law with parameters you care about, fit that law directly with nonlinear regression. Part 13 covers these exits.

Part 9

Zeros, negatives, proportions and offsets

Zeros

\(\ln 0\) and \(1/0\) do not exist. The usual patch is \(\ln(y + c)\), but the answer depends on \(c\). Adding 1 to counts in the hundreds changes almost nothing; adding 1 to counts near 0 changes the shape completely. If you must shift, choose \(c\) for a reason (half the detection limit, or the smallest non-zero value divided by two), report it, and check that the conclusions survive a different \(c\). Often the better answer is \(\sqrt{y}\), which handles zeros, or a Poisson or Tweedie model with a log link, which models \(\ln E(y)\) without ever logging a zero.

Negative values

The power ladder needs positive data. For variables that cross zero (deviations, residuals, temperature in °C), use the Yeo–Johnson transformation, which applies Box–Cox to \(y + 1\) for \(y \ge 0\) and a mirrored version for \(y < 0\):

Better still, ask whether the variable has a natural zero. Temperature converted to kelvin is positive, and Arrhenius needs kelvin anyway.

Proportions and percentages

A proportion is trapped between 0 and 1, so its relationship with x flattens at both ends into an S-curve. The logit \(\ln[p/(1-p)]\) stretches both ends to infinity and usually straightens it. Exact 0% and 100% readings are undefined on the logit scale; with counts behind them, use logistic regression on the counts instead of transforming percentages.

Offsets and start values

Many curves are power laws in \(x - x_0\), not in \(x\): a threshold stress, an ambient temperature, a minimum detectable level. When the data span a narrow ratio, no rung bends them, but subtracting the right \(x_0\) restores a wide ratio and the ladder works again. Load Recovery time vs set-point in the lab, try \(-1/x\) with no offset, then set the offset to 95.

Part 10

Check that the transformation worked

A transformation is a hypothesis about the data. Test it the same way you test any regression model:

1. Residuals vs fitted. A horizontal band with no curve and no funnel. This is the single most important plot.
2. Residuals vs each predictor. Catches curvature in one x that the fitted-value plot hides in a multiple regression.
3. Normal probability plot of residuals. Roughly straight. Heavy tails after transforming may mean outliers rather than the wrong power.
4. Residuals in time order. Transformations do nothing for autocorrelation or drift.
5. Lack-of-fit test. With replicate x values, Minitab's pure-error lack-of-fit test checks linearity formally.

Comparing models fairly

\(R^{2}\) measures the share of variance explained in the response you fitted. If one model predicts \(y\) and another predicts \(\ln y\), their \(R^{2}\) values describe different quantities and cannot be ranked. Three fair comparisons:

• Back-transform each model's predictions to the original scale and compute \(R^{2}_{\text{orig}} = 1 - \sum(y_i - \tilde{y}_i)^{2} / \sum(y_i - \bar{y})^{2}\), where \(\tilde{y}_i\) is the back-transformed prediction. The lab reports this.
• Compare the Box–Cox log-likelihood \(\ell(\lambda)\) at each candidate power. The geometric-mean scaling makes them comparable.
• Compare out-of-sample prediction error on the original scale with cross-validation.

Comparing models with different x transforms but the same y is simpler: the \(R^{2}\) values and \(S\) are directly comparable.

Part 11

Worked example: defect growth after a roll change

A finishing mill records surface defects per 100 pipes at six points after a roll change. Tonnage x is in thousands of tonnes.

Step 1: the raw fit looks good and is wrong. A straight line on raw y gives \(R^{2} = 0.943\). The residuals run plus, minus, minus, minus, minus, plus: a U-shape. The curve bulges down and right, and the jumps between points grow with y, which suggests a multiplicative process. Move y down.

Step 2: fit \(\ln y\) on x.

The residuals on the log scale are all within \(\pm 0.015\) with no pattern. Strictly, \(0.9998\) cannot be compared with \(0.943\) because the responses differ; the evidence is the residual pattern.

Step 3: interpret the slope. Each extra thousand tonnes multiplies the defect rate by \(e^{0.3415}\):

Step 4: back-transform a prediction. At 7 kt the model gives \(\widehat{\ln y} = -0.5430 + 0.3415 \times 7 = 1.8478\), so the predicted median rate is \(e^{1.8478} = 6.35\) defects per 100 pipes. The smearing correction here is \(\tfrac{1}{n}\sum e^{e_i} = 1.00003\) because the residuals are tiny, so the mean and median agree. With noisier data they would not; Part 12 shows how far apart they can drift.

Step 5: act. The rate doubles about every 2 kt. Seven kilotonnes is outside the observed range, so treat the 7 kt figure as a warning for planning the next roll change, not a forecast to bank on.

Part 17

Python and R equivalents

The same workflow in code: profile Box–Cox on the regression, run one Box–Tidwell step, and back-transform with smearing.

import numpy as np
import statsmodels.api as sm
from scipy import stats

def boxcox_profile(x, y, lambdas=np.arange(-2, 3.01, 0.01)):
    """Regression-aware Box-Cox: profile SSE of z(lambda) on x."""
    gm = np.exp(np.mean(np.log(y)))
    X = sm.add_constant(x)
    ll = []
    for lam in lambdas:
        z = gm * np.log(y) if abs(lam) < 1e-9 else (y**lam - 1) / (lam * gm**(lam - 1))
        sse = sm.OLS(z, X).fit().ssr
        ll.append(-len(y) / 2 * np.log(sse / len(y)))
    ll = np.array(ll)
    best = lambdas[ll.argmax()]
    ci = lambdas[ll >= ll.max() - stats.chi2.ppf(0.95, 1) / 2]
    return best, (ci.min(), ci.max())

def box_tidwell_step(x, y, alpha=1.0):
    """One Box-Tidwell update for the power on x (x > 0)."""
    z = x**alpha
    b1 = sm.OLS(y, sm.add_constant(z)).fit().params[1]
    gamma = sm.OLS(y, sm.add_constant(np.column_stack([z, z * np.log(z)]))).fit().params[2]
    return alpha * (1 + gamma / b1)

# Fit ln(y) on x, then back-transform with Duan's smearing estimator
fit = sm.OLS(np.log(y), sm.add_constant(x)).fit()
smear = np.mean(np.exp(fit.resid))
mean_pred = np.exp(fit.fittedvalues) * smear      # mean of y
median_pred = np.exp(fit.fittedvalues)            # median of y

# Counts: model the mean directly instead of transforming
poisson = sm.GLM(counts, sm.add_constant(length), family=sm.families.Poisson()).fit()

library(MASS)   # boxcox()
library(car)    # boxTidwell()
fit <- lm(y ~ x)
boxcox(fit, lambda = seq(-2, 2, 0.05))   # regression-aware profile with 95% CI
boxTidwell(y ~ x)                        # estimated power on x with a score test

scipy.stats.boxcox(y) and scipy.stats.yeojohnson(y) are univariate, like Minitab's control-chart Box–Cox. They are right for normalizing a single variable and wrong for picking a regression transform.

References

• Box, G. E. P., & Cox, D. R. (1964). An analysis of transformations. Journal of the Royal Statistical Society, Series B, 26(2), 211–252.
• Box, G. E. P., & Tidwell, P. W. (1962). Transformation of the independent variables. Technometrics, 4(4), 531–550.
• Duan, N. (1983). Smearing estimate: A nonparametric retransformation method. Journal of the American Statistical Association, 78(383), 605–610.
• Mosteller, F., & Tukey, J. W. (1977). Data Analysis and Regression. Addison-Wesley.
• Tukey, J. W. (1977). Exploratory Data Analysis. Addison-Wesley.
• Yeo, I.-K., & Johnson, R. A. (2000). A new family of power transformations to improve normality or symmetry. Biometrika, 87(4), 954–959.
• Montgomery, D. C., Peck, E. A., & Vining, G. G. (2021). Introduction to Linear Regression Analysis (6th ed.). Wiley.