R/vpc.R

In nlmixr2est: Nonlinear Mixed Effects Models in Population PK/PD, Estimation Routines

.lastPredSimulationInfo <- NULL # to get observation dataset with pred attached for pred_corr
#' VPC simulation
#'
#' @param object This is the nlmixr2 fit object
#' @param ... Other arguments sent to `rxSolve()`
#' @param keep Keep character vector
#' @param n Number of simulations
#' @param pred Should predictions be added to the simulation
#' @param seed Seed to set for the VPC simulation
#' @param nretry Number of times to retry the simulation if there is
#'   NA values in the simulation
#' @param normRelated should the VPC style simulation be for normal
#'   related variables only
#' @return data frame of the VPC simulation
#' @author Matthew L. Fidler
#' @export
#' @examples
#'
#' \donttest{
#' one.cmt <- function() {
#'  ini({
#'    ## You may label each parameter with a comment
#'    tka <- 0.45 # Log Ka
#'    tcl <- log(c(0, 2.7, 100)) # Log Cl
#'    ## This works with interactive models
#'    ## You may also label the preceding line with label("label text")
#'    tv <- 3.45; label("log V")
#'    ## the label("Label name") works with all models
#'    eta.ka ~ 0.6
#'    eta.cl ~ 0.3
#'    eta.v ~ 0.1
#'    add.sd <- 0.7
#'  })
#'  model({
#'    ka <- exp(tka + eta.ka)
#'    cl <- exp(tcl + eta.cl)
#'    v <- exp(tv + eta.v)
#'    linCmt() ~ add(add.sd)
#'  })
#' }
#'
#' fit <- nlmixr(one.cmt, theo_sd, est="focei")
#'
#' head(vpcSim(fit, pred=TRUE))
#'
#' }
vpcSim <- function(object, ..., keep=NULL, n=300,
                   pred=FALSE, seed=1009, nretry=50,
                   normRelated=TRUE) {
  assignInMyNamespace(".finalUiCompressed", FALSE)
  on.exit(assignInMyNamespace(".finalUiCompressed", TRUE))
  set.seed(seed)
  .si <- object$simInfo
  .si$object <- eval(.getSimModel(object, hideIpred=FALSE))
  .w <- which(names(.si) == "rx")
  .si <- .si[-.w]
  .si$nsim <- n
  .si <- c(.si, list(...))
  .pt <- proc.time()
  .si$keep <- unique(c(keep, "nlmixrRowNums"))
  .data <- .si$events
  .data$nlmixrRowNums <- seq_along(.data[, 1])
  if (normRelated) {
    .ui <- object$ui
    .predDf <-.ui$predDf
    if (all(.predDf$dist %in% c("norm", "dnorm","t", "cauchy"))) {
    } else {
      if (is.null(.data$CMT)) {
        .ds <- object$dataSav
        .data$nlmixrRowNums <- seq_along(.data[,1])
        .ds <- .ds[, c("CMT", "nlmixrRowNums")]
        .data <- merge(.data, .ds, by ="nlmixrRowNums")
        .data <- .data[order(.data$nlmixrRowNums),]
      }
      .lst <- .Call(`_nlmixr2est_filterNormalLikeAndDoses`,
                    .data$CMT, .predDf$distribution, .predDf$cmt)
      .lst$nlmixrRowNums <- .data[.lst$filter, "nlmixrRowNums"]
      if (.lst$nnorm == 0L) {
        stop("need normal data for vpcSim (or use normRelated=FALSE)")
      }
      .data <- .data[.lst$filter, ]
    }
  }
  .si$events <- .data
  .si$thetaMat <- NULL
  .si$dfSub <- NULL
  .si$dfObs <- NULL
  .si$returnType <- "data.frame.TBS"
  .sim <- do.call(rxode2::rxSolve, .si)
  # now look for how many have missing values
  .w <- which(is.na(.sim$ipred))
  .nretry <- 0
  while (length(.w) > 0 && .nretry < nretry) {
    .w <- which(is.na(.sim$ipred))
    .simIds <- unique(.sim$sim.id[.w])
    .sim <- .sim[!(.sim$sim.id %in% .simIds), ]
    .sim$sim.id <- as.integer(factor(.sim$sim.id))
    .mx <- max(.sim$sim.id)
    .si$nsim <- n - .mx
    .sim2 <- do.call(rxode2::rxSolve, .si)
    .sim2$sim.id <- .sim2$sim.id + .mx
    .sim <- rbind(.sim, .sim2)
    .w <- which(is.na(.sim$ipred))
    .nretry <- .nretry + 1
  }
  if (.nretry != 0) {
    if (length(.w) == 0) {
      warning("'NA' values in vpc or npde simulation, re-simulated until all simulations were successful",
              call.=FALSE)
    } else {
      warning("'NA' values in vpc or npde simulation",
              call.=FALSE)
    }
  }
  if (pred) {
    .si$nsim <- n # restore for pred
    .si2 <- .si
    .si2$params <- c(
      .si$params, setNames(rep(0, dim(.si$omega)[1]), dimnames(.si$omega)[[2]]),
      setNames(rep(0, dim(.si$sigma)[1]), dimnames(.si$sigma)[[2]]))
    .si2$omega <- NULL
    .si2$sigma <- NULL
    .si2$returnType <- "data.frame"
    .si2$nStud <- 1
    .si2$nsim <- NULL
    assignInMyNamespace(".lastPredSimulationInfo", .si2)
    .sim2 <- do.call(rxode2::rxSolve, .si2)
    .sim$pred <- .sim2$sim
  }
  .sim <- vpcNameDataCmts(object, .sim)
  .cls <- c("nlmixr2vpcSim", class(.sim))
  .fit <- object
  .cls0 <- c("rxHidden", class(.fit))
  attr(.cls0, ".foceiEnv") <- attr(class(.fit), ".foceiEnv")
  class(.fit) <- .cls0
  attr(.cls, "fit") <- .fit
  class(.sim) <- .cls
  return(.sim)
}
#' Name the data and compartments
#'
#' @param object nlmixr2 fit object
#' @param data dataset to name `dvid` and `cmt` columns to correspond with the model
#' @return Updated object/data
#' @author Matthew L. Fidler
#' @keywords internal
#' @export
vpcNameDataCmts <- function(object, data) {
  assignInMyNamespace(".finalUiCompressed", FALSE)
  on.exit(assignInMyNamespace(".finalUiCompressed", TRUE))
  .wdvid <- which(tolower(names(data)) == "dvid")
  .wcmt <- which(tolower(names(data)) == "cmt")
  .info <- get("predDf", object$ui)
  if (is.null(.info)) {
    return(invisible(data))
  }
  if (length(.info$cond) == 1L) {
    return(invisible(data))
  }
  .state0 <- rxode2::rxModelVars(object$ui)$state
  .maxCmt <- max(.info$cmt)
  .cmtF <- character(max(length(.state0), .maxCmt))
  .dvidF <- character(length(.info$cmt))
  for (.i in seq_along(.state0)) {
    .cmtF[.i] <- .state0[.i]
  }
  for (.i in seq_along(.info$cmt)) {
    .cmtF[.info$cmt[.i]] <- paste(.info$cond[.i])
    .dvidF[.i] <- paste(.info$cond[.i])
  }
  if (length(.wdvid) == 1) {
    if (inherits(data[[.wdvid]], "numeric")) {
      data[[.wdvid]] <- as.integer(data[[.wdvid]])
    }
    if (!inherits(data[[.wdvid]], "factor") &&
          inherits(data[[.wdvid]], "integer")) {
      .tmp <- data[[.wdvid]]
      attr(.tmp, "levels") <- .dvidF
      attr(.tmp, "class") <- "factor"
      data[[.wdvid]] <- .tmp
    } else if (inherits(data[[.wdvid]], "character")) {
      data[[.wdvid]] <- factor(data[[.wdvid]], .dvidF)
    }
  }
  if (length(.wcmt) == 1) {
    if (inherits(data[[.wcmt]], "numeric")) {
      data[[.wcmt]] <- as.integer(data[[.wcmt]])
    }
    if (!inherits(data[[.wcmt]], "factor") &&
          inherits(data[[.wcmt]], "integer")) {
      .tmp <- data[[.wcmt]]
      attr(.tmp, "levels") <- .cmtF
      attr(.tmp, "class") <- "factor"
      data[[.wcmt]] <- .tmp
    } else if (inherits(data[[.wcmt]], "character")) {
      data[[.wcmt]] <- factor(data[[.wcmt]], .cmtF)
    }
  }
  data
}

