From 3529f5ff498d7d054c7b1911ddfc4b242902b85d Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 28 Sep 2022 16:34:57 +0200 Subject: Fix handling of multistart fits with failures --- NAMESPACE | 3 + R/llhist.R | 11 +- R/mhmkin.R | 2 +- R/multistart.R | 48 +++++++- R/parhist.R | 8 +- R/parms.R | 82 ++++++++++++++ R/parms.mkinfit.R | 69 ------------ R/saem.R | 125 +++++++++++++-------- log/check.log | 16 +-- man/multistart.Rd | 3 - man/parms.Rd | 8 +- vignettes/web_only/multistart.R | 26 ----- .../figure-html/unnamed-chunk-2-1.png | Bin 0 -> 76647 bytes .../figure-html/unnamed-chunk-3-1.png | Bin 0 -> 70831 bytes 14 files changed, 228 insertions(+), 173 deletions(-) create mode 100644 R/parms.R delete mode 100644 R/parms.mkinfit.R delete mode 100644 vignettes/web_only/multistart.R create mode 100644 vignettes/web_only/multistart_files/figure-html/unnamed-chunk-2-1.png create mode 100644 vignettes/web_only/multistart_files/figure-html/unnamed-chunk-3-1.png diff --git a/NAMESPACE b/NAMESPACE index bd12d0db..388a16fd 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -13,6 +13,8 @@ S3method(aw,multistart) S3method(confint,mkinfit) S3method(convergence,mhmkin) S3method(convergence,mmkin) +S3method(convergence,multistart) +S3method(convergence,multistart.saem.mmkin) S3method(f_time_norm_focus,mkindsg) S3method(f_time_norm_focus,numeric) S3method(illparms,mhmkin) @@ -42,6 +44,7 @@ S3method(plot,mmkin) S3method(plot,nafta) S3method(print,convergence.mhmkin) S3method(print,convergence.mmkin) +S3method(print,convergence.multistart) S3method(print,illparms.mhmkin) S3method(print,illparms.mmkin) S3method(print,mhmkin) diff --git a/R/llhist.R b/R/llhist.R index b7f6de21..9ddf5b10 100644 --- a/R/llhist.R +++ b/R/llhist.R @@ -16,7 +16,16 @@ llhist <- function(object, breaks = "Sturges", lpos = "topleft", main = "", ...) oldpar <- par(no.readonly = TRUE) on.exit(par(oldpar, no.readonly = TRUE)) - ll <- sapply(object, logLik) + if (inherits(object, "multistart.saem.mmkin")) { + llfunc <- function(object) { + if (inherits(object$so, "try-error")) return(NA) + else return(logLik(object$so)) + } + } else { + stop("llhist is only implemented for multistart.saem.mmkin objects") + } + + ll <- stats::na.omit(sapply(object, llfunc)) kde <- KernSmooth::bkde(ll) par(las = 1) diff --git a/R/mhmkin.R b/R/mhmkin.R index b72ae318..a1475ef9 100644 --- a/R/mhmkin.R +++ b/R/mhmkin.R @@ -179,7 +179,7 @@ AIC.mhmkin <- function(object, ..., k = 2) { BIC.mhmkin <- function(object, ...) { res <- sapply(object, function(x) { if (inherits(x, "try-error")) return(NA) - else return(BIC(x$so, k = k)) + else return(BIC(x$so)) }) dim(res) <- dim(object) dimnames(res) <- dimnames(object) diff --git a/R/multistart.R b/R/multistart.R index cc55feae..b65c0bee 100644 --- a/R/multistart.R +++ b/R/multistart.R @@ -92,15 +92,51 @@ multistart.saem.mmkin <- function(object, n = 50, cores = 1, return(res) } -#' @rdname multistart #' @export -print.multistart <- function(x, ...) { - cat("Multistart object with", length(x), "fits of the following type:\n\n") - print(x[[1]]) +convergence.multistart <- function(object, ...) { + all_summary_warnings <- character() + + result <- lapply(object, + function(fit) { + if (inherits(fit, "try-error")) return("E") + else { + return("OK") + } + }) + result <- unlist(result) + + class(result) <- "convergence.multistart" + return(result) +} + +#' @export +convergence.multistart.saem.mmkin <- function(object, ...) { + all_summary_warnings <- character() + + result <- lapply(object, + function(fit) { + if (inherits(fit$so, "try-error")) return("E") + else { + return("OK") + } + }) + result <- unlist(result) + + class(result) <- "convergence.multistart" + return(result) +} + +#' @export +print.convergence.multistart <- function(x, ...) { + class(x) <- NULL + print(table(x, dnn = NULL)) + if (any(x == "OK")) cat("OK: Fit terminated successfully\n") + if (any(x == "E")) cat("E: Error\n") } #' @rdname multistart #' @export -parms.multistart <- function(object, ...) { - t(sapply(object, parms)) +print.multistart <- function(x, ...) { + cat(" object with", length(x), "fits:\n") + print(convergence(x)) } diff --git a/R/parhist.R b/R/parhist.R index 69aafe02..10730873 100644 --- a/R/parhist.R +++ b/R/parhist.R @@ -29,20 +29,20 @@ parhist <- function(object, lpos = "bottomleft", main = "", ...) { degparm_index <- which(names(orig_parms) %in% degparm_names_transformed) orig_parms[degparm_names_transformed] <- backtransform_odeparms( orig_parms[degparm_names_transformed], - orig$mmkin$mkinmod, + orig$mmkin[[1]]$mkinmod, transform_rates = orig$mmkin[[1]]$transform_rates, transform_fractions = orig$mmkin[[1]]$transform_fractions) start_parms <- backtransform_odeparms(start_parms, - orig$mmkin$mkinmod, + orig$mmkin[[1]]$mkinmod, transform_rates = orig$mmkin[[1]]$transform_rates, transform_fractions = orig$mmkin[[1]]$transform_fractions) degparm_names <- names(start_parms) names(orig_parms) <- c(degparm_names, names(orig_parms[-degparm_index])) - + all_parms[, degparm_names_transformed] <- t(apply(all_parms[, degparm_names_transformed], 1, backtransform_odeparms, - orig$mmkin$mkinmod, + orig$mmkin[[1]]$mkinmod, transform_rates = orig$mmkin[[1]]$transform_rates, transform_fractions = orig$mmkin[[1]]$transform_fractions)) colnames(all_parms)[1:length(degparm_names)] <- degparm_names diff --git a/R/parms.R b/R/parms.R new file mode 100644 index 00000000..bd4e479b --- /dev/null +++ b/R/parms.R @@ -0,0 +1,82 @@ +#' Extract model parameters +#' +#' This function returns degradation model parameters as well as error +#' model parameters per default, in order to avoid working with a fitted model +#' without considering the error structure that was assumed for the fit. +#' +#' @param object A fitted model object. +#' @param \dots Not used +#' @return Depending on the object, a numeric vector of fitted model parameters, +#' a matrix (e.g. for mmkin row objects), or a list of matrices (e.g. for +#' mmkin objects with more than one row). +#' @seealso [saem], [multistart] +#' @examples +#' # mkinfit objects +#' fit <- mkinfit("SFO", FOCUS_2006_C, quiet = TRUE) +#' parms(fit) +#' parms(fit, transformed = TRUE) +#' +#' # mmkin objects +#' ds <- lapply(experimental_data_for_UBA_2019[6:10], +#' function(x) subset(x$data[c("name", "time", "value")])) +#' names(ds) <- paste("Dataset", 6:10) +#' \dontrun{ +#' fits <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE, cores = 1) +#' parms(fits["SFO", ]) +#' parms(fits[, 2]) +#' parms(fits) +#' parms(fits, transformed = TRUE) +#' } +#' @export +parms <- function(object, ...) +{ + UseMethod("parms", object) +} + +#' @param transformed Should the parameters be returned as used internally +#' during the optimisation? +#' @param errparms Should the error model parameters be returned +#' in addition to the degradation parameters? +#' @rdname parms +#' @export +parms.mkinfit <- function(object, transformed = FALSE, errparms = TRUE, ...) +{ + res <- if (transformed) object$par + else c(object$bparms.optim, object$errparms) + if (!errparms) { + res[setdiff(names(res), names(object$errparms))] + } + else return(res) +} + +#' @rdname parms +#' @export +parms.mmkin <- function(object, transformed = FALSE, errparms = TRUE, ...) +{ + if (nrow(object) == 1) { + res <- sapply(object, parms, transformed = transformed, + errparms = errparms, ...) + colnames(res) <- colnames(object) + } else { + res <- list() + for (i in 1:nrow(object)) { + res[[i]] <- parms(object[i, ], transformed = transformed, + errparms = errparms, ...) + } + names(res) <- rownames(object) + } + return(res) +} + +#' @param exclude_failed For [multistart] objects, should rows for failed fits +#' be removed from the returned parameter matrix? +#' @rdname parms +#' @export +parms.multistart <- function(object, exclude_failed = TRUE, ...) { + res <- t(sapply(object, parms)) + successful <- which(!is.na(res[, 1])) + first_success <- successful[1] + colnames(res) <- names(parms(object[[first_success]])) + if (exclude_failed) res <- res[successful, ] + return(res) +} diff --git a/R/parms.mkinfit.R b/R/parms.mkinfit.R deleted file mode 100644 index 83766355..00000000 --- a/R/parms.mkinfit.R +++ /dev/null @@ -1,69 +0,0 @@ -#' Extract model parameters -#' -#' This function returns degradation model parameters as well as error -#' model parameters per default, in order to avoid working with a fitted model -#' without considering the error structure that was assumed for the fit. -#' -#' @param object A fitted model object. -#' @param \dots Not used -#' @return Depending on the object, a numeric vector of fitted model parameters, -#' a matrix (e.g. for mmkin row objects), or a list of matrices (e.g. for -#' mmkin objects with more than one row). -#' @seealso [saem], [multistart] -#' @examples -#' # mkinfit objects -#' fit <- mkinfit("SFO", FOCUS_2006_C, quiet = TRUE) -#' parms(fit) -#' parms(fit, transformed = TRUE) -#' -#' # mmkin objects -#' ds <- lapply(experimental_data_for_UBA_2019[6:10], -#' function(x) subset(x$data[c("name", "time", "value")])) -#' names(ds) <- paste("Dataset", 6:10) -#' \dontrun{ -#' fits <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE, cores = 1) -#' parms(fits["SFO", ]) -#' parms(fits[, 2]) -#' parms(fits) -#' parms(fits, transformed = TRUE) -#' } -#' @export -parms <- function(object, ...) -{ - UseMethod("parms", object) -} - -#' @param transformed Should the parameters be returned as used internally -#' during the optimisation? -#' @param errparms Should the error model parameters be returned -#' in addition to the degradation parameters? -#' @rdname parms -#' @export -parms.mkinfit <- function(object, transformed = FALSE, errparms = TRUE, ...) -{ - res <- if (transformed) object$par - else c(object$bparms.optim, object$errparms) - if (!errparms) { - res[setdiff(names(res), names(object$errparms))] - } - else return(res) -} - -#' @rdname parms -#' @export -parms.mmkin <- function(object, transformed = FALSE, errparms = TRUE, ...) -{ - if (nrow(object) == 1) { - res <- sapply(object, parms, transformed = transformed, - errparms = errparms, ...) - colnames(res) <- colnames(object) - } else { - res <- list() - for (i in 1:nrow(object)) { - res[[i]] <- parms(object[i, ], transformed = transformed, - errparms = errparms, ...) - } - names(res) <- rownames(object) - } - return(res) -} diff --git a/R/saem.R b/R/saem.R index 0d0d9b8a..8929cf6f 100644 --- a/R/saem.R +++ b/R/saem.R @@ -137,9 +137,16 @@ saem.mmkin <- function(object, transformations = transformations, ...) d_saemix <- saemix_data(object, verbose = verbose) + fit_failed <- FALSE + FIM_failed <- NULL fit_time <- system.time({ - utils::capture.output(f_saemix <- saemix::saemix(m_saemix, d_saemix, control), split = !quiet) - FIM_failed <- NULL + utils::capture.output(f_saemix <- try(saemix::saemix(m_saemix, d_saemix, control)), split = !quiet) + if (inherits(f_saemix, "try-error")) fit_failed <- TRUE + }) + + return_data <- nlme_data(object) + + if (!fit_failed) { if (any(is.na(f_saemix@results@se.fixed))) FIM_failed <- c(FIM_failed, "fixed effects") if (any(is.na(c(f_saemix@results@se.omega, f_saemix@results@se.respar)))) { FIM_failed <- c(FIM_failed, "random effects and residual error parameters") @@ -147,38 +154,37 @@ saem.mmkin <- function(object, if (!is.null(FIM_failed) & fail_with_errors) { stop("Could not invert FIM for ", paste(FIM_failed, collapse = " and ")) } - }) - transparms_optim <- f_saemix@results@fixed.effects - names(transparms_optim) <- f_saemix@results@name.fixed + transparms_optim <- f_saemix@results@fixed.effects + names(transparms_optim) <- f_saemix@results@name.fixed - if (transformations == "mkin") { - bparms_optim <- backtransform_odeparms(transparms_optim, - object[[1]]$mkinmod, - object[[1]]$transform_rates, - object[[1]]$transform_fractions) - } else { - bparms_optim <- transparms_optim - } + if (transformations == "mkin") { + bparms_optim <- backtransform_odeparms(transparms_optim, + object[[1]]$mkinmod, + object[[1]]$transform_rates, + object[[1]]$transform_fractions) + } else { + bparms_optim <- transparms_optim + } - return_data <- nlme_data(object) - saemix_data_ds <- f_saemix@data@data$ds - mkin_ds_order <- as.character(unique(return_data$ds)) - saemix_ds_order <- unique(saemix_data_ds) - - psi <- saemix::psi(f_saemix) - rownames(psi) <- saemix_ds_order - return_data$predicted <- f_saemix@model@model( - psi = psi[mkin_ds_order, ], - id = as.numeric(return_data$ds), - xidep = return_data[c("time", "name")]) - - return_data <- transform(return_data, - residual = value - predicted, - std = sigma_twocomp(predicted, - f_saemix@results@respar[1], f_saemix@results@respar[2])) - return_data <- transform(return_data, - standardized = residual / std) + saemix_data_ds <- f_saemix@data@data$ds + mkin_ds_order <- as.character(unique(return_data$ds)) + saemix_ds_order <- unique(saemix_data_ds) + + psi <- saemix::psi(f_saemix) + rownames(psi) <- saemix_ds_order + return_data$predicted <- f_saemix@model@model( + psi = psi[mkin_ds_order, ], + id = as.numeric(return_data$ds), + xidep = return_data[c("time", "name")]) + + return_data <- transform(return_data, + residual = value - predicted, + 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, @@ -187,15 +193,13 @@ saem.mmkin <- function(object, transformations = transformations, transform_rates = object[[1]]$transform_rates, transform_fractions = object[[1]]$transform_fractions, + sm = m_saemix, so = f_saemix, call = call, time = fit_time, mean_dp_start = attr(m_saemix, "mean_dp_start"), - bparms.optim = bparms_optim, bparms.fixed = object[[1]]$bparms.fixed, data = return_data, - mkin_ds_order = mkin_ds_order, - saemix_ds_order = saemix_ds_order, err_mod = object[[1]]$err_mod, date.fit = date(), saemixversion = as.character(utils::packageVersion("saemix")), @@ -203,6 +207,12 @@ saem.mmkin <- function(object, Rversion = paste(R.version$major, R.version$minor, sep=".") ) + if (!fit_failed) { + result$mkin_ds_order <- mkin_ds_order + result$saemix_ds_order <- saemix_ds_order + result$bparms.optim <- bparms_optim + } + class(result) <- c("saem.mmkin", "mixed.mmkin") return(result) } @@ -222,16 +232,20 @@ print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) { length(unique(x$data$name)), "variable(s) grouped in", length(unique(x$data$ds)), "datasets\n") - cat("\nLikelihood computed by importance sampling\n") - print(data.frame( - AIC = AIC(x$so, type = "is"), - BIC = BIC(x$so, type = "is"), - logLik = logLik(x$so, type = "is"), - row.names = " "), digits = digits) - - cat("\nFitted parameters:\n") - conf.int <- parms(x, ci = TRUE) - print(conf.int, digits = digits) + if (inherits(x$so, "try-error")) { + cat("\nFit did not terminate successfully\n") + } else { + cat("\nLikelihood computed by importance sampling\n") + print(data.frame( + AIC = AIC(x$so, type = "is"), + BIC = BIC(x$so, type = "is"), + logLik = logLik(x$so, type = "is"), + row.names = " "), digits = digits) + + cat("\nFitted parameters:\n") + conf.int <- parms(x, ci = TRUE) + print(conf.int, digits = digits) + } invisible(x) } @@ -618,12 +632,27 @@ update.saem.mmkin <- function(object, ..., evaluate = TRUE) { #' @param ci Should a matrix with estimates and confidence interval boundaries #' be returned? If FALSE (default), a vector of estimates is returned. parms.saem.mmkin <- function(object, ci = FALSE, ...) { - conf.int <- object$so@results@conf.int[c("estimate", "lower", "upper")] - rownames(conf.int) <- object$so@results@conf.int[["name"]] - conf.int.var <- grepl("^Var\\.", rownames(conf.int)) - conf.int <- conf.int[!conf.int.var, ] + cov.mod <- object$sm@covariance.model + n_cov_mod_parms <- sum(cov.mod[upper.tri(cov.mod, diag = TRUE)]) + n_parms <- length(object$sm@name.modpar) + + n_cov_mod_parms + + length(object$sm@name.sigma) + + if (inherits(object$so, "try-error")) { + conf.int <- matrix(rep(NA, 3 * n_parms), ncol = 3) + colnames(conf.int) <- c("estimate", "lower", "upper") + } else { + conf.int <- object$so@results@conf.int[c("estimate", "lower", "upper")] + rownames(conf.int) <- object$so@results@conf.int[["name"]] + conf.int.var <- grepl("^Var\\.", rownames(conf.int)) + conf.int <- conf.int[!conf.int.var, ] + conf.int.cov <- grepl("^Cov\\.", rownames(conf.int)) + conf.int <- conf.int[!conf.int.cov, ] + } estimate <- conf.int[, "estimate"] + names(estimate) <- rownames(conf.int) + if (ci) return(conf.int) else return(estimate) } diff --git a/log/check.log b/log/check.log index 60298c09..8fce3fd1 100644 --- a/log/check.log +++ b/log/check.log @@ -48,15 +48,7 @@ Maintainer: ‘Johannes Ranke ’ * checking Rd cross-references ... OK * checking for missing documentation entries ... OK * checking for code/documentation mismatches ... OK -* checking Rd \usage sections ... WARNING -Undocumented arguments in documentation object 'multistart' - ‘x’ - -Functions with \usage entries need to have the appropriate \alias -entries, and all their arguments documented. -The \usage entries must correspond to syntactically valid R code. -See chapter ‘Writing R documentation files’ in the ‘Writing R -Extensions’ manual. +* checking Rd \usage sections ... OK * checking Rd contents ... OK * checking for unstated dependencies in examples ... OK * checking contents of ‘data’ directory ... OK @@ -77,9 +69,5 @@ Extensions’ manual. * checking for detritus in the temp directory ... OK * DONE -Status: 1 WARNING -See - ‘/home/jranke/git/mkin/mkin.Rcheck/00check.log’ -for details. - +Status: OK diff --git a/man/multistart.Rd b/man/multistart.Rd index aa352267..78ff4614 100644 --- a/man/multistart.Rd +++ b/man/multistart.Rd @@ -4,7 +4,6 @@ \alias{multistart} \alias{multistart.saem.mmkin} \alias{print.multistart} -\alias{parms.multistart} \title{Perform a hierarchical model fit with multiple starting values} \usage{ multistart( @@ -18,8 +17,6 @@ multistart( \method{multistart}{saem.mmkin}(object, n = 50, cores = 1, cluster = NULL, ...) \method{print}{multistart}(x, ...) - -\method{parms}{multistart}(object, ...) } \arguments{ \item{object}{The fit object to work with} diff --git a/man/parms.Rd b/man/parms.Rd index acae2d91..5c0e8895 100644 --- a/man/parms.Rd +++ b/man/parms.Rd @@ -1,9 +1,10 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/parms.mkinfit.R +% Please edit documentation in R/parms.R \name{parms} \alias{parms} \alias{parms.mkinfit} \alias{parms.mmkin} +\alias{parms.multistart} \title{Extract model parameters} \usage{ parms(object, ...) @@ -11,6 +12,8 @@ parms(object, ...) \method{parms}{mkinfit}(object, transformed = FALSE, errparms = TRUE, ...) \method{parms}{mmkin}(object, transformed = FALSE, errparms = TRUE, ...) + +\method{parms}{multistart}(object, exclude_failed = TRUE, ...) } \arguments{ \item{object}{A fitted model object.} @@ -22,6 +25,9 @@ during the optimisation?} \item{errparms}{Should the error model parameters be returned in addition to the degradation parameters?} + +\item{exclude_failed}{For \link{multistart} objects, should rows for failed fits +be removed from the returned parameter matrix?} } \value{ Depending on the object, a numeric vector of fitted model parameters, diff --git a/vignettes/web_only/multistart.R b/vignettes/web_only/multistart.R deleted file mode 100644 index b995c545..00000000 --- a/vignettes/web_only/multistart.R +++ /dev/null @@ -1,26 +0,0 @@ -## ----------------------------------------------------------------------------- -library(mkin) -dmta_ds <- lapply(1:7, function(i) { - ds_i <- dimethenamid_2018$ds[[i]]$data - ds_i[ds_i$name == "DMTAP", "name"] <- "DMTA" - ds_i$time <- ds_i$time * dimethenamid_2018$f_time_norm[i] - ds_i -}) -names(dmta_ds) <- sapply(dimethenamid_2018$ds, function(ds) ds$title) -dmta_ds[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]]) -dmta_ds[["Elliot 1"]] <- dmta_ds[["Elliot 2"]] <- NULL - -## ----------------------------------------------------------------------------- -f_mmkin <- mmkin("DFOP", dmta_ds, error_model = "tc", cores = 7, quiet = TRUE) -f_saem_full <- saem(f_mmkin) -f_saem_full_multi <- multistart(f_saem_full, n = 16, cores = 16) -parhist(f_saem_full_multi) - -## ----------------------------------------------------------------------------- -f_saem_reduced <- update(f_saem_full, covariance.model = diag(c(1, 1, 0, 1))) -f_saem_reduced_multi <- multistart(f_saem_reduced, n = 16, cores = 16) -parhist(f_saem_reduced_multi, lpos = "topright") - -## ----------------------------------------------------------------------------- -llhist(f_saem_reduced_multi) - diff --git a/vignettes/web_only/multistart_files/figure-html/unnamed-chunk-2-1.png b/vignettes/web_only/multistart_files/figure-html/unnamed-chunk-2-1.png new file mode 100644 index 00000000..17168f74 Binary files /dev/null and b/vignettes/web_only/multistart_files/figure-html/unnamed-chunk-2-1.png differ diff --git a/vignettes/web_only/multistart_files/figure-html/unnamed-chunk-3-1.png b/vignettes/web_only/multistart_files/figure-html/unnamed-chunk-3-1.png new file mode 100644 index 00000000..17d99e4e Binary files /dev/null and b/vignettes/web_only/multistart_files/figure-html/unnamed-chunk-3-1.png differ -- cgit v1.2.1