R: Shadindg under several normals with ggplot

Question

I would like to plot several normal distributions and shade the area under every curve to the right of x=6 and to the left of the symetrical value...

I've defined two functions to do it easily.

mifun <- function(x,mm,ss,alt) {
  y <- alt+dnorm(x, mean = mm, sd = ss)
  return(y)
}

mifun2 <- function(x,mm,ss,alt) {
  y <- alt+dnorm(x, mean = mm, sd = ss)
  y[x<6 & (2*mm-x)<6 ] <- NA
  return(y)
}


ggplot(data.frame(x = c(-5, 11)), aes(x)) +
  stat_function(fun = function(x) mifun(x,0,2,0), geom = "line") +
  stat_function(fun = function(x) mifun2(x,0,2,0), geom = "area") +
  stat_function(fun = function(x) mifun(x,2,2,-0.4), geom = "line") +
  stat_function(fun = function(x) mifun2(x,2,2,-0.4), geom = "area") +
  stat_function(fun = function(x) mifun(x,4,2,-0.8), geom = "line") + 
  stat_function(fun = function(x) mifun2(x,4,2,-0.8), geom = "area") +
  stat_function(fun = function(x) mifun(x,6,2,-1.2), geom = "line") +
  stat_function(fun = function(x) mifun2(x,6,2,-1.2), geom = "area") +theme_bw()

This is the picture without shading:

This is the picture with the shading, it doesn't work as expected:

If I only plot one normal I get this

As you can see the shade appears above the normals instead of below.

PD: picture for Mark As you can see the first plots are aligned on the right shadow limit, and the latter ones are aligned on the left.

Answer 1

I think that the comments above are right: faceting is going to be the way to go. Unfortunately, facet_wrap (and facet_grid) don't play well with stat_function because there is no good way to pass the faceting variable into stat_function (that I can find anyway).

So, instead, you may need to generate your density curves and area to fill first. Note that geom_area fills in the gaps between points, so you need to have mifun2 output 0 instead of NA where there is supposed to be no fill:

mifun2 <- function(x,mm,ss,alt) {
  y <- alt+dnorm(x, mean = mm, sd = ss)
  y[x<6 & (2*mm-x)<6 ] <- 0
  return(y)
}

Then, generate the gridded data.frame you want. Here I am using 4 values for the means (from your question) and 1001 points in the range you want to plot:

normCurves <-
  data.frame(x = rep(seq(-5, 11, length.out = 1001)
                     , times = 4)
             , myMean = rep(c(0,2,4,6)
                            , each = 1001))

Then, use your two functions to generate plottable columns for the density and area filling:

normCurves$density <-
  mifun(normCurves$x, normCurves$myMean, 2, 0)

normCurves$toHighlight <-
  mifun2(normCurves$x, normCurves$myMean, 2, 0)

And then plot those directly, using facet_wrap to separate out the various means:

ggplot(normCurves
       , aes(x)) +
  geom_line(aes(y = density)) +
  geom_area(aes(y = toHighlight)) +
  facet_wrap(~myMean, ncol = 1, labeller = label_both)

gives:

This approach also lends itself to plotting everything on top of each other with different colors:

ggplot(normCurves
       , aes(x
             , col = factor(myMean)
             , fill = factor(myMean)
             , group = myMean)) +
  geom_line(aes(y = density)) +
  geom_area(aes(y = toHighlight)
            , alpha = 0.2
            , col = NA
            , position = "identity")

gives

If you want it to highlight things that are more extreme than 6 in the other direction (WHY??), you just need to change your mifun2 definition:

mifun2 <- function(x,mm,ss,alt) {
  y <- alt+dnorm(x, mean = mm, sd = ss)

  # Make changes for small means
  y[mm <= 6 & x<6 & (2*mm-x)<6 ] <- 0  

  # Make changes for large means
  y[mm > 6 & x>6 & (2*mm-x)>6 ] <- 0  

  return(y)
}

Then, add the addtional means:

normCurves <-
  data.frame(x = rep(seq(-6, 18, length.out = 1001)
                     , times = 7)
             , myMean = rep(seq(0, 12, 2)
                            , each = 1001))

Now, if you re-run the same plotting code as above, you get:

You now also mention you want to animate, which is relatively straightforward:

toAnimate <-
  ggplot(normCurves
         , aes(x
               , group = myMean
               , frame = myMean)) +
  geom_line(aes(y = density)) +
  geom_area(aes(y = toHighlight)
            , position = "identity") +
  ggtitle("Mean = ")

gganimate::gg_animate(toAnimate)

gives