diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/anova.saem.mmkin.R | 106 | ||||
| -rw-r--r-- | R/aw.R | 36 | ||||
| -rw-r--r-- | R/illparms.R | 137 | ||||
| -rw-r--r-- | R/intervals.R | 5 | ||||
| -rw-r--r-- | R/llhist.R | 43 | ||||
| -rw-r--r-- | R/mhmkin.R | 192 | ||||
| -rw-r--r-- | R/mixed.mmkin.R | 3 | ||||
| -rw-r--r-- | R/mkinfit.R | 9 | ||||
| -rw-r--r-- | R/mmkin.R | 51 | ||||
| -rw-r--r-- | R/multistart.R | 224 | ||||
| -rw-r--r-- | R/nlme.R | 4 | ||||
| -rw-r--r-- | R/parms.R | 82 | ||||
| -rw-r--r-- | R/parms.mkinfit.R | 62 | ||||
| -rw-r--r-- | R/parplot.R | 105 | ||||
| -rw-r--r-- | R/plot.mixed.mmkin.R | 5 | ||||
| -rw-r--r-- | R/saem.R | 402 | ||||
| -rw-r--r-- | R/set_nd_nq.R | 164 | ||||
| -rw-r--r-- | R/status.R | 113 | ||||
| -rw-r--r-- | R/summary.mkinfit.R | 4 | ||||
| -rw-r--r-- | R/summary.mmkin.R | 56 | ||||
| -rw-r--r-- | R/summary.saem.mmkin.R | 27 |
21 files changed, 1617 insertions, 213 deletions
diff --git a/R/anova.saem.mmkin.R b/R/anova.saem.mmkin.R new file mode 100644 index 00000000..9937a919 --- /dev/null +++ b/R/anova.saem.mmkin.R @@ -0,0 +1,106 @@ +#' Anova method for saem.mmkin objects +#' +#' Generate an anova object. The method to calculate the BIC is that from +#' the saemix package. As in other prominent anova methods, models are sorted +## by likelihood, and the test is relative to the model on the previous line. +#' +#' @param object An [saem.mmkin] object +#' @param ... further such objects +#' @param method Method for likelihood calculation: "is" (importance sampling), +#' "lin" (linear approximation), or "gq" (Gaussian quadrature). Passed +#' to [saemix::logLik.SaemixObject] +#' @param test Should a likelihood ratio test be performed? If TRUE, +#' the alternative models are tested against the first model. Should +#' only be done for nested models. +#' @param model.names Optional character vector of model names +#' @importFrom stats anova logLik update pchisq terms +#' @importFrom methods is +#' @importFrom utils capture.output +#' @export +#' @return an "anova" data frame; the traditional (S3) result of anova() +anova.saem.mmkin <- function(object, ..., + method = c("is", "lin", "gq"), test = FALSE, model.names = NULL) +{ + # The following code is heavily inspired by anova.lmer in the lme4 package + mCall <- match.call(expand.dots = TRUE) + dots <- list(...) + method <- match.arg(method) + + is_model <- sapply(dots, is, "saem.mmkin") + if (any(is_model)) { + mods <- c(list(object), dots[is_model]) + + # Ensure same data, ignoring covariates + same_data <- sapply(dots[is_model], function(x) { + identical(object$data[c("ds", "name", "time", "value")], + x$data[c("ds", "name", "time", "value")]) + }) + if (!all(same_data)) { + stop(sum(!same_data), " objects have not been fitted to the same data") + } + + if (is.null(model.names)) + model.names <- vapply(as.list(mCall)[c(FALSE, TRUE, is_model)], deparse1, "") + + # Sanitize model names + if (length(model.names) != length(mods)) + stop("Model names vector and model list have different lengths") + + if (any(duplicated(model.names))) + stop("Duplicate model names are not allowed") + + if (max(nchar(model.names)) > 200) { + warning("Model names longer than 200 characters, assigning generic names") + model.names <- paste0("MODEL",seq_along(model.names)) + } + names(mods) <- model.names + + llks <- lapply(model.names, function(x) { + llk <- try(logLik(mods[[x]], method = method)) + if (inherits(llk, "try-error")) + stop("Could not obtain log likelihood with '", method, "' method for ", x) + return(llk) + }) + + # Order models by increasing degrees of freedom: + npar <- vapply(llks, attr, FUN.VALUE=numeric(1), "df") + ii <- order(npar) + mods <- mods[ii] + llks <- llks[ii] + npar <- npar[ii] + + # Describe data for the header, as in summary.saem.mmkin + header <- paste("Data:", nrow(object$data), "observations of", + length(unique(object$data$name)), "variable(s) grouped in", + length(unique(object$data$ds)), "datasets\n") + + # Calculate statistics + llk <- unlist(llks) + chisq <- 2 * pmax(0, c(NA, diff(llk))) + dfChisq <- c(NA, diff(npar)) + + bic <- function(x, method) { + BIC(x$so, method = method) + } + val <- data.frame( + npar = npar, + AIC = sapply(llks, AIC), + BIC = sapply(mods, bic, method = method), # We use the saemix method here + Lik = llk, + row.names = names(mods), check.names = FALSE) + + if (test) { + testval <- data.frame( + Chisq = chisq, + Df = dfChisq, + "Pr(>Chisq)" = ifelse(dfChisq == 0, NA, + pchisq(chisq, dfChisq, lower.tail = FALSE)), + row.names = names(mods), check.names = FALSE) + val <- cbind(val, testval) + } + class(val) <- c("anova", class(val)) + structure(val, heading = c(header)) + } else { + stop("Currently, no anova method is implemented for the case of a single saem.mmkin object") + } +} @@ -30,6 +30,14 @@ #' @export aw <- function(object, ...) UseMethod("aw") +.aw <- function(all_objects) { + AIC_all <- sapply(all_objects, AIC) + delta_i <- AIC_all - min(AIC_all) + denom <- sum(exp(-delta_i/2)) + w_i <- exp(-delta_i/2) / denom + return(w_i) +} + #' @export #' @rdname aw aw.mkinfit <- function(object, ...) { @@ -43,11 +51,7 @@ aw.mkinfit <- function(object, ...) { } } all_objects <- list(object, ...) - AIC_all <- sapply(all_objects, AIC) - delta_i <- AIC_all - min(AIC_all) - denom <- sum(exp(-delta_i/2)) - w_i <- exp(-delta_i/2) / denom - return(w_i) + .aw(all_objects) } #' @export @@ -56,3 +60,25 @@ aw.mmkin <- function(object, ...) { if (ncol(object) > 1) stop("Please supply an mmkin column object") do.call(aw, object) } + +#' @export +#' @rdname aw +aw.mixed.mmkin <- function(object, ...) { + oo <- list(...) + data_object <- object$data[c("ds", "name", "time", "value")] + for (i in seq_along(oo)) { + if (!inherits(oo[[i]], "mixed.mmkin")) stop("Please supply objects inheriting from mixed.mmkin") + data_other_object <- oo[[i]]$data[c("ds", "name", "time", "value")] + if (!identical(data_object, data_other_object)) { + stop("It seems that the mixed.mmkin objects have not all been fitted to the same data") + } + } + all_objects <- list(object, ...) + .aw(all_objects) +} + +#' @export +#' @rdname aw +aw.multistart <- function(object, ...) { + do.call(aw, object) +} diff --git a/R/illparms.R b/R/illparms.R new file mode 100644 index 00000000..931d8f05 --- /dev/null +++ b/R/illparms.R @@ -0,0 +1,137 @@ +#' 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 +#' @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, ...) +{ + UseMethod("illparms", object) +} + +#' @rdname illparms +#' @export +#' @examples +#' fit <- mkinfit("FOMC", FOCUS_2006_A, quiet = TRUE) +#' illparms(fit) +illparms.mkinfit <- function(object, conf.level = 0.95, ...) { + p_values <- suppressWarnings(summary(object)$bpar[, "Pr(>t)"]) + na <- is.na(p_values) + failed <- p_values > 1 - conf.level + names(parms(object))[na | failed] +} + +#' @rdname illparms +#' @export +#' @examples +#' \dontrun{ +#' fits <- mmkin( +#' c("SFO", "FOMC"), +#' list("FOCUS A" = FOCUS_2006_A, +#' "FOCUS C" = FOCUS_2006_C), +#' quiet = TRUE) +#' illparms(fits) +#' } +illparms.mmkin <- function(object, conf.level = 0.95, ...) { + result <- lapply(object, + function(fit) { + if (inherits(fit, "try-error")) return("E") + ill <- illparms(fit, conf.level = conf.level) + if (length(ill) > 0) { + return(paste(ill, collapse = ", ")) + } else { + return("") + } + }) + result <- unlist(result) + dim(result) <- dim(object) + dimnames(result) <- dimnames(object) + class(result) <- "illparms.mmkin" + return(result) +} + +#' @rdname illparms +#' @export +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, ...) { + if (inherits(object[[1]], "saem.mmkin")) { + check_failed <- function(x) if (inherits(x$so, "try-error")) TRUE else FALSE + } + result <- lapply(object, + function(fit) { + if (check_failed(fit)) { + return("E") + } else { + if (!is.null(fit$FIM_failed) && + "random effects and error model parameters" %in% fit$FIM_failed) { + return("NA") + } else { + 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/intervals.R b/R/intervals.R index 258eb4ad..705ef6eb 100644 --- a/R/intervals.R +++ b/R/intervals.R @@ -77,8 +77,9 @@ intervals.saem.mmkin <- function(object, level = 0.95, backtransform = TRUE, ... attr(confint_ret, "label") <- "Fixed effects:" # Random effects - ranef_ret <- as.matrix(conf.int[paste0("SD.", pnames), c("lower", "est.", "upper")]) - rownames(ranef_ret) <- paste0(gsub("SD\\.", "sd(", rownames(ranef_ret)), ")") + sdnames <- intersect(rownames(conf.int), paste("SD", pnames, sep = ".")) + ranef_ret <- as.matrix(conf.int[sdnames, c("lower", "est.", "upper")]) + rownames(ranef_ret) <- paste0(gsub("SD\\.", "sd(", sdnames), ")") attr(ranef_ret, "label") <- "Random effects:" diff --git a/R/llhist.R b/R/llhist.R new file mode 100644 index 00000000..e158495d --- /dev/null +++ b/R/llhist.R @@ -0,0 +1,43 @@ +#' Plot the distribution of log likelihoods from multistart objects +#' +#' Produces a histogram of log-likelihoods. In addition, the likelihood of the +#' original fit is shown as a red vertical line. +#' +#' @param object The [multistart] object +#' @param breaks Passed to [hist] +#' @param lpos Positioning of the legend. +#' @param main Title of the plot +#' @param \dots Passed to [hist] +#' @seealso [multistart] +#' @export +llhist <- function(object, breaks = "Sturges", lpos = "topleft", main = "", + ...) +{ + oldpar <- par(no.readonly = TRUE) + on.exit(par(oldpar, no.readonly = TRUE)) + + 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_orig <- logLik(attr(object, "orig")) + ll <- stats::na.omit(sapply(object, llfunc)) + + par(las = 1) + h <- hist(ll, freq = TRUE, + xlab = "", main = main, + ylab = "Frequency of log likelihoods", breaks = breaks, ...) + + freq_factor <- h$counts[1] / h$density[1] + + abline(v = ll_orig, col = 2) + + legend(lpos, inset = c(0.05, 0.05), bty = "n", + lty = 1, col = c(2), + legend = "original fit") +} 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 +} diff --git a/R/mixed.mmkin.R b/R/mixed.mmkin.R index 682a9a34..fd6d975f 100644 --- a/R/mixed.mmkin.R +++ b/R/mixed.mmkin.R @@ -1,5 +1,6 @@ #' Create a mixed effects model from an mmkin row object #' +#' @importFrom rlang !!! #' @param object An [mmkin] row object #' @param method The method to be used #' @param \dots Currently not used @@ -65,7 +66,7 @@ mixed.mmkin <- function(object, method = c("none"), ...) { function(x) x$data[c("variable", "time", "observed", "predicted", "residual")]) names(ds_list) <- ds_names - res$data <- purrr::map_dfr(ds_list, function(x) x, .id = "ds") + res$data <- vctrs::vec_rbind(!!!ds_list, .names_to = "ds") names(res$data)[1:4] <- c("ds", "name", "time", "value") res$data$name <- as.character(res$data$name) res$data$ds <- ordered(res$data$ds, levels = unique(res$data$ds)) diff --git a/R/mkinfit.R b/R/mkinfit.R index 1c2f8b5e..693778fd 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -844,8 +844,13 @@ mkinfit <- function(mkinmod, observed, } } + dim_hessian <- length(c(degparms, errparms)) + fit$hessian <- try(numDeriv::hessian(cost_function, c(degparms, errparms), OLS = FALSE, update_data = FALSE), silent = TRUE) + if (inherits(fit$hessian, "try-error")) { + fit$hessian <- matrix(NA, nrow = dim_hessian, ncol = dim_hessian) + } dimnames(fit$hessian) <- list(names(c(degparms, errparms)), names(c(degparms, errparms))) @@ -858,7 +863,9 @@ mkinfit <- function(mkinmod, observed, fit$hessian_notrans <- try(numDeriv::hessian(cost_function, c(bparms.optim, errparms), OLS = FALSE, trans = FALSE, update_data = FALSE), silent = TRUE) - + if (inherits(fit$hessian_notrans, "try-error")) { + fit$hessian_notrans <- matrix(NA, nrow = dim_hessian, ncol = dim_hessian) + } dimnames(fit$hessian_notrans) <- list(names(c(bparms.optim, errparms)), names(c(bparms.optim, errparms))) }) @@ -114,15 +114,18 @@ mmkin <- function(models = c("SFO", "FOMC", "DFOP"), datasets, return(res) } - 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) - } + 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) <- "mmkin" return(results) } @@ -168,41 +171,7 @@ mmkin <- function(models = c("SFO", "FOMC", "DFOP"), datasets, print.mmkin <- function(x, ...) { cat("<mmkin> object\n") cat("Status of individual fits:\n\n") - all_summary_warnings <- character() - sww <- 0 # Counter for Shapiro-Wilks warnings - - display <- lapply(x, - function(fit) { - if (inherits(fit, "try-error")) return("E") - sw <- fit$summary_warnings - swn <- names(sw) - if (length(sw) > 0) { - if (any(grepl("S", swn))) { - sww <<- sww + 1 - swn <- gsub("S", paste0("S", sww), swn) - } - warnstring <- paste(swn, collapse = ", ") - names(sw) <- swn - all_summary_warnings <<- c(all_summary_warnings, sw) - return(warnstring) - } else { - return("OK") - } - }) - display <- unlist(display) - dim(display) <- dim(x) - dimnames(display) <- dimnames(x) - print(display, quote = FALSE) - - cat("\n") - if (any(display == "OK")) cat("OK: No warnings\n") - if (any(display == "E")) cat("E: Error\n") - u_swn <- unique(names(all_summary_warnings)) - u_w <- all_summary_warnings[u_swn] - for (i in seq_along(u_w)) { - cat(names(u_w)[i], ": ", u_w[i], "\n", sep = "") - } - + print(status(x)) } #' @export diff --git a/R/multistart.R b/R/multistart.R new file mode 100644 index 00000000..61ef43dc --- /dev/null +++ b/R/multistart.R @@ -0,0 +1,224 @@ +#' Perform a hierarchical model fit with multiple starting values +#' +#' The purpose of this method is to check if a certain algorithm for fitting +#' nonlinear hierarchical models (also known as nonlinear mixed-effects models) +#' will reliably yield results that are sufficiently similar to each other, if +#' started with a certain range of reasonable starting parameters. It is +#' inspired by the article on practical identifiabiliy in the frame of nonlinear +#' mixed-effects models by Duchesne et al (2021). +#' +#' @param object The fit object to work with +#' @param n How many different combinations of starting parameters should be +#' used? +#' @param cores How many fits should be run in parallel (only on posix platforms)? +#' @param cluster A cluster as returned by [parallel::makeCluster] to be used +#' for parallel execution. +#' @param \dots Passed to the update function. +#' @param x The multistart object to print +#' @return A list of [saem.mmkin] objects, with class attributes +#' 'multistart.saem.mmkin' and 'multistart'. +#' @seealso [parplot], [llhist] +#' +#' @references Duchesne R, Guillemin A, Gandrillon O, Crauste F. Practical +#' identifiability in the frame of nonlinear mixed effects models: the example +#' of the in vitro erythropoiesis. BMC Bioinformatics. 2021 Oct 4;22(1):478. +#' doi: 10.1186/s12859-021-04373-4. +#' @export +#' @examples +#' \dontrun{ +#' 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) +#' parplot(f_saem_full_multi, lpos = "topleft") +#' illparms(f_saem_full) +#' +#' f_saem_reduced <- update(f_saem_full, no_random_effect = "log_k2") +#' illparms(f_saem_reduced) +#' # On Windows, we need to create a cluster first. When working with +#' # such a cluster, we need to export the mmkin object to the cluster +#' # nodes, as it is referred to when updating the saem object on the nodes. +#' library(parallel) +#' cl <- makePSOCKcluster(12) +#' f_saem_reduced_multi <- multistart(f_saem_reduced, n = 16, cluster = cl) +#' parplot(f_saem_reduced_multi, lpos = "topright") +#' stopCluster(cl) +#' } +multistart <- function(object, n = 50, + cores = if (Sys.info()["sysname"] == "Windows") 1 else parallel::detectCores(), + cluster = NULL, ...) +{ + UseMethod("multistart", object) +} + +#' @rdname multistart +#' @export +multistart.saem.mmkin <- function(object, n = 50, cores = 1, + cluster = NULL, ...) { + call <- match.call() + if (n <= 1) stop("Please specify an n of at least 2") + + mmkin_object <- object$mmkin + + mmkin_parms <- parms(mmkin_object, errparms = FALSE, + transformed = object$transformations == "mkin") + start_parms <- apply( + mmkin_parms, 1, + function(x) stats::runif(n, min(x), max(x))) + + saem_call <- object$call + saem_call[[1]] <- saem + saem_call[[2]] <- mmkin_object + i_startparms <- which(names(saem_call) == "degparms_start") + + fit_function <- function(x) { + + new_startparms <- str2lang( + paste0(capture.output(dput(start_parms[x, ])), + collapse = "")) + + if (length(i_startparms) == 0) { + saem_call <- c(as.list(saem_call), degparms_start = new_startparms) + saem_call <- as.call(saem_call) + } else { + saem_call[i_startparms] <- new_startparms + } + + ret <- eval(saem_call) + + return(ret) + } + + if (is.null(cluster)) { + res <- parallel::mclapply(1:n, fit_function, mc.cores = cores) + } else { + res <- parallel::parLapply(cluster, 1:n, fit_function) + } + attr(res, "orig") <- object + attr(res, "start_parms") <- start_parms + attr(res, "call") <- call + class(res) <- c("multistart.saem.mmkin", "multistart") + return(res) +} + +#' @export +status.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) <- "status.multistart" + return(result) +} + +#' @export +status.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) <- "status.multistart" + return(result) +} + +#' @export +print.status.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 +print.multistart <- function(x, ...) { + cat("<multistart> object with", length(x), "fits:\n") + print(status(x)) +} + +#' @rdname multistart +#' @export +best <- function(object, ...) +{ + UseMethod("best", object) +} + +#' @export +#' @return The object with the highest likelihood +#' @rdname multistart +best.default <- function(object, ...) +{ + return(object[[which.best(object)]]) +} + +#' @return The index of the object with the highest likelihood +#' @rdname multistart +#' @export +which.best <- function(object, ...) +{ + UseMethod("which.best", object) +} + +#' @rdname multistart +#' @export +which.best.default <- function(object, ...) +{ + llfunc <- function(object) { + ret <- try(logLik(object)) + if (inherits(ret, "try-error")) return(NA) + else return(ret) + } + ll <- sapply(object, llfunc) + return(which.max(ll)) +} + +#' @export +update.multistart <- function(object, ..., evaluate = TRUE) { + call <- attr(object, "call") + # For some reason we get multistart.saem.mmkin in call[[1]] when using multistart + # from the loaded package so we need to fix this so we do not have to export + # multistart.saem.mmkin + call[[1]] <- multistart + + 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 +} @@ -125,7 +125,7 @@ nlme_function <- function(object) { } #' @rdname nlme -#' @importFrom purrr map_dfr +#' @importFrom rlang !!! #' @return A \code{\link{groupedData}} object #' @export nlme_data <- function(object) { @@ -134,7 +134,7 @@ nlme_data <- function(object) { ds_list <- lapply(object, function(x) x$data[c("time", "variable", "observed")]) names(ds_list) <- ds_names - ds_nlme <- purrr::map_dfr(ds_list, function(x) x, .id = "ds") + ds_nlme <- vctrs::vec_rbind(!!!ds_list, .names_to = "ds") ds_nlme$variable <- as.character(ds_nlme$variable) ds_nlme$ds <- ordered(ds_nlme$ds, levels = unique(ds_nlme$ds)) ds_nlme_renamed <- data.frame(ds = ds_nlme$ds, name = ds_nlme$variable, diff --git a/R/parms.R b/R/parms.R new file mode 100644 index 00000000..bd4e479b --- /dev/null +++ b/R/parms.R @@ -0,0 +1,82 @@ +#' Extract model parameters +#' +#' This function returns degradation model parameters as well as error +#' model parameters per default, in order to avoid working with a fitted model +#' without considering the error structure that was assumed for the fit. +#' +#' @param object A fitted model object. +#' @param \dots Not used +#' @return Depending on the object, a numeric vector of fitted model parameters, +#' a matrix (e.g. for mmkin row objects), or a list of matrices (e.g. for +#' mmkin objects with more than one row). +#' @seealso [saem], [multistart] +#' @examples +#' # mkinfit objects +#' fit <- mkinfit("SFO", FOCUS_2006_C, quiet = TRUE) +#' parms(fit) +#' parms(fit, transformed = TRUE) +#' +#' # mmkin objects +#' ds <- lapply(experimental_data_for_UBA_2019[6:10], +#' function(x) subset(x$data[c("name", "time", "value")])) +#' names(ds) <- paste("Dataset", 6:10) +#' \dontrun{ +#' fits <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE, cores = 1) +#' parms(fits["SFO", ]) +#' parms(fits[, 2]) +#' parms(fits) +#' parms(fits, transformed = TRUE) +#' } +#' @export +parms <- function(object, ...) +{ + UseMethod("parms", object) +} + +#' @param transformed Should the parameters be returned as used internally +#' during the optimisation? +#' @param errparms Should the error model parameters be returned +#' in addition to the degradation parameters? +#' @rdname parms +#' @export +parms.mkinfit <- function(object, transformed = FALSE, errparms = TRUE, ...) +{ + res <- if (transformed) object$par + else c(object$bparms.optim, object$errparms) + if (!errparms) { + res[setdiff(names(res), names(object$errparms))] + } + else return(res) +} + +#' @rdname parms +#' @export +parms.mmkin <- function(object, transformed = FALSE, errparms = TRUE, ...) +{ + if (nrow(object) == 1) { + res <- sapply(object, parms, transformed = transformed, + errparms = errparms, ...) + colnames(res) <- colnames(object) + } else { + res <- list() + for (i in 1:nrow(object)) { + res[[i]] <- parms(object[i, ], transformed = transformed, + errparms = errparms, ...) + } + names(res) <- rownames(object) + } + 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/parms.mkinfit.R b/R/parms.mkinfit.R deleted file mode 100644 index a1f2d209..00000000 --- a/R/parms.mkinfit.R +++ /dev/null @@ -1,62 +0,0 @@ -#' Extract model parameters from mkinfit models -#' -#' This function always returns degradation model parameters as well as error -#' model parameters, in order to avoid working with a fitted model without -#' considering the error structure that was assumed for the fit. -#' -#' @param object A fitted model object. Methods are implemented for -#' [mkinfit()] objects and for [mmkin()] objects. -#' @param \dots Not used -#' @return For mkinfit objects, a numeric vector of fitted model parameters. -#' For mmkin row objects, a matrix with the parameters with a -#' row for each dataset. If the mmkin object has more than one row, a list of -#' such matrices is returned. -#' @examples -#' # mkinfit objects -#' fit <- mkinfit("SFO", FOCUS_2006_C, quiet = TRUE) -#' parms(fit) -#' parms(fit, transformed = TRUE) -#' -#' # mmkin objects -#' ds <- lapply(experimental_data_for_UBA_2019[6:10], -#' function(x) subset(x$data[c("name", "time", "value")])) -#' names(ds) <- paste("Dataset", 6:10) -#' \dontrun{ -#' fits <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE, cores = 1) -#' parms(fits["SFO", ]) -#' parms(fits[, 2]) -#' parms(fits) -#' parms(fits, transformed = TRUE) -#' } -#' @export -parms <- function(object, ...) -{ - UseMethod("parms", object) -} - -#' @param transformed Should the parameters be returned -#' as used internally during the optimisation? -#' @rdname parms -#' @export -parms.mkinfit <- function(object, transformed = FALSE, ...) -{ - if (transformed) object$par - else c(object$bparms.optim, object$errparms) -} - -#' @rdname parms -#' @export -parms.mmkin <- function(object, transformed = FALSE, ...) -{ - if (nrow(object) == 1) { - res <- sapply(object, parms, transformed = transformed, ...) - colnames(res) <- colnames(object) - } else { - res <- list() - for (i in 1:nrow(object)) { - res[[i]] <- parms(object[i, ], transformed = transformed, ...) - } - names(res) <- rownames(object) - } - return(res) -} diff --git a/R/parplot.R b/R/parplot.R new file mode 100644 index 00000000..627a4eb8 --- /dev/null +++ b/R/parplot.R @@ -0,0 +1,105 @@ +#' Plot parameter variability of multistart objects +#' +#' Produces a boxplot with all parameters from the multiple runs, scaled +#' either by the parameters of the run with the highest likelihood, +#' or by their medians as proposed in the paper by Duchesne et al. (2021). +#' +#' @param object The [multistart] object +#' @param llmin The minimum likelihood of objects to be shown +#' @param scale By default, scale parameters using the best available fit. +#' If 'median', parameters are scaled using the median parameters from all fits. +#' @param main Title of the plot +#' @param lpos Positioning of the legend. +#' @param \dots Passed to [boxplot] +#' @references Duchesne R, Guillemin A, Gandrillon O, Crauste F. Practical +#' identifiability in the frame of nonlinear mixed effects models: the example +#' of the in vitro erythropoiesis. BMC Bioinformatics. 2021 Oct 4;22(1):478. +#' doi: 10.1186/s12859-021-04373-4. +#' @seealso [multistart] +#' @importFrom stats median +#' @export +parplot <- function(object, ...) { + UseMethod("parplot") +} + +#' @rdname parplot +#' @export +parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best", "median"), + lpos = "bottomleft", main = "", ...) +{ + oldpar <- par(no.readonly = TRUE) + on.exit(par(oldpar, no.readonly = TRUE)) + + orig <- attr(object, "orig") + orig_parms <- parms(orig) + start_parms <- orig$mean_dp_start + all_parms <- parms(object) + + if (inherits(object, "multistart.saem.mmkin")) { + llfunc <- function(object) { + if (inherits(object$so, "try-error")) return(NA) + else return(logLik(object$so)) + } + } else { + stop("parplot is only implemented for multistart.saem.mmkin objects") + } + ll <- sapply(object, llfunc) + selected <- which(ll > llmin) + selected_parms <- all_parms[selected, ] + + par(las = 1) + if (orig$transformations == "mkin") { + degparm_names_transformed <- names(start_parms) + degparm_index <- which(names(orig_parms) %in% degparm_names_transformed) + orig_parms[degparm_names_transformed] <- backtransform_odeparms( + orig_parms[degparm_names_transformed], + 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[[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])) + + selected_parms[, degparm_names_transformed] <- + t(apply(selected_parms[, degparm_names_transformed], 1, backtransform_odeparms, + orig$mmkin[[1]]$mkinmod, + transform_rates = orig$mmkin[[1]]$transform_rates, + transform_fractions = orig$mmkin[[1]]$transform_fractions)) + colnames(selected_parms)[1:length(degparm_names)] <- degparm_names + } + + scale <- match.arg(scale) + parm_scale <- switch(scale, + best = selected_parms[which.best(object[selected]), ], + median = apply(selected_parms, 2, median) + ) + + # Boxplots of all scaled parameters + selected_scaled_parms <- t(apply(selected_parms, 1, function(x) x / parm_scale)) + boxplot(selected_scaled_parms, log = "y", main = main, , + ylab = "Normalised parameters", ...) + + # Show starting parameters + start_scaled_parms <- rep(NA_real_, length(orig_parms)) + names(start_scaled_parms) <- names(orig_parms) + start_scaled_parms[names(start_parms)] <- + start_parms / parm_scale[names(start_parms)] + points(start_scaled_parms, col = 3, cex = 3) + + # Show parameters of original run + orig_scaled_parms <- orig_parms / parm_scale + points(orig_scaled_parms, col = 2, cex = 2) + + abline(h = 1, lty = 2) + + legend(lpos, inset = c(0.05, 0.05), bty = "n", + pch = 1, col = 3:1, lty = c(NA, NA, 1), + legend = c( + "Starting parameters", + "Original run", + "Multistart runs")) +} diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R index 3a444253..6815bfd2 100644 --- a/R/plot.mixed.mmkin.R +++ b/R/plot.mixed.mmkin.R @@ -153,7 +153,7 @@ plot.mixed.mmkin <- function(x, outtimes <- sort(unique(c(x$data$time, seq(xlim[1], xlim[2], length.out = 50)))) - pred_ds <- purrr::map_dfr(i, function(ds_i) { + pred_list <- lapply(i, function(ds_i) { odeparms_trans <- degparms_all[ds_i, odeparms_names] names(odeparms_trans) <- odeparms_names # needed if only one odeparm if (backtransform) { @@ -171,8 +171,9 @@ plot.mixed.mmkin <- function(x, out <- mkinpredict(x$mkinmod, odeparms, odeini, outtimes, solution_type = solution_type, atol = fit_1$atol, rtol = fit_1$rtol) - return(cbind(as.data.frame(out), ds = ds_names[ds_i])) }) + names(pred_list) <- ds_names[i] + pred_ds <- vctrs::vec_rbind(!!!pred_list, .names_to = "ds") odeparms_pop_trans <- degparms_all_pop[odeparms_names] @@ -15,35 +15,46 @@ utils::globalVariables(c("predicted", "std")) #' #' @importFrom utils packageVersion #' @param object An [mmkin] row object containing several fits of the same -#' [mkinmod] model to different datasets +#' [mkinmod] model to different datasets #' @param verbose Should we print information about created objects of -#' type [saemix::SaemixModel] and [saemix::SaemixData]? +#' type [saemix::SaemixModel] and [saemix::SaemixData]? #' @param transformations Per default, all parameter transformations are done -#' in mkin. If this argument is set to 'saemix', parameter transformations -#' are done in 'saemix' for the supported cases. Currently this is only -#' supported in cases where the initial concentration of the parent is not fixed, -#' SFO or DFOP is used for the parent and there is either no metabolite or one. +#' in mkin. If this argument is set to 'saemix', parameter transformations +#' are done in 'saemix' for the supported cases, i.e. (as of version 1.1.2) +#' SFO, FOMC, DFOP and HS without fixing `parent_0`, and SFO or DFOP with +#' one SFO metabolite. #' @param degparms_start Parameter values given as a named numeric vector will -#' be used to override the starting values obtained from the 'mmkin' object. +#' be used to override the starting values obtained from the 'mmkin' object. #' @param test_log_parms If TRUE, an attempt is made to use more robust starting -#' values for population parameters fitted as log parameters in mkin (like -#' rate constants) by only considering rate constants that pass the t-test -#' when calculating mean degradation parameters using [mean_degparms]. +#' values for population parameters fitted as log parameters in mkin (like +#' rate constants) by only considering rate constants that pass the t-test +#' when calculating mean degradation parameters using [mean_degparms]. #' @param conf.level Possibility to adjust the required confidence level -#' for parameter that are tested if requested by 'test_log_parms'. +#' for parameter that are tested if requested by 'test_log_parms'. #' @param solution_type Possibility to specify the solution type in case the -#' automatic choice is not desired +#' automatic choice is not desired +#' @param no_random_effect Character vector of degradation parameters for +#' which there should be no variability over the groups. Only used +#' if the covariance model is not explicitly specified. +#' @param covariance.model Will be passed to [saemix::SaemixModel()]. Per +#' default, uncorrelated random effects are specified for all degradation +#' parameters. +#' @param covariates A data frame with covariate data for use in +#' 'covariate_models', with dataset names as row names. +#' @param covariate_models A list containing linear model formulas with one explanatory +#' variable, i.e. of the type 'parameter ~ covariate'. Covariates must be available +#' in the 'covariates' data frame. #' @param fail_with_errors Should a failure to compute standard errors -#' from the inverse of the Fisher Information Matrix be a failure? +#' from the inverse of the Fisher Information Matrix be a failure? #' @param quiet Should we suppress the messages saemix prints at the beginning -#' and the end of the optimisation process? +#' and the end of the optimisation process? #' @param nbiter.saemix Convenience option to increase the number of -#' iterations +#' iterations #' @param control Passed to [saemix::saemix]. #' @param \dots Further parameters passed to [saemix::saemixModel]. #' @return An S3 object of class 'saem.mmkin', containing the fitted -#' [saemix::SaemixObject] as a list component named 'so'. The -#' object also inherits from 'mixed.mmkin'. +#' [saemix::SaemixObject] as a list component named 'so'. The +#' object also inherits from 'mixed.mmkin'. #' @seealso [summary.saem.mmkin] [plot.mixed.mmkin] #' @examples #' \dontrun{ @@ -58,7 +69,13 @@ 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", ]) +#' anova(f_saem_sfo, f_saem_fomc, f_saem_dfop) +#' anova(f_saem_sfo, f_saem_dfop, test = TRUE) +#' illparms(f_saem_dfop) +#' f_saem_dfop_red <- update(f_saem_dfop, no_random_effect = "g_qlogis") +#' anova(f_saem_dfop, f_saem_dfop_red, test = TRUE) #' +#' anova(f_saem_sfo, f_saem_fomc, f_saem_dfop) #' # The returned saem.mmkin object contains an SaemixObject, therefore we can use #' # functions from saemix #' library(saemix) @@ -70,7 +87,7 @@ utils::globalVariables(c("predicted", "std")) #' #' f_mmkin_parent_tc <- update(f_mmkin_parent, error_model = "tc") #' f_saem_fomc_tc <- saem(f_mmkin_parent_tc["FOMC", ]) -#' compare.saemix(f_saem_fomc$so, f_saem_fomc_tc$so) +#' anova(f_saem_fomc, f_saem_fomc_tc, test = TRUE) #' #' sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"), #' A1 = mkinsub("SFO")) @@ -118,6 +135,10 @@ saem.mmkin <- function(object, test_log_parms = TRUE, conf.level = 0.6, solution_type = "auto", + covariance.model = "auto", + covariates = NULL, + covariate_models = NULL, + no_random_effect = NULL, nbiter.saemix = c(300, 100), control = list(displayProgress = FALSE, print = FALSE, nbiter.saemix = nbiter.saemix, @@ -125,56 +146,68 @@ 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, test_log_parms = test_log_parms, conf.level = conf.level, solution_type = solution_type, - transformations = transformations, ...) - d_saemix <- saemix_data(object, verbose = verbose) + transformations = transformations, + covariance.model = covariance.model, + covariates = covariates, + covariate_models = covariate_models, + no_random_effect = no_random_effect, + ...) + d_saemix <- saemix_data(object, covariates = covariates, 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") + FIM_failed <- c(FIM_failed, "random effects and error model parameters") } 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, @@ -183,14 +216,14 @@ 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, + FIM_failed = FIM_failed, 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")), @@ -198,6 +231,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) } @@ -217,18 +256,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 <- x$so@results@conf.int[c("estimate", "lower", "upper")] - rownames(conf.int) <- x$so@results@conf.int[["name"]] - conf.int.var <- grepl("^Var\\.", rownames(conf.int)) - print(conf.int[!conf.int.var, ], 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) } @@ -236,14 +277,23 @@ print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) { #' @rdname saem #' @return An [saemix::SaemixModel] object. #' @export -saemix_model <- function(object, solution_type = "auto", transformations = c("mkin", "saemix"), - degparms_start = numeric(), test_log_parms = FALSE, conf.level = 0.6, verbose = FALSE, ...) +saemix_model <- function(object, solution_type = "auto", + transformations = c("mkin", "saemix"), degparms_start = numeric(), + covariance.model = "auto", no_random_effect = NULL, + covariates = NULL, covariate_models = NULL, + test_log_parms = FALSE, conf.level = 0.6, verbose = FALSE, ...) { if (nrow(object) > 1) stop("Only row objects allowed") mkin_model <- object[[1]]$mkinmod degparms_optim <- mean_degparms(object, test_log_parms = test_log_parms) + na_degparms <- names(which(is.na(degparms_optim))) + if (length(na_degparms) > 0) { + message("Did not find valid starting values for ", paste(na_degparms, collapse = ", "), "\n", + "Now trying with test_log_parms = FALSE") + degparms_optim <- mean_degparms(object, test_log_parms = FALSE) + } if (transformations == "saemix") { degparms_optim <- backtransform_odeparms(degparms_optim, object[[1]]$mkinmod, @@ -252,7 +302,8 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk } degparms_fixed <- object[[1]]$bparms.fixed - # Transformations are done in the degradation function + # Transformations are done in the degradation function by default + # (transformations = "mkin") transform.par = rep(0, length(degparms_optim)) odeini_optim_parm_names <- grep('_0$', names(degparms_optim), value = TRUE) @@ -275,7 +326,10 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk # Parent only if (length(mkin_model$spec) == 1) { parent_type <- mkin_model$spec[[1]]$type - if (length(odeini_fixed) == 1) { + if (length(odeini_fixed) == 1 && !grepl("_bound$", names(odeini_fixed))) { + if (transformations == "saemix") { + stop("saemix transformations are not supported for parent fits with fixed initial parent value") + } if (parent_type == "SFO") { stop("saemix needs at least two parameters to work on.") } @@ -303,6 +357,9 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk } } } else { + if (length(odeini_fixed) == 2) { + stop("SFORB with fixed initial parent value is not supported") + } if (parent_type == "SFO") { if (transformations == "mkin") { model_function <- function(psi, id, xidep) { @@ -316,8 +373,15 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk } } if (parent_type == "FOMC") { - model_function <- function(psi, id, xidep) { - psi[id, 1] / (xidep[, "time"]/exp(psi[id, 3]) + 1)^exp(psi[id, 2]) + if (transformations == "mkin") { + model_function <- function(psi, id, xidep) { + psi[id, 1] / (xidep[, "time"]/exp(psi[id, 3]) + 1)^exp(psi[id, 2]) + } + } else { + model_function <- function(psi, id, xidep) { + psi[id, 1] / (xidep[, "time"]/psi[id, 3] + 1)^psi[id, 2] + } + transform.par = c(0, 1, 1) } } if (parent_type == "DFOP") { @@ -338,14 +402,57 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk transform.par = c(0, 1, 1, 3) } } + if (parent_type == "SFORB") { + if (transformations == "mkin") { + model_function <- function(psi, id, xidep) { + k_12 <- exp(psi[id, 3]) + k_21 <- exp(psi[id, 4]) + k_1output <- exp(psi[id, 2]) + t <- xidep[, "time"] + + sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21) + b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp + b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp + + psi[id, 1] * (((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) + + ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t)) + } + } else { + model_function <- function(psi, id, xidep) { + k_12 <- psi[id, 3] + k_21 <- psi[id, 4] + k_1output <- psi[id, 2] + t <- xidep[, "time"] + + sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21) + b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp + b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp + + psi[id, 1] * (((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) + + ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t)) + } + transform.par = c(0, 1, 1, 1) + } + } if (parent_type == "HS") { - model_function <- function(psi, id, xidep) { - tb <- exp(psi[id, 4]) - t <- xidep[, "time"] - k1 = exp(psi[id, 2]) - psi[id, 1] * ifelse(t <= tb, - exp(- k1 * t), - exp(- k1 * tb) * exp(- exp(psi[id, 3]) * (t - tb))) + if (transformations == "mkin") { + model_function <- function(psi, id, xidep) { + tb <- exp(psi[id, 4]) + t <- xidep[, "time"] + k1 <- exp(psi[id, 2]) + psi[id, 1] * ifelse(t <= tb, + exp(- k1 * t), + exp(- k1 * tb) * exp(- exp(psi[id, 3]) * (t - tb))) + } + } else { + model_function <- function(psi, id, xidep) { + tb <- psi[id, 4] + t <- xidep[, "time"] + psi[id, 1] * ifelse(t <= tb, + exp(- psi[id, 2] * t), + exp(- psi[id, 2] * tb) * exp(- psi[id, 3] * (t - tb))) + } + transform.par = c(0, 1, 1, 1) } } } @@ -518,8 +625,52 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk degparms_psi0 <- degparms_optim degparms_psi0[names(degparms_start)] <- degparms_start - psi0_matrix <- matrix(degparms_psi0, nrow = 1) - colnames(psi0_matrix) <- names(degparms_psi0) + psi0_matrix <- matrix(degparms_psi0, nrow = 1, + dimnames = list("(Intercept)", names(degparms_psi0))) + + if (covariance.model[1] == "auto") { + covariance_diagonal <- rep(1, length(degparms_optim)) + if (!is.null(no_random_effect)) { + degparms_no_random <- which(names(degparms_psi0) %in% no_random_effect) + covariance_diagonal[degparms_no_random] <- 0 + } + covariance.model = diag(covariance_diagonal) + } + + if (is.null(covariate_models)) { + covariate.model <- matrix(nrow = 0, ncol = 0) # default in saemixModel() + } else { + degparms_dependent <- sapply(covariate_models, function(x) as.character(x[[2]])) + covariates_in_models = unique(unlist(lapply( + covariate_models, function(x) + colnames(attr(terms(x), "factors")) + ))) + covariates_not_available <- setdiff(covariates_in_models, names(covariates)) + if (length(covariates_not_available) > 0) { + stop("Covariate(s) ", paste(covariates_not_available, collapse = ", "), + " used in the covariate models not available in the covariate data") + } + psi0_matrix <- rbind(psi0_matrix, + matrix(0, nrow = length(covariates), ncol = ncol(psi0_matrix), + dimnames = list(names(covariates), colnames(psi0_matrix)))) + covariate.model <- matrix(0, nrow = length(covariates), + ncol = ncol(psi0_matrix), + dimnames = list( + covariates = names(covariates), + degparms = colnames(psi0_matrix))) + if (transformations == "saemix") { + stop("Covariate models with saemix transformations currently not supported") + } + parms_trans <- as.data.frame(t(sapply(object, parms, transformed = TRUE))) + for (covariate_model in covariate_models) { + covariate_name <- as.character(covariate_model[[2]]) + model_data <- cbind(parms_trans, covariates) + ini_model <- lm(covariate_model, data = model_data) + ini_coef <- coef(ini_model) + psi0_matrix[names(ini_coef), covariate_name] <- ini_coef + covariate.model[names(ini_coef)[-1], covariate_name] <- as.numeric(as.logical(ini_coef[-1])) + } + } res <- saemix::saemixModel(model_function, psi0 = psi0_matrix, @@ -527,6 +678,8 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk transform.par = transform.par, error.model = error.model, verbose = verbose, + covariance.model = covariance.model, + covariate.model = covariate.model, ... ) attr(res, "mean_dp_start") <- degparms_optim @@ -534,26 +687,107 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk } #' @rdname saem +#' @importFrom rlang !!! #' @return An [saemix::SaemixData] object. #' @export -saemix_data <- function(object, verbose = FALSE, ...) { +saemix_data <- function(object, covariates = NULL, verbose = FALSE, ...) { if (nrow(object) > 1) stop("Only row objects allowed") ds_names <- colnames(object) ds_list <- lapply(object, function(x) x$data[c("time", "variable", "observed")]) names(ds_list) <- ds_names - ds_saemix_all <- purrr::map_dfr(ds_list, function(x) x, .id = "ds") + ds_saemix_all <- vctrs::vec_rbind(!!!ds_list, .names_to = "ds") ds_saemix <- data.frame(ds = ds_saemix_all$ds, name = as.character(ds_saemix_all$variable), time = ds_saemix_all$time, value = ds_saemix_all$observed, stringsAsFactors = FALSE) + if (!is.null(covariates)) { + name.covariates <- names(covariates) + covariates$ds <- rownames(covariates) + ds_saemix <- merge(ds_saemix, covariates, sort = FALSE) + } else { + name.covariates <- character(0) + } res <- saemix::saemixData(ds_saemix, name.group = "ds", name.predictors = c("time", "name"), name.response = "value", + name.covariates = name.covariates, verbose = verbose, ...) return(res) } + +#' logLik method for saem.mmkin objects +#' +#' @param object The fitted [saem.mmkin] object +#' @param \dots Passed to [saemix::logLik.SaemixObject] +#' @param method Passed to [saemix::logLik.SaemixObject] +#' @export +logLik.saem.mmkin <- function(object, ..., method = c("is", "lin", "gq")) { + method <- match.arg(method) + return(logLik(object$so, method = method)) +} + +#' @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 so we do not have to export saem.mmkin in + # addition to the S3 method + call[[1]] <- saem + + # We also need to provide the mmkin object in the call, so it + # will also be found when called by testthat or pkgdown + call[[2]] <- object$mmkin + + 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 +} + +#' @export +#' @rdname saem +#' @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, ...) { + 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/R/set_nd_nq.R b/R/set_nd_nq.R new file mode 100644 index 00000000..37b9a894 --- /dev/null +++ b/R/set_nd_nq.R @@ -0,0 +1,164 @@ +#' Set non-detects and unquantified values in residue series without replicates +#' +#' This function automates replacing unquantified values in residue time and +#' depth series. For time series, the function performs part of the residue +#' processing proposed in the FOCUS kinetics guidance for parent compounds +#' and metabolites. For two-dimensional residue series over time and depth, +#' it automates the proposal of Boesten et al (2015). +#' +#' @param res_raw Character vector of a residue time series, or matrix of +#' residue values with rows representing depth profiles for a specific sampling +#' time, and columns representing time series of residues at the same depth. +#' Values below the limit of detection (lod) have to be coded as "nd", values +#' between the limit of detection and the limit of quantification, if any, have +#' to be coded as "nq". Samples not analysed have to be coded as "na". All +#' values that are not "na", "nd" or "nq" have to be coercible to numeric +#' @param lod Limit of detection (numeric) +#' @param loq Limit of quantification(numeric). Must be specified if the FOCUS rule to +#' stop after the first non-detection is to be applied +#' @param time_zero_presence Do we assume that residues occur at time zero? +#' This only affects samples from the first sampling time that have been +#' reported as "nd" (not detected). +#' @references Boesten, J. J. T. I., van der Linden, A. M. A., Beltman, W. H. +#' J. and Pol, J. W. (2015). Leaching of plant protection products and their +#' transformation products; Proposals for improving the assessment of leaching +#' to groundwater in the Netherlands — Version 2. Alterra report 2630, Alterra +#' Wageningen UR (University & Research centre) +#' @references FOCUS (2014) Generic Guidance for Estimating Persistence and Degradation +#' Kinetics from Environmental Fate Studies on Pesticides in EU Registration, Version 1.1, +#' 18 December 2014, p. 251 +#' @return A numeric vector, if a vector was supplied, or a numeric matrix otherwise +#' @export +#' @examples +#' # FOCUS (2014) p. 75/76 and 131/132 +#' parent_1 <- c(.12, .09, .05, .03, "nd", "nd", "nd", "nd", "nd", "nd") +#' set_nd_nq(parent_1, 0.02) +#' parent_2 <- c(.12, .09, .05, .03, "nd", "nd", .03, "nd", "nd", "nd") +#' set_nd_nq(parent_2, 0.02) +#' set_nd_nq_focus(parent_2, 0.02, loq = 0.05) +#' parent_3 <- c(.12, .09, .05, .03, "nd", "nd", .06, "nd", "nd", "nd") +#' set_nd_nq(parent_3, 0.02) +#' set_nd_nq_focus(parent_3, 0.02, loq = 0.05) +#' metabolite <- c("nd", "nd", "nd", 0.03, 0.06, 0.10, 0.11, 0.10, 0.09, 0.05, 0.03, "nd", "nd") +#' set_nd_nq(metabolite, 0.02) +#' set_nd_nq_focus(metabolite, 0.02, 0.05) +#' # +#' # Boesten et al. (2015), p. 57/58 +#' table_8 <- matrix( +#' c(10, 10, rep("nd", 4), +#' 10, 10, rep("nq", 2), rep("nd", 2), +#' 10, 10, 10, "nq", "nd", "nd", +#' "nq", 10, "nq", rep("nd", 3), +#' "nd", "nq", "nq", rep("nd", 3), +#' rep("nd", 6), rep("nd", 6)), +#' ncol = 6, byrow = TRUE) +#' set_nd_nq(table_8, 0.5, 1.5, time_zero_presence = TRUE) +#' table_10 <- matrix( +#' c(10, 10, rep("nd", 4), +#' 10, 10, rep("nd", 4), +#' 10, 10, 10, rep("nd", 3), +#' "nd", 10, rep("nd", 4), +#' rep("nd", 18)), +#' ncol = 6, byrow = TRUE) +#' set_nd_nq(table_10, 0.5, time_zero_presence = TRUE) +set_nd_nq <- function(res_raw, lod, loq = NA, time_zero_presence = FALSE) { + if (!is.character(res_raw)) { + stop("Please supply a vector or a matrix of character values") + } + if (is.vector(res_raw)) { + was_vector <- TRUE + res_raw <- as.matrix(res_raw) + } else { + was_vector <- FALSE + if (!is.matrix(res_raw)) { + stop("Please supply a vector or a matrix of character values") + } + } + nq <- 0.5 * (loq + lod) + nda <- 0.5 * lod # not detected but adjacent to detection + res_raw[res_raw == "nq"] <- nq + + if (!time_zero_presence) { + for (j in 1:ncol(res_raw)) { + if (res_raw[1, j] == "nd") res_raw[1, j] <- "na" + } + } + res_raw[res_raw == "na"] <- NA + + not_nd_na <- function(value) !(grepl("nd", value) | is.na(value)) + + for (i in 1:nrow(res_raw)) { + for (j in 1:ncol(res_raw)) { + if (!is.na(res_raw[i, j]) && res_raw[i, j] == "nd") { + if (i > 1) { # check earlier sample in same layer + if (not_nd_na(res_raw[i - 1, j])) res_raw[i, j] <- "nda" + } + if (i < nrow(res_raw)) { # check later sample + if (not_nd_na(res_raw[i + 1, j])) res_raw[i, j] <- "nda" + } + if (j > 1) { # check above sample at the same time + if (not_nd_na(res_raw[i, j - 1])) res_raw[i, j] <- "nda" + } + if (j < ncol(res_raw)) { # check sample below at the same time + if (not_nd_na(res_raw[i, j + 1])) res_raw[i, j] <- "nda" + } + } + } + } + res_raw[res_raw == "nda"] <- nda + res_raw[res_raw == "nd"] <- NA + + result <- as.numeric(res_raw) + dim(result) <- dim(res_raw) + dimnames(result) <- dimnames(res_raw) + if (was_vector) result <- as.vector(result) + return(result) +} + +#' @describeIn set_nd_nq Set non-detects in residue time series according to FOCUS rules +#' @param set_first_sample_nd Should the first sample be set to "first_sample_nd_value" +#' in case it is a non-detection? +#' @param first_sample_nd_value Value to be used for the first sample if it is a non-detection +#' @param ignore_below_loq_after_first_nd Should we ignore values below the LOQ after the first +#' non-detection that occurs after the quantified values? +#' @export +set_nd_nq_focus <- function(res_raw, lod, loq = NA, + set_first_sample_nd = TRUE, first_sample_nd_value = 0, + ignore_below_loq_after_first_nd = TRUE) +{ + + if (!is.vector(res_raw)) stop("FOCUS rules are only specified for one-dimensional time series") + + if (ignore_below_loq_after_first_nd & is.na(loq)) { + stop("You need to specify an LOQ") + } + + n <- length(res_raw) + if (ignore_below_loq_after_first_nd) { + for (i in 3:n) { + if (!res_raw[i - 2] %in% c("na", "nd")) { + if (res_raw[i - 1] == "nd") { + res_remaining <- res_raw[i:n] + res_remaining_unquantified <- ifelse(res_remaining == "na", TRUE, + ifelse(res_remaining == "nd", TRUE, + ifelse(res_remaining == "nq", TRUE, + ifelse(suppressWarnings(as.numeric(res_remaining)) < loq, TRUE, FALSE)))) + res_remaining_numeric <- suppressWarnings(as.numeric(res_remaining)) + res_remaining_below_loq <- ifelse(res_remaining == "nq", TRUE, + ifelse(!is.na(res_remaining_numeric) & res_remaining_numeric < loq, TRUE, FALSE)) + if (all(res_remaining_unquantified)) { + res_raw[i:n] <- ifelse(res_remaining_below_loq, "nd", res_remaining) + } + } + } + } + } + + result <- set_nd_nq(res_raw, lod = lod, loq = loq) + + if (set_first_sample_nd) { + if (res_raw[1] == "nd") result[1] <- first_sample_nd_value + } + + return(result) +} diff --git a/R/status.R b/R/status.R new file mode 100644 index 00000000..8bcd3a16 --- /dev/null +++ b/R/status.R @@ -0,0 +1,113 @@ +#' Method to get status information for fit array objects +#' +#' @param object The object to investigate +#' @param x The object to be printed +#' @param \dots For potential future extensions +#' @return An object with the same dimensions as the fit array +#' suitable printing method. +#' @export +status <- function(object, ...) +{ + UseMethod("status", object) +} + +#' @rdname status +#' @export +#' @examples +#' \dontrun{ +#' fits <- mmkin( +#' c("SFO", "FOMC"), +#' list("FOCUS A" = FOCUS_2006_A, +#' "FOCUS B" = FOCUS_2006_C), +#' quiet = TRUE) +#' status(fits) +#' } +status.mmkin <- function(object, ...) { + all_summary_warnings <- character() + sww <- 0 # Counter for Shapiro-Wilks warnings + + result <- lapply(object, + function(fit) { + if (inherits(fit, "try-error")) return("E") + sw <- fit$summary_warnings + swn <- names(sw) + if (length(sw) > 0) { + if (any(grepl("S", swn))) { + sww <<- sww + 1 + swn <- gsub("S", paste0("S", sww), swn) + } + warnstring <- paste(swn, collapse = ", ") + names(sw) <- swn + all_summary_warnings <<- c(all_summary_warnings, sw) + return(warnstring) + } else { + return("OK") + } + }) + result <- unlist(result) + dim(result) <- dim(object) + dimnames(result) <- dimnames(object) + + u_swn <- unique(names(all_summary_warnings)) + attr(result, "unique_warnings") <- all_summary_warnings[u_swn] + class(result) <- "status.mmkin" + return(result) +} + +#' @rdname status +#' @export +print.status.mmkin <- function(x, ...) { + u_w <- attr(x, "unique_warnings") + attr(x, "unique_warnings") <- NULL + class(x) <- NULL + print(x, quote = FALSE) + cat("\n") + for (i in seq_along(u_w)) { + cat(names(u_w)[i], ": ", u_w[i], "\n", sep = "") + } + if (any(x == "OK")) cat("OK: No warnings\n") + if (any(x == "E")) cat("E: Error\n") +} + +#' @rdname status +#' @export +status.mhmkin <- function(object, ...) { + if (inherits(object[[1]], "saem.mmkin")) { + test_func <- function(fit) { + if (inherits(fit$so, "try-error")) { + return("E") + } else { + if (!is.null(fit$FIM_failed)) { + return_values <- c("fixed effects" = "Fth", + "random effects and error model parameters" = "FO") + return(paste(return_values[fit$FIM_failed], collapse = ", ")) + } else { + return("OK") + } + } + } + } else { + stop("Only mhmkin objects containing saem.mmkin objects currently supported") + } + result <- lapply(object, test_func) + result <- unlist(result) + dim(result) <- dim(object) + dimnames(result) <- dimnames(object) + + class(result) <- "status.mhmkin" + return(result) +} + +#' @rdname status +#' @export +print.status.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 == "Fth")) cat("Fth: Could not invert FIM for fixed effects\n") + if (any(x == "FO")) cat("FO: Could not invert FIM for random effects and error model parameters\n") + if (any(x == "Fth, FO")) cat("Fth, FO: Could not invert FIM for fixed effects, nor for random effects and error model parameters\n") + if (any(x == "E")) cat("E: Error\n") +} + diff --git a/R/summary.mkinfit.R b/R/summary.mkinfit.R index f9858c32..4122873f 100644 --- a/R/summary.mkinfit.R +++ b/R/summary.mkinfit.R @@ -47,7 +47,7 @@ #' \url{http://esdac.jrc.ec.europa.eu/projects/degradation-kinetics} #' @examples #' -#' summary(mkinfit(mkinmod(parent = mkinsub("SFO")), FOCUS_2006_A, quiet = TRUE)) +#' summary(mkinfit("SFO", FOCUS_2006_A, quiet = TRUE)) #' #' @export summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, ...) { @@ -172,7 +172,7 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, if (distimes) ans$distimes <- ep$distimes if (length(ep$SFORB) != 0) ans$SFORB <- ep$SFORB if (!is.null(object$d_3_message)) ans$d_3_message <- object$d_3_message - class(ans) <- c("summary.mkinfit", "summary.modFit") + class(ans) <- "summary.mkinfit" return(ans) } diff --git a/R/summary.mmkin.R b/R/summary.mmkin.R new file mode 100644 index 00000000..06472e18 --- /dev/null +++ b/R/summary.mmkin.R @@ -0,0 +1,56 @@ +#' Summary method for class "mmkin" +#' +#' Shows status information on the [mkinfit] objects contained in the object +#' and gives an overview of ill-defined parameters calculated by [illparms]. +#' +#' @param object an object of class [mmkin] +#' @param x an object of class \code{summary.mmkin}. +#' @param conf.level confidence level for testing parameters +#' @param digits number of digits to use for printing +#' @param \dots optional arguments passed to methods like \code{print}. +#' @examples +#' +#' fits <- mmkin( +#' c("SFO", "FOMC"), +#' list("FOCUS A" = FOCUS_2006_A, +#' "FOCUS C" = FOCUS_2006_C), +#' quiet = TRUE, cores = 1) +#' summary(fits) +#' +#' @export +summary.mmkin <- function(object, conf.level = 0.95, ...) { + + ans <- list( + err_mod = object[[1, 1]]$err_mod, + time = attr(object, "time"), + illparms = illparms(object), + status = status(object) + ) + + class(ans) <- c("summary.mmkin") + return(ans) +} + +#' @rdname summary.mmkin +#' @export +print.summary.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) { + if (!is.null(x$err_mod)) { + cat("Error model: ") + cat(switch(x$err_mod, + const = "Constant variance", + obs = "Variance unique to each observed variable", + tc = "Two-component variance function"), "\n") + } + cat("Fitted in", x$time[["elapsed"]], "s\n") + + cat("\nStatus:\n") + print(x$status) + + if (any(x$illparms != "")) { + cat("\nIll-defined parameters:\n") + print(x$illparms) + } + + invisible(x) +} + diff --git a/R/summary.saem.mmkin.R b/R/summary.saem.mmkin.R index fa52a579..651cc29e 100644 --- a/R/summary.saem.mmkin.R +++ b/R/summary.saem.mmkin.R @@ -73,7 +73,12 @@ #' f_mmkin_dfop_sfo <- mmkin(list(dfop_sfo), ds_syn_dfop_sfo, #' quiet = TRUE, error_model = "tc", cores = 5) #' f_saem_dfop_sfo <- saem(f_mmkin_dfop_sfo) -#' summary(f_saem_dfop_sfo, data = TRUE) +#' print(f_saem_dfop_sfo) +#' illparms(f_saem_dfop_sfo) +#' f_saem_dfop_sfo_2 <- update(f_saem_dfop_sfo, covariance.model = diag(c(0, 0, 1, 1, 1, 0))) +#' illparms(f_saem_dfop_sfo_2) +#' intervals(f_saem_dfop_sfo_2) +#' summary(f_saem_dfop_sfo_2, data = TRUE) #' } #' #' @export @@ -82,18 +87,19 @@ summary.saem.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes = mod_vars <- names(object$mkinmod$diffs) pnames <- names(object$mean_dp_start) - np <- length(pnames) + names_fixed_effects <- object$so@results@name.fixed + n_fixed <- length(names_fixed_effects) conf.int <- object$so@results@conf.int rownames(conf.int) <- conf.int$name - confint_trans <- as.matrix(conf.int[pnames, c("estimate", "lower", "upper")]) + confint_trans <- as.matrix(parms(object, ci = TRUE)) colnames(confint_trans)[1] <- "est." # In case objects were produced by earlier versions of saem if (is.null(object$transformations)) object$transformations <- "mkin" if (object$transformations == "mkin") { - bp <- backtransform_odeparms(confint_trans[, "est."], object$mkinmod, + bp <- backtransform_odeparms(confint_trans[pnames, "est."], object$mkinmod, object$transform_rates, object$transform_fractions) bpnames <- names(bp) @@ -126,20 +132,20 @@ summary.saem.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes = } } } else { - confint_back <- confint_trans + confint_back <- confint_trans[names_fixed_effects, ] } # Correlation of fixed effects (inspired by summary.nlme) - varFix <- vcov(object$so)[1:np, 1:np] + varFix <- vcov(object$so)[1:n_fixed, 1:n_fixed] stdFix <- sqrt(diag(varFix)) object$corFixed <- array( t(varFix/stdFix)/stdFix, dim(varFix), - list(pnames, pnames)) + list(names_fixed_effects, names_fixed_effects)) # Random effects - rnames <- paste0("SD.", pnames) - confint_ranef <- as.matrix(conf.int[rnames, c("estimate", "lower", "upper")]) + sdnames <- intersect(rownames(conf.int), paste0("SD.", pnames)) + confint_ranef <- as.matrix(conf.int[sdnames, c("estimate", "lower", "upper")]) colnames(confint_ranef)[1] <- "est." # Error model @@ -147,7 +153,6 @@ summary.saem.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes = confint_errmod <- as.matrix(conf.int[enames, c("estimate", "lower", "upper")]) colnames(confint_errmod)[1] <- "est." - object$confint_trans <- confint_trans object$confint_ranef <- confint_ranef object$confint_errmod <- confint_errmod @@ -202,7 +207,7 @@ print.summary.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3) cat("\nModel predictions using solution type", x$solution_type, "\n") cat("\nFitted in", x$time[["elapsed"]], "s\n") - cat("Using", paste(x$so@options$nbiter.saemix, collapse = ", "), + cat("Using", paste(x$so@options$nbiter.saemix, collapse = ", "), "iterations and", x$so@options$nb.chains, "chains\n") cat("\nVariance model: ") |