diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/add_err.R | 2 | ||||
| -rw-r--r-- | R/nlme.mmkin.R | 62 | ||||
| -rw-r--r-- | R/sigma_twocomp.R | 33 |
3 files changed, 86 insertions, 11 deletions
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) +} |