0

I'm trying to plot kernel density estimates for my data, but I'm getting an error ("bandwidths must be strictly positive") because the 0.25 and 0.75 quantiles are the same. A solution for this in MASS::kde2d was presented here, but I'm unsure how to write a function for h in stat_density_2d() that would work similarly.

The solution for kde2d is:

kde2d(df$s1x, df$s1y,
      h = c(ifelse(bandwidth.nrd(df$s1x) == 0, 0.1, bandwidth.nrd(df$s1x)),
            ifelse(bandwidth.nrd(df$s1y) == 0, 0.1, bandwidth.nrd(df$s1y))))

Because I'm faceting the plot, I need "h" to take a function that will work dynamically with facets or else I would add static values to "h".

Update: I've tried

h = c(function(x) ifelse(bandwidth.nrd(x)==0,0.01,bandwidth.nrd(x)),
      function(y) ifelse(bandwidth.nrd(y)==0,0.01,bandwidth.nrd(y)))

and

h= function(x) ifelse(x==c(0,0),c(0.01,0.01),x)

as well as a few variations, but without knowing exactly what's happening under the hood, I'm unsure how to code the function.

Here is some data I generated that may help:

xx <- rbind(data.frame(x = c(rep(30.00000,40),rep(30.00000,10)+rnorm(10)),
                 y = c(rep(-90.00000,40),rep(-90.00000,10)+rnorm(10)),
                 grp = 'A'),
            data.frame(x = c(rep(30.00000,50)+rnorm(50)),
                       y = c(rep(-90.00000,50)+rnorm(50)),
                       grp = 'B'))
ggplot()+
  stat_density_2d(data = xx, aes(x,y, fill = ..level..), geom = 'polygon')+
  facet_wrap(~grp)

example.png

jared_mamrot
  • 22,354
  • 4
  • 21
  • 46
user100351
  • 11
  • 4
  • 1
    Can you provide data and code that produces the error? Dummy data is fine as long as it generates the error. – teunbrand Jan 11 '22 at 22:22
  • It looks like `h` is an argument to `stat_density_2d()`, why can't you just supply it in the same way? – DaveArmstrong Jan 11 '22 at 23:44
  • I posted an answer suggesting that, but it doesn't work with the requirement to work dynamically with facets, i.e. you'd like a solution where the `x` and `y` referenced by `h` were *facet-specific*. I haven't figured out how to do it yet. – Ben Bolker Jan 11 '22 at 23:49

0 Answers0