aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
Diffstat (limited to 'R')
-rw-r--r--R/llhist.R11
-rw-r--r--R/mhmkin.R2
-rw-r--r--R/multistart.R48
-rw-r--r--R/parhist.R8
-rw-r--r--R/parms.R (renamed from R/parms.mkinfit.R)13
-rw-r--r--R/saem.R125
6 files changed, 147 insertions, 60 deletions
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("<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)
+}
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)
}

Contact - Imprint