#' 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 show_errmin Should the chi2 error level be shown on top of the plots
#' to the left?
#' @param errmin_var The variable for which the FOCUS chi2 error value should
#' be shown.
#' @param errmin_digits The number of significant digits for rounding the FOCUS
#' chi2 error percentage.
#' @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
#' plot(f[, 3:4])
#' library(nlme)
#' f_nlme <- nlme(f)
#'
#' #plot(f_nlme) # too many panels for pkgdown
#' plot(f_nlme, 3:4)
#' @export
plot.nlme.mmkin <- function(x, i = 1:ncol(x$mmkin_orig),
main = "auto", legends = 1,
resplot = c("time", "errmod"),
standardized = FALSE,
show_errmin = TRUE,
errmin_var = "All data", errmin_digits = 3,
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_optim)] <- odeparms_optim
mkinfit_call[["observed"]] <- ds[[a]]
mkinfit_call[["parms.ini"]] <- odeparms
mkinfit_call[["state.ini"]] <- state_ini
mkinfit_call[["control"]] <- list(iter.max = 0)
mkinfit_call[["quiet"]] <- TRUE
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, main = main, legends = legends,
resplot = resplot, standardized = standardized,
show_errmin = show_errmin,
errmin_var = errmin_var, errmin_digits = errmin_digits,
cex = cex,
rel.height.middle = rel.height.middle,
ymax = ymax, ...)
}
