From 28197d5fcbaf85b39f4c032b8180d68b6f6a01b3 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 16 Jun 2021 18:03:22 +0200 Subject: 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. --- R/dimethenamid_2018.R | 11 ++++- R/nlmixr.R | 79 ++++++++++++++++++++++++++++++----- R/tffm0.R | 46 ++++++++++++++++++++ tests/testthat/test_nlmixr.r | 99 ++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 223 insertions(+), 12 deletions(-) create mode 100644 R/tffm0.R create mode 100644 tests/testthat/test_nlmixr.r 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) + -- cgit v1.2.1