diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/nlme.R | 88 | ||||
| -rw-r--r-- | R/nlme.mmkin.R | 100 | ||||
| -rw-r--r-- | R/plot.nlme.mmkin.R | 22 |
3 files changed, 131 insertions, 79 deletions
@@ -8,6 +8,7 @@ #' @param object An mmkin row object containing several fits of the same model to different datasets #' @import nlme #' @rdname nlme +#' @seealso \code{\link{nlme.mmkin}} #' @examples #' sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120) #' m_SFO <- mkinmod(parent = mkinsub("SFO")) @@ -47,73 +48,9 @@ #' start = mean_dp) #' summary(m_nlme) #' plot(augPred(m_nlme, level = 0:1), layout = c(3, 1)) +#' # augPred does not seem to work on fits with more than one state +#' # variable #' -#' \dontrun{ -#' # Test on some real data -#' ds_2 <- lapply(experimental_data_for_UBA_2019[6:10], -#' function(x) x$data[c("name", "time", "value")]) -#' m_sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"), -#' A1 = mkinsub("SFO"), use_of_ff = "min") -#' m_sfo_sfo_ff <- mkinmod(parent = mkinsub("SFO", "A1"), -#' A1 = mkinsub("SFO"), use_of_ff = "max") -#' m_fomc_sfo <- mkinmod(parent = mkinsub("FOMC", "A1"), -#' A1 = mkinsub("SFO")) -#' m_dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"), -#' A1 = mkinsub("SFO")) -#' m_sforb_sfo <- mkinmod(parent = mkinsub("SFORB", "A1"), -#' A1 = mkinsub("SFO")) -#' -#' f_2 <- mmkin(list("SFO-SFO" = m_sfo_sfo, -#' "SFO-SFO-ff" = m_sfo_sfo_ff, -#' "FOMC-SFO" = m_fomc_sfo, -#' "DFOP-SFO" = m_dfop_sfo, -#' "SFORB-SFO" = m_sforb_sfo), -#' ds_2) -#' -#' grouped_data_2 <- nlme_data(f_2["SFO-SFO", ]) -#' -#' mean_dp_sfo_sfo <- mean_degparms(f_2["SFO-SFO", ]) -#' mean_dp_sfo_sfo_ff <- mean_degparms(f_2["SFO-SFO-ff", ]) -#' mean_dp_fomc_sfo <- mean_degparms(f_2["FOMC-SFO", ]) -#' mean_dp_dfop_sfo <- mean_degparms(f_2["DFOP-SFO", ]) -#' mean_dp_sforb_sfo <- mean_degparms(f_2["SFORB-SFO", ]) -#' -#' nlme_f_sfo_sfo <- nlme_function(f_2["SFO-SFO", ]) -#' nlme_f_sfo_sfo_ff <- nlme_function(f_2["SFO-SFO-ff", ]) -#' nlme_f_fomc_sfo <- nlme_function(f_2["FOMC-SFO", ]) -#' assign("nlme_f_sfo_sfo", nlme_f_sfo_sfo, globalenv()) -#' assign("nlme_f_sfo_sfo_ff", nlme_f_sfo_sfo_ff, globalenv()) -#' assign("nlme_f_fomc_sfo", nlme_f_fomc_sfo, globalenv()) -#' -#' # Allowing for correlations between random effects (not shown) -#' # leads to non-convergence -#' f_nlme_sfo_sfo <- nlme(value ~ nlme_f_sfo_sfo(name, time, -#' parent_0, log_k_parent_sink, log_k_parent_A1, log_k_A1_sink), -#' data = grouped_data_2, -#' fixed = parent_0 + log_k_parent_sink + log_k_parent_A1 + log_k_A1_sink ~ 1, -#' random = pdDiag(parent_0 + log_k_parent_sink + log_k_parent_A1 + log_k_A1_sink ~ 1), -#' start = mean_dp_sfo_sfo) -#' # augPred does not see to work on this object, so no plot is shown -#' -#' # The same model fitted with transformed formation fractions does not converge -#' f_nlme_sfo_sfo_ff <- nlme(value ~ nlme_f_sfo_sfo_ff(name, time, -#' parent_0, log_k_parent, log_k_A1, f_parent_ilr_1), -#' data = grouped_data_2, -#' fixed = parent_0 + log_k_parent + log_k_A1 + f_parent_ilr_1 ~ 1, -#' random = pdDiag(parent_0 + log_k_parent + log_k_A1 + f_parent_ilr_1 ~ 1), -#' start = mean_dp_sfo_sfo_ff) -#' -#' f_nlme_fomc_sfo <- nlme(value ~ nlme_f_fomc_sfo(name, time, -#' parent_0, log_alpha, log_beta, log_k_A1, f_parent_ilr_1), -#' data = grouped_data_2, -#' fixed = parent_0 + log_alpha + log_beta + log_k_A1 + f_parent_ilr_1 ~ 1, -#' random = pdDiag(parent_0 + log_alpha + log_beta + log_k_A1 + f_parent_ilr_1 ~ 1), -#' start = mean_dp_fomc_sfo) -#' -#' # DFOP-SFO and SFORB-SFO did not converge with full random effects -#' -#' anova(f_nlme_fomc_sfo, f_nlme_sfo_sfo) -#' } #' @return A function that can be used with nlme #' @export nlme_function <- function(object) { @@ -185,14 +122,25 @@ nlme_function <- function(object) { } #' @rdname nlme -#' @return A named vector containing mean values of the fitted degradation model parameters +#' @return If random is FALSE (default), a named vector containing mean values +#' of the fitted degradation model parameters. If random is TRUE, a list with +#' fixed and random effects, in the format required by the start argument of +#' nlme for the case of a single grouping variable ds? +#' @param random Should a list with fixed and random effects be returned? #' @export -mean_degparms <- function(object) { +mean_degparms <- function(object, random = FALSE) { if (nrow(object) > 1) stop("Only row objects allowed") degparm_mat_trans <- sapply(object, parms, transformed = TRUE) mean_degparm_names <- setdiff(rownames(degparm_mat_trans), names(object[[1]]$errparms)) - res <- apply(degparm_mat_trans[mean_degparm_names, ], 1, mean) - return(res) + fixed <- apply(degparm_mat_trans[mean_degparm_names, ], 1, mean) + if (random) { + degparm_mat_trans[mean_degparm_names, ] + random <- t(apply(degparm_mat_trans[mean_degparm_names, ], 2, function(column) column - fixed)) + rownames(random) <- as.character(1:nrow(random)) + return(list(fixed = fixed, random = list(ds = random))) + } else { + return(fixed) + } } #' @rdname nlme diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R index 2ee46f33..6e3467ed 100644 --- a/R/nlme.mmkin.R +++ b/R/nlme.mmkin.R @@ -32,9 +32,52 @@ #' f <- mmkin("SFO", ds, quiet = TRUE, cores = 1) #' library(nlme) #' f_nlme <- nlme(f) -#' nlme(f, random = parent_0 ~ 1) -#' f_nlme <- nlme(f, start = c(parent_0 = 100, log_k_parent_sink = 0.1)) -#' update(f_nlme, random = parent_0 ~ 1) +#' print(f_nlme) +#' f_nlme_2 <- nlme(f, start = c(parent_0 = 100, log_k_parent_sink = 0.1)) +#' update(f_nlme_2, random = parent_0 ~ 1) +#' \dontrun{ +#' # Test on some real data +#' ds_2 <- lapply(experimental_data_for_UBA_2019[6:10], +#' function(x) x$data[c("name", "time", "value")]) +#' m_sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"), +#' A1 = mkinsub("SFO"), use_of_ff = "min", quiet = TRUE) +#' m_sfo_sfo_ff <- mkinmod(parent = mkinsub("SFO", "A1"), +#' A1 = mkinsub("SFO"), use_of_ff = "max", quiet = TRUE) +#' m_fomc_sfo <- mkinmod(parent = mkinsub("FOMC", "A1"), +#' A1 = mkinsub("SFO"), quiet = TRUE) +#' m_dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"), +#' A1 = mkinsub("SFO"), quiet = TRUE) +#' +#' f_2 <- mmkin(list("SFO-SFO" = m_sfo_sfo, +#' "SFO-SFO-ff" = m_sfo_sfo_ff, +#' "FOMC-SFO" = m_fomc_sfo, +#' "DFOP-SFO" = m_dfop_sfo), +#' ds_2, quiet = TRUE) +#' plot(f_2["SFO-SFO", 3:4]) # Separate fits for datasets 3 and 4 +#' +#' f_nlme_sfo_sfo <- nlme(f_2["SFO-SFO", ]) +#' # plot(f_nlme_sfo_sfo) # not feasible with pkgdown figures +#' plot(f_nlme_sfo_sfo, 3:4) # Global mixed model: Fits for datasets 3 and 4 +#' +#' # With formation fractions +#' f_nlme_sfo_sfo_ff <- nlme(f_2["SFO-SFO-ff", ]) +#' plot(f_nlme_sfo_sfo_ff, 3:4) # chi2 different due to different df attribution +#' +#' # For more parameters, we need to increase pnlsMaxIter and the tolerance +#' # to get convergence +#' f_nlme_fomc_sfo <- nlme(f_2["FOMC-SFO", ], +#' control = list(pnlsMaxIter = 100, tolerance = 1e-4), verbose = TRUE) +#' f_nlme_dfop_sfo <- nlme(f_2["DFOP-SFO", ], +#' control = list(pnlsMaxIter = 120, tolerance = 5e-4), verbose = TRUE) +#' plot(f_2["FOMC-SFO", 3:4]) +#' plot(f_nlme_fomc_sfo, 3:4) +#' +#' plot(f_2["DFOP-SFO", 3:4]) +#' plot(f_nlme_dfop_sfo, 3:4) +#' +#' anova(f_nlme_dfop_sfo, f_nlme_fomc_sfo, f_nlme_sfo_sfo) +#' anova(f_nlme_dfop_sfo, f_nlme_sfo_sfo) # if we ignore FOMC +#' } # Code inspired by nlme.nlsList nlme.mmkin <- function(model, data = sys.frame(sys.parent()), fixed, random = fixed, @@ -68,7 +111,7 @@ nlme.mmkin <- function(model, data = sys.frame(sys.parent()), thisCall[["data"]] <- nlme_data(model) if (missing(start)) { - thisCall[["start"]] <- mean_dp + thisCall[["start"]] <- mean_degparms(model, random = TRUE) } thisCall[["fixed"]] <- lapply(as.list(dp_names), function(el) @@ -84,3 +127,52 @@ nlme.mmkin <- function(model, data = sys.frame(sys.parent()), return(val) } +#' @export +#' @rdname nlme.mmkin +#' @param x An nlme.mmkin object to print +#' @param data Should the data be printed? +#' @param ... Further arguments as in the generic +print.nlme.mmkin <- function(x, ...) { + x$call$data <- "Not shown" + NextMethod("print", x) +} + +#' @export +#' @rdname nlme.mmkin +#' @param object An nlme.mmkin object to update +#' @param ... Update specifications passed to update.nlme +update.nlme.mmkin <- function(object, ...) { + res <- NextMethod() + res$mmkin_orig <- object$mmkin_orig + class(res) <- c("nlme.mmkin", "nlme", "lme") + return(res) +} + +# The following is necessary as long as R bug 17761 is not fixed +# https://bugs.r-project.org/bugzilla/show_bug.cgi?id=17761 +#' @export +anova.nlme.mmkin <- function(object, ...) { + thisCall <- as.list(match.call())[-1] + object_name <- as.character(thisCall[[1]]) + other_object_names <- sapply(thisCall[-1], as.character) + + remove_class <- function(object, classname) { + old_class <- class(object) + class(object) <- setdiff(old_class, classname) + return(object) + } + object <- remove_class(object, "nlme.mmkin") + other_objects <- list(...) + other_objects <- lapply(other_objects, remove_class, "nlme.mmkin") + + env <- new.env() + assign(object_name, object, env) + for (i in seq_along(other_objects)) { + assign(other_object_names[i], other_objects[[i]], env) + } + res <- eval(parse(text = paste0("anova.lme(", object_name, ", ", + paste(other_object_names, collapse = ", "), ")")), + envir = env) + + return(res) +} diff --git a/R/plot.nlme.mmkin.R b/R/plot.nlme.mmkin.R index ef6d100a..0f3ad715 100644 --- a/R/plot.nlme.mmkin.R +++ b/R/plot.nlme.mmkin.R @@ -11,6 +11,12 @@ #' @param standardized Should the residuals be standardized? This option #' is passed to \code{\link{mkinresplot}}, it only takes effect if #' `resplot = "time"`. +#' @param show_errmin Should the chi2 error level be shown on top of the plots +#' to the left? +#' @param errmin_var The variable for which the FOCUS chi2 error value should +#' be shown. +#' @param errmin_digits The number of significant digits for rounding the FOCUS +#' chi2 error percentage. #' @param cex Passed to the plot functions and \code{\link{mtext}}. #' @param rel.height.middle The relative height of the middle plot, if more #' than two rows of plots are shown. @@ -25,16 +31,19 @@ #' function(x) subset(x$data[c("name", "time", "value")], name == "parent")) #' f <- mmkin("SFO", ds, quiet = TRUE, cores = 1) #' #plot(f) # too many panels for pkgdown +#' plot(f[, 3:4]) #' library(nlme) #' f_nlme <- nlme(f) #' #' #plot(f_nlme) # too many panels for pkgdown -#' plot(f_nlme, 1:2) +#' plot(f_nlme, 3:4) #' @export plot.nlme.mmkin <- function(x, i = 1:ncol(x$mmkin_orig), main = "auto", legends = 1, resplot = c("time", "errmod"), standardized = FALSE, + show_errmin = TRUE, + errmin_var = "All data", errmin_digits = 3, cex = 0.7, rel.height.middle = 0.9, ymax = "auto", ...) { @@ -79,13 +88,14 @@ plot.nlme.mmkin <- function(x, i = 1:ncol(x$mmkin_orig), state_ini[names(odeini_optim)] <- odeini_optim odeparms <- fit_a$bparms.ode - odeparms[names(odeparms)] <- odeparms_optim + odeparms[names(odeparms_optim)] <- odeparms_optim mkinfit_call[["observed"]] <- ds[[a]] mkinfit_call[["parms.ini"]] <- odeparms mkinfit_call[["state.ini"]] <- state_ini - mkinfit_call[["control"]] <- list(iter.max = 1) + mkinfit_call[["control"]] <- list(iter.max = 0) + mkinfit_call[["quiet"]] <- TRUE res <- suppressWarnings(do.call("mkinfit", mkinfit_call)) return(res) @@ -94,9 +104,11 @@ plot.nlme.mmkin <- function(x, i = 1:ncol(x$mmkin_orig), # Set dimensions with names and the class (mmkin) attributes(mmkin_nlme) <- attributes(x$mmkin_orig[, i]) - plot(mmkin_nlme[, i], main = main, legends = legends, + plot(mmkin_nlme, main = main, legends = legends, resplot = resplot, standardized = standardized, - show_errmin = FALSE, cex = cex, + show_errmin = show_errmin, + errmin_var = errmin_var, errmin_digits = errmin_digits, + cex = cex, rel.height.middle = rel.height.middle, ymax = ymax, ...) |