Translate formation fractions to nlmixr language

Works for the dimethenamid data, at least for FOCEI. Very little testing
yet. The summary function does not yet handle the new transformations
of formation fractions (that are in fact very old, as they were used
in the very first version of mkin). The test file has no tests yet, just
some code that may be used for testing.

4 files changed, 223 insertions, 12 deletions

diff --git a/R/dimethenamid_2018.R b/R/dimethenamid_2018.R
index 79018c11..76b98efe 100644
--- a/R/dimethenamid_2018.R
+++ b/R/dimethenamid_2018.R
@@ -41,6 +41,13 @@
 #' f_dmta_mkin_tc <- mmkin(
 #'   list("DFOP-SFO3+" = dfop_sfo3_plus),
 #'   dmta_ds, quiet = TRUE, error_model = "tc")
-#' nlmixr_model(f_dmta_mkin_tc) # incomplete
-#' # nlmixr(f_dmta_mkin_tc, est = "saem") # not supported (yet)
+#' nlmixr_model(f_dmta_mkin_tc)
+#' f_dmta_nlmixr_saem <- nlmixr(f_dmta_mkin_tc, est = "saem",
+#'   control = saemControl(print = 500))
+#' summary(f_dmta_nlmixr_saem)
+#' plot(f_dmta_nlmixr_saem)
+#' f_dmta_nlmixr_focei <- nlmixr(f_dmta_mkin_tc, est = "focei",
+#'   control = foceiControl(print = 500))
+#' summary(f_dmta_nlmixr_focei)
+#' plot(f_dmta_nlmixr_focei)
 "dimethenamid_2018"
diff --git a/R/nlmixr.R b/R/nlmixr.R
index 6e0b5128..9c364c4f 100644
--- a/R/nlmixr.R
+++ b/R/nlmixr.R
@@ -20,6 +20,9 @@ nlmixr::nlmixr
 #' @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
