Rudimentary support for setting up nlmixr models

- All degradation models are specified as ODE models. This appears to be
fast enough
- Error models are being translated to nlmixr as close to the mkin error
model as possible. When using the 'saem' backend, it appears not to be
possible to use the same error model for more than one observed variable
- No support yet for models with parallel formation of metabolites, where
the ilr transformation is used in mkin per default
- There is a bug in nlmixr which appears to be triggered if the data are
not balanced, see nlmixrdevelopment/nlmixr#530
- There is a print and a plot method, the summary method is not finished

1 files changed, 467 insertions, 0 deletions

diff --git a/R/nlmixr.R b/R/nlmixr.R
new file mode 100644
index 00000000..223b23a1
--- /dev/null
+++ b/R/nlmixr.R
@@ -0,0 +1,467 @@
+utils::globalVariables(c("predicted", "std"))
+
+#' 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
+#' @param object An [mmkin] row object containing several fits of the same
+#'   [mkinmod] model to different datasets
+#' @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 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 control Passed to [nlmixr::nlmixr].
+#' @param \dots Passed to [nlmixr_model]
+#' @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
+#' 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)
+#'
+#' \dontrun{
+#' 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,
+  ...,
+  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)
+  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)
+  })
+
+  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)
+
+  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
+#' @return An function defining a model suitable for fitting with [nlmixr::nlmixr].
+#' @export
+nlmixr_model <- function(object,
+  est = c("saem", "focei"),
+  degparms_start = "auto",
+  test_log_parms = FALSE, 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 <- degparms_mmkin$fixed
+
+  if (degparms_start[1] == "auto") {
+    degparms_start <- degparms_optim
+  }
+  degparms_fixed <- object[[1]]$bparms.fixed
+
+  degparms_optim_back_names <- names(backtransform_odeparms(degparms_optim,
+      object[[1]]$mkinmod,
+      object[[1]]$transform_rates,
+      object[[1]]$transform_fractions))
+  names(degparms_optim_back_names) <- names(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)
+
+  # Definition of the model function
+  f <- function(){}
+
+  ini_block <- "ini({"
+
+  # Initial values for all degradation parameters
+  for (parm_name in names(degparms_optim)) {
+    # 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(degparms_mmkin$eta[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_optim), 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_optim), value = TRUE)) {
+    model_block <- paste0(
+      model_block,
+      degparms_optim_back_names[parm_name], " = ",
+      "expit(", parm_name, " + eta.", parm_name, ")\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, "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)
+}