aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
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