Why does stat_summary plot multiple/single lines depending on the variable?

We can pre-compute the median using group_by and mutate, which I often find more reliable and easy to understand in its behavior, and then just use geom_vline. Can't answer on the stat_summary side, but interested to know the answer.

mpg %>%
  group_by(cyl) %>%
  mutate(cty_med = median(cty)) %>%
  ggplot(aes(x=cty))+
  geom_histogram()+
  facet_grid(~cyl)+
  geom_vline(aes(xintercept=cty_med))

If you want to generalize this, you can just create a wrapper function that does your calculation and faceting together.

f <- function(df, fct, var) {
  df %>%
    group_by({{fct}}) %>%
    mutate(med = median({{var}})) %>%
    ggplot(aes(x={{var}}))+
    geom_histogram() +
    facet_grid(cols = vars({{fct}})) +
    geom_vline(aes(xintercept=med))
}

f(mpg, cyl, cty)
f(mpg, cyl, displ)

Demetri, here is the R code that will give you what you need:

library(tidyverse) 

g <- mpg %>% 
     ggplot(aes(x=cty)) + 
     geom_histogram() + 
     stat_summary(aes(x = 0, xintercept = stat(y), y = cty), 
                     fun.y = median, geom = "vline", colour = "red") + 
     facet_grid(~ cyl)


g 

The stat_summary() function is set up to compute a summary (in this case, the median) for the variable specified in its y argument. In contrast, the geom_histogram() function creates a histogram for the variable specified in its x argument. So you have to be careful with how you specify the y argument for the stat_summary() function, as seen in the code above.

Note that you dont't need to use group = cyl in your ggplot() call if you are using facet_grid() or facet_wrap() to produce multiple graphical panels. Grouping and facetting are totally different plotting operations: grouping will show different data groups in the same panel; facetting will show different data groups in different panels.

Addendum 1

To check that the summary statistics were computed correctly for each panel, the command below will come in handy:

ggplot_build(g)$data

Scroll to the bottom of the output produced by this command to find the xintercept values used by R - these should be the medians plotted in the various panels. Alternatively, extract these values directly with:

ggplot_build(g)$data[[2]]

The xintercept values can be compared with independently computed median values of cty for each cyl level to ensure agreement.

Addendum 2

The default choice of binwidth for geom_histogram() needs some attention. You can try something like this to allow variable binwidth choice across your different panels:

theme_set(theme_bw())

g <- mpg %>% 
  ggplot(aes(x=cty)) + 
  geom_histogram(binwidth = function(x) 2 * IQR(x) / (length(x)^(1/3)), 
             fill = "lightblue3", colour = "white") + 
  stat_summary(aes(x = 0, xintercept = stat(y), y = cty), 
        fun.y = median, geom = "vline", colour = "red2") + 
  facet_wrap(~ cyl, scales = "free_x")


g 

See this link for other possibilities of binwidth choice: https://github.com/tidyverse/ggplot2/issues/2312.

I have made another post because I believe my exploration and solution are different enough from the original to justify it, but I can merge if anyone disagrees. So, I believe I have figured out what is the cause of the issue with stat_summary and your current solution.

I believe that stat_summary calculates its summary statistic for each unique value of x, when the x variable takes on integer values.

library(tidyverse)

sapply(mpg, class)
#> manufacturer        model        displ         year          cyl        trans 
#>  "character"  "character"    "numeric"    "integer"    "integer"  "character" 
#>          drv          cty          hwy           fl        class 
#>  "character"    "integer"    "integer"  "character"  "character"

See below the same before when using hwy and cty, even when both explicitly converted to numeric rather than integer vectors.

mpg2 <- mpg %>%
  mutate(hwy = as.numeric(hwy),
         cty = as.numeric(cty))

sapply(mpg2, class)
#> manufacturer        model        displ         year          cyl        trans 
#>  "character"  "character"    "numeric"    "integer"    "integer"  "character" 
#>          drv          cty          hwy           fl        class 
#>  "character"    "numeric"    "numeric"  "character"  "character"

mpg2 %>%
  ggplot(aes(x=hwy, group=cyl))+
  geom_histogram()+
  facet_grid(~cyl)+
  stat_summary(aes(xintercept=stat(x), y=0), fun = median, geom = 'vline')

And example with cty:

mpg2 %>%
  ggplot(aes(x=cty, group=cyl))+
  geom_histogram()+
  facet_grid(~cyl)+
  stat_summary(aes(xintercept=stat(x), y=0), fun = median, geom = 'vline')

However, if we make a slight adjustment to cty prior to plotting, adding a minute decimal point, we get the desired behavior.

mpg %>%
  mutate(cty = cty + .000001) %>%
  ggplot(aes(x=cty, group=cyl))+
  geom_histogram()+
  facet_grid(~cyl)+
  stat_summary(aes(xintercept=stat(x), y=0), fun = median, geom = 'vline')

And we see the same behavior with hwy.

mpg %>%
  mutate(hwy = hwy + .000001) %>%
  ggplot(aes(x=hwy, group=cyl))+
  geom_histogram()+
  facet_grid(~cyl)+
  stat_summary(aes(xintercept=stat(x), y=0), fun = median, geom = 'vline')

Of course, this isn't necessarily a desirable solution. Since we are mapping vertical lines, we can instead create a new aes where we instead plot our xintercept as a function of y, and provide a single dummy variable to x within our data range. This then tricks the system into plotting only one median from our single x value, and gives us the desired graph.

mpg %>%
  ggplot(aes(x=cty, group = cyl))+
  geom_histogram()+
  facet_grid(~cyl)+
  stat_summary(aes(x = 3, y = cty, xintercept = stat(y)), fun = median, geom = 'vline')

And there we go! Quite convoluted, and don't really like it as a solution, but I believe this is the way you have to go if using stat_summary.