R/nmlist.R

In kylebaron/nmlist:

extract_matrix <- function(x,prefix) {
  nrow <- ncol <- length(x)
  m <- matrix(0, nrow = nrow, ncol = ncol)
  for(row in seq_along(x)) {
    for(col in seq_along(x[[row]])) {
      value <- as.numeric(x[[row]][[col]])
      m[row,col] <- value
      m[col,row] <- value
    }
  }
  #dimnames(m) <- list(paste0(prefix, seq_len(nrow)), NULL)
  m
}

extract_theta <- function(x) {
  if(is.null(x)) return(list())
  labels <- sapply(x, attr, which = "name")
  values <- as.numeric(unlist(x, use.names=FALSE))
  setNames(as.list(values), paste0("THETA",labels))
}

extract_shrink <- function(x) {
  if(is.null(x)) return(list())
  x <- x$row
  out <- vector("numeric", length(x))
  name <- vector("character", length(x))
  for(i in seq_along(x)) {
    out[i] <- as.numeric(x[[i]])
    name[i] <- attr(x[[i]], which = "cname")
  }
  setNames(out,name)
}

#' Get NONMEM run results
#'
#' @param run the run number
#' @param project the project directory
#' @param rundir the run directory
#' @param raw if \code{TRUE}, the raw list is returned
#'
#' @details
#' \itemize{
#'   \item status: misc run information
#'   \item theta: list of theta values
#'   \item omega: omega matrix
#'   \item sigma: sigma matrix
#'   \item par: theta estimates with standard errors
#'   \item etashk: eta shrinkage values (percent)
#'   \item epsshk: epsilon shrinkage values (percent)
#'   \item etabar: mean eta
#'   \item etabarp: etabar p-value
#' }
#'
#' @examples
#'
#' project <- system.file("nonmem", package="nmlist")
#'
#' x <- nmlist(1005, project)
#'
#' names(x)
#'
#' x$theta
#' x$omega
#' x$cov
#'
#' @export
nmlist <- function(run, project, rundir = file.path(project,run), raw = FALSE) {
  xmlfile <- file.path(rundir, paste0(run, ".xml"))
  xml <- xml2::read_xml(xmlfile)
  x <- xml2::as_list(xml)
  if(raw) return(x)
  dt <- x$stop_datetime[[1]]
  prob <- x$nonmem$problem$problem_information[[1]]
  prob <- strsplit(prob, "\n")[[1]]
  recs <- grep("TOT. NO. OF", prob)
  drecs <- grep("NO. OF DATA RECS", prob)
  recs <- prob[c(recs,drecs)]
  recs <- strsplit(recs, ":", fixed = TRUE)
  recs <- data.frame(name = sapply(recs, "[[", 1L),
                     value = as.integer(sapply(recs, "[[", 2L)))
  x <- x$nonmem$problem$estimation
  th <- extract_theta(x$theta)
  thse <- extract_theta(x$thetase)
  par <- data.frame(THETA = as.numeric(th))
  se <- as.numeric(thse)
  if(length(se)==nrow(par)) par$SE <- se
  om <- extract_matrix(x$omega, "OMEGA")
  omc <- extract_matrix(x$omegac, "OMEGA")
  sg <- extract_matrix(x$sigma,"SIGMA")
  sgc <- extract_matrix(x$sigmac, "SIGMA")

  covar <- extract_matrix(x$covariance, "COVAR")
  if(nrow(covar) > 0) {
    covar <- covar[seq_along(th),seq_along(th)]
  }
  corr <- extract_matrix(x$correlation, "CORR")
  if(nrow(covar) > 0) {
    corr <- corr[seq_along(th),seq_along(th)]
  }
  cor_flag <- any(corr[lower.tri(corr)] > 0.9)
  etash <- extract_shrink(x$etashrink)
  epssh <- extract_shrink(x$epsshrink)
  etab <- extract_shrink(x$etabar)
  etabp <- extract_shrink(x$etabarpvalue)
  eigen <- extract_theta(x$eigenvalues)
  eigen <- as.numeric(unname(eigen))
  number <- NA
  if(length(eigen) > 0) {
     number <- max(eigen)/min(eigen)
  }

  status <- list(OFV = as.numeric(x$final_objective_function),
                 COV = as.numeric(x$covariance_status),
                 TERM = as.numeric(x$termination_status),
                 TIME = paste0(x$estimation_elapsed_time,"/",x$covariance_elapsed_time),
                 COR_FLAG = cor_flag, RECS = recs, DATE = dt)

  x <- list(status = status,
            theta = th, omega = om, sigma = sg,
            par = par,
            etashk = etash, epsshk = epssh,
            etabar = etab, etabarp = etabp,
            eigen = eigen, condition = number,
            cov = covar, cor = corr)

  structure(x, class="nmlist")

}

#' @export
print.nmlist <- function(x,...) {
  print(x$status)
}