Single figure with multiple plots - R ggplot2

So in my efforts to revise some PhD chapters to submit for publication, I've had to figure out how to do some nice graphics in R to reduced, condense and improve the figures in my chapters.

Some simple steps are described below, as well as the example outputs.

I've tested the layout() function in R, as well as facet.grid, but the latter is for when axes on multiple graphs are the same. This is the case with me, and so I've put together each figure independently and then simple added them together and given them the same legends where applicable.

R code as follows:

library("ggplot2")

• Create the first figure called popPlot - how population changes over time

popPlot <- ggplot(modelPop, aes(x=factor(year), y=currentPop)) + geom_boxplot()
popPlot <- popPlot + theme(panel.background=element_rect(fill='white', colour='black')) + ylab("Population Size") + xlab("Year of Simulation")
***** the following overrides the automatic scale and sets the x axis labels *******
popPlot <- popPlot + scale_x_discrete(breaks=c(11, 20, 30, 40, 50, 60, 70, 80, 90, 100), labels=c("11", "20", "30", "40", "50", "60", "70", "80", "90", "100"))
***** leaving the x axis title as blank() means it doesn't get printed on this figure ******
popPlot <- popPlot + theme(axis.title.x=element_blank(), axis.title.y=element_text(size=30))
popPlot <- popPlot + theme(axis.text.x=element_text(size=25), axis.text.y=element_text(size=25))
****** the following adds the line indicating the mean value for each year of the simulation *****
popPlot <- popPlot + stat_summary(fun.y=mean, geom="line", aes(group=1), colour='red', size=2)
popPlot

• Create the second figure called popPlot.L - how population changes over a longer time period

popPlot.L <- ggplot(longModelPop, aes(x=factor(year), y=currentPop)) + geom_boxplot()
******* if the axis titles are not set out here, the data column names will be applied automatically ****
popPlot.L <- popPlot.L + theme(panel.background=element_rect(fill='white', colour='black')) + ylab("Population Size") + xlab("Year of Simulation")
popPlot.L <- popPlot.L + scale_x_discrete(breaks=c(11, 100, 200, 300, 400, 500), labels=c("11", "100", "200", "300", "400", "500"))
****note both axis titles are removed in this figure ******
popPlot.L <- popPlot.L + theme(axis.title.x=element_blank(), axis.title.y=element_blank())
popPlot.L <- popPlot.L + theme(axis.text.x=element_text(size=25), axis.text.y=element_text(size=25))
popPlot.L <- popPlot.L + stat_summary(fun.y=mean, geom="line", aes(group=1), colour='red', size=2)
popPlot.L

• Now plot both figures together

library("grid")
library("gridExtra")

******* this a single line codes for the plotting of both graphs in the same figure *******
grid.arrange(popPlot, popPlot.L, ncol=2, sub = textGrob("Year of simulation", gp=gpar(fontsize=30)))

This produces the figure below: