aboutsummaryrefslogtreecommitdiff
path: root/R/saemix.R
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2020-11-07 12:54:58 +0100
committerJohannes Ranke <jranke@uni-bremen.de>2020-11-07 12:54:58 +0100
commit2e37905a7bf0d72751ee7326ef05a268989e49b2 (patch)
treecd65d87f8ff258304b26f3aa64e5a8c7216400d6 /R/saemix.R
parentcda47972e2b6a9610e3118dcd2270d7a1c76de3d (diff)
Create saem generic for fitting saemix models
The reasons for this decision were - Creating an saemix generic in the saemix package caused problems with roxygen, because functions like saemix.plot.xy were documented in their help files as S3 methods, although explicitly exported with @export - Creating an saemix generic in this package is possible, but would make it necessary to load samix with exclude = "saemix" in order to avoid overwriting the generic when loading saemix. - The return object of such an saemix generic in this package cannot be an S3 class with class attribute c("saemix.mmkin", "SaemixObject") similar to nlme.mmkin, as saemix returns an S4 class. - Extending the S4 class SaemixObject using simple inheritance to a class SaemixMmkinObject with additional slots did not work as expected. When the initialize method was left untouched, it prevented creation of an SaemixMmkinObject even if it was based on an initialised SaemixObject, as the initialize method seems to always be called by new(). This could potentially be circumvented by a coerce method. If an alternative initialize method was used, an SaemixMmkinObject could be created. However, the methods written for SaemixObjects only worked in some instances, either because they checked for the class, and not for class inheritance (like compare.saemix), or because the initialize method was called for some reason. Therefore, the idea of creating a derived S4 class was abandoned. - A side effect of this decision is that the introduction of the saem generic opens the possibility to use the same generic also for other backends like nlmixr with the SAEM algorithm.
Diffstat (limited to 'R/saemix.R')
-rw-r--r--R/saemix.R85
1 files changed, 45 insertions, 40 deletions
diff --git a/R/saemix.R b/R/saemix.R
index 280490a0..1db8b011 100644
--- a/R/saemix.R
+++ b/R/saemix.R
@@ -1,20 +1,16 @@
-#' Create saemix models
+#' Fit nonlinear mixed models with SAEM
#'
-#' The saemix function defined in this package is an S3 generic function
-#' using [saemix::saemix()] as its method for [saemix::SaemixModel] objects.
+#' This function uses [saemix::saemix()] as a backend for fitting nonlinear mixed
+#' effects models created from [mmkin] row objects using the stochastic approximation
+#' to the expectation maximisation algorithm (SAEM).
#'
-#' The method for mmkin row objects sets up a nonlinear mixed effects model for
-#' use with the saemix package. An mmkin row object is essentially a list of
-#' mkinfit objects that have been obtained by fitting the same model to
-#' a list of datasets.
+#' An mmkin row object is essentially a list of mkinfit objects that have been
+#' obtained by fitting the same model to a list of datasets using [mkinfit].
#'
-#' Starting values for the fixed effects (population mean parameters, argument psi0 of
-#' [saemix::saemixModel()] are the mean values of the parameters found using
-#' [mmkin].
+#' Starting values for the fixed effects (population mean parameters, argument
+#' psi0 of [saemix::saemixModel()] are the mean values of the parameters found
+#' using [mmkin].
#'
-#' @param model For the default method, this is an [saemix::saemixModel] object.
-#' If this is an [mmkin] row object, the [saemix::saemixModel] is created
-#' internally from the [mmkin] object.
#' @param object An [mmkin] row object containing several fits of the same
#' [mkinmod] model to different datasets
#' @param verbose Should we print information about created objects?
@@ -22,75 +18,84 @@
#' [parallel::mclapply()]. Using more than 1 core is experimental and may
#' lead to uncontrolled forking, apparently depending on the BLAS version
#' used.
+#' @param suppressPlot Should we suppress any plotting that is done
+#' by the saemix function?
+#' @param control Passed to [saemix::saemix]
#' @param \dots Further parameters passed to [saemix::saemixData]
#' and [saemix::saemixModel].
#' @return An [saemix::SaemixObject].
#' @examples
#' \dontrun{
-#' # We can load saemix, but should exclude the saemix function
-#' # as it would mask our generic version of it
-#' library(saemix, exclude = "saemix")
#' ds <- lapply(experimental_data_for_UBA_2019[6:10],
#' function(x) subset(x$data[c("name", "time", "value")]))
#' names(ds) <- paste("Dataset", 6:10)
#' f_mmkin_parent_p0_fixed <- mmkin("FOMC", ds, cores = 1,
#' state.ini = c(parent = 100), fixed_initials = "parent", quiet = TRUE)
-#' f_saemix_p0_fixed <- saemix(f_mmkin_parent_p0_fixed)
+#' f_saem_p0_fixed <- saem(f_mmkin_parent_p0_fixed)
#'
#' f_mmkin_parent <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE)
-#' f_saemix_sfo <- saemix(f_mmkin_parent["SFO", ])
-#' f_saemix_fomc <- saemix(f_mmkin_parent["FOMC", ])
-#' f_saemix_dfop <- saemix(f_mmkin_parent["DFOP", ])
+#' f_saem_sfo <- saem(f_mmkin_parent["SFO", ])
+#' f_saem_fomc <- saem(f_mmkin_parent["FOMC", ])
+#' f_saem_dfop <- saem(f_mmkin_parent["DFOP", ])
#'
-#' # As this returns an SaemixObject, we can use functions from saemix
-#' compare.saemix(list(f_saemix_sfo, f_saemix_fomc, f_saemix_dfop))
+#' # The returned saem.mmkin object contains an SaemixObject, we can use
+#' # functions from saemix
+#' library(saemix)
+#' compare.saemix(list(f_saem_sfo$so, f_saem_fomc$so, f_saem_dfop$so))
#'
#' f_mmkin_parent_tc <- update(f_mmkin_parent, error_model = "tc")
-#' f_saemix_fomc_tc <- saemix(f_mmkin_parent_tc["FOMC", ])
-#' compare.saemix(list(f_saemix_fomc, f_saemix_fomc_tc))
+#' f_saem_fomc_tc <- saem(f_mmkin_parent_tc["FOMC", ])
+#' compare.saemix(list(f_saem_fomc$so, f_saem_fomc_tc$so))
#'
#' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
#' A1 = mkinsub("SFO"))
#' f_mmkin <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, solution_type = "analytical")
#' # This takes about 4 minutes on my system
-#' f_saemix <- saemix(f_mmkin)
+#' f_saem <- saem(f_mmkin)
#'
-#' # Using a single core, it takes about 6 minutes, using 10 cores it is slower
-#' # instead of faster
#' f_mmkin_des <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, solution_type = "deSolve")
-#' f_saemix_des <- saemix(f_mmkin_des, cores = 1)
-#' compare.saemix(list(f_saemix, f_saemix_des))
-#'
+#' # Using a single core, the following takes about 6 minutes, using 10 cores
+#' # it is slower instead of faster
+#' f_saem_des <- saem(f_mmkin_des, cores = 1)
+#' compare.saemix(list(f_saemix$so, f_saemix_des$so))
#' }
#' @export
-saemix <- function(model, data, control, ...) UseMethod("saemix")
+saem <- function(object, control, ...) UseMethod("saem")
-#' @rdname saemix
+#' @rdname saem
#' @export
-saemix.mmkin <- function(model, data,
+saem.mmkin <- function(object,
control = list(displayProgress = FALSE, print = FALSE,
save = FALSE, save.graphs = FALSE),
cores = 1,
verbose = FALSE, suppressPlot = TRUE, ...)
{
- m_saemix <- saemix_model(model, cores = cores, verbose = verbose)
- d_saemix <- saemix_data(model, verbose = verbose)
+ m_saemix <- saemix_model(object, cores = cores, verbose = verbose)
+ 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
# that we discard afterwards
tmp <- tempfile()
- png(tmp)
+ grDevices::png(tmp)
}
- result <- saemix::saemix(m_saemix, d_saemix, control)
+ f_saemix <- saemix::saemix(m_saemix, d_saemix, control)
if (suppressPlot) {
- dev.off()
+ grDevices::dev.off()
unlink(tmp)
}
+ result <- list(
+ mkinmod = object[[1]]$mkinmod,
+ mmkin = object,
+ solution_type = object[[1]]$solution_type,
+ transform_rates = object[[1]]$transform_rates,
+ transform_fractions = object[[1]]$transform_fractions,
+ so = f_saemix)
+ class(result) <- "saem.mmkin"
return(result)
}
-#' @rdname saemix
+#' @rdname saem
#' @return An [saemix::SaemixModel] object.
#' @export
saemix_model <- function(object, cores = 1, verbose = FALSE, ...) {
@@ -254,7 +259,7 @@ saemix_model <- function(object, cores = 1, verbose = FALSE, ...) {
return(res)
}
-#' @rdname saemix
+#' @rdname saem
#' @return An [saemix::SaemixData] object.
#' @export
saemix_data <- function(object, verbose = FALSE, ...) {

Contact - Imprint