Merge branch 'main' into dev

10 files changed, 88 insertions, 43 deletions

diff --git a/R/create_deg_func.R b/R/create_deg_func.R
index 5794c65c..d3e11f78 100644
--- a/R/create_deg_func.R
+++ b/R/create_deg_func.R
@@ -61,16 +61,21 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) {
     ),
   ")")
 
-  if (length(obs_vars) >= 2) {
-    supported <- FALSE # except for the following cases
+  if (length(obs_vars) >= 2) supported <- FALSE
+  # Except for the following cases:
+
+  if (length(obs_vars) == 2) {
     n1 <- obs_vars[1]
     n2 <- obs_vars[2]
     n10 <- paste0("odeini['", parent, "']")
     n20 <- paste0("odeini['", n2, "']")
 
     # sfo_sfo
-    if (all(spec[[1]]$sink == FALSE, length(obs_vars) == 2,
-        spec[[1]]$type == "SFO", spec[[2]]$type == "SFO")) {
+    if (all(
+        spec[[1]]$sink == FALSE,
+        spec[[1]]$type == "SFO",
+        spec[[2]]$type == "SFO",
+        is.null(spec[[2]]$to))) {
       supported <- TRUE
       k1 <- k_parent
       k2 <- paste0("k_", n2, if(min_ff) "_sink" else "")
@@ -80,8 +85,12 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) {
     }
 
     # sfo_f12_sfo
-    if (all(use_of_ff == "max", spec[[1]]$sink == TRUE, length(obs_vars) == 2,
-        spec[[1]]$type == "SFO", spec[[2]]$type == "SFO")) {
+    if (all(
+        use_of_ff == "max",
+        spec[[1]]$sink == TRUE,
+        spec[[1]]$type == "SFO",
+        spec[[2]]$type == "SFO",
+        is.null(spec[[2]]$to))) {
       supported <- TRUE
       k1 <- k_parent
       k2 <- paste0("k_", n2)
@@ -92,8 +101,12 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) {
     }
 
     # sfo_k120_sfo
-    if (all(use_of_ff == "min", spec[[1]]$sink == TRUE, length(obs_vars) == 2,
-        spec[[1]]$type == "SFO", spec[[2]]$type == "SFO")) {
+    if (all(
+        use_of_ff == "min",
+        spec[[1]]$sink == TRUE,
+        spec[[1]]$type == "SFO",
+        spec[[2]]$type == "SFO",
+        is.null(spec[[2]]$to))) {
       supported <- TRUE
       k12 <- paste0("k_", n1, "_", n2)
       k10 <- paste0("k_", n1, "_sink")
@@ -104,8 +117,12 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) {
     }
 
     # dfop_f12_sfo
-    if (all(use_of_ff == "max", spec[[1]]$sink == TRUE, length(obs_vars) == 2,
-        spec[[1]]$type == "DFOP", spec[[2]]$type == "SFO")) {
+    if (all(
+        use_of_ff == "max",
+        spec[[1]]$sink == TRUE,
+        spec[[1]]$type == "DFOP",
+        spec[[2]]$type == "SFO",
+        is.null(spec[[2]]$to))) {
       supported <- TRUE
       f12 <- paste0("f_", n1, "_to_", n2)
       k2 <- paste0("k_", n2)
@@ -119,7 +136,6 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) {
 
   }
 
-
   if (supported) {
     deg_func <- function(observed, odeini, odeparms) {}
 
diff --git a/R/illparms.R b/R/illparms.R
index 68a6bff6..b4b37fbb 100644
--- a/R/illparms.R
+++ b/R/illparms.R
@@ -102,12 +102,14 @@ illparms.saem.mmkin <- function(object, conf.level = 0.95, random = TRUE, errmod
     ints <- intervals(object, conf.level = conf.level)
     ill_parms <- character(0)
     if (random) {
-      ill_parms_random <- ints$random[, "lower"] < 0
-      ill_parms <- c(ill_parms, names(which(ill_parms_random)))
+      ill_parms_random_i <- which(ints$random[, "lower"] < 0)
+      ill_parms_random <- rownames(ints$random)[ill_parms_random_i]
+      ill_parms <- c(ill_parms, ill_parms_random)
     }
     if (errmod) {
-      ill_parms_errmod <- ints$errmod[, "lower"] < 0 & ints$errmod[, "upper"] > 0
-      ill_parms <- c(ill_parms, names(which(ill_parms_errmod)))
+      ill_parms_errmod_i <- which(ints$errmod[, "lower"] < 0 & ints$errmod[, "upper"] > 0)
+      ill_parms_errmod <- rownames(ints$errmod)[ill_parms_errmod_i]
+      ill_parms <- c(ill_parms, ill_parms_errmod)
     }
     if (slopes) {
       if (is.null(object$so)) stop("Slope testing is only implemented for the saemix backend")
diff --git a/R/mhmkin.R b/R/mhmkin.R
index 6265a59e..14a7ac29 100644
--- a/R/mhmkin.R
+++ b/R/mhmkin.R
@@ -219,11 +219,22 @@ print.mhmkin <- function(x, ...) {
   print(status(x))
 }
 
+#' Check if fit within an mhmkin object failed
+#' @param x The object to be checked
+check_failed <- function(x) {
+  if (inherits(x, "try-error")) {
+    return(TRUE)
+  } else {
+    if (inherits(x$so, "try-error")) {
+      return(TRUE)
+    } else {
+      return(FALSE)
+    }
+  }
+}
+
 #' @export
 AIC.mhmkin <- function(object, ..., k = 2) {
-  if (inherits(object[[1]], "saem.mmkin")) {
-    check_failed <- function(x) if (inherits(x$so, "try-error")) TRUE else FALSE
-  }
   res <- sapply(object, function(x) {
     if (check_failed(x)) return(NA)
     else return(AIC(x$so, k = k))
@@ -235,9 +246,6 @@ AIC.mhmkin <- function(object, ..., k = 2) {
 
 #' @export
 BIC.mhmkin <- function(object, ...) {
-  if (inherits(object[[1]], "saem.mmkin")) {
-    check_failed <- function(x) if (inherits(x$so, "try-error")) TRUE else FALSE
-  }
   res <- sapply(object, function(x) {
     if (check_failed(x)) return(NA)
     else return(BIC(x$so))
@@ -280,7 +288,13 @@ anova.mhmkin <- function(object, ...,
   if (identical(model.names, "auto")) {
     model.names <- outer(rownames(object), colnames(object), paste)
   }
-  rlang::inject(anova(!!!(object), method = method, test = test,
+  failed_index <- which(sapply(object, check_failed), arr.ind = TRUE)
+  if (length(failed_index > 0)) {
+    rlang::inject(anova(!!!(object[-failed_index]), method = method, test = test,
+      model.names = model.names[-failed_index]))
+  } else {
+    rlang::inject(anova(!!!(object), method = method, test = test,
       model.names = model.names))
+  }
 }
 
diff --git a/R/mkinfit.R b/R/mkinfit.R
index c851fddb..52053685 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -21,7 +21,8 @@ utils::globalVariables(c("name", "time", "value"))
 #'   "FOMC", "DFOP", "HS", "SFORB", "IORE"). If a shorthand name is given, a
 #'   parent only degradation model is generated for the variable with the
 #'   highest value in \code{observed}.
-#' @param observed A dataframe with the observed data.  The first column called
+#' @param observed A dataframe or an object coercible to a dataframe
+#'   (e.g. a \code{tibble}) with the observed data.  The first column called
 #'   "name" must contain the name of the observed variable for each data point.
 #'   The second column must contain the times of observation, named "time".
 #'   The third column must be named "value" and contain the observed values.
@@ -292,6 +293,9 @@ mkinfit <- function(mkinmod, observed,
   # Get the names of observed variables
   obs_vars <- names(mkinmod$spec)
 
+  # Coerce observed data to a dataframe
+  observed <- as.data.frame(observed)
+
   # Subset observed data with names of observed data in the model and remove NA values
   observed <- subset(observed, name %in% obs_vars)
   observed <- subset(observed, !is.na(value))
diff --git a/R/mkinpredict.R b/R/mkinpredict.R
index 60456fb2..4c6d7862 100644
--- a/R/mkinpredict.R
+++ b/R/mkinpredict.R
@@ -22,7 +22,7 @@
 #' variables. The third possibility "eigen" is fast in comparison to uncompiled
 #' ODE models, but not applicable to some models, e.g. using FOMC for the
 #' parent compound.
-#' @param method.ode The solution method passed via [mkinpredict] to [ode]] in
+#' @param method.ode The solution method passed via [mkinpredict] to `deSolve::ode()` in
 #' case the solution type is "deSolve" and we are not using compiled code.
 #' When using compiled code, only lsoda is supported.
 #' @param use_compiled If set to \code{FALSE}, no compiled version of the
@@ -36,7 +36,7 @@
 #'   the observed variables (default) or for all state variables (if set to
 #'   FALSE). Setting this to FALSE has no effect for analytical solutions,
 #'   as these always return mapped output.
-#' @param na_stop Should it be an error if [ode] returns NaN values
+#' @param na_stop Should it be an error if `deSolve::ode()` returns NaN values
 #' @param \dots Further arguments passed to the ode solver in case such a
 #'   solver is used.
 #' @return A matrix with the numeric solution in wide format
diff --git a/R/multistart.R b/R/multistart.R
index aeea2d81..4ba82a43 100644
--- a/R/multistart.R
+++ b/R/multistart.R
@@ -40,7 +40,7 @@
 #' f_mmkin <- mmkin("DFOP", dmta_ds, error_model = "tc", cores = 7, quiet = TRUE)
 #' f_saem_full <- saem(f_mmkin)
 #' f_saem_full_multi <- multistart(f_saem_full, n = 16, cores = 16)
-#' parplot(f_saem_full_multi, lpos = "topleft")
+#' parplot(f_saem_full_multi, lpos = "topleft", las = 2)
 #' illparms(f_saem_full)
 #'
 #' f_saem_reduced <- update(f_saem_full, no_random_effect = "log_k2")
@@ -50,7 +50,7 @@
 #' library(parallel)
 #' cl <- makePSOCKcluster(12)
 #' f_saem_reduced_multi <- multistart(f_saem_reduced, n = 16, cluster = cl)
-#' parplot(f_saem_reduced_multi, lpos = "topright", ylim = c(0.5, 2))
+#' parplot(f_saem_reduced_multi, lpos = "topright", ylim = c(0.5, 2), las = 2)
 #' stopCluster(cl)
 #' }
 multistart <- function(object, n = 50,
diff --git a/R/nlme.R b/R/nlme.R
index 6b2d06d0..de82d99c 100644
--- a/R/nlme.R
+++ b/R/nlme.R
@@ -126,7 +126,7 @@ nlme_function <- function(object) {
 
 #' @rdname nlme
 #' @importFrom rlang !!!
-#' @return A \code{\link{groupedData}} object
+#' @return A `nlme::groupedData` object
 #' @export
 nlme_data <- function(object) {
   if (nrow(object) > 1) stop("Only row objects allowed")
diff --git a/R/parplot.R b/R/parplot.R
index 3da4b51a..a33112a5 100644
--- a/R/parplot.R
+++ b/R/parplot.R
@@ -35,9 +35,6 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, llquant = NA,
   scale = c("best", "median"),
   lpos = "bottomleft", main = "", ...)
 {
-  oldpar <- par(no.readonly = TRUE)
-  on.exit(par(oldpar, no.readonly = TRUE))
-
   orig <- attr(object, "orig")
   orig_parms <- parms(orig)
   start_degparms <- orig$mean_dp_start
@@ -59,11 +56,10 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, llquant = NA,
   selected <- which(ll > llmin)
   selected_parms <- all_parms[selected, ]
 
-  par(las = 1)
   if (orig$transformations == "mkin") {
     degparm_names_transformed <- names(start_degparms)
     degparm_index <- which(names(orig_parms) %in% degparm_names_transformed)
-    orig_parms[degparm_names_transformed] <- backtransform_odeparms(
+    orig_degparms <- backtransform_odeparms(
       orig_parms[degparm_names_transformed],
       orig$mmkin[[1]]$mkinmod,
       transform_rates = orig$mmkin[[1]]$transform_rates,
@@ -74,14 +70,17 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, llquant = NA,
       transform_fractions = orig$mmkin[[1]]$transform_fractions)
     degparm_names <- names(start_degparms)
 
-    names(orig_parms) <- c(degparm_names, names(orig_parms[-degparm_index]))
+    orig_parms_back <- orig_parms
+    orig_parms_back[degparm_index] <- orig_degparms
+    names(orig_parms_back)[degparm_index] <- degparm_names
+    orig_parms <- orig_parms_back
 
     selected_parms[, degparm_names_transformed] <-
       t(apply(selected_parms[, degparm_names_transformed], 1, backtransform_odeparms,
       orig$mmkin[[1]]$mkinmod,
       transform_rates = orig$mmkin[[1]]$transform_rates,
       transform_fractions = orig$mmkin[[1]]$transform_fractions))
-    colnames(selected_parms)[1:length(degparm_names)] <- degparm_names
+    colnames(selected_parms)[degparm_index] <- degparm_names
   }
 
   start_errparms <- orig$so@model@error.init
@@ -99,6 +98,12 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, llquant = NA,
 
   # Boxplots of all scaled parameters
   selected_scaled_parms <- t(apply(selected_parms, 1, function(x) x / parm_scale))
+  i_negative <- selected_scaled_parms <= 0
+  parms_with_negative_scaled_values <- paste(names(which(apply(i_negative, 2, any))), collapse = ", ")
+  if (any(i_negative)) {
+    warning("There are negative values for ", parms_with_negative_scaled_values, " which are set to NA for plotting")
+  }
+  selected_scaled_parms[i_negative] <- NA
   boxplot(selected_scaled_parms, log = "y", main = main, ,
     ylab = "Normalised parameters", ...)
 
diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R
index d6c3d0de..f05f1110 100644
--- a/R/plot.mixed.mmkin.R
+++ b/R/plot.mixed.mmkin.R
@@ -93,7 +93,7 @@ plot.mixed.mmkin <- function(x,
   nrow.legend = ceiling((length(i) + 1) / ncol.legend),
   rel.height.legend = 0.02 + 0.07 * nrow.legend,
   rel.height.bottom = 1.1,
-  pch_ds = 1:length(i),
+  pch_ds = c(1:25, 33, 35:38, 40:41, 47:57, 60:90)[1:length(i)],
   col_ds = pch_ds + 1,
   lty_ds = col_ds,
   frame = TRUE, ...
diff --git a/R/status.R b/R/status.R
index 8bcd3a16..f9d79e7d 100644
--- a/R/status.R
+++ b/R/status.R
@@ -74,15 +74,19 @@ print.status.mmkin <- function(x, ...) {
 status.mhmkin <- function(object, ...) {
   if (inherits(object[[1]], "saem.mmkin")) {
     test_func <- function(fit) {
-      if (inherits(fit$so, "try-error")) {
-        return("E")
+      if (inherits(fit, "try-error")) {
+          return("E")
       } else {
-        if (!is.null(fit$FIM_failed)) {
-          return_values <- c("fixed effects" = "Fth",
-            "random effects and error model parameters" = "FO")
-          return(paste(return_values[fit$FIM_failed], collapse = ", "))
+        if (inherits(fit$so, "try-error")) {
+          return("E")
         } else {
-          return("OK")
+          if (!is.null(fit$FIM_failed)) {
+            return_values <- c("fixed effects" = "Fth",
+              "random effects and error model parameters" = "FO")
+            return(paste(return_values[fit$FIM_failed], collapse = ", "))
+          } else {
+            return("OK")
+          }
         }
       }
     }