From 8d21c9966305627f9ea8599503fb3709ca916f9b Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Fri, 13 Nov 2020 10:27:00 +0100 Subject: Fixed plotting of non-analytical solutions in saem --- R/plot.mixed.mmkin.R | 8 +++----- R/saem.R | 31 ++++++++++++++++++++----------- R/summary.saem.mmkin.R | 9 +++------ 3 files changed, 26 insertions(+), 22 deletions(-) (limited to 'R') diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R index 0d6b15c0..d8a1c2ac 100644 --- a/R/plot.mixed.mmkin.R +++ b/R/plot.mixed.mmkin.R @@ -75,8 +75,8 @@ plot.mixed.mmkin <- function(x, rownames(degparms_optim) <- ds_names degparms_optim_names <- setdiff(names(fit_1$par), names(fit_1$errparms)) colnames(degparms_optim) <- degparms_optim_names - residual_type = ifelse(standardized, "iwres", "ires") - residuals <- x$so@results@predictions[[residual_type]] + residual_type = ifelse(standardized, "residual", "standardized") + residuals <- x$data[[residual_type]] degparms_pop <- x$so@results@fixed.effects names(degparms_pop) <- degparms_optim_names } @@ -95,9 +95,7 @@ plot.mixed.mmkin <- function(x, odeini_names <- grep("_0$", degparms_all_names, value = TRUE) odeparms_names <- setdiff(degparms_all_names, odeini_names) - residual_type = ifelse(standardized, "iwres", "ires") - - observed <- cbind(x$data, + observed <- cbind(x$data[c("ds", "name", "time", "value")], residual = residuals) solution_type = fit_1$solution_type diff --git a/R/saem.R b/R/saem.R index 377da903..18152baf 100644 --- a/R/saem.R +++ b/R/saem.R @@ -63,7 +63,7 @@ #' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"), #' A1 = mkinsub("SFO")) #' # The following fit uses analytical solutions for SFO-SFO and DFOP-SFO, -#' # and compiled ODEs for FOMC, both are fast +#' # and compiled ODEs for FOMC that are much slower #' f_mmkin <- mmkin(list( #' "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo), #' ds, quiet = TRUE) @@ -75,11 +75,12 @@ #' # We can use print, plot and summary methods to check the results #' print(f_saem_dfop_sfo) #' plot(f_saem_dfop_sfo) -#' summary(f_saem_dfop_sfo, data = FALSE) +#' summary(f_saem_dfop_sfo, data = TRUE) #' #' # Using a single core, the following takes about 6 minutes as we do not have an #' # analytical solution. Using 10 cores it is slower instead of faster -#' #f_saem_fomc <- saem(f_mmkin["FOMC-SFO", ], cores = 1) +#' f_saem_fomc <- saem(f_mmkin["FOMC-SFO", ], cores = 1) +#' plot(f_saem_fomc) #' } #' @export saem <- function(object, control, ...) UseMethod("saem") @@ -103,12 +104,6 @@ saem.mmkin <- function(object, } fit_time <- system.time({ utils::capture.output(f_saemix <- saemix::saemix(m_saemix, d_saemix, control), split = !quiet) - f_pred <- try(saemix::saemix.predict(f_saemix), silent = TRUE) - if (!inherits(f_pred, "try-error")) { - f_saemix <- f_pred - } else { - warning("Creating predictions from the saemix model failed") - } }) if (suppressPlot) { grDevices::dev.off() @@ -121,6 +116,20 @@ saem.mmkin <- function(object, object[[1]]$transform_rates, object[[1]]$transform_fractions) + return_data <- nlme_data(object) + + return_data$predicted <- f_saemix@model@model( + psi = saemix::psi(f_saemix), + id = as.numeric(return_data$ds), + xidep = return_data[c("time", "name")]) + + return_data <- transform(return_data, + residual = predicted - value, + std = sigma_twocomp(predicted, + f_saemix@results@respar[1], f_saemix@results@respar[2])) + return_data <- transform(return_data, + standardized = residual / std) + result <- list( mkinmod = object[[1]]$mkinmod, mmkin = object, @@ -132,7 +141,7 @@ saem.mmkin <- function(object, mean_dp_start = attr(m_saemix, "mean_dp_start"), bparms.optim = bparms_optim, bparms.fixed = object[[1]]$bparms.fixed, - data = nlme_data(object), + data = return_data, err_mod = object[[1]]$err_mod, date.fit = date(), saemixversion = as.character(utils::packageVersion("saemix")), @@ -323,7 +332,7 @@ saemix_model <- function(object, cores = 1, verbose = FALSE, ...) { uid <- unique(id) res_list <- parallel::mclapply(uid, function(i) { - transparms_optim <- psi[i, ] + transparms_optim <- as.numeric(psi[i, ]) names(transparms_optim) <- names(degparms_optim) odeini_optim <- transparms_optim[odeini_optim_parm_names] diff --git a/R/summary.saem.mmkin.R b/R/summary.saem.mmkin.R index 536693a3..97f9f2da 100644 --- a/R/summary.saem.mmkin.R +++ b/R/summary.saem.mmkin.R @@ -151,14 +151,11 @@ summary.saem.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes = err_mod = object$mmkin_orig[[1]]$err_mod object$diffs <- object$mkinmod$diffs - object$print_data <- data + object$print_data <- data # boolean: Should we print the data? so_pred <- object$so@results@predictions - object$data[["observed"]] <- object$data[["value"]] - object$data[["value"]] <- NULL - object$data[["predicted"]] <- so_pred$ipred - object$data[["residual"]] <- so_pred$ires - object$data[["standardized"]] <- so_pred$iwres + names(object$data)[4] <- "observed" # rename value to observed + object$verbose <- verbose object$fixed <- object$mmkin_orig[[1]]$fixed -- cgit v1.2.1