diff options
| -rw-r--r-- | .travis.yml | 5 | ||||
| -rw-r--r-- | DESCRIPTION | 2 | ||||
| -rw-r--r-- | R/add_err.R | 2 | ||||
| -rw-r--r-- | R/nlme.mmkin.R | 62 | ||||
| -rw-r--r-- | R/sigma_twocomp.R | 33 | ||||
| -rw-r--r-- | check.log | 4 | ||||
| -rw-r--r-- | man/add_err.Rd | 2 | ||||
| -rw-r--r-- | man/nlme.mmkin.Rd | 44 | ||||
| -rw-r--r-- | man/sigma_twocomp.Rd | 28 |
9 files changed, 158 insertions, 24 deletions
diff --git a/.travis.yml b/.travis.yml index b6982328..081a0f83 100644 --- a/.travis.yml +++ b/.travis.yml @@ -9,8 +9,9 @@ addons: apt: packages: - gcc -r_packages: - - nlme +repos: + CRAN: https://cloud.r-project.org + jranke: https://jranke.github.io/drat/ r_github_packages: - jranke/saemixextension script: diff --git a/DESCRIPTION b/DESCRIPTION index 80ba23d4..b56c77c7 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -17,7 +17,7 @@ Description: Calculation routines based on the FOCUS Kinetics Report (2006, note that no warranty is implied for correctness of results or fitness for a particular purpose. Imports: stats, graphics, methods, deSolve, R6, inline, parallel, numDeriv, - lmtest, pkgbuild, nlme (>= 3.1-149), purrr, saemix (>= 3.1.9000) + lmtest, pkgbuild, nlme (>= 3.1-150.1), purrr, saemix (>= 3.1.9000) Suggests: knitr, rbenchmark, tikzDevice, testthat, rmarkdown, covr, vdiffr, benchmarkme, tibble, stats4 License: GPL diff --git a/R/add_err.R b/R/add_err.R index d2092a84..8931a281 100644 --- a/R/add_err.R +++ b/R/add_err.R @@ -72,7 +72,7 @@ #' #' @export add_err <- function(prediction, sdfunc, secondary = c("M1", "M2"), - n = 1000, LOD = 0.1, reps = 2, digits = 1, seed = NA) + n = 10, LOD = 0.1, reps = 2, digits = 1, seed = NA) { if (!is.na(seed)) set.seed(seed) diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R index a94a26f7..7f7e34e9 100644 --- a/R/nlme.mmkin.R +++ b/R/nlme.mmkin.R @@ -47,16 +47,16 @@ get_deg_func <- function() { #' @examples #' ds <- lapply(experimental_data_for_UBA_2019[6:10], #' function(x) subset(x$data[c("name", "time", "value")], name == "parent")) -#' f <- mmkin("SFO", ds, quiet = TRUE, cores = 1) +#' f <- mmkin(c("SFO", "DFOP"), ds, quiet = TRUE, cores = 1) #' library(nlme) -#' endpoints(f[[1]]) -#' f_nlme <- nlme(f) -#' print(f_nlme) -#' endpoints(f_nlme) +#' f_nlme_sfo <- nlme(f["SFO", ]) +#' f_nlme_dfop <- nlme(f["DFOP", ]) +#' AIC(f_nlme_sfo, f_nlme_dfop) +#' print(f_nlme_dfop) +#' endpoints(f_nlme_dfop) #' \dontrun{ -#' f_nlme_2 <- nlme(f, start = c(parent_0 = 100, log_k_parent_sink = 0.1)) +#' f_nlme_2 <- nlme(f["SFO", ], start = c(parent_0 = 100, log_k_parent = 0.1)) #' update(f_nlme_2, random = parent_0 ~ 1) -#' # Test on some real data #' ds_2 <- lapply(experimental_data_for_UBA_2019[6:10], #' function(x) x$data[c("name", "time", "value")]) #' m_sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"), @@ -100,6 +100,36 @@ get_deg_func <- function() { #' #' endpoints(f_nlme_sfo_sfo) #' endpoints(f_nlme_dfop_sfo) +#' +#' if (findFunction("varConstProp")) { # tc error model for nlme available +#' # Attempts to fit metabolite kinetics with the tc error model +#' #f_2_tc <- mmkin(list("SFO-SFO" = m_sfo_sfo, +#' # "SFO-SFO-ff" = m_sfo_sfo_ff, +#' # "FOMC-SFO" = m_fomc_sfo, +#' # "DFOP-SFO" = m_dfop_sfo), +#' # ds_2, quiet = TRUE, +#' # error_model = "tc") +#' #f_nlme_sfo_sfo_tc <- nlme(f_2_tc["SFO-SFO", ], control = list(maxIter = 100)) +#' #f_nlme_dfop_sfo_tc <- nlme(f_2_tc["DFOP-SFO", ]) +#' #f_nlme_dfop_sfo_tc <- update(f_nlme_dfop_sfo, weights = varConstProp(), +#' # control = list(sigma = 1, msMaxIter = 100, pnlsMaxIter = 15)) +#' # Fitting metabolite kinetics with nlme.mmkin and the two-component +#' # error model currently does not work, at least not with these data. +#' +#' f_tc <- mmkin(c("SFO", "DFOP"), ds, quiet = TRUE, error_model = "tc") +#' f_nlme_sfo_tc <- nlme(f_tc["SFO", ]) +#' f_nlme_dfop_tc <- nlme(f_tc["DFOP", ]) +#' AIC(f_nlme_sfo, f_nlme_sfo_tc, f_nlme_dfop, f_nlme_dfop_tc) +#' print(f_nlme_dfop_tc) +#' } +#' f_2_obs <- mmkin(list("SFO-SFO" = m_sfo_sfo, +#' "DFOP-SFO" = m_dfop_sfo), +#' ds_2, quiet = TRUE, error_model = "obs") +#' f_nlme_sfo_sfo_obs <- nlme(f_2_obs["SFO-SFO", ]) +#' # The same with DFOP-SFO does not converge, apparently the variances of +#' # parent and A1 are too similar in this case, so that the model is +#' # overparameterised +#' #f_nlme_dfop_sfo_obs <- nlme(f_2_obs["DFOP-SFO", ], control = list(maxIter = 100)) #' } nlme.mmkin <- function(model, data = sys.frame(sys.parent()), fixed, random = fixed, @@ -145,6 +175,24 @@ nlme.mmkin <- function(model, data = sys.frame(sys.parent()), thisCall[["random"]] <- pdDiag(thisCall[["fixed"]]) } + error_model <- model[[1]]$err_mod + + if (missing(weights)) { + thisCall[["weights"]] <- switch(error_model, + const = NULL, + obs = varIdent(form = ~ 1 | name), + tc = varConstProp()) + sigma <- switch(error_model, + tc = 1, + NULL) + } + + control <- thisCall[["control"]] + if (error_model == "tc") { + control$sigma = 1 + thisCall[["control"]] <- control + } + val <- do.call("nlme.formula", thisCall) val$mmkin_orig <- model class(val) <- c("nlme.mmkin", "nlme", "lme") diff --git a/R/sigma_twocomp.R b/R/sigma_twocomp.R index 1e012d15..e8a92ced 100644 --- a/R/sigma_twocomp.R +++ b/R/sigma_twocomp.R @@ -23,7 +23,34 @@ #' #' Rocke, David M. and Lorenzato, Stefan (1995) A two-component model for #' measurement error in analytical chemistry. Technometrics 37(2), 176-184. +#' @examples +#' times <- c(0, 1, 3, 7, 14, 28, 60, 90, 120) +#' d_pred <- data.frame(time = times, parent = 100 * exp(- 0.03 * times)) +#' set.seed(123456) +#' d_syn <- add_err(d_pred, function(y) sigma_twocomp(y, 1, 0.07), +#' reps = 2, n = 1)[[1]] +#' f_nls <- nls(value ~ SSasymp(time, 0, parent_0, lrc), data = d_syn, +#' start = list(parent_0 = 100, lrc = -3)) +#' library(nlme) +#' f_gnls <- gnls(value ~ SSasymp(time, 0, parent_0, lrc), +#' data = d_syn, na.action = na.omit, +#' start = list(parent_0 = 100, lrc = -3)) +#' if (length(findFunction("varConstProp")) > 0) { +#' f_gnls_tc <- gnls(value ~ SSasymp(time, 0, parent_0, lrc), +#' data = d_syn, na.action = na.omit, +#' start = list(parent_0 = 100, lrc = -3), +#' weights = varConstProp()) +#' f_gnls_tc_sf <- gnls(value ~ SSasymp(time, 0, parent_0, lrc), +#' data = d_syn, na.action = na.omit, +#' start = list(parent_0 = 100, lrc = -3), +#' control = list(sigma = 1), +#' weights = varConstProp()) +#' } +#' f_mkin <- mkinfit("SFO", d_syn, error_model = "const", quiet = TRUE) +#' f_mkin_tc <- mkinfit("SFO", d_syn, error_model = "tc", quiet = TRUE) +#' plot_res(f_mkin_tc, standardized = TRUE) +#' AIC(f_nls, f_gnls, f_gnls_tc, f_gnls_tc_sf, f_mkin, f_mkin_tc) #' @export - sigma_twocomp <- function(y, sigma_low, rsd_high) { - sqrt(sigma_low^2 + y^2 * rsd_high^2) - } +sigma_twocomp <- function(y, sigma_low, rsd_high) { + sqrt(sigma_low^2 + y^2 * rsd_high^2) +} @@ -1,11 +1,11 @@ * using log directory ‘/home/jranke/git/mkin/mkin.Rcheck’ -* using R version 4.0.2 (2020-06-22) +* using R version 4.0.3 (2020-10-10) * using platform: x86_64-pc-linux-gnu (64-bit) * using session charset: UTF-8 * using options ‘--no-tests --as-cran’ * checking for file ‘mkin/DESCRIPTION’ ... OK * checking extension type ... Package -* this is package ‘mkin’ version ‘0.9.50.3’ +* this is package ‘mkin’ version ‘0.9.50.4’ * package encoding: UTF-8 * checking CRAN incoming feasibility ... Note_to_CRAN_maintainers Maintainer: ‘Johannes Ranke <jranke@uni-bremen.de>’ diff --git a/man/add_err.Rd b/man/add_err.Rd index 9527d508..fe7002ee 100644 --- a/man/add_err.Rd +++ b/man/add_err.Rd @@ -8,7 +8,7 @@ add_err( prediction, sdfunc, secondary = c("M1", "M2"), - n = 1000, + n = 10, LOD = 0.1, reps = 2, digits = 1, diff --git a/man/nlme.mmkin.Rd b/man/nlme.mmkin.Rd index 10c3ec78..0af670a0 100644 --- a/man/nlme.mmkin.Rd +++ b/man/nlme.mmkin.Rd @@ -77,16 +77,16 @@ have been obtained by fitting the same model to a list of datasets. \examples{ ds <- lapply(experimental_data_for_UBA_2019[6:10], function(x) subset(x$data[c("name", "time", "value")], name == "parent")) -f <- mmkin("SFO", ds, quiet = TRUE, cores = 1) +f <- mmkin(c("SFO", "DFOP"), ds, quiet = TRUE, cores = 1) library(nlme) -endpoints(f[[1]]) -f_nlme <- nlme(f) -print(f_nlme) -endpoints(f_nlme) +f_nlme_sfo <- nlme(f["SFO", ]) +f_nlme_dfop <- nlme(f["DFOP", ]) +AIC(f_nlme_sfo, f_nlme_dfop) +print(f_nlme_dfop) +endpoints(f_nlme_dfop) \dontrun{ - f_nlme_2 <- nlme(f, start = c(parent_0 = 100, log_k_parent_sink = 0.1)) + f_nlme_2 <- nlme(f["SFO", ], start = c(parent_0 = 100, log_k_parent = 0.1)) update(f_nlme_2, random = parent_0 ~ 1) - # Test on some real data ds_2 <- lapply(experimental_data_for_UBA_2019[6:10], function(x) x$data[c("name", "time", "value")]) m_sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"), @@ -130,6 +130,36 @@ endpoints(f_nlme) endpoints(f_nlme_sfo_sfo) endpoints(f_nlme_dfop_sfo) + + if (findFunction("varConstProp")) { # tc error model for nlme available + # Attempts to fit metabolite kinetics with the tc error model + #f_2_tc <- mmkin(list("SFO-SFO" = m_sfo_sfo, + # "SFO-SFO-ff" = m_sfo_sfo_ff, + # "FOMC-SFO" = m_fomc_sfo, + # "DFOP-SFO" = m_dfop_sfo), + # ds_2, quiet = TRUE, + # error_model = "tc") + #f_nlme_sfo_sfo_tc <- nlme(f_2_tc["SFO-SFO", ], control = list(maxIter = 100)) + #f_nlme_dfop_sfo_tc <- nlme(f_2_tc["DFOP-SFO", ]) + #f_nlme_dfop_sfo_tc <- update(f_nlme_dfop_sfo, weights = varConstProp(), + # control = list(sigma = 1, msMaxIter = 100, pnlsMaxIter = 15)) + # Fitting metabolite kinetics with nlme.mmkin and the two-component + # error model currently does not work, at least not with these data. + + f_tc <- mmkin(c("SFO", "DFOP"), ds, quiet = TRUE, error_model = "tc") + f_nlme_sfo_tc <- nlme(f_tc["SFO", ]) + f_nlme_dfop_tc <- nlme(f_tc["DFOP", ]) + AIC(f_nlme_sfo, f_nlme_sfo_tc, f_nlme_dfop, f_nlme_dfop_tc) + print(f_nlme_dfop_tc) + } + f_2_obs <- mmkin(list("SFO-SFO" = m_sfo_sfo, + "DFOP-SFO" = m_dfop_sfo), + ds_2, quiet = TRUE, error_model = "obs") + f_nlme_sfo_sfo_obs <- nlme(f_2_obs["SFO-SFO", ]) + # The same with DFOP-SFO does not converge, apparently the variances of + # parent and A1 are too similar in this case, so that the model is + # overparameterised + #f_nlme_dfop_sfo_obs <- nlme(f_2_obs["DFOP-SFO", ], control = list(maxIter = 100)) } } \seealso{ diff --git a/man/sigma_twocomp.Rd b/man/sigma_twocomp.Rd index 4e1f7c38..ed79d493 100644 --- a/man/sigma_twocomp.Rd +++ b/man/sigma_twocomp.Rd @@ -31,6 +31,34 @@ proposed by Rocke and Lorenzato (1995) can be written in this form as well, but assumes approximate lognormal distribution of errors for high values of y. } +\examples{ +times <- c(0, 1, 3, 7, 14, 28, 60, 90, 120) +d_pred <- data.frame(time = times, parent = 100 * exp(- 0.03 * times)) +set.seed(123456) +d_syn <- add_err(d_pred, function(y) sigma_twocomp(y, 1, 0.07), + reps = 2, n = 1)[[1]] +f_nls <- nls(value ~ SSasymp(time, 0, parent_0, lrc), data = d_syn, + start = list(parent_0 = 100, lrc = -3)) +library(nlme) +f_gnls <- gnls(value ~ SSasymp(time, 0, parent_0, lrc), + data = d_syn, na.action = na.omit, + start = list(parent_0 = 100, lrc = -3)) +if (length(findFunction("varConstProp")) > 0) { + f_gnls_tc <- gnls(value ~ SSasymp(time, 0, parent_0, lrc), + data = d_syn, na.action = na.omit, + start = list(parent_0 = 100, lrc = -3), + weights = varConstProp()) + f_gnls_tc_sf <- gnls(value ~ SSasymp(time, 0, parent_0, lrc), + data = d_syn, na.action = na.omit, + start = list(parent_0 = 100, lrc = -3), + control = list(sigma = 1), + weights = varConstProp()) +} +f_mkin <- mkinfit("SFO", d_syn, error_model = "const", quiet = TRUE) +f_mkin_tc <- mkinfit("SFO", d_syn, error_model = "tc", quiet = TRUE) +plot_res(f_mkin_tc, standardized = TRUE) +AIC(f_nls, f_gnls, f_gnls_tc, f_gnls_tc_sf, f_mkin, f_mkin_tc) +} \references{ Werner, Mario, Brooks, Samuel H., and Knott, Lancaster B. (1978) Additive, Multiplicative, and Mixed Analytical Errors. Clinical Chemistry |