Creating Multiplots with NADA / NADA2 Package in R

Question

I am trying to create multipannel plots using the NADA2 package. I would like to end with the three plots side-by-side in one figure. The benefit of the NADA2 package is that it creates a maximum liklihood estimation for values below the detection limit and runs statistics on and plots the estimations. This dataset is half nondetects, so I really prefer to keep using the NADA2 plots and simply arrange the layout. Any help would be appreciated.

Example (abbreviated) data and code are listed below. Please note that:

• Data is originally in one dataframe but has been split because I had difficulty calling an individual "name" factors to create individual plots. I am open to solutions using one or multiple dataframes.
• Errors received onsample dataset are slightly different than those received on real dataset, so I have included both when appropriate.
• Real dataset consists of ~22850 observations with 14 "names", 12 "years", and 24 parameters (e.g., "mc") to be plotted

Previous Attempts: Multipannel Plots: ggarange and layout.show The NADA2 package does not appear compatible with ggarrange

Error in gList(...) : only 'grobs' allowed in "gList" and I have been having difficulty getting the plots to appear using the layout.show function (plots show up with correct dimensions individually but not in the layout).

Groups within One Plot: xgroups argument

NADA plots include an xgroups argument, but I receive an error Error in levels<-(*tmp*, value = as.character(levels)) : factor level [2] is duplicated. When I attempt to remedy the levels error, it does not appear to work, as I receive the same error but with factor level [4]. On the test dataset I receive a different error Error in if (sum(1 - as.integer(nd[y.grp == levels(grp)[i]])) <= 2) { : missing value where TRUE/FALSE needed)

# Setup -----
## load packages -----
library(NADA2)
## create a data frames -----
df1 <- data.frame(
    name = c('Ja'), 
    year = c('2008', '2008', '2008', '2008', '2008', '2008', '2008', '2008', '2009', '2009', '2009', '2009', '2009', '2009', '2009'),
    mc = c(0.075, 0.075, 0.075, 3.544, 0.559, 0.508, 0.469, 0.46, 0.075, 0.075, 0.075, 3.402, 0.294, 0.293, 0.075),
    mc_cens = c(TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE), stringsAsFactors = FALSE)

df2 <- data.frame(
    name = c('Os'), 
    year = c('2008', '2008', '2008', '2008', '2008', '2008', '2008', '2008', '2009', '2009', '2009', '2009', '2009', '2009', '2009'),
    mc = c(0.075, 0.075, 0.075, 0.063, 1.505, 1.818, 1.312, 0.9, 0.075, 0.075, 0.075, 0.823, 0.807, 1.818, 1.312),
    mc_cens = c(TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE), stringsAsFactors = FALSE)

df3 <- data.frame(
    name = c('Tip'),
    year = c('2008', '2008', '2008', '2008', '2008', '2008', '2008', '2008', '2009', '2009', '2009', '2009', '2009', '2009', '2009'),
    mc = c(0.075, 0.075, 0.075, 4.515, 1.321, 0.82, 0.687, 0.687, 0.075, 0.075, 0.075, 1.56, 1.04, 1.027, 0.82),
    mc_cens = c(TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE), stringsAsFactors = FALSE)

### Change column data types
df1$mc <- sapply(df1$mc, as.numeric)
df2$mc <- sapply(df2$mc, as.numeric)
df3$mc <- sapply(df3$mc, as.numeric)

df1$mc_cens <- as.logical(df1$mc_cens)
df2$mc_cens <- as.logical(df2$mc_cens)
df3$mc_cens <- as.logical(df3$mc_cens)

### view dataframes' structure
str(df1)
str(df2)
str(df3)

# Plotting -----
## Create NADA2 plots -----
plot_j <- cboxplot(df1$mc, df1$mc_cens, df1$year, show=TRUE, LOG = FALSE, minmax = FALSE, Title = "Ja MC Per Year", Ylab = "MC", Xlab = "Year")
plot_o <- cboxplot(df2$mc, df2$mc_cens, df2$year, show=TRUE, LOG = FALSE, minmax = FALSE, Title = "Os MC Per Year", Ylab = "MC", Xlab = "Year")
plot_t <- cboxplot(df3$mc, df3$mc_cens, df3$year, show=TRUE, LOG = FALSE, minmax = FALSE, Title = "Tip MC Per Year", Ylab = "MC", Xlab = "Year")

## Merged Plots -----
### Create Merged Dataframe -----
df_all <- rbind(df1, df2, df3)

## Attempts to make unique factors
df_all$name <- factor(name[1:3], levels = df_all$name[1:3])
within(acc_routine_open_df_plotting,lake_name<- ave(as.character(lake_name),FUN=make.unique))

### Create NADA2 Plots Using xgroup -----
plot_all <- cboxplot(df_all$mc, df_all$mc_cens, df_all$year, xgroup=df_all$name, show=TRUE, LOG = FALSE, minmax = FALSE, Title = "MC Per Year", Ylab = "Microcystin Concentration (ppb)", Xlab = "Year")