Finding Optimal Normalizing Transformations via bestNormalize

The bestNormalize R package was designed to help users find a transformation that can effectively normalize a vector regardless of its actual distribution. Each of the many normalization techniques that have been developed has its own strengths and weaknesses, and deciding which to use until data are fully observed is difficult or impossible. This package facilitates choosing between a range of possible transformations and will automatically return the best one, i.e., the one that makes data look the most normal. To evaluate and compare the normalization efficacy across a suite of possible transformations, we developed a statistic based on a goodness of fit test divided by its degrees of freedom. Transformations can be seamlessly trained and applied to newly observed data and can be implemented in conjunction with caret and recipes for data preprocessing in machine learning workflows. Custom transformations and normalization statistics are supported.


Introduction
The bestNormalize package contains a suite of transformation-estimating functions that can be used to normalize data. The function of the same name attempts to find and execute the best of all of these potential normalizing transformations. In this package, we define "normalize" as in "to render data Gaussian", rather than transform them to a specific scale.
There are many instances where researchers may want to normalize a variable. First, there is the (often problematic) assumption of normality of the outcome (conditional on the covariates) in the classical linear regression problem. Over the years, many methods have been used to relax this assumption: generalized linear models, quantile regression, survival models, etc. One technique that is still somewhat popular in this context is to "beat the data" to look normal via some kind of normalizing transformation. This could be something as simple as a log transformation or something as complex as a Yeo-Johnson transformation (Yeo and Johnson, 2000). In fact, many complex normalization methods were designed expressly to find a transformation that could render regression residuals Gaussian. While perhaps not the most elegant solution to the problem, often, this technique works well as a quick solution. Another increasingly popular application of normalization occurs in applied regression settings with highly skewed distributions of the covariates (Kuhn and Johnson, 2013). In these settings, there exists the tendency to have high leverage points (and highly influential points), even when one centers and scales the covariates. When examining interactions, these influential points can become especially problematic since the leverage of that point gets amplified for every interaction in which it is involved. Normalization of such covariates can mitigate their leverage and influence, thereby allowing for easier model selection and more robust downstream predictor manipulations (such as principal components analysis), which can otherwise be sensitive to skew or outliers. As a result, popular model selection packages such as caret (Kuhn, 2017) and recipes (Kuhn and Wickham, 2018) have built-in mechanisms to normalize the predictor variables (they call this "preprocessing"). This concept is unique in that it forgoes the assumption of linearity between the outcome (Y) and the covariate, opting instead for a linear relationship between Y and the transformed value of the covariate (which in many cases may be more plausible).
This package is designed to make normalization effortless and consistent. We have also introduced Ordered Quantile (ORQ) normalization via the orderNorm function, which uses a rank mapping of the observed data to the normal distribution in order to guarantee normally distributed transformed data (if ties are not present). We have shown how ORQ normalization performs very consistently across different distributions, successfully normalizing left-or right-skewed data, multi-modal data, and even data generated from a Cauchy distribution (Peterson and Cavanaugh, 2019).
In this paper, we describe our R package bestNormalize, which is available via the Comprehensive R Archive Network (CRAN). First, we describe normalization methods that have been developed and that we implement in the package. Second, we describe the novel cross-validation-based estimation procedure, which we utilize to judge the normalization efficacy of our suite of normalization transformations. Third, we go through some basic examples of bestNormalize functionality and a simple implementation of our methods within the recipes package. We illustrate a more in-depth usecase in a car pricing application, performing a transform-both-sides regression as well as comparing the performance of several predictive models fit via caret. Finally, we conclude by discussing the pros and cons of normalization in general and future directions for the package.

Normalization methods
Many normalization transformation functions exist, and though some can be implemented well in existing R packages, bestNormalize puts them all under the same umbrella syntax. This section describes each transformation contained in the bestNormalize suite.

The Box-Cox transformation
The Box-Cox transformation was famously proposed in Box and Cox (1964) and can be implemented with differing syntax and methods in many existing packages in R (e.g., caret, MASS (Venables and Ripley, 2002), and more). It is a straightforward transformation that typically only involves one parameter, λ: where x refers to the datum in its original unit (pre-transformation). Given multiple observations, the λ parameter can be estimated via maximum likelihood, and x must be greater than zero.

