aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--NAMESPACE3
-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
-rw-r--r--log/check.log16
-rw-r--r--man/multistart.Rd3
-rw-r--r--man/parms.Rd8
-rw-r--r--vignettes/web_only/multistart.R26
-rw-r--r--vignettes/web_only/multistart_files/figure-html/unnamed-chunk-2-1.pngbin0 -> 76647 bytes
-rw-r--r--vignettes/web_only/multistart_files/figure-html/unnamed-chunk-3-1.pngbin0 -> 70831 bytes
13 files changed, 159 insertions, 104 deletions
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("<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)
}
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 <johannes.ranke@jrwb.de>’
* 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
--- /dev/null
+++ b/vignettes/web_only/multistart_files/figure-html/unnamed-chunk-2-1.png
Binary files 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
--- /dev/null
+++ b/vignettes/web_only/multistart_files/figure-html/unnamed-chunk-3-1.png
Binary files differ

Contact - Imprint