diff options
| -rw-r--r-- | NAMESPACE | 2 | ||||
| -rw-r--r-- | R/dimethenamid_2018.R | 14 | ||||
| -rw-r--r-- | R/summary.nlmixr.mmkin.R | 14 | ||||
| -rw-r--r-- | R/tffm0.R | 6 | ||||
| -rw-r--r-- | R/transform_odeparms.R | 13 | ||||
| -rw-r--r-- | check.log | 11 | ||||
| -rw-r--r-- | man/dimethenamid_2018.Rd | 36 | ||||
| -rw-r--r-- | man/nlmixr.mmkin.Rd | 15 | ||||
| -rw-r--r-- | man/tffm0.Rd | 42 | ||||
| -rw-r--r-- | test.log | 45 | ||||
| -rw-r--r-- | tests/testthat/test_nlmixr.r | 194 |
11 files changed, 243 insertions, 149 deletions
@@ -62,6 +62,7 @@ export(f_time_norm_focus) export(get_deg_func) export(ilr) export(invilr) +export(invtffm0) export(loftest) export(logistic.solution) export(lrtest) @@ -101,6 +102,7 @@ export(saem) export(saemix_data) export(saemix_model) export(sigma_twocomp) +export(tffm0) export(transform_odeparms) import(deSolve) import(graphics) diff --git a/R/dimethenamid_2018.R b/R/dimethenamid_2018.R index 76b98efe..6e0bda0c 100644 --- a/R/dimethenamid_2018.R +++ b/R/dimethenamid_2018.R @@ -31,6 +31,7 @@ #' dmta_ds[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]]) #' dmta_ds[["Elliot 1"]] <- NULL #' dmta_ds[["Elliot 2"]] <- NULL +#' \dontrun{ #' dfop_sfo3_plus <- mkinmod( #' DMTA = mkinsub("DFOP", c("M23", "M27", "M31")), #' M23 = mkinsub("SFO"), @@ -42,12 +43,15 @@ #' list("DFOP-SFO3+" = dfop_sfo3_plus), #' dmta_ds, quiet = TRUE, error_model = "tc") #' 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)) +#' control = nlmixr::foceiControl(print = 500)) #' summary(f_dmta_nlmixr_focei) #' plot(f_dmta_nlmixr_focei) +#' # saem has a problem with this model/data combination, maybe because of the +#' # overparameterised error model, to be investigated +#' #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) +#' } "dimethenamid_2018" diff --git a/R/summary.nlmixr.mmkin.R b/R/summary.nlmixr.mmkin.R index f2d7c607..a023f319 100644 --- a/R/summary.nlmixr.mmkin.R +++ b/R/summary.nlmixr.mmkin.R @@ -85,11 +85,11 @@ summary.nlmixr.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes mod_vars <- names(object$mkinmod$diffs) - pnames <- names(object$mean_dp_start) - np <- length(pnames) - conf.int <- confint(object$nm) - confint_trans <- as.matrix(conf.int[pnames, c(1, 3, 4)]) + dpnames <- setdiff(rownames(conf.int), names(object$mean_ep_start)) + ndp <- length(dpnames) + + confint_trans <- as.matrix(conf.int[dpnames, c(1, 3, 4)]) colnames(confint_trans) <- c("est.", "lower", "upper") bp <- backtransform_odeparms(confint_trans[, "est."], object$mkinmod, @@ -100,7 +100,7 @@ summary.nlmixr.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes # with the exception of sets of formation fractions (single fractions are OK). f_names_skip <- character(0) for (box in mod_vars) { # Figure out sets of fractions to skip - f_names <- grep(paste("^f", box, sep = "_"), pnames, value = TRUE) + f_names <- grep(paste("^f", box, sep = "_"), dpnames, value = TRUE) n_paths <- length(f_names) if (n_paths > 1) f_names_skip <- c(f_names_skip, f_names) } @@ -109,7 +109,7 @@ summary.nlmixr.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes dimnames = list(bpnames, colnames(confint_trans))) confint_back[, "est."] <- bp - for (pname in pnames) { + for (pname in dpnames) { if (!pname %in% f_names_skip) { par.lower <- confint_trans[pname, "lower"] par.upper <- confint_trans[pname, "upper"] @@ -131,7 +131,7 @@ summary.nlmixr.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes object$corFixed <- array( t(varFix/stdFix)/stdFix, dim(varFix), - list(pnames, pnames)) + list(dpnames, dpnames)) object$confint_trans <- confint_trans object$confint_back <- confint_back @@ -13,7 +13,8 @@ #' #' @param ff Vector of untransformed formation fractions. The sum #' must be smaller or equal to one -#' @param ff_trans +#' @param ff_trans Vector of transformed formation fractions that can be +#' restricted to the interval from 0 to 1 #' @return A vector of the transformed formation fractions #' @export #' @examples @@ -33,7 +34,8 @@ tffm0 <- function(ff) { return(res) } #' @rdname tffm0 -#' @return +#' @export +#' @return A vector of backtransformed formation fractions for natural use in degradation models invtffm0 <- function(ff_trans) { n <- length(ff_trans) res <- numeric(n) diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R index 4fe4e5c2..174e7c2d 100644 --- a/R/transform_odeparms.R +++ b/R/transform_odeparms.R @@ -229,13 +229,18 @@ backtransform_odeparms <- function(transparms, mkinmod, if (length(trans_f) > 0) { if(transform_fractions) { if (any(grepl("qlogis", names(trans_f)))) { - parms[f_names] <- plogis(trans_f) + f_tmp <- plogis(trans_f) + if (any(grepl("_tffm0_.*_qlogis$", names(f_tmp)))) { + parms[f_names] <- invtffm0(f_tmp) + } else { + parms[f_names] <- f_tmp + } } else { - f <- invilr(trans_f) + f_tmp <- invilr(trans_f) if (spec[[box]]$sink) { - parms[f_names] <- f[1:length(f)-1] + parms[f_names] <- f_tmp[1:length(f_tmp)-1] } else { - parms[f_names] <- f + parms[f_names] <- f_tmp } } } else { @@ -57,10 +57,7 @@ Maintainer: ‘Johannes Ranke <jranke@uni-bremen.de>’ * checking data for ASCII and uncompressed saves ... OK * checking installed files from ‘inst/doc’ ... OK * checking files in ‘vignettes’ ... OK -* checking examples ... NOTE -Examples with CPU (user + system) or elapsed time > 5s - user system elapsed -nlmixr.mmkin 8.129 0.375 5.384 +* checking examples ... OK * checking for unstated dependencies in ‘tests’ ... OK * checking tests ... SKIPPED * checking for unstated dependencies in vignettes ... OK @@ -71,9 +68,5 @@ nlmixr.mmkin 8.129 0.375 5.384 * checking for detritus in the temp directory ... OK * DONE -Status: 1 NOTE -See - ‘/home/jranke/git/mkin/mkin.Rcheck/00check.log’ -for details. - +Status: OK diff --git a/man/dimethenamid_2018.Rd b/man/dimethenamid_2018.Rd index b6f761e8..93fcad26 100644 --- a/man/dimethenamid_2018.Rd +++ b/man/dimethenamid_2018.Rd @@ -31,5 +31,41 @@ specific pieces of information in the comments. } \examples{ print(dimethenamid_2018) +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 +\dontrun{ +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") +nlmixr_model(f_dmta_mkin_tc) +f_dmta_nlmixr_focei <- nlmixr(f_dmta_mkin_tc, est = "focei", + control = nlmixr::foceiControl(print = 500)) +summary(f_dmta_nlmixr_focei) +plot(f_dmta_nlmixr_focei) +# saem has a problem with this model/data combination, maybe because of the +# overparameterised error model, to be investigated +#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) +} } \keyword{datasets} diff --git a/man/nlmixr.mmkin.Rd b/man/nlmixr.mmkin.Rd index 4ab30272..0f4f41a2 100644 --- a/man/nlmixr.mmkin.Rd +++ b/man/nlmixr.mmkin.Rd @@ -16,6 +16,8 @@ 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() @@ -27,7 +29,8 @@ nlmixr_model( object, est = c("saem", "focei"), degparms_start = "auto", - test_log_parms = FALSE, + eta_start = "auto", + test_log_parms = TRUE, conf.level = 0.6, error_model = object[[1]]$err_mod, add_attributes = FALSE @@ -58,6 +61,13 @@ when calculating mean degradation parameters using \link{mean_degparms}.} \item{conf.level}{Possibility to adjust the required confidence level for parameter that are tested if requested by 'test_log_parms'.} +\item{degparms_start}{Parameter values given as a named numeric vector will +be used to override the starting values obtained from the 'mmkin' object.} + +\item{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.} + \item{\dots}{Passed to \link{nlmixr_model}} \item{save}{Passed to \link[nlmixr:nlmixr]{nlmixr::nlmixr}} @@ -68,9 +78,6 @@ for parameter that are tested if requested by 'test_log_parms'.} \item{digits}{Number of digits to use for printing} -\item{degparms_start}{Parameter values given as a named numeric vector will -be used to override the starting values obtained from the 'mmkin' object.} - \item{add_attributes}{Should the starting values used for degradation model parameters and their distribution and for the error model parameters be returned as attributes?} diff --git a/man/tffm0.Rd b/man/tffm0.Rd new file mode 100644 index 00000000..46978d5e --- /dev/null +++ b/man/tffm0.Rd @@ -0,0 +1,42 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/tffm0.R +\name{tffm0} +\alias{tffm0} +\alias{invtffm0} +\title{Transform formation fractions as in the first published mkin version} +\usage{ +tffm0(ff) + +invtffm0(ff_trans) +} +\arguments{ +\item{ff}{Vector of untransformed formation fractions. The sum +must be smaller or equal to one} + +\item{ff_trans}{Vector of transformed formation fractions that can be +restricted to the interval from 0 to 1} +} +\value{ +A vector of the transformed formation fractions + +A vector of backtransformed formation fractions for natural use in degradation models +} +\description{ +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 \link{ilr} transformation because the ilr transformed +fractions can assumed to follow normal distribution. As the ilr +transformation is not available in \link{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 \link{nlmixr_model}. However, +this implementation cannot be used there, as it is not accessible +from RxODE. +} +\examples{ +ff_example <- c( + 0.10983681, 0.09035905, 0.08399383 +) +ff_example_trans <- tffm0(ff_example) +invtffm0(ff_example_trans) +} @@ -3,17 +3,17 @@ Loading required package: parallel ℹ Testing mkin ✔ | OK F W S | Context ✔ | 5 | AIC calculation -✔ | 5 | Analytical solutions for coupled models [3.3 s] +✔ | 5 | Analytical solutions for coupled models [3.5 s] ✔ | 5 | Calculation of Akaike weights ✔ | 2 | Export dataset for reading into CAKE -✔ | 12 | Confidence intervals and p-values [1.3 s] +✔ | 12 | Confidence intervals and p-values [1.0 s] ✔ | 14 | Error model fitting [4.7 s] ✔ | 5 | Time step normalisation -✔ | 4 | Calculation of FOCUS chi2 error levels [0.5 s] +✔ | 4 | Calculation of FOCUS chi2 error levels [0.6 s] ✔ | 14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.8 s] ✔ | 4 | Test fitting the decline of metabolites from their maximum [0.3 s] ✔ | 1 | Fitting the logistic model [0.2 s] -✔ | 35 1 | Nonlinear mixed-effects models [27.1 s] +✔ | 35 1 | Nonlinear mixed-effects models [26.9 s] ──────────────────────────────────────────────────────────────────────────────── Skip (test_mixed.R:161:3): saem results are reproducible for biphasic fits Reason: Fitting with saemix takes around 10 minutes when using deSolve @@ -21,33 +21,36 @@ Reason: Fitting with saemix takes around 10 minutes when using deSolve ✔ | 2 | Test dataset classes mkinds and mkindsg ✔ | 10 | Special cases of mkinfit calls [0.4 s] ✔ | 1 | mkinfit features [0.3 s] -✔ | 8 | mkinmod model generation and printing [0.3 s] +✔ | 8 | mkinmod model generation and printing [0.2 s] ✔ | 3 | Model predictions with mkinpredict [0.3 s] -✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.7 s] -✔ | 9 | Nonlinear mixed-effects models with nlme [8.1 s] -✖ | 14 2 | Plotting [1.9 s] +✔ | 14 2 | Evaluations according to 2015 NAFTA guidance [1.3 s] ──────────────────────────────────────────────────────────────────────────────── -Failure (test_plot.R:40:5): Plotting mkinfit, mmkin and mixed model objects is reproducible -Figures don't match: mixed-model-fit-for-saem-object-with-saemix-transformations.svg - - -Failure (test_plot.R:55:5): Plotting mkinfit, mmkin and mixed model objects is reproducible -Figures don't match: mixed-model-fit-for-saem-object-with-mkin-transformations.svg +Skip (test_nafta.R:25:5): Test data from Appendix B are correctly evaluated +Reason: getRversion() >= "4.1.0" is TRUE +Skip (test_nafta.R:53:5): Test data from Appendix D are correctly evaluated +Reason: getRversion() >= "4.1.0" is TRUE +──────────────────────────────────────────────────────────────────────────────── +✔ | 9 | Nonlinear mixed-effects models with nlme [8.1 s] +✔ | 0 1 | Plotting [0.8 s] +──────────────────────────────────────────────────────────────────────────────── +Skip (test_plot.R:18:3): Plotting mkinfit, mmkin and mixed model objects is reproducible +Reason: getRversion() >= "4.1.0" is TRUE ──────────────────────────────────────────────────────────────────────────────── ✔ | 4 | Residuals extracted from mkinfit models ✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.5 s] -✔ | 7 | Fitting the SFORB model [3.9 s] +✔ | 7 | Fitting the SFORB model [3.8 s] ✔ | 1 | Summaries of old mkinfit objects ✔ | 4 | Summary [0.1 s] -✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2 s] -✔ | 9 | Hypothesis tests [8.2 s] -✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.4 s] +✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.3 s] +✔ | 9 | Hypothesis tests [8.5 s] +✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2 s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 70.0 s +Duration: 68.2 s ── Skipped tests ────────────────────────────────────────────────────────────── -● Fitting with saemix takes around 10 minutes when using deSolve (1) +• Fitting with saemix takes around 10 minutes when using deSolve (1) +• getRversion() >= "4.1.0" is TRUE (3) -[ FAIL 2 | WARN 0 | SKIP 1 | PASS 204 ] +[ FAIL 0 | WARN 0 | SKIP 4 | PASS 188 ] diff --git a/tests/testthat/test_nlmixr.r b/tests/testthat/test_nlmixr.r index e3bd3d66..dcbb50ac 100644 --- a/tests/testthat/test_nlmixr.r +++ b/tests/testthat/test_nlmixr.r @@ -1,99 +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) - +# 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) +# |