mixed.mmkin and test coverage

5 files changed, 188 insertions, 57 deletions

diff --git a/R/mixed.mmkin.R b/R/mixed.mmkin.R
new file mode 100644
index 00000000..6fe5130d
--- /dev/null
+++ b/R/mixed.mmkin.R
@@ -0,0 +1,101 @@
+#' Create a mixed effects model from an mmkin row object
+#'
+#' @param object An [mmkin] row object
+#' @param method The method to be used
+#' @param \dots Currently not used
+#' @examples
+#' sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
+#' n_biphasic <- 8
+#' err_1 = list(const = 1, prop = 0.07)
+#'
+#' DFOP_SFO <- mkinmod(
+#'   parent = mkinsub("DFOP", "m1"),
+#'   m1 = mkinsub("SFO"),
+#'   quiet = TRUE)
+#'
+#' set.seed(123456)
+#' log_sd <- 0.3
+#' syn_biphasic_parms <- as.matrix(data.frame(
+#'   k1 = rlnorm(n_biphasic, log(0.05), log_sd),
+#'   k2 = rlnorm(n_biphasic, log(0.01), log_sd),
+#'   g = plogis(rnorm(n_biphasic, 0, log_sd)),
+#'   f_parent_to_m1 = plogis(rnorm(n_biphasic, 0, log_sd)),
+#'   k_m1 = rlnorm(n_biphasic, log(0.002), log_sd)))
+#'
+#' ds_biphasic_mean <- lapply(1:n_biphasic,
+#'   function(i) {
+#'     mkinpredict(DFOP_SFO, syn_biphasic_parms[i, ],
+#'       c(parent = 100, m1 = 0), sampling_times)
+#'   }
+#' )
+#'
+#' set.seed(123456L)
+#' ds_biphasic <- lapply(ds_biphasic_mean, function(ds) {
+#'   add_err(ds,
+#'     sdfunc = function(value) sqrt(err_1$const^2 + value^2 * err_1$prop^2),
+#'     n = 1, secondary = "m1")[[1]]
+#' })
+#'
+#' \dontrun{
+#' f_mmkin <- mmkin(list("DFOP-SFO" = DFOP_SFO), ds_biphasic, error_model = "tc", quiet = TRUE)
+#'
+#' f_mixed <- mixed(f_mmkin)
+#' print(f_mixed)
+#' plot(f_mixed)
+#' }
+#' @export
+mixed <- function(object, ...) {
+  UseMethod("mixed")
+}
+
+#' @export
+#' @rdname mixed
+mixed.mmkin <- function(object, method = c("none"), ...) {
+  if (nrow(object) > 1) stop("Only row objects allowed")
+
+  method <- match.arg(method)
+  if (method == "default") method = c("naive", "nlme")
+
+  ds_names <- colnames(object)
+  res <- list(mmkin = object, mkinmod = object[[1]]$mkinmod)
+
+  if (method[1] == "none") {
+    ds_list <- lapply(object,
+      function(x) x$data[c("variable", "time", "observed", "predicted", "residual")])
+
+    names(ds_list) <- ds_names
+    res$data <- purrr::map_dfr(ds_list, function(x) x, .id = "ds")
+    names(res$data)[1:4] <- c("ds", "name", "time", "value")
+    res$data$name <- as.character(res$data$name)
+    res$data$ds <- ordered(res$data$ds, levels = unique(res$data$ds))
+    standardized <- unlist(lapply(object, residuals, standardized = TRUE))
+    res$data$std <- res$data$residual / standardized
+    res$data$standardized <- standardized
+
+    class(res) <- c("mixed.mmkin")
+    return(res)
+  }
+}
+
+#' @export
+#' @rdname mixed
+#' @param x A mixed.mmkin object to print
+#' @param digits Number of digits to use for printing.
+print.mixed.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) {
+  cat("Kinetic model fitted by nonlinear regression to each dataset" )
+  cat("\nStructural model:\n")
+  diffs <- x$mmkin[[1]]$mkinmod$diffs
+  nice_diffs <- gsub("^(d.*) =", "\\1/dt =", diffs)
+  writeLines(strwrap(nice_diffs, exdent = 11))
+  cat("\nData:\n")
+  cat(nrow(x$data), "observations of",
+    length(unique(x$data$name)), "variable(s) grouped in",
+    length(unique(x$data$ds)), "datasets\n\n")
+
+  print(x$mmkin)
+
+  cat("\nMean fitted parameters:\n")
+  print(mean_degparms(x$mmkin))
+
+  invisible(x)
+}
diff --git a/R/mkinfit.R b/R/mkinfit.R
index a6efc858..c659e446 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -885,8 +885,8 @@ mkinfit <- function(mkinmod, observed,
   fit$rss <- function(P) cost_function(P, OLS = TRUE, update_data = FALSE)
 
   # Log-likelihood with possibility to fix degparms or errparms
-  fit$ll <- function(P, fixed_degparms = FALSE, fixed_errparms = FALSE) {
-    - cost_function(P, trans = FALSE, fixed_degparms = fixed_degparms,
+  fit$ll <- function(P, fixed_degparms = FALSE, fixed_errparms = FALSE, trans = FALSE) {
+    - cost_function(P, trans = trans, fixed_degparms = fixed_degparms,
       fixed_errparms = fixed_errparms, OLS = FALSE, update_data = FALSE)
   }
 
diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R
index d8a1c2ac..db29376e 100644
--- a/R/plot.mixed.mmkin.R
+++ b/R/plot.mixed.mmkin.R
@@ -2,7 +2,7 @@ utils::globalVariables("ds")
 
 #' Plot predictions from a fitted nonlinear mixed model obtained via an mmkin row object
 #'
-#' @param x An object of class [saem.mmkin] or [nlme.mmkin]
+#' @param x An object of class [mixed.mmkin], [saem.mmkin] or [nlme.mmkin]
 #' @param i A numeric index to select datasets for which to plot the individual predictions,
 #'   in case plots get too large
 #' @inheritParams plot.mkinfit
@@ -63,31 +63,43 @@ plot.mixed.mmkin <- function(x,
   fit_1 <- x$mmkin[[1]]
   ds_names <- colnames(x$mmkin)
 
+  backtransform = TRUE
+
+  if (identical(class(x), "mixed.mmkin")) {
+    degparms_pop <- mean_degparms(x$mmkin)
+
+    degparms_tmp <- parms(x$mmkin, transformed = TRUE)
+    degparms_i <- as.data.frame(t(degparms_tmp[setdiff(rownames(degparms_tmp), names(fit_1$errparms)), ]))
+    residual_type = ifelse(standardized, "standardized", "residual")
+    residuals <- x$data[[residual_type]]
+  }
+
   if (inherits(x, "nlme.mmkin")) {
-    degparms_optim <- coefficients(x)
+    degparms_i <- coefficients(x)
     degparms_pop <- nlme::fixef(x)
     residuals <- residuals(x,
       type = ifelse(standardized, "pearson", "response"))
   }
 
   if (inherits(x, "saem.mmkin")) {
-    degparms_optim <- saemix::psi(x$so)
-    rownames(degparms_optim) <- ds_names
-    degparms_optim_names <- setdiff(names(fit_1$par), names(fit_1$errparms))
-    colnames(degparms_optim) <- degparms_optim_names
-    residual_type = ifelse(standardized, "residual", "standardized")
+    if (x$transformations == "saemix") backtransform = FALSE
+    degparms_i <- saemix::psi(x$so)
+    rownames(degparms_i) <- ds_names
+    degparms_i_names <- setdiff(x$so@results@name.fixed, names(fit_1$errparms))
+    colnames(degparms_i) <- degparms_i_names
+    residual_type = ifelse(standardized, "standardized", "residual")
     residuals <- x$data[[residual_type]]
     degparms_pop <- x$so@results@fixed.effects
-    names(degparms_pop) <- degparms_optim_names
+    names(degparms_pop) <- degparms_i_names
   }
 
   degparms_fixed <- fit_1$fixed$value
   names(degparms_fixed) <- rownames(fit_1$fixed)
-  degparms_all <- cbind(as.matrix(degparms_optim),
-    matrix(rep(degparms_fixed, nrow(degparms_optim)),
+  degparms_all <- cbind(as.matrix(degparms_i),
+    matrix(rep(degparms_fixed, nrow(degparms_i)),
       ncol = length(degparms_fixed),
-      nrow = nrow(degparms_optim), byrow = TRUE))
-  degparms_all_names <- c(names(degparms_optim), names(degparms_fixed))
+      nrow = nrow(degparms_i), byrow = TRUE))
+  degparms_all_names <- c(names(degparms_i), names(degparms_fixed))
   colnames(degparms_all) <- degparms_all_names
 
   degparms_all_pop <- c(degparms_pop, degparms_fixed)
@@ -106,10 +118,14 @@ plot.mixed.mmkin <- function(x,
   pred_ds <- purrr::map_dfr(i, function(ds_i)   {
     odeparms_trans <- degparms_all[ds_i, odeparms_names]
     names(odeparms_trans) <- odeparms_names # needed if only one odeparm
-    odeparms <- backtransform_odeparms(odeparms_trans,
-      x$mkinmod,
-      transform_rates = fit_1$transform_rates,
-      transform_fractions = fit_1$transform_fractions)
+    if (backtransform) {
+      odeparms <- backtransform_odeparms(odeparms_trans,
+        x$mkinmod,
+        transform_rates = fit_1$transform_rates,
+        transform_fractions = fit_1$transform_fractions)
+    } else {
+      odeparms <- odeparms_trans
+    }
 
     odeini <- degparms_all[ds_i, odeini_names]
     names(odeini) <- gsub("_0", "", odeini_names)
@@ -121,10 +137,15 @@ plot.mixed.mmkin <- function(x,
   })
 
   odeparms_pop_trans <- degparms_all_pop[odeparms_names]
-  odeparms_pop <- backtransform_odeparms(odeparms_pop_trans,
-    x$mkinmod,
-    transform_rates = fit_1$transform_rates,
-    transform_fractions = fit_1$transform_fractions)
+
+  if (backtransform) {
+    odeparms_pop <- backtransform_odeparms(odeparms_pop_trans,
+      x$mkinmod,
+      transform_rates = fit_1$transform_rates,
+      transform_fractions = fit_1$transform_fractions)
+  } else {
+    odeparms_pop <- odeparms_pop_trans
+  }
 
   odeini_pop <- degparms_all_pop[odeini_names]
   names(odeini_pop) <- gsub("_0", "", odeini_names)
diff --git a/R/saem.R b/R/saem.R
index e5634ac7..88b8e172 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -92,7 +92,7 @@ utils::globalVariables(c("predicted", "std"))
 #' plot(f_saem_fomc)
 #' }
 #' @export
-saem <- function(object, control, ...) UseMethod("saem")
+saem <- function(object, ...) UseMethod("saem")
 
 #' @rdname saem
 #' @export
@@ -104,9 +104,11 @@ saem.mmkin <- function(object,
   cores = 1,
   verbose = FALSE, suppressPlot = TRUE, quiet = FALSE, ...)
 {
+  transformations <- match.arg(transformations)
   m_saemix <- saemix_model(object, cores = cores, verbose = verbose,
     solution_type = solution_type, transformations = transformations, ...)
   d_saemix <- saemix_data(object, verbose = verbose)
+
   if (suppressPlot) {
     # We suppress the log-likelihood curve that saemix currently
     # produces at the end of the fit by plotting to a file
diff --git a/R/summary.saem.mmkin.R b/R/summary.saem.mmkin.R
index 97f9f2da..27c2ce6c 100644
--- a/R/summary.saem.mmkin.R
+++ b/R/summary.saem.mmkin.R
@@ -89,9 +89,42 @@ summary.saem.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes =
   confint_trans <- as.matrix(conf.int[pnames, c("estimate", "lower", "upper")])
   colnames(confint_trans)[1] <- "est."
 
-  bp <- backtransform_odeparms(confint_trans[, "est."], object$mkinmod,
-    object$transform_rates, object$transform_fractions)
-  bpnames <- names(bp)
+  if (object$transformations == "mkin") {
+    bp <- backtransform_odeparms(confint_trans[, "est."], object$mkinmod,
+      object$transform_rates, object$transform_fractions)
+    bpnames <- names(bp)
+
+    # Transform boundaries of CI for one parameter at a time,
+    # with the exception of sets of formation fractions (single fractions are OK).
+    f_names_skip <- character(0)
+    for (box in mod_vars) { # Figure out sets of fractions to skip
+      f_names <- grep(paste("^f", box, sep = "_"), pnames, value = TRUE)
+      n_paths <- length(f_names)
+      if (n_paths > 1) f_names_skip <- c(f_names_skip, f_names)
+    }
+
+    confint_back <- matrix(NA, nrow = length(bp), ncol = 3,
+      dimnames = list(bpnames, colnames(confint_trans)))
+    confint_back[, "est."] <- bp
+
+    for (pname in pnames) {
+      if (!pname %in% f_names_skip) {
+        par.lower <- confint_trans[pname, "lower"]
+        par.upper <- confint_trans[pname, "upper"]
+        names(par.lower) <- names(par.upper) <- pname
+        bpl <- backtransform_odeparms(par.lower, object$mkinmod,
+                                              object$transform_rates,
+                                              object$transform_fractions)
+        bpu <- backtransform_odeparms(par.upper, object$mkinmod,
+                                              object$transform_rates,
+                                              object$transform_fractions)
+        confint_back[names(bpl), "lower"] <- bpl
+        confint_back[names(bpu), "upper"] <- bpu
+      }
+    }
+  } else {
+    confint_back <- confint_trans
+  }
 
   #  Correlation of fixed effects (inspired by summary.nlme)
   varFix <- vcov(object$so)[1:np, 1:np]
@@ -111,34 +144,6 @@ summary.saem.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes =
   confint_errmod <- as.matrix(conf.int[enames, c("estimate", "lower", "upper")])
   colnames(confint_errmod)[1] <- "est."
 
-  # Transform boundaries of CI for one parameter at a time,
-  # with the exception of sets of formation fractions (single fractions are OK).
-  f_names_skip <- character(0)
-  for (box in mod_vars) { # Figure out sets of fractions to skip
-    f_names <- grep(paste("^f", box, sep = "_"), pnames, value = TRUE)
-    n_paths <- length(f_names)
-    if (n_paths > 1) f_names_skip <- c(f_names_skip, f_names)
-  }
-
-  confint_back <- matrix(NA, nrow = length(bp), ncol = 3,
-    dimnames = list(bpnames, colnames(confint_trans)))
-  confint_back[, "est."] <- bp
-
-  for (pname in pnames) {
-    if (!pname %in% f_names_skip) {
-      par.lower <- confint_trans[pname, "lower"]
-      par.upper <- confint_trans[pname, "upper"]
-      names(par.lower) <- names(par.upper) <- pname
-      bpl <- backtransform_odeparms(par.lower, object$mkinmod,
-                                            object$transform_rates,
-                                            object$transform_fractions)
-      bpu <- backtransform_odeparms(par.upper, object$mkinmod,
-                                            object$transform_rates,
-                                            object$transform_fractions)
-      confint_back[names(bpl), "lower"] <- bpl
-      confint_back[names(bpu), "upper"] <- bpu
-    }
-  }
 
   object$confint_trans <- confint_trans
   object$confint_ranef <- confint_ranef
@@ -213,7 +218,7 @@ print.summary.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3)
   print(data.frame(AIC = x$AIC, BIC = x$BIC, logLik = x$logLik,
       row.names = " "), digits = digits)
 
-  cat("\nOptimised, transformed parameters with symmetric confidence intervals:\n")
+  cat("\nOptimised parameters:\n")
   print(x$confint_trans, digits = digits)
 
   if (nrow(x$confint_trans) > 1) {
@@ -228,8 +233,10 @@ print.summary.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3)
   cat("\nVariance model:\n")
   print(x$confint_errmod, digits = digits)
 
-  cat("\nBacktransformed parameters with asymmetric confidence intervals:\n")
-  print(x$confint_back, digits = digits)
+  if (x$transformations == "mkin") {
+    cat("\nBacktransformed parameters:\n")
+    print(x$confint_back, digits = digits)
+  }
 
   printSFORB <- !is.null(x$SFORB)
   if(printSFORB){