Choosing c such that log(x + c) would remove skew from the population

Question

I have data for which I would like to take the log transformation before doing OLS. The data include zeros. Thus, I want to do a log(x + c). I know a traditional c to choose is 1. I am wondering whether there is a way to have the data choose c such that there would be no skew anymore using features like the sample mean or variance? Is there a formaula?

Answer 1

Because the intention is to do OLS, the choice of $c$ should be made in this context.

In general, we ought to fit $c$ simultaneously with the rest of the regression. A quick and dirty way to do this recognizes that the regression $R^2$ is proportional to the log likelihood, so we could seek a value of $c$ that maximizes $R^2$.

This is a special example of the problem of choosing among a parameterized family of transformations $y \to f(y; \theta)$ to achieve the best possible fit of $y$ to explanatory values $x$. This can be solved in R rather simply and directly:

xform <- function(f, theta, x, y, ...) {
  g <- function(theta) -summary(lm(f(y, theta) ~ x))$r.squared
  nlm(g, theta, ...)
}

(I am glossing over a somewhat delicate matter of choosing good starting values for the parameter: it is possible to obtain bad solutions with nlm otherwise. Standard methods of exploratory data analysis will produce decent starting values, but that's a subject for another day.)

As an example of the use of xform, let's generate some highly skewed data for $y$ for which the "started logarithm" $\log(y+c)$ will produce an unskewed distribution:

set.seed(17)
y <- sort(exp(rnorm(32, 4, 1))) + 100

Evidently $\log(y-100)$ is drawn iid from a Normal distribution.

I will apply xform to three choices of $x$:

1. Values from which $y$ differs by additive, homoscedastic residuals. In this case it would be a mistake to take the logarithm of $y$: it is a grossly incorrect model of the relationship between $x$ and $y$.

2. Values from which $y$ differs by multiplicative lognormal residuals (more or less). In this case, taking the logarithm of $y$ is a good idea because it leads to a model for which OLS regression is appropriate.

3. Constant values of $x$, so that in effect we are looking at $y$ outside of the regression context altogether.

In cases (1) and (2) I will plot the histogram of $y$ (to show it is highly skewed), the scatterplot of $y$ against $x$ (to exhibit the data), and the scatterplot of the transformed $y$ against $x$ with the OLS line superimposed, to see the result of the transformation. In the third case those scatterplots are meaningless, so I only report the value of $c$ found by xform.

set.seed(17)
par(mfrow=c(2,3))

y <- sort(exp(rnorm(32, 4, 1))) + 100
x <- y - rnorm(length(y), 0, 50)
hist(y)
plot(x,y)
const <- xform(function(y,c) log(y+c), 1-min(y), x, y)$estimate
plot(x, log(y + const), ylab=sprintf("log(y + %0.1f)", const))
abline(coef(fit1<-lm(log(y+const)~x)), col="Gray")

x <- log((1:length(y) + rnorm(length(y), 10, 3)))
hist(y)
plot(x,y)
const <- xform(function(y,c) log(y+c), 1-min(y), x, y)$estimate
plot(x, log(y + const), ylab=sprintf("log(y + %0.1f)", const))
abline(coef(fit2<-lm(log(y+const)~x)), col="Gray")

x <- rep(1, length(y))
const <- xform(function(y,c) log(y+c), 1-min(y), x, y)$estimate

The top row is the first case and the second row of plots are for the second case.

Please observe:

• The $y$ values are identical in all three instances.

• The $y$ values are constructed from a model in which $c=-100$.

• The fitted value of $c$ in the first case, $1408392.5$, is essentially infinite. This indicates it's bad to be taking the logarithm at all for these $x$ values. (Adding this huge value of $c$ to $y$ before taking the logarithm basically does not change the shape of the data: that's why the two scatterplots in the top row look the same.)

• The fitted value of $c$ in the second case, $-104.0$, is close to the value of $-100$ used to generate the data. (Repeated simulations indicate that the fitted value in the second case will be biased slightly low, averaging around $-105$.)

• The fitted value of $c$ in the third case is $-115.8$, still close to the value used to generate the data.

If we were to use the "universal" value of $c$ found in the third case (essentially by ignoring the $x$ values), here is what the scatterplot would look like in conjunction with the $x$ values from case 1:

For these particular $x$ values, the OLS line fit to the transformed $y$ values is a terrible description of the relationship between $y$ and $x$. Notice how it underestimates most values of $y$ but grossly overestimates a few of them for $x$ between $140$ and $200$.

In summary, if you want to transform the response variable for a regression (to achieve symmetry or linearity), you must account for the regression itself. This is because the regression only "cares" about the residuals, not the raw values of $y$. As the extraordinarily bad value of $c$ in the first case shows, ignoring this advice could produce awful results.

Answer 2

You could try to minimize the square of the sample skewness with respect to the constant $c$. Here is a quick R program to get the job done, though I'm not sure if this is the best way to go about solving your problem.

log.skewness computes the skewness of the transformed data and f returns the square of the skewness of the transformed data.

log.skewness <- function(c, x) {
    y <- log(x+c)
    mean((y-mean(y))^3)/(mean((y-mean(y))^2))^(3/2)
}

f <- function(c, x) {
    log.skewness(c, x)^2
}

Here I generate some data, minimize the squared skewness of the transformed data with respect to $c$, and plot the histograms of the original data and the transformed data.

set.seed(seed = 1)
data <- rexp(100)
opt <- optimize(f = f, interval = c(min(data), 1E8), x = data)
par(mfrow=c(1,2))
hist(x = data, main = 'Original Data', xlab = '')
hist(x = log(data + opt$minimum), main = 'Transformed data', xlab = '')
par(mfrow=c(1,1))

