From 26085403289e29259e500282e8e88a5ab00c07a0 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 14 Apr 2020 17:22:00 +0200 Subject: Add a nlme method for mmkin row objects --- NAMESPACE | 2 ++ R/nlme.R | 17 ++++----- R/nlme.mmkin.R | 85 +++++++++++++++++++++++++++++++++++++++++++++ check.log | 102 +++++++++++++++++++++++++++++++++++++++++++++++++++--- man/nlme.Rd | 15 ++++---- man/nlme.mmkin.Rd | 72 ++++++++++++++++++++++++++++++++++++++ 6 files changed, 270 insertions(+), 23 deletions(-) create mode 100644 R/nlme.mmkin.R create mode 100644 man/nlme.mmkin.Rd 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 ’ @@ -41,7 +41,14 @@ Maintainer: ‘Johannes Ranke ’ * 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 ’ * 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 ’ * 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) +} -- cgit v1.2.1