diff --git a/NEWS.md b/NEWS.md
index a01fd00..71a322b 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,5 @@
 # tidyvpc 1.5.0
-* Support for generating percentage blq/alq plots using `plot.tidyvpcobj`. For VPC with `censoring()`, users can supply arguments `censoring.type` (options are `'none'`, `'blq'`, `'alq'`, or `'both'`, defaults to 'none') and `censoring.output` (options are `'grid'` or `'list'`, defaults to `'grid'`).[#21](https://github.com/certara/issues/21)
+* Support for generating percentage blq/alq plots using `plot.tidyvpcobj`. For VPC with `censoring()`, users can supply arguments `censoring.type` (options are `'none'`, `'blq'`, `'alq'`, or `'both'`, defaults to 'none') and `censoring.output` (options are `'grid'` or `'list'`, defaults to `'grid'`).[#21](https://github.com/certara/tidyvpc/issues/21)
 * Plotting updates were made for ggplot2 version 3.4.0 to use `linewidth` instead of `size` for lines[#39](https://github.com/certara/tidyvpc/issues/39).
 * `simulated.tidyvpcobj()` detects if the number of simulated rows is not an integer multiple of the number of observed rows and adds the new `xsim` argument to test that x values match between replicated simulations.  It will suggest that MDV filtering may not have occurred if either of these fails [#35](https://github.com/certara/tidyvpc/issues/35).
 * Prevent division by zero in `predcorrect()` transformation [#31](https://github.com/certara/tidyvpc/issues/31).
diff --git a/R/plot.R b/R/plot.R
index 7ddf450..988231c 100644
--- a/R/plot.R
+++ b/R/plot.R
@@ -178,13 +178,15 @@ plot.tidyvpcobj <- function(x,
     if (censoring.output == "grid") {
       #Return egg
       g <- do.call(egg::ggarrange, grid_list)
+      return(invisible(g))
     } else {
       #Return list
       g <- setdiff(grid_list, grid_args)
+      return(g)
     }
   }
 
-  g
+  return(g)
 }
 
 #' Expand single-value vpc groups to a finite width so that they show up with `geom_ribbon()`
diff --git a/vignettes/tidyvpc_cont.Rmd b/vignettes/tidyvpc_cont.Rmd
index 72e8c5f..39cf8f7 100644
--- a/vignettes/tidyvpc_cont.Rmd
+++ b/vignettes/tidyvpc_cont.Rmd
@@ -229,7 +229,7 @@ vpc <- observed(obs_data, x=TIME, y=DV) %>%
     binning(bin = "jenks", nbins = 5) %>%
     vpcstats()
 
-plot(vpc)
+plot(vpc, censoring.type = "blq")
 ```
 
 Censoring using `LLOQ` variable in the data.
@@ -248,11 +248,11 @@ vpc <- observed(obs_data, x=TIME, y=DV) %>%
     binless(optimize = FALSE, lambda = c(1.5, 2.5, 1.7)) %>%
     vpcstats()
 
-plot(vpc)
+plot(vpc, censoring.type = "blq")
 ```
 
 
-The `tidyvpc` package also allows you to use different LLOQ for each level of stratification variable. We'll set an `LLOQ` value of `50` for `Study A` and `25` for `Study B` and calculate statistics at 5%, 50%, 90% quantiles.
+The `tidyvpc` package also allows you to use different LLOQ for each level of stratification variable. We'll set an `LLOQ` value of `50` for `Study A` and `25` for `Study B` and calculate statistics at 10%, 50%, 90% quantiles.
 
 ```{r  fig.width = 9, fig.height = 6, out.width=640, warning = FALSE}
 obs_data$LLOQ <- obs_data[, ifelse(STUDY == "Study A", 50, 25)]
@@ -264,9 +264,43 @@ vpc <- observed(obs_data, x=TIME, y=DV) %>%
     binning(bin = "pam", nbins = 4) %>%
     vpcstats(qpred = c(0.1, 0.5, 0.9))
 
-plot(vpc)
+plot(vpc, censoring.type = "blq", facet.scales = "fixed")
+```
+
+The `tidyvpc` package also supports usage of `censoring()` with ALQ data, similar to above usage with BLQ data.
+
+```{r  fig.width = 9, fig.height = 6, out.width=640, warning = FALSE}
+obs_data$ULOQ <- obs_data[, ifelse(STUDY == "Study A", 125, 100)]
+
+vpc <- observed(obs_data, x = TIME, y = DV) |>
+  simulated(sim_data, y = DV) |>
+  censoring(alq = DV > ULOQ, uloq = ULOQ) |>
+  stratify(~ STUDY) |>
+  binning(bin = NTIME) |>
+  vpcstats(qpred = c(0.1, 0.5, 0.9))
+
+plot(vpc, censoring.type = "alq")
+
+```
+
+If using `censoring()` with both ALQ and BLQ data, set `censoring.type = "both"` in the `plot()` function to display both percentage of BLQ and ALQ plots as a grid in the resulting VPC plot.
+
+```{r  fig.width = 9, fig.height = 6, out.width=640, warning = FALSE}
+obs_data$LLOQ <- obs_data[, ifelse(STUDY == "Study A", 50, 25)]
+obs_data$ULOQ <- obs_data[, ifelse(STUDY == "Study A", 125, 100)]
+
+vpc <- observed(obs_data, x = TIME, y = DV) |>
+  simulated(sim_data, y = DV) |>
+  censoring(blq = DV < LLOQ, lloq = LLOQ,  alq = DV > ULOQ, uloq = ULOQ) |>
+  stratify(~ STUDY) |>
+  binning(bin = NTIME) |>
+  vpcstats(qpred = c(0.1, 0.5, 0.9))
+
+plot(vpc, censoring.type = "both", facet.scales = "fixed")
 ```
 
+If you want to return the percentage of BLQ and/or ALQ plots individually as elements in a list, instead of arranged in a grid, use the `censoring.output` argument e.g.,  `plot_list <- plot(vpc, censoring.type = "both", censoring.output = "list")`.
+
 ## `predcorrect()`
 
 *Note: The `predcorrect()` function is only applicable for continuous VPC.*
diff --git a/vignettes/tidyvpc_whats_new.Rmd b/vignettes/tidyvpc_whats_new.Rmd
index 1723bb9..0013a20 100644
--- a/vignettes/tidyvpc_whats_new.Rmd
+++ b/vignettes/tidyvpc_whats_new.Rmd
@@ -9,6 +9,56 @@ vignette: >
   \usepackage[utf8]{inputenc}
 ---
 
+```{r, warning = FALSE, echo = FALSE, message = FALSE}
+knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 9, fig.height = 6, out.width = 640, warning = FALSE)
+```
+
+## tidyvpc 1.5.0
+
+### Plot Percentage of BLQ and/or ALQ
+
+For VPC with censoring, users can supply additional arguments to `plot.tidyvpcobj` e.g., `censoring.type` (options are `'none'`, `'blq'`, `'alq'`, or `'both'`, defaults to `'none'`) and `censoring.output` (options are `'grid'` or `'list'`, defaults to `'grid'`).
+
+If `censoring.output = 'grid'`, the plots will be arranged into single grid plot. Users may pass additional arguments via ellipsis to `egg::ggarrange` e.g., `nrow = 1`, `ncol = 2` in order to customize plots in grid arrangement.
+
+If `censoring.output = 'list'`, the resulting plots will be returned individually as elements in list.
+
+Example usage is below:
+
+```{r}
+library(tidyvpc)
+obs_data <- obs_data[MDV == 0]
+sim_data <- sim_data[MDV == 0]
+obs_data$LLOQ <- obs_data[, ifelse(STUDY == "Study A", 50, 25)]
+obs_data$ULOQ <- obs_data[, ifelse(STUDY == "Study A", 125, 100)]
+
+vpc <- observed(obs_data, x = TIME, y = DV) |>
+  simulated(sim_data, y = DV) |>
+  censoring(blq = DV < LLOQ, lloq = LLOQ,  alq = DV > ULOQ, uloq = ULOQ) |>
+  stratify(~ STUDY) |>
+  binning(bin = NTIME) |>
+  vpcstats(qpred = c(0.1, 0.5, 0.9))
+```
+
+If blq data, users may supply `censoring.type = "blq"`:
+
+```{r}
+plot(vpc, censoring.type = "blq", censoring.output = "grid", facet.scales = "fixed")
+```
+
+If alq data, users may supply `censoring.type = "alq"`:
+
+```{r}
+plot(vpc, censoring.type = "alq", censoring.output = "grid", ncol = 2, nrow = 1)
+```
+
+If both blq and alq data, users may supply `censoring.type = "both"`
+
+```{r}
+vpc_plots <- plot(vpc, censoring.type = "both", censoring.output = "list")
+```
+
+By default, when `censoring.tidyvpcobj` is used, no percentage blq/alq plots will be returned e.g., default for `censoring.type = 'none'`. If users specify `censoring.type='both'` and only blq censoring was performed, for example, they will receive an error stating e.g., `pctalq data.frame was not found in tidyvpcobj. Use censoring() to create censored data for plotting alq`.
 
 ## tidyvpc 1.4.0