diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/llhist.R | 11 | ||||
| -rw-r--r-- | R/mhmkin.R | 2 | ||||
| -rw-r--r-- | R/multistart.R | 48 | ||||
| -rw-r--r-- | R/parhist.R | 8 | ||||
| -rw-r--r-- | R/parms.R (renamed from R/parms.mkinfit.R) | 13 | ||||
| -rw-r--r-- | R/saem.R | 125 |
6 files changed, 147 insertions, 60 deletions
@@ -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) @@ -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("<multistart> 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.mkinfit.R b/R/parms.R index 83766355..bd4e479b 100644 --- a/R/parms.mkinfit.R +++ b/R/parms.R @@ -67,3 +67,16 @@ parms.mmkin <- function(object, transformed = FALSE, errparms = TRUE, ...) } 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) +} @@ -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) } |