61
$\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$
4
  • 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, 2012 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$ Sep 5, 2012 at 14:19
  • $\begingroup$ I think (from memory ...) there is some implementation of this in the car package (R). But you'd also investigate gam's like in package mgcv. $\endgroup$ Sep 10, 2012 at 23:22
  • 1
    $\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, 2013 at 17:06

4 Answers 4

87
+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$
16
  • 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, 2014 at 18:09
  • 3
    $\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, 2014 at 23:15
  • 3
    $\begingroup$ I don't think Tukey's approach will come close to preserving confidence interval coverage. It is meant for a more predictive/curve fitting model rather than for use with formal statistical inference. $\endgroup$ Mar 24, 2014 at 20:57
  • 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, 2014 at 21:08
  • 6
    $\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$ Sep 26, 2014 at 18:01
12
$\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$
10
$\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$
0
$\begingroup$

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.

$\endgroup$
1
  • 1
    $\begingroup$ There is an R package on CRAN, mfp. There are also functions in some other packages. $\endgroup$ Nov 15, 2022 at 23:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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