diff --git a/NEWS.md b/NEWS.md
index 8c060db..71a322b 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,8 +1,10 @@
 # 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/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).
 * Usability enhancements for prediction corrected VPC (pcVPC), which include support for `binning.tidyvpcobj()` either before or after usage of `predcorrect.tidyvpcobj()`, and automatically performing LOESS pcVPC when `binless.tidyvpcobj()` is used. As a result, the `loess.ypc` argument is no longer required[#43](https://github.com/certara/tidyvpc/issues/43).
+* The `binless.tidyvpcobj` function is now compatible with usage of `censoring.tidyvpcobj` with ALQ data, in addition to BLQ data.[#49](https://github.com/certara/tidyvpc/issues/49)
 * VPC can work with a single value in a group [#51](https://github.com/certara/tidyvpc/issues/51)
 
 # tidyvpc 1.4.0
diff --git a/R/binless.R b/R/binless.R
index a1c237c..9654b0e 100644
--- a/R/binless.R
+++ b/R/binless.R
@@ -276,7 +276,7 @@ binlessfit <- function(o, conf.level = .95, llam.quant = NULL, span = NULL, ...)
 }
 
 .binlessvpcstats <-  function(o, qpred=c(0.05, 0.5, 0.95), conf.level=0.95, quantile.type=7, vpc.type, ...){
-  y <- x <- blq <- fit <- . <- repl <- cprop <- rqssmed <- llam.med <- c.rqssmed <-  NULL
+  y <- x <- blq <- fit <- . <- repl <- cprop <- rqssmed <- llam.med <- c.rqssmed <- lloq <- uloq <- alq <- NULL
 
   obs.fits <- o$rqss.obs.fits
   sim.fits <- o$rqss.sim.fits
diff --git a/R/plot.R b/R/plot.R
index e4e56c7..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()`
@@ -194,6 +196,7 @@ plot.tidyvpcobj <- function(x,
 #'   single-value groups
 #' @noRd
 expand_vpc_stats_single_value <- function(vpc, xvar, width = 0.0001) {
+  n_xvar <- NULL
   d_vpc_stats <- vpc$stats
   if (!is.null(vpc$strat)) {
     d_vpc_stats[, n_xvar := length(unique(get(xvar))), by = names(vpc$strat)]
@@ -508,7 +511,7 @@ plot_censored <-
            show.binning) {
 
     stopifnot(inherits(vpc, "tidyvpcobj"))
-    hi <- lo <- md <- xbin <- y <- NULL
+    hi <- lo <- md <- xbin <- y <- x <- xleft <- xright <- blq <- alq <- NULL
     . <- list
 
     method <- vpc$vpc.method$method
diff --git a/R/vpcstats.R b/R/vpcstats.R
index c461410..73fd805 100644
--- a/R/vpcstats.R
+++ b/R/vpcstats.R
@@ -129,6 +129,7 @@ simulated.tidyvpcobj <- function(o, data, ysim, ..., xsim) {
 #'  values for lower/upper limit of quantification. Logicals for \code{blq} and \code{alq} are returned that indicate whether the DV value lies below/above limit
 #'  of quantification.
 #' @examples
+#' \donttest{
 #' require(magrittr)
 #'
 #' vpc <- observed(obs_data, x=TIME, y=DV) %>%
@@ -147,7 +148,7 @@ simulated.tidyvpcobj <- function(o, data, ysim, ..., xsim) {
 #'     stratify(~ STUDY) %>%
 #'     binning(bin = "kmeans", nbins = 4) %>%
 #'     vpcstats()
-#'
+#' }
 #'
 #' @seealso \code{\link{observed}} \code{\link{simulated}} \code{\link{stratify}} \code{\link{predcorrect}} \code{\link{binning}} \code{\link{binless}} \code{\link{vpcstats}}
 
@@ -245,6 +246,7 @@ censoring.tidyvpcobj <- function(o, blq, lloq, alq, uloq, data=o$data, ...) {
 #'   is \code{obs} split by unique levels of stratification variable(s). Resulting datasets are of class object \code{data.frame}
 #'   and \code{data.table}.
 #' @examples
+#' \donttest{
 #' require(magrittr)
 #'
 #' vpc <- observed(obs_data, x=TIME, y=DV) %>%
@@ -259,6 +261,7 @@ censoring.tidyvpcobj <- function(o, blq, lloq, alq, uloq, data=o$data, ...) {
 #'     stratify(~ GENDER + STUDY) %>%
 #'     binning(bin = "centers", centers = c(1,3,5,7,10)) %>%
 #'     vpcstats()
+#' }
 #'
 #' @seealso \code{\link{observed}} \code{\link{simulated}} \code{\link{censoring}} \code{\link{predcorrect}} \code{\link{binning}} \code{\link{binless}} \code{\link{vpcstats}}
 
@@ -372,7 +375,7 @@ binning <- function(o, ...) UseMethod("binning")
 #' @rdname binning
 #' @export
 binning.tidyvpcobj <- function(o, bin, data=o$data, xbin="xmedian", centers, breaks, nbins, altx, stratum=NULL, by.strata=TRUE,  ...) {
-  keep <- i <- NULL
+  keep <- i <- ypc <- y <- NULL
   . <- list
 
   # If xbin is numeric, then that is the bin
@@ -557,6 +560,7 @@ binning.tidyvpcobj <- function(o, bin, data=o$data, xbin="xmedian", centers, bre
 #'   prediction corrected VPC is to be performed, the \code{predcor.log} logical indicating whether the DV is on a log-scale, and the \code{pred} prediction
 #'   column from the original data.
 #' @examples
+#' \donttest{
 #' require(magrittr)
 #'
 #' obs_data <- obs_data[MDV == 0]
@@ -573,7 +577,6 @@ binning.tidyvpcobj <- function(o, bin, data=o$data, xbin="xmedian", centers, bre
 #'        predcorrect(pred=PRED) %>%
 #'        vpcstats()
 #'
-#'  \donttest{
 #'  # For binless loess prediction corrected, use predcorrect() before
 #'  # binless() and set loess.ypc = TRUE
 #'
diff --git a/cran-comments.md b/cran-comments.md
index 67a85d3..0a6e6c5 100644
--- a/cran-comments.md
+++ b/cran-comments.md
@@ -1,11 +1,13 @@
 ## Release summary
 
-This is a minor release that provides support for the following additional binning methods available in `classInt v0.4.8`: `headtails`, `maximum`, and `box`. In addition, a [fix](https://github.com/certara/tidyvpc/pull/16) has been made for an issue found inside the `npde()` function.
+This is a minor release `1.5.0` that provides misc usability enhancements, including the ability to plot percentage of BLQ and/or ALQ in `plot.tidyvpcobj` if `censoring()` is used, support for `binning()` either before or after usage of `predcorrect()`, and additional validation checks within `predcorrect()` and `simulated()` functions to identify potential issues within data.
 
 ## Test environments
 
 * Windows 10 Enterprise, R 4.2.1
-* Ubuntu 20.04.5, R 4.2.1
+* Windows Server 2022, R 4.3.1
+* macOS 12.6.9, R 4.3.1
+* Ubuntu 22.04.3, R 4.3.1
 
 ## R CMD check results
 
diff --git a/man/censoring.Rd b/man/censoring.Rd
index 3d246de..6042efe 100644
--- a/man/censoring.Rd
+++ b/man/censoring.Rd
@@ -33,6 +33,7 @@ Updates \code{obs} \code{data.frame} in \code{tidypcobj} with censored values fo
 Specify censoring variable or censoring value for VPC.
 }
 \examples{
+\donttest{
 require(magrittr)
 
 vpc <- observed(obs_data, x=TIME, y=DV) \%>\%
@@ -51,7 +52,7 @@ vpc <- observed(obs_data, x=TIME, y=DV) \%>\%
     stratify(~ STUDY) \%>\%
     binning(bin = "kmeans", nbins = 4) \%>\%
     vpcstats()
-
+}
 
 }
 \seealso{
diff --git a/man/predcorrect.Rd b/man/predcorrect.Rd
index f20d162..b3470a8 100644
--- a/man/predcorrect.Rd
+++ b/man/predcorrect.Rd
@@ -29,6 +29,7 @@ Updates \code{tidyvpcobj} with required information to performing prediction cor
 Specify prediction variable for pcVPC.
 }
 \examples{
+\donttest{
 require(magrittr)
 
 obs_data <- obs_data[MDV == 0]
@@ -45,7 +46,6 @@ obs_data$PRED <- sim_data[REP == 1, PRED]
        predcorrect(pred=PRED) \%>\%
        vpcstats()
 
- \donttest{
  # For binless loess prediction corrected, use predcorrect() before
  # binless() and set loess.ypc = TRUE
 
diff --git a/man/stratify.Rd b/man/stratify.Rd
index 254e4c6..268a255 100644
--- a/man/stratify.Rd
+++ b/man/stratify.Rd
@@ -27,6 +27,7 @@ Returns updated \code{tidyvpcobj} with stratification formula, stratification co
 Use to specify stratification variables for VPC.
 }
 \examples{
+\donttest{
 require(magrittr)
 
 vpc <- observed(obs_data, x=TIME, y=DV) \%>\%
@@ -41,6 +42,7 @@ vpc <- vpc \%>\%
     stratify(~ GENDER + STUDY) \%>\%
     binning(bin = "centers", centers = c(1,3,5,7,10)) \%>\%
     vpcstats()
+}
 
 }
 \seealso{
diff --git a/tests/testthat.R b/tests/testthat.R
index 190005c..330e711 100644
--- a/tests/testthat.R
+++ b/tests/testthat.R
@@ -1,4 +1,6 @@
 library(testthat)
 library(tidyvpc)
-
+registered_threads <- data.table::getDTthreads()
+data.table::setDTthreads(2)
+on.exit(data.table::setDTthreads(registered_threads))
 test_check("tidyvpc")
diff --git a/tests/testthat/test-binless.R b/tests/testthat/test-binless.R
index 676c07b..dbf43ab 100644
--- a/tests/testthat/test-binless.R
+++ b/tests/testthat/test-binless.R
@@ -114,6 +114,7 @@ test_that("cont vpc predcorrect binless vpcstats are correct", {
 })
 
 test_that("cont vpc binless censoring vpcstats are correct", {
+  skip_on_cran()
   obs_data <- obs_data[MDV == 0]
   sim_data <- sim_data[MDV == 0]
 
@@ -200,6 +201,7 @@ test_that("binless.tidyvpcobj returns correct errors and warnings", {
 
 test_that("binless.tidyvpcobj uses supplied lambda and span if optimize = FALSE", {
   # continuous VPC
+  skip_on_cran()
   obs_data <- obs_data[MDV == 0]
   sim_data <- sim_data[MDV == 0]
   obs_data$PRED <- sim_data[REP == 1, PRED]
diff --git a/vignettes/tidyvpc_cat.Rmd b/vignettes/tidyvpc_cat.Rmd
index 1bfe7f2..86ee0b4 100644
--- a/vignettes/tidyvpc_cat.Rmd
+++ b/vignettes/tidyvpc_cat.Rmd
@@ -18,6 +18,10 @@ library(data.table)
 set.seed(1014)
 ```
 
+```{r, echo=FALSE}
+data.table::setDTthreads(2)
+```
+
 
 ## Introduction
 
diff --git a/vignettes/tidyvpc_cont.Rmd b/vignettes/tidyvpc_cont.Rmd
index 72e8c5f..436d007 100644
--- a/vignettes/tidyvpc_cont.Rmd
+++ b/vignettes/tidyvpc_cont.Rmd
@@ -16,6 +16,11 @@ library(magrittr)
 set.seed(1014)
 ```
 
+```{r, echo=FALSE}
+data.table::setDTthreads(2)
+```
+
+
 ## Introduction 
 
 **When deriving a Visual Predictive Check (VPC) you must:**
@@ -229,44 +234,57 @@ vpc <- observed(obs_data, x=TIME, y=DV) %>%
     binning(bin = "jenks", nbins = 5) %>%
     vpcstats()
 
-plot(vpc)
-```
-
-Censoring using `LLOQ` variable in the data.
-
-First, let's add an LLOQ variable to the observed data.
-```{r}
-obs_data$LLOQ <- 50
+plot(vpc, censoring.type = "blq")
 ```
 
-Then, we'll specify lower limit of quantification as the unquoted variable name `LLOQ` in our data. Let's also provide our own lambda values.
+The `tidyvpc` package also allows you to use LLOQ values within your data and 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)]
+
 vpc <- observed(obs_data, x=TIME, y=DV) %>%
     simulated(sim_data, y=DV) %>%
     censoring(blq=(DV < LLOQ), lloq=LLOQ) %>%
-    binless(optimize = FALSE, lambda = c(1.5, 2.5, 1.7)) %>%
-    vpcstats()
+    stratify(~ STUDY) %>%
+    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.
 
-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.
+```{r  fig.width = 9, fig.height = 6, out.width=640, warning = FALSE, eval = 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))
+```
+
+To plot the above `vpc` and include the percentage of ALQ plot, you'd run `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) %>%
-    stratify(~ STUDY) %>%
-    binning(bin = "pam", nbins = 4) %>%
-    vpcstats(qpred = c(0.1, 0.5, 0.9))
+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)
+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_nlmixr2.Rmd b/vignettes/tidyvpc_nlmixr2.Rmd
index 3511cff..404ac55 100644
--- a/vignettes/tidyvpc_nlmixr2.Rmd
+++ b/vignettes/tidyvpc_nlmixr2.Rmd
@@ -14,6 +14,10 @@ knitr::opts_chunk$set(
 )
 ```
 
+```{r, echo=FALSE}
+data.table::setDTthreads(2)
+```
+
 # Introduction
 
 `tidyvpc` and `nlmixr2` can work together seamlessly.  The information below
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