Good methods for density plots of non-negative variables in R?

Question

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

Answer 1

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

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

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

Answer 2

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:

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)

Answer 3

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

Answer 4

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