Merge pull request #323 from nlmixr2/311-condition-numbers

311 condition numbers
nlmixr2 · Mar 10, 2023 · d4a11da · d4a11da
2 parents 5b63a4d + d2ed943
commit d4a11da
Show file tree

Hide file tree

Showing 10 changed files with 127 additions and 25 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -1,8 +1,14 @@
 # nlmixr2est (development version)
 
+- Breaking change, now calculate condition number based on covariance
+  and correlation, the names have changed to be more explicit.
+  `conditionNumber` changed to `conditionNumberCov` and a new metric
+  `conditionNumberCor` has been added.
+
 - A bug in boundary value detection prevented automatic covariance calculation
   with FOCEi estimation (#318)
 
+
 - Fix `vpcSim` so that it will be a bit more robust when it is
   difficult to simulate.
 

diff --git a/R/assert.R b/R/assert.R
@@ -60,9 +60,9 @@ assertNlmixrObjDataFrameRow <- function(df, allowNa=FALSE) {
          paste(.diff, collapse="', '"), "'",
          call.=FALSE)
   }
-  .w <- which(.name == "Condition Number")
+  .w <- which(.name == "Condition#(Cov)")
   if (length(.w) == 1) {
-    .cn <- df[["Condition Number"]]
+    .cn <- df[["Condition#(Cov)"]]
     if (inherits(.cn, "numeric")) {
       if (!is.na(.cn)) {
         .cn <- setNames(.cn, NULL)
@@ -75,6 +75,21 @@ assertNlmixrObjDataFrameRow <- function(df, allowNa=FALSE) {
   } else {
     .cn <- NA
   }
+  .w <- which(.name == "Condition#(Cor)")
+  if (length(.w) == 1) {
+    .cnr <- df[["Condition#(Cor)"]]
+    if (inherits(.cnr, "numeric")) {
+      if (!is.na(.cnr)) {
+        .cnr <- setNames(.cnr, NULL)
+      } else {
+        .cnr <- NA
+      }
+    } else {
+      .cnr <- NA
+    }
+  } else {
+    .cnr <- NA
+  }
   .df1 <- df[, .needed]
   if (length(.df1[, 1]) == 0) {
     stop("missing objective function in objective function data frame", call.=FALSE)
@@ -87,5 +102,5 @@ assertNlmixrObjDataFrameRow <- function(df, allowNa=FALSE) {
       checkmate::assertNumeric(df[[x]], len=1, any.missing=FALSE, .var.name=x)
     })
   }
-  list(.df1, .cn)
+  list(.df1, .cn, .cnr)
 }
diff --git a/R/broom.R b/R/broom.R
@@ -254,7 +254,8 @@ glance.nlmixr2FitCore <- function(x, ...) {
   .aic <- AIC(x) ## To calculate AIC if needed
   .df <- x$objDf[which(tolower(rownames(x$objDf)) == tolower(x$ofvType)), ]
   names(.df) <- gsub("Log-likelihood", "logLik", names(.df))
-  names(.df) <- gsub("Condition Number", "conditionNumber", names(.df))
+  names(.df) <- gsub("Condition#(Cor)", "conditionNumberCor", names(.df))
+  names(.df) <- gsub("Condition#(Cov)", "conditionNumberCov", names(.df))
   tibble::as_tibble(.df)
 }
 

diff --git a/R/complete.R b/R/complete.R
@@ -1,7 +1,8 @@
 
 .nmObjGetEnvInfo <- list(
   ui="rxode2 user interface",
-  conditionNumber="Condition Number",
+  conditionNumberCor="Condition Number (Correlation)",
+  conditionNumberCov="Condition Number (Covariance)",
   cov="Covariance of fixed effects",
   covMethod="Covariance Method for fixed effects",
   etaObf="ETAs and their individual objective function contribution (if applicable)",

diff --git a/R/cov.R b/R/cov.R
@@ -161,3 +161,12 @@ getVarCov.nlmixr2FitCore <- function(obj, ...) {
 
 ##' @export
 getVarCov.nlmixr2FitCoreSilent <- getVarCov.nlmixr2FitCore
+
+
+.cov2cor <- function(cov) {
+  .sd2 <- sqrt(diag(cov))
+  .cor <- stats::cov2cor(cov)
+  dimnames(.cor) <- dimnames(cov)
+  diag(.cor) <- .sd2
+  .cor
+}
diff --git a/R/focei.R b/R/focei.R
@@ -1750,13 +1750,22 @@ nlmixr2Est.posthoc <- function(env, ...) {
     ret <- attr(class(ret), ".foceiEnv")
   }
   .objDf1 <- get("objDf", ret)
-  if (any(names(.objDf1) == "Condition Number")) {
-    if (!any(names(objDf) == "Condition Number")) {
-      objDf[["Condition Number"]] <- NA_real_
+  if (any(names(.objDf1) == "Condition#(Cov)")) {
+    if (!any(names(objDf) == "Condition#(Cov)")) {
+      objDf[["Condition#(Cov)"]] <- NA_real_
     }
-  } else if (any(names(objDf) == "Condition Number")) {
-    if (!any(names(.objDf1) == "Condition Number")) {
-      .objDf1[["Condition Number"]] <- NA_real_
+  } else if (any(names(objDf) == "Condition#(Cov)")) {
+    if (!any(names(.objDf1) == "Condition#(Cov)")) {
+      .objDf1[["Condition#(Cov)"]] <- NA_real_
+    }
+  }
+  if (any(names(.objDf1) == "Condition#(Cor)")) {
+    if (!any(names(objDf) == "Condition#(Cor)")) {
+      objDf[["Condition#(Cor)"]] <- NA_real_
+    }
+  } else if (any(names(objDf) == "Condition#(Cor)")) {
+    if (!any(names(.objDf1) == "Condition#(Cor)")) {
+      .objDf1[["Condition#(Cor)"]] <- NA_real_
     }
   }
   assign("objDf", rbind(.objDf1, objDf), envir=ret)

diff --git a/R/nlmixrAddObjectiveFunction.R b/R/nlmixrAddObjectiveFunction.R
@@ -3,7 +3,7 @@
 #' @param fit nlmixr fit object
 #' @param objDf nlmixr objective function data frame which has column
 #'   names "OBJF", "AIC", "BIC", "Log-likelihood" and
-#'   "Condition Number"
+#'   "Condition#(Cov)" "Condition#(Cor)"
 #' @param type Objective Function Type
 #' @param etaObf Eta objective function table to add (with focei) to
 #'   give focei objective function
@@ -27,32 +27,48 @@ nlmixrAddObjectiveFunctionDataFrame <- function(fit, objDf, type, etaObf=NULL) {
     } else if (!is.na(.inRow[[2]])) {
       .cn <- .inRow[[2]]
     }
+    .cnr <- NA_real_
+    if (!is.na(.inRow2[[3]])) {
+      .cnr <- .inRow2[[3]]
+    } else if (!is.na(.inRow[[3]])) {
+      .cnr <- .inRow[[3]]
+    }
     if (is.na(.inRow2[[1]][[1]])) {
       # Here the original data frame is NA, that is the objective function has not been calculated
-      .tmp <- cbind(.inRow[[1]], data.frame("Condition Number"=.cn, check.names=FALSE))
+      .tmp <- cbind(.inRow[[1]], data.frame("Condition#(Cov)"=.cn, "Condition#(Cor)"=.cnr, check.names=FALSE))
       row.names(.tmp) <- type
       assign("objDf", .tmp, envir=fit$env)
       setOfv(fit, type)
     } else {
       if (any(.rownames == type)) stop("objective function '", type, "' already present", call.=FALSE)
       # Now the original data frame is not NA.
       .tmp <- rbind(.inRow[[1]], .inRow2[[1]])
-      .tmp[["Condition Number"]] <- .cn
+      .tmp[["Condition#(Cov)"]] <- .cn
+      .tmp[["Condition#(Cor)"]] <- .cnr
       row.names(.tmp) <- c(type, .rownames)
       assign("objDf", .tmp, envir=fit$env)
       setOfv(fit, type)
     }
   } else {
     if (any(.rownames == type)) stop("objective function '", type, "' already present", call.=FALSE)
     ## Now there is at least one interesting objective function
-    .cn <- .cur[["Condition Number"]][1]
+    .cn <- .cur[["Condition#(Cov)"]][1]
     if (is.null(.cn)) {
       .cn <- NA_real_
     }
     if (is.na(.cn) & !is.na(.inRow[[2]])) {
       .cn <- .inRow[[2]]
     }
-    .cur <- rbind(.cur[, names(.cur) != "Condition Number"], .inRow[[1]])
+    .cnr <- .cur[["Condition#(Cor)"]][1]
+    if (is.null(.cnr)) {
+      .cnr <- NA_real_
+    }
+    if (is.na(.cnr) & !is.na(.inRow[[3]])) {
+      .cnr <- .inRow[[3]]
+    }
+    .cur <- rbind(.cur[, !(names(.cur) %in% c("Condition#(Cor)", "Condition#(Cov)"))], .inRow[[1]])
+    .cur[["Condition#(Cov)"]] <- .cn
+    .cur[["Condition#(Cor)"]] <- .cnr
     row.names(.cur) <- c(.rownames, type)
     assign("objDf", .cur, envir=fit$env)
     setOfv(fit, type)

diff --git a/R/nmObjGet.R b/R/nmObjGet.R
@@ -494,8 +494,8 @@ nmObjGet.env <- function(x, ...) {
 nmObjGet.condition <- function(x, ...) {
   .env <- x[[1]]
   .objDf <- .env$objDf
-  if (any(names(.objDf) == "Condition Number")) {
-    .cn <- .objDf[, "Condition Number"]
+  if (any(names(.objDf) == "Condition#(Cor)")) {
+    .cn <- .objDf[, "Condition#(Cor)"]
     .cn <- .cn[!is.na(.cn)]
     return(.cn)
   }

diff --git a/src/inner.cpp b/src/inner.cpp
@@ -5877,13 +5877,44 @@ void foceiFinalizeTables(Environment e){
   }
   LogicalVector skipCov = e["skipCov"];
 
+  if (covExists) {
+    Function loadNamespace("loadNamespace", R_BaseNamespace);    
+    Environment nlmixr2 = loadNamespace("nlmixr2est");
+    Function getCor = nlmixr2[".cov2cor"];
+    e["fullCor"] = getCor(e["cov"]);
+    arma::mat cor = as<arma::mat>(e["fullCor"]);
+    cor.diag().ones();
+    arma::vec eigval;
+    arma::mat eigvec;
+    eig_sym(eigval, eigvec, cor);
+    e["eigenCor"] = eigval;
+    e["eigenVecCor"] = eigvec;
+    unsigned int k=0;
+    if (eigval.size() > 0){
+      double mx=std::fabs(eigval[0]), mn, cur;
+      mn=mx;
+      for (k = eigval.size(); k--;){
+        cur = std::fabs(eigval[k]);
+        if (cur > mx){
+          mx=cur;
+        }
+        if (cur < mn){
+          mn=cur;
+        }
+      }
+      e["conditionNumberCor"] = mx/mn;
+    } else {
+      e["conditionNumberCor"] = NA_REAL;
+    }
+  }
+
   if (covExists && op_focei.eigen){
     arma::vec eigval;
     arma::mat eigvec;
 
     eig_sym(eigval, eigvec, cov);
-    e["eigen"] = eigval;
-    e["eigenVec"] = eigvec;
+    e["eigenCov"] = eigval;
+    e["eigenVecCov"] = eigvec;
     unsigned int k=0;
     if (eigval.size() > 0){
       double mx=std::fabs(eigval[0]), mn, cur;
@@ -5897,9 +5928,9 @@ void foceiFinalizeTables(Environment e){
           mn=cur;
         }
       }
-      e["conditionNumber"] = mx/mn;
+      e["conditionNumberCov"] = mx/mn;
     } else {
-      e["conditionNumber"] = NA_REAL;
+      e["conditionNumberCov"] = NA_REAL;
     }
   }
   arma::vec se1;
@@ -6255,10 +6286,11 @@ void foceiFinalizeTables(Environment e){
     }
   }
   List objDf;
-  if (e.exists("conditionNumber")){
+  if (e.exists("conditionNumberCov")){
     objDf = List::create(_["OBJF"] = as<double>(e["objective"]), _["AIC"]=as<double>(e["AIC"]),
                          _["BIC"] = as<double>(e["BIC"]), _["Log-likelihood"]=as<double>(e["logLik"]),
-                         _["Condition Number"]=as<double>(e["conditionNumber"]));
+                         _["Condition#(Cov)"]=as<double>(e["conditionNumberCov"]),
+                         _["Condition#(Cor)"]=as<double>(e["conditionNumberCor"]));
   } else {
     objDf = List::create(_["OBJF"] = as<double>(e["objective"]), _["AIC"]=as<double>(e["AIC"]),
                          _["BIC"] = as<double>(e["BIC"]), _["Log-likelihood"]=as<double>(e["logLik"]));

diff --git a/tests/testthat/test-saem-aic.R b/tests/testthat/test-saem-aic.R
@@ -21,18 +21,31 @@ nmTest({
 
     fit <-
       suppressMessages(
-        nlmixr(one.cmt, theo_sd, est = "saem", control = list(calcTables = FALSE, print = 0))
+        nlmixr(one.cmt, theo_sd, est = "saem",
+               control = saemControl(calcTables = FALSE, print = 0, nBurn=10, nEm=10))
       )
 
     expect_s3_class(fit, "nlmixr2FitCore")
     expect_false(inherits(fit, "data.frame"))
     expect_false(inherits(fit, "nlmixrFitData"))
 
     expect_error(suppressMessages(setOfv(fit, "focei")), NA)
+    expect_true(any(names(fit$objDf) == "Condition#(Cov)"))
+    expect_true(any(names(fit$objDf) == "Condition#(Cor)"))
     expect_error(suppressMessages(setOfv(fit, "foce")), NA)
+    expect_true(any(names(fit$objDf) == "Condition#(Cov)"))
+    expect_true(any(names(fit$objDf) == "Condition#(Cor)"))
     expect_error(suppressMessages(setOfv(fit, "fo")), NA)
+    expect_true(any(names(fit$objDf) == "Condition#(Cov)"))
+    expect_true(any(names(fit$objDf) == "Condition#(Cor)"))
     expect_error(suppressMessages(setOfv(fit, "gauss3_1.6")), NA)
+    expect_true(any(names(fit$objDf) == "Condition#(Cov)"))
+    expect_true(any(names(fit$objDf) == "Condition#(Cor)"))
     expect_error(suppressMessages(setOfv(fit, "laplace2")), NA)
+    expect_true(any(names(fit$objDf) == "Condition#(Cov)"))
+    expect_true(any(names(fit$objDf) == "Condition#(Cor)"))
     expect_error(suppressMessages(setOfv(fit, "laplace3")), NA)
+    expect_true(any(names(fit$objDf) == "Condition#(Cov)"))
+    expect_true(any(names(fit$objDf) == "Condition#(Cor)"))
   })
 })