using stat_function and facet_wrap together in ggplot2 in R

Viewed 17086

I am trying to plot lattice type data with ggplot2 and then superimpose a normal distribution over the sample data to illustrate how far off normal the underlying data is. I would like to have the normal dist on top to have the same mean and stdev as the panel.

here's an example:

library(ggplot2)

#make some example data
dd<-data.frame(matrix(rnorm(144, mean=2, sd=2),72,2),c(rep("A",24),rep("B",24),rep("C",24)))
colnames(dd) <- c("x_value", "Predicted_value",  "State_CD")

#This works
pg <- ggplot(dd) + geom_density(aes(x=Predicted_value)) +  facet_wrap(~State_CD)
print(pg)

That all works great and produces a nice three panel graph of the data. How do I add the normal dist on top? It seems I would use stat_function, but this fails:

#this fails
pg <- ggplot(dd) + geom_density(aes(x=Predicted_value)) + stat_function(fun=dnorm) +  facet_wrap(~State_CD)
print(pg)

It appears that the stat_function is not getting along with the facet_wrap feature. How do I get these two to play nicely?

------------EDIT---------

I tried to integrate ideas from two of the answers below and I am still not there:

using a combination of both answers I can hack together this:

library(ggplot)
library(plyr)

#make some example data
dd<-data.frame(matrix(rnorm(108, mean=2, sd=2),36,2),c(rep("A",24),rep("B",24),rep("C",24)))
colnames(dd) <- c("x_value", "Predicted_value",  "State_CD")

DevMeanSt <- ddply(dd, c("State_CD"), function(df)mean(df$Predicted_value)) 
colnames(DevMeanSt) <- c("State_CD", "mean")
DevSdSt <- ddply(dd, c("State_CD"), function(df)sd(df$Predicted_value) )
colnames(DevSdSt) <- c("State_CD", "sd")
DevStatsSt <- merge(DevMeanSt, DevSdSt)

pg <- ggplot(dd, aes(x=Predicted_value))
pg <- pg + geom_density()
pg <- pg + stat_function(fun=dnorm, colour='red', args=list(mean=DevStatsSt$mean, sd=DevStatsSt$sd))
pg <- pg + facet_wrap(~State_CD)
print(pg)

which is really close... except something is wrong with the normal dist plotting:

enter image description here

what am I doing wrong here?

6 Answers

Orginally posted as an answer to this question, I was encouraged to share my solution here too.

I too became frustrated with overlaying theoretical densities over empirical data, so I wrote a function that automated this process. Since 2009 when this question was first posed, ggplot2 has greatly expanded the extensibility, so I've put it in a extension package on github (EDIT: you can find it on CRAN now).

library(ggplot2)
library(ggh4x)

set.seed(0)

# Make the example data
dd <- data.frame(matrix(rnorm(144, mean=2, sd=2),72,2),
                 c(rep("A",24),rep("B",24),rep("C",24)))
colnames(dd) <- c("x_value", "Predicted_value",  "State_CD")

ggplot(dd, aes(Predicted_value)) +
  geom_density() +
  stat_theodensity(colour = "red") +
  facet_wrap(~ State_CD)

Created on 2021-01-28 by the reprex package (v0.3.0)

If you are willing to use ggformula, then this is pretty easy. (It is also possible to mix and match and use ggformula just for the distribution overlay, but I'll illustrate the full on ggformula approach.)

library(ggformula)
theme_set(theme_bw())

gf_dens( ~ Sepal.Length | Species, data = iris) %>%
  gf_fitdistr(color = "red") %>% 
  gf_fitdistr(dist = "gamma", color = "blue")

Created on 2019-01-15 by the reprex package (v0.2.1)

Related