How small a quantity should be added to x to avoid taking the log of zero?

Question

I have analysed my data as they are. Now I want to look at my analyses after taking the log of all variables. Many variables contain many zeros. Therefore I add a small quantity to avoid taking the log of zero.

So far I've added 10^-10, without any rationale really, just because I felt like adding a very small quantity would be advisable to minimize the effect of my arbitrarily chosen quantity. But some variables contain mostly zeros, and therefore when logged mostly -23.02. The range of the ranges of my variables is 1.33-8819.21, and the frequency of zeroes also varies dramatically. Therefore my personal choice of "small quantity" affects the variables very differently. It is clear now that 10^-10 is a completely unacceptable choice, as most of the variance in all the variables then comes from this arbitrary "small quantity".

I wonder what would be a more correct way of doing this.

Maybe it's better to derive the quantity from each variables individual distribution? Are there any guidelines about how big this "small quantity" should be?

My analyses are mostly simple cox models with each variable and age/sex as IVs. The variables are the concentrations of various blood lipids, with often considerable coefficients of variation.

Edit: Adding the smallest non-zero value of the variable seems practical for my data. But maybe there is a general solution?

Edit 2: As the zeros merely indicate concentrations below the detection limit, maybe setting them to (detection limit)/2 would be appropriate?

Answer 1

As the zeros merely indicate concentrations below the detection limit, maybe setting them to (detection limit)/2 would be appropriate

I was just typing that the thing that comes to my mind where log does (frequently) make sense and 0 may occur are concentrations when you did the 2nd edit. As you say, for measured concentrations the 0 just means "I couldn't measure that low concentrations".

Side note: do you mean LOQ instead of LOD?

Whether setting the 0 to $\frac{1}{2}$LOQ is a good idea or not depends:

• from the point of view that $\frac{1}{2}\mathrm{LOQ}$ is your "guess" expressing that c is anywhere between 0 and LOQ, it does make sense.
  But consider the corresponding calibration function:
  
  On the left, the calibration function yields c = 0 below the LOQ. On the right, $\frac{1}{2}\mathrm{LOQ}$ is used instead of 0.

• However, if the original measured value is available, that may provide a better guess. After all, LOQ usually just means that the relative error is 10%. Below that the measurement still carries information, but the relative error becomes huge.
  
  (blue: LOD, red: LOQ)

• An alternative would be to exclude these measurements. That can be reasonable, too
  e.g. think of a calibration curve. In practice you often observe a sigmoid shape: for low c, signal ≈ constant, intermediate linear behaviour, then detector saturation.
  In that situation you may want to restrict yourself to statements about concentrations that are clearly in the linear range as both below and above other processes heavily influence the result.
  Make sure you explain that the data was selected that way and why.

---

edit: What is sensible or acceptable, depends of course on the problem. Hopefully, we're talking here about a small part of the data that does not influence the analyis.

Maybe a quick and dirty check is: run your data analysis with and without excluding the data (or whatever treatment you propose) and see whether anything changes substantially.

If you see changes, then of course you're in trouble. However, from the analytical chemistry point of view, I'd say your trouble does not primarily lie in which method you use to deal with the data, but the underlying problem is that the analytical method (or its working range) was not appropriate for the problem at hand. There is of course a zone where the better statistical approach can save your day, but in the end the approximation "garbage in, garbage out" usually holds also for the more fancy methods.

Quotations for the topic:

• A statistician once told me:

  The problem with you (chemists/spectroscopists) is that your problems are either so hard that they cannot be solved or so easy that there is no fun in solving them.

• Fisher about the statistical post-mortem of experiments

Answer 2

Chemical concentration data often have zeros, but these do not represent zero values: they are codes that variously (and confusingly) represent both nondetects (the measurement indicated, with a high degree of likelihood, that the analyte was not present) and "unquantified" values (the measurement detected the analyte but could not produce a reliable numeric value). Let's just vaguely call these "NDs" here.

Typically, there is a limit associated with an ND variously known as a "detection limit," "quantitation limit," or (much more honestly) a "reporting limit," because the laboratory chooses not to provide a numerical value (often for legal reasons). About all we really know of an ND is that the true value is likely less than the associated limit: it's almost (but not quite) a form of left censoring. (Well, that's not really true either: it's a convenient fiction. These limits are determined via calibrations which, in most cases, have poor to terrible statistical properties. They may be grossly over- or under-estimated. This is important to know when you're looking at a set of concentration data which appear to have a lognormal right tail which is cut off (say) at $1.33$, plus a "spike" at $0$ representing all the NDs. That would strongly suggest the reporting limit is just a little less than $1.33$, but the lab data might try to tell you it is $0.5$ or $0.1$ or something like that.)

Extensive research has been done over the last 30 years or so concerning how best to summarize and evaluate such datasets. Dennis Helsel published a book on this, Nondetects and Data Analysis (Wiley, 2005), teaches a course, and released an R package based on some of the techniques he favors. His website is comprehensive.

This field is fraught with error and misconception. Helsel is frank about this: on the first page of chapter 1 of his book he writes,

