Plotting boxplots in one graph of rasterlayers of different extents

Question

I have 4 rasterstacks of different zones containing each 12 layers (one for every month) containing data on atmospheric NO2 concentration. I made a layer for each rasterstack containing the mean values for the 12 layers. Now I created a boxplot for each of these mean layers but now I want to combine these boxplots in one plot to easily compare.

I tried to stack the different mean layers but that was not possible because they have a different extent. Here is the code

#import rasters
NS <- stack("north_seca_2019_stack.tiff")
BY <- stack("Biscay_2019_stack.tiff")
BP <- stack("BQPZJR_2019_stack.tiff")
EC <- stack("EngChan_2019_stack.tiff")


# mean of year

meanIgnoringZeroes <- function(x) {
  mean(x[x>0],na.rm=T)
}

meanBY<- overlay(BY,fun=meanIgnoringZeroes)
meanNS<- overlay(NS,fun=meanIgnoringZeroes)
meanBP<- overlay(BP,fun=meanIgnoringZeroes)
meanEC<- overlay(EC,fun=meanIgnoringZeroes)

t<- stack(meanBY,meanNS,meanBP,meanEC)

> t<- stack(meanBY,meanNS)
Error in compareRaster(x) : different extent

B<-boxplot(meanBY)
N<-boxplot(meanNS)
BZ<- boxplot(meanBP)
E<-boxplot(meanEC)

This creates boxplots for each zone but I can't find a way to plot them all together in one plot.