aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2020-04-15 01:27:02 +0200
committerJohannes Ranke <jranke@uni-bremen.de>2020-04-15 01:51:55 +0200
commit082be7baa35d7e8131a70ca8cc061d90e08fa584 (patch)
tree29774198fae54ac13f2a0c2d2fb385a5fc44ab92 /R
parent26085403289e29259e500282e8e88a5ab00c07a0 (diff)
A plot method for nlme.mmkin fits
Diffstat (limited to 'R')
-rw-r--r--R/nlme.mmkin.R9
-rw-r--r--R/plot.nlme.mmkin.R103
2 files changed, 108 insertions, 4 deletions
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, ...)
+
+}

Contact - Imprint