Merge pull request #60 from certara-mtomashevskiy/Variance_Correction…

…_implementation added varcorr argument to predcorrect method
certara · Mar 22, 2024 · 36d329f · 36d329f
2 parents cff5ec6 + ed2a394
commit 36d329f
Show file tree

Hide file tree

Showing 9 changed files with 279 additions and 119 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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]