The Yeo-Johnson transformation
The Yeo-Johnson transformation (Yeo and Johnson, 2000) attempts to find the value of λ in the following equation that minimizes the Kullback-Leibler distance between the normal distribution and the transformed distribution.
This method has the advantage of working without having to worry about the domain of x. As with the Box-Cox λ, this λ parameter can be estimated via maximum likelihood.

The Lambert W x F transformation
The Lambert W x F transformation, proposed in Goerg (2011) and implemented in the LambertW package, is essentially a mechanism that de-skews a random variable X using moments. The method is motivated by a system theory and is alleged to be able to transform any random variable into any other kind of random variable, thus being applicable to a large number of contexts. One of the package's main functions is Gaussianize, which is similar in spirit to the purpose of this package. However, this method may not perform as well on certain shapes of distributions as other candidate transformations; see Peterson and Cavanaugh (2019) for some examples.
The Gaussianize transformation can handle three types of transformations: skewed, heavytailed, and skewed heavy-tailed. For more details on this transformation, consult the LambertW documentation. 1 While the transformations contained and implemented by bestNormalize are reversible (i.e., 1-1), in rare circumstances, we have observed that the lambert function can yield non-reversible transformations.

The Ordered Quantile technique
The ORQ normalization technique (orderNorm) is based on the following transformation (originally discussed, as far as we can find, in Bartlett (1947) and further developed in Van der Waerden (1952)): Let x refer to the original data. Then the transformation is: This nonparametric transformation as defined works well on the observed data, but it is not trivial to implement in modern settings where the transformation needs to be applied on new data; we discussed this issue and our solution to it in Peterson and Cavanaugh (2019). Basically, on new data within the range of the original data, ORQ normalization will linearly interpolate between two of the original data points. On new data outside the range of the original data, the transformation extrapolates using a shifted logit approximation of the ranks to the original data. This is visualized below via the iris data set on the Petal.Width variable.

ORQ Normalization
Original Value

Transformed Value
Original data Transformed values for new data Approximation for extrapolation The shifted logit extrapolation ensures that the function is 1-1 and can handle data outside the original (observed) domain. The effects of the approximation will usually be relatively minimal since we should not expect to see many observations outside the observed range if the training set sample size is large relative to the test set. The ORQ technique will not guarantee a normal distribution in the presence of ties, but it still could yield the best normalizing transformation when compared to the other possible approaches. More information on ORQ normalization can be found in Peterson and Cavanaugh (2019) or in the bestNormalize documentation.

Other included transformations
In addition to the techniques above, the bestNormalize package performs and evaluates:

Other not-included transformations
A range of other normalization techniques has been proposed that are not included in this package (at the time of writing). These include (but are not limited to): Modified Box-Cox (Box and Cox, 1964), Manly's Exponential (Manly, 1976), John/Draper's Modulus (John and Draper, 1980), and Bickel/Doksum's Modified Box-Cox (Bickel and Doksum, 1981). However, it is straightforward to add new transformations into the same framework as other included transformations; each one is treated as its own S3 class, so in order to add other transformations, all one must do is define a new S3 class and provide the requisite S3 methods. To this end, we encourage readers to submit a pull request to the package's GitHub page with new transformation techniques that could be then added as a default in bestNormalize. Otherwise, in a later section, we show how users can implement custom transformations alongside the default ones described above.
Which transformation "best normalizes" the data?
The bestNormalize function selects the best transformation according to an extra-sample estimate of the Pearson P statistic divided by its degrees of freedom (DF ). This P statistic is defined as where O i is the number observed, and E i is the number of expected (under the hypothesis of normality) to fall into "bin" i. The bins (or "classes") are built such that observations will fall into each one with equal probability under the hypothesis of normality. A variety of alternative normality tests exist, but this particular one is relatively interpretable as a goodness of fit test, and the ratio P /DF can be compared between transformations as an absolute measure of departure from normality. Specifically, if the data in question follow a normal distribution, this ratio will be close to 1 or lower. The transformation which produces data with the lowest normality statistic is thus the most effective at normalizing the data, and gets selected by bestNormalize. The bestNormalize package utilizes nortest (Gross and Ligges, 2015) to compute this statistic; more information on its computation and degrees of freedom can be found in D'Agostino (1986) and Thode (2002).
Normality statistics for all candidate transformations can be estimated and compared with one simple call to bestNormalize, whose output makes it easy to see which transformations are viable and which are not. We have found that while complicated transformations are often most effective and therefore selected automatically, sometimes a simple transformation (e.g., the log or identity transforms) may be almost as effective, and ultimately the latter type will yield more interpretable results.
It is worth noting that when the normality statistic is estimated on in-sample data, the ORQ technique is predestined to be most effective since it is forcing its transformed data to follow a normal distribution exactly (Peterson and Cavanaugh, 2019). For this reason, by default, the bestNormalize function calculates an out-of-sample estimate for the P /DF statistic. Since this method necessitates cross-validation, it can be computationally frustrating for three reasons: (1) the results and the chosen transformation can depend on the seed, (2) it takes considerably longer to estimate than the in-sample statistic, and (3) it is unclear how to choose the number of folds and repeats.
In order to mediate these issues, we have built several features into bestNormalize. Issue (1) is only important for small sample sizes, and when it is a concern, the best transformations should look similar to one another. We address two solutions to (2) in the next section. In short, we have methods to parallelize or simplify the estimation of the statistic. For (3), we recommend 10-fold cross-validation with 5 repeats as the default, but if the sample is small, we suggest using 5 (or fewer) folds instead with more repeats; accurate estimation of P /DF requires a relatively large fold size (as a rule of thumb, 20 observations per fold seems to be enough for most cases, but this unfortunately depends on the distribution of the observed data).

