1 files changed, 137 insertions, 0 deletions
diff --git a/R/illparms.R b/R/illparms.R
new file mode 100644
index 00000000..931d8f05
--- /dev/null
+++ b/R/illparms.R
@@ -0,0 +1,137 @@
+#' Method to get the names of ill-defined parameters
+#'
+#' The method for generalised nonlinear regression fits as obtained
+#' with [mkinfit] and [mmkin] checks if the degradation parameters
+#' pass the Wald test (in degradation kinetics often simply called t-test) for
+#' significant difference from zero. For this test, the parameterisation
+#' without parameter transformations is used.
+#'
+#' The method for hierarchical model fits, also known as nonlinear
+#' mixed-effects model fits as obtained with [saem] and [mhmkin]
+#' checks if any of the confidence intervals for the random
+#' effects expressed as standard deviations include zero, and if
+#' the confidence intervals for the error model parameters include
+#' zero.
+#'
+#' @param object The object to investigate
+#' @param x The object to be printed
+#' @param conf.level The confidence level for checking p values
+#' @param \dots For potential future extensions
+#' @param random For hierarchical fits, should random effects be tested?
+#' @param errmod For hierarchical fits, should error model parameters be
+#' tested?
+#' @return For [mkinfit] or [saem] objects, a character vector of parameter
+#' names. For [mmkin] or [mhmkin] objects, a matrix like object of class
+#' 'illparms.mmkin' or 'illparms.mhmkin'. The latter objects have a suitable
+#' printing method.
+#' @export
+illparms <- function(object, ...)
+{
+  UseMethod("illparms", object)
+}
+
+#' @rdname illparms
+#' @export
+#' @examples
+#' fit <- mkinfit("FOMC", FOCUS_2006_A, quiet = TRUE)
+#' illparms(fit)
+illparms.mkinfit <- function(object, conf.level = 0.95, ...) {
+  p_values <- suppressWarnings(summary(object)$bpar[, "Pr(>t)"])
+  na <- is.na(p_values)
+  failed <- p_values > 1 - conf.level
+  names(parms(object))[na | failed]
+}
+
+#' @rdname illparms
+#' @export
+#' @examples
+#' \dontrun{
+#' fits <- mmkin(
+#'   c("SFO", "FOMC"),
+#'   list("FOCUS A" = FOCUS_2006_A,
+#'        "FOCUS C" = FOCUS_2006_C),
+#'   quiet = TRUE)
+#' illparms(fits)
+#' }
+illparms.mmkin <- function(object, conf.level = 0.95, ...) {
+  result <- lapply(object,
+    function(fit) {
+      if (inherits(fit, "try-error")) return("E")
+      ill <- illparms(fit, conf.level = conf.level)
+      if (length(ill) > 0) {
+        return(paste(ill, collapse = ", "))
+      } else {
+        return("")
+      }
+    })
+  result <- unlist(result)
+  dim(result) <- dim(object)
+  dimnames(result) <- dimnames(object)
+  class(result) <- "illparms.mmkin"
+  return(result)
+}
+
+#' @rdname illparms
+#' @export
+print.illparms.mmkin <- function(x, ...) {
+  class(x) <- NULL
+  print(x, quote = FALSE)
+}
+
+#' @rdname illparms
+#' @export
+illparms.saem.mmkin <- function(object, conf.level = 0.95, random = TRUE, errmod = TRUE, ...) {
+  if (inherits(object, "try-error")) {
+    return(NA)
+  } else {
+    ints <- intervals(object, conf.level = conf.level)
+    failed <- NULL
+    if (random) {
+      failed_random <- ints$random[, "lower"] < 0
+      failed <- c(failed, names(which(failed_random)))
+    }
+    if (errmod) {
+      failed_errmod <- ints$errmod[, "lower"] < 0 & ints$errmod[, "est."] > 0
+      failed <- c(failed, names(which(failed_errmod)))
+    }
+  }
+  return(failed)
+}
+
+#' @rdname illparms
+#' @export
+illparms.mhmkin <- function(object, conf.level = 0.95, random = TRUE, errmod = TRUE, ...) {
+  if (inherits(object[[1]], "saem.mmkin")) {
+    check_failed <- function(x) if (inherits(x$so, "try-error")) TRUE else FALSE
+  }
+  result <- lapply(object,
+    function(fit) {
+      if (check_failed(fit)) {
+        return("E")
+      } else {
+        if (!is.null(fit$FIM_failed) &&
+          "random effects and error model parameters" %in% fit$FIM_failed) {
+          return("NA")
+        } else {
+          ill <- illparms(fit, conf.level = conf.level, random = random, errmod = errmod)
+          if (length(ill) > 0) {
+            return(paste(ill, collapse = ", "))
+          } else {
+            return("")
+          }
+        }
+      }
+    })
+  result <- unlist(result)
+  dim(result) <- dim(object)
+  dimnames(result) <- dimnames(object)
+  class(result) <- "illparms.mhmkin"
+  return(result)
+}
+
+#' @rdname illparms
+#' @export
+print.illparms.mhmkin <- function(x, ...) {
+  class(x) <- NULL
+  print(x, quote = FALSE)
+}