diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/illparms.R | 70 | ||||
| -rw-r--r-- | R/mhmkin.R | 184 | ||||
| -rw-r--r-- | R/saem.R | 34 |
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) +} @@ -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 +} |