8
$\begingroup$

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?

$\endgroup$
  • 5
    $\begingroup$ The question is good, but it is predicated on several misunderstandings. The first is that removing the skewness in the values is desirable for OLS. It is not: the objective would be to reduce the skewness of the OLS residuals, not of the data themselves. The second misunderstanding is that taking a logarithm is a good approach to analyzing data that include zeros. It might or might not be, depending on how those zeros arise. (This question has been addressed in several past threads.) Yes, of course there is a formula to find such a $c$--but applying it may make things worse, not better. $\endgroup$ – whuber Oct 28 '12 at 15:05
  • $\begingroup$ Yes, log might not be the best approach, but I am going to go with it. What is the formula? $\endgroup$ – user1690130 Oct 28 '12 at 15:13
  • $\begingroup$ I posted a separate answer providing the formula and software to implement it. $\endgroup$ – whuber Oct 29 '12 at 10:51
10
$\begingroup$

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

Plots

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:

Scatterplot

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.

| cite | improve this answer | |
$\endgroup$
8
$\begingroup$

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))

Distribution plots

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.

| cite | improve this answer | |
$\endgroup$
7
$\begingroup$

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))
| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you for this code! So, in words, this is saying to choose c such that it minimizes the skewness formula, where skewness is based on the third moment divided by (the second moment to the 1.5 power)? $\endgroup$ – user1690130 Oct 28 '12 at 15:28
  • $\begingroup$ @JG I just edited the post. I meant to write that we minimize the squared skewness of the transformed data because the best possible skewness is zero. For more information on skewness, see the wikipedia article. $\endgroup$ – assumednormal Oct 28 '12 at 15:32

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.