aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2020-04-14 17:22:00 +0200
committerJohannes Ranke <jranke@uni-bremen.de>2020-04-14 17:22:00 +0200
commit26085403289e29259e500282e8e88a5ab00c07a0 (patch)
treee86d5195c725ee0c08a45e4393ff9c1dfad64314
parentd4a49b4837de347d34b2c198de7342c34b0fab63 (diff)
Add a nlme method for mmkin row objects
-rw-r--r--NAMESPACE2
-rw-r--r--R/nlme.R17
-rw-r--r--R/nlme.mmkin.R85
-rw-r--r--check.log102
-rw-r--r--man/nlme.Rd15
-rw-r--r--man/nlme.mmkin.Rd72
6 files changed, 270 insertions, 23 deletions
diff --git a/NAMESPACE b/NAMESPACE
index 1b53f101..5f305a3d 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -12,6 +12,7 @@ S3method(lrtest,mkinfit)
S3method(lrtest,mmkin)
S3method(mkinpredict,mkinfit)
S3method(mkinpredict,mkinmod)
+S3method(nlme,mmkin)
S3method(nobs,mkinfit)
S3method(parms,mkinfit)
S3method(plot,mkinfit)
@@ -91,6 +92,7 @@ importFrom(stats,dist)
importFrom(stats,dnorm)
importFrom(stats,lm)
importFrom(stats,logLik)
+importFrom(stats,na.fail)
importFrom(stats,nlminb)
importFrom(stats,nobs)
importFrom(stats,optimize)
diff --git a/R/nlme.R b/R/nlme.R
index 12a3104c..fafaa7b6 100644
--- a/R/nlme.R
+++ b/R/nlme.R
@@ -1,4 +1,4 @@
-#' Estimation of parameter distributions from mmkin row objects
+#' Helper functions to create nlme models from mmkin row objects
#'
#' These functions facilitate setting up a nonlinear mixed effects model for
#' an mmkin row object. An mmkin row object is essentially a list of mkinfit
@@ -81,14 +81,19 @@
#' nlme_f_sfo_sfo <- nlme_function(f_2["SFO-SFO", ])
#' nlme_f_sfo_sfo_ff <- nlme_function(f_2["SFO-SFO-ff", ])
#' nlme_f_fomc_sfo <- nlme_function(f_2["FOMC-SFO", ])
+#' assign("nlme_f_sfo_sfo", nlme_f_sfo_sfo, globalenv())
+#' assign("nlme_f_sfo_sfo_ff", nlme_f_sfo_sfo_ff, globalenv())
+#' assign("nlme_f_fomc_sfo", nlme_f_fomc_sfo, globalenv())
#'
-#' # Allowing for correlations between random effects leads to non-convergence
+#' # Allowing for correlations between random effects (not shown)
+#' # leads to non-convergence
#' f_nlme_sfo_sfo <- nlme(value ~ nlme_f_sfo_sfo(name, time,
#' parent_0, log_k_parent_sink, log_k_parent_A1, log_k_A1_sink),
#' data = grouped_data_2,
#' fixed = parent_0 + log_k_parent_sink + log_k_parent_A1 + log_k_A1_sink ~ 1,
#' random = pdDiag(parent_0 + log_k_parent_sink + log_k_parent_A1 + log_k_A1_sink ~ 1),
#' start = mean_dp_sfo_sfo)
+#' # augPred does not see to work on this object, so no plot is shown
#'
#' # The same model fitted with transformed formation fractions does not converge
#' f_nlme_sfo_sfo_ff <- nlme(value ~ nlme_f_sfo_sfo_ff(name, time,
@@ -98,14 +103,6 @@
#' random = pdDiag(parent_0 + log_k_parent + log_k_A1 + f_parent_ilr_1 ~ 1),
#' start = mean_dp_sfo_sfo_ff)
#'
-#' # It does converge with this version of reduced random effects
-#' f_nlme_sfo_sfo_ff <- nlme(value ~ nlme_f_sfo_sfo_ff(name, time,
-#' parent_0, log_k_parent, log_k_A1, f_parent_ilr_1),
-#' data = grouped_data_2,
-#' fixed = parent_0 + log_k_parent + log_k_A1 + f_parent_ilr_1 ~ 1,
-#' random = pdDiag(parent_0 + log_k_parent ~ 1),
-#' start = mean_dp_sfo_sfo_ff)
-#'
#' f_nlme_fomc_sfo <- nlme(value ~ nlme_f_fomc_sfo(name, time,
#' parent_0, log_alpha, log_beta, log_k_A1, f_parent_ilr_1),
#' data = grouped_data_2,
diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R
new file mode 100644
index 00000000..784ba609
--- /dev/null
+++ b/R/nlme.mmkin.R
@@ -0,0 +1,85 @@
+#' Create an nlme model for an mmkin row object
+#'
+#' This functions sets up a nonlinear mixed effects model for an mmkin row
+#' object. An mmkin row object is essentially a list of mkinfit objects that
+#' have been obtained by fitting the same model to a list of datasets.
+#'
+#' @param model An \code{\link{mmkin}} row object.
+#' @param data Ignored, data are taken from the mmkin model
+#' @param fixed Ignored, all degradation parameters fitted in the
+#' mmkin model are used as fixed parameters
+#' @param random If not specified, all fixed effects are complemented
+#' with uncorrelated random effects
+#' @param groups See the documentation of nlme
+#' @param start If not specified, mean values of the fitted degradation
+#' parameters taken from the mmkin object are used
+#' @param correlation See the documentation of nlme
+#' @param weights passed to nlme
+#' @param subset passed to nlme
+#' @param method passed to nlme
+#' @param na.action passed to nlme
+#' @param naPattern passed to nlme
+#' @param control passed to nlme
+#' @param verbose passed to nlme
+#' @importFrom stats na.fail
+#' @return Upon success, a fitted nlme.mmkin object, which is
+#' an nlme object with additional elements
+#' @export
+#' @seealso \code{\link{nlme_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)
+#' library(nlme)
+#' f_nlme <- nlme(f)
+#' nlme(f, random = parent_0 ~ 1)
+#' f_nlme <- nlme(f, start = c(parent_0 = 100, log_k_parent_sink = 0.1))
+#' update(f_nlme, random = parent_0 ~ 1)
+# Code inspired by nlme.nlsList
+nlme.mmkin <- function(model, data = sys.frame(sys.parent()),
+ fixed, random = fixed,
+ groups, start, correlation = NULL, weights = NULL,
+ subset, method = c("ML", "REML"),
+ na.action = na.fail, naPattern,
+ control = list(), verbose= FALSE)
+{
+ if (nrow(model) > 1) stop("Only row objects allowed")
+
+ thisCall <- as.list(match.call())[-1]
+
+ # warn in case of use of arguments that are overriden
+ if (any(!is.na(match(names(thisCall),
+ c("fixed", "data"))))) {
+ warning("'nlme.mmkin' will redefine 'fixed' and 'data'")
+ }
+
+ deg_func <- nlme_function(model)
+ assign("deg_func", deg_func, parent.frame())
+
+ # specify the model formula
+ this_model_text <- paste0("value ~ deg_func(",
+ paste(names(formals(deg_func)), collapse = ", "), ")")
+ this_model <- eval(parse(text = this_model_text))
+ thisCall[["model"]] <- this_model
+
+ mean_dp <- mean_degparms(model)
+ dp_names <- names(mean_dp)
+
+ thisCall[["data"]] <- nlme_data(model)
+
+ if (missing(start)) {
+ thisCall[["start"]] <- mean_dp
+ }
+
+ thisCall[["fixed"]] <- lapply(as.list(dp_names), function(el)
+ eval(parse(text = paste(el, 1, sep = "~"))))
+
+ if (missing(random)) {
+ thisCall[["random"]] <- pdDiag(thisCall[["fixed"]])
+ }
+
+ val <- do.call("nlme.formula", thisCall)
+ return(val)
+ class(val) <- c("nlme.mmkin", "nlme", "lme")
+}
+
diff --git a/check.log b/check.log
index 0a75a713..fc06bc0d 100644
--- a/check.log
+++ b/check.log
@@ -5,7 +5,7 @@
* using options ‘--no-tests --as-cran’
* checking for file ‘mkin/DESCRIPTION’ ... OK
* checking extension type ... Package
-* this is package ‘mkin’ version ‘0.9.49.9’
+* this is package ‘mkin’ version ‘0.9.49.10’
* package encoding: UTF-8
* checking CRAN incoming feasibility ... Note_to_CRAN_maintainers
Maintainer: ‘Johannes Ranke <jranke@uni-bremen.de>’
@@ -41,7 +41,14 @@ Maintainer: ‘Johannes Ranke <jranke@uni-bremen.de>’
* checking S3 generic/method consistency ... OK
* checking replacement functions ... OK
* checking foreign function calls ... OK
-* checking R code for possible problems ... OK
+* checking R code for possible problems ... NOTE
+nlme.mmkin: no visible binding for global variable ‘na.fail’
+nlme.mmkin: no visible binding for global variable ‘f’
+Undefined global functions or variables:
+ f na.fail
+Consider adding
+ importFrom("stats", "na.fail")
+to your NAMESPACE file.
* checking Rd files ... OK
* checking Rd metadata ... OK
* checking Rd line widths ... OK
@@ -56,7 +63,91 @@ 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 ... OK
+* checking examples ... ERROR
+Running examples in ‘mkin-Ex.R’ failed
+The error most likely occurred in:
+
+> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
+> ### Name: nlme_function
+> ### Title: Estimation of parameter distributions from mmkin row objects
+> ### Aliases: nlme_function mean_degparms nlme_data
+>
+> ### ** Examples
+>
+> sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120)
+> m_SFO <- mkinmod(parent = mkinsub("SFO"))
+> d_SFO_1 <- mkinpredict(m_SFO,
++ c(k_parent_sink = 0.1),
++ c(parent = 98), sampling_times)
+> d_SFO_1_long <- mkin_wide_to_long(d_SFO_1, time = "time")
+> d_SFO_2 <- mkinpredict(m_SFO,
++ c(k_parent_sink = 0.05),
++ c(parent = 102), sampling_times)
+> d_SFO_2_long <- mkin_wide_to_long(d_SFO_2, time = "time")
+> d_SFO_3 <- mkinpredict(m_SFO,
++ c(k_parent_sink = 0.02),
++ c(parent = 103), sampling_times)
+> d_SFO_3_long <- mkin_wide_to_long(d_SFO_3, time = "time")
+>
+> d1 <- add_err(d_SFO_1, function(value) 3, n = 1)
+> d2 <- add_err(d_SFO_2, function(value) 2, n = 1)
+> d3 <- add_err(d_SFO_3, function(value) 4, n = 1)
+> ds <- c(d1 = d1, d2 = d2, d3 = d3)
+>
+> f <- mmkin("SFO", ds, cores = 1, quiet = TRUE)
+> mean_dp <- mean_degparms(f)
+> grouped_data <- nlme_data(f)
+> nlme_f <- nlme_function(f)
+> # These assignments are necessary for these objects to be
+> # visible to nlme and augPred when evaluation is done by
+> # pkgdown to generated the html docs.
+> assign("nlme_f", nlme_f, globalenv())
+> assign("grouped_data", grouped_data, globalenv())
+>
+> library(nlme)
+> m_nlme <- nlme(value ~ nlme_f(name, time, parent_0, log_k_parent_sink),
++ data = grouped_data,
++ fixed = parent_0 + log_k_parent_sink ~ 1,
++ random = pdDiag(parent_0 + log_k_parent_sink ~ 1),
++ start = mean_dp)
+> summary(m_nlme)
+Nonlinear mixed-effects model fit by maximum likelihood
+ Model: value ~ nlme_f(name, time, parent_0, log_k_parent_sink)
+ Data: grouped_data
+ AIC BIC logLik
+ 253.6377 262.9937 -121.8188
+
+Random effects:
+ Formula: list(parent_0 ~ 1, log_k_parent_sink ~ 1)
+ Level: ds
+ Structure: Diagonal
+ parent_0 log_k_parent_sink Residual
+StdDev: 2.486822 0.6481355 2.368137
+
+Fixed effects: parent_0 + log_k_parent_sink ~ 1
+ Value Std.Error DF t-value p-value
+parent_0 101.43096 1.5960157 44 63.55261 0
+log_k_parent_sink -3.07614 0.3827294 44 -8.03737 0
+ Correlation:
+ prnt_0
+log_k_parent_sink 0.011
+
+Standardized Within-Group Residuals:
+ Min Q1 Med Q3 Max
+-1.9643731 -0.5430765 0.1659516 0.6160408 1.8082871
+
+Number of Observations: 48
+Number of Groups: 3
+> plot(augPred(m_nlme, level = 0:1), layout = c(3, 1))
+>
+> m_nlme <- nlme(nlme_formula,
++ data = grouped_data,
++ fixed = parent_0 + log_k_parent_sink ~ 1,
++ random = pdDiag(parent_0 + log_k_parent_sink ~ 1),
++ start = mean_dp)
+Error in nlme(nlme_formula, data = grouped_data, fixed = parent_0 + log_k_parent_sink ~ :
+ object 'nlme_formula' not found
+Execution halted
* checking for unstated dependencies in ‘tests’ ... OK
* checking tests ... SKIPPED
* checking for unstated dependencies in vignettes ... OK
@@ -66,5 +157,8 @@ Maintainer: ‘Johannes Ranke <jranke@uni-bremen.de>’
* checking for detritus in the temp directory ... OK
* DONE
-Status: OK
+Status: 1 ERROR, 1 NOTE
+See
+ ‘/home/jranke/git/mkin/mkin.Rcheck/00check.log’
+for details.
diff --git a/man/nlme.Rd b/man/nlme.Rd
index 8e5c2aa0..971ba3f5 100644
--- a/man/nlme.Rd
+++ b/man/nlme.Rd
@@ -101,14 +101,19 @@ plot(augPred(m_nlme, level = 0:1), layout = c(3, 1))
nlme_f_sfo_sfo <- nlme_function(f_2["SFO-SFO", ])
nlme_f_sfo_sfo_ff <- nlme_function(f_2["SFO-SFO-ff", ])
nlme_f_fomc_sfo <- nlme_function(f_2["FOMC-SFO", ])
+ assign("nlme_f_sfo_sfo", nlme_f_sfo_sfo, globalenv())
+ assign("nlme_f_sfo_sfo_ff", nlme_f_sfo_sfo_ff, globalenv())
+ assign("nlme_f_fomc_sfo", nlme_f_fomc_sfo, globalenv())
- # Allowing for correlations between random effects leads to non-convergence
+ # Allowing for correlations between random effects (not shown)
+ # leads to non-convergence
f_nlme_sfo_sfo <- nlme(value ~ nlme_f_sfo_sfo(name, time,
parent_0, log_k_parent_sink, log_k_parent_A1, log_k_A1_sink),
data = grouped_data_2,
fixed = parent_0 + log_k_parent_sink + log_k_parent_A1 + log_k_A1_sink ~ 1,
random = pdDiag(parent_0 + log_k_parent_sink + log_k_parent_A1 + log_k_A1_sink ~ 1),
start = mean_dp_sfo_sfo)
+ # augPred does not see to work on this object, so no plot is shown
# The same model fitted with transformed formation fractions does not converge
f_nlme_sfo_sfo_ff <- nlme(value ~ nlme_f_sfo_sfo_ff(name, time,
@@ -118,14 +123,6 @@ plot(augPred(m_nlme, level = 0:1), layout = c(3, 1))
random = pdDiag(parent_0 + log_k_parent + log_k_A1 + f_parent_ilr_1 ~ 1),
start = mean_dp_sfo_sfo_ff)
- # It does converge with this version of reduced random effects
- f_nlme_sfo_sfo_ff <- nlme(value ~ nlme_f_sfo_sfo_ff(name, time,
- parent_0, log_k_parent, log_k_A1, f_parent_ilr_1),
- data = grouped_data_2,
- fixed = parent_0 + log_k_parent + log_k_A1 + f_parent_ilr_1 ~ 1,
- random = pdDiag(parent_0 + log_k_parent ~ 1),
- start = mean_dp_sfo_sfo_ff)
-
f_nlme_fomc_sfo <- nlme(value ~ nlme_f_fomc_sfo(name, time,
parent_0, log_alpha, log_beta, log_k_A1, f_parent_ilr_1),
data = grouped_data_2,
diff --git a/man/nlme.mmkin.Rd b/man/nlme.mmkin.Rd
new file mode 100644
index 00000000..5f937488
--- /dev/null
+++ b/man/nlme.mmkin.Rd
@@ -0,0 +1,72 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/nlme.mmkin.R
+\name{nlme.mmkin}
+\alias{nlme.mmkin}
+\title{Create an nlme model for an mmkin row object}
+\usage{
+\method{nlme}{mmkin}(
+ model,
+ data = sys.frame(sys.parent()),
+ fixed,
+ random = fixed,
+ groups,
+ start,
+ correlation = NULL,
+ weights = NULL,
+ subset,
+ method = c("ML", "REML"),
+ na.action = na.fail,
+ naPattern,
+ control = list(),
+ verbose = FALSE
+)
+}
+\arguments{
+\item{model}{An \code{\link{mmkin}} row object.}
+
+\item{data}{Ignored, data are taken from the mmkin model}
+
+\item{fixed}{Ignored, all degradation parameters fitted in the
+mmkin model are used as fixed parameters}
+
+\item{random}{If not specified, all fixed effects are complemented
+with uncorrelated random effects}
+
+\item{groups}{See the documentation of nlme}
+
+\item{start}{If not specified, mean values of the fitted degradation
+parameters taken from the mmkin object are used}
+
+\item{correlation}{See the documentation of nlme}
+
+\item{weights}{passed to nlme}
+
+\item{subset}{passed to nlme}
+
+\item{method}{passed to nlme}
+
+\item{na.action}{passed to nlme}
+
+\item{naPattern}{passed to nlme}
+
+\item{control}{passed to nlme}
+
+\item{verbose}{passed to nlme}
+}
+\value{
+Upon success, a fitted nlme.mmkin object, which is
+ an nlme object with additional elements
+}
+\description{
+Create an nlme model for an mmkin row object
+}
+\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)
+library(nlme)
+f_nlme <- nlme(f)
+nlme(f, random = parent_0 ~ 1)
+f_nlme <- nlme(f, start = c(parent_0 = 100, log_k_parent_sink = 0.1))
+update(f_nlme, random = parent_0 ~ 1)
+}

Contact - Imprint