How to better plot and compare overlapping histograms?

I want to compare the distribution of 3 different time spans:

So I plot the histograms together, along with the model curve.

But I'm afraid that overlapping histograms makes each hard to see. The histogram are also set to half-transparent so the overlapping could be seen. But it also makes the color overlap, making it hard to discern one from another.

An additional problem is that, I also want to compare the bootstrapped result of the histogram, so I re-sample and plot for a lot of times:

I believe this is also very hard to see.

I'm wondering, what would be a good way for plotting this comparison? What can I do to make the plot more discernible?

• I think you want this, but you get this. Commented May 7, 2016 at 14:49
• @AntoniParellada, 3d plot may be a way, and I've seen this kind plot before. But in my case, I don't care about the addtional axis and the comparision along that axis, so I didn't use it. Commented May 9, 2016 at 1:27
• OK Nice question in any case. Commented May 9, 2016 at 1:44

The usua alternatives to display "overlapping" histograms are to:

• place the bar side by side (but I don't think that it is working well visually in most of the situations):

• connect the heights of the bars with a line (and drop the bar itself - there exists alternatives where the outline of the histogram is plotted, like a skyline):

I am adding R code used to make the figures:

dataf <- bind_rows(lapply(1:10,
function(x) {
data.frame(grp=x,
value=rnorm(100,
mean=runif(1)))
}))

ggplot(dataf) +
geom_histogram(aes(x=value, fill=factor(grp)),
position="dodge", binwidth=.5)

ggplot(dataf) +
geom_freqpoly(aes(x=value, color=factor(grp)), binwidth=.5)

• One could take the second plot a step farther and do a density plot instead of using the heights of the histogram bars, which depend heavily on the choices of cutoffs. I often find a density plot, which is essentially a smoothed histogram, to be more informative about the underlying data distributions. For histograms you're making an arbitrary choice about the cutoff values for the bars,for density plots an arbitrary choice about the kernel for smoothing. Don't really see a good reason to prefer histograms.
– EdM
Commented May 7, 2016 at 14:56
• @EdM : I think that the question author already has density estimates. Beside that histograms may have their merit when low number of values or when the binning has physical meaning. Commented May 7, 2016 at 15:04

Plotting histograms together can be fine, but it breaks down when you have more than two histograms, or the more they overlap, both of which apply in your case. I would suggest you start by making a plot matrix (so long as you don't have so many groups the plots become unusable).

Likewise, plots with too many, and too different, objects can become difficult to interpret. You want to compare histograms, and you want to compare kernel density plots. Note that a plot matrix has a main diagonal for each group, and then the upper and lower triangles are symmetrical. For a given plot in the upper triangle that compares two groups, there is a corresponding plot in the lower triangle that compares the same two groups. Thus, I would suggest you compare histograms in the upper triangle plots, and kernel density plots in the lower triangle plots.

Because it might still be difficult to compare two overlapping histograms in the subplots, I would suggest you make back to back histograms instead of overlapping histograms.

Putting these suggestions together, you could get something like this:

This was coded using R. The double histogram code was adapted from here. I suspect the code won't be interpretable to people who aren't already very familiar with R, but for those who do want to see it, it is displayed below:

data(mtcars)
d  = mtcars[,c("qsec","cyl")]
ud = unstack(d)
ud = data.frame(four  = c(ud[[1]], rep(NA,3)),
six   = c(ud[[2]], rep(NA,7)),
eight =   ud[[3]]             )

upper = function(x, y){
usr = par("usr"); on.exit(par(usr)); par(usr = c(0, 1, 0, 1), new=TRUE)
hx        = hist(x, plot=FALSE)
hy        = hist(y, plot=FALSE)
lim       = ifelse(max(hy$counts)>max(hx$counts), max(hy$counts), max(hx$counts))
hy$counts = - hy$counts
plot(hy, ylim=c(-lim, lim), col="red", xlim=c(14,23), axes=FALSE, main="")
lines(hx, col="blue")
abline(h=0)
}
diag.hist = function(x, ...){
usr = par("usr"); on.exit(par(usr)); par(usr=c(usr[1:2], 0, 1.5), new=TRUE)
hist(x, freq=FALSE, xlim=c(14,23), ylim=c(0,.8), main="", axes=FALSE)
lines(density(na.omit(x)))
}
lower = function(x, y){
usr = par("usr"); on.exit(par(usr)); par(usr = c(0, 1, 0, 1), new=TRUE)
plot( density(na.omit(x)), xlim=c(14,23),ylim=c(0,.5),main="",axes=FALSE, col="blue")
lines(density(na.omit(y)), col="red")
}

windows()
pairs(ud, upper.panel=upper, diag.panel=diag.hist, lower.panel=lower)

• This would be fine if you want to compare the histogram one by one. But in my case, I only want to know the variability in this the 3 overlapping histograms (year 2000, 2005, 2010). So for me, overlapping them would be just enough. Commented May 9, 2016 at 1:30