Simple examples
In this section, we illustrate a simple use-case of the functions provided in bestNormalize.

Basic implementation
First, we will generate and plot some skewed data: x <-rgamma(250, 1, 1) Finally, we can execute the best performing normalization on new data with predict(BNobject, new_x) or reverse the transformation with predict(BNobject, new_x_t,inverse = TRUE). Note that normalized values can either be obtained using predict or by extracting x.t from the object. The best transformation, in this case, is plotted in Figure 4. Normalizing transformations

In-sample normalization efficacy
To examine how each of the normalization methods performed (in-sample), we can visualize the transformed values in histograms (Figure 6), which plot the transformed data, x.t, stored in the transformation objects we created previously.  Evidently, ORQ normalization appears to have worked perfectly to normalize the data (as expected), and the Box-Cox method seemed to do quite well too.

Out-of-sample normalization efficacy
The bestNormalize function performs repeated (r=5) 10-fold cross-validation (CV) by default and stores the estimated normality statistic for each left-out fold/repeat into oos_preds. Users can access and visualize these results via a boxplot (see below), which may give some insight into whether the transformation is truly preferred by the normality statistic or if another (possibly simpler) transformation can be applied that would achieve the approximately the same results. In this example, Box-Cox, square-root, Yeo-Johnson, and ORQ seem to do similarly well, whereas the identity transform 2 , hyperbolic arc-sine, logging, and exponentiation are performing worse. Leave-one-out CV can be optionally performed in bestNormalize via the loo argument, which, if set to TRUE, will compute the leave-one-out CV transformations for each observation and method. Specifically, bestNormalize will be run n separate times where each observation is individually left out of the fitting process and subsequently plugged back in to get a "leave-one-out transformed value". Instead of taking the mean across repeats and folds, in this case, we estimate normalization efficacy using the full distribution of leave-one-out transformed values. This option is computationally intensive. Note that as with the "in-sample" normality statistics, the leave-one-out CV approach tends to select the ORQ transformation since ORQ's performance improves as the number of points in the training set relative to the testing set increases.

