aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2020-10-22 12:34:40 +0200
committerJohannes Ranke <jranke@uni-bremen.de>2020-10-22 12:34:40 +0200
commit4a6beafe6ca119500232ecda4b5672dd4a1877c2 (patch)
treeade255f256a2cebf6262f12f816925ca3ce9944c /R
parenta9c7a1a8322567e9406a59ba0a4f910b89bd05e6 (diff)
Improve interface to experimental version of nlme
The experimental nlme version in my drat repository contains the variance function structure varConstProp which makes it possible to use the two-component error model in generalized nonlinear models using nlme::gnls() and in mixed effects models using nlme::nlme().
Diffstat (limited to 'R')
-rw-r--r--R/add_err.R2
-rw-r--r--R/nlme.mmkin.R62
-rw-r--r--R/sigma_twocomp.R33
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)
+}

Contact - Imprint