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 --- 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 +++++++++++++++++++++++++++++++++--------------------- 7 files changed, 216 insertions(+), 129 deletions(-) create mode 100644 R/parms.R delete mode 100644 R/parms.mkinfit.R (limited to 'R') 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) } -- cgit v1.2.1