diff --git a/R/binless.R b/R/binless.R
index 9654b0e..77bfc5c 100644
--- a/R/binless.R
+++ b/R/binless.R
@@ -147,7 +147,12 @@ binlessfit <- function(o, conf.level = .95, llam.quant = NULL, span = NULL, ...)
   obs <- o$obs
   sim <- o$sim
 
-  if(isTRUE(o$loess.ypc)) {
+  if (isTRUE(o$loess.ypc)) {
+    if (isTRUE(o$varcorr)) {
+      warning("Variability correction is not supported by binless VPC and won't be applied")
+      o$varcorr <- FALSE
+    }
+    
     pred <- o$pred
     obs <- cbind(obs, pred)
     sim[, pred := rep(pred, len=.N), by = .(repl)]
@@ -218,13 +223,13 @@ binlessfit <- function(o, conf.level = .95, llam.quant = NULL, span = NULL, ...)
     strat.split.sim <- strat.split.sim[lapply(strat.split.sim,NROW)>0]
   }
 
-  if(isTRUE(o$loess.ypc) && !is.null(o$strat)) {
-    if(isTRUE(o$predcor.log)) {
-      for(i in seq_along(strat.split)) {
+  if (isTRUE(o$loess.ypc) && !is.null(o$strat)) {
+    if (isTRUE(o$predcor.log)) {
+      for (i in seq_along(strat.split)) {
         strat.split[[i]][, l.ypc := y +  (fitted(loess(pred ~ x, span = span[[i]], na.action = na.exclude, .SD)) - pred)]
       }
     } else {
-      for(i in seq_along(strat.split)) {
+      for (i in seq_along(strat.split)) {
         strat.split[[i]][, l.ypc := y * (fitted(loess(pred ~ x, span = span[[i]], na.action = na.exclude, .SD)) / pred)]
       }
     }
diff --git a/R/plot.R b/R/plot.R
index 988231c..3ce7eae 100644
--- a/R/plot.R
+++ b/R/plot.R
@@ -122,9 +122,15 @@ plot.tidyvpcobj <- function(x,
       ylab <- NULL
     }
     if (isTRUE(vpc$predcor)) {
-      ylab <- ifelse(length(ylab) == 0,
-                     "Prediction Corrected",
-                     paste0(ylab, "\nPrediction Corrected"))
+      if (isTRUE(vpc$varcorr)) {
+        ylab <- ifelse(length(ylab) == 0,
+                       "Prediction and Variability Corrected",
+                       paste0(ylab, "\nPrediction and Variability Corrected"))        
+      } else {
+        ylab <- ifelse(length(ylab) == 0,
+                       "Prediction Corrected",
+                       paste0(ylab, "\nPrediction Corrected"))
+      }
     }
   }
 
@@ -233,7 +239,8 @@ plot_continuous <-
            point.shape,
            point.stroke,
            point.alpha) {
-    alq <-bin <- blq <- hi <-l.ypc <-lo <- md <- pname <- qname <- x <- xleft <- xright <- y <- ypc <- NULL
+    alq <- bin <- blq <- hi <- l.ypc <- lo <- md <- pname <- qname <- NULL
+    x <- xleft <- xright <- y <- ypc <- ypcvc <- NULL
     . <- list
     method <- vpc$vpc.method$method
     qlvls <- levels(vpc$stats$qname)
@@ -306,8 +313,13 @@ plot_continuous <-
       if (isTRUE(vpc$predcor) && method == "binless") {
         points.dat[, y := l.ypc]
       } else if (isTRUE(vpc$predcor)) {
-        points.dat[, y := ypc]
+        if (isTRUE(vpc$varcorr)) {
+          points.dat[, y := ypcvc]
+        } else {
+          points.dat[, y := ypc]
+        }
       }
+      
       if (show.binning) {
         reorder2 <- function(y, x) {
           y <- stats::reorder(y, x)
diff --git a/R/vpcstats.R b/R/vpcstats.R
index 4ff1501..d68d0ff 100644
--- a/R/vpcstats.R
+++ b/R/vpcstats.R
@@ -21,19 +21,25 @@ NULL
 
 #' Specify observed dataset and variables for VPC
 #'
-#' The observed function is the first function in the vpc piping chain and is used for specifying observed data and variables for VPC. Note: Observed
-#' data must not contain missing DV and may require filtering \code{MDV == 0} before generating VPC.
+#' The observed function is the first function in the vpc piping chain and is
+#' used for specifying observed data and variables for VPC. Note: Observed data
+#' must not contain missing DV and may require filtering \code{MDV == 0} before
+#' generating VPC. Also observed data must be ordered by: Subject (ID), IVAR
+#' (Time)
 #'
 #' @param o A \code{data.frame} of observation data.
 #' @param x Numeric x-variable, typically named TIME.
 #' @param yobs Numeric y-variable, typically named DV.
 #' @param pred Population prediction variable, typically named PRED.
 #' @param blq Logical variable indicating below limit of quantification.
-#' @param lloq Number or numeric variable in data indicating the lower limit of quantification.
+#' @param lloq Number or numeric variable in data indicating the lower limit of
+#'   quantification.
 #' @param alq Logical variable indicating above limit of quantification .
-#' @param uloq Number or numeric variable in data indicating the upper limit of quantification.
+#' @param uloq Number or numeric variable in data indicating the upper limit of
+#'   quantification.
 #' @param ... Other arguments.
-#' @return A \code{tidyvpcobj} containing both original data and observed data formatted with \code{x} and \code{y} variables as specified in function.
+#' @return A \code{tidyvpcobj} containing both original data and observed data
+#'   formatted with \code{x} and \code{y} variables as specified in function.
 #'   Resulting data is of class \code{data.frame} and \code{data.table}.
 #' @examples
 #'
@@ -67,10 +73,11 @@ observed.data.frame <- function(o, x, yobs, pred=NULL, blq=NULL, lloq=-Inf, alq=
 
 #' Specify simulated dataset and variables for VPC
 #'
-#' The simulated function is used for specifying simulated input data and variables for VPC. Note: Simulated data must not
-#' contain missing DV and may require filtering \code{MDV == 0} before generating VPC. The ordering of observed and simulated
-#' data must also be consistent, with replicates in simulated data stacked on top of each other.
-#'
+#' The simulated function is used for specifying simulated input data and
+#' variables for VPC. Note: Simulated data must not contain missing DV and may
+#' require filtering \code{MDV == 0} before generating VPC. Simulated data must
+#' be ordered by: Replicate, Subject (ID), IVAR (Time).
+#' 
 #' @param o A \code{tidyvpcobj}.
 #' @param data A \code{data.frame} of simulated data.
 #' @param ysim Numeric y-variable, typically named DV.
@@ -375,7 +382,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 <- ypc <- y <- NULL
+  keep <- i <- NULL
   . <- list
 
   # If xbin is numeric, then that is the bin
@@ -484,7 +491,7 @@ binning.tidyvpcobj <- function(o, bin, data=o$data, xbin="xmedian", centers, bre
   }
 
   if (is.function(bin)) {
-    xdat <- data.table(i=1:nrow(o$obs), x=x)
+    xdat <- data.table(i = 1:nrow(o$obs), x = x)
     if (any(is.na(xdat[filter]$x))) {
       warning("x contains missing values, which could affect binning")
     }
@@ -493,7 +500,7 @@ binning.tidyvpcobj <- function(o, bin, data=o$data, xbin="xmedian", centers, bre
     }
     if (by.strata && !is.null(o$strat)) {
       sdat <- copy(o$strat)
-      temp <- xdat[filter, .(i=i, j=do.call(bin, c(list(x), args, .BY))), by=sdat[filter]]
+      temp <- xdat[filter, .(i = i, j = do.call(bin, c(list(x), args, .BY))), by = sdat[filter]]
       j <- temp[order(i), j]
     } else {
       j <- xdat[filter, do.call(bin, c(list(x), args))]
@@ -527,38 +534,48 @@ binning.tidyvpcobj <- function(o, bin, data=o$data, xbin="xmedian", centers, bre
   } else {
     stop("Invalid xbin")
   }
+  
   vpc.method <- list(method = "binning")
 
   # check if user supplied predcorrect before binning
-  if (!is.null(o$predcor) && o$predcor) {
-    pred <- o$pred
-    log <- o$predcor.log
-    mpred <- data.table(stratbin, pred)[, mpred := median(pred), by = stratbin]$mpred
-
-    if (log) {
-      o$obs[, ypc := (mpred - pred) + y]
-      o$sim[, ypc := (mpred - pred) + y]
-    } else {
-      o$obs[, ypc := ifelse(pred == 0, 0, (mpred / pred) * y)]
-      o$sim[, ypc := ifelse(rep(pred, times = nrow(o$sim) / nrow(o$obs)) == 0, 0, (mpred / pred) * y)]
-    }
-  }
-
-  update(o, xbin=xbin, vpc.method = vpc.method)
+  if (isTRUE(o$predcor)) {
+    o <-
+      get_predcorrect(
+        o,
+        stratbin,
+        pred = o$pred,
+        log = o$predcor.log,
+        varcorr = o$varcorr
+      )
+  } else if (isTRUE(o$varcorr)) {
+    warning("Cannot apply variance prediction correction when predcor flag is off.")
+  }
+
+  update(o, xbin = xbin, vpc.method = vpc.method)
 }
 
+
 #' Prediction corrected Visual Predictive Check (pcVPC)
 #'
 #' Specify prediction variable for pcVPC.
 #'
-#' @param o A \code{tidyvpcobj}.
+#' @param o A `tidyvpcobj`.
 #' @param pred Prediction variable in observed data.
-#' @param data Observed data supplied in \code{observed()} function.
+#' @param data Observed data supplied in `observed()` function.
 #' @param ... Other arguments to include.
 #' @param log Logical indicating whether DV was modeled in logarithmic scale.
-#' @return Updates \code{tidyvpcobj} with required information to performing prediction correction, which includes the \code{predcor} logical indicating whether
-#'   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.
+#' @param varcorr Logical indicating whether variability correction should be
+#'   applied for prediction corrected dependent variable
+#' @return Updates `tidyvpcobj` with required information to perform prediction
+#'   correction, which includes the `predcor` logical indicating whether
+#'   prediction corrected VPC is to be performed, the `predcor.log` logical
+#'   indicating whether the DV is on a log-scale, the `varcorr` logical
+#'   indicating whether variability correction for prediction corrected
+#'   dependent variable is applied and the `pred` prediction column from the
+#'   original data. Both `obs` and `sim` data tables in the returned
+#'   `tidyvpcobj` object have additional `ypc` column with the results of
+#'   prediction correction and `ypcvc` column if variability correction is
+#'   requested.
 #' @examples
 #' \donttest{
 #' require(magrittr)
@@ -571,62 +588,92 @@ binning.tidyvpcobj <- function(o, bin, data=o$data, xbin="xmedian", centers, bre
 #'
 #' obs_data$PRED <- sim_data[REP == 1, PRED]
 #'
-#'   vpc <- observed(obs_data, x=TIME, y=DV) %>%
-#'        simulated(sim_data, y=DV) %>%
+#'   vpc <- observed(obs_data, x=TIME, yobs=DV) %>%
+#'        simulated(sim_data, ysim=DV) %>%
 #'        binning(bin = NTIME) %>%
-#'        predcorrect(pred=PRED) %>%
+#'        predcorrect(pred=PRED, varcorr = TRUE) %>%
 #'        vpcstats()
 #'
 #'  # For binless loess prediction corrected, use predcorrect() before
 #'  # binless() and set loess.ypc = TRUE
 #'
-#'   vpc <- observed(obs_data, x=TIME, y=DV) %>%
-#'        simulated(sim_data, y=DV) %>%
+#'   vpc <- observed(obs_data, x=TIME, yobs=DV) %>%
+#'        simulated(sim_data, ysim=DV) %>%
 #'        predcorrect(pred=PRED) %>%
-#'        binless(loess.ypc = TRUE) %>%
+#'        binless() %>%
 #'        vpcstats()
 #'        }
 #'
-#' @seealso \code{\link{observed}} \code{\link{simulated}} \code{\link{censoring}} \code{\link{stratify}} \code{\link{binning}} \code{\link{binless}} \code{\link{vpcstats}}
+#' @seealso \code{\link{observed}} \code{\link{simulated}}
+#'   \code{\link{censoring}} \code{\link{stratify}} \code{\link{binning}}
+#'   \code{\link{binless}} \code{\link{vpcstats}}
 
 #' @export
 predcorrect <- function(o, ...) UseMethod("predcorrect")
 
 #' @rdname predcorrect
 #' @export
-predcorrect.tidyvpcobj <- function(o, pred, data=o$data, ..., log=FALSE) {
-
-  ypc <- y <- NULL
-
-  if (missing(pred)) {
-    pred <- o$pred
-  } else {
-    pred <- rlang::eval_tidy(rlang::enquo(pred), data)
-  }
-  if (is.null(pred)) {
-    stop("No pred specified")
-  }
-
-  stratbin <- o$.stratbin
-  # predcorrect after binning, check if binning/binless has already been specified
-
-  if (!is.null(o$vpc.method)) {
-    if(o$vpc.method$method == "binless") {
-      o$vpc.method$loess.ypc <- TRUE
-    } else { #binning specified, perform ypc calculcation
-      mpred <- data.table(stratbin, pred)[, mpred := median(pred), by = stratbin]$mpred
-
-      if (log) {
-        o$obs[, ypc := (mpred - pred) + y]
-        o$sim[, ypc := (mpred - pred) + y]
+predcorrect.tidyvpcobj <-
+  function(o,
+           pred,
+           data = o$data,
+           ...,
+           log = FALSE,
+           varcorr = FALSE) {
+    if (missing(pred)) {
+      pred <- o$pred
+    } else {
+      pred <- rlang::eval_tidy(rlang::enquo(pred), data)
+    }
+    if (is.null(pred)) {
+      stop("No pred specified")
+    }
+    
+    stratbin <- o$.stratbin
+    # predcorrect after binning, check if binning/binless has already been specified
+    
+    if (!is.null(o$vpc.method)) {
+      if (o$vpc.method$method == "binless") {
+        o$vpc.method$loess.ypc <- TRUE
       } else {
-        o$obs[, ypc := ifelse(pred == 0, 0, (mpred / pred) * y)]
-        o$sim[, ypc := ifelse(rep(pred, times = nrow(o$sim) / nrow(o$obs)) == 0, 0, (mpred / pred) * y)]
+        #binning specified, perform ypc calculcation
+        o <- get_predcorrect(o, stratbin, pred, log, varcorr)
       }
     }
-  }
-
-  update(o, predcor=TRUE, predcor.log=log, pred=pred)
+    
+    update(
+      o,
+      predcor = TRUE,
+      predcor.log = log,
+      varcorr = varcorr,
+      pred = pred
+    )
+  }
+
+get_predcorrect <- function(o, stratbin, pred, log, varcorr) {
+  ypcvc <- ypc <- y <- ij <- NULL
+  . <- list
+  mpred <-
+    data.table(stratbin, pred)[, mpred := median(pred), by = stratbin]$mpred
+  
+  if (log) {
+    o$obs[, ypc := (mpred - pred) + y]
+    o$sim[, ypc := (mpred - pred) + y]
+  } else {
+    o$obs[, ypc := ifelse(pred == 0, 0, (mpred / pred) * y)]
+    o$sim[, ypc := ifelse(rep(pred, times = nrow(o$sim) / nrow(o$obs)) == 0, 0, (mpred / pred) * y)]
+  }
+  
+  if (varcorr) {
+    ypcsdij <-
+      data.table(ypc = o$sim$ypc, ij = row.names(o$obs))[, .(ypcsdij = stats::sd(ypc)), by = ij]$ypcsdij
+    mypcsdij <-
+      data.table(stratbin, ypcsdij)[, mypcsdij := median(ypcsdij), by = stratbin]$mypcsdij
+    o$obs[, ypcvc := mpred + (y - mpred) * mypcsdij / ypcsdij]
+    o$sim[, ypcvc := mpred + (y - mpred) * mypcsdij / ypcsdij]
+  }
+  
+  o
 }
 
 #' Remove prediction correction for Visual Predictive Check (VPC)
@@ -641,7 +688,7 @@ nopredcorrect <- function(o, ...) UseMethod("nopredcorrect")
 #' @rdname nopredcorrect
 #' @export
 nopredcorrect.tidyvpcobj <- function(o, ...) {
-  update(o, predcor=FALSE)
+  update(o, predcor = FALSE)
 }
 
 #' @export
diff --git a/R/vpcstats_fun.R b/R/vpcstats_fun.R
index 5ad37ef..ce515b6 100644
--- a/R/vpcstats_fun.R
+++ b/R/vpcstats_fun.R
@@ -29,16 +29,22 @@ vpcstats.tidyvpcobj <- function(o, vpc.type =c("continuous", "categorical"), qpr
   type <- match.arg(vpc.type)
   method <- o$vpc.method
 
-  stopifnot(method$method %in% c("binless", "binning"))
+  if (!isTRUE(method$method %in% c("binless", "binning"))) {
+    stop(
+      "A binning method should be specified through binning() or binless() execution before vpcstats() run."
+    )
+  }
+  
   stopifnot(length(qpred) == 3)
 
-  repl <- ypc <- y <- x <- blq <- lloq <- alq <- uloq <- NULL
+  repl <- ypc <- ypcvc <- y <- x <- blq <- lloq <- alq <- uloq <- NULL
   . <- list
   qconf <- c(0, 0.5, 1) + c(1, 0, -1)*(1 - conf.level)/2
 
   obs      <- o$obs
   sim      <- o$sim
   predcor  <- o$predcor
+  varcorr  <- o$varcorr
   stratbin <- o$.stratbin
   xbin     <- o$xbin
   strat    <- o$strat
@@ -197,8 +203,13 @@ vpcstats.tidyvpcobj <- function(o, vpc.type =c("continuous", "categorical"), qpr
       .stratbinrepl <- data.table(stratbin, sim[, .(repl)])
 
       if (isTRUE(predcor)) {
-        qobs <- obs[, quant_loq(ypc, probs=qpred, blq=blq,   alq=alq,   type = quantile.type),   by=stratbin]
-        qsim <- sim[, quant_loq(ypc, probs=qpred, blq=FALSE, alq=FALSE, type = quantile.type), by=.stratbinrepl]
+        if (isTRUE(varcorr)) {
+          qobs <- obs[, quant_loq(ypcvc, probs=qpred, blq=blq,   alq=alq,   type = quantile.type),   by=stratbin]
+          qsim <- sim[, quant_loq(ypcvc, probs=qpred, blq=FALSE, alq=FALSE, type = quantile.type), by=.stratbinrepl]
+        } else {
+          qobs <- obs[, quant_loq(ypc, probs=qpred, blq=blq,   alq=alq,   type = quantile.type),   by=stratbin]
+          qsim <- sim[, quant_loq(ypc, probs=qpred, blq=FALSE, alq=FALSE, type = quantile.type), by=.stratbinrepl]
+        }
       } else {
         qobs <- obs[, quant_loq(y, probs=qpred, blq=blq,   alq=alq  , type = quantile.type),   by=stratbin]
         qsim <- sim[, quant_loq(y, probs=qpred, blq=FALSE, alq=FALSE, type = quantile.type), by=.stratbinrepl]
diff --git a/man/npde.Rd b/man/npde.Rd
index b97b109..00234eb 100644
--- a/man/npde.Rd
+++ b/man/npde.Rd
@@ -42,7 +42,7 @@ vpc <- observed(npde$npdeobs, x=epred, y=npde) \%>\%
     binning("eqcut", nbins=10) \%>\%
     vpcstats()
 
-plot(vpc) + 
+plot(vpc) +
 labs(x="Simulation-based Population Prediction", y="Normalized Prediction Distribution Error")
 }
 
diff --git a/man/observed.Rd b/man/observed.Rd
index 92ff931..7765953 100644
--- a/man/observed.Rd
+++ b/man/observed.Rd
@@ -32,19 +32,25 @@ observed(o, ...)
 
 \item{blq}{Logical variable indicating below limit of quantification.}
 
-\item{lloq}{Number or numeric variable in data indicating the lower limit of quantification.}
+\item{lloq}{Number or numeric variable in data indicating the lower limit of
+quantification.}
 
 \item{alq}{Logical variable indicating above limit of quantification .}
 
-\item{uloq}{Number or numeric variable in data indicating the upper limit of quantification.}
+\item{uloq}{Number or numeric variable in data indicating the upper limit of
+quantification.}
 }
 \value{
-A \code{tidyvpcobj} containing both original data and observed data formatted with \code{x} and \code{y} variables as specified in function.
+A \code{tidyvpcobj} containing both original data and observed data
+  formatted with \code{x} and \code{y} variables as specified in function.
   Resulting data is of class \code{data.frame} and \code{data.table}.
 }
 \description{
-The observed function is the first function in the vpc piping chain and is used for specifying observed data and variables for VPC. Note: Observed
-data must not contain missing DV and may require filtering \code{MDV == 0} before generating VPC.
+The observed function is the first function in the vpc piping chain and is
+used for specifying observed data and variables for VPC. Note: Observed data
+must not contain missing DV and may require filtering \code{MDV == 0} before
+generating VPC. Also observed data must be ordered by: Subject (ID), IVAR
+(Time)
 }
 \examples{
 
diff --git a/man/predcorrect.Rd b/man/predcorrect.Rd
index b3470a8..d9744f0 100644
--- a/man/predcorrect.Rd
+++ b/man/predcorrect.Rd
@@ -7,23 +7,33 @@
 \usage{
 predcorrect(o, ...)
 
-\method{predcorrect}{tidyvpcobj}(o, pred, data = o$data, ..., log = FALSE)
+\method{predcorrect}{tidyvpcobj}(o, pred, data = o$data, ..., log = FALSE, varcorr = FALSE)
 }
 \arguments{
-\item{o}{A \code{tidyvpcobj}.}
+\item{o}{A `tidyvpcobj`.}
 
 \item{...}{Other arguments to include.}
 
 \item{pred}{Prediction variable in observed data.}
 
-\item{data}{Observed data supplied in \code{observed()} function.}
+\item{data}{Observed data supplied in `observed()` function.}
 
 \item{log}{Logical indicating whether DV was modeled in logarithmic scale.}
+
+\item{varcorr}{Logical indicating whether variability correction should be
+applied for prediction corrected dependent variable}
 }
 \value{
-Updates \code{tidyvpcobj} with required information to performing prediction correction, which includes the \code{predcor} logical indicating whether
-  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.
+Updates `tidyvpcobj` with required information to perform prediction
+  correction, which includes the `predcor` logical indicating whether
+  prediction corrected VPC is to be performed, the `predcor.log` logical
+  indicating whether the DV is on a log-scale, the `varcorr` logical
+  indicating whether variability correction for prediction corrected
+  dependent variable is applied and the `pred` prediction column from the
+  original data. Both `obs` and `sim` data tables in the returned
+  `tidyvpcobj` object have additional `ypc` column with the results of
+  prediction correction and `ypcvc` column if variability correction is
+  requested.
 }
 \description{
 Specify prediction variable for pcVPC.
@@ -40,23 +50,25 @@ sim_data <- sim_data[MDV == 0]
 
 obs_data$PRED <- sim_data[REP == 1, PRED]
 
-  vpc <- observed(obs_data, x=TIME, y=DV) \%>\%
-       simulated(sim_data, y=DV) \%>\%
+  vpc <- observed(obs_data, x=TIME, yobs=DV) \%>\%
+       simulated(sim_data, ysim=DV) \%>\%
        binning(bin = NTIME) \%>\%
-       predcorrect(pred=PRED) \%>\%
+       predcorrect(pred=PRED, varcorr = TRUE) \%>\%
        vpcstats()
 
  # For binless loess prediction corrected, use predcorrect() before
  # binless() and set loess.ypc = TRUE
 
-  vpc <- observed(obs_data, x=TIME, y=DV) \%>\%
-       simulated(sim_data, y=DV) \%>\%
+  vpc <- observed(obs_data, x=TIME, yobs=DV) \%>\%
+       simulated(sim_data, ysim=DV) \%>\%
        predcorrect(pred=PRED) \%>\%
-       binless(loess.ypc = TRUE) \%>\%
+       binless() \%>\%
        vpcstats()
        }
 
 }
 \seealso{
-\code{\link{observed}} \code{\link{simulated}} \code{\link{censoring}} \code{\link{stratify}} \code{\link{binning}} \code{\link{binless}} \code{\link{vpcstats}}
+\code{\link{observed}} \code{\link{simulated}}
+  \code{\link{censoring}} \code{\link{stratify}} \code{\link{binning}}
+  \code{\link{binless}} \code{\link{vpcstats}}
 }
diff --git a/man/simulated.Rd b/man/simulated.Rd
index b43ccc7..7d0991a 100644
--- a/man/simulated.Rd
+++ b/man/simulated.Rd
@@ -25,9 +25,10 @@ A \code{tidyvpcobj} containing simulated dataset \code{sim} formatted with colum
  The column \code{x} is used from the \code{observed()} function. Resulting dataset is of class \code{data.frame} and \code{data.table}.
 }
 \description{
-The simulated function is used for specifying simulated input data and variables for VPC. Note: Simulated data must not
-contain missing DV and may require filtering \code{MDV == 0} before generating VPC. The ordering of observed and simulated
-data must also be consistent, with replicates in simulated data stacked on top of each other.
+The simulated function is used for specifying simulated input data and
+variables for VPC. Note: Simulated data must not contain missing DV and may
+require filtering \code{MDV == 0} before generating VPC. Simulated data must
+be ordered by: Replicate, Subject (ID), IVAR (Time).
 }
 \examples{
 require(magrittr)
diff --git a/tests/testthat/test-vpcstats.R b/tests/testthat/test-vpcstats.R
index 4cec26b..a6c40d2 100644
--- a/tests/testthat/test-vpcstats.R
+++ b/tests/testthat/test-vpcstats.R
@@ -1,22 +1,52 @@
 test_that("simulated.tidyvpcobj detects row count mismatches", {
-  vpcobj_o <- observed(o = data.frame(x = 0:1, y = c(0, 2), pred = c(0, 2.5)), x = x, yobs = y)
-  expect_silent(
-    simulated(vpcobj_o, data = data.frame(x = 0:1, y = c(0, 3)), x = x, y = y)
-  )
+  vpcobj_o <-
+    observed(o = data.frame(
+      x = 0:1,
+      y = c(0, 2),
+      pred = c(0, 2.5)
+    ),
+    x = x,
+    yobs = y)
+  expect_silent(simulated(
+    vpcobj_o,
+    data = data.frame(x = 0:1, y = c(0, 3)),
+    xsim = x,
+    ysim = y
+  ))
   expect_error(
-    simulated(vpcobj_o, data = data.frame(x = 0:2, y = c(0, 3, 4)), x = x, y = y),
+    simulated(
+      vpcobj_o,
+      data = data.frame(x = 0:2, y = c(0, 3, 4)),
+      x = x,
+      ysim = y
+    ),
     regexp = "The number of simulated rows is not a multiple of the number of observed rows.  Ensure that you filtered your observed data to remove MDV rows.",
     fixed = TRUE
   )
 })
 
 test_that("simulated.tidyvpcobj detects x-variable mismatches", {
-  vpcobj_o <- observed(o = data.frame(x = 0:1, y = c(0, 2), pred = c(0, 2.5)), x = x, yobs = y)
-  expect_silent(
-    simulated(vpcobj_o, data = data.frame(x = c(0:1, 0:1), y = c(0, 3, 0, 4)), xsim = x, y = y)
-  )
+  vpcobj_o <-
+    observed(o = data.frame(
+      x = 0:1,
+      y = c(0, 2),
+      pred = c(0, 2.5)
+    ),
+    x = x,
+    yobs = y)
+  expect_silent(simulated(
+    vpcobj_o,
+    data = data.frame(x = c(0:1, 0:1), y = c(0, 3, 0, 4)),
+    xsim = x,
+    ysim = y
+  ))
   expect_error(
-    simulated(vpcobj_o, data = data.frame(x = c(0:2, 1), y = c(0, 3, 0, 4)), xsim = x, y = y),
+    simulated(
+      vpcobj_o,
+      data = data.frame(x = c(0:2, 1), y = c(0, 3, 0, 4)),
+      xsim = x,
+      ysim = y
+    ),
     regexp = "Values of `xsim` do not match observed data x-values.  Ensure that you filtered your observed data to remove MDV rows.",
     fixed = TRUE
   )
@@ -24,9 +54,45 @@ test_that("simulated.tidyvpcobj detects x-variable mismatches", {
 
 test_that("predcorrect.tidyvpcobj", {
   # Prevent division by zero for prediction correction
-  vpcobj_o <- observed(o = data.frame(x = 0:1, y = c(0, 2), pred = c(0, 2.5)), x = x, yobs = y)
-  vpcobj_s <- simulated(vpcobj_o, data = data.frame(x = 0:1, y = c(0, 3)), x = x, y = y)
+  vpcobj_o <-
+    observed(o = data.frame(
+      x = 0:1,
+      y = c(0, 2),
+      pred = c(0, 2.5)
+    ),
+    x = x,
+    yobs = y)
+  vpcobj_s <-
+    simulated(
+      vpcobj_o,
+      data = data.frame(x = 0:1, y = c(0, 3)),
+      xsim = x,
+      ysim = y
+    )
   vpcobj_b <- binning(vpcobj_s, bin = x)
-  vpcobj_predcorr <- predcorrect(vpcobj_b, pred = pred)
+  vpcobj_predcorr <- predcorrect(vpcobj_b, pred = pred, varcorr = TRUE)
   expect_equal(vpcobj_predcorr$obs$ypc, c(0, 2))
 })
+
+test_that("predcorrect.tidyvpcobj with variance correction", {
+  vpcobj_o <-
+    observed(o = data.frame(
+      x = 0:1,
+      y = c(0, 2, 1, 3),
+      pred = c(1, 2, 0, 4)
+    ),
+    x = x,
+    yobs = y)
+  vpcobj_s <-
+    simulated(
+      vpcobj_o,
+      data = data.frame(x = 0:1, y = c(0, 1, 1, 2, 1, 3, 0, 4)),
+      xsim = x,
+      ysim = y
+    )
+  vpcobj_b <- binning(vpcobj_s, bin = x)
+  vpcobj_bp <- predcorrect(vpcobj_b, pred = pred, varcorr = TRUE)
+  vpcobj_p <- predcorrect(vpcobj_s, pred = pred, varcorr = TRUE)
+  vpcobj_pb <- binning(vpcobj_p, bin = x)
+  expect_equal(vpcobj_bp$sim$ypcvc, vpcobj_pb$sim$ypcvc)
+})
\ No newline at end of file