#' Expand a VPC simulation
#'
#'
#' @param object nlmixr fit object
#' @param sim vpc simulation object
#' @param extra extra data from original fit to add
#' @param fullData is the full data (possibly modified); This is used
#'   for the vpc tad calculation
#' @return Expanded data frame with extra pieces added
#' @author Matthew L. Fidler
#' @export
#' @keywords internal
vpcSimExpand <- function(object, sim, extra, fullData=NULL) {
  if (is.null(extra)) return(sim)
  if (is.null(fullData)) {
    .fullData <- object$origData
  } else {
    .fullData <- fullData
  }
  .fullData$nlmixrRowNums <- seq_along(.fullData[, 1])
  .extra <- extra[extra %in% names(.fullData)]
  .extra <- extra[!(extra %in% names(sim))]
  if (length(.extra) == 0) return(sim)
  .wid <- which(tolower(names(.fullData)) == "id")
  names(.fullData)[.wid] <- "ID"
  .sim <- sim
  .wid <- which(tolower(names(.sim)) == "id")
  names(.sim)[.wid] <- "ID"
  .ret <- merge(.fullData, .sim, by=c("ID", "nlmixrRowNums"))
  .w <- which(names(.ret) == "nlmixrRowNums")
  vpcNameDataCmts(object, .ret[, -.w])
}

#' Get the least prediction simulation information for VPC
#'
#' @return The last prediction simulation from the `vpcSim` function (data.frame)
#'
#' @export
#' @keywords internal
.nlmixr2estLastPredSimulationInfo <- function() {
  .lastPredSimulationInfo
}