Answer 3

In a comment, I asserted there is a formula for a constant $c$ which eliminates the "skewness" of a dataset $(y_i)$ upon applying the "started logarithm" $y \to \log(y+c)$. There are many formulas. This answer describes a family of formulas that are extremely simple, fast to compute, robust, and flexible. (Note, though, that another answer I posted explains why you might not want to apply such a formula to the dependent variable in a regression context and gives an alternative for such cases.)

Motivation

In practice the term "skewness" in the questions usually refers qualitatively to the shape of dataset, not to the sample skewness. The sample skewness (a standardized central third moment) is highly sensitive to even moderate outliers. Procedures that are sensitive to outliers are good procedures only for--letting outliers ruin your analysis. Thus, sample skewness is usually a poor choice to use for exploring or characterizing data. (Yes, there do exist some statistic tests based on skewness, but they tend to be inferior to more robust tests or of limited application.)

Instead, pick some low percent $\alpha$--typically between $.05$ and $0.25$--and examine the $\alpha$ and $1-\alpha$ percentiles of the dataset $(y_i)$. Let's call these $q_{-}$ and $q_{+}$, respectively. Let the median of the dataset be $m$. A good measure of symmetry of the data is the difference between the upper spread $q_{+}-m$ and the lower spread $m - q_{-}$: in a perfectly symmetric data set, these two spreads around the median are equal no matter what the value of $\alpha$. Accommodating the reality that almost no dataset is, or can readily be made, symmetric, we anticipate that achieving symmetry of the spreads at the $\alpha$ and $1-\alpha$ percentiles will likely do a good job of symmetrizing the whole dataset, except possibly for some outlying values beyond those percentiles.

The formula

Here is the signal beauty of this approach: because we seek a monotonic re-expression of the data, say $y_i \to f(y_i,c)$ for a given function $f$ (the logarithm) and parameter $c$ yet to be found, the $\alpha$ percentile of the re-expressed data $(f(y_i,c))$ is (up to a tiny difference related to the discreteness of the data) equal to $f$ applied to the $\alpha$ percentile of $(y_i)$. Therefore, it will suffice to find $c$ so that

$$\log(q_{+} + c) - \log(m + c) = \log(m + c) - \log(q_{-} + c).$$

By means of elementary properties of logarithms, this equation implies

$$(q_{+} + c)(q_{-} + c) = (m + c)^2.$$

Provided the data are not already symmetric (that is, there actually is a difference in the original spreads), this has at most one solution

$$c = \frac{m^2 - q_{+}q_{-}}{q_{+} + q_{-} - 2m}.$$

Because a fixed quantile of $n$ data can be found with $O(n)$ computation, this is a fast procedure. For large datasets, applying this algorithm to a small random sample (as small as a few dozen values) will work fine, leading to an extremely fast constant time algorithm.

Implementation

There are some messy practical details to deal with if one wants a general-purpose software solution. In order to apply logarithms, we will need all of the $(y_i+c)$ to be positive. Thus, when the solution $c$ is less than or equal to $-\min(\{y_i\})$, some alternative must be found. One reasonable approach includes two things:

1. Return a code to indicate the original solution will not be applicable to all data values and

2. Adjust $c$ to be slightly greater than $-\min(\{y_i\})$.

The adjustment should be relative to the range of values in $(y_i)$. That is, we should select some other percent $\beta$ and extrapolate from the $\beta$ percentile of the $(y_i)$ down below $\min(\{y_i\}$ by a tiny proportion $\gamma \gt 0$. That is the solution adopted by the following R function.

log.start <- function(y, alpha=0.25, beta=1, gamma=0.01) {
  #
  # Find a constant `const` for which log(y+c) is as symmetric
  # as possible at the `alpha` and 1-alpha percentiles.
  # Returns a list containing the solution and a code:
  #   Code=0: Constant is good.
  #   Code=1: Data are already symmetric.
  #   Code=2: The best solution is not large enough and has been adjusted.
  #           In this case, the constant is obtained by extrapolating to the
  #           left from the `beta` percentile beyond the minimum value of `y`;
  #           the extrapolation is by a (small, positive) amount `gamma`.
  #
  stats <- quantile(y, c(alpha, 1/2, 1-alpha))
  code <- 0
  if (diff(diff(stats)==0)) {
    const <- Inf
    code <- 1
  }
  else {
    const <- (stats[2]^2 - stats[1]*stats[3]) / (stats[1] + stats[3] - 2*stats[2])
    y.min <- min(y)
    if (const < -y.min) {
      const <- gamma * quantile(y, beta) - (1+gamma) * y.min
      code <- 2
    }
  }
  list(offset=as.numeric(const), code=code)
}

Example

As an example of its use, let's generate data $(y_i)$ for which $\log(y_i+100)$ are iid Normal:

set.seed(17)
y <- exp(rnorm(32, 4, 1)) + 100

Apply the formula. Here I chose $\alpha=0.10$, in an effort to symmetrize the middle $1-2\alpha$ = $80$% of the data, allowing for as much as $10$% in the upper or lower tails (or both) to exhibit outlying behavior:

const <- log.start(y, .10)$offset

Display the distributions of the data (in the left column) and the transformed data (in the right column):

par(mfrow=c(2,2))
hist(y); hist(log(y + const))
qqnorm(y); qqnorm(log(y + const))

The almost symmetric histogram and nearly linear q-q plot in the right column of graphs show this choice of $c$ (which was $-97.7$, quite close to canceling the value of $100$ used to generate the data) works very well in this instance.