R/npde.R

In tidyvpc: VPC Percentiles and Prediction Intervals

#' Normalized Prediction Distribution Errors
#'
#' @param o A \code{tidyvpcobj}.
#' @param id A vector of IDs. Used to associate observations (\code{y}) that
#' originate from the same individual. Evaluated in the \code{data.frame}
#' \code{data}.
#' @param data A \code{data.frame}.
#' @param smooth Should a uniform random perturbation be used to smooth the pd/pde values?
#' @param ... Additional arguments.
#'
#' @references
#'
#' Brendel, K., Comets, E., Laffont, C., Laveille, C. & Mentrée, F. Metrics
#' for external model evaluation with an application to the population
#' pharmacokinetics of gliclazide.  Pharm. Res. (2006) 23(9), 2036–2049.
#'
#' Nguyen, T.H.T., et al. Model evaluation of continuous data pharmacometric
#' models: metrics and graphics.  CPT Pharmacometrics Syst. Pharmacol. (2017)
#' 6(2), 87–109; doi:10.1002/psp4.12161.
#'
#' @examples
#' \donttest{
#' require(magrittr)
#' require(ggplot2)
#' 
#' obs <- obs_data[MDV==0]
#' sim <- sim_data[MDV==0]
#'
#' npde <- observed(obs, x=NULL, y=DV) %>%
#'     simulated(sim, y=DV) %>%
#'     npde(id=ID)
#'
#' vpc <- observed(npde$npdeobs, x=epred, y=npde) %>%
#'     simulated(npde$npdesim, y=npde) %>%
#'     binning("eqcut", nbins=10) %>%
#'     vpcstats()
#'
#' plot(vpc) + 
#' labs(x="Simulation-based Population Prediction", y="Normalized Prediction Distribution Error")
#' }
#'
#' @export
npde <- function(o, ...) UseMethod("npde")

#' @rdname npde
#' @export
npde.tidyvpcobj <- function(o, id, data=o$data, smooth=FALSE, ...) {

  obs <- rn <- y <- iter <- NULL
  
  if (missing(id)) {
    id <- o$id
  } else {
    id <- rlang::eval_tidy(rlang::enquo(id), data)
  }
  if (is.null(id)) {
    stop("No id specified")
  }

  niter <- nrow(o$sim)/nrow(o$obs)

  calc.npde <- function(dv, iter, smooth=FALSE) {
    yobs <- dv[iter==0]
    nobs <- length(yobs)
    niter <- max(iter) # = length(dv)/nobs - 1
    ysim <- matrix(dv[iter != 0], nrow=nobs, ncol=niter)

    # -- Non-decorrelated --

    pd.obs <- apply(ysim < yobs, 1, mean)

    # Handle extremes
    if (isTRUE(smooth)) {
      u <- runif(nobs, 0, 1/niter)
      pd.obs[pd.obs == 1] <- rep(1 - (1/niter), len=sum(pd.obs == 1))
      pd.obs <- pd.obs + u
    } else {
      pd.obs[pd.obs == 0] <- 1/(2*niter)
      pd.obs[pd.obs == 1] <- 1 - 1/(2*niter)
    }

    # Normalize
    npd.obs <- qnorm(pd.obs)

    # Compute for simulated data
    r <- apply(ysim, 1, rank)
    pd.sim <- r/(niter + 1)   # +1 prevents Inf after normalization
    npd.sim <- qnorm(pd.sim)

    # -- Decorrelated --

    epred <- apply(ysim, 1, mean)
    ec <- cov(t(ysim))
    ec.inv.chol <- chol(solve(ec))  # In future, should we consider different decorrelation methods?

    eres <- yobs - epred
    esres <- ysim - epred

    ewres <- drop(ec.inv.chol %*% eres)
    eswres <- ec.inv.chol %*% esres

    pde.obs <- apply(eswres < ewres, 1, mean)

    # Handle extremes
    if (isTRUE(smooth)) {
      # Can use the same uniform as for pd, no need to call runif again
      pde.obs[pde.obs == 1] <- rep(1 - (1/niter), len=sum(pde.obs == 1))
      pde.obs <- pde.obs + u
    } else {
      pde.obs[pde.obs == 0] <- 1/(2*niter)
      pde.obs[pde.obs == 1] <- 1 - 1/(2*niter)
    }

    # Normalize
    npde.obs <- qnorm(pde.obs)

    # Compute for simulated data
    r <- apply(eswres, 1, rank)
    pde.sim <- r/(niter + 1)   # +1 prevents Inf after normalization
    npde.sim <- qnorm(pde.sim)

    data.table(
      iter=iter,
      epred=epred,
      eres=c(eres, esres),
      ewres=c(ewres, eswres),
      npd=c(npd.obs, npd.sim),
      npde=c(npde.obs, npde.sim))
  }

  obssim <- data.table(
    y    = c(o$obs$y, o$sim$y),
    id   = rep(id, length.out=(niter + 1)*nrow(o$obs)),
    iter = rep(0:niter, each=nrow(o$obs)))

  obssim <- obssim[, rn := (1:.N)][, cbind(rn, calc.npde(y, iter, smooth=smooth)), by=id][order(rn)][, rn := NULL]

  npdeobs <- obssim[iter==0]
  npdesim <- obssim[iter!=0]

  update(o, npdeobs=npdeobs, npdesim=npdesim)
}

# vim: ts=2 sw=2 et