...the most commonly used method in environmental studies today, substitution of one-half the detection limit, is NOT a reasonable method for interpreting censored data.

So, what to do? Options include ignoring this good advice, applying some of the methods in Helsel's book, and using some alternative methods. That's right, the book is not comprehensive and valid alternatives do exist. Adding a constant to all values in the dataset ("starting" them) is one. But consider:

• Adding $1$ is not a good place to start, because this recipe depends on the units of measurement. Adding $1$ microgram per deciliter will not have the same result as adding $1$ millimole per liter.

• After starting all the values, you will still have a spike at the smallest value, representing that collection of NDs. Your hope is that this spike is consistent with the quantified data in the sense that its total mass is approximately equal to the mass of a lognormal distribution between $0$ and the start value.

  An excellent tool for determining the start value is a lognormal probability plot: apart from the NDs, the data should be approximately linear.

• The collection of NDs can also be described with a so-called "delta lognormal" distribution. This is a mixture of a point mass and a lognormal.

As is evident in the following histograms of simulated values, the censored and delta distributions are not the same. The delta approach is most useful for explanatory variables in regression: you can create a "dummy" variable to indicate the NDs, take logarithms of the detected values (or otherwise transform them as needed), and not worry about the replacement values for the NDs.

In these histograms, approximately 20% of the lowest values have been replaced by zeros. For comparability, they are all based on the same 1000 simulated underlying lognormal values (upper left). The delta distribution was created by replacing 200 of the values by zeros at random. The censored distribution was created by replacing the 200 smallest values by zeros. The "realistic" distribution conforms to my experience, which is that the reporting limits actually vary in practice (even when that is not indicated by the laboratory!): I made them vary randomly (by just a little bit, rarely more than 30 in either direction) and replaced all simulated values less than their reporting limits by zeros.

To show the utility of the probability plot and to explain its interpretation, the next figure displays normal probability plots related to the logarithms of the preceding data.

The upper left shows all the data (before any censoring or replacement). It's a good fit to the ideal diagonal line (we expect some deviations in the extreme tails). This is what we are aiming to achieve in all the subsequent plots (but, due to the NDs, we will inevitably fall short of this ideal.) The upper right is a probability plot for the censored dataset, using a start value of 1. It's a terrible fit, because all the NDs (plotted at 0, because $\log(1+0)=0$) are plotted much too low. The lower left is a probability plot for the censored dataset with a start value of 120, which is close to a typical reporting limit. The fit in the bottom left is now decent--we only hope that all these values come somewhere near to, but to the right of, the fitted line--but the curvature in the upper tail shows that adding 120 is starting to alter the shape of the distribution. The bottom right shows what happens to the delta-lognormal data: there's a good fit to the upper tail, but some pronounced curvature near the reporting limit (at the middle of the plot).

Finally, let's explore some of the more realistic scenarios:

The upper left shows the censored dataset with the zeros set to one-half the reporting limit. It's a pretty good fit. On the upper right is the more realistic dataset (with randomly varying reporting limits). A start value of 1 does not help, but--on the lower left--for a start value of 120 (near the upper range of the reporting limits) the fit is quite good. Interestingly, the curvature near the middle as the points rise up from the NDs to the quantified values is reminiscent of the delta lognormal distribution (even though these data were not generated from such a mixture). On the lower right is the probability plot you get when the realistic data have their NDs replaced by one-half the (typical) reporting limit. This is the best fit, even though it shows some delta-lognormal-like behavior in the middle.

What you ought to do, then, is to use probability plots to explore the distributions as various constants are used in place of the NDs. Start the search with one-half the nominal, average, reporting limit, then vary it up and down from there. Choose a plot that looks like the bottom right: roughly a diagonal straight line for the quantified values, a quick drop-off to a low plateau, and a plateau of values that (just barely) meets the extension of the diagonal. However, following Helsel's advice (which is strongly supported in the literature), for actual statistical summaries, avoid any method that replaces the NDs by any constant. For regression, consider adding in a dummy variable to indicate the NDs. For some graphical displays, the constant replacement of NDs by the value found with the probability plot exercise will work well. For other graphical displays it may be important to depict the actual reporting limits, so replace the NDs by their reporting limits instead. You need to be flexible!

Answer 3

@miura

I came across this article by Bill Gould on the Stata blog (I think he actually founded Stata) which I think could provide help with your analysis. Near the end of the article he cautions against the use of arbitrary numbers that are close to zero, such as 0.01, 0.0001, 0.0000001, and 0 since in logs they are -4.61, -9.21, -16.12, and $-\infty$. In this situation they are not arbitrary at all. He advises the use of a Poisson regression since it recognizes that the above number are actually close together.

Answer 4

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

Answer 5

You can set the zeros of the $i^{th}$ variable to the ${\rm mean}(x_i) - n\times{\rm stddev}(x_i)$ where $n$ is large enough to distinguish these cases from the rest (e.g., 6 or 10).

Note that any such artificial setup will affect your analyses so you should be careful with your interpretation and in some cases discard these cases to avoid artifacts.

Using the detection limit is also a reasonable idea.