Is there a Box-Cox like transformation for independent variables? That is, a transformation that optimizes the $x$ variable so that the y~f(x) will make a more reasonable fit for a linear model?
If so, is there a function to perform this with R?
John Tukey advocated his "three point method" for finding re-expressions of variables to linearize relationships.
I will illustrate with an exercise from his book, Exploratory Data Analysis. These are mercury vapor pressure data from an experiment in which temperature was varied and vapor pressure was measured.
pressure <- c(0.0004, 0.0013, 0.006, 0.03, 0.09, 0.28, 0.8, 1.85, 4.4,
9.2, 18.3, 33.7, 59, 98, 156, 246, 371, 548, 790) # mm Hg
temperature <- seq(0, 360, 20) # Degrees C
The relation is strongly nonlinear: see the left panel in the illustration.
Because this is an exploratory exercise, we expect it to be interactive. The analyst is asked to begin by identifying three "typical" points in the plot: one near each end and one in the middle. I have done so here and marked them in red. (When I first did this exercise long ago, I used a different set of points but arrived at the same results.)
In the three point method, one searches--by brute force or otherwise--for a Box-Cox transformation that when applied to one of the coordinates--either y or x--will (a) place the typical points approximately on a line and (b) uses a "nice" power, usually chosen from a "ladder" of powers that might be interpretable by the analyst.
For reasons that will become apparent later, I have extended the Box-Cox family by allowing an "offset" so that the transformations are in the form
$$x \to \frac{(x + \alpha)^\lambda - 1}{\lambda}.$$
Here's a quick and dirty R implementation. It first finds an optimal $(\lambda,\alpha)$ solution, then rounds $\lambda$ to the nearest value on the ladder and, subject to that restriction, optimizes $\alpha$ (within reasonable limits). It's incredibly quick because all the calculations are based on just those three typical points out of the original dataset. (You could do them with pencil and paper, even, which is exactly what Tukey did.)
box.cox <- function(x, parms=c(1,0)) {
lambda <- parms[1]
offset <- parms[2]
if (lambda==0) log(x+offset) else ((x+offset)^lambda - 1)/lambda
}
threepoint <- function(x, y, ladder=c(1, 1/2, 1/3, 0, -1/2, -1)) {
# x and y are length-three samples from a dataset.
dx <- diff(x)
f <- function(parms) (diff(diff(box.cox(y, parms)) / dx))^2
fit <- nlm(f, c(1,0))
parms <- fit$estimate #$
lambda <- ladder[which.min(abs(parms[1] - ladder))]
if (lambda==0) offset = 0 else {
do <- diff(range(y))
offset <- optimize(function(x) f(c(lambda, x)),
c(max(-min(x), parms[2]-do), parms[2]+do))$minimum
}
c(lambda, offset)
}
When the three-point method is applied to the pressure (y) values in the mercury vapor dataset, we obtain the middle panel of the plots.
data <- cbind(temperature, pressure)
n <- dim(data)[1]
i3 <- c(2, floor((n+1)/2), n-1)
parms <- threepoint(temperature[i3], pressure[i3])
y <- box.cox(pressure, parms)
In this case, parms turns out to equal $(0,0)$: the method elects to log-transform the pressure.
We have reached a point analogous to the context of the question: for whatever reason (usually to stabilize residual variance), we have re-expressed the dependent variable, but we find that the relation with an independent variable is nonlinear. So now we turn to re-expressing the independent variable in an effort to linearize the relation. This is done in the same way, merely reversing the roles of x and y:
parms <- threepoint(y[i3], temperature[i3])
x <- box.cox(temperature, parms)
The values of parms for the independent variable (temperature) are found to be $(-1, 253.75)$: in other words, we should express the temperature as degrees Celsius above $-254$C and use its reciprocal (the $-1$ power). (For technical reasons, the Box-Cox transformation further adds $1$ to the result.) The resulting relation is shown in the right panel.
By now, anybody with the least science background has recognized that the data are "telling" us to use absolute temperatures--where the offset is $273$ instead of $254$--because those will be physically meaningful. (When the last plot is re-drawn using an offset of $273$ instead of $254$, there is little visible change. A physicist would then label the x-axis with $1/(1-x)$: that is, reciprocal absolute temperature.)
This is a nice example of how statistical exploration needs to interact with understanding of the subject of investigation. In fact, reciprocal absolute temperatures show up all the time in physical laws. Consequently, using simple EDA methods alone to explore this century-old, simple, dataset, we have rediscovered the Clausius-Clapeyron relation: the logarithm of the vapor pressure is a linear function of the reciprocal absolute temperature. Not only that, we have a not very bad estimate of absolute zero ($-254$ degrees C), from the slope of the righthand plot we can calculate the specific enthalpy of vaporization, and--as it turns out--a careful analysis of the residuals identifies an outlier (the value at a temperature of $0$ degrees C), shows us how the enthalphy of vaporization varies (very slightly) with temperature (thereby violating the Ideal Gas Law), and ultimately can give us accurate information about the effective radius of the mercury gas molecules! All that from 19 data points and some basic skills in EDA.
data <- cbind(temperature, pressure). (This was done in order to match the way in which (x,y) data are usually stored in R, either as a matrix or data frame.) I apologize for omitting that step and have edited the answer to include it--see whether that works.
$\endgroup$
Take a look at these slides on "Regression diagnostics" by John Fox (available from here, complete with references), which briefly discuss the issue of transforming nonlinearity. It covers Tukey's "bulging rule" for selecting power transformations (addressed by the accepted answer), but also mentions the Box-Cox and Yeo-Johnson families of transformations. See Section 3.6 of the slides. For a more formal take by the same author see J. Fox, Applied Regression Analysis and Generalized Linear Models, Second Edition (Sage, 2008).
As for actual R packages that help with this, absolutely take a look at the car package, authored by J. Fox and S. Weisberg. This package accompanies J. Fox and S. Weisberg, An R Companion to Applied Regression, Second Edition, (Sage, 2011), another must-read. Using that package you can start off from basicPower() (simple power transformations), bcPower() (Box-Cox transformations) and yjPower() (Yeo-Johnson transformations). There is also powerTransform():
The function powerTransform is used to estimate normalizing transformations of a univariate or a multivariate random variable.
Check both books for more details on the theory behind these transformations and on computational approaches.
There are many advantages to making estimation of covariate transformations a formal part of the estimation process. This will recognize the number of parameters involved and produced good confidence interval coverage and type I error preservation. Regression splines are some of the best approaches. And splines will work with zero and negative values of $X$ unlike logarithmic approaches.
The method of fractional polynomials due to Royston and Altman (1994) https://doi.org/10.2307/2986270 (paywalled) may be just such a method you might use to handle the case of an optimal power law relating the X to the Y. In this case, you believe that the $Y$ (or an appropriate transformation) is truly normally distributed about the conditional mean response given an $X$ (one or more predictors), but the exact functional form of the mean is unknown:
$$ Y = p(X) + \epsilon $$
with $p(X)$ an unknown mapping of a linear combination of $X$, and $\epsilon \sim \mathcal{N}(0, \sigma^2$). Note in this particular case it does not suffice to transform the $Y$ variable according to Box-Cox, because the error term loses normality, and the error becomes correlated with the response.
In the spectral theory of analysis, polynomials give us good local approximations to smooth functions, that's why we are often interested (as in the case of splines, or in the Box Cox transformation) expressing the functional relationship according to a power law.
The advantage of fractional polynomials over splines is that the procedure gives one an estimate of the optimal power law that may be used generally and with relative simplicity. I consider this an optimal process when the point of the analysis is to communicate a simple approximation, as would be the case in statistical mechanics when we might estimate simple functional relationships (like stopping distance relates to the square root of velocity at pre-braking speed). Splines, however, are far more flexible, and the extra degrees of freedom are well spent in terms of predictiveness to allow for breakpoints, and higher order terms where necessary.
Fractional polynominals are an interesting and underutilized procedure. Stata implemented fracpoly but to my knowledge nobody has done this in R.
mfp. There are also functions in some other packages.
$\endgroup$
Commented
Nov 15, 2022 at 23:27
Rand, thinking about it for a moment, I'm not sure exactly how one would do this at all. What criteria would you optimize to ensure the "most linear" transformation? $R^2$ is tempting but, as seen in my answer here, $R^2$ alone cannot be used to see whether the linearity assumption of a model is satisfied. Did you have some criteria in mind? $\endgroup$mgcv. $\endgroup$