diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/dimethenamid_2018.R | 11 | ||||
| -rw-r--r-- | R/nlmixr.R | 79 | ||||
| -rw-r--r-- | R/tffm0.R | 46 |
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" @@ -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) +} |