From 082be7baa35d7e8131a70ca8cc061d90e08fa584 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 15 Apr 2020 01:27:02 +0200 Subject: A plot method for nlme.mmkin fits --- NAMESPACE | 2 + NEWS.md | 2 + R/nlme.mmkin.R | 9 +- R/plot.nlme.mmkin.R | 103 ++++++++++++ _pkgdown.yml | 5 + check.log | 102 +----------- docs/news/index.html | 1 + docs/reference/index.html | 23 ++- docs/reference/nlme.html | 32 ++-- docs/reference/nlme.mmkin.html | 301 +++++++++++++++++++++++++++++++++++ docs/reference/plot.nlme.mmkin-1.png | Bin 0 -> 32297 bytes docs/reference/plot.nlme.mmkin.html | 252 +++++++++++++++++++++++++++++ docs/sitemap.xml | 6 + man/nlme.Rd | 2 +- man/nlme.mmkin.Rd | 11 +- man/plot.nlme.mmkin.Rd | 67 ++++++++ 16 files changed, 800 insertions(+), 118 deletions(-) create mode 100644 R/plot.nlme.mmkin.R create mode 100644 docs/reference/nlme.mmkin.html create mode 100644 docs/reference/plot.nlme.mmkin-1.png create mode 100644 docs/reference/plot.nlme.mmkin.html create mode 100644 man/plot.nlme.mmkin.Rd diff --git a/NAMESPACE b/NAMESPACE index 5f305a3d..fce1c5f7 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -18,6 +18,7 @@ S3method(parms,mkinfit) S3method(plot,mkinfit) S3method(plot,mmkin) S3method(plot,nafta) +S3method(plot,nlme.mmkin) S3method(print,mkinds) S3method(print,mkinmod) S3method(print,nafta) @@ -87,6 +88,7 @@ importFrom(stats,AIC) importFrom(stats,BIC) importFrom(stats,aggregate) importFrom(stats,coef) +importFrom(stats,coefficients) importFrom(stats,cov2cor) importFrom(stats,dist) importFrom(stats,dnorm) diff --git a/NEWS.md b/NEWS.md index 46d2d711..a4976488 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,7 @@ # mkin 0.9.49.10 (unreleased) +- 'nlme.mmkin': An nlme method for mmkin row objects and an associated class with plot method + - 'mean_degparms, nlme_data, nlme_function': Three new functions to facilitate building nlme models from mmkin row objects - 'endpoints': Don't return the SFORB list component if it's empty. This reduces distraction and complies with the documentation diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R index 784ba609..2ee46f33 100644 --- a/R/nlme.mmkin.R +++ b/R/nlme.mmkin.R @@ -22,8 +22,8 @@ #' @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 +#' @return Upon success, a fitted nlme.mmkin object, which is an nlme object +#' with additional elements #' @export #' @seealso \code{\link{nlme_function}} #' @examples @@ -54,7 +54,7 @@ nlme.mmkin <- function(model, data = sys.frame(sys.parent()), } deg_func <- nlme_function(model) - assign("deg_func", deg_func, parent.frame()) + assign("deg_func", deg_func, globalenv()) # specify the model formula this_model_text <- paste0("value ~ deg_func(", @@ -79,7 +79,8 @@ nlme.mmkin <- function(model, data = sys.frame(sys.parent()), } val <- do.call("nlme.formula", thisCall) - return(val) + val$mmkin_orig <- model class(val) <- c("nlme.mmkin", "nlme", "lme") + return(val) } diff --git a/R/plot.nlme.mmkin.R b/R/plot.nlme.mmkin.R new file mode 100644 index 00000000..ef6d100a --- /dev/null +++ b/R/plot.nlme.mmkin.R @@ -0,0 +1,103 @@ +#' Plot a fitted nonlinear mixed model obtained via an mmkin row object +#' +#' @param x An object of class \code{\link{nlme.mmkin}} +#' @param i A numeric index to select datasets for which to plot the nlme fit, +#' in case plots get too large +#' @param main The main title placed on the outer margin of the plot. +#' @param legends An index for the fits for which legends should be shown. +#' @param resplot Should the residuals plotted against time, using +#' \code{\link{mkinresplot}}, or as squared residuals against predicted +#' values, with the error model, using \code{\link{mkinerrplot}}. +#' @param standardized Should the residuals be standardized? This option +#' is passed to \code{\link{mkinresplot}}, it only takes effect if +#' `resplot = "time"`. +#' @param cex Passed to the plot functions and \code{\link{mtext}}. +#' @param rel.height.middle The relative height of the middle plot, if more +#' than two rows of plots are shown. +#' @param ymax Maximum y axis value for \code{\link{plot.mkinfit}}. +#' @param \dots Further arguments passed to \code{\link{plot.mkinfit}} and +#' \code{\link{mkinresplot}}. +#' @importFrom stats coefficients +#' @return The function is called for its side effect. +#' @author Johannes Ranke +#' @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) +#' #plot(f) # too many panels for pkgdown +#' library(nlme) +#' f_nlme <- nlme(f) +#' +#' #plot(f_nlme) # too many panels for pkgdown +#' plot(f_nlme, 1:2) +#' @export +plot.nlme.mmkin <- function(x, i = 1:ncol(x$mmkin_orig), + main = "auto", legends = 1, + resplot = c("time", "errmod"), + standardized = FALSE, + cex = 0.7, rel.height.middle = 0.9, + ymax = "auto", ...) +{ + + degparms_optim_nlme <- coefficients(x) + degparms_optim_names <- names(degparms_optim_nlme) + + odeini_optim_names <- grep("_0$", degparms_optim_names, value = TRUE) + odeparms_optim_names <- setdiff(degparms_optim_names, odeini_optim_names) + + fit_1 <- x$mmkin_orig[[1]] + + mkinfit_call <- as.list(fit_1$call)[-1] + mkinfit_call[["mkinmod"]] <- fit_1$mkinmod + + ds <- lapply(x$mmkin_orig, function(x) { + data.frame(name = x$data$variable, + time = x$data$time, + value = x$data$observed) + }) + + # This takes quite some time. This could be greatly reduced + # if the plot.mkinfit code would be imported and adapted, + # allowing also to overly plots of mmkin fits and nlme fits + mmkin_nlme <- lapply(i, function(a) { + + degparms_optim <- as.numeric(degparms_optim_nlme[a, ]) + names(degparms_optim) <- degparms_optim_names + + odeini_optim <- degparms_optim[odeini_optim_names] + names(odeini_optim) <- gsub("_0$", "", names(odeini_optim)) + + odeparms_optim_trans <- degparms_optim[odeparms_optim_names] + odeparms_optim <- backtransform_odeparms(odeparms_optim_trans, + fit_1$mkinmod, + transform_rates = fit_1$transform_rates, + transform_fractions = fit_1$transform_fractions) + + fit_a <- x$mmkin_orig[[a]] + + state_ini <- fit_a$bparms.state + state_ini[names(odeini_optim)] <- odeini_optim + + odeparms <- fit_a$bparms.ode + odeparms[names(odeparms)] <- odeparms_optim + + mkinfit_call[["observed"]] <- ds[[a]] + mkinfit_call[["parms.ini"]] <- odeparms + mkinfit_call[["state.ini"]] <- state_ini + + mkinfit_call[["control"]] <- list(iter.max = 1) + + res <- suppressWarnings(do.call("mkinfit", mkinfit_call)) + return(res) + }) + + # Set dimensions with names and the class (mmkin) + attributes(mmkin_nlme) <- attributes(x$mmkin_orig[, i]) + + plot(mmkin_nlme[, i], main = main, legends = legends, + resplot = resplot, standardized = standardized, + show_errmin = FALSE, cex = cex, + rel.height.middle = rel.height.middle, + ymax = ymax, ...) + +} diff --git a/_pkgdown.yml b/_pkgdown.yml index 0cecdc30..2ede6857 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -33,6 +33,11 @@ reference: - "`[.mmkin`" - plot.mmkin - AIC.mmkin + - title: Mixed models + desc: Create and work with nonlinear mixed models + contents: + - nlme.mmkin + - plot.nlme.mmkin - nlme_function - title: Datasets and known results contents: diff --git a/check.log b/check.log index fc06bc0d..43cc5f64 100644 --- a/check.log +++ b/check.log @@ -41,14 +41,7 @@ Maintainer: ‘Johannes Ranke ’ * checking S3 generic/method consistency ... OK * checking replacement functions ... OK * checking foreign function calls ... 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 R code for possible problems ... OK * checking Rd files ... OK * checking Rd metadata ... OK * checking Rd line widths ... OK @@ -63,91 +56,11 @@ to your NAMESPACE file. * checking data for ASCII and uncompressed saves ... OK * checking installed files from ‘inst/doc’ ... OK * checking files in ‘vignettes’ ... 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 examples ... NOTE +Examples with CPU or elapsed time > 5s + user system elapsed +nlme.mmkin 4.801 1.040 4.497 +plot.nlme.mmkin 4.751 0.285 4.680 * checking for unstated dependencies in ‘tests’ ... OK * checking tests ... SKIPPED * checking for unstated dependencies in vignettes ... OK @@ -157,8 +70,9 @@ Execution halted * checking for detritus in the temp directory ... OK * DONE -Status: 1 ERROR, 1 NOTE +Status: 1 NOTE See ‘/home/jranke/git/mkin/mkin.Rcheck/00check.log’ for details. + diff --git a/docs/news/index.html b/docs/news/index.html index c8d7ce3d..5fd97344 100644 --- a/docs/news/index.html +++ b/docs/news/index.html @@ -134,6 +134,7 @@ mkin 0.9.49.10 (unreleased) Unreleased