diff options
Diffstat (limited to 'R/plot.mixed.mmkin.R')
| -rw-r--r-- | R/plot.mixed.mmkin.R | 48 |
1 files changed, 46 insertions, 2 deletions
diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R index 4c1f1531..1ac62b07 100644 --- a/R/plot.mixed.mmkin.R +++ b/R/plot.mixed.mmkin.R @@ -2,7 +2,7 @@ utils::globalVariables("ds") #' Plot predictions from a fitted nonlinear mixed model obtained via an mmkin row object #' -#' @param x An object of class [mixed.mmkin], [nlme.mmkin] +#' @param x An object of class [mixed.mmkin], [saem.mmkin] or [nlme.mmkin] #' @param i A numeric index to select datasets for which to plot the individual predictions, #' in case plots get too large #' @inheritParams plot.mkinfit @@ -10,6 +10,10 @@ utils::globalVariables("ds") #' `resplot = "time"`. #' @param pred_over Named list of alternative predictions as obtained #' from [mkinpredict] with a compatible [mkinmod]. +#' @param test_log_parms Passed to [mean_degparms] in the case of an +#' [mixed.mmkin] object +#' @param conf.level Passed to [mean_degparms] in the case of an +#' [mixed.mmkin] object #' @param rel.height.legend The relative height of the legend shown on top #' @param rel.height.bottom The relative height of the bottom plot row #' @param ymax Vector of maximum y axis values @@ -36,9 +40,23 @@ utils::globalVariables("ds") #' #' # For this fit we need to increase pnlsMaxiter, and we increase the #' # tolerance in order to speed up the fit for this example evaluation +#' # It still takes 20 seconds to run #' f_nlme <- nlme(f, control = list(pnlsMaxIter = 120, tolerance = 1e-3)) #' plot(f_nlme) #' +#' f_saem <- saem(f, transformations = "saemix") +#' plot(f_saem) +#' +#' f_obs <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, error_model = "obs") +#' f_nlmix <- nlmix(f_obs) +#' plot(f_nlmix) +#' +#' # We can overlay the two variants if we generate predictions +#' pred_nlme <- mkinpredict(dfop_sfo, +#' f_nlme$bparms.optim[-1], +#' c(parent = f_nlme$bparms.optim[[1]], A1 = 0), +#' seq(0, 180, by = 0.2)) +#' plot(f_saem, pred_over = list(nlme = pred_nlme)) #' } #' @export plot.mixed.mmkin <- function(x, @@ -49,6 +67,8 @@ plot.mixed.mmkin <- function(x, xlim = range(x$data$time), resplot = c("predicted", "time"), pred_over = NULL, + test_log_parms = FALSE, + conf.level = 0.6, ymax = "auto", maxabs = "auto", ncol.legend = ifelse(length(i) <= 3, length(i) + 1, ifelse(length(i) <= 8, 3, 4)), nrow.legend = ceiling((length(i) + 1) / ncol.legend), @@ -67,7 +87,7 @@ plot.mixed.mmkin <- function(x, backtransform = TRUE if (identical(class(x), "mixed.mmkin")) { - degparms_pop <- mean_degparms(x$mmkin) + degparms_pop <- mean_degparms(x$mmkin, test_log_parms = test_log_parms, conf.level = conf.level) degparms_tmp <- parms(x$mmkin, transformed = TRUE) degparms_i <- as.data.frame(t(degparms_tmp[setdiff(rownames(degparms_tmp), names(fit_1$errparms)), ])) @@ -82,6 +102,30 @@ plot.mixed.mmkin <- function(x, type = ifelse(standardized, "pearson", "response")) } + if (inherits(x, "saem.mmkin")) { + if (x$transformations == "saemix") backtransform = FALSE + degparms_i <- saemix::psi(x$so) + rownames(degparms_i) <- ds_names + degparms_i_names <- setdiff(x$so@results@name.fixed, names(fit_1$errparms)) + colnames(degparms_i) <- degparms_i_names + residual_type = ifelse(standardized, "standardized", "residual") + residuals <- x$data[[residual_type]] + degparms_pop <- x$so@results@fixed.effects + names(degparms_pop) <- degparms_i_names + } + + if (inherits(x, "nlmixr.mmkin")) { + eta_i <- random.effects(x$nm)[-1] + names(eta_i) <- gsub("^eta.", "", names(eta_i)) + degparms_i <- eta_i + degparms_pop <- x$nm$theta + for (parm_name in names(degparms_i)) { + degparms_i[parm_name] <- eta_i[parm_name] + degparms_pop[parm_name] + } + residual_type = ifelse(standardized, "standardized", "residual") + residuals <- x$data[[residual_type]] + } + degparms_fixed <- fit_1$fixed$value names(degparms_fixed) <- rownames(fit_1$fixed) degparms_all <- cbind(as.matrix(degparms_i), |