aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
Diffstat (limited to 'R')
-rw-r--r--R/dimethenamid_2018.R53
-rw-r--r--R/endpoints.R8
-rw-r--r--R/mean_degparms.R62
-rw-r--r--R/mixed.mmkin.R3
-rw-r--r--R/mkinsub.R5
-rw-r--r--R/nlme.R26
-rw-r--r--R/nlme.mmkin.R11
-rw-r--r--R/nlmixr.R539
-rw-r--r--R/plot.mixed.mmkin.R48
-rw-r--r--R/saem.R543
-rw-r--r--R/summary.nlmixr.mmkin.R250
-rw-r--r--R/summary.saem.mmkin.R268
-rw-r--r--R/tffm0.R48
-rw-r--r--R/transform_odeparms.R13
14 files changed, 1829 insertions, 48 deletions
diff --git a/R/dimethenamid_2018.R b/R/dimethenamid_2018.R
index 189da618..770649e2 100644
--- a/R/dimethenamid_2018.R
+++ b/R/dimethenamid_2018.R
@@ -15,7 +15,58 @@
#' @source Rapporteur Member State Germany, Co-Rapporteur Member State Bulgaria (2018)
#' Renewal Assessment Report Dimethenamid-P Volume 3 - B.8 Environmental fate and behaviour
#' Rev. 2 - November 2017
-#' \url{http://registerofquestions.efsa.europa.eu/roqFrontend/outputLoader?output=ON-5211}
+#' \url{https://open.efsa.europa.eu/study-inventory/EFSA-Q-2014-00716}
#' @examples
#' print(dimethenamid_2018)
+#' dmta_ds <- lapply(1:8, function(i) {
+#' ds_i <- dimethenamid_2018$ds[[i]]$data
+#' ds_i[ds_i$name == "DMTAP", "name"] <- "DMTA"
+#' ds_i$time <- ds_i$time * dimethenamid_2018$f_time_norm[i]
+#' ds_i
+#' })
+#' names(dmta_ds) <- sapply(dimethenamid_2018$ds, function(ds) ds$title)
+#' dmta_ds[["Borstel"]] <- rbind(dmta_ds[["Borstel 1"]], dmta_ds[["Borstel 2"]])
+#' dmta_ds[["Borstel 1"]] <- NULL
+#' dmta_ds[["Borstel 2"]] <- NULL
+#' dmta_ds[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]])
+#' dmta_ds[["Elliot 1"]] <- NULL
+#' dmta_ds[["Elliot 2"]] <- NULL
+#' \dontrun{
+#' dfop_sfo3_plus <- mkinmod(
+#' DMTA = mkinsub("DFOP", c("M23", "M27", "M31")),
+#' M23 = mkinsub("SFO"),
+#' M27 = mkinsub("SFO"),
+#' M31 = mkinsub("SFO", "M27", sink = FALSE),
+#' quiet = TRUE
+#' )
+#' f_dmta_mkin_tc <- mmkin(
+#' list("DFOP-SFO3+" = dfop_sfo3_plus),
+#' dmta_ds, quiet = TRUE, error_model = "tc")
+#' nlmixr_model(f_dmta_mkin_tc)
+#' # The focei fit takes about four minutes on my system
+#' system.time(
+#' f_dmta_nlmixr_focei <- nlmixr(f_dmta_mkin_tc, est = "focei",
+#' control = nlmixr::foceiControl(print = 500))
+#' )
+#' summary(f_dmta_nlmixr_focei)
+#' plot(f_dmta_nlmixr_focei)
+#' # Using saemix takes about 18 minutes
+#' system.time(
+#' f_dmta_saemix <- saem(f_dmta_mkin_tc, test_log_parms = TRUE)
+#' )
+#'
+#' # nlmixr with est = "saem" is pretty fast with default iteration numbers, most
+#' # of the time (about 2.5 minutes) is spent for calculating the log likelihood at the end
+#' # The likelihood calculated for the nlmixr fit is much lower than that found by saemix
+#' # Also, the trace plot and the plot of the individual predictions is not
+#' # convincing for the parent. It seems we are fitting an overparameterised
+#' # model, so the result we get strongly depends on starting parameters and control settings.
+#' system.time(
+#' f_dmta_nlmixr_saem <- nlmixr(f_dmta_mkin_tc, est = "saem",
+#' control = nlmixr::saemControl(print = 500, logLik = TRUE, nmc = 9))
+#' )
+#' traceplot(f_dmta_nlmixr_saem$nm)
+#' summary(f_dmta_nlmixr_saem)
+#' plot(f_dmta_nlmixr_saem)
+#' }
"dimethenamid_2018"
diff --git a/R/endpoints.R b/R/endpoints.R
index b5872e68..6bf52f07 100644
--- a/R/endpoints.R
+++ b/R/endpoints.R
@@ -10,8 +10,8 @@
#' Additional DT50 values are calculated from the FOMC DT90 and k1 and k2 from
#' HS and DFOP, as well as from Eigenvalues b1 and b2 of any SFORB models
#'
-#' @param fit An object of class [mkinfit] or [nlme.mmkin]
-#' or another object that has list components
+#' @param fit An object of class [mkinfit], [nlme.mmkin], [saem.mmkin] or
+#' [nlmixr.mmkin]. Or another object that has list components
#' mkinmod containing an [mkinmod] degradation model, and two numeric vectors,
#' bparms.optim and bparms.fixed, that contain parameter values
#' for that model.
@@ -20,8 +20,8 @@
#' and, if applicable, a vector of formation fractions named ff
#' and, if the SFORB model was in use, a vector of eigenvalues
#' of these SFORB models, equivalent to DFOP rate constants
-#' @note The function is used internally by [summary.mkinfit]
-#' and [summary.nlme.mmkin]
+#' @note The function is used internally by [summary.mkinfit],
+#' [summary.nlme.mmkin] and [summary.saem.mmkin].
#' @author Johannes Ranke
#' @examples
#'
diff --git a/R/mean_degparms.R b/R/mean_degparms.R
new file mode 100644
index 00000000..ec20c068
--- /dev/null
+++ b/R/mean_degparms.R
@@ -0,0 +1,62 @@
+#' Calculate mean degradation parameters for an mmkin row object
+#'
+#' @return If random is FALSE (default), a named vector containing mean values
+#' of the fitted degradation model parameters. If random is TRUE, a list with
+#' fixed and random effects, in the format required by the start argument of
+#' nlme for the case of a single grouping variable ds.
+#' @param object An mmkin row object containing several fits of the same model to different datasets
+#' @param random Should a list with fixed and random effects be returned?
+#' @param test_log_parms If TRUE, log parameters are only considered in
+#' the mean calculations if their untransformed counterparts (most likely
+#' rate constants) pass the t-test for significant difference from zero.
+#' @param conf.level Possibility to adjust the required confidence level
+#' for parameter that are tested if requested by 'test_log_parms'.
+#' @export
+mean_degparms <- function(object, random = FALSE, test_log_parms = FALSE, conf.level = 0.6)
+{
+ if (nrow(object) > 1) stop("Only row objects allowed")
+ parm_mat_trans <- sapply(object, parms, transformed = TRUE)
+
+ if (test_log_parms) {
+ parm_mat_dim <- dim(parm_mat_trans)
+ parm_mat_dimnames <- dimnames(parm_mat_trans)
+
+ log_parm_trans_names <- grep("^log_", rownames(parm_mat_trans), value = TRUE)
+ log_parm_names <- gsub("^log_", "", log_parm_trans_names)
+
+ t_test_back_OK <- matrix(
+ sapply(object, function(o) {
+ suppressWarnings(summary(o)$bpar[log_parm_names, "Pr(>t)"] < (1 - conf.level))
+ }), nrow = length(log_parm_names))
+ rownames(t_test_back_OK) <- log_parm_trans_names
+
+ parm_mat_trans_OK <- parm_mat_trans
+ for (trans_parm in log_parm_trans_names) {
+ parm_mat_trans_OK[trans_parm, ] <- ifelse(t_test_back_OK[trans_parm, ],
+ parm_mat_trans[trans_parm, ], NA)
+ }
+ } else {
+ parm_mat_trans_OK <- parm_mat_trans
+ }
+
+ mean_degparm_names <- setdiff(rownames(parm_mat_trans), names(object[[1]]$errparms))
+ degparm_mat_trans <- parm_mat_trans[mean_degparm_names, , drop = FALSE]
+ degparm_mat_trans_OK <- parm_mat_trans_OK[mean_degparm_names, , drop = FALSE]
+
+ fixed <- apply(degparm_mat_trans_OK, 1, mean, na.rm = TRUE)
+ if (random) {
+ random <- t(apply(degparm_mat_trans[mean_degparm_names, , drop = FALSE], 2, function(column) column - fixed))
+ # If we only have one parameter, apply returns a vector so we get a single row
+ if (nrow(degparm_mat_trans) == 1) random <- t(random)
+ rownames(random) <- levels(nlme_data(object)$ds)
+
+ # For nlmixr we can specify starting values for standard deviations eta, and
+ # we ignore uncertain parameters if test_log_parms is FALSE
+ eta <- apply(degparm_mat_trans_OK, 1, stats::sd, na.rm = TRUE)
+
+ return(list(fixed = fixed, random = list(ds = random), eta = eta))
+ } else {
+ return(fixed)
+ }
+}
+
diff --git a/R/mixed.mmkin.R b/R/mixed.mmkin.R
index 7aa5edd5..682a9a34 100644
--- a/R/mixed.mmkin.R
+++ b/R/mixed.mmkin.R
@@ -3,6 +3,8 @@
#' @param object An [mmkin] row object
#' @param method The method to be used
#' @param \dots Currently not used
+#' @return An object of class 'mixed.mmkin' which has the observed data in a
+#' single dataframe which is convenient for plotting
#' @examples
#' sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
#' n_biphasic <- 8
@@ -54,7 +56,6 @@ 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)
diff --git a/R/mkinsub.R b/R/mkinsub.R
index 886f712c..93af3f16 100644
--- a/R/mkinsub.R
+++ b/R/mkinsub.R
@@ -1,8 +1,3 @@
-#' Function to set up a kinetic submodel for one state variable
-#'
-#' This is a convenience function to set up the lists used as arguments for
-#' \code{\link{mkinmod}}.
-#'
#' @rdname mkinmod
#' @param submodel Character vector of length one to specify the submodel type.
#' See \code{\link{mkinmod}} for the list of allowed submodel names.
diff --git a/R/nlme.R b/R/nlme.R
index 9215aab0..8762f137 100644
--- a/R/nlme.R
+++ b/R/nlme.R
@@ -36,7 +36,7 @@
#' nlme_f <- nlme_function(f)
#' # These assignments are necessary for these objects to be
#' # visible to nlme and augPred when evaluation is done by
-#' # pkgdown to generated the html docs.
+#' # pkgdown to generate the html docs.
#' assign("nlme_f", nlme_f, globalenv())
#' assign("grouped_data", grouped_data, globalenv())
#'
@@ -125,30 +125,6 @@ nlme_function <- function(object) {
}
#' @rdname nlme
-#' @return If random is FALSE (default), a named vector containing mean values
-#' of the fitted degradation model parameters. If random is TRUE, a list with
-#' fixed and random effects, in the format required by the start argument of
-#' nlme for the case of a single grouping variable ds.
-#' @param random Should a list with fixed and random effects be returned?
-#' @export
-mean_degparms <- function(object, random = FALSE) {
- if (nrow(object) > 1) stop("Only row objects allowed")
- parm_mat_trans <- sapply(object, parms, transformed = TRUE)
- mean_degparm_names <- setdiff(rownames(parm_mat_trans), names(object[[1]]$errparms))
- degparm_mat_trans <- parm_mat_trans[mean_degparm_names, , drop = FALSE]
- fixed <- apply(degparm_mat_trans, 1, mean)
- if (random) {
- random <- t(apply(degparm_mat_trans[mean_degparm_names, , drop = FALSE], 2, function(column) column - fixed))
- # If we only have one parameter, apply returns a vector so we get a single row
- if (nrow(degparm_mat_trans) == 1) random <- t(random)
- rownames(random) <- levels(nlme_data(object)$ds)
- return(list(fixed = fixed, random = list(ds = random)))
- } else {
- return(fixed)
- }
-}
-
-#' @rdname nlme
#' @importFrom purrr map_dfr
#' @return A \code{\link{groupedData}} object
#' @export
diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R
index ff1f2fff..7049a9a1 100644
--- a/R/nlme.mmkin.R
+++ b/R/nlme.mmkin.R
@@ -24,7 +24,7 @@ get_deg_func <- function() {
#' This functions sets up a nonlinear mixed effects model for an mmkin row
#' object. 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.
-#'
+#'
#' Note that the convergence of the nlme algorithms depends on the quality
#' of the data. In degradation kinetics, we often only have few datasets
#' (e.g. data for few soils) and complicated degradation models, which may
@@ -34,10 +34,9 @@ get_deg_func <- function() {
#' @param data Ignored, data are taken from the mmkin model
#' @param fixed Ignored, all degradation parameters fitted in the
#' mmkin model are used as fixed parameters
-#' @param random If not specified, correlated random effects are set up
-#' for all optimised degradation model parameters using the log-Cholesky
-#' parameterization [nlme::pdLogChol] that is also the default of
-#' the generic [nlme] method.
+#' @param random If not specified, no correlations between random effects are
+#' set up for the optimised degradation model parameters. This is
+#' achieved by using the [nlme::pdDiag] method.
#' @param groups See the documentation of nlme
#' @param start If not specified, mean values of the fitted degradation
#' parameters taken from the mmkin object are used
@@ -135,7 +134,7 @@ nlme.mmkin <- function(model, data = "auto",
function(el) eval(parse(text = paste(el, 1, sep = "~")))),
random = pdDiag(fixed),
groups,
- start = mean_degparms(model, random = TRUE),
+ start = mean_degparms(model, random = TRUE, test_log_parms = TRUE),
correlation = NULL, weights = NULL,
subset, method = c("ML", "REML"),
na.action = na.fail, naPattern,
diff --git a/R/nlmixr.R b/R/nlmixr.R
new file mode 100644
index 00000000..9c364c4f
--- /dev/null
+++ b/R/nlmixr.R
@@ -0,0 +1,539 @@
+utils::globalVariables(c("predicted", "std", "ID", "TIME", "CMT", "DV", "IPRED", "IRES", "IWRES"))
+
+#' @export
+nlmixr::nlmixr
+
+#' Fit nonlinear mixed models using nlmixr
+#'
+#' This function uses [nlmixr::nlmixr()] as a backend for fitting nonlinear mixed
+#' effects models created from [mmkin] row objects using the Stochastic Approximation
+#' Expectation Maximisation algorithm (SAEM).
+#'
+#' 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].
+#'
+#' @importFrom nlmixr nlmixr tableControl
+#' @importFrom dplyr %>%
+#' @param object An [mmkin] row object containing several fits of the same
+#' [mkinmod] model to different datasets
+#' @param data Not used, the data are extracted from the mmkin row object
+#' @param est Estimation method passed to [nlmixr::nlmixr]
+#' @param degparms_start Parameter values given as a named numeric vector will
+#' be used to override the starting values obtained from the 'mmkin' object.
+#' @param eta_start Standard deviations on the transformed scale given as a
+#' named numeric vector will be used to override the starting values obtained
+#' from the 'mmkin' object.
+#' @param test_log_parms If TRUE, an attempt is made to use more robust starting
+#' values for population parameters fitted as log parameters in mkin (like
+#' rate constants) by only considering rate constants that pass the t-test
+#' when calculating mean degradation parameters using [mean_degparms].
+#' @param conf.level Possibility to adjust the required confidence level
+#' for parameter that are tested if requested by 'test_log_parms'.
+#' @param data Not used, as the data are extracted from the mmkin row object
+#' @param table Passed to [nlmixr::nlmixr]
+#' @param error_model Possibility to override the error model which is being
+#' set based on the error model used in the mmkin row object.
+#' @param control Passed to [nlmixr::nlmixr]
+#' @param \dots Passed to [nlmixr_model]
+#' @param save Passed to [nlmixr::nlmixr]
+#' @param envir Passed to [nlmixr::nlmixr]
+#' @return An S3 object of class 'nlmixr.mmkin', containing the fitted
+#' [nlmixr::nlmixr] object as a list component named 'nm'. The
+#' object also inherits from 'mixed.mmkin'.
+#' @seealso [summary.nlmixr.mmkin] [plot.mixed.mmkin]
+#' @examples
+#' \dontrun{
+#' 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 <- mmkin(c("SFO", "FOMC", "DFOP", "HS"), ds, quiet = TRUE, cores = 1)
+#' f_mmkin_parent_tc <- mmkin(c("SFO", "FOMC", "DFOP"), ds, error_model = "tc",
+#' cores = 1, quiet = TRUE)
+#'
+#' f_nlmixr_sfo_saem <- nlmixr(f_mmkin_parent["SFO", ], est = "saem")
+#' f_nlmixr_sfo_focei <- nlmixr(f_mmkin_parent["SFO", ], est = "focei")
+#'
+#' f_nlmixr_fomc_saem <- nlmixr(f_mmkin_parent["FOMC", ], est = "saem")
+#' f_nlmixr_fomc_focei <- nlmixr(f_mmkin_parent["FOMC", ], est = "focei")
+#'
+#' f_nlmixr_dfop_saem <- nlmixr(f_mmkin_parent["DFOP", ], est = "saem")
+#' f_nlmixr_dfop_focei <- nlmixr(f_mmkin_parent["DFOP", ], est = "focei")
+#'
+#' f_nlmixr_hs_saem <- nlmixr(f_mmkin_parent["HS", ], est = "saem")
+#' f_nlmixr_hs_focei <- nlmixr(f_mmkin_parent["HS", ], est = "focei")
+#'
+#' f_nlmixr_fomc_saem_tc <- nlmixr(f_mmkin_parent_tc["FOMC", ], est = "saem")
+#' f_nlmixr_fomc_focei_tc <- nlmixr(f_mmkin_parent_tc["FOMC", ], est = "focei")
+#'
+#' AIC(
+#' f_nlmixr_sfo_saem$nm, f_nlmixr_sfo_focei$nm,
+#' f_nlmixr_fomc_saem$nm, f_nlmixr_fomc_focei$nm,
+#' f_nlmixr_dfop_saem$nm, f_nlmixr_dfop_focei$nm,
+#' f_nlmixr_hs_saem$nm, f_nlmixr_hs_focei$nm,
+#' f_nlmixr_fomc_saem_tc$nm, f_nlmixr_fomc_focei_tc$nm)
+#'
+#' AIC(nlme(f_mmkin_parent["FOMC", ]))
+#' AIC(nlme(f_mmkin_parent["HS", ]))
+#'
+#' # nlme is comparable to nlmixr with focei, saem finds a better
+#' # solution, the two-component error model does not improve it
+#' plot(f_nlmixr_fomc_saem)
+#'
+#' sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"),
+#' A1 = mkinsub("SFO"))
+#' fomc_sfo <- mkinmod(parent = mkinsub("FOMC", "A1"),
+#' A1 = mkinsub("SFO"))
+#' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
+#' A1 = mkinsub("SFO"))
+#'
+#' f_mmkin_const <- mmkin(list(
+#' "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo),
+#' ds, quiet = TRUE, error_model = "const")
+#' f_mmkin_obs <- mmkin(list(
+#' "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo),
+#' ds, quiet = TRUE, error_model = "obs")
+#' f_mmkin_tc <- mmkin(list(
+#' "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo),
+#' ds, quiet = TRUE, error_model = "tc")
+#'
+#' # A single constant variance is currently only possible with est = 'focei' in nlmixr
+#' f_nlmixr_sfo_sfo_focei_const <- nlmixr(f_mmkin_const["SFO-SFO", ], est = "focei")
+#' f_nlmixr_fomc_sfo_focei_const <- nlmixr(f_mmkin_const["FOMC-SFO", ], est = "focei")
+#' f_nlmixr_dfop_sfo_focei_const <- nlmixr(f_mmkin_const["DFOP-SFO", ], est = "focei")
+#'
+#' # Variance by variable is supported by 'saem' and 'focei'
+#' f_nlmixr_fomc_sfo_saem_obs <- nlmixr(f_mmkin_obs["FOMC-SFO", ], est = "saem")
+#' f_nlmixr_fomc_sfo_focei_obs <- nlmixr(f_mmkin_obs["FOMC-SFO", ], est = "focei")
+#' f_nlmixr_dfop_sfo_saem_obs <- nlmixr(f_mmkin_obs["DFOP-SFO", ], est = "saem")
+#' f_nlmixr_dfop_sfo_focei_obs <- nlmixr(f_mmkin_obs["DFOP-SFO", ], est = "focei")
+#'
+#' # Identical two-component error for all variables is only possible with
+#' # est = 'focei' in nlmixr
+#' f_nlmixr_fomc_sfo_focei_tc <- nlmixr(f_mmkin_tc["FOMC-SFO", ], est = "focei")
+#' f_nlmixr_dfop_sfo_focei_tc <- nlmixr(f_mmkin_tc["DFOP-SFO", ], est = "focei")
+#'
+#' # Two-component error by variable is possible with both estimation methods
+#' # Variance by variable is supported by 'saem' and 'focei'
+#' f_nlmixr_fomc_sfo_saem_obs_tc <- nlmixr(f_mmkin_tc["FOMC-SFO", ], est = "saem",
+#' error_model = "obs_tc")
+#' f_nlmixr_fomc_sfo_focei_obs_tc <- nlmixr(f_mmkin_tc["FOMC-SFO", ], est = "focei",
+#' error_model = "obs_tc")
+#' f_nlmixr_dfop_sfo_saem_obs_tc <- nlmixr(f_mmkin_tc["DFOP-SFO", ], est = "saem",
+#' error_model = "obs_tc")
+#' f_nlmixr_dfop_sfo_focei_obs_tc <- nlmixr(f_mmkin_tc["DFOP-SFO", ], est = "focei",
+#' error_model = "obs_tc")
+#'
+#' AIC(
+#' f_nlmixr_sfo_sfo_focei_const$nm,
+#' f_nlmixr_fomc_sfo_focei_const$nm,
+#' f_nlmixr_dfop_sfo_focei_const$nm,
+#' f_nlmixr_fomc_sfo_saem_obs$nm,
+#' f_nlmixr_fomc_sfo_focei_obs$nm,
+#' f_nlmixr_dfop_sfo_saem_obs$nm,
+#' f_nlmixr_dfop_sfo_focei_obs$nm,
+#' f_nlmixr_fomc_sfo_focei_tc$nm,
+#' f_nlmixr_dfop_sfo_focei_tc$nm,
+#' f_nlmixr_fomc_sfo_saem_obs_tc$nm,
+#' f_nlmixr_fomc_sfo_focei_obs_tc$nm,
+#' f_nlmixr_dfop_sfo_saem_obs_tc$nm,
+#' f_nlmixr_dfop_sfo_focei_obs_tc$nm
+#' )
+#' # Currently, FOMC-SFO with two-component error by variable fitted by focei gives the
+#' # lowest AIC
+#' plot(f_nlmixr_fomc_sfo_focei_obs_tc)
+#' summary(f_nlmixr_fomc_sfo_focei_obs_tc)
+#' }
+#' @export
+nlmixr.mmkin <- function(object, data = NULL,
+ est = NULL, control = list(),
+ table = tableControl(),
+ error_model = object[[1]]$err_mod,
+ test_log_parms = TRUE,
+ conf.level = 0.6,
+ degparms_start = "auto",
+ eta_start = "auto",
+ ...,
+ save = NULL,
+ envir = parent.frame()
+)
+{
+ m_nlmixr <- nlmixr_model(object, est = est,
+ error_model = error_model, add_attributes = TRUE,
+ test_log_parms = test_log_parms, conf.level = conf.level,
+ degparms_start = degparms_start, eta_start = eta_start
+ )
+ d_nlmixr <- nlmixr_data(object)
+
+ mean_dp_start <- attr(m_nlmixr, "mean_dp_start")
+ mean_ep_start <- attr(m_nlmixr, "mean_ep_start")
+
+ attributes(m_nlmixr) <- NULL
+
+ fit_time <- system.time({
+ f_nlmixr <- nlmixr(m_nlmixr, d_nlmixr, est = est, control = control)
+ })
+
+ if (is.null(f_nlmixr$CMT)) {
+ nlmixr_df <- as.data.frame(f_nlmixr[c("ID", "TIME", "DV", "IPRED", "IRES", "IWRES")])
+ nlmixr_df$CMT <- as.character(object[[1]]$data$variable[1])
+ } else {
+ nlmixr_df <- as.data.frame(f_nlmixr[c("ID", "TIME", "DV", "CMT", "IPRED", "IRES", "IWRES")])
+ }
+
+ return_data <- nlmixr_df %>%
+ dplyr::transmute(ds = ID, name = CMT, time = TIME, value = DV,
+ predicted = IPRED, residual = IRES,
+ std = IRES/IWRES, standardized = IWRES) %>%
+ dplyr::arrange(ds, name, time)
+
+ bparms_optim <- backtransform_odeparms(f_nlmixr$theta,
+ object[[1]]$mkinmod,
+ object[[1]]$transform_rates,
+ object[[1]]$transform_fractions)
+
+ result <- list(
+ mkinmod = object[[1]]$mkinmod,
+ mmkin = object,
+ transform_rates = object[[1]]$transform_rates,
+ transform_fractions = object[[1]]$transform_fractions,
+ nm = f_nlmixr,
+ est = est,
+ time = fit_time,
+ mean_dp_start = mean_dp_start,
+ mean_ep_start = mean_ep_start,
+ bparms.optim = bparms_optim,
+ bparms.fixed = object[[1]]$bparms.fixed,
+ data = return_data,
+ err_mod = error_model,
+ date.fit = date(),
+ nlmixrversion = as.character(utils::packageVersion("nlmixr")),
+ mkinversion = as.character(utils::packageVersion("mkin")),
+ Rversion = paste(R.version$major, R.version$minor, sep=".")
+ )
+
+ class(result) <- c("nlmixr.mmkin", "mixed.mmkin")
+ return(result)
+}
+
+#' @export
+#' @rdname nlmixr.mmkin
+#' @param x An nlmixr.mmkin object to print
+#' @param digits Number of digits to use for printing
+print.nlmixr.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) {
+ cat("Kinetic nonlinear mixed-effects model fit by", x$est, "using nlmixr")
+ 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")
+
+ cat("\nLikelihood:\n")
+ print(data.frame(
+ AIC = AIC(x$nm),
+ BIC = BIC(x$nm),
+ logLik = logLik(x$nm),
+ row.names = " "), digits = digits)
+
+ cat("\nFitted parameters:\n")
+ print(x$nm$parFixed, digits = digits)
+
+ invisible(x)
+}
+
+#' @rdname nlmixr.mmkin
+#' @param add_attributes Should the starting values used for degradation model
+#' parameters and their distribution and for the error model parameters
+#' be returned as attributes?
+#' @return An function defining a model suitable for fitting with [nlmixr::nlmixr].
+#' @export
+nlmixr_model <- function(object,
+ est = c("saem", "focei"),
+ degparms_start = "auto",
+ eta_start = "auto",
+ test_log_parms = TRUE, conf.level = 0.6,
+ error_model = object[[1]]$err_mod, add_attributes = FALSE)
+{
+ if (nrow(object) > 1) stop("Only row objects allowed")
+ est = match.arg(est)
+
+ mkin_model <- object[[1]]$mkinmod
+ obs_vars <- names(mkin_model$spec)
+
+ if (error_model == object[[1]]$err_mod) {
+
+ if (length(object[[1]]$mkinmod$spec) > 1 & est == "saem") {
+ if (error_model == "const") {
+ message(
+ "Constant variance for more than one variable is not supported for est = 'saem'\n",
+ "Changing the error model to 'obs' (variance by observed variable)")
+ error_model <- "obs"
+ }
+ if (error_model =="tc") {
+ message(
+ "With est = 'saem', a different error model is required for each observed variable",
+ "Changing the error model to 'obs_tc' (Two-component error for each observed variable)")
+ error_model <- "obs_tc"
+ }
+ }
+ }
+
+ degparms_mmkin <- mean_degparms(object,
+ test_log_parms = test_log_parms,
+ conf.level = conf.level, random = TRUE)
+
+ degparms_optim <- degparms_mmkin$fixed
+
+ degparms_optim_ilr_names <- grep("^f_.*_ilr", names(degparms_optim), value = TRUE)
+ obs_vars_ilr <- unique(gsub("f_(.*)_ilr.*$", "\\1", degparms_optim_ilr_names))
+ degparms_optim_noilr <- degparms_optim[setdiff(names(degparms_optim),
+ degparms_optim_ilr_names)]
+
+ degparms_optim_back <- backtransform_odeparms(degparms_optim,
+ object[[1]]$mkinmod,
+ object[[1]]$transform_rates,
+ object[[1]]$transform_fractions)
+
+ if (degparms_start[1] == "auto") {
+ degparms_start <- degparms_optim_noilr
+ for (obs_var_ilr in obs_vars_ilr) {
+ ff_names <- grep(paste0("^f_", obs_var_ilr, "_"),
+ names(degparms_optim_back), value = TRUE)
+ f_tffm0 <- tffm0(degparms_optim_back[ff_names])
+ f_tffm0_qlogis <- qlogis(f_tffm0)
+ names(f_tffm0_qlogis) <- paste0("f_", obs_var_ilr,
+ "_tffm0_", 1:length(f_tffm0), "_qlogis")
+ degparms_start <- c(degparms_start, f_tffm0_qlogis)
+ }
+ }
+
+ if (eta_start[1] == "auto") {
+ eta_start <- degparms_mmkin$eta[setdiff(names(degparms_optim),
+ degparms_optim_ilr_names)]
+ for (obs_var_ilr in obs_vars_ilr) {
+ ff_n <- length(grep(paste0("^f_", obs_var_ilr, "_"),
+ names(degparms_optim_back), value = TRUE))
+ eta_start_ff <- rep(0.3, ff_n)
+ names(eta_start_ff) <- paste0("f_", obs_var_ilr,
+ "_tffm0_", 1:ff_n, "_qlogis")
+ eta_start <- c(eta_start, eta_start_ff)
+ }
+ }
+
+
+ degparms_fixed <- object[[1]]$bparms.fixed
+
+ odeini_optim_parm_names <- grep('_0$', names(degparms_optim), value = TRUE)
+ odeini_fixed_parm_names <- grep('_0$', names(degparms_fixed), value = TRUE)
+
+ odeparms_fixed_names <- setdiff(names(degparms_fixed), odeini_fixed_parm_names)
+ odeparms_fixed <- degparms_fixed[odeparms_fixed_names]
+
+ odeini_fixed <- degparms_fixed[odeini_fixed_parm_names]
+ names(odeini_fixed) <- gsub('_0$', '', odeini_fixed_parm_names)
+
+ # Definition of the model function
+ f <- function(){}
+
+ ini_block <- "ini({"
+
+ # Initial values for all degradation parameters
+ for (parm_name in names(degparms_start)) {
+ # As initials for state variables are not transformed,
+ # we need to modify the name here as we want to
+ # use the original name in the model block
+ ini_block <- paste0(
+ ini_block,
+ parm_name, " = ",
+ as.character(degparms_start[parm_name]),
+ "\n",
+ "eta.", parm_name, " ~ ",
+ as.character(eta_start[parm_name]),
+ "\n"
+ )
+ }
+
+ # Error model parameters
+ error_model_mkin <- object[[1]]$err_mod
+
+ errparm_names_mkin <- names(object[[1]]$errparms)
+ errparms_mkin <- sapply(errparm_names_mkin, function(parm_name) {
+ mean(sapply(object, function(x) x$errparms[parm_name]))
+ })
+
+ sigma_tc_mkin <- errparms_ini <- errparms_mkin[1] +
+ mean(unlist(sapply(object, function(x) x$data$observed)), na.rm = TRUE) *
+ errparms_mkin[2]
+
+ if (error_model == "const") {
+ if (error_model_mkin == "tc") {
+ errparms_ini <- sigma_tc_mkin
+ } else {
+ errparms_ini <- mean(errparms_mkin)
+ }
+ names(errparms_ini) <- "sigma"
+ }
+
+ if (error_model == "obs") {
+ errparms_ini <- switch(error_model_mkin,
+ const = rep(errparms_mkin["sigma"], length(obs_vars)),
+ obs = errparms_mkin,
+ tc = sigma_tc_mkin)
+ names(errparms_ini) <- paste0("sigma_", obs_vars)
+ }
+
+ if (error_model == "tc") {
+ if (error_model_mkin != "tc") {
+ stop("Not supported")
+ } else {
+ errparms_ini <- errparms_mkin
+ }
+ }
+
+ if (error_model == "obs_tc") {
+ if (error_model_mkin != "tc") {
+ stop("Not supported")
+ } else {
+ errparms_ini <- rep(errparms_mkin, length(obs_vars))
+ names(errparms_ini) <- paste0(
+ rep(names(errparms_mkin), length(obs_vars)),
+ "_",
+ rep(obs_vars, each = 2))
+ }
+ }
+
+ for (parm_name in names(errparms_ini)) {
+ ini_block <- paste0(
+ ini_block,
+ parm_name, " = ",
+ as.character(errparms_ini[parm_name]),
+ "\n"
+ )
+ }
+
+ ini_block <- paste0(ini_block, "})")
+
+ body(f)[2] <- parse(text = ini_block)
+
+ model_block <- "model({"
+
+ # Population initial values for the ODE state variables
+ for (parm_name in odeini_optim_parm_names) {
+ model_block <- paste0(
+ model_block,
+ parm_name, "_model = ",
+ parm_name, " + eta.", parm_name, "\n",
+ gsub("(.*)_0", "\\1(0)", parm_name), " = ", parm_name, "_model\n")
+ }
+
+ # Population initial values for log rate constants
+ for (parm_name in grep("^log_", names(degparms_start), value = TRUE)) {
+ model_block <- paste0(
+ model_block,
+ gsub("^log_", "", parm_name), " = ",
+ "exp(", parm_name, " + eta.", parm_name, ")\n")
+ }
+
+ # Population initial values for logit transformed parameters
+ for (parm_name in grep("_qlogis$", names(degparms_start), value = TRUE)) {
+ model_block <- paste0(
+ model_block,
+ gsub("_qlogis$", "", parm_name), " = ",
+ "expit(", parm_name, " + eta.", parm_name, ")\n")
+ }
+
+ # Calculate formation fractions from tffm0 transformed values
+ for (obs_var_ilr in obs_vars_ilr) {
+ ff_names <- grep(paste0("^f_", obs_var_ilr, "_"),
+ names(degparms_optim_back), value = TRUE)
+ pattern <- paste0("^f_", obs_var_ilr, "_to_(.*)$")
+ target_vars <- gsub(pattern, "\\1",
+ grep(paste0("^f_", obs_var_ilr, "_to_"), names(degparms_optim_back), value = TRUE))
+ for (i in 1:length(target_vars)) {
+ ff_name <- ff_names[i]
+ ff_line <- paste0(ff_name, " = f_", obs_var_ilr, "_tffm0_", i)
+ if (i > 1) {
+ for (j in (i - 1):1) {
+ ff_line <- paste0(ff_line, " * (1 - f_", obs_var_ilr, "_tffm0_", j , ")")
+ }
+ }
+ model_block <- paste0(
+ model_block,
+ ff_line,
+ "\n"
+ )
+ }
+ }
+
+ # Differential equations
+ model_block <- paste0(
+ model_block,
+ paste(
+ gsub("d_(.*) =", "d/dt(\\1) =", mkin_model$diffs),
+ collapse = "\n"),
+ "\n"
+ )
+
+ # Error model
+ if (error_model == "const") {
+ model_block <- paste0(model_block,
+ paste(paste0(obs_vars, " ~ add(sigma)"), collapse = "\n"))
+ }
+ if (error_model == "obs") {
+ model_block <- paste0(model_block,
+ paste(paste0(obs_vars, " ~ add(sigma_", obs_vars, ")"), collapse = "\n"),
+ "\n")
+ }
+ if (error_model == "tc") {
+ model_block <- paste0(model_block,
+ paste(paste0(obs_vars, " ~ add(sigma_low) + prop(rsd_high)"), collapse = "\n"),
+ "\n")
+ }
+ if (error_model == "obs_tc") {
+ model_block <- paste0(model_block,
+ paste(
+ paste0(obs_vars, " ~ add(sigma_low_", obs_vars, ") + ",
+ "prop(rsd_high_", obs_vars, ")"), collapse = "\n"),
+ "\n")
+ }
+
+ model_block <- paste0(model_block, "})")
+
+ body(f)[3] <- parse(text = model_block)
+
+ if (add_attributes) {
+ attr(f, "mean_dp_start") <- degparms_optim
+ attr(f, "eta_start") <- degparms_mmkin$eta
+ attr(f, "mean_ep_start") <- errparms_ini
+ }
+
+ return(f)
+}
+
+#' @rdname nlmixr.mmkin
+#' @return An dataframe suitable for use with [nlmixr::nlmixr]
+#' @export
+nlmixr_data <- function(object, ...) {
+ if (nrow(object) > 1) stop("Only row objects allowed")
+ d <- lapply(object, function(x) x$data)
+ compartment_map <- 1:length(object[[1]]$mkinmod$spec)
+ names(compartment_map) <- names(object[[1]]$mkinmod$spec)
+ ds_names <- colnames(object)
+
+ ds_list <- lapply(object, function(x) x$data[c("time", "variable", "observed")])
+ names(ds_list) <- ds_names
+ ds_nlmixr <- purrr::map_dfr(ds_list, function(x) x, .id = "ds")
+ ds_nlmixr$variable <- as.character(ds_nlmixr$variable)
+ ds_nlmixr_renamed <- data.frame(
+ ID = ds_nlmixr$ds,
+ TIME = ds_nlmixr$time,
+ AMT = 0, EVID = 0,
+ CMT = ds_nlmixr$variable,
+ DV = ds_nlmixr$observed,
+ stringsAsFactors = FALSE)
+
+ return(ds_nlmixr_renamed)
+}
diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R
index 4c1f1531..1ac62b07 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 [mixed.mmkin], [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
@@ -10,6 +10,10 @@ utils::globalVariables("ds")
#' `resplot = "time"`.
#' @param pred_over Named list of alternative predictions as obtained
#' from [mkinpredict] with a compatible [mkinmod].
+#' @param test_log_parms Passed to [mean_degparms] in the case of an
+#' [mixed.mmkin] object
+#' @param conf.level Passed to [mean_degparms] in the case of an
+#' [mixed.mmkin] object
#' @param rel.height.legend The relative height of the legend shown on top
#' @param rel.height.bottom The relative height of the bottom plot row
#' @param ymax Vector of maximum y axis values
@@ -36,9 +40,23 @@ utils::globalVariables("ds")
#'
#' # For this fit we need to increase pnlsMaxiter, and we increase the
#' # tolerance in order to speed up the fit for this example evaluation
+#' # It still takes 20 seconds to run
#' f_nlme <- nlme(f, control = list(pnlsMaxIter = 120, tolerance = 1e-3))
#' plot(f_nlme)
#'
+#' f_saem <- saem(f, transformations = "saemix")
+#' plot(f_saem)
+#'
+#' f_obs <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, error_model = "obs")
+#' f_nlmix <- nlmix(f_obs)
+#' plot(f_nlmix)
+#'
+#' # We can overlay the two variants if we generate predictions
+#' pred_nlme <- mkinpredict(dfop_sfo,
+#' f_nlme$bparms.optim[-1],
+#' c(parent = f_nlme$bparms.optim[[1]], A1 = 0),
+#' seq(0, 180, by = 0.2))
+#' plot(f_saem, pred_over = list(nlme = pred_nlme))
#' }
#' @export
plot.mixed.mmkin <- function(x,
@@ -49,6 +67,8 @@ plot.mixed.mmkin <- function(x,
xlim = range(x$data$time),
resplot = c("predicted", "time"),
pred_over = NULL,
+ test_log_parms = FALSE,
+ conf.level = 0.6,
ymax = "auto", maxabs = "auto",
ncol.legend = ifelse(length(i) <= 3, length(i) + 1, ifelse(length(i) <= 8, 3, 4)),
nrow.legend = ceiling((length(i) + 1) / ncol.legend),
@@ -67,7 +87,7 @@ plot.mixed.mmkin <- function(x,
backtransform = TRUE
if (identical(class(x), "mixed.mmkin")) {
- degparms_pop <- mean_degparms(x$mmkin)
+ degparms_pop <- mean_degparms(x$mmkin, test_log_parms = test_log_parms, conf.level = conf.level)
degparms_tmp <- parms(x$mmkin, transformed = TRUE)
degparms_i <- as.data.frame(t(degparms_tmp[setdiff(rownames(degparms_tmp), names(fit_1$errparms)), ]))
@@ -82,6 +102,30 @@ plot.mixed.mmkin <- function(x,
type = ifelse(standardized, "pearson", "response"))
}
+ if (inherits(x, "saem.mmkin")) {
+ 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_i_names
+ }
+
+ if (inherits(x, "nlmixr.mmkin")) {
+ eta_i <- random.effects(x$nm)[-1]
+ names(eta_i) <- gsub("^eta.", "", names(eta_i))
+ degparms_i <- eta_i
+ degparms_pop <- x$nm$theta
+ for (parm_name in names(degparms_i)) {
+ degparms_i[parm_name] <- eta_i[parm_name] + degparms_pop[parm_name]
+ }
+ residual_type = ifelse(standardized, "standardized", "residual")
+ residuals <- x$data[[residual_type]]
+ }
+
degparms_fixed <- fit_1$fixed$value
names(degparms_fixed) <- rownames(fit_1$fixed)
degparms_all <- cbind(as.matrix(degparms_i),
diff --git a/R/saem.R b/R/saem.R
new file mode 100644
index 00000000..9db2c04a
--- /dev/null
+++ b/R/saem.R
@@ -0,0 +1,543 @@
+utils::globalVariables(c("predicted", "std"))
+
+#' Fit nonlinear mixed models with SAEM
+#'
+#' This function uses [saemix::saemix()] as a backend for fitting nonlinear mixed
+#' effects models created from [mmkin] row objects using the Stochastic Approximation
+#' Expectation Maximisation algorithm (SAEM).
+#'
+#' 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].
+#'
+#' @importFrom utils packageVersion
+#' @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 of
+#' type [saemix::SaemixModel] and [saemix::SaemixData]?
+#' @param transformations Per default, all parameter transformations are done
+#' in mkin. If this argument is set to 'saemix', parameter transformations
+#' are done in 'saemix' for the supported cases. Currently this is only
+#' supported in cases where the initial concentration of the parent is not fixed,
+#' SFO or DFOP is used for the parent and there is either no metabolite or one.
+#' @param degparms_start Parameter values given as a named numeric vector will
+#' be used to override the starting values obtained from the 'mmkin' object.
+#' @param test_log_parms If TRUE, an attempt is made to use more robust starting
+#' values for population parameters fitted as log parameters in mkin (like
+#' rate constants) by only considering rate constants that pass the t-test
+#' when calculating mean degradation parameters using [mean_degparms].
+#' @param conf.level Possibility to adjust the required confidence level
+#' for parameter that are tested if requested by 'test_log_parms'.
+#' @param solution_type Possibility to specify the solution type in case the
+#' automatic choice is not desired
+#' @param fail_with_errors Should a failure to compute standard errors
+#' from the inverse of the Fisher Information Matrix be a failure?
+#' @param quiet Should we suppress the messages saemix prints at the beginning
+#' and the end of the optimisation process?
+#' @param nbiter.saemix Convenience option to increase the number of
+#' iterations
+#' @param control Passed to [saemix::saemix].
+#' @param \dots Further parameters passed to [saemix::saemixModel].
+#' @return An S3 object of class 'saem.mmkin', containing the fitted
+#' [saemix::SaemixObject] as a list component named 'so'. The
+#' object also inherits from 'mixed.mmkin'.
+#' @seealso [summary.saem.mmkin] [plot.mixed.mmkin]
+#' @examples
+#' \dontrun{
+#' 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,
+#' state.ini = c(parent = 100), fixed_initials = "parent", quiet = TRUE)
+#' f_saem_p0_fixed <- saem(f_mmkin_parent_p0_fixed)
+#'
+#' f_mmkin_parent <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE)
+#' f_saem_sfo <- saem(f_mmkin_parent["SFO", ])
+#' f_saem_fomc <- saem(f_mmkin_parent["FOMC", ])
+#' f_saem_dfop <- saem(f_mmkin_parent["DFOP", ])
+#'
+#' # The returned saem.mmkin object contains an SaemixObject, therefore we can use
+#' # functions from saemix
+#' library(saemix)
+#' compare.saemix(f_saem_sfo$so, f_saem_fomc$so, f_saem_dfop$so)
+#' plot(f_saem_fomc$so, plot.type = "convergence")
+#' plot(f_saem_fomc$so, plot.type = "individual.fit")
+#' plot(f_saem_fomc$so, plot.type = "npde")
+#' plot(f_saem_fomc$so, plot.type = "vpc")
+#'
+#' f_mmkin_parent_tc <- update(f_mmkin_parent, error_model = "tc")
+#' f_saem_fomc_tc <- saem(f_mmkin_parent_tc["FOMC", ])
+#' compare.saemix(f_saem_fomc$so, f_saem_fomc_tc$so)
+#'
+#' sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"),
+#' A1 = mkinsub("SFO"))
+#' fomc_sfo <- mkinmod(parent = mkinsub("FOMC", "A1"),
+#' A1 = mkinsub("SFO"))
+#' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"),
+#' A1 = mkinsub("SFO"))
+#' # The following fit uses analytical solutions for SFO-SFO and DFOP-SFO,
+#' # and compiled ODEs for FOMC that are much slower
+#' f_mmkin <- mmkin(list(
+#' "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo),
+#' ds, quiet = TRUE)
+#' # saem fits of SFO-SFO and DFOP-SFO to these data take about five seconds
+#' # each on this system, as we use analytical solutions written for saemix.
+#' # When using the analytical solutions written for mkin this took around
+#' # four minutes
+#' f_saem_sfo_sfo <- saem(f_mmkin["SFO-SFO", ])
+#' f_saem_dfop_sfo <- saem(f_mmkin["DFOP-SFO", ])
+#' # We can use print, plot and summary methods to check the results
+#' print(f_saem_dfop_sfo)
+#' plot(f_saem_dfop_sfo)
+#' summary(f_saem_dfop_sfo, data = TRUE)
+#'
+#' # The following takes about 6 minutes
+#' #f_saem_dfop_sfo_deSolve <- saem(f_mmkin["DFOP-SFO", ], solution_type = "deSolve",
+#' # control = list(nbiter.saemix = c(200, 80), nbdisplay = 10))
+#'
+#' #saemix::compare.saemix(list(
+#' # f_saem_dfop_sfo$so,
+#' # f_saem_dfop_sfo_deSolve$so))
+#'
+#' # If the model supports it, we can also use eigenvalue based solutions, which
+#' # take a similar amount of time
+#' #f_saem_sfo_sfo_eigen <- saem(f_mmkin["SFO-SFO", ], solution_type = "eigen",
+#' # control = list(nbiter.saemix = c(200, 80), nbdisplay = 10))
+#' }
+#' @export
+saem <- function(object, ...) UseMethod("saem")
+
+#' @rdname saem
+#' @export
+saem.mmkin <- function(object,
+ transformations = c("mkin", "saemix"),
+ degparms_start = numeric(),
+ test_log_parms = TRUE,
+ conf.level = 0.6,
+ solution_type = "auto",
+ nbiter.saemix = c(300, 100),
+ control = list(displayProgress = FALSE, print = FALSE,
+ nbiter.saemix = nbiter.saemix,
+ save = FALSE, save.graphs = FALSE),
+ fail_with_errors = TRUE,
+ verbose = FALSE, quiet = FALSE, ...)
+{
+ transformations <- match.arg(transformations)
+ m_saemix <- saemix_model(object, verbose = verbose,
+ degparms_start = degparms_start,
+ test_log_parms = test_log_parms, conf.level = conf.level,
+ solution_type = solution_type,
+ transformations = transformations, ...)
+ d_saemix <- saemix_data(object, verbose = verbose)
+
+ fit_time <- system.time({
+ utils::capture.output(f_saemix <- saemix::saemix(m_saemix, d_saemix, control), split = !quiet)
+ FIM_failed <- NULL
+ if (any(is.na(f_saemix@results@se.fixed))) FIM_failed <- c(FIM_failed, "fixed effects")
+ if (any(is.na(c(f_saemix@results@se.omega, f_saemix@results@se.respar)))) {
+ FIM_failed <- c(FIM_failed, "random effects and residual error parameters")
+ }
+ if (!is.null(FIM_failed) & fail_with_errors) {
+ stop("Could not invert FIM for ", paste(FIM_failed, collapse = " and "))
+ }
+ })
+
+ transparms_optim <- f_saemix@results@fixed.effects
+ names(transparms_optim) <- f_saemix@results@name.fixed
+
+ if (transformations == "mkin") {
+ bparms_optim <- backtransform_odeparms(transparms_optim,
+ object[[1]]$mkinmod,
+ object[[1]]$transform_rates,
+ object[[1]]$transform_fractions)
+ } else {
+ bparms_optim <- transparms_optim
+ }
+
+ return_data <- nlme_data(object)
+
+ return_data$predicted <- f_saemix@model@model(
+ psi = saemix::psi(f_saemix),
+ id = as.numeric(return_data$ds),
+ xidep = return_data[c("time", "name")])
+
+ return_data <- transform(return_data,
+ residual = value - predicted,
+ std = sigma_twocomp(predicted,
+ f_saemix@results@respar[1], f_saemix@results@respar[2]))
+ return_data <- transform(return_data,
+ standardized = residual / std)
+
+ result <- list(
+ mkinmod = object[[1]]$mkinmod,
+ mmkin = object,
+ solution_type = object[[1]]$solution_type,
+ transformations = transformations,
+ transform_rates = object[[1]]$transform_rates,
+ transform_fractions = object[[1]]$transform_fractions,
+ so = f_saemix,
+ time = fit_time,
+ mean_dp_start = attr(m_saemix, "mean_dp_start"),
+ bparms.optim = bparms_optim,
+ bparms.fixed = object[[1]]$bparms.fixed,
+ data = return_data,
+ err_mod = object[[1]]$err_mod,
+ date.fit = date(),
+ saemixversion = as.character(utils::packageVersion("saemix")),
+ mkinversion = as.character(utils::packageVersion("mkin")),
+ Rversion = paste(R.version$major, R.version$minor, sep=".")
+ )
+
+ class(result) <- c("saem.mmkin", "mixed.mmkin")
+ return(result)
+}
+
+#' @export
+#' @rdname saem
+#' @param x An saem.mmkin object to print
+#' @param digits Number of digits to use for printing
+print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) {
+ cat( "Kinetic nonlinear mixed-effects model fit by SAEM" )
+ 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")
+
+ cat("\nLikelihood computed by importance sampling\n")
+ print(data.frame(
+ AIC = AIC(x$so, type = "is"),
+ BIC = BIC(x$so, type = "is"),
+ logLik = logLik(x$so, type = "is"),
+ row.names = " "), digits = digits)
+
+ cat("\nFitted parameters:\n")
+ conf.int <- x$so@results@conf.int[c("estimate", "lower", "upper")]
+ rownames(conf.int) <- x$so@results@conf.int[["name"]]
+ print(conf.int, digits = digits)
+
+ invisible(x)
+}
+
+#' @rdname saem
+#' @return An [saemix::SaemixModel] object.
+#' @export
+saemix_model <- function(object, solution_type = "auto", transformations = c("mkin", "saemix"),
+ degparms_start = numeric(), test_log_parms = FALSE, verbose = FALSE, ...)
+{
+ if (packageVersion("saemix") < "3.1.9000") {
+ stop("To use the interface to saemix, you need to install a development version\n",
+ "preferably https://github.com/jranke/saemixextension@warp_combined")
+ }
+
+ if (nrow(object) > 1) stop("Only row objects allowed")
+
+ mkin_model <- object[[1]]$mkinmod
+
+ degparms_optim <- mean_degparms(object, test_log_parms = test_log_parms)
+ if (transformations == "saemix") {
+ degparms_optim <- backtransform_odeparms(degparms_optim,
+ object[[1]]$mkinmod,
+ object[[1]]$transform_rates,
+ object[[1]]$transform_fractions)
+ }
+ degparms_fixed <- object[[1]]$bparms.fixed
+
+ # Transformations are done in the degradation function
+ transform.par = rep(0, length(degparms_optim))
+
+ odeini_optim_parm_names <- grep('_0$', names(degparms_optim), value = TRUE)
+ odeini_fixed_parm_names <- grep('_0$', names(degparms_fixed), value = TRUE)
+
+ odeparms_fixed_names <- setdiff(names(degparms_fixed), odeini_fixed_parm_names)
+ odeparms_fixed <- degparms_fixed[odeparms_fixed_names]
+
+ odeini_fixed <- degparms_fixed[odeini_fixed_parm_names]
+ names(odeini_fixed) <- gsub('_0$', '', odeini_fixed_parm_names)
+
+ model_function <- FALSE
+
+ # Model functions with analytical solutions
+ # Fixed parameters, use_of_ff = "min" and turning off sinks currently not supported here
+ # In general, we need to consider exactly how the parameters in mkinfit were specified,
+ # as the parameters are currently mapped by position in these solutions
+ sinks <- sapply(mkin_model$spec, function(x) x$sink)
+ if (length(odeparms_fixed) == 0 & mkin_model$use_of_ff == "max" & all(sinks)) {
+ # Parent only
+ if (length(mkin_model$spec) == 1) {
+ parent_type <- mkin_model$spec[[1]]$type
+ if (length(odeini_fixed) == 1) {
+ if (parent_type == "SFO") {
+ stop("saemix needs at least two parameters to work on.")
+ }
+ if (parent_type == "FOMC") {
+ model_function <- function(psi, id, xidep) {
+ odeini_fixed / (xidep[, "time"]/exp(psi[id, 2]) + 1)^exp(psi[id, 1])
+ }
+ }
+ if (parent_type == "DFOP") {
+ model_function <- function(psi, id, xidep) {
+ g <- plogis(psi[id, 3])
+ t <- xidep[, "time"]
+ odeini_fixed * (g * exp(- exp(psi[id, 1]) * t) +
+ (1 - g) * exp(- exp(psi[id, 2]) * t))
+ }
+ }
+ if (parent_type == "HS") {
+ model_function <- function(psi, id, xidep) {
+ tb <- exp(psi[id, 3])
+ t <- xidep[, "time"]
+ k1 = exp(psi[id, 1])
+ odeini_fixed * ifelse(t <= tb,
+ exp(- k1 * t),
+ exp(- k1 * tb) * exp(- exp(psi[id, 2]) * (t - tb)))
+ }
+ }
+ } else {
+ if (parent_type == "SFO") {
+ if (transformations == "mkin") {
+ model_function <- function(psi, id, xidep) {
+ psi[id, 1] * exp( - exp(psi[id, 2]) * xidep[, "time"])
+ }
+ } else {
+ model_function <- function(psi, id, xidep) {
+ psi[id, 1] * exp( - psi[id, 2] * xidep[, "time"])
+ }
+ transform.par = c(0, 1)
+ }
+ }
+ if (parent_type == "FOMC") {
+ model_function <- function(psi, id, xidep) {
+ psi[id, 1] / (xidep[, "time"]/exp(psi[id, 3]) + 1)^exp(psi[id, 2])
+ }
+ }
+ if (parent_type == "DFOP") {
+ if (transformations == "mkin") {
+ model_function <- function(psi, id, xidep) {
+ g <- plogis(psi[id, 4])
+ t <- xidep[, "time"]
+ psi[id, 1] * (g * exp(- exp(psi[id, 2]) * t) +
+ (1 - g) * exp(- exp(psi[id, 3]) * t))
+ }
+ } else {
+ model_function <- function(psi, id, xidep) {
+ g <- psi[id, 4]
+ t <- xidep[, "time"]
+ psi[id, 1] * (g * exp(- psi[id, 2] * t) +
+ (1 - g) * exp(- psi[id, 3] * t))
+ }
+ transform.par = c(0, 1, 1, 3)
+ }
+ }
+ if (parent_type == "HS") {
+ model_function <- function(psi, id, xidep) {
+ tb <- exp(psi[id, 4])
+ t <- xidep[, "time"]
+ k1 = exp(psi[id, 2])
+ psi[id, 1] * ifelse(t <= tb,
+ exp(- k1 * t),
+ exp(- k1 * tb) * exp(- exp(psi[id, 3]) * (t - tb)))
+ }
+ }
+ }
+ }
+
+ # Parent with one metabolite
+ # Parameter names used in the model functions are as in
+ # https://nbviewer.jupyter.org/urls/jrwb.de/nb/Symbolic%20ODE%20solutions%20for%20mkin.ipynb
+ types <- unname(sapply(mkin_model$spec, function(x) x$type))
+ if (length(mkin_model$spec) == 2 &! "SFORB" %in% types ) {
+ # Initial value for the metabolite (n20) must be fixed
+ if (names(odeini_fixed) == names(mkin_model$spec)[2]) {
+ n20 <- odeini_fixed
+ parent_name <- names(mkin_model$spec)[1]
+ if (identical(types, c("SFO", "SFO"))) {
+ if (transformations == "mkin") {
+ model_function <- function(psi, id, xidep) {
+ t <- xidep[, "time"]
+ n10 <- psi[id, 1]
+ k1 <- exp(psi[id, 2])
+ k2 <- exp(psi[id, 3])
+ f12 <- plogis(psi[id, 4])
+ ifelse(xidep[, "name"] == parent_name,
+ n10 * exp(- k1 * t),
+ (((k2 - k1) * n20 - f12 * k1 * n10) * exp(- k2 * t)) / (k2 - k1) +
+ (f12 * k1 * n10 * exp(- k1 * t)) / (k2 - k1)
+ )
+ }
+ } else {
+ model_function <- function(psi, id, xidep) {
+ t <- xidep[, "time"]
+ n10 <- psi[id, 1]
+ k1 <- psi[id, 2]
+ k2 <- psi[id, 3]
+ f12 <- psi[id, 4]
+ ifelse(xidep[, "name"] == parent_name,
+ n10 * exp(- k1 * t),
+ (((k2 - k1) * n20 - f12 * k1 * n10) * exp(- k2 * t)) / (k2 - k1) +
+ (f12 * k1 * n10 * exp(- k1 * t)) / (k2 - k1)
+ )
+ }
+ transform.par = c(0, 1, 1, 3)
+ }
+ }
+ if (identical(types, c("DFOP", "SFO"))) {
+ if (transformations == "mkin") {
+ model_function <- function(psi, id, xidep) {
+ t <- xidep[, "time"]
+ n10 <- psi[id, 1]
+ k2 <- exp(psi[id, 2])
+ f12 <- plogis(psi[id, 3])
+ l1 <- exp(psi[id, 4])
+ l2 <- exp(psi[id, 5])
+ g <- plogis(psi[id, 6])
+ ifelse(xidep[, "name"] == parent_name,
+ n10 * (g * exp(- l1 * t) + (1 - g) * exp(- l2 * t)),
+ ((f12 * g - f12) * l2 * n10 * exp(- l2 * t)) / (l2 - k2) -
+ (f12 * g * l1 * n10 * exp(- l1 * t)) / (l1 - k2) +
+ ((((l1 - k2) * l2 - k2 * l1 + k2^2) * n20 +
+ ((f12 * l1 + (f12 * g - f12) * k2) * l2 -
+ f12 * g * k2 * l1) * n10) * exp( - k2 * t)) /
+ ((l1 - k2) * l2 - k2 * l1 + k2^2)
+ )
+ }
+ } else {
+ model_function <- function(psi, id, xidep) {
+ t <- xidep[, "time"]
+ n10 <- psi[id, 1]
+ k2 <- psi[id, 2]
+ f12 <- psi[id, 3]
+ l1 <- psi[id, 4]
+ l2 <- psi[id, 5]
+ g <- psi[id, 6]
+ ifelse(xidep[, "name"] == parent_name,
+ n10 * (g * exp(- l1 * t) + (1 - g) * exp(- l2 * t)),
+ ((f12 * g - f12) * l2 * n10 * exp(- l2 * t)) / (l2 - k2) -
+ (f12 * g * l1 * n10 * exp(- l1 * t)) / (l1 - k2) +
+ ((((l1 - k2) * l2 - k2 * l1 + k2^2) * n20 +
+ ((f12 * l1 + (f12 * g - f12) * k2) * l2 -
+ f12 * g * k2 * l1) * n10) * exp( - k2 * t)) /
+ ((l1 - k2) * l2 - k2 * l1 + k2^2)
+ )
+ }
+ transform.par = c(0, 1, 3, 1, 1, 3)
+ }
+ }
+ }
+ }
+ }
+
+ if (is.function(model_function) & solution_type == "auto") {
+ solution_type = "analytical saemix"
+ } else {
+
+ if (solution_type == "auto")
+ solution_type <- object[[1]]$solution_type
+
+ model_function <- function(psi, id, xidep) {
+
+ uid <- unique(id)
+
+ res_list <- lapply(uid, function(i) {
+
+ transparms_optim <- as.numeric(psi[i, ]) # psi[i, ] is a dataframe when called in saemix.predict
+ names(transparms_optim) <- names(degparms_optim)
+
+ odeini_optim <- transparms_optim[odeini_optim_parm_names]
+ names(odeini_optim) <- gsub('_0$', '', odeini_optim_parm_names)
+
+ odeini <- c(odeini_optim, odeini_fixed)[names(mkin_model$diffs)]
+
+ ode_transparms_optim_names <- setdiff(names(transparms_optim), odeini_optim_parm_names)
+ odeparms_optim <- backtransform_odeparms(transparms_optim[ode_transparms_optim_names], mkin_model,
+ transform_rates = object[[1]]$transform_rates,
+ transform_fractions = object[[1]]$transform_fractions)
+ odeparms <- c(odeparms_optim, odeparms_fixed)
+
+ xidep_i <- subset(xidep, id == i)
+
+ if (solution_type == "analytical") {
+ out_values <- mkin_model$deg_func(xidep_i, odeini, odeparms)
+ } else {
+
+ i_time <- xidep_i$time
+ i_name <- xidep_i$name
+
+ out_wide <- mkinpredict(mkin_model,
+ odeparms = odeparms, odeini = odeini,
+ solution_type = solution_type,
+ outtimes = sort(unique(i_time)),
+ na_stop = FALSE
+ )
+
+ out_index <- cbind(as.character(i_time), as.character(i_name))
+ out_values <- out_wide[out_index]
+ }
+ return(out_values)
+ })
+ res <- unlist(res_list)
+ return(res)
+ }
+ }
+
+ error.model <- switch(object[[1]]$err_mod,
+ const = "constant",
+ tc = "combined",
+ obs = "constant")
+
+ if (object[[1]]$err_mod == "obs") {
+ warning("The error model 'obs' (variance by variable) can currently not be transferred to an saemix model")
+ }
+
+ error.init <- switch(object[[1]]$err_mod,
+ const = c(a = mean(sapply(object, function(x) x$errparms)), b = 1),
+ tc = c(a = mean(sapply(object, function(x) x$errparms[1])),
+ b = mean(sapply(object, function(x) x$errparms[2]))),
+ obs = c(a = mean(sapply(object, function(x) x$errparms)), b = 1))
+
+ degparms_psi0 <- degparms_optim
+ degparms_psi0[names(degparms_start)] <- degparms_start
+ psi0_matrix <- matrix(degparms_psi0, nrow = 1)
+ colnames(psi0_matrix) <- names(degparms_psi0)
+
+ res <- saemix::saemixModel(model_function,
+ psi0 = psi0_matrix,
+ "Mixed model generated from mmkin object",
+ transform.par = transform.par,
+ error.model = error.model,
+ verbose = verbose
+ )
+ attr(res, "mean_dp_start") <- degparms_optim
+ return(res)
+}
+
+#' @rdname saem
+#' @return An [saemix::SaemixData] object.
+#' @export
+saemix_data <- function(object, verbose = FALSE, ...) {
+ if (nrow(object) > 1) stop("Only row objects allowed")
+ ds_names <- colnames(object)
+
+ ds_list <- lapply(object, function(x) x$data[c("time", "variable", "observed")])
+ names(ds_list) <- ds_names
+ ds_saemix_all <- purrr::map_dfr(ds_list, function(x) x, .id = "ds")
+ ds_saemix <- data.frame(ds = ds_saemix_all$ds,
+ name = as.character(ds_saemix_all$variable),
+ time = ds_saemix_all$time,
+ value = ds_saemix_all$observed,
+ stringsAsFactors = FALSE)
+
+ res <- saemix::saemixData(ds_saemix,
+ name.group = "ds",
+ name.predictors = c("time", "name"),
+ name.response = "value",
+ verbose = verbose,
+ ...)
+ return(res)
+}
diff --git a/R/summary.nlmixr.mmkin.R b/R/summary.nlmixr.mmkin.R
new file mode 100644
index 00000000..a023f319
--- /dev/null
+++ b/R/summary.nlmixr.mmkin.R
@@ -0,0 +1,250 @@
+#' Summary method for class "nlmixr.mmkin"
+#'
+#' Lists model equations, initial parameter values, optimised parameters
+#' for fixed effects (population), random effects (deviations from the
+#' population mean) and residual error model, as well as the resulting
+#' endpoints such as formation fractions and DT50 values. Optionally
+#' (default is FALSE), the data are listed in full.
+#'
+#' @importFrom stats confint sd
+#' @param object an object of class [nlmixr.mmkin]
+#' @param x an object of class [summary.nlmixr.mmkin]
+#' @param data logical, indicating whether the full data should be included in
+#' the summary.
+#' @param verbose Should the summary be verbose?
+#' @param distimes logical, indicating whether DT50 and DT90 values should be
+#' included.
+#' @param digits Number of digits to use for printing
+#' @param \dots optional arguments passed to methods like \code{print}.
+#' @return The summary function returns a list obtained in the fit, with at
+#' least the following additional components
+#' \item{nlmixrversion, mkinversion, Rversion}{The nlmixr, mkin and R versions used}
+#' \item{date.fit, date.summary}{The dates where the fit and the summary were
+#' produced}
+#' \item{diffs}{The differential equations used in the degradation model}
+#' \item{use_of_ff}{Was maximum or minimum use made of formation fractions}
+#' \item{data}{The data}
+#' \item{confint_back}{Backtransformed parameters, with confidence intervals if available}
+#' \item{ff}{The estimated formation fractions derived from the fitted
+#' model.}
+#' \item{distimes}{The DT50 and DT90 values for each observed variable.}
+#' \item{SFORB}{If applicable, eigenvalues of SFORB components of the model.}
+#' The print method is called for its side effect, i.e. printing the summary.
+#' @importFrom stats predict vcov
+#' @author Johannes Ranke for the mkin specific parts
+#' nlmixr authors for the parts inherited from nlmixr.
+#' @examples
+#' # Generate five datasets following DFOP-SFO kinetics
+#' sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
+#' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "m1"),
+#' m1 = mkinsub("SFO"), quiet = TRUE)
+#' set.seed(1234)
+#' k1_in <- rlnorm(5, log(0.1), 0.3)
+#' k2_in <- rlnorm(5, log(0.02), 0.3)
+#' g_in <- plogis(rnorm(5, qlogis(0.5), 0.3))
+#' f_parent_to_m1_in <- plogis(rnorm(5, qlogis(0.3), 0.3))
+#' k_m1_in <- rlnorm(5, log(0.02), 0.3)
+#'
+#' pred_dfop_sfo <- function(k1, k2, g, f_parent_to_m1, k_m1) {
+#' mkinpredict(dfop_sfo,
+#' c(k1 = k1, k2 = k2, g = g, f_parent_to_m1 = f_parent_to_m1, k_m1 = k_m1),
+#' c(parent = 100, m1 = 0),
+#' sampling_times)
+#' }
+#'
+#' ds_mean_dfop_sfo <- lapply(1:5, function(i) {
+#' mkinpredict(dfop_sfo,
+#' c(k1 = k1_in[i], k2 = k2_in[i], g = g_in[i],
+#' f_parent_to_m1 = f_parent_to_m1_in[i], k_m1 = k_m1_in[i]),
+#' c(parent = 100, m1 = 0),
+#' sampling_times)
+#' })
+#' names(ds_mean_dfop_sfo) <- paste("ds", 1:5)
+#'
+#' ds_syn_dfop_sfo <- lapply(ds_mean_dfop_sfo, function(ds) {
+#' add_err(ds,
+#' sdfunc = function(value) sqrt(1^2 + value^2 * 0.07^2),
+#' n = 1)[[1]]
+#' })
+#'
+#' \dontrun{
+#' # Evaluate using mmkin and nlmixr
+#' f_mmkin_dfop_sfo <- mmkin(list(dfop_sfo), ds_syn_dfop_sfo,
+#' quiet = TRUE, error_model = "tc", cores = 5)
+#' f_saemix_dfop_sfo <- mkin::saem(f_mmkin_dfop_sfo)
+#' f_nlme_dfop_sfo <- mkin::nlme(f_mmkin_dfop_sfo)
+#' f_nlmixr_dfop_sfo_saem <- nlmixr(f_mmkin_dfop_sfo, est = "saem")
+#' # The following takes a very long time but gives
+#' f_nlmixr_dfop_sfo_focei <- nlmixr(f_mmkin_dfop_sfo, est = "focei")
+#' AIC(f_nlmixr_dfop_sfo_saem$nm, f_nlmixr_dfop_sfo_focei$nm)
+#' summary(f_nlmixr_dfop_sfo_sfo, data = TRUE)
+#' }
+#'
+#' @export
+summary.nlmixr.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes = TRUE, ...) {
+
+ mod_vars <- names(object$mkinmod$diffs)
+
+ conf.int <- confint(object$nm)
+ dpnames <- setdiff(rownames(conf.int), names(object$mean_ep_start))
+ ndp <- length(dpnames)
+
+ confint_trans <- as.matrix(conf.int[dpnames, c(1, 3, 4)])
+ colnames(confint_trans) <- c("est.", "lower", "upper")
+
+ 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 = "_"), dpnames, 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 dpnames) {
+ 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
+ }
+ }
+
+ # Correlation of fixed effects (inspired by summary.nlme)
+ varFix <- vcov(object$nm)
+ stdFix <- sqrt(diag(varFix))
+ object$corFixed <- array(
+ t(varFix/stdFix)/stdFix,
+ dim(varFix),
+ list(dpnames, dpnames))
+
+ object$confint_trans <- confint_trans
+ object$confint_back <- confint_back
+
+ object$date.summary = date()
+ object$use_of_ff = object$mkinmod$use_of_ff
+
+ object$diffs <- object$mkinmod$diffs
+ object$print_data <- data # boolean: Should we print the data?
+
+ names(object$data)[4] <- "observed" # rename value to observed
+
+ object$verbose <- verbose
+
+ object$fixed <- object$mmkin_orig[[1]]$fixed
+ object$AIC = AIC(object$nm)
+ object$BIC = BIC(object$nm)
+ object$logLik = logLik(object$nm)
+
+ ep <- endpoints(object)
+ if (length(ep$ff) != 0)
+ object$ff <- ep$ff
+ if (distimes) object$distimes <- ep$distimes
+ if (length(ep$SFORB) != 0) object$SFORB <- ep$SFORB
+ class(object) <- c("summary.nlmixr.mmkin")
+ return(object)
+}
+
+#' @rdname summary.nlmixr.mmkin
+#' @export
+print.summary.nlmixr.mmkin <- function(x, digits = max(3, getOption("digits") - 3), verbose = x$verbose, ...) {
+ cat("nlmixr version used for fitting: ", x$nlmixrversion, "\n")
+ cat("mkin version used for pre-fitting: ", x$mkinversion, "\n")
+ cat("R version used for fitting: ", x$Rversion, "\n")
+
+ cat("Date of fit: ", x$date.fit, "\n")
+ cat("Date of summary:", x$date.summary, "\n")
+
+ cat("\nEquations:\n")
+ nice_diffs <- gsub("^(d.*) =", "\\1/dt =", x[["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")
+
+ cat("\nDegradation model predictions using RxODE\n")
+
+ cat("\nFitted in", x$time[["elapsed"]], "s\n")
+
+ cat("\nVariance model: ")
+ cat(switch(x$err_mod,
+ const = "Constant variance",
+ obs = "Variance unique to each observed variable",
+ tc = "Two-component variance function",
+ obs_tc = "Two-component variance unique to each observed variable"), "\n")
+
+ cat("\nMean of starting values for individual parameters:\n")
+ print(x$mean_dp_start, digits = digits)
+
+ cat("\nMean of starting values for error model parameters:\n")
+ print(x$mean_ep_start, digits = digits)
+
+ cat("\nFixed degradation parameter values:\n")
+ if(length(x$fixed$value) == 0) cat("None\n")
+ else print(x$fixed, digits = digits)
+
+ cat("\nResults:\n\n")
+ cat("Likelihood calculated by", nlmixr::getOfvType(x$nm), " \n")
+ print(data.frame(AIC = x$AIC, BIC = x$BIC, logLik = x$logLik,
+ row.names = " "), digits = digits)
+
+ cat("\nOptimised parameters:\n")
+ print(x$confint_trans, digits = digits)
+
+ if (nrow(x$confint_trans) > 1) {
+ corr <- x$corFixed
+ class(corr) <- "correlation"
+ print(corr, title = "\nCorrelation:", ...)
+ }
+
+ cat("\nRandom effects (omega):\n")
+ print(x$nm$omega, digits = digits)
+
+ cat("\nVariance model:\n")
+ print(x$nm$sigma, digits = digits)
+
+ cat("\nBacktransformed parameters:\n")
+ print(x$confint_back, digits = digits)
+
+ printSFORB <- !is.null(x$SFORB)
+ if(printSFORB){
+ cat("\nEstimated Eigenvalues of SFORB model(s):\n")
+ print(x$SFORB, digits = digits,...)
+ }
+
+ printff <- !is.null(x$ff)
+ if(printff){
+ cat("\nResulting formation fractions:\n")
+ print(data.frame(ff = x$ff), digits = digits,...)
+ }
+
+ printdistimes <- !is.null(x$distimes)
+ if(printdistimes){
+ cat("\nEstimated disappearance times:\n")
+ print(x$distimes, digits = digits,...)
+ }
+
+ if (x$print_data){
+ cat("\nData:\n")
+ print(format(x$data, digits = digits, ...), row.names = FALSE)
+ }
+
+ invisible(x)
+}
diff --git a/R/summary.saem.mmkin.R b/R/summary.saem.mmkin.R
new file mode 100644
index 00000000..e92c561c
--- /dev/null
+++ b/R/summary.saem.mmkin.R
@@ -0,0 +1,268 @@
+#' Summary method for class "saem.mmkin"
+#'
+#' Lists model equations, initial parameter values, optimised parameters
+#' for fixed effects (population), random effects (deviations from the
+#' population mean) and residual error model, as well as the resulting
+#' endpoints such as formation fractions and DT50 values. Optionally
+#' (default is FALSE), the data are listed in full.
+#'
+#' @param object an object of class [saem.mmkin]
+#' @param x an object of class [summary.saem.mmkin]
+#' @param data logical, indicating whether the full data should be included in
+#' the summary.
+#' @param verbose Should the summary be verbose?
+#' @param distimes logical, indicating whether DT50 and DT90 values should be
+#' included.
+#' @param digits Number of digits to use for printing
+#' @param \dots optional arguments passed to methods like \code{print}.
+#' @return The summary function returns a list based on the [saemix::SaemixObject]
+#' obtained in the fit, with at least the following additional components
+#' \item{saemixversion, mkinversion, Rversion}{The saemix, mkin and R versions used}
+#' \item{date.fit, date.summary}{The dates where the fit and the summary were
+#' produced}
+#' \item{diffs}{The differential equations used in the degradation model}
+#' \item{use_of_ff}{Was maximum or minimum use made of formation fractions}
+#' \item{data}{The data}
+#' \item{confint_trans}{Transformed parameters as used in the optimisation, with confidence intervals}
+#' \item{confint_back}{Backtransformed parameters, with confidence intervals if available}
+#' \item{confint_errmod}{Error model parameters with confidence intervals}
+#' \item{ff}{The estimated formation fractions derived from the fitted
+#' model.}
+#' \item{distimes}{The DT50 and DT90 values for each observed variable.}
+#' \item{SFORB}{If applicable, eigenvalues of SFORB components of the model.}
+#' The print method is called for its side effect, i.e. printing the summary.
+#' @importFrom stats predict vcov
+#' @author Johannes Ranke for the mkin specific parts
+#' saemix authors for the parts inherited from saemix.
+#' @examples
+#' # Generate five datasets following DFOP-SFO kinetics
+#' sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
+#' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "m1"),
+#' m1 = mkinsub("SFO"), quiet = TRUE)
+#' set.seed(1234)
+#' k1_in <- rlnorm(5, log(0.1), 0.3)
+#' k2_in <- rlnorm(5, log(0.02), 0.3)
+#' g_in <- plogis(rnorm(5, qlogis(0.5), 0.3))
+#' f_parent_to_m1_in <- plogis(rnorm(5, qlogis(0.3), 0.3))
+#' k_m1_in <- rlnorm(5, log(0.02), 0.3)
+#'
+#' pred_dfop_sfo <- function(k1, k2, g, f_parent_to_m1, k_m1) {
+#' mkinpredict(dfop_sfo,
+#' c(k1 = k1, k2 = k2, g = g, f_parent_to_m1 = f_parent_to_m1, k_m1 = k_m1),
+#' c(parent = 100, m1 = 0),
+#' sampling_times)
+#' }
+#'
+#' ds_mean_dfop_sfo <- lapply(1:5, function(i) {
+#' mkinpredict(dfop_sfo,
+#' c(k1 = k1_in[i], k2 = k2_in[i], g = g_in[i],
+#' f_parent_to_m1 = f_parent_to_m1_in[i], k_m1 = k_m1_in[i]),
+#' c(parent = 100, m1 = 0),
+#' sampling_times)
+#' })
+#' names(ds_mean_dfop_sfo) <- paste("ds", 1:5)
+#'
+#' ds_syn_dfop_sfo <- lapply(ds_mean_dfop_sfo, function(ds) {
+#' add_err(ds,
+#' sdfunc = function(value) sqrt(1^2 + value^2 * 0.07^2),
+#' n = 1)[[1]]
+#' })
+#'
+#' \dontrun{
+#' # Evaluate using mmkin and saem
+#' f_mmkin_dfop_sfo <- mmkin(list(dfop_sfo), ds_syn_dfop_sfo,
+#' quiet = TRUE, error_model = "tc", cores = 5)
+#' f_saem_dfop_sfo <- saem(f_mmkin_dfop_sfo)
+#' summary(f_saem_dfop_sfo, data = TRUE)
+#' }
+#'
+#' @export
+summary.saem.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes = TRUE, ...) {
+
+ mod_vars <- names(object$mkinmod$diffs)
+
+ pnames <- names(object$mean_dp_start)
+ np <- length(pnames)
+
+ conf.int <- object$so@results@conf.int
+ rownames(conf.int) <- conf.int$name
+ confint_trans <- as.matrix(conf.int[pnames, c("estimate", "lower", "upper")])
+ colnames(confint_trans)[1] <- "est."
+
+ # In case objects were produced by earlier versions of saem
+ if (is.null(object$transformations)) object$transformations <- "mkin"
+
+ 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]
+ stdFix <- sqrt(diag(varFix))
+ object$corFixed <- array(
+ t(varFix/stdFix)/stdFix,
+ dim(varFix),
+ list(pnames, pnames))
+
+ # Random effects
+ rnames <- paste0("SD.", pnames)
+ confint_ranef <- as.matrix(conf.int[rnames, c("estimate", "lower", "upper")])
+ colnames(confint_ranef)[1] <- "est."
+
+ # Error model
+ enames <- if (object$err_mod == "const") "a.1" else c("a.1", "b.1")
+ confint_errmod <- as.matrix(conf.int[enames, c("estimate", "lower", "upper")])
+ colnames(confint_errmod)[1] <- "est."
+
+
+ object$confint_trans <- confint_trans
+ object$confint_ranef <- confint_ranef
+ object$confint_errmod <- confint_errmod
+ object$confint_back <- confint_back
+
+ object$date.summary = date()
+ object$use_of_ff = object$mkinmod$use_of_ff
+ object$error_model_algorithm = object$mmkin_orig[[1]]$error_model_algorithm
+ err_mod = object$mmkin_orig[[1]]$err_mod
+
+ object$diffs <- object$mkinmod$diffs
+ object$print_data <- data # boolean: Should we print the data?
+ so_pred <- object$so@results@predictions
+
+ names(object$data)[4] <- "observed" # rename value to observed
+
+ object$verbose <- verbose
+
+ object$fixed <- object$mmkin_orig[[1]]$fixed
+ object$AIC = AIC(object$so)
+ object$BIC = BIC(object$so)
+ object$logLik = logLik(object$so, method = "is")
+
+ ep <- endpoints(object)
+ if (length(ep$ff) != 0)
+ object$ff <- ep$ff
+ if (distimes) object$distimes <- ep$distimes
+ if (length(ep$SFORB) != 0) object$SFORB <- ep$SFORB
+ class(object) <- c("summary.saem.mmkin")
+ return(object)
+}
+
+#' @rdname summary.saem.mmkin
+#' @export
+print.summary.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), verbose = x$verbose, ...) {
+ cat("saemix version used for fitting: ", x$saemixversion, "\n")
+ cat("mkin version used for pre-fitting: ", x$mkinversion, "\n")
+ cat("R version used for fitting: ", x$Rversion, "\n")
+
+ cat("Date of fit: ", x$date.fit, "\n")
+ cat("Date of summary:", x$date.summary, "\n")
+
+ cat("\nEquations:\n")
+ nice_diffs <- gsub("^(d.*) =", "\\1/dt =", x[["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")
+
+ cat("\nModel predictions using solution type", x$solution_type, "\n")
+
+ cat("\nFitted in", x$time[["elapsed"]], "s using", paste(x$so@options$nbiter.saemix, collapse = ", "), "iterations\n")
+
+ cat("\nVariance model: ")
+ cat(switch(x$err_mod,
+ const = "Constant variance",
+ obs = "Variance unique to each observed variable",
+ tc = "Two-component variance function"), "\n")
+
+ cat("\nMean of starting values for individual parameters:\n")
+ print(x$mean_dp_start, digits = digits)
+
+ cat("\nFixed degradation parameter values:\n")
+ if(length(x$fixed$value) == 0) cat("None\n")
+ else print(x$fixed, digits = digits)
+
+ cat("\nResults:\n\n")
+ cat("Likelihood computed by importance sampling\n")
+ print(data.frame(AIC = x$AIC, BIC = x$BIC, logLik = x$logLik,
+ row.names = " "), digits = digits)
+
+ cat("\nOptimised parameters:\n")
+ print(x$confint_trans, digits = digits)
+
+ if (nrow(x$confint_trans) > 1) {
+ corr <- x$corFixed
+ class(corr) <- "correlation"
+ print(corr, title = "\nCorrelation:", ...)
+ }
+
+ cat("\nRandom effects:\n")
+ print(x$confint_ranef, digits = digits)
+
+ cat("\nVariance model:\n")
+ print(x$confint_errmod, digits = digits)
+
+ if (x$transformations == "mkin") {
+ cat("\nBacktransformed parameters:\n")
+ print(x$confint_back, digits = digits)
+ }
+
+ printSFORB <- !is.null(x$SFORB)
+ if(printSFORB){
+ cat("\nEstimated Eigenvalues of SFORB model(s):\n")
+ print(x$SFORB, digits = digits,...)
+ }
+
+ printff <- !is.null(x$ff)
+ if(printff){
+ cat("\nResulting formation fractions:\n")
+ print(data.frame(ff = x$ff), digits = digits,...)
+ }
+
+ printdistimes <- !is.null(x$distimes)
+ if(printdistimes){
+ cat("\nEstimated disappearance times:\n")
+ print(x$distimes, digits = digits,...)
+ }
+
+ if (x$print_data){
+ cat("\nData:\n")
+ print(format(x$data, digits = digits, ...), row.names = FALSE)
+ }
+
+ invisible(x)
+}
diff --git a/R/tffm0.R b/R/tffm0.R
new file mode 100644
index 00000000..bb5f4cf5
--- /dev/null
+++ b/R/tffm0.R
@@ -0,0 +1,48 @@
+#' Transform formation fractions as in the first published mkin version
+#'
+#' The transformed fractions can be restricted between 0 and 1 in model
+#' optimisations. Therefore this transformation was used originally in mkin. It
+#' was later replaced by the [ilr] transformation because the ilr transformed
+#' fractions can assumed to follow normal distribution. As the ilr
+#' transformation is not available in [RxODE] and can therefore not be used in
+#' the nlmixr modelling language, this transformation is currently used for
+#' translating mkin models with formation fractions to more than one target
+#' compartment for fitting with nlmixr in [nlmixr_model]. However,
+#' this implementation cannot be used there, as it is not accessible
+#' from RxODE.
+#'
+#' @param ff Vector of untransformed formation fractions. The sum
+#' must be smaller or equal to one
+#' @param ff_trans Vector of transformed formation fractions that can be
+#' restricted to the interval from 0 to 1
+#' @return A vector of the transformed formation fractions
+#' @export
+#' @examples
+#' ff_example <- c(
+#' 0.10983681, 0.09035905, 0.08399383
+#' )
+#' ff_example_trans <- tffm0(ff_example)
+#' invtffm0(ff_example_trans)
+tffm0 <- function(ff) {
+ n <- length(ff)
+ res <- numeric(n)
+ f_remaining <- 1
+ for (i in 1:n) {
+ res[i] <- ff[i]/f_remaining
+ f_remaining <- f_remaining - ff[i]
+ }
+ return(res)
+}
+#' @rdname tffm0
+#' @export
+#' @return A vector of backtransformed formation fractions for natural use in degradation models
+invtffm0 <- function(ff_trans) {
+ n <- length(ff_trans)
+ res <- numeric(n)
+ f_remaining <- 1
+ for (i in 1:n) {
+ res[i] <- ff_trans[i] * f_remaining
+ f_remaining <- f_remaining - res[i]
+ }
+ return(res)
+}
diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R
index 4fe4e5c2..174e7c2d 100644
--- a/R/transform_odeparms.R
+++ b/R/transform_odeparms.R
@@ -229,13 +229,18 @@ backtransform_odeparms <- function(transparms, mkinmod,
if (length(trans_f) > 0) {
if(transform_fractions) {
if (any(grepl("qlogis", names(trans_f)))) {
- parms[f_names] <- plogis(trans_f)
+ f_tmp <- plogis(trans_f)
+ if (any(grepl("_tffm0_.*_qlogis$", names(f_tmp)))) {
+ parms[f_names] <- invtffm0(f_tmp)
+ } else {
+ parms[f_names] <- f_tmp
+ }
} else {
- f <- invilr(trans_f)
+ f_tmp <- invilr(trans_f)
if (spec[[box]]$sink) {
- parms[f_names] <- f[1:length(f)-1]
+ parms[f_names] <- f_tmp[1:length(f_tmp)-1]
} else {
- parms[f_names] <- f
+ parms[f_names] <- f_tmp
}
}
} else {

Contact - Imprint