aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2021-06-16 18:03:22 +0200
committerJohannes Ranke <jranke@uni-bremen.de>2021-06-16 18:03:22 +0200
commit28197d5fcbaf85b39f4c032b8180d68b6f6a01b3 (patch)
treeba4cc748b6109c86e97292774d1ebb0ec3066803 /R
parent88cf130615a6cde0c4e65d14db32fed7f6e43085 (diff)
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.
Diffstat (limited to 'R')
-rw-r--r--R/dimethenamid_2018.R11
-rw-r--r--R/nlmixr.R79
-rw-r--r--R/tffm0.R46
3 files changed, 124 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)
+}

Contact - Imprint