36
$\begingroup$
plot(density(rexp(100))

Obviously all density to the left of zero represents bias.

I'm looking to summarize some data for non-statisticians, and I want to avoid questions about why non-negative data has density to the left of zero. The plots are for randomization checking; I want to show the distributions of variables by treatment and control groups. The distributions are often exponential-ish. Histograms are tricky for various reasons.

A quick google search gives me work by statisticians on non-negative kernels, e.g.: this.

But has any of it been implemented in R? Of implemented methods, are any of them "best" in some way for descriptive statistics?

EDIT: even if the from command can solve my current problem, it'd be nice to know whether anyone has implemented kernels based on literature on non-negative density estimation

$\endgroup$
  • 3
    $\begingroup$ Not what you are asking, but I wouldn't apply kernel density estimation to something that should be exponential, especially for presentation to non-statistical audiences. I would use a quantile-quantile plot and explain that the plot should be straight if the distribution were exponential. $\endgroup$ – Nick Cox Jul 29 '13 at 7:28
  • 6
    $\begingroup$ plot(density(rexp(100), from=0)) ? $\endgroup$ – Stéphane Laurent Jul 29 '13 at 7:59
  • 4
    $\begingroup$ One thing that I have sometimes done fairly successfully is to get a kde on the logs, and then transform the density estimate (not forgetting the Jacobian). Another possibility would be to use a log-spline density estimate set up so it knows about the bound. $\endgroup$ – Glen_b -Reinstate Monica Jul 29 '13 at 9:52
  • 1
    $\begingroup$ possible duplicate of How can I estimate the density of a zero-inflated parameter in R? $\endgroup$ – Andy W Jul 29 '13 at 12:07
  • 1
    $\begingroup$ I discussed the transformation method mentioned by @Glen_b in stata-journal.com/sjpdf.html?articlenum=gr0003 (see pp.76-78). Zeros might be accommodated by using log(x + 1) rather than log and modifying the Jacobian. $\endgroup$ – Nick Cox Jul 29 '13 at 13:44
21
$\begingroup$

One solution, borrowed from approaches to edge-weighting of spatial statistics, is to truncate the density on the left at zero but to up-weight the data that are closest to zero. The idea is that each value $x$ is "spread" into a kernel of unit total area centered at $x$; any part of the kernel that would spill over into negative territory is removed and the kernel is renormalized to unit area.

For instance, with a Gaussian kernel $K_h(y,x) = \exp(-\frac{1}{2}((y-x)/h)^2) / \sqrt{2\pi}$, the renormalization weight is

$$w(x) = 1 / \int_0^\infty K(y,x) dy = \frac{1}{1 - \Phi_{x, h}(0)}$$

where $\Phi$ is the cumulative distribution function of a Normal variate of mean $x$ and standard deviation $h$. Comparable formulas are available for other kernels.

This is simpler--and much faster in computation--than trying to narrow the bandwidths near $0$. It is difficult to prescribe exactly how the bandwidths should be changed near $0$, anyway. Nevertheless, this method is also ad hoc: there will still be some bias near $0$. It appears to work better than the default density estimate. Here is a comparison using a largish dataset:

Figure

The blue shows the default density while the red shows the density adjusted for the edge at $0$. The true underlying distribution is traced as a dotted line for reference.


R code

The density function in R will complain that the sum of weights is not unity, because it wants the integral over all real numbers to be unity, whereas this approach makes the integral over positive numbers equal to unity. As a check, the latter integral is estimated as a Riemann sum.

set.seed(17)
x <- rexp(1000)
#
# Compute a bandwidth.
#
h <- density(x, kernel="gaussian")$bw # $
#
# Compute edge weights.
#
w <- 1 / pnorm(0, mean=x, sd=h, lower.tail=FALSE)
#
# The truncated weighted density is what we want.
#
d <- density(x, bw=h, kernel="gaussian", weights=w / length(x))
d$y[d$x < 0] <- 0
#
# Check: the integral ought to be close to 1:
#
sum(d$y * diff(d$x)[1])
#
# Plot the two density estimates.
#
par(mfrow=c(1,1))
plot(d, type="n", main="Default and truncated densities", xlim=c(-1, 5))
polygon(density(x, kernel="gaussian", bw=h), col="#6060ff80", border=NA)
polygon(d, col="#ff606080", border=NA)
curve(exp(-x), from=0, to=max(x), lty=2, add=TRUE)
$\endgroup$
21
$\begingroup$

An alternative is the approach of Kooperberg and colleagues, based on estimating the density using splines to approximate the log-density of the data. I'll show an example using the data from @whuber's answer, which will allow for a comparison of approaches.

set.seed(17)
x <- rexp(1000)

You'll need the logspline package installed for this; install it if it is not:

install.packages("logspline")

Load the package and estimate the density using the logspline() function:

require("logspline")
m <- logspline(x)

In the following, I assume that the object d from @whuber's answer is present in the workspace.

plot(d, type="n", main="Default, truncated, and logspline densities", 
     xlim=c(-1, 5), ylim = c(0, 1))
polygon(density(x, kernel="gaussian", bw=h), col="#6060ff80", border=NA)
polygon(d, col="#ff606080", border=NA)
plot(m, add = TRUE, col = "red", lwd = 3, xlim = c(-0.001, max(x)))
curve(exp(-x), from=0, to=max(x), lty=2, add=TRUE)
rug(x, side = 3)

The resulting plot is shown below, with the logspline density shown by the red line

Default, truncated, and logspline densities

Additionally, the support for the density can be specified via arguments lbound and ubound. If we wish to assume that the density is 0 to the left of 0 and there is a discontinuity at 0, we could use lbound = 0 in the call to logspline(), for example

m2 <- logspline(x, lbound = 0)

Yielding the following density estimate (shown here with the original m logspline fit as the previous figure was already getting busy).

plot.new()
plot.window(xlim = c(-1, max(x)), ylim = c(0, 1.2))
title(main = "Logspline densities with & without a lower bound",
      ylab = "Density", xlab = "x")
plot(m,  col = "red",  xlim = c(0, max(x)), lwd = 3, add = TRUE)
plot(m2, col = "blue", xlim = c(0, max(x)), lwd = 2, add = TRUE)
curve(exp(-x), from=0, to=max(x), lty=2, add=TRUE)
rug(x, side = 3)
axis(1)
axis(2)
box()

The resulting plot is shown below

Comparison of logspline density estimates with and without a lower bound on the support

In this case, exploiting knowledge of x results in a density estimate that doesn't tend to 0 at $x = 0$, but is similar to the standard logspline fit elsewhere over x

$\endgroup$
  • 1
    $\begingroup$ +1 I like this idea but am surprised it works so well with the example data (from an exponential distribution). Do you have some intuition as to why that is? In some sense it does great near $0$ but it is missing the "lump" of values near $1$ in the actual data, so I wonder whether there's not some kind of trade-off between good accuracy at low values and low accuracy (or, equivalently, large bandwidths) at high values. $\endgroup$ – whuber Sep 27 '13 at 19:58
  • $\begingroup$ @whuber Good question. I only came across this approach recently myself. I suspect a good question to ask here is, as the truncated and logspline methods are just estimates of the true density, are the differences in fit significant, statistically? I'm not sure exactly why it does so well at zero, though. I'd appreciate knowing why too. $\endgroup$ – Reinstate Monica - G. Simpson Sep 27 '13 at 20:09
  • $\begingroup$ @GavinSimpson, Thanks for this nice answer. Can you reproduce the last plot with the lastest version of logspline? For me, the density of both, the bounded and the unbounded version go to zero at x = 0. $\endgroup$ – cel Nov 21 '14 at 16:27
4
$\begingroup$

To compare distributions by groups (which you say is the goal in one of your comments) why not something simpler? Parallel box plots work nicely if N is large; parallel strip plots work if N is small (and both show outliers well, which you say is a problem in your data).

$\endgroup$
  • 1
    $\begingroup$ Yeah, thanks, that works. But I like density plots. They show more about the data than boxplots do. I guess I am kind of surprised that nothing seems to have been already implemented. Maybe I'll implement one of these things myself one day. People would probably find it useful. $\endgroup$ – generic_user Jul 29 '13 at 10:56
  • 1
    $\begingroup$ I like density plots, too; but you have to consider your audience. $\endgroup$ – Peter Flom - Reinstate Monica Jul 29 '13 at 11:27
  • 1
    $\begingroup$ Have to agree with @PeterFlom on this one. Don't get too complicated if your audience isn't statistically knowledgeable. You could also do comparative/parallel box-plots with an overlay of butterfly plots on top. That way the box-plot summary is visible as well as all the data. $\endgroup$ – doug.numbers Jul 29 '13 at 21:45
  • $\begingroup$ The suggestion that different people comprehend aggregated plots differently is certainly correct. Despite understanding what a density plot is (and understanding that it is not a probability) I have no understanding of what a "parallel boxplot" might be. It suggests a parallel coordinates plot but I suspect that is not correct. $\endgroup$ – DWin Oct 3 '13 at 17:14
2
$\begingroup$

As Stéphane comments you can use from = 0 and, additionally, you can represent your values under the density curve with rug (x)

$\endgroup$
  • 4
    $\begingroup$ Correct me if I am wrong but from=0 looks as if it just suppresses plotting for values below 0; it doesn't correct the calculation for the fact that some of the distribution has been smeared below 0. $\endgroup$ – Nick Cox Jul 29 '13 at 9:36
  • 1
    $\begingroup$ That is correct. Using the from command yields a plot that looks like it has the peak just right of zero. But if you look at histograms with continually smaller bins, a lot of data will show the peak AT zero. The from is just a graphical trick. $\endgroup$ – generic_user Jul 29 '13 at 10:26
  • $\begingroup$ @NickCox I'm not sure but I don't think from=0 suppresses anything. It just starts the "grid" at zero. $\endgroup$ – Stéphane Laurent Jul 29 '13 at 10:41
  • $\begingroup$ The difference is whether the estimated density is non-zero for negative values, not whether or not it is plotted. Researchers may decide not to worry about this if all they want is a visualization. $\endgroup$ – Nick Cox Jul 29 '13 at 13:36
  • $\begingroup$ @NickCox The command density(rexp(100), from=0) has nothing to do with the graphic $\endgroup$ – Stéphane Laurent Jul 29 '13 at 16:55

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.