diff options
Diffstat (limited to 'R/mhmkin.R')
| -rw-r--r-- | R/mhmkin.R | 192 |
1 files changed, 192 insertions, 0 deletions
diff --git a/R/mhmkin.R b/R/mhmkin.R new file mode 100644 index 00000000..7f3ff9fa --- /dev/null +++ b/R/mhmkin.R @@ -0,0 +1,192 @@ +#' 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(status(x)) +} + +#' @export +AIC.mhmkin <- function(object, ..., k = 2) { + if (inherits(object[[1]], "saem.mmkin")) { + check_failed <- function(x) if (inherits(x$so, "try-error")) TRUE else FALSE + } + res <- sapply(object, function(x) { + if (check_failed(x)) return(NA) + else return(AIC(x$so, k = k)) + }) + dim(res) <- dim(object) + dimnames(res) <- dimnames(object) + return(res) +} + +#' @export +BIC.mhmkin <- function(object, ...) { + if (inherits(object[[1]], "saem.mmkin")) { + check_failed <- function(x) if (inherits(x$so, "try-error")) TRUE else FALSE + } + res <- sapply(object, function(x) { + if (check_failed(x)) return(NA) + else return(BIC(x$so)) + }) + dim(res) <- dim(object) + dimnames(res) <- dimnames(object) + return(res) +} + +#' @export +update.mhmkin <- function(object, ..., evaluate = TRUE) { + call <- attr(object, "call") + # For some reason we get mhkin.list in call[[1]] when using mhmkin from the + # loaded package so we need to fix this so we do not have to export + # mhmkin.list in addition to the S3 method mhmkin + call[[1]] <- mhmkin + + 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 +} |