53
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ I don't know of any facilities for doing this in R and, 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$ – Macro Sep 5 '12 at 12:18
  • 1
    $\begingroup$ I am under the impression I've already seen a paper talking about that. Maybe googling with "covariate" instead of "independent variable" is more judicious. $\endgroup$ – Stéphane Laurent Sep 5 '12 at 14:19
  • $\begingroup$ I think (from memory ...) there is some implementation of this in the car package (R). But youd also investigate gam's like in package gmcv. $\endgroup$ – kjetil b halvorsen Sep 10 '12 at 23:22
  • $\begingroup$ A thread that discusses the automatic inference of Box-Cox parameters (by simultaneously transforming all independent variables along with the dependent variable) has appeared at stats.stackexchange.com/questions/60431/…. $\endgroup$ – whuber May 30 '13 at 17:06
72
+50
$\begingroup$

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.

Plots

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.

$\endgroup$
  • 2
    $\begingroup$ Hello dear Whuber. What an interesting reply, I read it with pleasure, thank you! (and I will also play with it a bit to see how it might fit the issue I am working with) $\endgroup$ – Tal Galili Sep 6 '12 at 4:03
  • 3
    $\begingroup$ @Andre To use the last part of this code you need to bind the x and y variables into an $n$ by $2$ matrix, as in 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$ – whuber Feb 25 '14 at 18:09
  • 2
    $\begingroup$ @landroni Integral and small fractional powers often arise in physical, chemical, and biological theories, as well as through geometric considerations. (For example, when a variable is a volume, its cube root is a length--which is interpretable--whereas, say, its seventh root has no simple geometric interpretation.) Other powers rarely have any such interpretation. $\endgroup$ – whuber Mar 6 '14 at 23:15
  • 3
    $\begingroup$ @Frank That's right; it is explicitly and unabashedly an exploratory technique. Note that it does not even claim to be predictive. Exploration can only suggest ways to proceed. It is possible to imagine allocating four df of your modeling budget to estimating these transformations, though--and the estimation could be incorporated automatically into the fitting algorithm using Tukey's approach or otherwise (ML is an obvious possibility). $\endgroup$ – whuber Mar 24 '14 at 21:08
  • 5
    $\begingroup$ They don't take into account the (very significant amount of) uncertainty in the $Y$-transform unless the centering of $Y$ and estimation of $\lambda$ are explicit parts of the information matrix. $\endgroup$ – Frank Harrell Sep 26 '14 at 18:01
11
$\begingroup$

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.

$\endgroup$
9
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.