Q:

R - Histograms with shared/same x and y axes

I want to plot two histograms where the x and y ranges are the same for both. After reading some posts, my solution is to use ggplot2, geom_histogram twice. The first time I am creating the plots without plotting for each dataset of interest with the intention to get the maximum y/count and x axes values among all plots of interest. For example, having two plots, if for the first one ymax_1 = 10 for the other ymax_2 = 15, then both plots will have a y axis range from 0 to 15 at least. Similarly holds for the x axis.
After this plot, I take the ymax/xmax values and plot the histograms like before with the addition of xlim( 0, xmax) and ylim( 0, ymax). However, when I do this, the amount of counts changes. More specifically, in the first plots where I don't have any xlim/ylim specified I get from ggplot_build( ggplot(...) + geom_histogram(...)) ymax = 2000 but when I use xlim the second time I get ymax = 4000. Nevertheless, from the 1st plot I have ymax = 2000 and hence the second time the histograms are not plotted properly. When I remove the xlim option I get the same result.
How and why the xlim option affects the amount of counts? I hope that was clear.
df = read.table( paste( path, f, sep = "/"), header = TRUE, fill = TRUE, sep = ",", stringsAsFactors = TRUE)
measure = colnames( df)[ 7]
combs = unique( df[, c( 'A', 'B', 'C')])
# order combs in specific order to get a specific sequence of plots
combs = combs[ with( combs, order( B, C, A)), ]
bns = lst()
xmxs = lst()
ymxs = lst()
for( j in seq( 1, length( combs[ , 1]), 2)) {
 if( combs[ j, 2] == combs[ j, 3]) {
  next
 }
 tmp = subset( df, A == combs[ j, 1] & B == combs[ j, 2] & C == combs[ j, 3], select = c( measure))
 # Freedman – Diaconis rule, "On the histogram as a density estimator: L2 theory"
 bw = 2 * IQR( tmp[ , 1]) / ( length( tmp[ , 1])^(1/3))
 bns[[ j]] = ceiling( ( max( tmp[ , 1]) - min( tmp[ , 1])) / bw)

 plots[[ j]] = ggplot( tmp, aes_string( measure)) + geom_histogram( bins = bns[[ j]], aes( fill = ..count..))
 histg = ggplot_build( plots[[ j]])$data[[ 1]]
 ymxs[[ j]] = max( histg$count)
 xmxs[[ j]] = max( histg$x)

 tmp = subset( df, A == combs[ j + 1, 1] & B == combs[ j + 1, 2] & C == combs[ j + 1, 3], select = c( measure))
 # Freedman – Diaconis rule, "On the histogram as a density estimator: L2 theory"
 bw = 2 * IQR( tmp[ , 1]) / ( length( tmp[ , 1])^(1/3))
 bns[[ j + 1]] = ceiling( ( max( tmp[ , 1]) - min( tmp[ , 1])) / bw)

 plots[[ j + 1]] = ggplot( tmp, aes_string( measure)) + geom_histogram( bins = bns[[ j + 1]], aes( fill = ..count..))
 histg = ggplot_build( plots[[ j + 1]])$data[[ 1]]
 ymxs[[ j + 1]] = max( histg$count)
 xmxs[[ j + 1]] = max( histg$x)
 if( ymxs[[ j]] > ymxs[[ j + 1]]) {
  ymxs[[ j + 1]] = ymxs[[ j]]
 }
 else {
  ymxs[[ j]] = ymxs[[ j + 1]]
 }
 if( xmxs[[ j]] > xmxs[[ j + 1]]) {
  xmxs[[ j + 1]] = xmxs[[ j]]
 }
 else {
  xmxs[[ j]] = xmxs[[ j + 1]]
 }
}
pplots = lst()
for( j in 1 : length( combs[ , 1])) {
 if( combs[ j, 2] == combs[ j, 3]) {
  next
 }
 tmp = subset( df, A == combs[ j, 1] & B == combs[ j, 2] & C == combs[ j, 3], select = c( measure))
 avg = sprintf( "%.2f", mean( tmp[ , 1]))
 stdv = sprintf( "%.2f", std( tmp[ , 1]))
 count = length( tmp[ , 1])
 entities[[ j]] = paste( combs[ j, 1], " ", combs[ j, 2], " vs ", combs[ j, 3])
pplots[[ j]] = ggplot( tmp, aes_string( measure)) +
 geom_histogram( bins = bns[[ j]], aes( fill = ..count..)) +
 # xlim( 0, 1.2*xmxs[[ j]]) +
 # ylim( 0, 1.2*ymxs[[ j]]) +
 ggtitle( bquote( atop( paste( .(entities[[ j]])), paste( mu, " = ", .( avg), ", ", sigma, " = ", .( stdv), ", #cells = ", .( count), sep = " ")))) +
 theme( plot.title = element_text( size = 20), axis.text = element_text( size = 12), axis.title = element_text( size = 15))
 }

 # plot every two plots because the Reference.Population is the same
 for( j in seq( 1, length( plots), 2)) {
 fileext = str_remove_all( entities[[ j]], 'N')
 filename_hi = paste( gsub( '.{4}$', '', f), "_distribution_", fileext, ".png", sep = "")
 png( filename = paste( path, filename_hi, sep = "/"))
 grid.draw( rbind( ggplotGrob( pplots[[ j]]), ggplotGrob( pplots[[ j + 1]]), size = "last"))
 dev.off()
 }

So, in the code above, plots contains the initial plots from which I get the min and max values for y,x axes and pplots contains the plots that I plot finally using the xlim/ylim options. However, for example,
max( plots[[ 8]]$data[[ 1]]$count) != max( plots[[ 8]]$data[[ 1]]$count)
when I use the xlim option. The first one gives 1947 and the other one gives 4529 for my data.
Thanks

A:

As an alternative to the other posts you read, I suggest combining the datasets into one, and faceting them. To do that, you need to select the column you want histogrammed, and add a column that indicates the dataset from which the data is extracted.
For this example, I'll combine iris$Sepal.Length and mtcars$disp.
range(mtcars$disp)
# [1]  71.1 472.0
range(iris$Sepal.Length)
# [1] 4.3 7.9

Since these example data are so different, I'll scale one so that plot looks more comparable ... but different enough so that you can see that the axes are shared.
400 * (range(iris$Sepal.Length) - 4)
# [1]  120 1560

over to you if you need anything like this for your data.
From here, combine the relevant fields:
combined_dat <- rbind(
  cbind.data.frame(src = "iris Sepal.Length", val = 400 * (iris[, c("Sepal.Length")] - 4)),
  cbind.data.frame(src = "mtcars disp*", val = mtcars[, c("disp")])
)

head(combined_dat)
#                 src val
# 1 iris Sepal.Length 440
# 2 iris Sepal.Length 360
# 3 iris Sepal.Length 280
# 4 iris Sepal.Length 240
# 5 iris Sepal.Length 400
# 6 iris Sepal.Length 560

tail(combined_dat)
#              src   val
# 177 mtcars disp* 120.3
# 178 mtcars disp*  95.1
# 179 mtcars disp* 351.0
# 180 mtcars disp* 145.0
# 181 mtcars disp* 301.0
# 182 mtcars disp* 121.0

And then plot.
ggplot(combined_dat, aes(val)) +
  geom_histogram() +
  facet_wrap(~ src, ncol = 1)
# `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

