aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
Diffstat (limited to 'R')
-rw-r--r--R/illparms.R70
-rw-r--r--R/mhmkin.R184
-rw-r--r--R/saem.R34
3 files changed, 285 insertions, 3 deletions
diff --git a/R/illparms.R b/R/illparms.R
index f23f1cae..e4b28c56 100644
--- a/R/illparms.R
+++ b/R/illparms.R
@@ -1,12 +1,29 @@
#' Method to get the names of ill-defined parameters
#'
+#' The method for generalised nonlinear regression fits as obtained
+#' with [mkinfit] and [mmkin] checks if the degradation parameters
+#' pass the Wald test (in degradation kinetics often simply called t-test) for
+#' significant difference from zero. For this test, the parameterisation
+#' without parameter transformations is used.
+#'
+#' The method for hierarchical model fits, also known as nonlinear
+#' mixed-effects model fits as obtained with [saem] and [mhmkin]
+#' checks if any of the confidence intervals for the random
+#' effects expressed as standard deviations include zero, and if
+#' the confidence intervals for the error model parameters include
+#' zero.
+#'
#' @param object The object to investigate
#' @param x The object to be printed
#' @param conf.level The confidence level for checking p values
#' @param \dots For potential future extensions
-#' @return For [mkinfit] objects, a character vector of parameter names
-#' For [mmkin] objects, an object of class 'illparms.mmkin' with a
-#' suitable printing method.
+#' @param random For hierarchical fits, should random effects be tested?
+#' @param errmod For hierarchical fits, should error model parameters be
+#' tested?
+#' @return For [mkinfit] or [saem] objects, a character vector of parameter
+#' names. For [mmkin] or [mhmkin] objects, a matrix like object of class
+#' 'illparms.mmkin' or 'illparms.mhmkin'. The latter objects have a suitable
+#' printing method.
#' @export
illparms <- function(object, ...)
{
@@ -60,3 +77,50 @@ print.illparms.mmkin <- function(x, ...) {
class(x) <- NULL
print(x, quote = FALSE)
}
+
+#' @rdname illparms
+#' @export
+illparms.saem.mmkin <- function(object, conf.level = 0.95, random = TRUE, errmod = TRUE, ...) {
+ if (inherits(object, "try-error")) {
+ return(NA)
+ } else {
+ ints <- intervals(object, conf.level = conf.level)
+ failed <- NULL
+ if (random) {
+ failed_random <- ints$random[, "lower"] < 0
+ failed <- c(failed, names(which(failed_random)))
+ }
+ if (errmod) {
+ failed_errmod <- ints$errmod[, "lower"] < 0 & ints$errmod[, "est."] > 0
+ failed <- c(failed, names(which(failed_errmod)))
+ }
+ }
+ return(failed)
+}
+
+#' @rdname illparms
+#' @export
+illparms.mhmkin <- function(object, conf.level = 0.95, random = TRUE, errmod = TRUE, ...) {
+ result <- lapply(object,
+ function(fit) {
+ if (inherits(fit, "try-error")) return("E")
+ ill <- illparms(fit, conf.level = conf.level, random = random, errmod = errmod)
+ if (length(ill) > 0) {
+ return(paste(ill, collapse = ", "))
+ } else {
+ return("")
+ }
+ })
+ result <- unlist(result)
+ dim(result) <- dim(object)
+ dimnames(result) <- dimnames(object)
+ class(result) <- "illparms.mhmkin"
+ return(result)
+}
+
+#' @rdname illparms
+#' @export
+print.illparms.mhmkin <- function(x, ...) {
+ class(x) <- NULL
+ print(x, quote = FALSE)
+}
diff --git a/R/mhmkin.R b/R/mhmkin.R
new file mode 100644
index 00000000..61d67204
--- /dev/null
+++ b/R/mhmkin.R
@@ -0,0 +1,184 @@
+#' Fit nonlinear mixed-effects models built from one or more kinetic
+#' degradation models and one or more error models
+#'
+#' The name of the methods expresses that (**m**ultiple) **h**ierarchichal
+#' (also known as multilevel) **m**ulticompartment **kin**etic models are
+#' fitted. Our kinetic models are nonlinear, so we can use various nonlinear
+#' mixed-effects model fitting functions.
+#'
+#' @param objects A list of [mmkin] objects containing fits of the same
+#' degradation models to the same data, but using different error models.
+#' @param backend The backend to be used for fitting. Currently, only saemix is
+#' supported
+#' @param algorithm The algorithm to be used for fitting (currently not used)
+#' @param \dots Further arguments that will be passed to the nonlinear mixed-effects
+#' model fitting function.
+#' @param cores The number of cores to be used for multicore processing. This
+#' is only used when the \code{cluster} argument is \code{NULL}. On Windows
+#' machines, cores > 1 is not supported, you need to use the \code{cluster}
+#' argument to use multiple logical processors. Per default, all cores detected
+#' by [parallel::detectCores()] are used, except on Windows where the default
+#' is 1.
+#' @param cluster A cluster as returned by [makeCluster] to be used for
+#' parallel execution.
+#' @importFrom parallel mclapply parLapply detectCores
+#' @return A two-dimensional [array] of fit objects and/or try-errors that can
+#' be indexed using the degradation model names for the first index (row index)
+#' and the error model names for the second index (column index), with class
+#' attribute 'mhmkin'.
+#' @author Johannes Ranke
+#' @seealso \code{\link{[.mhmkin}} for subsetting [mhmkin] objects
+#' @export
+mhmkin <- function(objects, backend = "saemix", algorithm = "saem", ...) {
+ UseMethod("mhmkin")
+}
+
+#' @export
+#' @rdname mhmkin
+mhmkin.list <- function(objects, backend = "saemix",
+ ...,
+ cores = if (Sys.info()["sysname"] == "Windows") 1 else parallel::detectCores(), cluster = NULL)
+{
+ call <- match.call()
+ dot_args <- list(...)
+ backend_function <- switch(backend,
+ saemix = "saem"
+ )
+
+ deg_models <- lapply(objects[[1]][, 1], function(x) x$mkinmod)
+ names(deg_models) <- dimnames(objects[[1]])$model
+ n.deg <- length(deg_models)
+
+ ds <- lapply(objects[[1]][1, ], function(x) x$data)
+
+ for (other in objects[-1]) {
+ # Check if the degradation models in all objects are the same
+ for (deg_model_name in names(deg_models)) {
+ if (!all.equal(other[[deg_model_name, 1]]$mkinmod$spec,
+ deg_models[[deg_model_name]]$spec))
+ {
+ stop("The mmkin objects have to be based on the same degradation models")
+ }
+ }
+ # Check if they have been fitted to the same dataset
+ other_object_ds <- lapply(other[1, ], function(x) x$data)
+ for (i in 1:length(ds)) {
+ if (!all.equal(ds[[i]][c("time", "variable", "observed")],
+ other_object_ds[[i]][c("time", "variable", "observed")]))
+ {
+ stop("The mmkin objects have to be fitted to the same datasets")
+ }
+ }
+ }
+
+ n.o <- length(objects)
+
+ error_models = sapply(objects, function(x) x[[1]]$err_mod)
+ n.e <- length(error_models)
+
+ n.fits <- n.deg * n.e
+ fit_indices <- matrix(1:n.fits, ncol = n.e)
+ dimnames(fit_indices) <- list(degradation = names(deg_models),
+ error = error_models)
+
+ fit_function <- function(fit_index) {
+ w <- which(fit_indices == fit_index, arr.ind = TRUE)
+ deg_index <- w[1]
+ error_index <- w[2]
+ mmkin_row <- objects[[error_index]][deg_index, ]
+ res <- try(do.call(backend_function, args = c(list(mmkin_row), dot_args)))
+ return(res)
+ }
+
+ fit_time <- system.time({
+ if (is.null(cluster)) {
+ results <- parallel::mclapply(as.list(1:n.fits), fit_function,
+ mc.cores = cores, mc.preschedule = FALSE)
+ } else {
+ results <- parallel::parLapply(cluster, as.list(1:n.fits), fit_function)
+ }
+ })
+
+ attributes(results) <- attributes(fit_indices)
+ attr(results, "call") <- call
+ attr(results, "time") <- fit_time
+ class(results) <- "mhmkin"
+ return(results)
+}
+
+#' Subsetting method for mhmkin objects
+#'
+#' @param x An [mhmkin] object.
+#' @param i Row index selecting the fits for specific models
+#' @param j Column index selecting the fits to specific datasets
+#' @param drop If FALSE, the method always returns an mhmkin object, otherwise
+#' either a list of fit objects or a single fit object.
+#' @return An object of class \code{\link{mhmkin}}.
+#' @rdname mhmkin
+#' @export
+`[.mhmkin` <- function(x, i, j, ..., drop = FALSE) {
+ class(x) <- NULL
+ x_sub <- x[i, j, drop = drop]
+
+ if (!drop) {
+ class(x_sub) <- "mhmkin"
+ }
+ return(x_sub)
+}
+
+#' Print method for mhmkin objects
+#'
+#' @rdname mhmkin
+#' @export
+print.mhmkin <- function(x, ...) {
+ cat("<mhmkin> object\n")
+ cat("Status of individual fits:\n\n")
+ print(convergence(x))
+}
+
+#' @export
+convergence.mhmkin <- function(object, ...) {
+ all_summary_warnings <- character()
+
+ result <- lapply(object,
+ function(fit) {
+ if (inherits(fit, "try-error")) return("E")
+ else {
+ return("OK")
+ }
+ })
+ result <- unlist(result)
+ dim(result) <- dim(object)
+ dimnames(result) <- dimnames(object)
+
+ class(result) <- "convergence.mhmkin"
+ return(result)
+}
+
+#' @export
+print.convergence.mhmkin <- function(x, ...) {
+ class(x) <- NULL
+ print(x, quote = FALSE)
+ cat("\n")
+ if (any(x == "OK")) cat("OK: Fit terminated successfully\n")
+ if (any(x == "E")) cat("E: Error\n")
+}
+
+#' @export
+AIC.mhmkin <- function(object, ..., k = 2) {
+ res <- sapply(object, function(x) {
+ if (inherits(x, "try-error")) return(NA)
+ else return(AIC(x$so, k = k))
+ })
+ dim(res) <- dim(object)
+ dimnames(res) <- dimnames(object)
+ return(res)
+}
+
+#' @export
+BIC.mhmkin <- function(object, ...) {
+ res <- sapply(object, function(x) BIC(x$so))
+ dim(res) <- dim(object)
+ dimnames(res) <- dimnames(object)
+ return(res)
+}
diff --git a/R/saem.R b/R/saem.R
index e2b4d21c..370de3d8 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -58,6 +58,9 @@ utils::globalVariables(c("predicted", "std"))
#' f_saem_sfo <- saem(f_mmkin_parent["SFO", ])
#' f_saem_fomc <- saem(f_mmkin_parent["FOMC", ])
#' f_saem_dfop <- saem(f_mmkin_parent["DFOP", ])
+#' illparms(f_saem_dfop)
+#' update(f_saem_dfop, covariance.model = diag(c(1, 1, 1, 0)))
+#' AIC(f_saem_dfop)
#'
#' # The returned saem.mmkin object contains an SaemixObject, therefore we can use
#' # functions from saemix
@@ -125,6 +128,7 @@ saem.mmkin <- function(object,
fail_with_errors = TRUE,
verbose = FALSE, quiet = FALSE, ...)
{
+ call <- match.call()
transformations <- match.arg(transformations)
m_saemix <- saemix_model(object, verbose = verbose,
degparms_start = degparms_start,
@@ -184,6 +188,7 @@ saem.mmkin <- function(object,
transform_rates = object[[1]]$transform_rates,
transform_fractions = object[[1]]$transform_fractions,
so = f_saemix,
+ call = call,
time = fit_time,
mean_dp_start = attr(m_saemix, "mean_dp_start"),
bparms.optim = bparms_optim,
@@ -579,3 +584,32 @@ saemix_data <- function(object, verbose = FALSE, ...) {
...)
return(res)
}
+
+#' @export
+logLik.saem.mmkin <- function(object, ...) return(logLik(object$so))
+
+#' @export
+update.saem.mmkin <- function(object, ..., evaluate = TRUE) {
+ call <- object$call
+ # For some reason we get saem.mmkin in the call when using mhmkin
+ # so we need to fix this in order to avoid exporting saem.mmkin
+ # in addition to the S3 method
+ call[[1]] <- saem
+ update_arguments <- match.call(expand.dots = FALSE)$...
+
+ if (length(update_arguments) > 0) {
+ update_arguments_in_call <- !is.na(match(names(update_arguments), names(call)))
+ }
+
+ for (a in names(update_arguments)[update_arguments_in_call]) {
+ call[[a]] <- update_arguments[[a]]
+ }
+
+ update_arguments_not_in_call <- !update_arguments_in_call
+ if(any(update_arguments_not_in_call)) {
+ call <- c(as.list(call), update_arguments[update_arguments_not_in_call])
+ call <- as.call(call)
+ }
+ if(evaluate) eval(call, parent.frame())
+ else call
+}

Contact - Imprint