Improving speed of estimation
Because bestNormalize uses repeated CV by default to estimate the out-of-sample normalization efficacy, it can be quite slow for larger objects. There are several means of speeding up the process. Each comes with some pros and cons. The first option is to specify out_of_sample = FALSE. This will highly speed up the process. However, for reasons previously discussed, ORQ normalization will always be chosen unless allow_orderNorm = FALSE. Therefore, a user might as well use the orderNorm function directly as opposed to only setting out_of_sample = FALSE since the end result will be the same (and orderNorm will run much faster). Note below that the in-sample normality results may differ slightly from the leave-one-out even when this may be unexpected (i.e., for the log transformation); this is due to slight differences in the standardization statistics. Another option to improve estimation efficiency is to use the built-in parallelization functionality. The repeated CV process can be parallelized via the cluster argument and the parallel and doRNG (Gaujoux, 2020) packages. A cluster can be set up with makeCluster and passed to bestNormalize via the cluster = argument.
cl <-parallel::makeCluster(5) b <-bestNormalize(x, cluster = cl, r = 10, quiet = TRUE) parallel::stopCluster(cl) The amount by which this parallelization will speed up the estimation of out-of-sample estimates depends (for the most part) on the number of repeats, the number of cores, and the sample size of the vector to be normalized. The plot below shows the estimation time for a run of bestNormalize with 15 repeats of 10-fold CV on a gamma-distributed random variable with various sample sizes and numbers of cores.

Implementation with caret, recipes
The step_best_normalize and the step_orderNorm functions can be utilized in conjunction with the recipes package to preprocess data in machine learning workflows with tidymodels (Kuhn and Wickham, 2020) or in combination with caret. The basic usage within recipes is shown below; for implementation with caret, refer to this paper's application.

Additional customization
Two important means of customization are available: 1) users may add custom transformation functions to be assessed alongside the default suite of normalization methods, and 2) users may change the statistic used "under the hood" by bestNormalize to estimate the departure from normality of the transformed data. This section contains examples and guidance for both extensions.

1) Adding user-defined functions
Via the new_transforms argument, users can use bestNormalize's machinery to compare custom, user-defined transformation functions to those included in the package. Below, I consider an example where a user may wish to compare the cube-root function with those provided in the package. bestNormalize requires two functions to implement this: the transformation function and an associated predict method. The custom cube-root transformation shown below is simple, but its skeleton can readily be made arbitrarily more complex. Evidently, the cube-root was the best normalizing transformation for this gamma-distributed random variable, performing comparably to the Box-Cox transformation.

2) Re-defining normality
The question "what is normal?" outside of a statistical discussion is quite loaded and subjective. Even in statistical discussions, many authors have contributed to the question of how to best detect departures from normality; these solutions are diverse, and several have been implemented well in nortest already. In order to accommodate those with varying opinions on the best definition of normality, we have included a feature that allows users to specify a custom definition of a normality statistic. This customization can be accomplished via the norm_stat_fn argument, which takes a function that will then be applied in lieu of the Pearson test statistic divided by its degree of freedom to assess normality.
The user-defined function must take an argument x, which indicates the data on which a user wants to evaluate the statistic.
Here is an example using the Lilliefors (Kolmogorov-Smirnov) normality test statistic: Here is an example using the Lillifors (Kolmogorov-Smirnov) normality test's p-value: (dont_do_this <-bestNormalize(x, norm_stat_fn = function(x) nortest::lillie.test(x)$p)) #> Best Normalizing transformation with 100 Observations #> Estimated Normality Statistics (using custom normalization statistic) #> -arcsinh ( Note: bestNormalize will attempt to minimize this statistic by default, which is definitely not what you want to do when calculating the p-value. This is seen in the example above, where the worst normalization transformation, exponentiation, is chosen. In this case, a user is advised to either manually select the best one or reverse their defined normalization statistic (in this case by subtracting it from 1): best_transform <-names(which.max(dont_do_this$norm_stats)) do_this <-dont_do_this$other_transforms[[best_transform]] or_this <-bestNormalize(x, norm_stat_fn = function(x) 1-nortest::lillie.test(x)$p) A p-value for normality should not be routinely used as the sole selector of a normalizing transformation. A normality test's p-value, as a measure of the departure from normality, is confounded by the sample size (a high sample size may yield strong evidence of a practically insignificant departure from normality). Therefore, we suggest the statistic used should estimate the departure from normality rather the strength of evidence against normality (e.g., Royston, 1991).

Background
The autotrader data set was scraped from the autotrader website as part of this package (and because at the time of data collection in 2017, the package author needed to purchase a car). We apply the bestNormalize functionality to de-skew mileage, age, and price in a pricing model. See ?autotrader for more information on this data set. data("autotrader") autotrader$yearsold <-2017 -autotrader$Year