@@ -148,6 +151,8 @@ nlmixr.mmkin <- function(object, data = NULL,
   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()
@@ -155,7 +160,9 @@ nlmixr.mmkin <- function(object, data = NULL,
 {
   m_nlmixr <- nlmixr_model(object, est = est,
     error_model = error_model, add_attributes = TRUE,
-    test_log_parms = test_log_parms, conf.level = conf.level)
+    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")
@@ -164,7 +171,7 @@ nlmixr.mmkin <- function(object, data = NULL,
   attributes(m_nlmixr) <- NULL
 
   fit_time <- system.time({
-    f_nlmixr <- nlmixr(m_nlmixr, d_nlmixr, est = est)
+    f_nlmixr <- nlmixr(m_nlmixr, d_nlmixr, est = est, control = control)
   })
 
   if (is.null(f_nlmixr$CMT)) {
@@ -246,7 +253,8 @@ print.nlmixr.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...)
 nlmixr_model <- function(object,
   est = c("saem", "focei"),
   degparms_start = "auto",
-  test_log_parms = FALSE, conf.level = 0.6,
+  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")
@@ -278,16 +286,44 @@ nlmixr_model <- function(object,
     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)
-  degparms_optim_back_names <- names(degparms_optim_back)
-  names(degparms_optim_back_names) <- names(degparms_optim)
 
   if (degparms_start[1] == "auto") {
-    degparms_start <- degparms_optim
+    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)
@@ -315,7 +351,7 @@ nlmixr_model <- function(object,
       as.character(degparms_start[parm_name]),
       "\n",
       "eta.", parm_name, " ~ ",
-      as.character(degparms_mmkin$eta[parm_name]),
+      as.character(eta_start[parm_name]),
       "\n"
     )
   }
@@ -394,7 +430,7 @@ nlmixr_model <- function(object,
   }
 
   # Population initial values for log rate constants
-  for (parm_name in grep("^log_", names(degparms_optim), value = TRUE)) {
+  for (parm_name in grep("^log_", names(degparms_start), value = TRUE)) {
     model_block <- paste0(
       model_block,
       gsub("^log_", "", parm_name), " = ",
@@ -402,13 +438,36 @@ nlmixr_model <- function(object,
   }
 
   # Population initial values for logit transformed parameters
-  for (parm_name in grep("_qlogis$", names(degparms_optim), value = TRUE)) {
+  for (parm_name in grep("_qlogis$", names(degparms_start), value = TRUE)) {
     model_block <- paste0(
       model_block,
-      degparms_optim_back_names[parm_name], " = ",
+      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,
diff --git a/R/tffm0.R b/R/tffm0.R
new file mode 100644
index 00000000..25787962
--- /dev/null
+++ b/R/tffm0.R
@@ -0,0 +1,46 @@
+#' 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
+#' @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
+#' @return
+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/tests/testthat/test_nlmixr.r b/tests/testthat/test_nlmixr.r
new file mode 100644
index 00000000..e3bd3d66
--- /dev/null
+++ b/tests/testthat/test_nlmixr.r
@@ -0,0 +1,99 @@
+
+
+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
+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")
+
+d_dmta_nlmixr <- nlmixr_data(f_dmta_mkin_tc)
+m_dmta_nlmixr <- function ()
+{
+    ini({
+        DMTA_0 = 98.7697627680706
+        eta.DMTA_0 ~ 2.35171765917765
+        log_k_M23 = -3.92162409637283
+        eta.log_k_M23 ~ 0.549278519419884
+        log_k_M27 = -4.33774620773911
+        eta.log_k_M27 ~ 0.864474956685295
+        log_k_M31 = -4.24767627688461
+        eta.log_k_M31 ~ 0.750297149164171
+        f_DMTA_tffm0_1_qlogis = -2.092409
+        eta.f_DMTA_tffm0_1_qlogis ~ 0.3
+        f_DMTA_tffm0_2_qlogis = -2.180576
+        eta.f_DMTA_tffm0_2_qlogis ~ 0.3
+        f_DMTA_tffm0_3_qlogis = -2.142672
+        eta.f_DMTA_tffm0_3_qlogis ~ 0.3
+        log_k1 = -2.2341008812259
+        eta.log_k1 ~ 0.902976221565793
+        log_k2 = -3.7762779983269
+        eta.log_k2 ~ 1.57684519529298
+        g_qlogis = 0.450175725479389
+        eta.g_qlogis ~ 3.0851335687675
+        sigma_low_DMTA = 0.697933852349996
+        rsd_high_DMTA = 0.0257724286053519
+        sigma_low_M23 = 0.697933852349996
+        rsd_high_M23 = 0.0257724286053519
+        sigma_low_M27 = 0.697933852349996
+        rsd_high_M27 = 0.0257724286053519
+        sigma_low_M31 = 0.697933852349996
+        rsd_high_M31 = 0.0257724286053519
+    })
+    model({
+        DMTA_0_model = DMTA_0 + eta.DMTA_0
+        DMTA(0) = DMTA_0_model
+        k_M23 = exp(log_k_M23 + eta.log_k_M23)
+        k_M27 = exp(log_k_M27 + eta.log_k_M27)
+        k_M31 = exp(log_k_M31 + eta.log_k_M31)
+        k1 = exp(log_k1 + eta.log_k1)
+        k2 = exp(log_k2 + eta.log_k2)
+        g = expit(g_qlogis + eta.g_qlogis)
+        f_DMTA_tffm0_1 = expit(f_DMTA_tffm0_1_qlogis + eta.f_DMTA_tffm0_1_qlogis)
+        f_DMTA_tffm0_2 = expit(f_DMTA_tffm0_2_qlogis + eta.f_DMTA_tffm0_2_qlogis)
+        f_DMTA_tffm0_3 = expit(f_DMTA_tffm0_3_qlogis + eta.f_DMTA_tffm0_3_qlogis)
+        f_DMTA_to_M23 = f_DMTA_tffm0_1
+        f_DMTA_to_M27 = (1 - f_DMTA_tffm0_1) * f_DMTA_tffm0_2
+        f_DMTA_to_M31 = (1 - f_DMTA_tffm0_1) * (1 - f_DMTA_tffm0_2) * f_DMTA_tffm0_3
+        d/dt(DMTA) = -((k1 * g * exp(-k1 * time) + k2 * (1 -
+            g) * exp(-k2 * time))/(g * exp(-k1 * time) + (1 -
+            g) * exp(-k2 * time))) * DMTA
+        d/dt(M23) = +f_DMTA_to_M23 * ((k1 * g * exp(-k1 * time) +
+            k2 * (1 - g) * exp(-k2 * time))/(g * exp(-k1 * time) +
+            (1 - g) * exp(-k2 * time))) * DMTA - k_M23 * M23
+        d/dt(M27) = +f_DMTA_to_M27 * ((k1 * g * exp(-k1 * time) +
+            k2 * (1 - g) * exp(-k2 * time))/(g * exp(-k1 * time) +
+            (1 - g) * exp(-k2 * time))) * DMTA - k_M27 * M27 +
+            k_M31 * M31
+        d/dt(M31) = +f_DMTA_to_M31 * ((k1 * g * exp(-k1 * time) +
+            k2 * (1 - g) * exp(-k2 * time))/(g * exp(-k1 * time) +
+            (1 - g) * exp(-k2 * time))) * DMTA - k_M31 * M31
+        DMTA ~ add(sigma_low_DMTA) + prop(rsd_high_DMTA)
+        M23 ~ add(sigma_low_M23) + prop(rsd_high_M23)
+        M27 ~ add(sigma_low_M27) + prop(rsd_high_M27)
+        M31 ~ add(sigma_low_M31) + prop(rsd_high_M31)
+    })
+}
+m_dmta_nlmixr_mkin <- nlmixr_model(f_dmta_mkin_tc, test_log_parms = TRUE)
+f_dmta_nlmixr_saem <- nlmixr(f_dmta_mkin_tc, est = "saem", control = saemControl(print = 250))
+f_dmta_nlmixr_focei <- nlmixr(f_dmta_mkin_tc, est = "focei", control = foceiControl(print = 250))
+plot(f_dmta_nlmixr_saem)
+plot(f_dmta_nlmixr_focei)
+