10
$\begingroup$

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?

enter image description here

$\endgroup$
14
$\begingroup$

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.

enter image description here

$\endgroup$
  • $\begingroup$ Thank you for your answer, but I am confused- you say 'estimate the spike using the proportion of zeros' but plot it with no bounds. does the spike have a discrete height or is it infinite, if discrete, is it $P(X=0)$? $\endgroup$ – Abe Jan 27 '11 at 17:45
  • $\begingroup$ This is a mixture of a discrete distribution and a continuous distribution. When plotted as a density, the spike is infinite (actually a Dirac delta function). Sometimes people plot the discrete part as a probability mass function (so then the spike has height $P(X=0)$) and the continuous part as a density function. That probably makes a better visual, but it does involve two different scales. $\endgroup$ – Rob Hyndman Jan 27 '11 at 23:10
  • $\begingroup$ this is coming in handy. fyi: it appears that, although bw = "SJ" affects the density in untransformed space, the logdensity is the same using "SJ" and the default "nrd0"... I'm about to read the SJ reference: "Sheather and Jones (1991) A reliable data-based bandwidth selection method for kernel density estimation." jstor.org/stable/2345597 $\endgroup$ – Abe Mar 29 '11 at 20:54
4
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

You may try lowering bandwidth (blue line is for adjust=0.5), enter image description here

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

$\endgroup$
  • $\begingroup$ is there another method that you would recommend? $\endgroup$ – Abe Jan 26 '11 at 22:22
  • $\begingroup$ @Abe Well, this depends on what you want to do... $\endgroup$ – user88 Jan 26 '11 at 22:26

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.