Transform-both-sides regression
Transform-both-sides (TBS) regression has several benefits that have been explored thoroughly elsewhere (see Harrell (2015) for an overview). Importantly, TBS regression can often (though not always) yield models that better satisfy assumptions of linear regression and mitigate the influence of outliers/skew. This approach has been shown to be useful in shrinking the size of prediction intervals while maintaining closer to nominal coverage in this data set (Peterson and Cavanaugh, 2019).
First, we will normalize the outcome (price).
(priceBN <-bestNormalize(autotrader$price)) #> We can see that the estimated normality statistic for the ORQ transformation is close to 1, so we know it is performing quite well despite the ties in the data. It is also performing considerably better than all of the other transformations. For age, we see something peculiar; none of the normalizing transformations performed well according to the normality statistics. By plotting the data, it becomes evident that the frequency of ties in age makes it very difficult to find a normalizing transformation (see figure below). Even so, orderNorm is chosen as it has the lowest estimated P /DF statistic. Next, we will fit a linear model on the transformed values of each variable for our TBS regression. The reverse-transformation functions will allow us to visualize how these variables affect model predictions in terms of their original units. -priceBN$x.t; m.t <-mileageBN$x.t; yo.t <-yearsoldBN$x.t fit <-lm(p.t~m.t + yo.t)  Unsurprisingly, we find that there are very significant relationships between transformed car price, mileage, and age. However, to interpret these values, we must resort to visualizations since there is no inherent meaning of a "one-unit increase" in the ORQ normalized measurements. We utilize the visreg package (Breheny and Burchett, 2017) to perform our visualizations, using predict.bestNormalize in conjunction with visreg's trans and xtrans options to view the relationship in terms of the original unit for the response and covariate respectively (formatting omitted). 3 For the sake of illustration, we have also plotted the estimated effect of a generalized additive (spline) model fit with mgcv (Wood, 2011). Below, we visualize the age effect, demonstrating how one might visualize the effect outside of visreg (plot formatting is omitted).

Implementation with recipes
To build a predictive model for the price variable that uses each vehicle's model and make in addition to its mileage and age, we can utilize the caret and recipes functionality to do so. This section outlines how to use bestNormalize in conjunction with these other popular ML packages. Price is logged instead of ORQ transformed in order to facilitate the interpretation of measures for prediction accuracy.
# Out of sample prediction accuracy results <-lapply(fits, function(x) { p <-c(predict(x, newdata = bake(rec, df_test))) yardstick::metrics(data.frame(est = exp(p), truth = df_test$price), truth = truth, estimate = est) }) results After normalization of mileage and age, a random forest had the optimal predictive performance on car price given a car's make, model, age, and mileage compared to other ML models, achieving out-of-sample R-squared 0.853 on a left-out test data set. We conjecture that the random forest performs best because it can better capture differential depreciation by make and model than the other methods.

Discussion
We have shown how the bestNormalize package can effectively and efficiently find the best normalizing transformation for a vector or set of vectors. However, normalization is by no means something that should be applied universally and without motivation. In situations where units have meaning, normalizing prior to analysis can contaminate the relationships suspected in the data and/or reduce predictive accuracy. Further, depending on the type of transformations used, interpreting regression coefficients post-transformation can be difficult or impossible without using a figure since the transformation function itself will look completely different for different distributions. So, while normalization transformations may well be able to increase the robustness of results and mitigate violations to the classical linear regression assumption of Gaussian residuals, it is by no means a universal solution.
On the other hand, when hypotheses are exploratory or when data is of poor quality with high amounts of skew/outliers, normalization can be an effective means of mitigating downstream issues this can cause in the analyses. For example, in machine learning contexts, some predictor manipulations rely on second-order statistics (e.g., principal components analysis or partial least squares), for which the variance calculation can be sensitive to skew and outliers. Normalizing transformations can improve the quality and stability of these calculations. Similarly, predictor normalization reduces the tendency for high-leverage points to have their leverage propagated into engineered features such as interactions or polynomials. Ultimately, these benefits can often produce predictive models that are more robust and stable.
We focused on making this package useful in a variety of machine learning workflows. We are enthusiastic in our support of bestNormalize, and will continue to maintain the package while it is found to be useful by R users. We hope to continue to build up the repertoire of candidate transformations using the same infrastructure so that additional ones can be considered by default in the future.