In ggplot, how to fill area between two normal curves

Question

I have two normal curves and I want to fill the right area between both curves, so left curve is inferior y limit and right curve is superior y limit. To plot the curves I am using stat_function() so ggplot draws the curve without defining an y-column in aes(). I have drawn the fill area between the curve and the X axis, but I need the area between both curves and the trick of emptying the left curve with NA doesn't seem to work as I expected.

The code to generate the plot is in a function as I need to plot several different couples of normal curves.

How can I do that?

library(ggplot2)
library(ggthemes)

graf_normal <- function(Xmedia1, Xdt1, Xmedia2, Xdt2) {

  Xmin1 <- Xmedia1-4*Xdt1
  Xmax1 <- Xmedia1+4*Xdt1

  Xmin2 <- Xmedia2-4*Xdt2
  Xmax2 <- Xmedia2+4*Xdt2

  Ymax1 <- max(dnorm(Xmedia1, Xmedia1, Xdt1))
  Ymax2 <- max(dnorm(Xmedia2, Xmedia2, Xdt2))

  Xmin <- min(Xmin1, Xmin2)
  Xmax <- max(Xmax1, Xmax2)

  ggplot(data.frame(X = c(Xmin, Xmax)), aes(x = X)) +
  geom_hline(yintercept = 0, colour = "grey", linewidth = 1) +
  stat_function(fun = dnorm, 
              args = c(Xmedia1, Xdt1), 
              linewidth = 1, 
              colour = "grey") +
  stat_function(fun = dnorm, 
              args = c(Xmedia2, Xdt2), 
              linewidth = 1, 
              colour = "black") +
  geom_segment(aes(x = Xmedia1, y = 0, xend = Xmedia1, yend = Ymax1), 
             linetype = "dashed", 
             linewidth = 0, 
             colour = "grey") +
  geom_segment(aes(x = Xmedia2, y = 0, xend = Xmedia2, yend = Ymax2), 
             linetype = "dashed", 
             linewidth = 0, 
             colour = "black") +
  ####################################################################
  stat_function(fun = dnorm,
            args = c(Xmedia2, Xdt2),
            xlim = c(Xmedia2+1.5*Xdt2,Xmax2),
            geom = "area",
            fill = "red",
            alpha = 0.5) +
  stat_function(fun = dnorm,
              args = c(Xmedia1, Xdt1),
              xlim = c(Xmedia1,Xmax1),
              geom = "area",
              fill = NA,
              alpha = 0.01) +

  ##################################################################

  theme(
    line = element_blank(),
    axis.line.y = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis. Ticks = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    legend. Position = "none",
    panel. Grid = element_blank(),
    panel. Background = element_rect(fill = "lightgray", colour = NA),
 ) +
 xlim(c(Xmin, Xmax)) 
}

g1 <- graf_normal(250, 7, 253, 7)

g1

The plot of both curves I get is this:

Thanks,

EDIT: Using @stephan's code and playing with data filtering, I've been able to do this, easier using geom_ribbon():

Cool way of differencing overlapping zones! Complete code:

library(ggplot2)

graf_normal <- function(Xmedia1, Xdt1, Xmedia2, Xdt2, n = 1000) {
  x1 <- Xmedia1 + 4 * Xdt1 * seq(-1, 1, length. Out = n)
  x2 <- Xmedia2 + 4 * Xdt2 * seq(-1, 1, length. Out = n)

  dat <- data. Frame(
    x = union(x1, x2)
  )
  dat$y1 <- dnorm(dat$x, Xmedia1, Xdt1)
  dat$y2 <- dnorm(dat$x, Xmedia2, Xdt2)

  Ymax1 <- dnorm(Xmedia1, Xmedia1, Xdt1)
  Ymax2 <- dnorm(Xmedia2, Xmedia2, Xdt2)

  ggplot(dat, aes(x)) +
    geom_hline(yintercept = 0, colour = "grey", linewidth = 1) +
    geom_ribbon(
      data =  subset(dat, x >= Xmedia2 + 1.5 * Xdt2),
      aes(ymin = y1, ymax = y2),
      fill = "red", alpha = 0.8
    ) +
    geom_ribbon(
      data =  subset(dat, (x <= Xmedia2 + 1.5 * Xdt2) & (y2 > y1)),
      aes(ymin = y1, ymax = y2),
      fill = "red", alpha = 0.2
    ) +
    geom_ribbon(
      data = subset(dat, x <= Xmedia1 - 1.5 * Xdt2),
      aes(ymin = y1, ymax = y2),
      fill = "blue", alpha = 0.8
    ) +
    geom_ribbon(
      data = subset(dat, (x <= Xmedia2 ) & (y1 > y2)),
      aes(ymin = y1, ymax = y2),
      fill = "blue", alpha = 0.2
    ) +
    annotate(
      geom = "segment",
      x = c(Xmedia1, Xmedia2), y = 0,
      xend = c(Xmedia1, Xmedia2), yend = c(Ymax1, Ymax2),
      linetype = "dashed",
      linewidth = 1,
      colour = c("grey", "black")
    ) +
    geom_line(aes(y = y1), linewidth = 1, colour = "grey") +
    geom_line(aes(y = y2), linewidth = 1, colour = "black") +
    theme(
      line = element_blank(),
      axis.line.y = element_blank(),
      axis. Text = element_blank(),
      axis. Ticks = element_blank(),
      axis. Title = element_blank(),
      legend. Position = "none",
      panel. Grid = element_blank(),
      panel. Background = element_rect(fill = "lightgray", colour = NA),
    )
}

graf_normal(250, 7, 253, 7)

However, the code doesn't work for all curves, working on it!:

graf_normal(250, 7, 253, 3)

Answer 1

One option to fill the area between the normal curves would be to use ggh4x::stat_difference which however requires to compute the values for the densities manually and drawing via geom_line instead of relying on stat_function():

library(ggplot2)
library(ggh4x)
graf_normal <- function(Xmedia1, Xdt1, Xmedia2, Xdt2, n = 101) {
  x1 <- Xmedia1 + 4 * Xdt1 * seq(-1, 1, length.out = n)
  x2 <- Xmedia2 + 4 * Xdt2 * seq(-1, 1, length.out = n)
  
  dat <- data.frame(
    x = union(x1, x2)
  )
  dat$y1 <- dnorm(dat$x, Xmedia1, Xdt1)
  dat$y2 <- dnorm(dat$x, Xmedia2, Xdt2)

  Ymax1 <- dnorm(Xmedia1, Xmedia1, Xdt1)
  Ymax2 <- dnorm(Xmedia2, Xmedia2, Xdt2)

  ggplot(dat, aes(x)) +
    geom_hline(yintercept = 0, colour = "grey", linewidth = 1) +
    ggh4x::stat_difference(
      data = ~ subset(.x, x >= Xmedia2 + 1.5 * Xdt2),
      aes(ymin = y1, ymax = y2)
    ) +
    annotate(
      geom = "segment",
      x = c(Xmedia1, Xmedia2), y = 0,
      xend = c(Xmedia1, Xmedia2), yend = c(Ymax1, Ymax2),
      linetype = "dashed",
      linewidth = 1,
      colour = c("grey", "black")
    ) +
    geom_line(aes(y = y1), linewidth = 1, colour = "grey") +
    geom_line(aes(y = y2), linewidth = 1, colour = "black") +
    scale_fill_manual(values = c(scales::alpha("red", .5), "transparent")) +
    theme(
      line = element_blank(),
      axis.line.y = element_blank(),
      axis.text = element_blank(),
      axis.ticks = element_blank(),
      axis.title = element_blank(),
      legend.position = "none",
      panel.grid = element_blank(),
      panel.background = element_rect(fill = "lightgray", colour = NA),
    )
}

graf_normal(250, 7, 253, 7)

EDIT Actually stat_differnce is not really needed for your case. Was thinking too complicated. As @JuanRiera mentioned in his comment, we could fill the area using a geom_ribbon:

library(ggplot2)

graf_normal <- function(Xmedia1, Xdt1, Xmedia2, Xdt2, n = 101) {
  x1 <- Xmedia1 + 4 * Xdt1 * seq(-1, 1, length.out = n)
  x2 <- Xmedia2 + 4 * Xdt2 * seq(-1, 1, length.out = n)

  dat <- data.frame(
    x = union(x1, x2)
  )
  dat$y1 <- dnorm(dat$x, Xmedia1, Xdt1)
  dat$y2 <- dnorm(dat$x, Xmedia2, Xdt2)

  Ymax1 <- dnorm(Xmedia1, Xmedia1, Xdt1)
  Ymax2 <- dnorm(Xmedia2, Xmedia2, Xdt2)

  ggplot(dat, aes(x)) +
    geom_hline(yintercept = 0, colour = "grey", linewidth = 1) +
    geom_ribbon(aes(ymin = y1, ymax = y2),
      data = subset(dat, x >= Xmedia2 + 1.5 * Xdt2),
      fill = "red", alpha = 0.5
    ) +
    annotate(
      geom = "segment",
      x = c(Xmedia1, Xmedia2), y = 0,
      xend = c(Xmedia1, Xmedia2), yend = c(Ymax1, Ymax2),
      linetype = "dashed",
      linewidth = 1,
      colour = c("grey", "black")
    ) +
    geom_line(aes(y = y1), linewidth = 1, colour = "grey") +
    geom_line(aes(y = y2), linewidth = 1, colour = "black") +
    theme(
      line = element_blank(),
      axis.line.y = element_blank(),
      axis.text = element_blank(),
      axis.ticks = element_blank(),
      axis.title = element_blank(),
      legend.position = "none",
      panel.grid = element_blank(),
      panel.background = element_rect(fill = "lightgray", colour = NA),
    )
}

graf_normal(250, 7, 253, 7)