How can I estimate the density of a zero-inflated parameter in R?

Question

I have a data set with lots of zeros that looks like this:

set.seed(1)
x <- c(rlnorm(100),rep(0,50))
hist(x,probability=TRUE,breaks = 25)

I would like to draw a line for its density, but the density() function uses a moving window that calculates negative values of x.

lines(density(x), col = 'grey')

There is a density(... from, to) arguments, but these seem to only truncate the calculation, not alter the window so that the density at 0 is consistent with the data as can be seen by the following plot :

lines(density(x, from = 0), col = 'black')

(if the interpolation was changed, I would expect that the black line would have higher density at 0 than the grey line)

Are there alternatives to this function that would provide a better calculation of the density at zero?

Answer 1

The density is infinite at zero because it includes a discrete spike. You need to estimate the spike using the proportion of zeros, and then estimate the positive part of the density assuming it is smooth. KDE will cause problems at the left hand end because it will put some weight on negative values. One useful approach is to transform to logs, estimate the density using KDE, and then transform back. See Wand, Marron & Ruppert (JASA 1991) for a reference.

The following R function will do the transformed density:

logdensity <- function (x, bw = "SJ") 
{
    y <- log(x)
    g <- density(y, bw = bw, n = 1001)
    xgrid <- exp(g$x)
    g$y <- c(0, g$y/xgrid)
    g$x <- c(0, xgrid)
    return(g)
}

Then the following will give the plot you want:

set.seed(1)
x <- c(rlnorm(100),rep(0,50))
hist(x,probability=TRUE,breaks = 25)
fit <- logdensity(x[x>0]) # Only take density of positive part
lines(fit$x,fit$y*mean(x>0),col="red") # Scale density by proportion positive
abline(v=0,col="blue") # Add spike at zero.

Answer 2

I'd agree with Rob Hyndman that you need to deal with the zeroes separately. There are a few methods of dealing with a kernel density estimation of a variable with bounded support, including 'reflection', 'rernormalisation' and 'linear combination'. These don't appear to have been implemented in R's density function, but are available in Benn Jann's kdens package for Stata.

Answer 3

Another option when you have data with a logical lower bound (such as 0, but could be other values) that you know the data will not go below and the regular kernel density estimate places values below that bound (or if you have an upper bound, or both) is to use logspline estimates. The logspline package for R implements these and the functions have arguments for specifying the bounds so the estimate will go to the bound, but not beyond and still scale to 1.

There are also methods (the oldlogspline function) that will take into account interval censoring, so if those 0's are not exact 0's, but are rounded so that you know they represent values between 0 and some other number (a detection limit for example) then you can give that information to the fitting function.

If the extra 0's are true 0's (not rounded) then estimating the spike or point mass is the better approach, but can also be combined with logspline estimation.

Answer 4

You may try lowering bandwidth (blue line is for adjust=0.5),

but probably KDE is just not the best method to deal with such data.