Boxplot - Displaying the same data, by two groups in one box plot

Question

I have a set of data that is collected from seven different sites, across two different regions. I want to use a box plot to show the differences between the sites, and between regions. I currently have two separate box plot figures. I am trying to work out (but having no luck) how to display both, in one figure... i.e. I would like my figure to display results as: Site 1, Site 2, Site 3, Site 4.... Region 1, Region 2 <- all on one figure. I have supplied the example data and my current code below. You will notice that my axis are flipped - this is intentional as I have longish site names.

Thanks!

#Example script

example<-structure(list(Site = c("NS", "NS", "NS", "NE", "NE", "NE", "KO", 
                                 "KO", "KO", "MI", "MI", "MI", "KI", "KI", "KI", "HP", "HP", "HP", 
                                 "HA", "HA", "HA"), Var1 = c(30, 5, 4, 0.3, 0, 0, 48, 79, 43, 
                                                             30, 30, 50, 5, 3, 3, 20, 31, 50, 30, 40, 68), Region = c("A", 
                                                                                                                      "A", "A", "A", "A", "A", "A", "A", "A", "B", "B", "B", "B", "B", 
                                                                                                                      "B", "B", "B", "B", "B", "B", "B")), class = "data.frame", row.names = c(NA, 
                                                                                                                                                                                               -21L))
######## By Site
library(forcats)
p <- example %>%
  mutate(Site = fct_relevel(Site, "HA", "HP",
                            "KI", "MI","KO", "NE",
                            "NS"
  )) %>%
  ggplot( aes(x=Site, y=Var1, fill = Region)) +
  geom_boxplot()+xlab("Site") + ylab ("Variable 1") + ylim(0,100) +
  coord_flip() + theme_bw() 
plot1<- p+scale_fill_manual(values=c("#0072B2", "#D55E00")) + theme(axis.text=element_text(size=11),
                                                                    axis.title=element_text(size=12))

plot1
####################### By Regions

example$Region<-factor(example$Region, levels=c("B", "A"))


pplot2<-ggplot(example, aes(x = Region, y = Var1, fill = Region)) + geom_boxplot()
pplot2<-pplot2 + ylab("Variable 1") + ylim (0,100) + 
  scale_x_discrete(limits = rev(levels("Region"))) + coord_flip()+ theme_bw() + scale_fill_manual(values=c("#D55E00", "#0072B2")) + theme(axis.text=element_text(size=11),
                                                                                                                                          axis.title=element_text(size=12))

pplot2

Answer 1

If you want to only stack the plots, you can try library(patchwork). It has a simple syntax for combining multiple plots together. Here is a tutorial.

library(patchwork)
p3 <- plot1/(pplot2 + theme(legend.position="none"))

If you actually want them to share the same axes, I have a sort of hacky/clunky solution. If you stack the dataframe so that "Region" is also included in "Site" column then you can put Site and Region as the same axis. Downside is that this makes labeling the axis a little clunkier. So, I toyed with hjust and adding spaces into the lab(x=) to get it to line up right. I also added geom_vline() (since coordinates are flipped) to get visually separate out where the region-based and site-based boxplots are. Again, not elegant, but gets the job done.

example2 <- example %>% mutate(Site = Region)
dat <- rbind(example, example2)

p4 <- dat %>% mutate(Site = fct_relevel(Site, "HA", "HP",
                                  "KI", "MI","KO", "NE",
                                  "NS", "B", "A"
)) %>% ggplot(aes(x=Site, y=Var1, fill=Region)) + geom_boxplot() + 
  coord_flip() + theme_bw() + scale_fill_manual(values=c("#0072B2", "#D55E00")) + 
  theme(axis.text=element_text(size=11),
        axis.title=element_text(size=12)) + labs(y="Variable 1", 
                                                 x="Site                            Region") + 
  theme(axis.title.y=element_text(hjust=.90)) +
  geom_vline(xintercept = 7.5)