ggplot: fill area between confidence interval area and 0

user_4198688 0
Viewed 200

I would like to nicely fill the area between the bounds of confidence intervals (between the upper or lower bound of an area) and 0 using different colours for positive or negative (below or above 0).

I managed to get close to what I expect, but there are some artefacts at the bounds of each area. Any idea on how to address these issues appearing at the edges? Thanks!

library(tidyverse)
library(modelbased)

data <- data.frame(x = rnorm(100, 0, 3))
data$y <- cos(data$x) + rnorm(100, 0, 0.5)

newdata <- lm(y ~ poly(x, 6), data=data) %>% 
  modelbased::estimate_link(length=50) 

head(newdata)
#>           x Predicted    CI_low  CI_high
#> 1 -6.387080 0.9694282 0.1751638 1.763693
#> 2 -6.114378 1.3088500 0.7334309 1.884269
#> 3 -5.841676 1.3578983 0.8677253 1.848071
#> 4 -5.568974 1.2084960 0.7437400 1.673252
#> 5 -5.296272 0.9360337 0.4889862 1.383081
#> 6 -5.023570 0.6010678 0.1821074 1.020028

newdata <- newdata %>% 
  rowwise() %>% 
  mutate(Limit = min(abs(CI_low), abs(CI_high)),
         Positive = ifelse(CI_low > 0, Limit, NA),
         Negative = ifelse(CI_low < 0, -Limit, NA)) 
  
  
newdata %>% 
  ggplot(aes(x=x, y=Predicted)) +
  geom_ribbon(aes(ymin = 0, ymax = Positive), fill="green") +
  geom_ribbon(aes(ymin = Negative, ymax = 0), fill="red") +
  geom_ribbon(aes(ymin=CI_low, ymax=CI_high), alpha=0.2) +
  geom_line() +
  geom_hline(yintercept=0)

Created on 2020-07-03 by the reprex package (v0.3.0)

dplyr ggplot2 plot r
1 Answers

I think that those artefacts are caused by the definitions of Positive and Negative variables. More precisely, I think you should specify Negative = ifelse(CI_high < 0, -Limit, NA) instead of Negative = ifelse(CI_low < 0, -Limit, NA). In that case, you obtain:

# packages
library(tidyverse)
library(modelbased)

# data
set.seed(1)
x <- rnorm(100, 0, 3)
data <- data.frame(
  x = x, 
  y = cos(x) + rnorm(100, 0, 0.5)
)

# model
newdata <- lm(y ~ poly(x, 6), data = data) %>% 
  modelbased::estimate_link(length = 50) %>% 
  mutate(
    Limit = pmin(abs(CI_low), abs(CI_high)), 
    Positive = ifelse(CI_low  > 0,  Limit, NA), 
    Negative = ifelse(CI_high < 0, -Limit, NA)
  )

# plot
ggplot(newdata, aes(x = x, y = Predicted)) + 
  geom_line() + 
  geom_ribbon(aes(ymin = CI_low, ymax = CI_high), alpha = 0.2) + 
  geom_ribbon(aes(ymin = 0, ymax = Positive), fill = "green") +  
  geom_ribbon(aes(ymin = Negative, ymax = 0), fill = "red")

Created on 2020-07-03 by the reprex package (v0.3.0)

If you want to further improve the CI at the boundaries, you should increase the length of the grid used, for example:

# packages
library(tidyverse)
library(modelbased)

# data
set.seed(1)
x <- rnorm(100, 0, 3)
data <- data.frame(
  x = x, 
  y = cos(x) + rnorm(100, 0, 0.5)
)

# model
newdata <- lm(y ~ poly(x, 6), data = data) %>% 
  modelbased::estimate_link(length = 150) %>% 
  mutate(
    Limit = pmin(abs(CI_low), abs(CI_high)), 
    Positive = ifelse(CI_low  > 0,  Limit, NA), 
    Negative = ifelse(CI_high < 0, -Limit, NA)
  )

# plot
ggplot(newdata, aes(x = x, y = Predicted)) + 
  geom_line() + 
  geom_ribbon(aes(ymin = CI_low, ymax = CI_high), alpha = 0.2) + 
  geom_ribbon(aes(ymin = 0, ymax = Positive), fill = "green") +  
  geom_ribbon(aes(ymin = Negative, ymax = 0), fill = "red")

Created on 2020-07-03 by the reprex package (v0.3.0)

user_102804110
Related