If I have highly skewed positive data I often take logs. But what should I do with highly skewed non-negative data that include zeros? I have seen two transformations used:
Are there any other approaches? Are there any good reasons to prefer one approach over the others?
It seems to me that the most appropriate choice of transformation is contingent on the model and the context.
The '0' point can arise from several different reasons each of which may have to be treated differently:
I am not really offering an answer as I suspect there is no universal, 'correct' transformation when you have zeros.
No-one mentioned the inverse hyperbolic sine transformation. So for completeness I'm adding it here.
This is an alternative to the Box-Cox transformations and is defined by \begin{equation} f(y,\theta) = \text{sinh}^{-1}(\theta y)/\theta = \log[\theta y + (\theta^2y^2+1)^{1/2}]/\theta, \end{equation} where $\theta>0$. For any value of $\theta$, zero maps to zero. There is also a two parameter version allowing a shift, just as with the two-parameter BC transformation. Burbidge, Magee and Robb (1988) discuss the IHS transformation including estimation of $\theta$.
The IHS transformation works with data defined on the whole real line including negative values and zeros. For large values of $y$ it behaves like a log transformation, regardless of the value of $\theta$ (except 0). The limiting case as $\theta\rightarrow0$ gives $f(y,\theta)\rightarrow y$.
It looks to me like the IHS transformation should be a lot better known than it is.
A useful approach when the variable is used as an independent factor in regression is to replace it by two variables: one is a binary indicator of whether it is zero and the other is the value of the original variable or a re-expression of it, such as its logarithm. This technique is discussed in Hosmer & Lemeshow's book on logistic regression (and in other places, I'm sure). Truncated probability plots of the positive part of the original variable are useful for identifying an appropriate re-expression. (See the analysis at https://stats.stackexchange.com/a/30749/919 for examples.)
When the variable is the dependent one in a linear model, censored regression (like Tobit) can be useful, again obviating the need to produce a started logarithm. This technique is common among econometricians.
The log transforms with shifts are special cases of the Box-Cox transformations:
$y(\lambda_{1}, \lambda_{2}) = \begin{cases} \frac {(y+\lambda_{2})^{\lambda_1} - 1} {\lambda_{1}} & \mbox{when } \lambda_{1} \neq 0 \\ \log (y + \lambda_{2}) & \mbox{when } \lambda_{1} = 0 \end{cases}$
These are the extended form for negative values, but also applicable to data containing zeros. Box and Cox (1964) presents an algorithm to find appropriate values for the $\lambda$'s using maximum likelihood. This gives you the ultimate transformation.
A reason to prefer Box-Cox transformations is that they're developed to ensure assumptions for the linear model. There's some work done to show that even if your data cannot be transformed to normality, then the estimated $\lambda$ still lead to a symmetric distribution.
I'm not sure how well this addresses your data, since it could be that $\lambda = (0, 1)$ which is just the log transform you mentioned, but it may be worth estimating the requried $\lambda$'s to see if another transformation is appropriate.
In R, the boxcox.fit function in package geoR will compute the parameters for you.
boxcox function in MASS only estimates $\lambda_1$ and assumes $\lambda_2=0$.
$\endgroup$
Commented
Aug 10, 2010 at 1:27
lambda2=TRUE in the arguments to boxcox.fit. (Also updated the answer.)
$\endgroup$
boxcoxfit.
$\endgroup$
I'm presuming that zero != missing data, as that's an entirely different question.
When thinking about how to handle zeros in multiple linear regression, I tend to consider how many zeros do we actually have?
Only a couple of zeros
If I have a single zero in a reasonably large data set, I tend to:
Does the model fit change? What about the parameter values? If the model is fairly robust to the removal of the point, I'll go for quick and dirty approach of adding $c$.
You could make this procedure a bit less crude and use the boxcox method with shifts described in ars' answer.
Large number of zeros
If my data set contains a large number of zeros, then this suggests that simple linear regression isn't the best tool for the job. Instead I would use something like mixture modelling (as suggested by Srikant and Robin).
If you want something quick and dirty why not use the square root?
Comparing the answer provided in by @RobHyndman to a log-plus-one transformation extended to negative values with the form:
$$T(x) = \text{sign}(x) \cdot \log{\left(|x|+1\right)} $$
(As Nick Cox pointed out in the comments, this is known as the 'neglog' transformation)
r = -1000:1000
l = sign(r)*log1p(abs(r))
l = l/max(l)
plot(r, l, type = "l", xlab = "Original", ylab = "Transformed", col = adjustcolor("red", alpha = 0.5), lwd = 3)
#We scale both to fit (-1,1)
for(i in exp(seq(-10, 100, 10))){
s = asinh(i*r)
s = s / max(s)
lines(r, s, col = adjustcolor("blue", alpha = 0.2), lwd = 3)
}
legend("topleft", c("asinh(x)", "sign(x) log(abs(x)+1)"), col = c("blue", "red"), lty = 1)
As you can see, as $\theta$ increases more the transform looks like a step function. With $\theta \approx 1$ it looks a lot like the log-plus-one transformation. And when $\theta \rightarrow 0$ it approaches a line.
EDIT: Keep in mind the log transform can be similarly altered to arbitrary scale, with similar results. I just wanted to show what $\theta$ gives similar results based on the previous answer. The biggest difference between both approaches is the region near $x=0$, as we can see by their derivatives.
I assume you have continuous data.
If the data include zeros this means you have a spike on zero which may be due to some particular aspect of your data. It appears for example in wind energy, wind below 2 m/s produce zero power (it is called cut in) and wind over (something around) 25 m/s also produce zero power (for security reason, it is called cut off). While the distribution of produced wind energy seems continuous there is a spike in zero.
My solution: In this case, I suggest to treat the zeros separately by working with a mixture of the spike in zero and the model you planned to use for the part of the distribution that is continuous (wrt Lebesgue).
Since the two-parameter fit Box-Cox has been proposed, here's some R to fit input data, run an arbitrary function on it (e.g. time series forecasting), and then return the inverted output:
# Two-parameter Box-Cox function
boxcox.f <- function(x, lambda1, lambda2) {
if (lambda1!=0) {
return(((x + lambda2) ^ lambda1 - 1) / lambda1)
} else {
return(log(x + lambda2))
}
}
# Two-parameter inverse Box-Cox function
boxcox.inv <- function(x, lambda1, lambda2) {
if (lambda1!=0) {
return((lambda1 * x + 1) ^ (1 / lambda1) - lambda2)
} else {
return(exp(x) - lambda2)
}
}
# Function to Box-Cox transform x, apply function g,
# and return inverted Box-Cox output y
boxcox.fit.apply <- function(x, g) {
require(geoR)
require(plyr)
# Fit lambdas
t <- try(lambda.pair <- boxcoxfit(x, lambda2=T)$lambda)
# Estimating both lambdas sometimes fails; if so, estimate lambda1 only
if (inherits(t, "try-error")) {
lambda1 <- boxcoxfit(x)$lambda
lambda2 <- 0
} else {
lambda1 <- lambda.pair[1]
lambda2 <- lambda.pair[2]
}
x.boxcox <- boxcox.f(x, lambda1, lambda2)
# Apply function g to x.boxcox. This should return data similar to x (e.g. ts)
y <- aaply(x.boxcox, 1, g)
return(boxcox.inv(y, lambda1, lambda2))
}
The Yeo-Johnson power transformation discussed here has excellent properties designed to handle zeros and negatives while building on the strengths of Box Cox power transformation. This is what I typically go to when I am dealing with zeros or negative data.
Here is a summary of transformations with pros/cons to illustrate why Yeo-Johnson is preferable.
Log
Pros: Does well with positive data.
Cons: Does not handle zeros.
> log(0)
[1] -Inf
Log Plus 1
Pros: The plus 1 offset adds the ability to handle zeros in addition to positive data.
Cons: Fails with negative data
> log1p(-1)
[1] -Inf
> log1p(-2)
[1] NaN
Warning message:
In log1p(-2) : NaNs produced
Square Root
Pros: Uses a power transformation that can handle zeros and positive data.
Cons: Fails with negative data
> sqrt(-1)
[1] NaN
Warning message:
In sqrt(-1) : NaNs produced
Box Cox
R Code:
box_cox <- function(x, lambda) {
eps <- 0.00001
if (abs(lambda) < eps)
log(x)
else
(x ^ lambda - 1) / lambda
}
Pros: Enables scaled power transformations
Cons: Suffers from issues with zeros and negatives (i.e. can only handle positive data.
> box_cox(0, lambda = 0)
[1] -Inf
> box_cox(0, lambda = -0.5)
[1] -Inf
> box_cox(-1, lambda = 0.5)
[1] NaN
Yeo Johnson
R Code:
yeo_johnson <- function(x, lambda) {
eps <- .000001
not_neg <- which(x >= 0)
is_neg <- which(x < 0)
not_neg_trans <- function(x, lambda) {
if (abs(lambda) < eps) log(x + 1)
else ((x + 1) ^ lambda - 1) / lambda
}
neg_trans <- function(x, lambda) {
if (abs(lambda - 2) < eps) - log(-x + 1)
else - ((-x + 1) ^ (2 - lambda) - 1) / (2 - lambda)
}
x[not_neg] <- not_neg_trans(x[not_neg], lambda)
x[is_neg] <- neg_trans(x[is_neg], lambda)
return(x)
}
Pros: Can handle positive, zero, and negative data.
Cons: None that I can think of. Properties are very similar to Box-Cox but can handle zero and negative data.
> yeo_johnson(0, lambda = 0)
[1] 0
> yeo_johnson(0, lambda = -0.5)
[1] 0
> yeo_johnson(-1, lambda = 0.5)
[1] -1.218951
Suppose Y is the amount of money each American spends on a new car in a given year (total purchase price). Y will spike at 0; will have no values at all between 0 and about 12,000; and will take other values mostly in the teens, twenties and thirties of thousands. Predictors would be proxies for the level of need and/or interest in making such a purchase. Need or interest could hardly be said to be zero for individuals who made no purchase; on these scales non-purchasers would be much closer to purchasers than Y or even the log of Y would suggest. In a case much like this but in health care, I found that the most accurate predictions, judged by test-set/training-set crossvalidation, were obtained by, in increasing order,
Some will recoil at this categorization of a continuous dependent variable. But although it sacrifices some information, categorizing seems to help by restoring an important underlying aspect of the situation -- again, that the "zeroes" are much more similar to the rest than Y would indicate.
To clarify how to deal with the log of zero in regression models, we have written a pedagogical paper explaining the best solution and the common mistakes people make in practice. We also came out with a new solution to tackle this issue.
You can find the paper by clicking here: https://ssrn.com/abstract=3444996
First, we think that ones should wonder why using a log transformation. In regression models, a log-log relationship leads to the identification of an elasticity. Indeed, if $\log(y) = \beta \log(x) + \varepsilon$, then $\beta$ corresponds to the elasticity of $y$ to $x$. The log can also linearize a theoretical model. It can also be used to reduce heteroskedasticity. However, in practice, it often occurs that the variable taken in log contains non-positive values.
A solution that is often proposed consists in adding a positive constant c to all observations $Y$ so that $Y + c > 0$. However, contrary to linear regressions, log-linear regressions are not robust to linear transformation of the dependent variable. This is due to the non-linear nature of the log function. Log transformation expands low values and squeezes high values. Therefore, adding a constant will distort the (linear) relationship between zeros and other observations in the data. The magnitude of the bias generated by the constant actually depends on the range of observations in the data. For that reason, adding the smallest possible constant is not necessarily the best worst solution.
In our article, we actually provide an example where adding very small constants is actually providing the highest bias. We provide derive an expression of the bias.
Actually, Poisson Pseudo Maximum Likelihood (PPML) can be considered as a good solution to this issue. One has to consider the following process:
$y_i = a_i \exp(\alpha + x_i' \beta)$ with $E(a_i | x_i) = 1$
This process is motivated by several features. First, it provides the same interpretation to $\beta$ as a semi-log model. Second, this data generating process provides a logical rationalization of zero values in the dependent variable. This situation can arise when the multiplicative error term, $a_i$ , is equal to zero. Third, estimating this model with PPML does not encounter the computational difficulty when $y_i = 0$. Under the assumption that $E(a_i|x_i) = 1$, we have $E( y_i - \exp(\alpha + x_i' \beta) | x_i) = 0$. We want to minimize the quadratic error of this moment, leading to the following first-order conditions:
$\sum_{i=1}^N ( y_i - \exp(\alpha + x_i' \beta) )x_i' = 0$
These conditions are defined even when $y_i = 0$. These first-order conditions are numerically equivalent to those of a Poisson model, so it can be estimated with any standard statistical software.
Finally, we propose a new solution that is also easy to implement and that provides unbiased estimator of $\beta$. One simply need to estimate:
$\log( y_i + \exp (\alpha + x_i' \beta)) = x_i' \beta + \eta_i $
We show that this estimator is unbiased and that it can simply be estimated with GMM with any standard statistical software. For instance, it can be estimated by executing just one line of code with Stata.
We hope that this article can help and we'd love to get feedback from you.
Christophe Bellégo and Louis-Daniel Pape CREST - Ecole Polytechnique - ENSAE
I had the same problem with data and no transformation would give reasonable distribution. I came up with the following idea. I would appreciate if someone decide whether it is worth utilising as I am not a statistitian.
We may adopt the assumption that 0 is not equal to 0. There is a hidden continuous value which we observe as zeros but, the low sensitivity of the test gives any values more than 0 only after reaching the treshold. So what we observe is more like half-normal distribution where all the left side of normal distribution is shown as one rectangle (x=0) in histogram. So maybe we can just perform following steps:
Depending on the problem's context, it may be useful to apply quantile transformations.
The idea itself is simple*, given a sample $x_1, \dots, x_n$, compute for each $i \in \{1, \dots, n\}$ the respective empirical cumulative density function values $F(x_i) = c_i$, then map $c_i$ to another distribution via the quantile function $Q$ of that distribution, i.e., $Q(c_i)$.
*Assuming you don't apply any interpolation and bounding logic.
With the method out of the way, there are several caveats, features, and notes which I will list below (mostly caveats). Details can be found in the references at the end.
Related questions and references:
