From 48c463680b51fa767b4cd7bd62865f192d0354ac Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Sat, 6 Feb 2021 18:30:32 +0100 Subject: Reintroduce interface to saemix Also after the upgrade from buster to bullseye of my local system, some test results for saemix have changed. --- R/saem.R | 512 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 512 insertions(+) create mode 100644 R/saem.R (limited to 'R/saem.R') diff --git a/R/saem.R b/R/saem.R new file mode 100644 index 00000000..fd2a77b4 --- /dev/null +++ b/R/saem.R @@ -0,0 +1,512 @@ +utils::globalVariables(c("predicted", "std")) + +#' Fit nonlinear mixed models with SAEM +#' +#' This function uses [saemix::saemix()] as a backend for fitting nonlinear mixed +#' effects models created from [mmkin] row objects using the Stochastic Approximation +#' Expectation Maximisation algorithm (SAEM). +#' +#' An mmkin row object is essentially a list of mkinfit objects that have been +#' obtained by fitting the same model to a list of datasets using [mkinfit]. +#' +#' Starting values for the fixed effects (population mean parameters, argument +#' psi0 of [saemix::saemixModel()] are the mean values of the parameters found +#' using [mmkin]. +#' +#' @param object An [mmkin] row object containing several fits of the same +#' [mkinmod] model to different datasets +#' @param verbose Should we print information about created objects of +#' 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. +#' @param degparms_start Parameter values given as a named numeric vector will +#' be used to override the starting values obtained from the 'mmkin' object. +#' @param solution_type Possibility to specify the solution type in case the +#' automatic choice is not desired +#' @param quiet Should we suppress the messages saemix prints at the beginning +#' and the end of the optimisation process? +#' @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'. +#' @seealso [summary.saem.mmkin] [plot.mixed.mmkin] +#' @examples +#' \dontrun{ +#' ds <- lapply(experimental_data_for_UBA_2019[6:10], +#' function(x) subset(x$data[c("name", "time", "value")])) +#' names(ds) <- paste("Dataset", 6:10) +#' f_mmkin_parent_p0_fixed <- mmkin("FOMC", ds, +#' state.ini = c(parent = 100), fixed_initials = "parent", quiet = TRUE) +#' f_saem_p0_fixed <- saem(f_mmkin_parent_p0_fixed) +#' +#' f_mmkin_parent <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE) +#' f_saem_sfo <- saem(f_mmkin_parent["SFO", ]) +#' f_saem_fomc <- saem(f_mmkin_parent["FOMC", ]) +#' f_saem_dfop <- saem(f_mmkin_parent["DFOP", ]) +#' +#' # The returned saem.mmkin object contains an SaemixObject, therefore we can use +#' # functions from saemix +#' library(saemix) +#' compare.saemix(list(f_saem_sfo$so, f_saem_fomc$so, f_saem_dfop$so)) +#' plot(f_saem_fomc$so, plot.type = "convergence") +#' plot(f_saem_fomc$so, plot.type = "individual.fit") +#' plot(f_saem_fomc$so, plot.type = "npde") +#' plot(f_saem_fomc$so, plot.type = "vpc") +#' +#' f_mmkin_parent_tc <- update(f_mmkin_parent, error_model = "tc") +#' f_saem_fomc_tc <- saem(f_mmkin_parent_tc["FOMC", ]) +#' compare.saemix(list(f_saem_fomc$so, f_saem_fomc_tc$so)) +#' +#' sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"), +#' A1 = mkinsub("SFO")) +#' fomc_sfo <- mkinmod(parent = mkinsub("FOMC", "A1"), +#' A1 = mkinsub("SFO")) +#' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"), +#' A1 = mkinsub("SFO")) +#' # The following fit uses analytical solutions for SFO-SFO and DFOP-SFO, +#' # and compiled ODEs for FOMC that are much slower +#' f_mmkin <- mmkin(list( +#' "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo), +#' ds, quiet = TRUE) +#' # saem fits of SFO-SFO and DFOP-SFO to these data take about five seconds +#' # each on this system, as we use analytical solutions written for saemix. +#' # When using the analytical solutions written for mkin this took around +#' # four minutes +#' f_saem_sfo_sfo <- saem(f_mmkin["SFO-SFO", ]) +#' f_saem_dfop_sfo <- saem(f_mmkin["DFOP-SFO", ]) +#' # We can use print, plot and summary methods to check the results +#' print(f_saem_dfop_sfo) +#' plot(f_saem_dfop_sfo) +#' summary(f_saem_dfop_sfo, data = TRUE) +#' +#' # The following takes about 6 minutes +#' #f_saem_dfop_sfo_deSolve <- saem(f_mmkin["DFOP-SFO", ], solution_type = "deSolve", +#' # control = list(nbiter.saemix = c(200, 80), nbdisplay = 10)) +#' +#' #saemix::compare.saemix(list( +#' # f_saem_dfop_sfo$so, +#' # f_saem_dfop_sfo_deSolve$so)) +#' +#' # If the model supports it, we can also use eigenvalue based solutions, which +#' # take a similar amount of time +#' #f_saem_sfo_sfo_eigen <- saem(f_mmkin["SFO-SFO", ], solution_type = "eigen", +#' # control = list(nbiter.saemix = c(200, 80), nbdisplay = 10)) +#' } +#' @export +saem <- function(object, ...) UseMethod("saem") + +#' @rdname saem +#' @export +saem.mmkin <- function(object, + transformations = c("mkin", "saemix"), + degparms_start = numeric(), + solution_type = "auto", + control = list(displayProgress = FALSE, print = FALSE, + save = FALSE, save.graphs = FALSE), + verbose = FALSE, quiet = FALSE, ...) +{ + transformations <- match.arg(transformations) + m_saemix <- saemix_model(object, verbose = verbose, + degparms_start = degparms_start, solution_type = solution_type, + transformations = transformations, ...) + d_saemix <- saemix_data(object, verbose = verbose) + + fit_time <- system.time({ + utils::capture.output(f_saemix <- saemix::saemix(m_saemix, d_saemix, control), split = !quiet) + }) + + 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 + } + + return_data <- nlme_data(object) + + return_data$predicted <- f_saemix@model@model( + psi = saemix::psi(f_saemix), + id = as.numeric(return_data$ds), + xidep = return_data[c("time", "name")]) + + return_data <- transform(return_data, + residual = predicted - value, + 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, + mmkin = object, + solution_type = object[[1]]$solution_type, + transformations = transformations, + transform_rates = object[[1]]$transform_rates, + transform_fractions = object[[1]]$transform_fractions, + so = f_saemix, + time = fit_time, + mean_dp_start = attr(m_saemix, "mean_dp_start"), + bparms.optim = bparms_optim, + bparms.fixed = object[[1]]$bparms.fixed, + data = return_data, + err_mod = object[[1]]$err_mod, + date.fit = date(), + saemixversion = as.character(utils::packageVersion("saemix")), + mkinversion = as.character(utils::packageVersion("mkin")), + Rversion = paste(R.version$major, R.version$minor, sep=".") + ) + + class(result) <- c("saem.mmkin", "mixed.mmkin") + return(result) +} + +#' @export +#' @rdname saem +#' @param x An saem.mmkin object to print +#' @param digits Number of digits to use for printing +print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) { + cat( "Kinetic nonlinear mixed-effects model fit by SAEM" ) + cat("\nStructural model:\n") + diffs <- x$mmkin[[1]]$mkinmod$diffs + nice_diffs <- gsub("^(d.*) =", "\\1/dt =", diffs) + writeLines(strwrap(nice_diffs, exdent = 11)) + cat("\nData:\n") + cat(nrow(x$data), "observations of", + 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"]] + print(conf.int, digits = digits) + + invisible(x) +} + +#' @rdname saem +#' @return An [saemix::SaemixModel] object. +#' @export +saemix_model <- function(object, solution_type = "auto", transformations = c("mkin", "saemix"), + degparms_start = numeric(), verbose = FALSE, ...) +{ + if (nrow(object) > 1) stop("Only row objects allowed") + + mkin_model <- object[[1]]$mkinmod + + degparms_optim <- mean_degparms(object) + if (transformations == "saemix") { + degparms_optim <- backtransform_odeparms(degparms_optim, + object[[1]]$mkinmod, + object[[1]]$transform_rates, + object[[1]]$transform_fractions) + } + degparms_fixed <- object[[1]]$bparms.fixed + + # Transformations are done in the degradation function + transform.par = rep(0, length(degparms_optim)) + + odeini_optim_parm_names <- grep('_0$', names(degparms_optim), value = TRUE) + odeini_fixed_parm_names <- grep('_0$', names(degparms_fixed), value = TRUE) + + odeparms_fixed_names <- setdiff(names(degparms_fixed), odeini_fixed_parm_names) + odeparms_fixed <- degparms_fixed[odeparms_fixed_names] + + odeini_fixed <- degparms_fixed[odeini_fixed_parm_names] + names(odeini_fixed) <- gsub('_0$', '', odeini_fixed_parm_names) + + model_function <- FALSE + + # Model functions with analytical solutions + # Fixed parameters, use_of_ff = "min" and turning off sinks currently not supported here + # In general, we need to consider exactly how the parameters in mkinfit were specified, + # as the parameters are currently mapped by position in these solutions + sinks <- sapply(mkin_model$spec, function(x) x$sink) + if (length(odeparms_fixed) == 0 & mkin_model$use_of_ff == "max" & all(sinks)) { + # Parent only + if (length(mkin_model$spec) == 1) { + parent_type <- mkin_model$spec[[1]]$type + if (length(odeini_fixed) == 1) { + if (parent_type == "SFO") { + stop("saemix needs at least two parameters to work on.") + } + if (parent_type == "FOMC") { + model_function <- function(psi, id, xidep) { + odeini_fixed / (xidep[, "time"]/exp(psi[id, 2]) + 1)^exp(psi[id, 1]) + } + } + if (parent_type == "DFOP") { + model_function <- function(psi, id, xidep) { + g <- plogis(psi[id, 3]) + t <- xidep[, "time"] + odeini_fixed * (g * exp(- exp(psi[id, 1]) * t) + + (1 - g) * exp(- exp(psi[id, 2]) * t)) + } + } + if (parent_type == "HS") { + model_function <- function(psi, id, xidep) { + tb <- exp(psi[id, 3]) + t <- xidep[, "time"] + k1 = exp(psi[id, 1]) + odeini_fixed * ifelse(t <= tb, + exp(- k1 * t), + exp(- k1 * tb) * exp(- exp(psi[id, 2]) * (t - tb))) + } + } + } else { + if (parent_type == "SFO") { + if (transformations == "mkin") { + model_function <- function(psi, id, xidep) { + psi[id, 1] * exp( - exp(psi[id, 2]) * xidep[, "time"]) + } + } else { + model_function <- function(psi, id, xidep) { + psi[id, 1] * exp( - psi[id, 2] * xidep[, "time"]) + } + transform.par = c(0, 1) + } + } + 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 (parent_type == "DFOP") { + if (transformations == "mkin") { + model_function <- function(psi, id, xidep) { + g <- plogis(psi[id, 4]) + t <- xidep[, "time"] + psi[id, 1] * (g * exp(- exp(psi[id, 2]) * t) + + (1 - g) * exp(- exp(psi[id, 3]) * t)) + } + } else { + model_function <- function(psi, id, xidep) { + g <- psi[id, 4] + t <- xidep[, "time"] + psi[id, 1] * (g * exp(- psi[id, 2] * t) + + (1 - g) * exp(- psi[id, 3] * t)) + } + transform.par = c(0, 1, 1, 3) + } + } + 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))) + } + } + } + } + + # Parent with one metabolite + # Parameter names used in the model functions are as in + # https://nbviewer.jupyter.org/urls/jrwb.de/nb/Symbolic%20ODE%20solutions%20for%20mkin.ipynb + types <- unname(sapply(mkin_model$spec, function(x) x$type)) + if (length(mkin_model$spec) == 2 &! "SFORB" %in% types ) { + # Initial value for the metabolite (n20) must be fixed + if (names(odeini_fixed) == names(mkin_model$spec)[2]) { + n20 <- odeini_fixed + parent_name <- names(mkin_model$spec)[1] + if (identical(types, c("SFO", "SFO"))) { + if (transformations == "mkin") { + model_function <- function(psi, id, xidep) { + t <- xidep[, "time"] + n10 <- psi[id, 1] + k1 <- exp(psi[id, 2]) + k2 <- exp(psi[id, 3]) + f12 <- plogis(psi[id, 4]) + ifelse(xidep[, "name"] == parent_name, + n10 * exp(- k1 * t), + (((k2 - k1) * n20 - f12 * k1 * n10) * exp(- k2 * t)) / (k2 - k1) + + (f12 * k1 * n10 * exp(- k1 * t)) / (k2 - k1) + ) + } + } else { + model_function <- function(psi, id, xidep) { + t <- xidep[, "time"] + n10 <- psi[id, 1] + k1 <- psi[id, 2] + k2 <- psi[id, 3] + f12 <- psi[id, 4] + ifelse(xidep[, "name"] == parent_name, + n10 * exp(- k1 * t), + (((k2 - k1) * n20 - f12 * k1 * n10) * exp(- k2 * t)) / (k2 - k1) + + (f12 * k1 * n10 * exp(- k1 * t)) / (k2 - k1) + ) + } + transform.par = c(0, 1, 1, 3) + } + } + if (identical(types, c("DFOP", "SFO"))) { + if (transformations == "mkin") { + model_function <- function(psi, id, xidep) { + t <- xidep[, "time"] + n10 <- psi[id, 1] + k2 <- exp(psi[id, 2]) + f12 <- plogis(psi[id, 3]) + l1 <- exp(psi[id, 4]) + l2 <- exp(psi[id, 5]) + g <- plogis(psi[id, 6]) + ifelse(xidep[, "name"] == parent_name, + n10 * (g * exp(- l1 * t) + (1 - g) * exp(- l2 * t)), + ((f12 * g - f12) * l2 * n10 * exp(- l2 * t)) / (l2 - k2) - + (f12 * g * l1 * n10 * exp(- l1 * t)) / (l1 - k2) + + ((((l1 - k2) * l2 - k2 * l1 + k2^2) * n20 + + ((f12 * l1 + (f12 * g - f12) * k2) * l2 - + f12 * g * k2 * l1) * n10) * exp( - k2 * t)) / + ((l1 - k2) * l2 - k2 * l1 + k2^2) + ) + } + } else { + model_function <- function(psi, id, xidep) { + t <- xidep[, "time"] + n10 <- psi[id, 1] + k2 <- psi[id, 2] + f12 <- psi[id, 3] + l1 <- psi[id, 4] + l2 <- psi[id, 5] + g <- psi[id, 6] + ifelse(xidep[, "name"] == parent_name, + n10 * (g * exp(- l1 * t) + (1 - g) * exp(- l2 * t)), + ((f12 * g - f12) * l2 * n10 * exp(- l2 * t)) / (l2 - k2) - + (f12 * g * l1 * n10 * exp(- l1 * t)) / (l1 - k2) + + ((((l1 - k2) * l2 - k2 * l1 + k2^2) * n20 + + ((f12 * l1 + (f12 * g - f12) * k2) * l2 - + f12 * g * k2 * l1) * n10) * exp( - k2 * t)) / + ((l1 - k2) * l2 - k2 * l1 + k2^2) + ) + } + transform.par = c(0, 1, 3, 1, 1, 3) + } + } + } + } + } + + if (is.function(model_function) & solution_type == "auto") { + solution_type = "analytical saemix" + } else { + + if (solution_type == "auto") + solution_type <- object[[1]]$solution_type + + model_function <- function(psi, id, xidep) { + + uid <- unique(id) + + res_list <- lapply(uid, function(i) { + + transparms_optim <- as.numeric(psi[i, ]) # psi[i, ] is a dataframe when called in saemix.predict + names(transparms_optim) <- names(degparms_optim) + + odeini_optim <- transparms_optim[odeini_optim_parm_names] + names(odeini_optim) <- gsub('_0$', '', odeini_optim_parm_names) + + odeini <- c(odeini_optim, odeini_fixed)[names(mkin_model$diffs)] + + ode_transparms_optim_names <- setdiff(names(transparms_optim), odeini_optim_parm_names) + odeparms_optim <- backtransform_odeparms(transparms_optim[ode_transparms_optim_names], mkin_model, + transform_rates = object[[1]]$transform_rates, + transform_fractions = object[[1]]$transform_fractions) + odeparms <- c(odeparms_optim, odeparms_fixed) + + xidep_i <- subset(xidep, id == i) + + if (solution_type == "analytical") { + out_values <- mkin_model$deg_func(xidep_i, odeini, odeparms) + } else { + + i_time <- xidep_i$time + i_name <- xidep_i$name + + out_wide <- mkinpredict(mkin_model, + odeparms = odeparms, odeini = odeini, + solution_type = solution_type, + outtimes = sort(unique(i_time)), + na_stop = FALSE + ) + + out_index <- cbind(as.character(i_time), as.character(i_name)) + out_values <- out_wide[out_index] + } + return(out_values) + }) + res <- unlist(res_list) + return(res) + } + } + + error.model <- switch(object[[1]]$err_mod, + const = "constant", + tc = "combined", + obs = "constant") + + if (object[[1]]$err_mod == "obs") { + warning("The error model 'obs' (variance by variable) can currently not be transferred to an saemix model") + } + + error.init <- switch(object[[1]]$err_mod, + const = c(a = mean(sapply(object, function(x) x$errparms)), b = 1), + tc = c(a = mean(sapply(object, function(x) x$errparms[1])), + b = mean(sapply(object, function(x) x$errparms[2]))), + obs = c(a = mean(sapply(object, function(x) x$errparms)), b = 1)) + + degparms_psi0 <- degparms_optim + degparms_psi0[names(degparms_start)] <- degparms_start + psi0_matrix <- matrix(degparms_psi0, nrow = 1) + colnames(psi0_matrix) <- names(degparms_psi0) + + res <- saemix::saemixModel(model_function, + psi0 = psi0_matrix, + "Mixed model generated from mmkin object", + transform.par = transform.par, + error.model = error.model, + verbose = verbose + ) + attr(res, "mean_dp_start") <- degparms_optim + return(res) +} + +#' @rdname saem +#' @return An [saemix::SaemixData] object. +#' @export +saemix_data <- function(object, 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 <- 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) + + res <- saemix::saemixData(ds_saemix, + name.group = "ds", + name.predictors = c("time", "name"), + name.response = "value", + verbose = verbose, + ...) + return(res) +} -- cgit v1.2.3 From c73b2f30ec836c949885784ab576e814eb8070a9 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 9 Mar 2021 17:35:47 +0100 Subject: Some improvements for borderline cases - fit_with_errors for saem() - test_log_parms for mean_degparms() and saem() --- NEWS.md | 2 + R/nlme.R | 37 +- R/nlme.mmkin.R | 2 +- R/saem.R | 31 +- build.log | 2 +- check.log | 6 +- docs/dev/404.html | 2 +- docs/dev/articles/index.html | 2 +- docs/dev/authors.html | 2 +- docs/dev/index.html | 7 +- docs/dev/news/index.html | 161 +- docs/dev/pkgdown.yml | 2 +- docs/dev/reference/Rplot001.png | Bin 13995 -> 1011 bytes docs/dev/reference/Rplot002.png | Bin 13648 -> 16859 bytes docs/dev/reference/Rplot003.png | Bin 28745 -> 28844 bytes docs/dev/reference/Rplot004.png | Bin 49269 -> 49360 bytes docs/dev/reference/Rplot005.png | Bin 59143 -> 59216 bytes docs/dev/reference/endpoints.html | 2 +- docs/dev/reference/index.html | 2 +- docs/dev/reference/nlme-1.png | Bin 70133 -> 68233 bytes docs/dev/reference/nlme-2.png | Bin 94031 -> 90552 bytes docs/dev/reference/nlme.html | 33 +- docs/dev/reference/nlme.mmkin-1.png | Bin 124677 -> 124827 bytes docs/dev/reference/nlme.mmkin-2.png | Bin 169523 -> 169698 bytes docs/dev/reference/nlme.mmkin-3.png | Bin 172692 -> 172809 bytes docs/dev/reference/nlme.mmkin.html | 2 +- docs/dev/reference/plot.mixed.mmkin-1.png | Bin 84734 -> 85433 bytes docs/dev/reference/plot.mixed.mmkin-2.png | Bin 173916 -> 174061 bytes docs/dev/reference/plot.mixed.mmkin-3.png | Bin 172396 -> 172540 bytes docs/dev/reference/plot.mixed.mmkin-4.png | Bin 175502 -> 175594 bytes docs/dev/reference/plot.mixed.mmkin.html | 6 +- docs/dev/reference/saem-1.png | Bin 47315 -> 47342 bytes docs/dev/reference/saem-2.png | Bin 48720 -> 48819 bytes docs/dev/reference/saem-3.png | Bin 82107 -> 82202 bytes docs/dev/reference/saem-4.png | Bin 128231 -> 128213 bytes docs/dev/reference/saem-5.png | Bin 173288 -> 173665 bytes docs/dev/reference/saem.html | 72 +- docs/dev/reference/summary.saem.mmkin.html | 422 +- man/nlme.Rd | 11 +- man/saem.Rd | 19 +- test.log | 36 +- .../plotting/mixed-model-fit-for-nlme-object.svg | 2402 +++++------ ...t-for-saem-object-with-mkin-transformations.svg | 4527 ++++++++++---------- ...for-saem-object-with-saemix-transformations.svg | 5 + tests/testthat/print_mmkin_biphasic_mixed.txt | 6 +- tests/testthat/print_nlme_biphasic.txt | 10 +- tests/testthat/print_sfo_saem_1.txt | 16 +- tests/testthat/setup_script.R | 19 +- tests/testthat/summary_nlme_biphasic_s.txt | 46 +- tests/testthat/summary_saem_biphasic_s.txt | 48 +- tests/testthat/test_mixed.R | 24 +- tests/testthat/test_nlme.R | 2 +- 52 files changed, 4040 insertions(+), 3926 deletions(-) (limited to 'R/saem.R') diff --git a/NEWS.md b/NEWS.md index 38cac245..5d0ea69a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -12,6 +12,8 @@ - 'saem': generic function to fit saemix models using 'saemix_model' and 'saemix_data', with a generator 'saem.mmkin', summary and plot methods +- 'mean_degparms': New argument 'test_log_parms' that makes the function only consider log-transformed parameters where the untransformed parameters pass the t-test for a certain confidence level. This can be used to check more plausible parameters for 'saem' + # mkin 1.0.4 (Unreleased) - 'plot.mixed.mmkin': Reset graphical parameters on exit diff --git a/R/nlme.R b/R/nlme.R index 9215aab0..d235a094 100644 --- a/R/nlme.R +++ b/R/nlme.R @@ -36,7 +36,7 @@ #' nlme_f <- nlme_function(f) #' # These assignments are necessary for these objects to be #' # visible to nlme and augPred when evaluation is done by -#' # pkgdown to generated the html docs. +#' # pkgdown to generate the html docs. #' assign("nlme_f", nlme_f, globalenv()) #' assign("grouped_data", grouped_data, globalenv()) #' @@ -130,13 +130,44 @@ nlme_function <- function(object) { #' 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? +#' @param test_log_parms If TRUE, log parameters are only considered in +#' the mean calculations if their untransformed counterparts (most likely +#' rate constants) pass the t-test for significant difference from zero. +#' @param conf.level Possibility to adjust the required confidence level +#' for parameter that are tested if requested by 'test_log_parms'. #' @export -mean_degparms <- function(object, random = FALSE) { +mean_degparms <- function(object, random = FALSE, test_log_parms = FALSE, conf.level = 0.6) +{ if (nrow(object) > 1) stop("Only row objects allowed") parm_mat_trans <- sapply(object, parms, transformed = TRUE) + + if (test_log_parms) { + parm_mat_dim <- dim(parm_mat_trans) + parm_mat_dimnames <- dimnames(parm_mat_trans) + + log_parm_trans_names <- grep("^log_", rownames(parm_mat_trans), value = TRUE) + log_parm_names <- gsub("^log_", "", log_parm_trans_names) + + t_test_back_OK <- matrix( + sapply(object, function(o) { + suppressWarnings(summary(o)$bpar[log_parm_names, "Pr(>t)"] < (1 - conf.level)) + }), nrow = length(log_parm_names)) + rownames(t_test_back_OK) <- log_parm_trans_names + + parm_mat_trans_OK <- parm_mat_trans + for (trans_parm in log_parm_trans_names) { + parm_mat_trans_OK[trans_parm, ] <- ifelse(t_test_back_OK[trans_parm, ], + parm_mat_trans[trans_parm, ], NA) + } + } else { + parm_mat_trans_OK <- parm_mat_trans + } + mean_degparm_names <- setdiff(rownames(parm_mat_trans), names(object[[1]]$errparms)) degparm_mat_trans <- parm_mat_trans[mean_degparm_names, , drop = FALSE] - fixed <- apply(degparm_mat_trans, 1, mean) + degparm_mat_trans_OK <- parm_mat_trans_OK[mean_degparm_names, , drop = FALSE] + + fixed <- apply(degparm_mat_trans_OK, 1, mean, na.rm = TRUE) if (random) { random <- t(apply(degparm_mat_trans[mean_degparm_names, , drop = FALSE], 2, function(column) column - fixed)) # If we only have one parameter, apply returns a vector so we get a single row diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R index ff1f2fff..306600c6 100644 --- a/R/nlme.mmkin.R +++ b/R/nlme.mmkin.R @@ -24,7 +24,7 @@ get_deg_func <- function() { #' This functions sets up a nonlinear mixed effects model for an mmkin row #' object. An mmkin row object is essentially a list of mkinfit objects that #' have been obtained by fitting the same model to a list of datasets. -#' +#' #' Note that the convergence of the nlme algorithms depends on the quality #' of the data. In degradation kinetics, we often only have few datasets #' (e.g. data for few soils) and complicated degradation models, which may diff --git a/R/saem.R b/R/saem.R index fd2a77b4..460edede 100644 --- a/R/saem.R +++ b/R/saem.R @@ -24,8 +24,16 @@ utils::globalVariables(c("predicted", "std")) #' SFO or DFOP is used for the parent and there is either no metabolite or one. #' @param degparms_start Parameter values given as a named numeric vector will #' 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]. +#' @param conf.level Possibility to adjust the required confidence level +#' 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 +#' @param fail_with_errors Should a failure to compute standard errors +#' 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? #' @param control Passed to [saemix::saemix] @@ -51,7 +59,7 @@ utils::globalVariables(c("predicted", "std")) #' # The returned saem.mmkin object contains an SaemixObject, therefore we can use #' # functions from saemix #' library(saemix) -#' compare.saemix(list(f_saem_sfo$so, f_saem_fomc$so, f_saem_dfop$so)) +#' compare.saemix(f_saem_sfo$so, f_saem_fomc$so, f_saem_dfop$so) #' plot(f_saem_fomc$so, plot.type = "convergence") #' plot(f_saem_fomc$so, plot.type = "individual.fit") #' plot(f_saem_fomc$so, plot.type = "npde") @@ -59,7 +67,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(list(f_saem_fomc$so, f_saem_fomc_tc$so)) +#' compare.saemix(f_saem_fomc$so, f_saem_fomc_tc$so) #' #' sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"), #' A1 = mkinsub("SFO")) @@ -104,19 +112,32 @@ saem <- function(object, ...) UseMethod("saem") saem.mmkin <- function(object, transformations = c("mkin", "saemix"), degparms_start = numeric(), + test_log_parms = FALSE, + conf.level = 0.6, solution_type = "auto", control = list(displayProgress = FALSE, print = FALSE, save = FALSE, save.graphs = FALSE), + fail_with_errors = TRUE, verbose = FALSE, quiet = FALSE, ...) { transformations <- match.arg(transformations) m_saemix <- saemix_model(object, verbose = verbose, - degparms_start = degparms_start, solution_type = solution_type, + 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) fit_time <- system.time({ utils::capture.output(f_saemix <- saemix::saemix(m_saemix, d_saemix, control), split = !quiet) + FIM_failed <- NULL + 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") + } + 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 @@ -203,13 +224,13 @@ print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) { #' @return An [saemix::SaemixModel] object. #' @export saemix_model <- function(object, solution_type = "auto", transformations = c("mkin", "saemix"), - degparms_start = numeric(), verbose = FALSE, ...) + degparms_start = numeric(), test_log_parms = FALSE, verbose = FALSE, ...) { if (nrow(object) > 1) stop("Only row objects allowed") mkin_model <- object[[1]]$mkinmod - degparms_optim <- mean_degparms(object) + degparms_optim <- mean_degparms(object, test_log_parms = test_log_parms) if (transformations == "saemix") { degparms_optim <- backtransform_odeparms(degparms_optim, object[[1]]$mkinmod, diff --git a/build.log b/build.log index d50a4860..ca1c0481 100644 --- a/build.log +++ b/build.log @@ -6,5 +6,5 @@ * creating vignettes ... OK * checking for LF line-endings in source and make files and shell scripts * checking for empty or unneeded directories -* building ‘mkin_1.0.3.9000.tar.gz’ +* building ‘mkin_1.0.4.9000.tar.gz’ diff --git a/check.log b/check.log index ac59f6af..6e19f958 100644 --- a/check.log +++ b/check.log @@ -1,16 +1,16 @@ * using log directory ‘/home/jranke/git/mkin/mkin.Rcheck’ -* using R version 4.0.3 (2020-10-10) +* using R version 4.0.4 (2021-02-15) * using platform: x86_64-pc-linux-gnu (64-bit) * using session charset: UTF-8 * using options ‘--no-tests --as-cran’ * checking for file ‘mkin/DESCRIPTION’ ... OK * checking extension type ... Package -* this is package ‘mkin’ version ‘1.0.3.9000’ +* this is package ‘mkin’ version ‘1.0.4.9000’ * package encoding: UTF-8 * checking CRAN incoming feasibility ... NOTE Maintainer: ‘Johannes Ranke ’ -Version contains large components (1.0.3.9000) +Version contains large components (1.0.4.9000) Unknown, possibly mis-spelled, fields in DESCRIPTION: ‘Remotes’ diff --git a/docs/dev/404.html b/docs/dev/404.html index f9e51aa3..58591997 100644 --- a/docs/dev/404.html +++ b/docs/dev/404.html @@ -71,7 +71,7 @@ mkin - 1.0.3.9000 + 1.0.4.9000 diff --git a/docs/dev/articles/index.html b/docs/dev/articles/index.html index 17ee4a69..3c00526e 100644 --- a/docs/dev/articles/index.html +++ b/docs/dev/articles/index.html @@ -71,7 +71,7 @@ mkin - 1.0.3.9000 + 1.0.4.9000 diff --git a/docs/dev/authors.html b/docs/dev/authors.html index 63050c0d..45db18f2 100644 --- a/docs/dev/authors.html +++ b/docs/dev/authors.html @@ -71,7 +71,7 @@ mkin - 1.0.3.9000 + 1.0.4.9000 diff --git a/docs/dev/index.html b/docs/dev/index.html index 57328658..d1fa1a52 100644 --- a/docs/dev/index.html +++ b/docs/dev/index.html @@ -38,7 +38,7 @@ mkin - 1.0.3.9000 + 1.0.4.9000 @@ -192,11 +192,12 @@

Many inspirations for improvements of mkin resulted from doing kinetic evaluations of degradation data for my clients while working at Harlan Laboratories and at Eurofins Regulatory AG, and now as an independent consultant.

Funding was received from the Umweltbundesamt in the course of the projects

    -
  • Grant Number 112407 (Testing and validation of modelling software as an alternative to ModelMaker 4.0, 2014-2015)
    • +
  • Project Number 27452 (Testing and validation of modelling software as an alternative to ModelMaker 4.0, 2014-2015)
  • Project Number 56703 (Optimization of gmkin for routine use in the Umweltbundesamt, 2015)
    • +
  • Project Number 92570 (Update of Project Number 27452, 2017-2018)
  • Project Number 112407 (Testing the feasibility of using an error model according to Rocke and Lorenzato for more realistic parameter estimates in the kinetic evaluation of degradation data, 2018-2019)
  • Project Number 120667 (Development of objective criteria for the evaluation of the visual fit in the kinetic evaluation of degradation data, 2019-2020)
    • -
  • Project 146839 (Checking the feasibility of using mixed-effects models for the derivation of kinetic modelling parameters from degradation studies, 2020-2021)
    • +
  • Project Number 146839 (Checking the feasibility of using mixed-effects models for the derivation of kinetic modelling parameters from degradation studies, 2020-2021)
diff --git a/docs/dev/news/index.html b/docs/dev/news/index.html index 31c392f7..10585403 100644 --- a/docs/dev/news/index.html +++ b/docs/dev/news/index.html @@ -71,7 +71,7 @@ mkin - 1.0.3.9000 + 1.0.4.9000
@@ -141,10 +141,9 @@ Source: NEWS.md -
-

-mkin 1.0.3.9000 Unreleased -

+
+

+mkin 1.0.4.9000

General

@@ -159,29 +158,35 @@

  • Reintroduce the interface to the current development version of saemix, in particular:

  • ‘saemix_model’ and ‘saemix_data’: Helper functions to set up nonlinear mixed-effects models for mmkin row objects

  • ‘saem’: generic function to fit saemix models using ‘saemix_model’ and ‘saemix_data’, with a generator ‘saem.mmkin’, summary and plot methods

    • +

  • ‘mean_degparms’: New argument ‘test_log_parms’ that makes the function only consider log-transformed parameters where the untransformed parameters pass the t-test for a certain confidence level. This can be used to check more plausible parameters for ‘saem’

    • -
    +
    +

    +mkin 1.0.4 (Unreleased)

    +
      +

    • ‘plot.mixed.mmkin’: Reset graphical parameters on exit

      • +

    • All plotting functions setting graphical parameters: Use on.exit() for resetting graphical parameters

      • +
    +
    +

    -mkin 1.0.3 Unreleased -

    +mkin 1.0.3 (2021-02-15)
    • Review and update README, the ‘Introduction to mkin’ vignette and some of the help pages
    -
    +

    -mkin 1.0.2 Unreleased -

    +mkin 1.0.2 (Unreleased)
    • ‘mkinfit’: Keep model names stored in ‘mkinmod’ objects, avoiding their loss in ‘gmkin’
    -
    +

    -mkin 1.0.1 2021-02-10 -

    +mkin 1.0.1 (2021-02-10)

    • ‘confint.mmkin’, ‘nlme.mmkin’, ‘transform_odeparms’: Fix example code in dontrun sections that failed with current defaults

    • ‘logLik.mkinfit’: Improve example code to avoid warnings and show convenient syntax

      • @@ -190,10 +195,9 @@

    • Increase test tolerance for some parameter comparisons that also proved to be platform dependent

    -
    +

    -mkin 1.0.0 2021-02-03 -

    +mkin 1.0.0 (2021-02-03)

    General

    @@ -221,8 +225,7 @@

    -mkin 0.9.50.3 (2020-10-08) 2020-10-08 -

    +mkin 0.9.50.3 (2020-10-08)

    • ‘parms’: Add a method for mmkin objects

    • ‘mmkin’ and ‘confint(method = ’profile’): Use all cores detected by parallel::detectCores() per default

      • @@ -237,8 +240,7 @@

    -mkin 0.9.50.2 (2020-05-12) 2020-05-12 -

    +mkin 0.9.50.2 (2020-05-12)

    • Increase tolerance for a platform specific test results on the Solaris test machine on CRAN

    • Updates and corrections (using the spelling package) to the documentation

      • @@ -246,8 +248,7 @@

    -mkin 0.9.50.1 (2020-05-11) 2020-05-11 -

    +mkin 0.9.50.1 (2020-05-11)

    • Support SFORB with formation fractions

    • ‘mkinmod’: Make ‘use_of_ff’ = “max” the default

      • @@ -256,16 +257,14 @@

    -mkin 0.9.49.11 (2020-04-20) 2020-04-20 -

    +mkin 0.9.49.11 (2020-04-20)
    • Increase a test tolerance to make it pass on all CRAN check machines

    -mkin 0.9.49.10 (2020-04-18) 2020-04-18 -

    +mkin 0.9.49.10 (2020-04-18)

    • ‘nlme.mmkin’: An nlme method for mmkin row objects and an associated S3 class with print, plot, anova and endpoint methods

    • ‘mean_degparms, nlme_data, nlme_function’: Three new functions to facilitate building nlme models from mmkin row objects

      • @@ -277,8 +276,7 @@

    -mkin 0.9.49.9 (2020-03-31) 2020-03-31 -

    +mkin 0.9.49.9 (2020-03-31)

    • ‘mkinmod’: Use pkgbuild::has_compiler instead of Sys.which(‘gcc’), as the latter will often fail even if Rtools are installed

    • ‘mkinds’: Use roxygen for documenting fields and methods of this R6 class

      • @@ -286,8 +284,7 @@

    -mkin 0.9.49.8 (2020-01-09) 2020-01-09 -

    +mkin 0.9.49.8 (2020-01-09)

    • ‘aw’: Generic function for calculating Akaike weights, methods for mkinfit objects and mmkin columns

    • ‘loftest’: Add a lack-of-fit test

      • @@ -298,8 +295,7 @@

    -mkin 0.9.49.7 (2019-11-01) 2019-11-02 -

    +mkin 0.9.49.7 (2019-11-01)

    • Fix a bug introduced in 0.9.49.6 that occurred if the direct optimisation yielded a higher likelihood than the three-step optimisation in the d_3 algorithm, which caused the fitted parameters of the three-step optimisation to be returned instead of the parameters of the direct optimisation

    • Add a ‘nobs’ method for mkinfit objects, enabling the default ‘BIC’ method from the stats package. Also, add a ‘BIC’ method for mmkin column objects.

      • @@ -307,8 +303,7 @@

    -mkin 0.9.49.6 (2019-10-31) 2019-10-31 -

    +mkin 0.9.49.6 (2019-10-31)

    • Implement a likelihood ratio test as a method for ‘lrtest’ from the lmtest package

    • Add an ‘update’ method for mkinfit objects which remembers fitted parameters if appropriate

      • @@ -327,8 +322,7 @@

    -mkin 0.9.49.5 (2019-07-04) 2019-07-04 -

    +mkin 0.9.49.5 (2019-07-04)

    • Several algorithms for minimization of the negative log-likelihood for non-constant error models (two-component and variance by variable). In the case the error model is constant variance, least squares is used as this is more stable. The default algorithm ‘d_3’ tries direct minimization and a three-step procedure, and returns the model with the highest likelihood.

    • The argument ‘reweight.method’ to mkinfit and mmkin is now obsolete, use ‘error_model’ and ‘error_model_algorithm’ instead

      • @@ -346,8 +340,7 @@

    -mkin 0.9.48.1 (2019-03-04) 2019-03-04 -

    +mkin 0.9.48.1 (2019-03-04)

    • Add the function ‘logLik.mkinfit’ which makes it possible to calculate an AIC for mkinfit objects

    • Add the function ‘AIC.mmkin’ to make it easy to compare columns of mmkin objects

      • @@ -363,8 +356,7 @@

    -mkin 0.9.47.5 (2018-09-14) 2018-09-14 -

    +mkin 0.9.47.5 (2018-09-14)

    • Make the two-component error model stop in cases where it is inadequate to avoid nls crashes on windows

    • Move two vignettes to a location where they will not be built on CRAN (to avoid more NOTES from long execution times)

      • @@ -373,8 +365,7 @@

    -mkin 0.9.47.3 Unreleased -

    +mkin 0.9.47.3

    • ‘mkinfit’: Improve fitting the error model for reweight.method = ‘tc’. Add ‘manual’ to possible arguments for ‘weight’

    • Test that FOCUS_2006_C can be evaluated with DFOP and reweight.method = ‘tc’

      • @@ -382,8 +373,7 @@

    -mkin 0.9.47.2 (2018-07-19) 2018-07-19 -

    +mkin 0.9.47.2 (2018-07-19)

    • ‘sigma_twocomp’: Rename ‘sigma_rl’ to ‘sigma_twocomp’ as the Rocke and Lorenzato model assumes lognormal distribution for large y. Correct references to the Rocke and Lorenzato model accordingly.

    • ‘mkinfit’: Use 1.1 as starting value for N parameter of IORE models to obtain convergence in more difficult cases. Show parameter names when ‘trace_parms’ is ‘TRUE’.

      • @@ -391,8 +381,7 @@

    -mkin 0.9.47.1 (2018-02-06) 2018-02-06 -

    +mkin 0.9.47.1 (2018-02-06)

    • Skip some tests on CRAN and winbuilder to avoid timeouts

    • ‘test_data_from_UBA_2014’: Added this list of datasets containing experimental data used in the expertise from 2014

      • @@ -404,8 +393,7 @@

    -mkin 0.9.46.3 (2017-11-16) 2017-11-16 -

    +mkin 0.9.46.3 (2017-11-16)

    • README.md, vignettes/mkin.Rmd: URLs were updated

    • synthetic_data_for_UBA: Add the code used to generate the data in the interest of reproducibility

      • @@ -413,8 +401,7 @@

    -mkin 0.9.46.2 (2017-10-10) 2017-10-10 -

    +mkin 0.9.46.2 (2017-10-10)

    • Converted the vignette FOCUS_Z from tex/pdf to markdown/html

    • DESCRIPTION: Add ORCID

      • @@ -422,8 +409,7 @@

    -mkin 0.9.46.1 (2017-09-14) 2017-09-14 -

    +mkin 0.9.46.1 (2017-09-14)

    • plot.mkinfit: Fix scaling of residual plots for the case of separate plots for each observed variable

    • plot.mkinfit: Use all data points of the fitted curve for y axis scaling for the case of separate plots for each observed variable

      • @@ -432,16 +418,14 @@

    -mkin 0.9.46 (2017-07-24) 2017-07-29 -

    +mkin 0.9.46 (2017-07-24)
    • Remove test_FOMC_ill-defined.R as it is too platform dependent

    -mkin 0.9.45.2 (2017-07-24) 2017-07-22 -

    +mkin 0.9.45.2 (2017-07-24)

    • Rename twa to max_twa_parent to avoid conflict with twa from my pfm package

    • Update URLs in documentation

      • @@ -451,8 +435,7 @@

    -mkin 0.9.45.1 (2016-12-20) Unreleased -

    +mkin 0.9.45.1 (2016-12-20)

    New features

    @@ -463,8 +446,7 @@

    -mkin 0.9.45 (2016-12-08) 2016-12-08 -

    +mkin 0.9.45 (2016-12-08)

    Minor changes

    @@ -477,8 +459,7 @@

    -mkin 0.9.44 (2016-06-29) 2016-06-29 -

    +mkin 0.9.44 (2016-06-29)

    Bug fixes

    @@ -489,8 +470,7 @@

    -mkin 0.9.43 (2016-06-28) 2016-06-28 -

    +mkin 0.9.43 (2016-06-28)

    Major changes

    @@ -528,8 +508,7 @@

    -mkin 0.9.42 (2016-03-25) 2016-03-25 -

    +mkin 0.9.42 (2016-03-25)

    Major changes

    @@ -549,8 +528,7 @@

    -mkin 0.9-41 (2015-11-09) 2015-11-09 -

    +mkin 0.9-41 (2015-11-09)

    Minor changes

    @@ -572,8 +550,7 @@

    -mkin 0.9-40 (2015-07-21) 2015-07-21 -

    +mkin 0.9-40 (2015-07-21)

    Bug fixes

    @@ -593,8 +570,7 @@

    -mkin 0.9-39 (2015-06-26) 2015-06-26 -

    +mkin 0.9-39 (2015-06-26)

    Major changes

    @@ -614,8 +590,7 @@

    -mkin 0.9-38 (2015-06-24) 2015-06-23 -

    +mkin 0.9-38 (2015-06-24)

    Minor changes

    @@ -635,8 +610,7 @@

    -mkin 0.9-36 (2015-06-21) 2015-06-21 -

    +mkin 0.9-36 (2015-06-21)

    Major changes

    @@ -657,8 +631,7 @@

    -mkin 0.9-35 (2015-05-15) 2015-05-15 -

    +mkin 0.9-35 (2015-05-15)

    Major changes

    @@ -689,8 +662,7 @@

    -mkin 0.9-34 (2014-11-22) 2014-11-22 -

    +mkin 0.9-34 (2014-11-22)

    New features

    @@ -711,8 +683,7 @@

    -mkin 0.9-33 (2014-10-22) 2014-10-12 -

    +mkin 0.9-33 (2014-10-22)

    New features

    @@ -744,8 +715,7 @@

    -mkin 0.9-32 (2014-07-24) 2014-07-24 -

    +mkin 0.9-32 (2014-07-24)

    New features

    @@ -781,8 +751,7 @@

    -mkin 0.9-31 (2014-07-14) 2014-07-14 -

    +mkin 0.9-31 (2014-07-14)

    Bug fixes

    @@ -793,8 +762,7 @@

    -mkin 0.9-30 (2014-07-11) 2014-07-11 -

    +mkin 0.9-30 (2014-07-11)

    New features

    @@ -825,8 +793,7 @@

    -mkin 0.9-29 (2014-06-27) 2014-06-27 -

    +mkin 0.9-29 (2014-06-27)

    • R/mkinresplot.R: Make it possible to specify xlim

    • R/geometric_mean.R, man/geometric_mean.Rd: Add geometric mean function

      • @@ -835,8 +802,7 @@

    -mkin 0.9-28 (2014-05-20) 2014-05-20 -

    +mkin 0.9-28 (2014-05-20)

    • Do not backtransform confidence intervals for formation fractions if more than one compound is formed, as such parameters only define the pathways as a set

    • Add historical remarks and some background to the main package vignette

      • @@ -845,8 +811,7 @@

    -mkin 0.9-27 (2014-05-10) 2014-05-10 -

    +mkin 0.9-27 (2014-05-10)

    • Fork the GUI into a separate package gmkin

    • DESCRIPTION, NAMESPACE, TODO: Adapt and add copyright information

      • @@ -869,8 +834,7 @@

    -mkin 0.9-24 (2013-11-06) 2013-11-06 -

    +mkin 0.9-24 (2013-11-06)

    • Bugfix re-enabling the fixing of any combination of initial values for state variables

    • Default values for kinetic rate constants are not all 0.1 any more but are “salted” with a small increment to avoid numeric artefacts with the eigenvalue based solutions

      • @@ -879,8 +843,7 @@

    -mkin 0.9-22 (2013-10-26) 2013-10-26 -

    +mkin 0.9-22 (2013-10-26)

    • Get rid of the optimisation step in mkinerrmin - this was unnecessary. Thanks to KinGUII for the inspiration - actually this is equation 6-2 in FOCUS kinetics p. 91 that I had overlooked originally

    • Fix plot.mkinfit as it passed graphical arguments like main to the solver

      • diff --git a/docs/dev/pkgdown.yml b/docs/dev/pkgdown.yml index 4df60994..dbacd0ab 100644 --- a/docs/dev/pkgdown.yml +++ b/docs/dev/pkgdown.yml @@ -10,7 +10,7 @@ articles: web_only/NAFTA_examples: NAFTA_examples.html web_only/benchmarks: benchmarks.html web_only/compiled_models: compiled_models.html -last_built: 2021-02-15T16:08Z +last_built: 2021-03-09T16:32Z urls: reference: https://pkgdown.jrwb.de/mkin/reference article: https://pkgdown.jrwb.de/mkin/articles diff --git a/docs/dev/reference/Rplot001.png b/docs/dev/reference/Rplot001.png index 7f498242..17a35806 100644 Binary files a/docs/dev/reference/Rplot001.png and b/docs/dev/reference/Rplot001.png differ diff --git a/docs/dev/reference/Rplot002.png b/docs/dev/reference/Rplot002.png index 54c31a3f..a9a972e5 100644 Binary files a/docs/dev/reference/Rplot002.png and b/docs/dev/reference/Rplot002.png differ diff --git a/docs/dev/reference/Rplot003.png b/docs/dev/reference/Rplot003.png index 2b011ec1..d077f01c 100644 Binary files a/docs/dev/reference/Rplot003.png and b/docs/dev/reference/Rplot003.png differ diff --git a/docs/dev/reference/Rplot004.png b/docs/dev/reference/Rplot004.png index 98dd019e..ffcd2d96 100644 Binary files a/docs/dev/reference/Rplot004.png and b/docs/dev/reference/Rplot004.png differ diff --git a/docs/dev/reference/Rplot005.png b/docs/dev/reference/Rplot005.png index 8c91d61e..dfb5965b 100644 Binary files a/docs/dev/reference/Rplot005.png and b/docs/dev/reference/Rplot005.png differ diff --git a/docs/dev/reference/endpoints.html b/docs/dev/reference/endpoints.html index c9912f9c..63bec6a8 100644 --- a/docs/dev/reference/endpoints.html +++ b/docs/dev/reference/endpoints.html @@ -78,7 +78,7 @@ advantage that the SFORB model can also be used for metabolites." /> mkin - 1.0.3.9000 + 1.0.4.9000
    diff --git a/docs/dev/reference/index.html b/docs/dev/reference/index.html index 03a21517..5533a01f 100644 --- a/docs/dev/reference/index.html +++ b/docs/dev/reference/index.html @@ -71,7 +71,7 @@ mkin - 1.0.3.9000 + 1.0.4.9000
    diff --git a/docs/dev/reference/nlme-1.png b/docs/dev/reference/nlme-1.png index 728cc557..fd68ae43 100644 Binary files a/docs/dev/reference/nlme-1.png and b/docs/dev/reference/nlme-1.png differ diff --git a/docs/dev/reference/nlme-2.png b/docs/dev/reference/nlme-2.png index e8167455..853cae40 100644 Binary files a/docs/dev/reference/nlme-2.png and b/docs/dev/reference/nlme-2.png differ diff --git a/docs/dev/reference/nlme.html b/docs/dev/reference/nlme.html index b850eb3d..78d132e9 100644 --- a/docs/dev/reference/nlme.html +++ b/docs/dev/reference/nlme.html @@ -75,7 +75,7 @@ datasets. They are used internally by the nlme.mmkin() method." /> mkin - 1.0.3.9000 + 1.0.4.9000
    @@ -155,7 +155,7 @@ datasets. They are used internally by the nlme.m 
    nlme_function(object)
 
-mean_degparms(object, random = FALSE)
+mean_degparms(object, random = FALSE, test_log_parms = FALSE, conf.level = 0.6)
 
 nlme_data(object)
    @@ -170,6 +170,17 @@ datasets. They are used internally by the     nlme.m random

    Should a list with fixed and random effects be returned?

    + + test_log_parms +

    If TRUE, log parameters are only considered in +the mean calculations if their untransformed counterparts (most likely +rate constants) pass the t-test for significant difference from zero.

    + + + conf.level +

    Possibility to adjust the required confidence level +for parameter that are tested if requested by 'test_log_parms'.

    +

    Value

    @@ -211,7 +222,7 @@ nlme for the case of a single grouping variable ds.

    nlme_f <- nlme_function(f) # These assignments are necessary for these objects to be # visible to nlme and augPred when evaluation is done by -# pkgdown to generated the html docs. +# pkgdown to generate the html docs. assign("nlme_f", nlme_f, globalenv()) assign("grouped_data", grouped_data, globalenv()) @@ -226,28 +237,28 @@ nlme for the case of a single grouping variable ds.

    #> Model: value ~ nlme_f(name, time, parent_0, log_k_parent_sink) #> Data: grouped_data #> AIC BIC logLik -#> 300.6824 310.2426 -145.3412 +#> 298.2781 307.7372 -144.1391 #> #> Random effects: #> Formula: list(parent_0 ~ 1, log_k_parent_sink ~ 1) #> Level: ds #> Structure: Diagonal #> parent_0 log_k_parent_sink Residual -#> StdDev: 1.697361 0.6801209 3.666073 +#> StdDev: 0.937473 0.7098105 3.83543 #> #> Fixed effects: parent_0 + log_k_parent_sink ~ 1 #> Value Std.Error DF t-value p-value -#> parent_0 100.99378 1.3890416 46 72.70753 0 -#> log_k_parent_sink -3.07521 0.4018589 46 -7.65246 0 +#> parent_0 101.76838 1.1445443 45 88.91607 0 +#> log_k_parent_sink -3.05444 0.4195622 45 -7.28008 0 #> Correlation: #> prnt_0 -#> log_k_parent_sink 0.027 +#> log_k_parent_sink 0.034 #> #> Standardized Within-Group Residuals: -#> Min Q1 Med Q3 Max -#> -1.9942823 -0.5622565 0.1791579 0.7165038 2.0704781 +#> Min Q1 Med Q3 Max +#> -2.61693595 -0.21853231 0.05740682 0.57209372 3.04598764 #> -#> Number of Observations: 50 +#> Number of Observations: 49 #> Number of Groups: 3
    plot(augPred(m_nlme, level = 0:1), layout = c(3, 1))
    # augPred does not work on fits with more than one state # variable diff --git a/docs/dev/reference/nlme.mmkin-1.png b/docs/dev/reference/nlme.mmkin-1.png index 9186c135..90ede880 100644 Binary files a/docs/dev/reference/nlme.mmkin-1.png and b/docs/dev/reference/nlme.mmkin-1.png differ diff --git a/docs/dev/reference/nlme.mmkin-2.png b/docs/dev/reference/nlme.mmkin-2.png index d395fe02..0d140fd1 100644 Binary files a/docs/dev/reference/nlme.mmkin-2.png and b/docs/dev/reference/nlme.mmkin-2.png differ diff --git a/docs/dev/reference/nlme.mmkin-3.png b/docs/dev/reference/nlme.mmkin-3.png index 40518a59..8a60b52b 100644 Binary files a/docs/dev/reference/nlme.mmkin-3.png and b/docs/dev/reference/nlme.mmkin-3.png differ diff --git a/docs/dev/reference/nlme.mmkin.html b/docs/dev/reference/nlme.mmkin.html index 925cf7cf..f308d8b7 100644 --- a/docs/dev/reference/nlme.mmkin.html +++ b/docs/dev/reference/nlme.mmkin.html @@ -74,7 +74,7 @@ have been obtained by fitting the same model to a list of datasets." /> mkin - 1.0.3.9000 + 1.0.4.9000
    diff --git a/docs/dev/reference/plot.mixed.mmkin-1.png b/docs/dev/reference/plot.mixed.mmkin-1.png index 9c9a2211..2224d96e 100644 Binary files a/docs/dev/reference/plot.mixed.mmkin-1.png and b/docs/dev/reference/plot.mixed.mmkin-1.png differ diff --git a/docs/dev/reference/plot.mixed.mmkin-2.png b/docs/dev/reference/plot.mixed.mmkin-2.png index 0f66ff0f..28168495 100644 Binary files a/docs/dev/reference/plot.mixed.mmkin-2.png and b/docs/dev/reference/plot.mixed.mmkin-2.png differ diff --git a/docs/dev/reference/plot.mixed.mmkin-3.png b/docs/dev/reference/plot.mixed.mmkin-3.png index 34212f1c..d18275dd 100644 Binary files a/docs/dev/reference/plot.mixed.mmkin-3.png and b/docs/dev/reference/plot.mixed.mmkin-3.png differ diff --git a/docs/dev/reference/plot.mixed.mmkin-4.png b/docs/dev/reference/plot.mixed.mmkin-4.png index c1450d24..2fd52425 100644 Binary files a/docs/dev/reference/plot.mixed.mmkin-4.png and b/docs/dev/reference/plot.mixed.mmkin-4.png differ diff --git a/docs/dev/reference/plot.mixed.mmkin.html b/docs/dev/reference/plot.mixed.mmkin.html index 630e95a3..36796580 100644 --- a/docs/dev/reference/plot.mixed.mmkin.html +++ b/docs/dev/reference/plot.mixed.mmkin.html @@ -72,7 +72,7 @@ mkin - 1.0.3.9000 + 1.0.4.9000
    @@ -283,10 +283,10 @@ corresponding model prediction lines for the different datasets.

    f_saem <- saem(f, transformations = "saemix")
    #> Running main SAEM algorithm -#> [1] "Mon Feb 15 17:12:17 2021" +#> [1] "Tue Mar 9 17:34:35 2021" #> .... #> Minimisation finished -#> [1] "Mon Feb 15 17:12:24 2021"
    plot(f_saem) +#> [1] "Tue Mar 9 17:34:42 2021"
    plot(f_saem)
    # We can overlay the two variants if we generate predictions pred_nlme <- mkinpredict(dfop_sfo, diff --git a/docs/dev/reference/saem-1.png b/docs/dev/reference/saem-1.png index 2df248bb..0da31388 100644 Binary files a/docs/dev/reference/saem-1.png and b/docs/dev/reference/saem-1.png differ diff --git a/docs/dev/reference/saem-2.png b/docs/dev/reference/saem-2.png index d4a2c1be..010950ba 100644 Binary files a/docs/dev/reference/saem-2.png and b/docs/dev/reference/saem-2.png differ diff --git a/docs/dev/reference/saem-3.png b/docs/dev/reference/saem-3.png index 4474b1f1..829f22bf 100644 Binary files a/docs/dev/reference/saem-3.png and b/docs/dev/reference/saem-3.png differ diff --git a/docs/dev/reference/saem-4.png b/docs/dev/reference/saem-4.png index bf24d6b0..4e976fa2 100644 Binary files a/docs/dev/reference/saem-4.png and b/docs/dev/reference/saem-4.png differ diff --git a/docs/dev/reference/saem-5.png b/docs/dev/reference/saem-5.png index 27ed3f8f..f50969b4 100644 Binary files a/docs/dev/reference/saem-5.png and b/docs/dev/reference/saem-5.png differ diff --git a/docs/dev/reference/saem.html b/docs/dev/reference/saem.html index bdb1226e..23102df3 100644 --- a/docs/dev/reference/saem.html +++ b/docs/dev/reference/saem.html @@ -74,7 +74,7 @@ Expectation Maximisation algorithm (SAEM)." /> mkin - 1.0.3.9000 + 1.0.4.9000
    @@ -158,9 +158,12 @@ Expectation Maximisation algorithm (SAEM).

    object, transformations = c("mkin", "saemix"), degparms_start = numeric(), + test_log_parms = FALSE, + conf.level = 0.6, solution_type = "auto", control = list(displayProgress = FALSE, print = FALSE, save = FALSE, save.graphs = FALSE), + fail_with_errors = TRUE, verbose = FALSE, quiet = FALSE, ... @@ -174,6 +177,7 @@ Expectation Maximisation algorithm (SAEM).

    solution_type = "auto", transformations = c("mkin", "saemix"), degparms_start = numeric(), + test_log_parms = FALSE, verbose = FALSE, ... ) @@ -204,6 +208,18 @@ SFO or DFOP is used for the parent and there is either no metabolite or one.

    degparms_start

    Parameter values given as a named numeric vector will be used to override the starting values obtained from the 'mmkin' object.

    + + + 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.

    + + + conf.level +

    Possibility to adjust the required confidence level +for parameter that are tested if requested by 'test_log_parms'.

    solution_type @@ -214,6 +230,11 @@ automatic choice is not desired

    control

    Passed to saemix::saemix

    + + fail_with_errors +

    Should a failure to compute standard errors +from the inverse of the Fisher Information Matrix be a failure?

    + verbose

    Should we print information about created objects of @@ -261,33 +282,36 @@ using mmkin.

    state.ini = c(parent = 100), fixed_initials = "parent", quiet = TRUE) f_saem_p0_fixed <- saem(f_mmkin_parent_p0_fixed)
    #> Running main SAEM algorithm -#> [1] "Mon Feb 15 17:12:32 2021" +#> [1] "Tue Mar 9 17:34:44 2021" #> .... #> Minimisation finished -#> [1] "Mon Feb 15 17:12:34 2021"
    +#> [1] "Tue Mar 9 17:34:45 2021"
    f_mmkin_parent <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE) f_saem_sfo <- saem(f_mmkin_parent["SFO", ])
    #> Running main SAEM algorithm -#> [1] "Mon Feb 15 17:12:35 2021" +#> [1] "Tue Mar 9 17:34:46 2021" #> .... #> Minimisation finished -#> [1] "Mon Feb 15 17:12:36 2021"
    f_saem_fomc <- saem(f_mmkin_parent["FOMC", ]) +#> [1] "Tue Mar 9 17:34:48 2021"
    f_saem_fomc <- saem(f_mmkin_parent["FOMC", ])
    #> Running main SAEM algorithm -#> [1] "Mon Feb 15 17:12:36 2021" +#> [1] "Tue Mar 9 17:34:48 2021" #> .... #> Minimisation finished -#> [1] "Mon Feb 15 17:12:38 2021"
    f_saem_dfop <- saem(f_mmkin_parent["DFOP", ]) +#> [1] "Tue Mar 9 17:34:50 2021"
    f_saem_dfop <- saem(f_mmkin_parent["DFOP", ])
    #> Running main SAEM algorithm -#> [1] "Mon Feb 15 17:12:39 2021" +#> [1] "Tue Mar 9 17:34:51 2021" #> .... #> Minimisation finished -#> [1] "Mon Feb 15 17:12:42 2021"
    +#> [1] "Tue Mar 9 17:34:53 2021"
    # The returned saem.mmkin object contains an SaemixObject, therefore we can use # functions from saemix library(saemix)
    #> Package saemix, version 3.1.9000 -#> please direct bugs, questions and feedback to emmanuelle.comets@inserm.fr
    compare.saemix(list(f_saem_sfo$so, f_saem_fomc$so, f_saem_dfop$so)) -
    #> Error in compare.saemix(list(f_saem_sfo$so, f_saem_fomc$so, f_saem_dfop$so)): 'compare.saemix' requires at least two models.
    plot(f_saem_fomc$so, plot.type = "convergence") +#> please direct bugs, questions and feedback to emmanuelle.comets@inserm.fr
    compare.saemix(f_saem_sfo$so, f_saem_fomc$so, f_saem_dfop$so) +
    #> Likelihoods calculated by importance sampling
    #> AIC BIC +#> 1 624.2484 622.2956 +#> 2 467.7096 464.9757 +#> 3 495.4373 491.9222
    plot(f_saem_fomc$so, plot.type = "convergence")
    #> Plotting convergence plots
    plot(f_saem_fomc$so, plot.type = "individual.fit")
    #> Plotting individual fits
    plot(f_saem_fomc$so, plot.type = "npde")
    #> Simulating data using nsim = 1000 simulated datasets @@ -324,11 +348,13 @@ using mmkin.

    f_mmkin_parent_tc <- update(f_mmkin_parent, error_model = "tc") f_saem_fomc_tc <- saem(f_mmkin_parent_tc["FOMC", ])
    #> Running main SAEM algorithm -#> [1] "Mon Feb 15 17:12:44 2021" +#> [1] "Tue Mar 9 17:34:55 2021" #> .... #> Minimisation finished -#> [1] "Mon Feb 15 17:12:49 2021"
    compare.saemix(list(f_saem_fomc$so, f_saem_fomc_tc$so)) -
    #> Error in compare.saemix(list(f_saem_fomc$so, f_saem_fomc_tc$so)): 'compare.saemix' requires at least two models.
    +#> [1] "Tue Mar 9 17:35:00 2021"
    compare.saemix(f_saem_fomc$so, f_saem_fomc_tc$so) +
    #> Likelihoods calculated by importance sampling
    #> AIC BIC +#> 1 467.7096 464.9757 +#> 2 469.6831 466.5586
    sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"), A1 = mkinsub("SFO"))
    #> Temporary DLL for differentials generated and loaded
    fomc_sfo <- mkinmod(parent = mkinsub("FOMC", "A1"), @@ -346,15 +372,15 @@ using mmkin.

    # four minutes f_saem_sfo_sfo <- saem(f_mmkin["SFO-SFO", ])
    #> Running main SAEM algorithm -#> [1] "Mon Feb 15 17:12:51 2021" +#> [1] "Tue Mar 9 17:35:02 2021" #> .... #> Minimisation finished -#> [1] "Mon Feb 15 17:12:56 2021"
    f_saem_dfop_sfo <- saem(f_mmkin["DFOP-SFO", ]) +#> [1] "Tue Mar 9 17:35:07 2021"
    f_saem_dfop_sfo <- saem(f_mmkin["DFOP-SFO", ])
    #> Running main SAEM algorithm -#> [1] "Mon Feb 15 17:12:56 2021" +#> [1] "Tue Mar 9 17:35:07 2021" #> .... #> Minimisation finished -#> [1] "Mon Feb 15 17:13:05 2021"
    # We can use print, plot and summary methods to check the results +#> [1] "Tue Mar 9 17:35:15 2021"
    # We can use print, plot and summary methods to check the results print(f_saem_dfop_sfo)
    #> Kinetic nonlinear mixed-effects model fit by SAEM #> Structural model: @@ -395,10 +421,10 @@ using mmkin.

    #> SD.g_qlogis 0.44771 -0.86417 1.7596
    plot(f_saem_dfop_sfo)
    summary(f_saem_dfop_sfo, data = TRUE)
    #> saemix version used for fitting: 3.1.9000 -#> mkin version used for pre-fitting: 1.0.3.9000 -#> R version used for fitting: 4.0.3 -#> Date of fit: Mon Feb 15 17:13:05 2021 -#> Date of summary: Mon Feb 15 17:13:06 2021 +#> mkin version used for pre-fitting: 1.0.4.9000 +#> R version used for fitting: 4.0.4 +#> Date of fit: Tue Mar 9 17:35:16 2021 +#> Date of summary: Tue Mar 9 17:35:16 2021 #> #> Equations: #> d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * @@ -413,7 +439,7 @@ using mmkin.

    #> #> Model predictions using solution type analytical #> -#> Fitted in 8.985 s using 300, 100 iterations +#> Fitted in 8.668 s using 300, 100 iterations #> #> Variance model: Constant variance #> diff --git a/docs/dev/reference/summary.saem.mmkin.html b/docs/dev/reference/summary.saem.mmkin.html index 0d661ee9..1166abb1 100644 --- a/docs/dev/reference/summary.saem.mmkin.html +++ b/docs/dev/reference/summary.saem.mmkin.html @@ -76,7 +76,7 @@ endpoints such as formation fractions and DT50 values. Optionally mkin - 1.0.3.9000 + 1.0.4.9000
    @@ -260,15 +260,15 @@ saemix authors for the parts inherited from saemix.

    quiet = TRUE, error_model = "tc", cores = 5) f_saem_dfop_sfo <- saem(f_mmkin_dfop_sfo)
    #> Running main SAEM algorithm -#> [1] "Mon Feb 15 17:13:15 2021" +#> [1] "Tue Mar 9 17:35:19 2021" #> .... #> Minimisation finished -#> [1] "Mon Feb 15 17:13:26 2021"
    summary(f_saem_dfop_sfo, data = TRUE) +#> [1] "Tue Mar 9 17:35:30 2021"
    summary(f_saem_dfop_sfo, data = TRUE)
    #> saemix version used for fitting: 3.1.9000 -#> mkin version used for pre-fitting: 1.0.3.9000 -#> R version used for fitting: 4.0.3 -#> Date of fit: Mon Feb 15 17:13:27 2021 -#> Date of summary: Mon Feb 15 17:13:27 2021 +#> mkin version used for pre-fitting: 1.0.4.9000 +#> R version used for fitting: 4.0.4 +#> Date of fit: Tue Mar 9 17:35:31 2021 +#> Date of summary: Tue Mar 9 17:35:31 2021 #> #> Equations: #> d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * @@ -283,7 +283,7 @@ saemix authors for the parts inherited from saemix.

    #> #> Model predictions using solution type analytical #> -#> Fitted in 12.204 s using 300, 100 iterations +#> Fitted in 12.058 s using 300, 100 iterations #> #> Variance model: Two-component variance function #> @@ -300,231 +300,231 @@ saemix authors for the parts inherited from saemix.

    #> #> Likelihood computed by importance sampling #> AIC BIC logLik -#> 829.3 823.9 -400.7 +#> 825.9 820.4 -398.9 #> #> Optimised parameters: -#> est. lower upper -#> parent_0 101.29457 97.855 104.7344 -#> log_k_m1 -4.06337 -4.182 -3.9445 -#> f_parent_qlogis -0.94546 -1.307 -0.5841 -#> log_k1 -2.98794 -3.844 -2.1321 -#> log_k2 -3.47891 -4.253 -2.7050 -#> g_qlogis -0.03211 -1.157 1.0931 +#> est. lower upper +#> parent_0 101.118986 97.368 104.8695 +#> log_k_m1 -4.057591 -4.177 -3.9379 +#> f_parent_qlogis -0.933087 -1.290 -0.5763 +#> log_k1 -2.945520 -3.833 -2.0576 +#> log_k2 -3.531954 -4.310 -2.7542 +#> g_qlogis -0.009584 -1.688 1.6687 #> #> Correlation: #> prnt_0 lg_k_1 f_prn_ log_k1 log_k2 -#> log_k_m1 -0.202 -#> f_parent_qlogis -0.145 0.195 -#> log_k1 0.094 -0.099 -0.049 -#> log_k2 -0.042 0.056 0.024 -0.097 -#> g_qlogis -0.005 0.000 0.007 -0.160 -0.113 +#> log_k_m1 -0.198 +#> f_parent_qlogis -0.153 0.184 +#> log_k1 0.080 -0.077 -0.045 +#> log_k2 0.005 0.008 -0.003 -0.019 +#> g_qlogis -0.059 0.048 0.041 -0.334 -0.253 #> #> Random effects: #> est. lower upper -#> SD.parent_0 2.70085 -0.64980 6.0515 -#> SD.log_k_m1 0.08408 -0.04023 0.2084 -#> SD.f_parent_qlogis 0.39215 0.13695 0.6473 -#> SD.log_k1 0.89280 0.27466 1.5109 -#> SD.log_k2 0.82387 0.26388 1.3838 -#> SD.g_qlogis 0.36468 -0.86978 1.5991 +#> SD.parent_0 2.97797 -0.62927 6.5852 +#> SD.log_k_m1 0.09235 -0.02448 0.2092 +#> SD.f_parent_qlogis 0.38712 0.13469 0.6396 +#> SD.log_k1 0.88671 0.27052 1.5029 +#> SD.log_k2 0.80497 0.25587 1.3541 +#> SD.g_qlogis 0.36812 -3.56188 4.2981 #> #> Variance model: #> est. lower upper -#> a.1 0.65724 0.49361 0.82086 -#> b.1 0.06434 0.05034 0.07835 +#> a.1 0.85879 0.68143 1.03615 +#> b.1 0.07787 0.06288 0.09286 #> #> Backtransformed parameters: #> est. lower upper -#> parent_0 101.29457 97.85477 104.73437 -#> k_m1 0.01719 0.01526 0.01936 -#> f_parent_to_m1 0.27980 0.21302 0.35798 -#> k1 0.05039 0.02141 0.11859 -#> k2 0.03084 0.01422 0.06687 -#> g 0.49197 0.23916 0.74896 +#> parent_0 101.11899 97.36850 104.86947 +#> k_m1 0.01729 0.01534 0.01949 +#> f_parent_to_m1 0.28230 0.21587 0.35979 +#> k1 0.05257 0.02163 0.12776 +#> k2 0.02925 0.01344 0.06366 +#> g 0.49760 0.15606 0.84140 #> #> Resulting formation fractions: #> ff -#> parent_m1 0.2798 -#> parent_sink 0.7202 +#> parent_m1 0.2823 +#> parent_sink 0.7177 #> #> Estimated disappearance times: #> DT50 DT90 DT50back DT50_k1 DT50_k2 -#> parent 17.49 61.05 18.38 13.76 22.47 -#> m1 40.32 133.94 NA NA NA +#> parent 17.47 62.31 18.76 13.18 23.7 +#> m1 40.09 133.17 NA NA NA #> #> Data: -#> ds name time observed predicted residual std standardized -#> ds 1 parent 0 89.8 9.878e+01 8.98039 6.3899 1.40541 -#> ds 1 parent 0 104.1 9.878e+01 -5.31961 6.3899 -0.83251 -#> ds 1 parent 1 88.7 9.422e+01 5.52084 6.0981 0.90533 -#> ds 1 parent 1 95.5 9.422e+01 -1.27916 6.0981 -0.20976 -#> ds 1 parent 3 81.8 8.587e+01 4.06752 5.5641 0.73103 -#> ds 1 parent 3 94.5 8.587e+01 -8.63248 5.5641 -1.55147 -#> ds 1 parent 7 71.5 7.180e+01 0.29615 4.6662 0.06347 -#> ds 1 parent 7 70.3 7.180e+01 1.49615 4.6662 0.32063 -#> ds 1 parent 14 54.2 5.360e+01 -0.59602 3.5112 -0.16975 -#> ds 1 parent 14 49.6 5.360e+01 4.00398 3.5112 1.14035 -#> ds 1 parent 28 31.5 3.213e+01 0.62529 2.1691 0.28828 -#> ds 1 parent 28 28.8 3.213e+01 3.32529 2.1691 1.53306 -#> ds 1 parent 60 12.1 1.271e+01 0.60718 1.0490 0.57879 -#> ds 1 parent 60 13.6 1.271e+01 -0.89282 1.0490 -0.85108 -#> ds 1 parent 90 6.2 6.080e+00 -0.12020 0.7649 -0.15716 -#> ds 1 parent 90 8.3 6.080e+00 -2.22020 0.7649 -2.90279 -#> ds 1 parent 120 2.2 3.011e+00 0.81059 0.6852 1.18302 -#> ds 1 parent 120 2.4 3.011e+00 0.61059 0.6852 0.89113 -#> ds 1 m1 1 0.3 1.131e+00 0.83071 0.6613 1.25628 -#> ds 1 m1 1 0.2 1.131e+00 0.93071 0.6613 1.40750 -#> ds 1 m1 3 2.2 3.147e+00 0.94691 0.6877 1.37688 -#> ds 1 m1 3 3.0 3.147e+00 0.14691 0.6877 0.21361 -#> ds 1 m1 7 6.5 6.341e+00 -0.15949 0.7736 -0.20618 -#> ds 1 m1 7 5.0 6.341e+00 1.34051 0.7736 1.73290 -#> ds 1 m1 14 10.2 9.910e+00 -0.28991 0.9157 -0.31659 -#> ds 1 m1 14 9.5 9.910e+00 0.41009 0.9157 0.44783 -#> ds 1 m1 28 12.2 1.255e+01 0.34690 1.0410 0.33323 -#> ds 1 m1 28 13.4 1.255e+01 -0.85310 1.0410 -0.81949 -#> ds 1 m1 60 11.8 1.087e+01 -0.92713 0.9599 -0.96586 -#> ds 1 m1 60 13.2 1.087e+01 -2.32713 0.9599 -2.42434 -#> ds 1 m1 90 6.6 7.813e+00 1.21254 0.8274 1.46541 -#> ds 1 m1 90 9.3 7.813e+00 -1.48746 0.8274 -1.79766 -#> ds 1 m1 120 3.5 5.295e+00 1.79489 0.7403 2.42457 -#> ds 1 m1 120 5.4 5.295e+00 -0.10511 0.7403 -0.14198 -#> ds 2 parent 0 118.0 1.074e+02 -10.63436 6.9396 -1.53242 -#> ds 2 parent 0 99.8 1.074e+02 7.56564 6.9396 1.09021 -#> ds 2 parent 1 90.2 1.012e+02 10.96486 6.5425 1.67594 -#> ds 2 parent 1 94.6 1.012e+02 6.56486 6.5425 1.00342 -#> ds 2 parent 3 96.1 9.054e+01 -5.56014 5.8627 -0.94839 -#> ds 2 parent 3 78.4 9.054e+01 12.13986 5.8627 2.07069 -#> ds 2 parent 7 77.9 7.468e+01 -3.21805 4.8501 -0.66350 -#> ds 2 parent 7 77.7 7.468e+01 -3.01805 4.8501 -0.62226 -#> ds 2 parent 14 56.0 5.748e+01 1.47774 3.7563 0.39340 -#> ds 2 parent 14 54.7 5.748e+01 2.77774 3.7563 0.73948 -#> ds 2 parent 28 36.6 3.996e+01 3.36317 2.6541 1.26717 -#> ds 2 parent 28 36.8 3.996e+01 3.16317 2.6541 1.19182 -#> ds 2 parent 60 22.1 2.132e+01 -0.78225 1.5210 -0.51430 -#> ds 2 parent 60 24.7 2.132e+01 -3.38225 1.5210 -2.22369 -#> ds 2 parent 90 12.4 1.215e+01 -0.25010 1.0213 -0.24487 -#> ds 2 parent 90 10.8 1.215e+01 1.34990 1.0213 1.32169 -#> ds 2 parent 120 6.8 6.931e+00 0.13105 0.7943 0.16500 -#> ds 2 parent 120 7.9 6.931e+00 -0.96895 0.7943 -1.21994 -#> ds 2 m1 1 1.3 1.519e+00 0.21924 0.6645 0.32995 -#> ds 2 m1 3 3.7 4.049e+00 0.34891 0.7070 0.49351 -#> ds 2 m1 3 4.7 4.049e+00 -0.65109 0.7070 -0.92094 -#> ds 2 m1 7 8.1 7.565e+00 -0.53526 0.8179 -0.65448 -#> ds 2 m1 7 7.9 7.565e+00 -0.33526 0.8179 -0.40993 -#> ds 2 m1 14 10.1 1.071e+01 0.60614 0.9521 0.63663 -#> ds 2 m1 14 10.3 1.071e+01 0.40614 0.9521 0.42657 -#> ds 2 m1 28 10.7 1.224e+01 1.54440 1.0260 1.50526 -#> ds 2 m1 28 12.2 1.224e+01 0.04440 1.0260 0.04327 -#> ds 2 m1 60 10.7 1.056e+01 -0.14005 0.9453 -0.14815 -#> ds 2 m1 60 12.5 1.056e+01 -1.94005 0.9453 -2.05226 -#> ds 2 m1 90 9.1 8.089e+00 -1.01088 0.8384 -1.20577 -#> ds 2 m1 90 7.4 8.089e+00 0.68912 0.8384 0.82197 -#> ds 2 m1 120 6.1 5.855e+00 -0.24463 0.7576 -0.32292 -#> ds 2 m1 120 4.5 5.855e+00 1.35537 0.7576 1.78911 -#> ds 3 parent 0 106.2 1.095e+02 3.30335 7.0765 0.46680 -#> ds 3 parent 0 106.9 1.095e+02 2.60335 7.0765 0.36788 -#> ds 3 parent 1 107.4 9.939e+01 -8.01282 6.4287 -1.24641 -#> ds 3 parent 1 96.1 9.939e+01 3.28718 6.4287 0.51133 -#> ds 3 parent 3 79.4 8.365e+01 4.24698 5.4222 0.78326 -#> ds 3 parent 3 82.6 8.365e+01 1.04698 5.4222 0.19309 -#> ds 3 parent 7 63.9 6.405e+01 0.14704 4.1732 0.03523 -#> ds 3 parent 7 62.4 6.405e+01 1.64704 4.1732 0.39467 -#> ds 3 parent 14 51.0 4.795e+01 -3.04985 3.1546 -0.96681 -#> ds 3 parent 14 47.1 4.795e+01 0.85015 3.1546 0.26950 -#> ds 3 parent 28 36.1 3.501e+01 -1.09227 2.3465 -0.46549 -#> ds 3 parent 28 36.6 3.501e+01 -1.59227 2.3465 -0.67858 -#> ds 3 parent 60 20.1 2.012e+01 0.02116 1.4520 0.01457 -#> ds 3 parent 60 19.8 2.012e+01 0.32116 1.4520 0.22119 -#> ds 3 parent 90 11.3 1.206e+01 0.76096 1.0170 0.74826 -#> ds 3 parent 90 10.7 1.206e+01 1.36096 1.0170 1.33825 -#> ds 3 parent 120 8.2 7.230e+00 -0.97022 0.8052 -1.20493 -#> ds 3 parent 120 7.3 7.230e+00 -0.07022 0.8052 -0.08721 -#> ds 3 m1 0 0.8 -5.684e-13 -0.80000 0.6572 -1.21722 -#> ds 3 m1 1 1.8 2.045e+00 0.24538 0.6703 0.36608 -#> ds 3 m1 1 2.3 2.045e+00 -0.25462 0.6703 -0.37987 -#> ds 3 m1 3 4.2 5.136e+00 0.93594 0.7356 1.27228 -#> ds 3 m1 3 4.1 5.136e+00 1.03594 0.7356 1.40822 -#> ds 3 m1 7 6.8 8.674e+00 1.87438 0.8623 2.17381 -#> ds 3 m1 7 10.1 8.674e+00 -1.42562 0.8623 -1.65335 -#> ds 3 m1 14 11.4 1.083e+01 -0.56746 0.9580 -0.59233 -#> ds 3 m1 14 12.8 1.083e+01 -1.96746 0.9580 -2.05369 -#> ds 3 m1 28 11.5 1.098e+01 -0.51762 0.9651 -0.53637 -#> ds 3 m1 28 10.6 1.098e+01 0.38238 0.9651 0.39623 -#> ds 3 m1 60 7.5 8.889e+00 1.38911 0.8713 1.59436 -#> ds 3 m1 60 8.6 8.889e+00 0.28911 0.8713 0.33183 -#> ds 3 m1 90 7.3 6.774e+00 -0.52608 0.7886 -0.66708 -#> ds 3 m1 90 8.1 6.774e+00 -1.32608 0.7886 -1.68150 -#> ds 3 m1 120 5.3 4.954e+00 -0.34584 0.7305 -0.47345 -#> ds 3 m1 120 3.8 4.954e+00 1.15416 0.7305 1.58004 -#> ds 4 parent 0 104.7 9.957e+01 -5.13169 6.4403 -0.79681 -#> ds 4 parent 0 88.3 9.957e+01 11.26831 6.4403 1.74966 -#> ds 4 parent 1 94.2 9.644e+01 2.23888 6.2400 0.35879 -#> ds 4 parent 1 94.6 9.644e+01 1.83888 6.2400 0.29469 -#> ds 4 parent 3 78.1 9.054e+01 12.43946 5.8627 2.12180 -#> ds 4 parent 3 96.5 9.054e+01 -5.96054 5.8627 -1.01669 -#> ds 4 parent 7 76.2 8.004e+01 3.83771 5.1918 0.73919 -#> ds 4 parent 7 77.8 8.004e+01 2.23771 5.1918 0.43101 -#> ds 4 parent 14 70.8 6.511e+01 -5.69246 4.2406 -1.34238 -#> ds 4 parent 14 67.3 6.511e+01 -2.19246 4.2406 -0.51702 -#> ds 4 parent 28 43.1 4.454e+01 1.43744 2.9401 0.48890 -#> ds 4 parent 28 45.1 4.454e+01 -0.56256 2.9401 -0.19134 -#> ds 4 parent 60 21.3 2.132e+01 0.02005 1.5211 0.01318 -#> ds 4 parent 60 23.5 2.132e+01 -2.17995 1.5211 -1.43310 -#> ds 4 parent 90 11.8 1.182e+01 0.02167 1.0053 0.02156 -#> ds 4 parent 90 12.1 1.182e+01 -0.27833 1.0053 -0.27687 -#> ds 4 parent 120 7.0 6.852e+00 -0.14780 0.7914 -0.18675 -#> ds 4 parent 120 6.2 6.852e+00 0.65220 0.7914 0.82408 -#> ds 4 m1 0 1.6 -5.684e-14 -1.60000 0.6572 -2.43444 -#> ds 4 m1 1 0.9 6.918e-01 -0.20821 0.6587 -0.31607 -#> ds 4 m1 3 3.7 1.959e+00 -1.74131 0.6692 -2.60204 -#> ds 4 m1 3 2.0 1.959e+00 -0.04131 0.6692 -0.06173 -#> ds 4 m1 7 3.6 4.076e+00 0.47590 0.7076 0.67252 -#> ds 4 m1 7 3.8 4.076e+00 0.27590 0.7076 0.38989 -#> ds 4 m1 14 7.1 6.698e+00 -0.40189 0.7859 -0.51135 -#> ds 4 m1 14 6.6 6.698e+00 0.09811 0.7859 0.12483 -#> ds 4 m1 28 9.5 9.175e+00 -0.32492 0.8835 -0.36779 -#> ds 4 m1 28 9.3 9.175e+00 -0.12492 0.8835 -0.14141 -#> ds 4 m1 60 8.3 8.818e+00 0.51810 0.8683 0.59671 -#> ds 4 m1 60 9.0 8.818e+00 -0.18190 0.8683 -0.20949 -#> ds 4 m1 90 6.6 6.645e+00 0.04480 0.7841 0.05713 -#> ds 4 m1 90 7.7 6.645e+00 -1.05520 0.7841 -1.34581 -#> ds 4 m1 120 3.7 4.648e+00 0.94805 0.7221 1.31293 -#> ds 4 m1 120 3.5 4.648e+00 1.14805 0.7221 1.58991 -#> ds 5 parent 0 110.4 1.026e+02 -7.81752 6.6333 -1.17853 -#> ds 5 parent 0 112.1 1.026e+02 -9.51752 6.6333 -1.43482 -#> ds 5 parent 1 93.5 9.560e+01 2.10274 6.1865 0.33989 -#> ds 5 parent 1 91.0 9.560e+01 4.60274 6.1865 0.74399 -#> ds 5 parent 3 71.0 8.356e+01 12.55799 5.4165 2.31846 -#> ds 5 parent 3 89.7 8.356e+01 -6.14201 5.4165 -1.13394 -#> ds 5 parent 7 60.4 6.550e+01 5.09732 4.2653 1.19506 -#> ds 5 parent 7 59.1 6.550e+01 6.39732 4.2653 1.49984 -#> ds 5 parent 14 56.5 4.641e+01 -10.09145 3.0576 -3.30044 -#> ds 5 parent 14 47.0 4.641e+01 -0.59145 3.0576 -0.19344 -#> ds 5 parent 28 30.2 2.982e+01 -0.37647 2.0284 -0.18560 -#> ds 5 parent 28 23.9 2.982e+01 5.92353 2.0284 2.92028 -#> ds 5 parent 60 17.0 1.754e+01 0.53981 1.3060 0.41332 -#> ds 5 parent 60 18.7 1.754e+01 -1.16019 1.3060 -0.88834 -#> ds 5 parent 90 11.3 1.175e+01 0.45050 1.0018 0.44969 -#> ds 5 parent 90 11.9 1.175e+01 -0.14950 1.0018 -0.14923 -#> ds 5 parent 120 9.0 7.915e+00 -1.08476 0.8315 -1.30462 -#> ds 5 parent 120 8.1 7.915e+00 -0.18476 0.8315 -0.22220 -#> ds 5 m1 0 0.7 0.000e+00 -0.70000 0.6572 -1.06507 -#> ds 5 m1 1 3.0 3.062e+00 0.06170 0.6861 0.08992 -#> ds 5 m1 1 2.6 3.062e+00 0.46170 0.6861 0.67290 -#> ds 5 m1 3 5.1 8.209e+00 3.10938 0.8432 3.68760 -#> ds 5 m1 3 7.5 8.209e+00 0.70938 0.8432 0.84130 -#> ds 5 m1 7 16.5 1.544e+01 -1.05567 1.1914 -0.88605 -#> ds 5 m1 7 19.0 1.544e+01 -3.55567 1.1914 -2.98436 -#> ds 5 m1 14 22.9 2.181e+01 -1.08765 1.5498 -0.70181 -#> ds 5 m1 14 23.2 2.181e+01 -1.38765 1.5498 -0.89539 -#> ds 5 m1 28 22.2 2.404e+01 1.83624 1.6805 1.09270 -#> ds 5 m1 28 24.4 2.404e+01 -0.36376 1.6805 -0.21647 -#> ds 5 m1 60 15.5 1.875e+01 3.25390 1.3741 2.36805 -#> ds 5 m1 60 19.8 1.875e+01 -1.04610 1.3741 -0.76131 -#> ds 5 m1 90 14.9 1.380e+01 -1.09507 1.1050 -0.99102 -#> ds 5 m1 90 14.2 1.380e+01 -0.39507 1.1050 -0.35753 -#> ds 5 m1 120 10.9 1.002e+01 -0.88429 0.9205 -0.96069 -#> ds 5 m1 120 10.4 1.002e+01 -0.38429 0.9205 -0.41749
    # } +#> ds name time observed predicted residual std standardized +#> ds 1 parent 0 89.8 9.838e+01 8.584661 7.7094 1.113536 +#> ds 1 parent 0 104.1 9.838e+01 -5.715339 7.7094 -0.741350 +#> ds 1 parent 1 88.7 9.388e+01 5.182489 7.3611 0.704041 +#> ds 1 parent 1 95.5 9.388e+01 -1.617511 7.3611 -0.219739 +#> ds 1 parent 3 81.8 8.563e+01 3.825382 6.7229 0.569010 +#> ds 1 parent 3 94.5 8.563e+01 -8.874618 6.7229 -1.320062 +#> ds 1 parent 7 71.5 7.169e+01 0.188290 5.6482 0.033336 +#> ds 1 parent 7 70.3 7.169e+01 1.388290 5.6482 0.245795 +#> ds 1 parent 14 54.2 5.361e+01 -0.586595 4.2624 -0.137621 +#> ds 1 parent 14 49.6 5.361e+01 4.013405 4.2624 0.941587 +#> ds 1 parent 28 31.5 3.219e+01 0.688936 2.6496 0.260011 +#> ds 1 parent 28 28.8 3.219e+01 3.388936 2.6496 1.279016 +#> ds 1 parent 60 12.1 1.278e+01 0.678998 1.3145 0.516562 +#> ds 1 parent 60 13.6 1.278e+01 -0.821002 1.3145 -0.624595 +#> ds 1 parent 90 6.2 6.157e+00 -0.043461 0.9835 -0.044188 +#> ds 1 parent 90 8.3 6.157e+00 -2.143461 0.9835 -2.179316 +#> ds 1 parent 120 2.2 3.076e+00 0.876218 0.8916 0.982775 +#> ds 1 parent 120 2.4 3.076e+00 0.676218 0.8916 0.758453 +#> ds 1 m1 1 0.3 1.134e+00 0.833749 0.8633 0.965750 +#> ds 1 m1 1 0.2 1.134e+00 0.933749 0.8633 1.081583 +#> ds 1 m1 3 2.2 3.157e+00 0.957400 0.8933 1.071763 +#> ds 1 m1 3 3.0 3.157e+00 0.157400 0.8933 0.176202 +#> ds 1 m1 7 6.5 6.369e+00 -0.130995 0.9917 -0.132090 +#> ds 1 m1 7 5.0 6.369e+00 1.369005 0.9917 1.380438 +#> ds 1 m1 14 10.2 9.971e+00 -0.229362 1.1577 -0.198112 +#> ds 1 m1 14 9.5 9.971e+00 0.470638 1.1577 0.406513 +#> ds 1 m1 28 12.2 1.265e+01 0.447735 1.3067 0.342637 +#> ds 1 m1 28 13.4 1.265e+01 -0.752265 1.3067 -0.575683 +#> ds 1 m1 60 11.8 1.097e+01 -0.832027 1.2112 -0.686945 +#> ds 1 m1 60 13.2 1.097e+01 -2.232027 1.2112 -1.842825 +#> ds 1 m1 90 6.6 7.876e+00 1.275985 1.0553 1.209109 +#> ds 1 m1 90 9.3 7.876e+00 -1.424015 1.0553 -1.349381 +#> ds 1 m1 120 3.5 5.336e+00 1.835829 0.9540 1.924292 +#> ds 1 m1 120 5.4 5.336e+00 -0.064171 0.9540 -0.067263 +#> ds 2 parent 0 118.0 1.092e+02 -8.812058 8.5459 -1.031142 +#> ds 2 parent 0 99.8 1.092e+02 9.387942 8.5459 1.098529 +#> ds 2 parent 1 90.2 1.023e+02 12.114268 8.0135 1.511724 +#> ds 2 parent 1 94.6 1.023e+02 7.714268 8.0135 0.962654 +#> ds 2 parent 3 96.1 9.066e+01 -5.436165 7.1122 -0.764344 +#> ds 2 parent 3 78.4 9.066e+01 12.263835 7.1122 1.724339 +#> ds 2 parent 7 77.9 7.365e+01 -4.245773 5.7995 -0.732090 +#> ds 2 parent 7 77.7 7.365e+01 -4.045773 5.7995 -0.697604 +#> ds 2 parent 14 56.0 5.593e+01 -0.073803 4.4389 -0.016626 +#> ds 2 parent 14 54.7 5.593e+01 1.226197 4.4389 0.276236 +#> ds 2 parent 28 36.6 3.892e+01 2.320837 3.1502 0.736737 +#> ds 2 parent 28 36.8 3.892e+01 2.120837 3.1502 0.673248 +#> ds 2 parent 60 22.1 2.136e+01 -0.741020 1.8719 -0.395868 +#> ds 2 parent 60 24.7 2.136e+01 -3.341020 1.8719 -1.784841 +#> ds 2 parent 90 12.4 1.251e+01 0.113999 1.2989 0.087765 +#> ds 2 parent 90 10.8 1.251e+01 1.713999 1.2989 1.319575 +#> ds 2 parent 120 6.8 7.338e+00 0.537708 1.0315 0.521281 +#> ds 2 parent 120 7.9 7.338e+00 -0.562292 1.0315 -0.545113 +#> ds 2 m1 1 1.3 1.576e+00 0.276176 0.8675 0.318352 +#> ds 2 m1 3 3.7 4.177e+00 0.476741 0.9183 0.519146 +#> ds 2 m1 3 4.7 4.177e+00 -0.523259 0.9183 -0.569801 +#> ds 2 m1 7 8.1 7.724e+00 -0.376365 1.0485 -0.358970 +#> ds 2 m1 7 7.9 7.724e+00 -0.176365 1.0485 -0.168214 +#> ds 2 m1 14 10.1 1.077e+01 0.674433 1.2006 0.561738 +#> ds 2 m1 14 10.3 1.077e+01 0.474433 1.2006 0.395158 +#> ds 2 m1 28 10.7 1.212e+01 1.416179 1.2758 1.110010 +#> ds 2 m1 28 12.2 1.212e+01 -0.083821 1.2758 -0.065699 +#> ds 2 m1 60 10.7 1.041e+01 -0.294930 1.1807 -0.249793 +#> ds 2 m1 60 12.5 1.041e+01 -2.094930 1.1807 -1.774316 +#> ds 2 m1 90 9.1 8.079e+00 -1.020859 1.0646 -0.958929 +#> ds 2 m1 90 7.4 8.079e+00 0.679141 1.0646 0.637941 +#> ds 2 m1 120 6.1 5.968e+00 -0.131673 0.9765 -0.134843 +#> ds 2 m1 120 4.5 5.968e+00 1.468327 0.9765 1.503683 +#> ds 3 parent 0 106.2 1.036e+02 -2.638248 8.1101 -0.325303 +#> ds 3 parent 0 106.9 1.036e+02 -3.338248 8.1101 -0.411614 +#> ds 3 parent 1 107.4 9.580e+01 -11.600063 7.5094 -1.544743 +#> ds 3 parent 1 96.1 9.580e+01 -0.300063 7.5094 -0.039958 +#> ds 3 parent 3 79.4 8.297e+01 3.574516 6.5182 0.548391 +#> ds 3 parent 3 82.6 8.297e+01 0.374516 6.5182 0.057457 +#> ds 3 parent 7 63.9 6.517e+01 1.272397 5.1472 0.247200 +#> ds 3 parent 7 62.4 6.517e+01 2.772397 5.1472 0.538618 +#> ds 3 parent 14 51.0 4.821e+01 -2.790075 3.8512 -0.724475 +#> ds 3 parent 14 47.1 4.821e+01 1.109925 3.8512 0.288205 +#> ds 3 parent 28 36.1 3.385e+01 -2.250573 2.7723 -0.811811 +#> ds 3 parent 28 36.6 3.385e+01 -2.750573 2.7723 -0.992168 +#> ds 3 parent 60 20.1 1.964e+01 -0.455700 1.7543 -0.259760 +#> ds 3 parent 60 19.8 1.964e+01 -0.155700 1.7543 -0.088753 +#> ds 3 parent 90 11.3 1.210e+01 0.795458 1.2746 0.624068 +#> ds 3 parent 90 10.7 1.210e+01 1.395458 1.2746 1.094792 +#> ds 3 parent 120 8.2 7.451e+00 -0.749141 1.0364 -0.722816 +#> ds 3 parent 120 7.3 7.451e+00 0.150859 1.0364 0.145558 +#> ds 3 m1 0 0.8 3.695e-13 -0.800000 0.8588 -0.931542 +#> ds 3 m1 1 1.8 1.740e+00 -0.059741 0.8694 -0.068714 +#> ds 3 m1 1 2.3 1.740e+00 -0.559741 0.8694 -0.643812 +#> ds 3 m1 3 4.2 4.531e+00 0.331379 0.9285 0.356913 +#> ds 3 m1 3 4.1 4.531e+00 0.431379 0.9285 0.464618 +#> ds 3 m1 7 6.8 8.113e+00 1.312762 1.0661 1.231333 +#> ds 3 m1 7 10.1 8.113e+00 -1.987238 1.0661 -1.863971 +#> ds 3 m1 14 11.4 1.079e+01 -0.613266 1.2013 -0.510507 +#> ds 3 m1 14 12.8 1.079e+01 -2.013266 1.2013 -1.675923 +#> ds 3 m1 28 11.5 1.133e+01 -0.174252 1.2310 -0.141553 +#> ds 3 m1 28 10.6 1.133e+01 0.725748 1.2310 0.589558 +#> ds 3 m1 60 7.5 8.948e+00 1.448281 1.1059 1.309561 +#> ds 3 m1 60 8.6 8.948e+00 0.348281 1.1059 0.314922 +#> ds 3 m1 90 7.3 6.665e+00 -0.634932 1.0034 -0.632752 +#> ds 3 m1 90 8.1 6.665e+00 -1.434932 1.0034 -1.430004 +#> ds 3 m1 120 5.3 4.795e+00 -0.504936 0.9365 -0.539199 +#> ds 3 m1 120 3.8 4.795e+00 0.995064 0.9365 1.062586 +#> ds 4 parent 0 104.7 9.985e+01 -4.850494 7.8227 -0.620050 +#> ds 4 parent 0 88.3 9.985e+01 11.549506 7.8227 1.476402 +#> ds 4 parent 1 94.2 9.676e+01 2.556304 7.5834 0.337093 +#> ds 4 parent 1 94.6 9.676e+01 2.156304 7.5834 0.284346 +#> ds 4 parent 3 78.1 9.092e+01 12.817485 7.1318 1.797230 +#> ds 4 parent 3 96.5 9.092e+01 -5.582515 7.1318 -0.782764 +#> ds 4 parent 7 76.2 8.050e+01 4.297338 6.3270 0.679204 +#> ds 4 parent 7 77.8 8.050e+01 2.697338 6.3270 0.426320 +#> ds 4 parent 14 70.8 6.562e+01 -5.179989 5.1816 -0.999687 +#> ds 4 parent 14 67.3 6.562e+01 -1.679989 5.1816 -0.324222 +#> ds 4 parent 28 43.1 4.499e+01 1.886936 3.6069 0.523140 +#> ds 4 parent 28 45.1 4.499e+01 -0.113064 3.6069 -0.031346 +#> ds 4 parent 60 21.3 2.151e+01 0.214840 1.8827 0.114114 +#> ds 4 parent 60 23.5 2.151e+01 -1.985160 1.8827 -1.054433 +#> ds 4 parent 90 11.8 1.190e+01 0.098528 1.2633 0.077990 +#> ds 4 parent 90 12.1 1.190e+01 -0.201472 1.2633 -0.159475 +#> ds 4 parent 120 7.0 6.886e+00 -0.113832 1.0125 -0.112431 +#> ds 4 parent 120 6.2 6.886e+00 0.686168 1.0125 0.677724 +#> ds 4 m1 0 1.6 4.263e-14 -1.600000 0.8588 -1.863085 +#> ds 4 m1 1 0.9 7.140e-01 -0.185984 0.8606 -0.216112 +#> ds 4 m1 3 3.7 2.022e+00 -1.678243 0.8731 -1.922160 +#> ds 4 m1 3 2.0 2.022e+00 0.021757 0.8731 0.024919 +#> ds 4 m1 7 3.6 4.207e+00 0.607229 0.9192 0.660633 +#> ds 4 m1 7 3.8 4.207e+00 0.407229 0.9192 0.443044 +#> ds 4 m1 14 7.1 6.912e+00 -0.188339 1.0135 -0.185828 +#> ds 4 m1 14 6.6 6.912e+00 0.311661 1.0135 0.307506 +#> ds 4 m1 28 9.5 9.449e+00 -0.050714 1.1309 -0.044843 +#> ds 4 m1 28 9.3 9.449e+00 0.149286 1.1309 0.132004 +#> ds 4 m1 60 8.3 8.997e+00 0.697403 1.1083 0.629230 +#> ds 4 m1 60 9.0 8.997e+00 -0.002597 1.1083 -0.002343 +#> ds 4 m1 90 6.6 6.697e+00 0.096928 1.0047 0.096472 +#> ds 4 m1 90 7.7 6.697e+00 -1.003072 1.0047 -0.998348 +#> ds 4 m1 120 3.7 4.622e+00 0.921607 0.9312 0.989749 +#> ds 4 m1 120 3.5 4.622e+00 1.121607 0.9312 1.204537 +#> ds 5 parent 0 110.4 1.045e+02 -5.942426 8.1795 -0.726502 +#> ds 5 parent 0 112.1 1.045e+02 -7.642426 8.1795 -0.934338 +#> ds 5 parent 1 93.5 9.739e+01 3.893915 7.6327 0.510162 +#> ds 5 parent 1 91.0 9.739e+01 6.393915 7.6327 0.837700 +#> ds 5 parent 3 71.0 8.519e+01 14.188275 6.6891 2.121098 +#> ds 5 parent 3 89.7 8.519e+01 -4.511725 6.6891 -0.674487 +#> ds 5 parent 7 60.4 6.684e+01 6.439546 5.2753 1.220701 +#> ds 5 parent 7 59.1 6.684e+01 7.739546 5.2753 1.467133 +#> ds 5 parent 14 56.5 4.736e+01 -9.138979 3.7868 -2.413407 +#> ds 5 parent 14 47.0 4.736e+01 0.361021 3.7868 0.095338 +#> ds 5 parent 28 30.2 3.033e+01 0.131178 2.5132 0.052195 +#> ds 5 parent 28 23.9 3.033e+01 6.431178 2.5132 2.558936 +#> ds 5 parent 60 17.0 1.771e+01 0.705246 1.6243 0.434177 +#> ds 5 parent 60 18.7 1.771e+01 -0.994754 1.6243 -0.612409 +#> ds 5 parent 90 11.3 1.180e+01 0.504856 1.2580 0.401315 +#> ds 5 parent 90 11.9 1.180e+01 -0.095144 1.2580 -0.075631 +#> ds 5 parent 120 9.0 7.917e+00 -1.083499 1.0571 -1.024928 +#> ds 5 parent 120 8.1 7.917e+00 -0.183499 1.0571 -0.173579 +#> ds 5 m1 0 0.7 3.553e-15 -0.700000 0.8588 -0.815100 +#> ds 5 m1 1 3.0 3.204e+00 0.204414 0.8943 0.228572 +#> ds 5 m1 1 2.6 3.204e+00 0.604414 0.8943 0.675845 +#> ds 5 m1 3 5.1 8.586e+00 3.485889 1.0884 3.202858 +#> ds 5 m1 3 7.5 8.586e+00 1.085889 1.0884 0.997722 +#> ds 5 m1 7 16.5 1.612e+01 -0.376855 1.5211 -0.247743 +#> ds 5 m1 7 19.0 1.612e+01 -2.876855 1.5211 -1.891237 +#> ds 5 m1 14 22.9 2.267e+01 -0.228264 1.9633 -0.116267 +#> ds 5 m1 14 23.2 2.267e+01 -0.528264 1.9633 -0.269072 +#> ds 5 m1 28 22.2 2.468e+01 2.480178 2.1050 1.178211 +#> ds 5 m1 28 24.4 2.468e+01 0.280178 2.1050 0.133099 +#> ds 5 m1 60 15.5 1.860e+01 3.099615 1.6838 1.840794 +#> ds 5 m1 60 19.8 1.860e+01 -1.200385 1.6838 -0.712883 +#> ds 5 m1 90 14.9 1.326e+01 -1.636345 1.3433 -1.218195 +#> ds 5 m1 90 14.2 1.326e+01 -0.936345 1.3433 -0.697072 +#> ds 5 m1 120 10.9 9.348e+00 -1.551535 1.1258 -1.378133 +#> ds 5 m1 120 10.4 9.348e+00 -1.051535 1.1258 -0.934014
    # }
    diff --git a/man/nlme.Rd b/man/nlme.Rd index 307cca82..c367868b 100644 --- a/man/nlme.Rd +++ b/man/nlme.Rd @@ -8,7 +8,7 @@ \usage{ nlme_function(object) -mean_degparms(object, random = FALSE) +mean_degparms(object, random = FALSE, test_log_parms = FALSE, conf.level = 0.6) nlme_data(object) } @@ -16,6 +16,13 @@ nlme_data(object) \item{object}{An mmkin row object containing several fits of the same model to different datasets} \item{random}{Should a list with fixed and random effects be returned?} + +\item{test_log_parms}{If TRUE, log parameters are only considered in +the mean calculations if their untransformed counterparts (most likely +rate constants) pass the t-test for significant difference from zero.} + +\item{conf.level}{Possibility to adjust the required confidence level +for parameter that are tested if requested by 'test_log_parms'.} } \value{ A function that can be used with nlme @@ -60,7 +67,7 @@ grouped_data <- nlme_data(f) nlme_f <- nlme_function(f) # These assignments are necessary for these objects to be # visible to nlme and augPred when evaluation is done by -# pkgdown to generated the html docs. +# pkgdown to generate the html docs. assign("nlme_f", nlme_f, globalenv()) assign("grouped_data", grouped_data, globalenv()) diff --git a/man/saem.Rd b/man/saem.Rd index d5a8f17e..45f74e44 100644 --- a/man/saem.Rd +++ b/man/saem.Rd @@ -14,9 +14,12 @@ saem(object, ...) object, transformations = c("mkin", "saemix"), degparms_start = numeric(), + test_log_parms = FALSE, + conf.level = 0.6, solution_type = "auto", control = list(displayProgress = FALSE, print = FALSE, save = FALSE, save.graphs = FALSE), + fail_with_errors = TRUE, verbose = FALSE, quiet = FALSE, ... @@ -29,6 +32,7 @@ saemix_model( solution_type = "auto", transformations = c("mkin", "saemix"), degparms_start = numeric(), + test_log_parms = FALSE, verbose = FALSE, ... ) @@ -50,11 +54,22 @@ SFO or DFOP is used for the parent and there is either no metabolite or one.} \item{degparms_start}{Parameter values given as a named numeric vector will be used to override the starting values obtained from the 'mmkin' object.} +\item{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 \link{mean_degparms}.} + +\item{conf.level}{Possibility to adjust the required confidence level +for parameter that are tested if requested by 'test_log_parms'.} + \item{solution_type}{Possibility to specify the solution type in case the automatic choice is not desired} \item{control}{Passed to \link[saemix:saemix]{saemix::saemix}} +\item{fail_with_errors}{Should a failure to compute standard errors +from the inverse of the Fisher Information Matrix be a failure?} + \item{verbose}{Should we print information about created objects of type \link[saemix:SaemixModel-class]{saemix::SaemixModel} and \link[saemix:SaemixData-class]{saemix::SaemixData}?} @@ -104,7 +119,7 @@ f_saem_dfop <- saem(f_mmkin_parent["DFOP", ]) # The returned saem.mmkin object contains an SaemixObject, therefore we can use # functions from saemix library(saemix) -compare.saemix(list(f_saem_sfo$so, f_saem_fomc$so, f_saem_dfop$so)) +compare.saemix(f_saem_sfo$so, f_saem_fomc$so, f_saem_dfop$so) plot(f_saem_fomc$so, plot.type = "convergence") plot(f_saem_fomc$so, plot.type = "individual.fit") plot(f_saem_fomc$so, plot.type = "npde") @@ -112,7 +127,7 @@ plot(f_saem_fomc$so, plot.type = "vpc") f_mmkin_parent_tc <- update(f_mmkin_parent, error_model = "tc") f_saem_fomc_tc <- saem(f_mmkin_parent_tc["FOMC", ]) -compare.saemix(list(f_saem_fomc$so, f_saem_fomc_tc$so)) +compare.saemix(f_saem_fomc$so, f_saem_fomc_tc$so) sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"), A1 = mkinsub("SFO")) diff --git a/test.log b/test.log index 2c77a113..5f50c623 100644 --- a/test.log +++ b/test.log @@ -6,39 +6,39 @@ Testing mkin ✔ | 2 | Export dataset for reading into CAKE ✔ | 14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [1.0 s] ✔ | 4 | Calculation of FOCUS chi2 error levels [0.5 s] -✔ | 7 | Fitting the SFORB model [3.6 s] -✔ | 5 | Analytical solutions for coupled models [3.3 s] +✔ | 7 | Fitting the SFORB model [3.5 s] +✔ | 5 | Analytical solutions for coupled models [3.2 s] ✔ | 5 | Calculation of Akaike weights -✔ | 12 | Confidence intervals and p-values [1.0 s] -✔ | 14 | Error model fitting [4.6 s] +✔ | 12 | Confidence intervals and p-values [1.1 s] +✔ | 14 | Error model fitting [4.5 s] ✔ | 5 | Time step normalisation ✔ | 4 | Test fitting the decline of metabolites from their maximum [0.3 s] ✔ | 1 | Fitting the logistic model [0.2 s] -✔ | 34 1 | Nonlinear mixed-effects models [40.8 s] +✔ | 35 1 | Nonlinear mixed-effects models [27.1 s] ──────────────────────────────────────────────────────────────────────────────── -Skip (test_mixed.R:150:3): saem results are reproducible for biphasic fits +Skip (test_mixed.R:161:3): saem results are reproducible for biphasic fits Reason: Fitting with saemix takes around 10 minutes when using deSolve ──────────────────────────────────────────────────────────────────────────────── ✔ | 2 | Test dataset classes mkinds and mkindsg -✔ | 1 | mkinfit features [0.5 s] -✔ | 10 | Special cases of mkinfit calls [0.6 s] -✔ | 8 | mkinmod model generation and printing [0.4 s] -✔ | 3 | Model predictions with mkinpredict [0.7 s] -✔ | 16 | Evaluations according to 2015 NAFTA guidance [3.1 s] -✔ | 9 | Nonlinear mixed-effects models [14.5 s] -✔ | 16 | Plotting [2.1 s] +✔ | 1 | mkinfit features [0.3 s] +✔ | 10 | Special cases of mkinfit calls [0.3 s] +✔ | 8 | mkinmod model generation and printing [0.2 s] +✔ | 3 | Model predictions with mkinpredict [0.2 s] +✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.8 s] +✔ | 9 | Nonlinear mixed-effects models with nlme [8.1 s] +✔ | 16 | Plotting [2.0 s] ✔ | 4 | Residuals extracted from mkinfit models -✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.7 s] -✔ | 4 | Summary [0.2 s] +✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.5 s] +✔ | 4 | Summary [0.1 s] ✔ | 1 | Summaries of old mkinfit objects ✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.3 s] -✔ | 9 | Hypothesis tests [8.4 s] +✔ | 9 | Hypothesis tests [8.3 s] ✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.5 s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 92.7 s +Duration: 69.4 s ── Skipped tests ────────────────────────────────────────────────────────────── ● Fitting with saemix takes around 10 minutes when using deSolve (1) -[ FAIL 0 | WARN 0 | SKIP 1 | PASS 205 ] +[ FAIL 0 | WARN 0 | SKIP 1 | PASS 206 ] diff --git a/tests/figs/plotting/mixed-model-fit-for-nlme-object.svg b/tests/figs/plotting/mixed-model-fit-for-nlme-object.svg index 3eb2b0f8..db13b159 100644 --- a/tests/figs/plotting/mixed-model-fit-for-nlme-object.svg +++ b/tests/figs/plotting/mixed-model-fit-for-nlme-object.svg @@ -86,7 +86,7 @@ - + @@ -107,19 +107,19 @@ 80 100 120 - + - - - - - + + + + + 0 -20 -40 -60 -80 -100 +20 +40 +60 +80 +100 @@ -133,597 +133,597 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -734,30 +734,34 @@ - + - - - - - + + + + + 0 -20 -40 -60 -80 -100 - - - +20 +40 +60 +80 +100 + + + + - - --4 --2 + + + +-3 +-2 +-1 0 -2 -4 +1 +2 +3 @@ -772,582 +776,582 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/figs/plotting/mixed-model-fit-for-saem-object-with-mkin-transformations.svg b/tests/figs/plotting/mixed-model-fit-for-saem-object-with-mkin-transformations.svg index 0c2992d5..209b3dee 100644 --- a/tests/figs/plotting/mixed-model-fit-for-saem-object-with-mkin-transformations.svg +++ b/tests/figs/plotting/mixed-model-fit-for-saem-object-with-mkin-transformations.svg @@ -86,7 +86,7 @@ - + @@ -107,19 +107,19 @@ 80 100 120 - + - - - - - + + + + + 0 -20 -40 -60 -80 -100 +20 +40 +60 +80 +100 @@ -133,597 +133,597 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -734,30 +734,30 @@ - + - - - - - + + + + + 0 -20 -40 -60 -80 -100 - - - +20 +40 +60 +80 +100 + + + - - --4 --2 + + +-4 +-2 0 -2 -4 +2 +4 @@ -772,588 +772,588 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + @@ -1374,17 +1374,17 @@ 80 100 120 - + - - - - + + + + 0 -10 -20 -30 -40 +10 +20 +30 +40 @@ -1398,531 +1398,530 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -1933,28 +1932,28 @@ - + - - - - + + + + 0 -10 -20 -30 -40 - - - +10 +20 +30 +40 + + + - - --4 --2 + + +-4 +-2 0 -2 -4 +2 +4 @@ -1969,514 +1968,518 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/figs/plotting/mixed-model-fit-for-saem-object-with-saemix-transformations.svg b/tests/figs/plotting/mixed-model-fit-for-saem-object-with-saemix-transformations.svg index e65bc3bb..66d1d172 100644 --- a/tests/figs/plotting/mixed-model-fit-for-saem-object-with-saemix-transformations.svg +++ b/tests/figs/plotting/mixed-model-fit-for-saem-object-with-saemix-transformations.svg @@ -710,4 +710,9 @@ + + + + + diff --git a/tests/testthat/print_mmkin_biphasic_mixed.txt b/tests/testthat/print_mmkin_biphasic_mixed.txt index 11e11bfc..0b23fe58 100644 --- a/tests/testthat/print_mmkin_biphasic_mixed.txt +++ b/tests/testthat/print_mmkin_biphasic_mixed.txt @@ -8,7 +8,7 @@ d_m1/dt = + f_parent_to_m1 * ((k1 * g * exp(-k1 * time) + k2 * (1 - g) exp(-k2 * time))) * parent - k_m1 * m1 Data: -509 observations of 2 variable(s) grouped in 15 datasets +507 observations of 2 variable(s) grouped in 15 datasets object Status of individual fits: @@ -21,6 +21,6 @@ OK: No warnings Mean fitted parameters: parent_0 log_k_m1 f_parent_qlogis log_k1 log_k2 - 100.702 -5.347 -0.078 -2.681 -4.366 + 100.667 -5.378 -0.095 -2.743 -4.510 g_qlogis - -0.335 + -0.180 diff --git a/tests/testthat/print_nlme_biphasic.txt b/tests/testthat/print_nlme_biphasic.txt index f86bda76..f40d438d 100644 --- a/tests/testthat/print_nlme_biphasic.txt +++ b/tests/testthat/print_nlme_biphasic.txt @@ -9,21 +9,21 @@ d_m1/dt = + f_parent_to_m1 * ((k1 * g * exp(-k1 * time) + k2 * (1 - g) exp(-k2 * time))) * parent - k_m1 * m1 Data: -509 observations of 2 variable(s) grouped in 15 datasets +507 observations of 2 variable(s) grouped in 15 datasets -Log-likelihood: -1329 +Log-likelihood: -1326 Fixed effects: list(parent_0 ~ 1, log_k_m1 ~ 1, f_parent_qlogis ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_qlogis ~ 1) parent_0 log_k_m1 f_parent_qlogis log_k1 log_k2 - 100.43 -5.34 -0.08 -2.90 -4.34 + 100.7 -5.4 -0.1 -2.8 -4.5 g_qlogis - -0.19 + -0.1 Random effects: Formula: list(parent_0 ~ 1, log_k_m1 ~ 1, f_parent_qlogis ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_qlogis ~ 1) Level: ds Structure: Diagonal parent_0 log_k_m1 f_parent_qlogis log_k1 log_k2 g_qlogis Residual -StdDev: 1 0.1 0.3 0.6 0.5 0.3 3 +StdDev: 1 0.03 0.3 0.3 0.2 0.3 3 diff --git a/tests/testthat/print_sfo_saem_1.txt b/tests/testthat/print_sfo_saem_1.txt index d341e9e7..0c0e32ce 100644 --- a/tests/testthat/print_sfo_saem_1.txt +++ b/tests/testthat/print_sfo_saem_1.txt @@ -3,19 +3,19 @@ Structural model: d_parent/dt = - k_parent * parent Data: -264 observations of 1 variable(s) grouped in 15 datasets +262 observations of 1 variable(s) grouped in 15 datasets Likelihood computed by importance sampling AIC BIC logLik - 1320 1324 -654 + 1310 1315 -649 Fitted parameters: estimate lower upper -parent_0 1e+02 98.78 1e+02 +parent_0 1e+02 98.87 1e+02 k_parent 4e-02 0.03 4e-02 -Var.parent_0 8e-01 -1.76 3e+00 -Var.k_parent 9e-02 0.03 2e-01 -a.1 9e-01 0.70 1e+00 -b.1 4e-02 0.03 4e-02 -SD.parent_0 9e-01 -0.57 2e+00 +Var.parent_0 1e+00 -1.72 5e+00 +Var.k_parent 1e-01 0.03 2e-01 +a.1 9e-01 0.75 1e+00 +b.1 5e-02 0.04 5e-02 +SD.parent_0 1e+00 -0.12 3e+00 SD.k_parent 3e-01 0.20 4e-01 diff --git a/tests/testthat/setup_script.R b/tests/testthat/setup_script.R index 9229c198..96e865d4 100644 --- a/tests/testthat/setup_script.R +++ b/tests/testthat/setup_script.R @@ -106,6 +106,7 @@ const <- function(value) 2 set.seed(123456) SFO <- mkinmod(parent = mkinsub("SFO")) k_parent = rlnorm(n, log(0.03), log_sd) +set.seed(123456) ds_sfo <- lapply(1:n, function(i) { ds_mean <- mkinpredict(SFO, c(k_parent = k_parent[i]), c(parent = 100), sampling_times) @@ -118,6 +119,7 @@ fomc_pop <- list(parent_0 = 100, alpha = 2, beta = 8) fomc_parms <- as.matrix(data.frame( alpha = rlnorm(n, log(fomc_pop$alpha), 0.4), beta = rlnorm(n, log(fomc_pop$beta), 0.2))) +set.seed(123456) ds_fomc <- lapply(1:3, function(i) { ds_mean <- mkinpredict(FOMC, fomc_parms[i, ], c(parent = 100), sampling_times) @@ -131,6 +133,7 @@ dfop_parms <- as.matrix(data.frame( k1 = rlnorm(n, log(dfop_pop$k1), log_sd), k2 = rlnorm(n, log(dfop_pop$k2), log_sd), g = plogis(rnorm(n, qlogis(dfop_pop$g), log_sd)))) +set.seed(123456) ds_dfop <- lapply(1:n, function(i) { ds_mean <- mkinpredict(DFOP, dfop_parms[i, ], c(parent = dfop_pop$parent_0), sampling_times) @@ -144,6 +147,7 @@ hs_parms <- as.matrix(data.frame( k1 = rlnorm(n, log(hs_pop$k1), log_sd), k2 = rlnorm(n, log(hs_pop$k2), log_sd), tb = rlnorm(n, log(hs_pop$tb), 0.1))) +set.seed(123456) ds_hs <- lapply(1:10, function(i) { ds_mean <- mkinpredict(HS, hs_parms[i, ], c(parent = hs_pop$parent_0), sampling_times) @@ -171,6 +175,7 @@ ds_biphasic_mean <- lapply(1:n_biphasic, c(parent = 100, m1 = 0), sampling_times) } ) +set.seed(123456) ds_biphasic <- lapply(ds_biphasic_mean, function(ds) { add_err(ds, sdfunc = function(value) sqrt(err_1$const^2 + value^2 * err_1$prop^2), @@ -193,8 +198,18 @@ nlme_biphasic <- nlme(mmkin_biphasic) if (saemix_available) { sfo_saem_1 <- saem(mmkin_sfo_1, quiet = TRUE, transformations = "saemix") - dfop_saemix_1 <- saem(mmkin_dfop_1, quiet = TRUE, transformations = "mkin") - dfop_saemix_2 <- saem(mmkin_dfop_1, quiet = TRUE, transformations = "saemix") + # With default control parameters, we do not get good results with mkin + # transformations here + dfop_saemix_1 <- saem(mmkin_dfop_1, quiet = TRUE, transformations = "mkin", + control = list( + displayProgress = FALSE, print = FALSE, save = FALSE, save.graphs = FALSE, + rw.init = 1, nbiter.saemix = c(600, 100)) + ) + dfop_saemix_2 <- saem(mmkin_dfop_1, quiet = TRUE, transformations = "saemix", + control = list( + displayProgress = FALSE, print = FALSE, save = FALSE, save.graphs = FALSE, + rw.init = 0.5, nbiter.saemix = c(600, 100)) + ) saem_biphasic_m <- saem(mmkin_biphasic, transformations = "mkin", quiet = TRUE) saem_biphasic_s <- saem(mmkin_biphasic, transformations = "saemix", quiet = TRUE) diff --git a/tests/testthat/summary_nlme_biphasic_s.txt b/tests/testthat/summary_nlme_biphasic_s.txt index 65aead62..86049775 100644 --- a/tests/testthat/summary_nlme_biphasic_s.txt +++ b/tests/testthat/summary_nlme_biphasic_s.txt @@ -13,19 +13,19 @@ d_m1/dt = + f_parent_to_m1 * ((k1 * g * exp(-k1 * time) + k2 * (1 - g) exp(-k2 * time))) * parent - k_m1 * m1 Data: -509 observations of 2 variable(s) grouped in 15 datasets +507 observations of 2 variable(s) grouped in 15 datasets Model predictions using solution type analytical -Fitted in test time 0 s using 3 iterations +Fitted in test time 0 s using 4 iterations Variance model: Constant variance Mean of starting values for individual parameters: parent_0 log_k_m1 f_parent_qlogis log_k1 log_k2 - 100.70 -5.35 -0.08 -2.68 -4.37 + 100.67 -5.38 -0.09 -2.74 -4.51 g_qlogis - -0.33 + -0.18 Fixed degradation parameter values: value type @@ -34,40 +34,40 @@ m1_0 0 state Results: AIC BIC logLik - 2683 2738 -1329 + 2678 2733 -1326 Optimised, transformed parameters with symmetric confidence intervals: - lower est. upper -parent_0 99.6 100.43 101.26 -log_k_m1 -5.5 -5.34 -5.18 -f_parent_qlogis -0.3 -0.08 0.09 -log_k1 -3.2 -2.90 -2.60 -log_k2 -4.6 -4.34 -4.07 -g_qlogis -0.5 -0.19 0.08 + lower est. upper +parent_0 99.8 100.7 101.62 +log_k_m1 -5.6 -5.4 -5.25 +f_parent_qlogis -0.3 -0.1 0.06 +log_k1 -3.0 -2.8 -2.57 +log_k2 -4.7 -4.5 -4.35 +g_qlogis -0.4 -0.1 0.17 Correlation: prnt_0 lg_k_1 f_prn_ log_k1 log_k2 -log_k_m1 -0.177 -f_parent_qlogis -0.164 0.385 -log_k1 0.108 -0.017 -0.025 -log_k2 0.036 0.054 0.008 0.096 -g_qlogis -0.068 -0.110 -0.030 -0.269 -0.267 +log_k_m1 -0.167 +f_parent_qlogis -0.145 0.380 +log_k1 0.170 0.005 -0.026 +log_k2 0.083 0.168 0.032 0.365 +g_qlogis -0.088 -0.170 -0.044 -0.472 -0.631 Random effects: Formula: list(parent_0 ~ 1, log_k_m1 ~ 1, f_parent_qlogis ~ 1, log_k1 ~ 1, log_k2 ~ 1, g_qlogis ~ 1) Level: ds Structure: Diagonal parent_0 log_k_m1 f_parent_qlogis log_k1 log_k2 g_qlogis Residual -StdDev: 1 0.1 0.3 0.6 0.5 0.3 3 +StdDev: 1 0.03 0.3 0.3 0.2 0.3 3 Backtransformed parameters with asymmetric confidence intervals: lower est. upper parent_0 1e+02 1e+02 1e+02 -k_m1 4e-03 5e-03 6e-03 +k_m1 4e-03 4e-03 5e-03 f_parent_to_m1 4e-01 5e-01 5e-01 -k1 4e-02 6e-02 7e-02 -k2 1e-02 1e-02 2e-02 +k1 5e-02 6e-02 8e-02 +k2 9e-03 1e-02 1e-02 g 4e-01 5e-01 5e-01 Resulting formation fractions: @@ -77,5 +77,5 @@ parent_sink 0.5 Estimated disappearance times: DT50 DT90 DT50back DT50_k1 DT50_k2 -parent 26 131 39 13 53 -m1 144 479 NA NA NA +parent 25 150 45 11 63 +m1 154 512 NA NA NA diff --git a/tests/testthat/summary_saem_biphasic_s.txt b/tests/testthat/summary_saem_biphasic_s.txt index 1e0f1ccc..8dfae367 100644 --- a/tests/testthat/summary_saem_biphasic_s.txt +++ b/tests/testthat/summary_saem_biphasic_s.txt @@ -13,7 +13,7 @@ d_m1/dt = + f_parent_to_m1 * ((k1 * g * exp(-k1 * time) + k2 * (1 - g) exp(-k2 * time))) * parent - k_m1 * m1 Data: -509 observations of 2 variable(s) grouped in 15 datasets +507 observations of 2 variable(s) grouped in 15 datasets Model predictions using solution type analytical @@ -23,9 +23,9 @@ Variance model: Constant variance Mean of starting values for individual parameters: parent_0 k_m1 f_parent_to_m1 k1 k2 - 1.0e+02 4.8e-03 4.8e-01 6.8e-02 1.3e-02 + 1.0e+02 4.6e-03 4.8e-01 6.4e-02 1.1e-02 g - 4.2e-01 + 4.6e-01 Fixed degradation parameter values: None @@ -34,37 +34,37 @@ Results: Likelihood computed by importance sampling AIC BIC logLik - 2645 2654 -1310 + 2702 2711 -1338 Optimised parameters: est. lower upper -parent_0 1.0e+02 99.627 1.0e+02 -k_m1 4.8e-03 0.004 5.6e-03 -f_parent_to_m1 4.8e-01 0.437 5.2e-01 -k1 6.5e-02 0.051 8.0e-02 -k2 1.2e-02 0.010 1.4e-02 -g 4.3e-01 0.362 5.0e-01 +parent_0 1.0e+02 1.0e+02 1.0e+02 +k_m1 4.7e-03 3.9e-03 5.6e-03 +f_parent_to_m1 4.8e-01 4.3e-01 5.2e-01 +k1 4.8e-02 3.1e-02 6.5e-02 +k2 1.3e-02 8.7e-03 1.7e-02 +g 5.0e-01 4.1e-01 5.8e-01 Correlation: prnt_0 k_m1 f_p__1 k1 k2 -k_m1 -0.156 -f_parent_to_m1 -0.157 0.372 -k1 0.159 0.000 -0.029 -k2 0.074 0.145 0.032 0.332 -g -0.072 -0.142 -0.044 -0.422 -0.570 +k_m1 -0.152 +f_parent_to_m1 -0.143 0.366 +k1 0.097 -0.014 -0.021 +k2 0.022 0.083 0.023 0.101 +g -0.084 -0.144 -0.044 -0.303 -0.364 Random effects: est. lower upper -SD.parent_0 1.14 0.251 2.03 -SD.k_m1 0.14 -0.073 0.35 -SD.f_parent_to_m1 0.29 0.176 0.41 -SD.k1 0.36 0.211 0.52 -SD.k2 0.18 0.089 0.27 -SD.g 0.32 0.098 0.53 +SD.parent_0 1.22 0.316 2.12 +SD.k_m1 0.15 -0.079 0.38 +SD.f_parent_to_m1 0.32 0.191 0.44 +SD.k1 0.66 0.416 0.90 +SD.k2 0.59 0.368 0.80 +SD.g 0.16 -0.373 0.70 Variance model: est. lower upper -a.1 2.7 2.5 2.9 +a.1 2.9 2.7 3 Resulting formation fractions: ff @@ -73,5 +73,5 @@ parent_sink 0.52 Estimated disappearance times: DT50 DT90 DT50back DT50_k1 DT50_k2 -parent 25 145 44 11 58 -m1 145 481 NA NA NA +parent 26 127 38 14 54 +m1 146 485 NA NA NA diff --git a/tests/testthat/test_mixed.R b/tests/testthat/test_mixed.R index 0eb1f0d5..5d15530b 100644 --- a/tests/testthat/test_mixed.R +++ b/tests/testthat/test_mixed.R @@ -53,20 +53,26 @@ test_that("Parent fits using saemix are correctly implemented", { expect_true(all(s_dfop_s2$confint_back[, "lower"] < dfop_pop)) expect_true(all(s_dfop_s2$confint_back[, "upper"] > dfop_pop)) + dfop_mmkin_means_trans_tested <- mean_degparms(mmkin_dfop_1, test_log_parms = TRUE) dfop_mmkin_means_trans <- apply(parms(mmkin_dfop_1, transformed = TRUE), 1, mean) + + dfop_mmkin_means_tested <- backtransform_odeparms(dfop_mmkin_means_trans_tested, mmkin_dfop_1$mkinmod) dfop_mmkin_means <- backtransform_odeparms(dfop_mmkin_means_trans, mmkin_dfop_1$mkinmod) - # We get < 22% deviations by averaging the transformed parameters + # We get < 20% deviations for parent_0 and k1 by averaging the transformed parameters + # If we average only parameters passing the t-test, the deviation for k2 is also < 20% rel_diff_mmkin <- (dfop_mmkin_means - dfop_pop) / dfop_pop - expect_true(all(rel_diff_mmkin < 0.22)) + rel_diff_mmkin_tested <- (dfop_mmkin_means_tested - dfop_pop) / dfop_pop + expect_true(all(rel_diff_mmkin[c("parent_0", "k1")] < 0.20)) + expect_true(all(rel_diff_mmkin_tested[c("parent_0", "k1", "k2")] < 0.20)) - # We get < 50% deviations with transformations made in mkin + # We get < 30% deviations with transformations made in mkin rel_diff_1 <- (s_dfop_s1$confint_back[, "est."] - dfop_pop) / dfop_pop expect_true(all(rel_diff_1 < 0.5)) - # We get < 12% deviations with transformations made in saemix + # We get < 20% deviations with transformations made in saemix rel_diff_2 <- (s_dfop_s2$confint_back[, "est."] - dfop_pop) / dfop_pop - expect_true(all(rel_diff_2 < 0.12)) + expect_true(all(rel_diff_2 < 0.2)) mmkin_hs_1 <- mmkin("HS", ds_hs, quiet = TRUE, error_model = "const", cores = n_cores) hs_saem_1 <- saem(mmkin_hs_1, quiet = TRUE) @@ -107,9 +113,14 @@ test_that("nlme results are reproducible to some degree", { expect_known_output(print(test_summary, digits = 1), "summary_nlme_biphasic_s.txt") + # k1 just fails the first test (lower bound of the ci), so we need to excluded it + dfop_no_k1 <- c("parent_0", "k_m1", "f_parent_to_m1", "k2", "g") + dfop_sfo_pop_no_k1 <- as.numeric(dfop_sfo_pop[dfop_no_k1]) dfop_sfo_pop <- as.numeric(dfop_sfo_pop) + ci_dfop_sfo_n <- summary(nlme_biphasic)$confint_back - expect_true(all(ci_dfop_sfo_n[, "lower"] < dfop_sfo_pop)) + + expect_true(all(ci_dfop_sfo_n[dfop_no_k1, "lower"] < dfop_sfo_pop_no_k1)) expect_true(all(ci_dfop_sfo_n[, "upper"] > dfop_sfo_pop)) }) @@ -155,4 +166,3 @@ test_that("saem results are reproducible for biphasic fits", { expect_true(all(ci_dfop_sfo_s_d[no_k2, "lower"] < dfop_sfo_pop[no_k2])) expect_true(all(ci_dfop_sfo_s_d[no_k1, "upper"] > dfop_sfo_pop[no_k1])) }) - diff --git a/tests/testthat/test_nlme.R b/tests/testthat/test_nlme.R index 989914da..a3bc9413 100644 --- a/tests/testthat/test_nlme.R +++ b/tests/testthat/test_nlme.R @@ -1,4 +1,4 @@ -context("Nonlinear mixed-effects models") +context("Nonlinear mixed-effects models with nlme") library(nlme) -- cgit v1.2.3 From cb112e53163f9dc63d439dba50ca051877d67a79 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 16 Mar 2021 16:47:17 +0100 Subject: Convenience option to set nbiter.saemix --- R/saem.R | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) (limited to 'R/saem.R') diff --git a/R/saem.R b/R/saem.R index 460edede..184890f4 100644 --- a/R/saem.R +++ b/R/saem.R @@ -36,7 +36,9 @@ utils::globalVariables(c("predicted", "std")) #' 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? -#' @param control Passed to [saemix::saemix] +#' @param nbiter.saemix Convenience option to increase the number of +#' 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 @@ -115,7 +117,9 @@ saem.mmkin <- function(object, test_log_parms = FALSE, conf.level = 0.6, solution_type = "auto", + nbiter.saemix = c(300, 100), control = list(displayProgress = FALSE, print = FALSE, + nbiter.saemix = nbiter.saemix, save = FALSE, save.graphs = FALSE), fail_with_errors = TRUE, verbose = FALSE, quiet = FALSE, ...) -- cgit v1.2.3 From 34d1c5f23edfb60548bc5a9dd99c2f3af92acef1 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Sat, 20 Mar 2021 21:26:40 +0100 Subject: Fix mkin calculation of saemix residuals --- R/saem.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'R/saem.R') diff --git a/R/saem.R b/R/saem.R index 184890f4..6f28a47a 100644 --- a/R/saem.R +++ b/R/saem.R @@ -164,7 +164,7 @@ saem.mmkin <- function(object, xidep = return_data[c("time", "name")]) return_data <- transform(return_data, - residual = predicted - value, + residual = value - predicted, std = sigma_twocomp(predicted, f_saemix@results@respar[1], f_saemix@results@respar[2])) return_data <- transform(return_data, -- cgit v1.2.3 From c6eb6b2bb598002523c3d34d71b0e4a99671ccd6 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 9 Jun 2021 16:53:31 +0200 Subject: Rudimentary support for setting up nlmixr models - All degradation models are specified as ODE models. This appears to be fast enough - Error models are being translated to nlmixr as close to the mkin error model as possible. When using the 'saem' backend, it appears not to be possible to use the same error model for more than one observed variable - No support yet for models with parallel formation of metabolites, where the ilr transformation is used in mkin per default - There is a bug in nlmixr which appears to be triggered if the data are not balanced, see nlmixrdevelopment/nlmixr#530 - There is a print and a plot method, the summary method is not finished --- .travis.yml | 2 + DESCRIPTION | 9 +- NAMESPACE | 8 + NEWS.md | 14 +- R/mean_degparms.R | 61 ++++++ R/nlme.R | 55 ------ R/nlme.mmkin.R | 2 +- R/nlmixr.R | 467 ++++++++++++++++++++++++++++++++++++++++++++ R/plot.mixed.mmkin.R | 17 ++ R/saem.R | 6 + R/summary.nlmixr.mmkin.R | 250 ++++++++++++++++++++++++ build.log | 2 +- check.log | 68 +++++-- man/mean_degparms.Rd | 27 +++ man/nlme.Rd | 17 -- man/nlme.mmkin.Rd | 2 +- man/nlmixr.mmkin.Rd | 188 ++++++++++++++++++ man/plot.mixed.mmkin.Rd | 5 + man/summary.nlmixr.mmkin.Rd | 100 ++++++++++ man/summary.saem.mmkin.Rd | 24 +-- 20 files changed, 1209 insertions(+), 115 deletions(-) create mode 100644 R/mean_degparms.R create mode 100644 R/nlmixr.R create mode 100644 R/summary.nlmixr.mmkin.R create mode 100644 man/mean_degparms.Rd create mode 100644 man/nlmixr.mmkin.Rd create mode 100644 man/summary.nlmixr.mmkin.Rd (limited to 'R/saem.R') diff --git a/.travis.yml b/.travis.yml index 6c03b451..60e37230 100644 --- a/.travis.yml +++ b/.travis.yml @@ -11,6 +11,8 @@ addons: cache: packages repos: CRAN: https://cloud.r-project.org +r_github_packages: + - jranke/saemixextension@warp_combined script: - R CMD build . - R CMD check --no-tests mkin_*.tar.gz diff --git a/DESCRIPTION b/DESCRIPTION index 48aaf81f..5b90ef37 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: mkin Type: Package Title: Kinetic Evaluation of Chemical Degradation Data -Version: 1.0.4.9000 -Date: 2021-04-21 +Version: 1.0.5 +Date: 2021-06-03 Authors@R: c( person("Johannes", "Ranke", role = c("aut", "cre", "cph"), email = "jranke@uni-bremen.de", @@ -18,11 +18,10 @@ Description: Calculation routines based on the FOCUS Kinetics Report (2006, note that no warranty is implied for correctness of results or fitness for a particular purpose. Depends: R (>= 2.15.1), parallel -Imports: stats, graphics, methods, deSolve, R6, inline (>= 0.3.17), numDeriv, - lmtest, pkgbuild, nlme (>= 3.1-151), purrr, saemix (>= 3.1.9000) +Imports: stats, graphics, methods, deSolve, R6, inline (>= 0.3.19), numDeriv, + lmtest, pkgbuild, nlme (>= 3.1-151), purrr, saemix, nlmixr Suggests: knitr, rbenchmark, tikzDevice, testthat, rmarkdown, covr, vdiffr, benchmarkme, tibble, stats4 -Remotes: github::saemixdevelopment/saemixextension License: GPL LazyLoad: yes LazyData: yes diff --git a/NAMESPACE b/NAMESPACE index f2497283..bb4f5f92 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -16,6 +16,7 @@ S3method(mixed,mmkin) S3method(mkinpredict,mkinfit) S3method(mkinpredict,mkinmod) S3method(nlme,mmkin) +S3method(nlmixr,mmkin) S3method(nobs,mkinfit) S3method(parms,mkinfit) S3method(parms,mmkin) @@ -30,6 +31,7 @@ S3method(print,mkinmod) S3method(print,mmkin) S3method(print,nafta) S3method(print,nlme.mmkin) +S3method(print,nlmixr.mmkin) S3method(print,saem.mmkin) S3method(print,summary.mkinfit) S3method(print,summary.nlme.mmkin) @@ -38,6 +40,7 @@ S3method(residuals,mkinfit) S3method(saem,mmkin) S3method(summary,mkinfit) S3method(summary,nlme.mmkin) +S3method(summary,nlmixr.mmkin) S3method(summary,saem.mmkin) S3method(update,mkinfit) S3method(update,mmkin) @@ -86,6 +89,8 @@ export(nafta) export(nlme) export(nlme_data) export(nlme_function) +export(nlmixr_data) +export(nlmixr_model) export(parms) export(plot_err) export(plot_res) @@ -102,6 +107,8 @@ importFrom(R6,R6Class) importFrom(grDevices,dev.cur) importFrom(lmtest,lrtest) importFrom(methods,signature) +importFrom(nlmixr,nlmixr) +importFrom(nlmixr,tableControl) importFrom(parallel,detectCores) importFrom(parallel,mclapply) importFrom(parallel,parLapply) @@ -135,4 +142,5 @@ importFrom(stats,shapiro.test) importFrom(stats,update) importFrom(stats,vcov) importFrom(utils,getFromNamespace) +importFrom(utils,packageVersion) importFrom(utils,write.table) diff --git a/NEWS.md b/NEWS.md index d80e152c..03098106 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,18 +1,12 @@ -# mkin 1.0.4.9000 - -## General - -- Switch to a versioning scheme where the fourth version component indicates development versions +# mkin 1.0.5 ## Mixed-effects models -- Reintroduce the interface to the current development version of saemix, in particular: - -- 'saemix_model' and 'saemix_data': Helper functions to set up nonlinear mixed-effects models for mmkin row objects +- Introduce an interface to nlmixr, supporting estimation methods 'saem' and 'focei': S3 method 'nlmixr.mmkin' using the helper functions 'nlmixr_model' and 'nlmixr_data' to set up nlmixr models for mmkin row objects, with summary and plot methods. -- 'saem': generic function to fit saemix models using 'saemix_model' and 'saemix_data', with a generator 'saem.mmkin', summary and plot methods +- Reintroduce the interface to current development versions (not on CRAN) of saemix, in particular the generic function 'saem' with a generator 'saem.mmkin', currently using 'saemix_model' and 'saemix_data', summary and plot methods -- 'mean_degparms': New argument 'test_log_parms' that makes the function only consider log-transformed parameters where the untransformed parameters pass the t-test for a certain confidence level. This can be used to obtain more plausible starting parameters for 'saem' +- 'mean_degparms': New argument 'test_log_parms' that makes the function only consider log-transformed parameters where the untransformed parameters pass the t-test for a certain confidence level. This can be used to obtain more plausible starting parameters for the different mixed-effects model backends - 'plot.mixed.mmkin': Gains arguments 'test_log_parms' and 'conf.level' diff --git a/R/mean_degparms.R b/R/mean_degparms.R new file mode 100644 index 00000000..ec7f4342 --- /dev/null +++ b/R/mean_degparms.R @@ -0,0 +1,61 @@ +#' Calculate mean degradation parameters for an mmkin row object +#' +#' @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? +#' @param test_log_parms If TRUE, log parameters are only considered in +#' the mean calculations if their untransformed counterparts (most likely +#' rate constants) pass the t-test for significant difference from zero. +#' @param conf.level Possibility to adjust the required confidence level +#' for parameter that are tested if requested by 'test_log_parms'. +#' @export +mean_degparms <- function(object, random = FALSE, test_log_parms = FALSE, conf.level = 0.6) +{ + if (nrow(object) > 1) stop("Only row objects allowed") + parm_mat_trans <- sapply(object, parms, transformed = TRUE) + + if (test_log_parms) { + parm_mat_dim <- dim(parm_mat_trans) + parm_mat_dimnames <- dimnames(parm_mat_trans) + + log_parm_trans_names <- grep("^log_", rownames(parm_mat_trans), value = TRUE) + log_parm_names <- gsub("^log_", "", log_parm_trans_names) + + t_test_back_OK <- matrix( + sapply(object, function(o) { + suppressWarnings(summary(o)$bpar[log_parm_names, "Pr(>t)"] < (1 - conf.level)) + }), nrow = length(log_parm_names)) + rownames(t_test_back_OK) <- log_parm_trans_names + + parm_mat_trans_OK <- parm_mat_trans + for (trans_parm in log_parm_trans_names) { + parm_mat_trans_OK[trans_parm, ] <- ifelse(t_test_back_OK[trans_parm, ], + parm_mat_trans[trans_parm, ], NA) + } + } else { + parm_mat_trans_OK <- parm_mat_trans + } + + mean_degparm_names <- setdiff(rownames(parm_mat_trans), names(object[[1]]$errparms)) + degparm_mat_trans <- parm_mat_trans[mean_degparm_names, , drop = FALSE] + degparm_mat_trans_OK <- parm_mat_trans_OK[mean_degparm_names, , drop = FALSE] + + fixed <- apply(degparm_mat_trans_OK, 1, mean, na.rm = TRUE) + if (random) { + random <- t(apply(degparm_mat_trans[mean_degparm_names, , drop = FALSE], 2, function(column) column - fixed)) + # If we only have one parameter, apply returns a vector so we get a single row + if (nrow(degparm_mat_trans) == 1) random <- t(random) + rownames(random) <- levels(nlme_data(object)$ds) + + # For nlmixr we can specify starting values for standard deviations eta, and + # we ignore uncertain parameters if test_log_parms is FALSE + eta <- apply(degparm_mat_trans_OK, 1, sd, na.rm = TRUE) + + return(list(fixed = fixed, random = list(ds = random), eta = eta)) + } else { + return(fixed) + } +} + diff --git a/R/nlme.R b/R/nlme.R index d235a094..8762f137 100644 --- a/R/nlme.R +++ b/R/nlme.R @@ -124,61 +124,6 @@ nlme_function <- function(object) { return(model_function) } -#' @rdname nlme -#' @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? -#' @param test_log_parms If TRUE, log parameters are only considered in -#' the mean calculations if their untransformed counterparts (most likely -#' rate constants) pass the t-test for significant difference from zero. -#' @param conf.level Possibility to adjust the required confidence level -#' for parameter that are tested if requested by 'test_log_parms'. -#' @export -mean_degparms <- function(object, random = FALSE, test_log_parms = FALSE, conf.level = 0.6) -{ - if (nrow(object) > 1) stop("Only row objects allowed") - parm_mat_trans <- sapply(object, parms, transformed = TRUE) - - if (test_log_parms) { - parm_mat_dim <- dim(parm_mat_trans) - parm_mat_dimnames <- dimnames(parm_mat_trans) - - log_parm_trans_names <- grep("^log_", rownames(parm_mat_trans), value = TRUE) - log_parm_names <- gsub("^log_", "", log_parm_trans_names) - - t_test_back_OK <- matrix( - sapply(object, function(o) { - suppressWarnings(summary(o)$bpar[log_parm_names, "Pr(>t)"] < (1 - conf.level)) - }), nrow = length(log_parm_names)) - rownames(t_test_back_OK) <- log_parm_trans_names - - parm_mat_trans_OK <- parm_mat_trans - for (trans_parm in log_parm_trans_names) { - parm_mat_trans_OK[trans_parm, ] <- ifelse(t_test_back_OK[trans_parm, ], - parm_mat_trans[trans_parm, ], NA) - } - } else { - parm_mat_trans_OK <- parm_mat_trans - } - - mean_degparm_names <- setdiff(rownames(parm_mat_trans), names(object[[1]]$errparms)) - degparm_mat_trans <- parm_mat_trans[mean_degparm_names, , drop = FALSE] - degparm_mat_trans_OK <- parm_mat_trans_OK[mean_degparm_names, , drop = FALSE] - - fixed <- apply(degparm_mat_trans_OK, 1, mean, na.rm = TRUE) - if (random) { - random <- t(apply(degparm_mat_trans[mean_degparm_names, , drop = FALSE], 2, function(column) column - fixed)) - # If we only have one parameter, apply returns a vector so we get a single row - if (nrow(degparm_mat_trans) == 1) random <- t(random) - rownames(random) <- levels(nlme_data(object)$ds) - return(list(fixed = fixed, random = list(ds = random))) - } else { - return(fixed) - } -} - #' @rdname nlme #' @importFrom purrr map_dfr #' @return A \code{\link{groupedData}} object diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R index 306600c6..a1aa32e5 100644 --- a/R/nlme.mmkin.R +++ b/R/nlme.mmkin.R @@ -135,7 +135,7 @@ nlme.mmkin <- function(model, data = "auto", function(el) eval(parse(text = paste(el, 1, sep = "~")))), random = pdDiag(fixed), groups, - start = mean_degparms(model, random = TRUE), + start = mean_degparms(model, random = TRUE, test_log_parms = TRUE), correlation = NULL, weights = NULL, subset, method = c("ML", "REML"), na.action = na.fail, naPattern, diff --git a/R/nlmixr.R b/R/nlmixr.R new file mode 100644 index 00000000..223b23a1 --- /dev/null +++ b/R/nlmixr.R @@ -0,0 +1,467 @@ +utils::globalVariables(c("predicted", "std")) + +#' Fit nonlinear mixed models using nlmixr +#' +#' This function uses [nlmixr::nlmixr()] as a backend for fitting nonlinear mixed +#' effects models created from [mmkin] row objects using the Stochastic Approximation +#' Expectation Maximisation algorithm (SAEM). +#' +#' An mmkin row object is essentially a list of mkinfit objects that have been +#' obtained by fitting the same model to a list of datasets using [mkinfit]. +#' +#' @importFrom nlmixr nlmixr tableControl +#' @param object An [mmkin] row object containing several fits of the same +#' [mkinmod] model to different datasets +#' @param est Estimation method passed to [nlmixr::nlmixr] +#' @param degparms_start Parameter values given as a named numeric vector will +#' 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]. +#' @param conf.level Possibility to adjust the required confidence level +#' 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 +#' @param control Passed to [nlmixr::nlmixr]. +#' @param \dots Passed to [nlmixr_model] +#' @return An S3 object of class 'nlmixr.mmkin', containing the fitted +#' [nlmixr::nlmixr] object as a list component named 'nm'. The +#' object also inherits from 'mixed.mmkin'. +#' @seealso [summary.nlmixr.mmkin] [plot.mixed.mmkin] +#' @examples +#' ds <- lapply(experimental_data_for_UBA_2019[6:10], +#' function(x) subset(x$data[c("name", "time", "value")])) +#' names(ds) <- paste("Dataset", 6:10) +#' f_mmkin_parent <- mmkin(c("SFO", "FOMC", "DFOP", "HS"), ds, quiet = TRUE, cores = 1) +#' f_mmkin_parent_tc <- mmkin(c("SFO", "FOMC", "DFOP"), ds, error_model = "tc", +#' cores = 1, quiet = TRUE) +#' +#' f_nlmixr_sfo_saem <- nlmixr(f_mmkin_parent["SFO", ], est = "saem") +#' f_nlmixr_sfo_focei <- nlmixr(f_mmkin_parent["SFO", ], est = "focei") +#' +#' f_nlmixr_fomc_saem <- nlmixr(f_mmkin_parent["FOMC", ], est = "saem") +#' f_nlmixr_fomc_focei <- nlmixr(f_mmkin_parent["FOMC", ], est = "focei") +#' +#' f_nlmixr_dfop_saem <- nlmixr(f_mmkin_parent["DFOP", ], est = "saem") +#' f_nlmixr_dfop_focei <- nlmixr(f_mmkin_parent["DFOP", ], est = "focei") +#' +#' f_nlmixr_hs_saem <- nlmixr(f_mmkin_parent["HS", ], est = "saem") +#' f_nlmixr_hs_focei <- nlmixr(f_mmkin_parent["HS", ], est = "focei") +#' +#' f_nlmixr_fomc_saem_tc <- nlmixr(f_mmkin_parent_tc["FOMC", ], est = "saem") +#' f_nlmixr_fomc_focei_tc <- nlmixr(f_mmkin_parent_tc["FOMC", ], est = "focei") +#' +#' AIC( +#' f_nlmixr_sfo_saem$nm, f_nlmixr_sfo_focei$nm, +#' f_nlmixr_fomc_saem$nm, f_nlmixr_fomc_focei$nm, +#' f_nlmixr_dfop_saem$nm, f_nlmixr_dfop_focei$nm, +#' f_nlmixr_hs_saem$nm, f_nlmixr_hs_focei$nm, +#' f_nlmixr_fomc_saem_tc$nm, f_nlmixr_fomc_focei_tc$nm) +#' +#' AIC(nlme(f_mmkin_parent["FOMC", ])) +#' AIC(nlme(f_mmkin_parent["HS", ])) +#' +#' # nlme is comparable to nlmixr with focei, saem finds a better +#' # solution, the two-component error model does not improve it +#' plot(f_nlmixr_fomc_saem) +#' +#' \dontrun{ +#' sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"), +#' A1 = mkinsub("SFO")) +#' fomc_sfo <- mkinmod(parent = mkinsub("FOMC", "A1"), +#' A1 = mkinsub("SFO")) +#' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"), +#' A1 = mkinsub("SFO")) +#' +#' f_mmkin_const <- mmkin(list( +#' "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo), +#' ds, quiet = TRUE, error_model = "const") +#' f_mmkin_obs <- mmkin(list( +#' "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo), +#' ds, quiet = TRUE, error_model = "obs") +#' f_mmkin_tc <- mmkin(list( +#' "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo), +#' ds, quiet = TRUE, error_model = "tc") +#' +#' # A single constant variance is currently only possible with est = 'focei' in nlmixr +#' f_nlmixr_sfo_sfo_focei_const <- nlmixr(f_mmkin_const["SFO-SFO", ], est = "focei") +#' f_nlmixr_fomc_sfo_focei_const <- nlmixr(f_mmkin_const["FOMC-SFO", ], est = "focei") +#' f_nlmixr_dfop_sfo_focei_const <- nlmixr(f_mmkin_const["DFOP-SFO", ], est = "focei") +#' +#' # Variance by variable is supported by 'saem' and 'focei' +#' f_nlmixr_fomc_sfo_saem_obs <- nlmixr(f_mmkin_obs["FOMC-SFO", ], est = "saem") +#' f_nlmixr_fomc_sfo_focei_obs <- nlmixr(f_mmkin_obs["FOMC-SFO", ], est = "focei") +#' f_nlmixr_dfop_sfo_saem_obs <- nlmixr(f_mmkin_obs["DFOP-SFO", ], est = "saem") +#' f_nlmixr_dfop_sfo_focei_obs <- nlmixr(f_mmkin_obs["DFOP-SFO", ], est = "focei") +#' +#' # Identical two-component error for all variables is only possible with +#' # est = 'focei' in nlmixr +#' f_nlmixr_fomc_sfo_focei_tc <- nlmixr(f_mmkin_tc["FOMC-SFO", ], est = "focei") +#' f_nlmixr_dfop_sfo_focei_tc <- nlmixr(f_mmkin_tc["DFOP-SFO", ], est = "focei") +#' +#' # Two-component error by variable is possible with both estimation methods +#' # Variance by variable is supported by 'saem' and 'focei' +#' f_nlmixr_fomc_sfo_saem_obs_tc <- nlmixr(f_mmkin_tc["FOMC-SFO", ], est = "saem", +#' error_model = "obs_tc") +#' f_nlmixr_fomc_sfo_focei_obs_tc <- nlmixr(f_mmkin_tc["FOMC-SFO", ], est = "focei", +#' error_model = "obs_tc") +#' f_nlmixr_dfop_sfo_saem_obs_tc <- nlmixr(f_mmkin_tc["DFOP-SFO", ], est = "saem", +#' error_model = "obs_tc") +#' f_nlmixr_dfop_sfo_focei_obs_tc <- nlmixr(f_mmkin_tc["DFOP-SFO", ], est = "focei", +#' error_model = "obs_tc") +#' +#' AIC( +#' f_nlmixr_sfo_sfo_focei_const$nm, +#' f_nlmixr_fomc_sfo_focei_const$nm, +#' f_nlmixr_dfop_sfo_focei_const$nm, +#' f_nlmixr_fomc_sfo_saem_obs$nm, +#' f_nlmixr_fomc_sfo_focei_obs$nm, +#' f_nlmixr_dfop_sfo_saem_obs$nm, +#' f_nlmixr_dfop_sfo_focei_obs$nm, +#' f_nlmixr_fomc_sfo_focei_tc$nm, +#' f_nlmixr_dfop_sfo_focei_tc$nm, +#' f_nlmixr_fomc_sfo_saem_obs_tc$nm, +#' f_nlmixr_fomc_sfo_focei_obs_tc$nm, +#' f_nlmixr_dfop_sfo_saem_obs_tc$nm, +#' f_nlmixr_dfop_sfo_focei_obs_tc$nm +#' ) +#' # Currently, FOMC-SFO with two-component error by variable fitted by focei gives the +#' # lowest AIC +#' plot(f_nlmixr_fomc_sfo_focei_obs_tc) +#' summary(f_nlmixr_fomc_sfo_focei_obs_tc) +#' } +#' @export +nlmixr.mmkin <- function(object, data = NULL, + est = NULL, control = list(), + table = tableControl(), + error_model = object[[1]]$err_mod, + test_log_parms = TRUE, + conf.level = 0.6, + ..., + save = NULL, + envir = parent.frame() +) +{ + m_nlmixr <- nlmixr_model(object, est = est, + error_model = error_model, add_attributes = TRUE, + test_log_parms = test_log_parms, conf.level = conf.level) + d_nlmixr <- nlmixr_data(object) + + mean_dp_start <- attr(m_nlmixr, "mean_dp_start") + mean_ep_start <- attr(m_nlmixr, "mean_ep_start") + + attributes(m_nlmixr) <- NULL + + fit_time <- system.time({ + f_nlmixr <- nlmixr(m_nlmixr, d_nlmixr, est = est) + }) + + if (is.null(f_nlmixr$CMT)) { + nlmixr_df <- as.data.frame(f_nlmixr[c("ID", "TIME", "DV", "IPRED", "IRES", "IWRES")]) + nlmixr_df$CMT <- as.character(object[[1]]$data$variable[1]) + } else { + nlmixr_df <- as.data.frame(f_nlmixr[c("ID", "TIME", "DV", "CMT", "IPRED", "IRES", "IWRES")]) + } + + return_data <- nlmixr_df %>% + dplyr::transmute(ds = ID, name = CMT, time = TIME, value = DV, + predicted = IPRED, residual = IRES, + std = IRES/IWRES, standardized = IWRES) + + bparms_optim <- backtransform_odeparms(f_nlmixr$theta, + object[[1]]$mkinmod, + object[[1]]$transform_rates, + object[[1]]$transform_fractions) + + result <- list( + mkinmod = object[[1]]$mkinmod, + mmkin = object, + transform_rates = object[[1]]$transform_rates, + transform_fractions = object[[1]]$transform_fractions, + nm = f_nlmixr, + est = est, + time = fit_time, + mean_dp_start = mean_dp_start, + mean_ep_start = mean_ep_start, + bparms.optim = bparms_optim, + bparms.fixed = object[[1]]$bparms.fixed, + data = return_data, + err_mod = error_model, + date.fit = date(), + nlmixrversion = as.character(utils::packageVersion("nlmixr")), + mkinversion = as.character(utils::packageVersion("mkin")), + Rversion = paste(R.version$major, R.version$minor, sep=".") + ) + + class(result) <- c("nlmixr.mmkin", "mixed.mmkin") + return(result) +} + +#' @export +#' @rdname nlmixr.mmkin +#' @param x An nlmixr.mmkin object to print +#' @param digits Number of digits to use for printing +print.nlmixr.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) { + cat("Kinetic nonlinear mixed-effects model fit by", x$est, "using nlmixr") + cat("\nStructural model:\n") + diffs <- x$mmkin[[1]]$mkinmod$diffs + nice_diffs <- gsub("^(d.*) =", "\\1/dt =", diffs) + writeLines(strwrap(nice_diffs, exdent = 11)) + cat("\nData:\n") + cat(nrow(x$data), "observations of", + length(unique(x$data$name)), "variable(s) grouped in", + length(unique(x$data$ds)), "datasets\n") + + cat("\nLikelihood:\n") + print(data.frame( + AIC = AIC(x$nm), + BIC = BIC(x$nm), + logLik = logLik(x$nm), + row.names = " "), digits = digits) + + cat("\nFitted parameters:\n") + print(x$nm$parFixed, digits = digits) + + invisible(x) +} + +#' @rdname nlmixr.mmkin +#' @return An function defining a model suitable for fitting with [nlmixr::nlmixr]. +#' @export +nlmixr_model <- function(object, + est = c("saem", "focei"), + degparms_start = "auto", + test_log_parms = FALSE, conf.level = 0.6, + error_model = object[[1]]$err_mod, add_attributes = FALSE) +{ + if (nrow(object) > 1) stop("Only row objects allowed") + est = match.arg(est) + + mkin_model <- object[[1]]$mkinmod + obs_vars <- names(mkin_model$spec) + + if (error_model == object[[1]]$err_mod) { + + if (length(object[[1]]$mkinmod$spec) > 1 & est == "saem") { + if (error_model == "const") { + message( + "Constant variance for more than one variable is not supported for est = 'saem'\n", + "Changing the error model to 'obs' (variance by observed variable)") + error_model <- "obs" + } + if (error_model =="tc") { + message( + "With est = 'saem', a different error model is required for each observed variable", + "Changing the error model to 'obs_tc' (Two-component error for each observed variable)") + error_model <- "obs_tc" + } + } + } + + degparms_mmkin <- mean_degparms(object, + test_log_parms = test_log_parms, + conf.level = conf.level, random = TRUE) + + degparms_optim <- degparms_mmkin$fixed + + degparms_optim <- degparms_mmkin$fixed + + if (degparms_start[1] == "auto") { + degparms_start <- degparms_optim + } + degparms_fixed <- object[[1]]$bparms.fixed + + degparms_optim_back_names <- names(backtransform_odeparms(degparms_optim, + object[[1]]$mkinmod, + object[[1]]$transform_rates, + object[[1]]$transform_fractions)) + names(degparms_optim_back_names) <- names(degparms_optim) + + odeini_optim_parm_names <- grep('_0$', names(degparms_optim), value = TRUE) + odeini_fixed_parm_names <- grep('_0$', names(degparms_fixed), value = TRUE) + + odeparms_fixed_names <- setdiff(names(degparms_fixed), odeini_fixed_parm_names) + odeparms_fixed <- degparms_fixed[odeparms_fixed_names] + + odeini_fixed <- degparms_fixed[odeini_fixed_parm_names] + names(odeini_fixed) <- gsub('_0$', '', odeini_fixed_parm_names) + + # Definition of the model function + f <- function(){} + + ini_block <- "ini({" + + # Initial values for all degradation parameters + for (parm_name in names(degparms_optim)) { + # As initials for state variables are not transformed, + # we need to modify the name here as we want to + # use the original name in the model block + ini_block <- paste0( + ini_block, + parm_name, " = ", + as.character(degparms_start[parm_name]), + "\n", + "eta.", parm_name, " ~ ", + as.character(degparms_mmkin$eta[parm_name]), + "\n" + ) + } + + # Error model parameters + error_model_mkin <- object[[1]]$err_mod + + errparm_names_mkin <- names(object[[1]]$errparms) + errparms_mkin <- sapply(errparm_names_mkin, function(parm_name) { + mean(sapply(object, function(x) x$errparms[parm_name])) + }) + + sigma_tc_mkin <- errparms_ini <- errparms_mkin[1] + + mean(unlist(sapply(object, function(x) x$data$observed)), na.rm = TRUE) * + errparms_mkin[2] + + if (error_model == "const") { + if (error_model_mkin == "tc") { + errparms_ini <- sigma_tc_mkin + } else { + errparms_ini <- mean(errparms_mkin) + } + names(errparms_ini) <- "sigma" + } + + if (error_model == "obs") { + errparms_ini <- switch(error_model_mkin, + const = rep(errparms_mkin["sigma"], length(obs_vars)), + obs = errparms_mkin, + tc = sigma_tc_mkin) + names(errparms_ini) <- paste0("sigma_", obs_vars) + } + + if (error_model == "tc") { + if (error_model_mkin != "tc") { + stop("Not supported") + } else { + errparms_ini <- errparms_mkin + } + } + + if (error_model == "obs_tc") { + if (error_model_mkin != "tc") { + stop("Not supported") + } else { + errparms_ini <- rep(errparms_mkin, length(obs_vars)) + names(errparms_ini) <- paste0( + rep(names(errparms_mkin), length(obs_vars)), + "_", + rep(obs_vars, each = 2)) + } + } + + for (parm_name in names(errparms_ini)) { + ini_block <- paste0( + ini_block, + parm_name, " = ", + as.character(errparms_ini[parm_name]), + "\n" + ) + } + + ini_block <- paste0(ini_block, "})") + + body(f)[2] <- parse(text = ini_block) + + model_block <- "model({" + + # Population initial values for the ODE state variables + for (parm_name in odeini_optim_parm_names) { + model_block <- paste0( + model_block, + parm_name, "_model = ", + parm_name, " + eta.", parm_name, "\n", + gsub("(.*)_0", "\\1(0)", parm_name), " = ", parm_name, "_model\n") + } + + # Population initial values for log rate constants + for (parm_name in grep("^log_", names(degparms_optim), value = TRUE)) { + model_block <- paste0( + model_block, + gsub("^log_", "", parm_name), " = ", + "exp(", parm_name, " + eta.", parm_name, ")\n") + } + + # Population initial values for logit transformed parameters + for (parm_name in grep("_qlogis$", names(degparms_optim), value = TRUE)) { + model_block <- paste0( + model_block, + degparms_optim_back_names[parm_name], " = ", + "expit(", parm_name, " + eta.", parm_name, ")\n") + } + + # Differential equations + model_block <- paste0( + model_block, + paste( + gsub("d_(.*) =", "d/dt(\\1) =", mkin_model$diffs), + collapse = "\n"), + "\n" + ) + + # Error model + if (error_model == "const") { + model_block <- paste0(model_block, + paste(paste0(obs_vars, " ~ add(sigma)"), collapse = "\n")) + } + if (error_model == "obs") { + model_block <- paste0(model_block, + paste(paste0(obs_vars, " ~ add(sigma_", obs_vars, ")"), collapse = "\n"), + "\n") + } + if (error_model == "tc") { + model_block <- paste0(model_block, + paste(paste0(obs_vars, " ~ add(sigma_low) + prop(rsd_high)"), collapse = "\n"), + "\n") + } + if (error_model == "obs_tc") { + model_block <- paste0(model_block, + paste( + paste0(obs_vars, " ~ add(sigma_low_", obs_vars, ") + ", + "prop(rsd_high_", obs_vars, ")"), collapse = "\n"), + "\n") + } + + model_block <- paste0(model_block, "})") + + body(f)[3] <- parse(text = model_block) + + if (add_attributes) { + attr(f, "mean_dp_start") <- degparms_optim + attr(f, "mean_ep_start") <- errparms_ini + } + + return(f) +} + +#' @rdname nlmixr.mmkin +#' @return An dataframe suitable for use with [nlmixr::nlmixr] +#' @export +nlmixr_data <- function(object, ...) { + if (nrow(object) > 1) stop("Only row objects allowed") + d <- lapply(object, function(x) x$data) + compartment_map <- 1:length(object[[1]]$mkinmod$spec) + names(compartment_map) <- names(object[[1]]$mkinmod$spec) + ds_names <- colnames(object) + + ds_list <- lapply(object, function(x) x$data[c("time", "variable", "observed")]) + names(ds_list) <- ds_names + ds_nlmixr <- purrr::map_dfr(ds_list, function(x) x, .id = "ds") + ds_nlmixr$variable <- as.character(ds_nlmixr$variable) + ds_nlmixr_renamed <- data.frame( + ID = ds_nlmixr$ds, + TIME = ds_nlmixr$time, + AMT = 0, EVID = 0, + CMT = ds_nlmixr$variable, + DV = ds_nlmixr$observed, + stringsAsFactors = FALSE) + + return(ds_nlmixr_renamed) +} diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R index f0682c10..1ac62b07 100644 --- a/R/plot.mixed.mmkin.R +++ b/R/plot.mixed.mmkin.R @@ -40,12 +40,17 @@ utils::globalVariables("ds") #' #' # For this fit we need to increase pnlsMaxiter, and we increase the #' # tolerance in order to speed up the fit for this example evaluation +#' # It still takes 20 seconds to run #' f_nlme <- nlme(f, control = list(pnlsMaxIter = 120, tolerance = 1e-3)) #' plot(f_nlme) #' #' f_saem <- saem(f, transformations = "saemix") #' plot(f_saem) #' +#' f_obs <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, error_model = "obs") +#' f_nlmix <- nlmix(f_obs) +#' plot(f_nlmix) +#' #' # We can overlay the two variants if we generate predictions #' pred_nlme <- mkinpredict(dfop_sfo, #' f_nlme$bparms.optim[-1], @@ -109,6 +114,18 @@ plot.mixed.mmkin <- function(x, names(degparms_pop) <- degparms_i_names } + if (inherits(x, "nlmixr.mmkin")) { + eta_i <- random.effects(x$nm)[-1] + names(eta_i) <- gsub("^eta.", "", names(eta_i)) + degparms_i <- eta_i + degparms_pop <- x$nm$theta + for (parm_name in names(degparms_i)) { + degparms_i[parm_name] <- eta_i[parm_name] + degparms_pop[parm_name] + } + residual_type = ifelse(standardized, "standardized", "residual") + residuals <- x$data[[residual_type]] + } + degparms_fixed <- fit_1$fixed$value names(degparms_fixed) <- rownames(fit_1$fixed) degparms_all <- cbind(as.matrix(degparms_i), diff --git a/R/saem.R b/R/saem.R index 6f28a47a..5daf4be8 100644 --- a/R/saem.R +++ b/R/saem.R @@ -13,6 +13,7 @@ utils::globalVariables(c("predicted", "std")) #' psi0 of [saemix::saemixModel()] are the mean values of the parameters found #' using [mmkin]. #' +#' @importFrom utils packageVersion #' @param object An [mmkin] row object containing several fits of the same #' [mkinmod] model to different datasets #' @param verbose Should we print information about created objects of @@ -230,6 +231,11 @@ print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) { saemix_model <- function(object, solution_type = "auto", transformations = c("mkin", "saemix"), degparms_start = numeric(), test_log_parms = FALSE, verbose = FALSE, ...) { + if (packageVersion("saemix") < "3.1.9000") { + stop("To use the interface to saemix, you need to install a development version\n", + "preferably https://github.com/jranke/saemixextension@warp_combined") + } + if (nrow(object) > 1) stop("Only row objects allowed") mkin_model <- object[[1]]$mkinmod diff --git a/R/summary.nlmixr.mmkin.R b/R/summary.nlmixr.mmkin.R new file mode 100644 index 00000000..ae8e32cf --- /dev/null +++ b/R/summary.nlmixr.mmkin.R @@ -0,0 +1,250 @@ +#' Summary method for class "nlmixr.mmkin" +#' +#' Lists model equations, initial parameter values, optimised parameters +#' for fixed effects (population), random effects (deviations from the +#' population mean) and residual error model, as well as the resulting +#' endpoints such as formation fractions and DT50 values. Optionally +#' (default is FALSE), the data are listed in full. +#' +#' @param object an object of class [nlmix.mmkin] +#' @param x an object of class [summary.nlmix.mmkin] +#' @param data logical, indicating whether the full data should be included in +#' the summary. +#' @param verbose Should the summary be verbose? +#' @param distimes logical, indicating whether DT50 and DT90 values should be +#' included. +#' @param digits Number of digits to use for printing +#' @param \dots optional arguments passed to methods like \code{print}. +#' @return The summary function returns a list obtained in the fit, with at +#' least the following additional components +#' \item{nlmixrversion, mkinversion, Rversion}{The nlmixr, mkin and R versions used} +#' \item{date.fit, date.summary}{The dates where the fit and the summary were +#' produced} +#' \item{diffs}{The differential equations used in the degradation model} +#' \item{use_of_ff}{Was maximum or minimum use made of formation fractions} +#' \item{data}{The data} +#' \item{confint_trans}{Transformed parameters as used in the optimisation, with confidence intervals} +#' \item{confint_back}{Backtransformed parameters, with confidence intervals if available} +#' \item{confint_errmod}{Error model parameters with confidence intervals} +#' \item{ff}{The estimated formation fractions derived from the fitted +#' model.} +#' \item{distimes}{The DT50 and DT90 values for each observed variable.} +#' \item{SFORB}{If applicable, eigenvalues of SFORB components of the model.} +#' The print method is called for its side effect, i.e. printing the summary. +#' @importFrom stats predict vcov +#' @author Johannes Ranke for the mkin specific parts +#' nlmixr authors for the parts inherited from nlmixr. +#' @examples +#' # Generate five datasets following DFOP-SFO kinetics +#' sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120) +#' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "m1"), +#' m1 = mkinsub("SFO"), quiet = TRUE) +#' set.seed(1234) +#' k1_in <- rlnorm(5, log(0.1), 0.3) +#' k2_in <- rlnorm(5, log(0.02), 0.3) +#' g_in <- plogis(rnorm(5, qlogis(0.5), 0.3)) +#' f_parent_to_m1_in <- plogis(rnorm(5, qlogis(0.3), 0.3)) +#' k_m1_in <- rlnorm(5, log(0.02), 0.3) +#' +#' pred_dfop_sfo <- function(k1, k2, g, f_parent_to_m1, k_m1) { +#' mkinpredict(dfop_sfo, +#' c(k1 = k1, k2 = k2, g = g, f_parent_to_m1 = f_parent_to_m1, k_m1 = k_m1), +#' c(parent = 100, m1 = 0), +#' sampling_times) +#' } +#' +#' ds_mean_dfop_sfo <- lapply(1:5, function(i) { +#' mkinpredict(dfop_sfo, +#' c(k1 = k1_in[i], k2 = k2_in[i], g = g_in[i], +#' f_parent_to_m1 = f_parent_to_m1_in[i], k_m1 = k_m1_in[i]), +#' c(parent = 100, m1 = 0), +#' sampling_times) +#' }) +#' names(ds_mean_dfop_sfo) <- paste("ds", 1:5) +#' +#' ds_syn_dfop_sfo <- lapply(ds_mean_dfop_sfo, function(ds) { +#' add_err(ds, +#' sdfunc = function(value) sqrt(1^2 + value^2 * 0.07^2), +#' n = 1)[[1]] +#' }) +#' +#' \dontrun{ +#' # Evaluate using mmkin and nlmixr +#' f_mmkin_dfop_sfo <- mmkin(list(dfop_sfo), ds_syn_dfop_sfo, +#' quiet = TRUE, error_model = "tc", cores = 5) +#' f_saemix_dfop_sfo <- mkin::saem(f_mmkin_dfop_sfo) +#' f_nlme_dfop_sfo <- mkin::nlme(f_mmkin_dfop_sfo) +#' f_nlmixr_dfop_sfo_saem <- nlmixr(f_mmkin_dfop_sfo, est = "saem") +#' # The following takes a very long time but gives +#' f_nlmixr_dfop_sfo_focei <- nlmixr(f_mmkin_dfop_sfo, est = "focei") +#' AIC(f_nlmixr_dfop_sfo_saem$nm, f_nlmixr_dfop_sfo_focei$nm) +#' summary(f_nlmixr_dfop_sfo, data = TRUE) +#' } +#' +#' @export +summary.nlmixr.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes = TRUE, ...) { + + mod_vars <- names(object$mkinmod$diffs) + + pnames <- names(object$mean_dp_start) + np <- length(pnames) + + conf.int <- confint(object$nm) + confint_trans <- as.matrix(conf.int[pnames, c(1, 3, 4)]) + colnames(confint_trans) <- c("est.", "lower", "upper") + + bp <- backtransform_odeparms(confint_trans[, "est."], object$mkinmod, + object$transform_rates, object$transform_fractions) + bpnames <- names(bp) + + # Transform boundaries of CI for one parameter at a time, + # with the exception of sets of formation fractions (single fractions are OK). + f_names_skip <- character(0) + for (box in mod_vars) { # Figure out sets of fractions to skip + f_names <- grep(paste("^f", box, sep = "_"), pnames, value = TRUE) + n_paths <- length(f_names) + if (n_paths > 1) f_names_skip <- c(f_names_skip, f_names) + } + + confint_back <- matrix(NA, nrow = length(bp), ncol = 3, + dimnames = list(bpnames, colnames(confint_trans))) + confint_back[, "est."] <- bp + + for (pname in pnames) { + if (!pname %in% f_names_skip) { + par.lower <- confint_trans[pname, "lower"] + par.upper <- confint_trans[pname, "upper"] + names(par.lower) <- names(par.upper) <- pname + bpl <- backtransform_odeparms(par.lower, object$mkinmod, + object$transform_rates, + object$transform_fractions) + bpu <- backtransform_odeparms(par.upper, object$mkinmod, + object$transform_rates, + object$transform_fractions) + confint_back[names(bpl), "lower"] <- bpl + confint_back[names(bpu), "upper"] <- bpu + } + } + + # Correlation of fixed effects (inspired by summary.nlme) + varFix <- vcov(object$nm) + stdFix <- sqrt(diag(varFix)) + object$corFixed <- array( + t(varFix/stdFix)/stdFix, + dim(varFix), + list(pnames, pnames)) + + object$confint_back <- confint_back + + object$date.summary = date() + object$use_of_ff = object$mkinmod$use_of_ff + + object$diffs <- object$mkinmod$diffs + object$print_data <- data # boolean: Should we print the data? + predict(object$nm) + so_pred <- object$so@results@predictions + + names(object$data)[4] <- "observed" # rename value to observed + + object$verbose <- verbose + + object$fixed <- object$mmkin_orig[[1]]$fixed + object$AIC = AIC(object$so) + object$BIC = BIC(object$so) + object$logLik = logLik(object$so, method = "is") + + ep <- endpoints(object) + if (length(ep$ff) != 0) + object$ff <- ep$ff + if (distimes) object$distimes <- ep$distimes + if (length(ep$SFORB) != 0) object$SFORB <- ep$SFORB + class(object) <- c("summary.saem.mmkin") + return(object) +} + +#' @rdname summary.saem.mmkin +#' @export +print.summary.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), verbose = x$verbose, ...) { + cat("saemix version used for fitting: ", x$saemixversion, "\n") + cat("mkin version used for pre-fitting: ", x$mkinversion, "\n") + cat("R version used for fitting: ", x$Rversion, "\n") + + cat("Date of fit: ", x$date.fit, "\n") + cat("Date of summary:", x$date.summary, "\n") + + cat("\nEquations:\n") + nice_diffs <- gsub("^(d.*) =", "\\1/dt =", x[["diffs"]]) + writeLines(strwrap(nice_diffs, exdent = 11)) + + cat("\nData:\n") + cat(nrow(x$data), "observations of", + length(unique(x$data$name)), "variable(s) grouped in", + length(unique(x$data$ds)), "datasets\n") + + cat("\nModel predictions using solution type", x$solution_type, "\n") + + cat("\nFitted in", x$time[["elapsed"]], "s using", paste(x$so@options$nbiter.saemix, collapse = ", "), "iterations\n") + + cat("\nVariance model: ") + cat(switch(x$err_mod, + const = "Constant variance", + obs = "Variance unique to each observed variable", + tc = "Two-component variance function"), "\n") + + cat("\nMean of starting values for individual parameters:\n") + print(x$mean_dp_start, digits = digits) + + cat("\nFixed degradation parameter values:\n") + if(length(x$fixed$value) == 0) cat("None\n") + else print(x$fixed, digits = digits) + + cat("\nResults:\n\n") + cat("Likelihood computed by importance sampling\n") + print(data.frame(AIC = x$AIC, BIC = x$BIC, logLik = x$logLik, + row.names = " "), digits = digits) + + cat("\nOptimised parameters:\n") + print(x$confint_trans, digits = digits) + + if (nrow(x$confint_trans) > 1) { + corr <- x$corFixed + class(corr) <- "correlation" + print(corr, title = "\nCorrelation:", ...) + } + + cat("\nRandom effects:\n") + print(x$confint_ranef, digits = digits) + + cat("\nVariance model:\n") + print(x$confint_errmod, digits = digits) + + if (x$transformations == "mkin") { + cat("\nBacktransformed parameters:\n") + print(x$confint_back, digits = digits) + } + + printSFORB <- !is.null(x$SFORB) + if(printSFORB){ + cat("\nEstimated Eigenvalues of SFORB model(s):\n") + print(x$SFORB, digits = digits,...) + } + + printff <- !is.null(x$ff) + if(printff){ + cat("\nResulting formation fractions:\n") + print(data.frame(ff = x$ff), digits = digits,...) + } + + printdistimes <- !is.null(x$distimes) + if(printdistimes){ + cat("\nEstimated disappearance times:\n") + print(x$distimes, digits = digits,...) + } + + if (x$print_data){ + cat("\nData:\n") + print(format(x$data, digits = digits, ...), row.names = FALSE) + } + + invisible(x) +} diff --git a/build.log b/build.log index ca1c0481..13f76240 100644 --- a/build.log +++ b/build.log @@ -6,5 +6,5 @@ * creating vignettes ... OK * checking for LF line-endings in source and make files and shell scripts * checking for empty or unneeded directories -* building ‘mkin_1.0.4.9000.tar.gz’ +* building ‘mkin_1.0.5.tar.gz’ diff --git a/check.log b/check.log index 7de944a5..f6ee39db 100644 --- a/check.log +++ b/check.log @@ -1,19 +1,14 @@ * using log directory ‘/home/jranke/git/mkin/mkin.Rcheck’ -* using R version 4.0.5 (2021-03-31) +* using R version 4.1.0 (2021-05-18) * using platform: x86_64-pc-linux-gnu (64-bit) * using session charset: UTF-8 * using options ‘--no-tests --as-cran’ * checking for file ‘mkin/DESCRIPTION’ ... OK * checking extension type ... Package -* this is package ‘mkin’ version ‘1.0.4.9000’ +* this is package ‘mkin’ version ‘1.0.5’ * package encoding: UTF-8 -* checking CRAN incoming feasibility ... NOTE +* checking CRAN incoming feasibility ... Note_to_CRAN_maintainers Maintainer: ‘Johannes Ranke ’ - -Version contains large components (1.0.4.9000) - -Unknown, possibly mis-spelled, fields in DESCRIPTION: - ‘Remotes’ * checking package namespace information ... OK * checking package dependencies ... OK * checking if this is a source package ... OK @@ -46,22 +41,68 @@ Unknown, possibly mis-spelled, fields in DESCRIPTION: * checking S3 generic/method consistency ... OK * checking replacement functions ... OK * checking foreign function calls ... OK -* checking R code for possible problems ... OK +* checking R code for possible problems ... NOTE +saemix_model: no visible global function definition for + ‘packageVersion’ +Undefined global functions or variables: + packageVersion +Consider adding + importFrom("utils", "packageVersion") +to your NAMESPACE file. * checking Rd files ... OK * checking Rd metadata ... OK * checking Rd line widths ... OK -* checking Rd cross-references ... OK +* checking Rd cross-references ... WARNING +Missing link or links in documentation object 'nlmixr.mmkin.Rd': + ‘nlmix_model’ ‘summary.nlmixr.mmkin’ + +See section 'Cross-references' in the 'Writing R Extensions' manual. * checking for missing documentation entries ... OK * checking for code/documentation mismatches ... OK -* checking Rd \usage sections ... OK +* checking Rd \usage sections ... WARNING +Undocumented arguments in documentation object 'mean_degparms' + ‘object’ + +Undocumented arguments in documentation object 'nlmixr.mmkin' + ‘data’ ‘table’ ‘error_model’ ‘save’ ‘envir’ +Documented arguments not in \usage in documentation object 'nlmixr.mmkin': + ‘solution_type’ + +Functions with \usage entries need to have the appropriate \alias +entries, and all their arguments documented. +The \usage entries must correspond to syntactically valid R code. +See chapter ‘Writing R documentation files’ in the ‘Writing R +Extensions’ manual. * checking Rd contents ... OK * checking for unstated dependencies in examples ... OK * checking contents of ‘data’ directory ... OK * checking data for non-ASCII characters ... OK +* checking LazyData ... OK * checking data for ASCII and uncompressed saves ... OK * checking installed files from ‘inst/doc’ ... OK * checking files in ‘vignettes’ ... OK -* checking examples ... OK +* checking examples ... ERROR +Running examples in ‘mkin-Ex.R’ failed +The error most likely occurred in: + +> base::assign(".ptime", proc.time(), pos = "CheckExEnv") +> ### Name: nlmixr.mmkin +> ### Title: Fit nonlinear mixed models using nlmixr +> ### Aliases: nlmixr.mmkin print.nlmixr.mmkin nlmixr_model nlmixr_data +> +> ### ** Examples +> +> ds <- lapply(experimental_data_for_UBA_2019[6:10], ++ function(x) subset(x$data[c("name", "time", "value")])) +> names(ds) <- paste("Dataset", 6:10) +> f_mmkin_parent <- mmkin(c("SFO", "FOMC", "DFOP", "HS"), ds, quiet = TRUE, cores = 1) +> f_mmkin_parent_tc <- mmkin(c("SFO", "FOMC", "DFOP"), ds, error_model = "tc", ++ cores = 1, quiet = TRUE) +> +> f_nlmixr_sfo_saem <- nlmixr(f_mmkin_parent["SFO", ], est = "saem") +Error in nlmixr(f_mmkin_parent["SFO", ], est = "saem") : + could not find function "nlmixr" +Execution halted * checking for unstated dependencies in ‘tests’ ... OK * checking tests ... SKIPPED * checking for unstated dependencies in vignettes ... OK @@ -72,9 +113,8 @@ Unknown, possibly mis-spelled, fields in DESCRIPTION: * checking for detritus in the temp directory ... OK * DONE -Status: 1 NOTE +Status: 1 ERROR, 2 WARNINGs, 1 NOTE See ‘/home/jranke/git/mkin/mkin.Rcheck/00check.log’ for details. - diff --git a/man/mean_degparms.Rd b/man/mean_degparms.Rd new file mode 100644 index 00000000..92ed4c9d --- /dev/null +++ b/man/mean_degparms.Rd @@ -0,0 +1,27 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/mean_degparms.R +\name{mean_degparms} +\alias{mean_degparms} +\title{Calculate mean degradation parameters for an mmkin row object} +\usage{ +mean_degparms(object, random = FALSE, test_log_parms = FALSE, conf.level = 0.6) +} +\arguments{ +\item{random}{Should a list with fixed and random effects be returned?} + +\item{test_log_parms}{If TRUE, log parameters are only considered in +the mean calculations if their untransformed counterparts (most likely +rate constants) pass the t-test for significant difference from zero.} + +\item{conf.level}{Possibility to adjust the required confidence level +for parameter that are tested if requested by 'test_log_parms'.} +} +\value{ +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. +} +\description{ +Calculate mean degradation parameters for an mmkin row object +} diff --git a/man/nlme.Rd b/man/nlme.Rd index c367868b..e87b7a00 100644 --- a/man/nlme.Rd +++ b/man/nlme.Rd @@ -2,36 +2,19 @@ % Please edit documentation in R/nlme.R \name{nlme_function} \alias{nlme_function} -\alias{mean_degparms} \alias{nlme_data} \title{Helper functions to create nlme models from mmkin row objects} \usage{ nlme_function(object) -mean_degparms(object, random = FALSE, test_log_parms = FALSE, conf.level = 0.6) - nlme_data(object) } \arguments{ \item{object}{An mmkin row object containing several fits of the same model to different datasets} - -\item{random}{Should a list with fixed and random effects be returned?} - -\item{test_log_parms}{If TRUE, log parameters are only considered in -the mean calculations if their untransformed counterparts (most likely -rate constants) pass the t-test for significant difference from zero.} - -\item{conf.level}{Possibility to adjust the required confidence level -for parameter that are tested if requested by 'test_log_parms'.} } \value{ A function that can be used with nlme -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. - A \code{\link{groupedData}} object } \description{ diff --git a/man/nlme.mmkin.Rd b/man/nlme.mmkin.Rd index 2fb0488a..a2b45efa 100644 --- a/man/nlme.mmkin.Rd +++ b/man/nlme.mmkin.Rd @@ -13,7 +13,7 @@ paste(el, 1, sep = "~")))), random = pdDiag(fixed), groups, - start = mean_degparms(model, random = TRUE), + start = mean_degparms(model, random = TRUE, test_log_parms = TRUE), correlation = NULL, weights = NULL, subset, diff --git a/man/nlmixr.mmkin.Rd b/man/nlmixr.mmkin.Rd new file mode 100644 index 00000000..86bbdc9f --- /dev/null +++ b/man/nlmixr.mmkin.Rd @@ -0,0 +1,188 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/nlmixr.R +\name{nlmixr.mmkin} +\alias{nlmixr.mmkin} +\alias{print.nlmixr.mmkin} +\alias{nlmixr_model} +\alias{nlmixr_data} +\title{Fit nonlinear mixed models using nlmixr} +\usage{ +\method{nlmixr}{mmkin}( + object, + data = NULL, + est = NULL, + control = list(), + table = tableControl(), + error_model = object[[1]]$err_mod, + test_log_parms = TRUE, + conf.level = 0.6, + ..., + save = NULL, + envir = parent.frame() +) + +\method{print}{nlmixr.mmkin}(x, digits = max(3, getOption("digits") - 3), ...) + +nlmixr_model( + object, + est = c("saem", "focei"), + degparms_start = "auto", + test_log_parms = FALSE, + conf.level = 0.6, + error_model = object[[1]]$err_mod +) + +nlmixr_data(object, ...) +} +\arguments{ +\item{object}{An \link{mmkin} row object containing several fits of the same +\link{mkinmod} model to different datasets} + +\item{est}{Estimation method passed to \link[nlmixr:nlmixr]{nlmixr::nlmixr}} + +\item{control}{Passed to \link[nlmixr:nlmixr]{nlmixr::nlmixr}.} + +\item{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 \link{mean_degparms}.} + +\item{conf.level}{Possibility to adjust the required confidence level +for parameter that are tested if requested by 'test_log_parms'.} + +\item{\dots}{Passed to \link{nlmixr_model}} + +\item{x}{An nlmixr.mmkin object to print} + +\item{digits}{Number of digits to use for printing} + +\item{degparms_start}{Parameter values given as a named numeric vector will +be used to override the starting values obtained from the 'mmkin' object.} + +\item{solution_type}{Possibility to specify the solution type in case the +automatic choice is not desired} +} +\value{ +An S3 object of class 'nlmixr.mmkin', containing the fitted +\link[nlmixr:nlmixr]{nlmixr::nlmixr} object as a list component named 'nm'. The +object also inherits from 'mixed.mmkin'. + +An function defining a model suitable for fitting with \link[nlmixr:nlmixr]{nlmixr::nlmixr}. + +An dataframe suitable for use with \link[nlmixr:nlmixr]{nlmixr::nlmixr} +} +\description{ +This function uses \code{\link[nlmixr:nlmixr]{nlmixr::nlmixr()}} as a backend for fitting nonlinear mixed +effects models created from \link{mmkin} row objects using the Stochastic Approximation +Expectation Maximisation algorithm (SAEM). +} +\details{ +An mmkin row object is essentially a list of mkinfit objects that have been +obtained by fitting the same model to a list of datasets using \link{mkinfit}. +} +\examples{ +ds <- lapply(experimental_data_for_UBA_2019[6:10], + function(x) subset(x$data[c("name", "time", "value")])) +names(ds) <- paste("Dataset", 6:10) +f_mmkin_parent <- mmkin(c("SFO", "FOMC", "DFOP", "HS"), ds, quiet = TRUE, cores = 1) +f_mmkin_parent_tc <- mmkin(c("SFO", "FOMC", "DFOP"), ds, error_model = "tc", + cores = 1, quiet = TRUE) + +f_nlmixr_sfo_saem <- nlmixr(f_mmkin_parent["SFO", ], est = "saem") +f_nlmixr_sfo_focei <- nlmixr(f_mmkin_parent["SFO", ], est = "focei") + +f_nlmixr_fomc_saem <- nlmixr(f_mmkin_parent["FOMC", ], est = "saem") +f_nlmixr_fomc_focei <- nlmixr(f_mmkin_parent["FOMC", ], est = "focei") + +f_nlmixr_dfop_saem <- nlmixr(f_mmkin_parent["DFOP", ], est = "saem") +f_nlmixr_dfop_focei <- nlmixr(f_mmkin_parent["DFOP", ], est = "focei") + +f_nlmixr_hs_saem <- nlmixr(f_mmkin_parent["HS", ], est = "saem") +f_nlmixr_hs_focei <- nlmixr(f_mmkin_parent["HS", ], est = "focei") + +f_nlmixr_fomc_saem_tc <- nlmixr(f_mmkin_parent_tc["FOMC", ], est = "saem") +f_nlmixr_fomc_focei_tc <- nlmixr(f_mmkin_parent_tc["FOMC", ], est = "focei") + +AIC( + f_nlmixr_sfo_saem$nm, f_nlmixr_sfo_focei$nm, + f_nlmixr_fomc_saem$nm, f_nlmixr_fomc_focei$nm, + f_nlmixr_dfop_saem$nm, f_nlmixr_dfop_focei$nm, + f_nlmixr_hs_saem$nm, f_nlmixr_hs_focei$nm, + f_nlmixr_fomc_saem_tc$nm, f_nlmixr_fomc_focei_tc$nm) + +AIC(nlme(f_mmkin_parent["FOMC", ])) +AIC(nlme(f_mmkin_parent["HS", ])) + +# nlme is comparable to nlmixr with focei, saem finds a better +# solution, the two-component error model does not improve it +plot(f_nlmixr_fomc_saem) + +\dontrun{ +sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"), + A1 = mkinsub("SFO")) +fomc_sfo <- mkinmod(parent = mkinsub("FOMC", "A1"), + A1 = mkinsub("SFO")) +dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"), + A1 = mkinsub("SFO")) + +f_mmkin_const <- mmkin(list( + "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo), + ds, quiet = TRUE, error_model = "const") +f_mmkin_obs <- mmkin(list( + "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo), + ds, quiet = TRUE, error_model = "obs") +f_mmkin_tc <- mmkin(list( + "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo), + ds, quiet = TRUE, error_model = "tc") + +# A single constant variance is currently only possible with est = 'focei' in nlmixr +f_nlmixr_sfo_sfo_focei_const <- nlmixr(f_mmkin_const["SFO-SFO", ], est = "focei") +f_nlmixr_fomc_sfo_focei_const <- nlmixr(f_mmkin_const["FOMC-SFO", ], est = "focei") +f_nlmixr_dfop_sfo_focei_const <- nlmixr(f_mmkin_const["DFOP-SFO", ], est = "focei") + +# Variance by variable is supported by 'saem' and 'focei' +f_nlmixr_fomc_sfo_saem_obs <- nlmixr(f_mmkin_obs["FOMC-SFO", ], est = "saem") +f_nlmixr_fomc_sfo_focei_obs <- nlmixr(f_mmkin_obs["FOMC-SFO", ], est = "focei") +f_nlmixr_dfop_sfo_saem_obs <- nlmixr(f_mmkin_obs["DFOP-SFO", ], est = "saem") +f_nlmixr_dfop_sfo_focei_obs <- nlmixr(f_mmkin_obs["DFOP-SFO", ], est = "focei") + +# Identical two-component error for all variables is only possible with +# est = 'focei' in nlmixr +f_nlmixr_fomc_sfo_focei_tc <- nlmixr(f_mmkin_tc["FOMC-SFO", ], est = "focei") +f_nlmixr_dfop_sfo_focei_tc <- nlmixr(f_mmkin_tc["DFOP-SFO", ], est = "focei") + +# Two-component error by variable is possible with both estimation methods +# Variance by variable is supported by 'saem' and 'focei' +f_nlmixr_fomc_sfo_saem_obs_tc <- nlmixr(f_mmkin_tc["FOMC-SFO", ], est = "saem", + error_model = "obs_tc") +f_nlmixr_fomc_sfo_focei_obs_tc <- nlmixr(f_mmkin_tc["FOMC-SFO", ], est = "focei", + error_model = "obs_tc") +f_nlmixr_dfop_sfo_saem_obs_tc <- nlmixr(f_mmkin_tc["DFOP-SFO", ], est = "saem", + error_model = "obs_tc") +f_nlmixr_dfop_sfo_focei_obs_tc <- nlmixr(f_mmkin_tc["DFOP-SFO", ], est = "focei", + error_model = "obs_tc") + +AIC( + f_nlmixr_sfo_sfo_focei_const$nm, + f_nlmixr_fomc_sfo_focei_const$nm, + f_nlmixr_dfop_sfo_focei_const$nm, + f_nlmixr_fomc_sfo_saem_obs$nm, + f_nlmixr_fomc_sfo_focei_obs$nm, + f_nlmixr_dfop_sfo_saem_obs$nm, + f_nlmixr_dfop_sfo_focei_obs$nm, + f_nlmixr_fomc_sfo_focei_tc$nm, + f_nlmixr_dfop_sfo_focei_tc$nm, + f_nlmixr_fomc_sfo_saem_obs_tc$nm, + f_nlmixr_fomc_sfo_focei_obs_tc$nm, + f_nlmixr_dfop_sfo_saem_obs_tc$nm, + f_nlmixr_dfop_sfo_focei_obs_tc$nm +) +# Currently, FOMC-SFO with two-component error by variable fitted by focei gives the +# lowest AIC +plot(f_nlmixr_fomc_sfo_focei_obs_tc) +summary(f_nlmixr_fomc_sfo_focei_obs_tc) +} +} +\seealso{ +\link{summary.nlmixr.mmkin} \link{plot.mixed.mmkin} +} diff --git a/man/plot.mixed.mmkin.Rd b/man/plot.mixed.mmkin.Rd index bcab3e74..d87ca22c 100644 --- a/man/plot.mixed.mmkin.Rd +++ b/man/plot.mixed.mmkin.Rd @@ -99,12 +99,17 @@ plot(f[, 3:4], standardized = TRUE) # For this fit we need to increase pnlsMaxiter, and we increase the # tolerance in order to speed up the fit for this example evaluation +# It still takes 20 seconds to run f_nlme <- nlme(f, control = list(pnlsMaxIter = 120, tolerance = 1e-3)) plot(f_nlme) f_saem <- saem(f, transformations = "saemix") plot(f_saem) +f_obs <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, error_model = "obs") +f_nlmix <- nlmix(f_obs) +plot(f_nlmix) + # We can overlay the two variants if we generate predictions pred_nlme <- mkinpredict(dfop_sfo, f_nlme$bparms.optim[-1], diff --git a/man/summary.nlmixr.mmkin.Rd b/man/summary.nlmixr.mmkin.Rd new file mode 100644 index 00000000..03f0ffb2 --- /dev/null +++ b/man/summary.nlmixr.mmkin.Rd @@ -0,0 +1,100 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/summary.nlmixr.mmkin.R +\name{summary.nlmixr.mmkin} +\alias{summary.nlmixr.mmkin} +\title{Summary method for class "nlmixr.mmkin"} +\usage{ +\method{summary}{nlmixr.mmkin}(object, data = FALSE, verbose = FALSE, distimes = TRUE, ...) +} +\arguments{ +\item{object}{an object of class \link{nlmix.mmkin}} + +\item{data}{logical, indicating whether the full data should be included in +the summary.} + +\item{verbose}{Should the summary be verbose?} + +\item{distimes}{logical, indicating whether DT50 and DT90 values should be +included.} + +\item{\dots}{optional arguments passed to methods like \code{print}.} + +\item{x}{an object of class \link{summary.nlmix.mmkin}} + +\item{digits}{Number of digits to use for printing} +} +\value{ +The summary function returns a list obtained in the fit, with at +least the following additional components +\item{nlmixrversion, mkinversion, Rversion}{The nlmixr, mkin and R versions used} +\item{date.fit, date.summary}{The dates where the fit and the summary were +produced} +\item{diffs}{The differential equations used in the degradation model} +\item{use_of_ff}{Was maximum or minimum use made of formation fractions} +\item{data}{The data} +\item{confint_trans}{Transformed parameters as used in the optimisation, with confidence intervals} +\item{confint_back}{Backtransformed parameters, with confidence intervals if available} +\item{confint_errmod}{Error model parameters with confidence intervals} +\item{ff}{The estimated formation fractions derived from the fitted +model.} +\item{distimes}{The DT50 and DT90 values for each observed variable.} +\item{SFORB}{If applicable, eigenvalues of SFORB components of the model.} +The print method is called for its side effect, i.e. printing the summary. +} +\description{ +Lists model equations, initial parameter values, optimised parameters +for fixed effects (population), random effects (deviations from the +population mean) and residual error model, as well as the resulting +endpoints such as formation fractions and DT50 values. Optionally +(default is FALSE), the data are listed in full. +} +\examples{ +# Generate five datasets following DFOP-SFO kinetics +sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120) +dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "m1"), + m1 = mkinsub("SFO"), quiet = TRUE) +set.seed(1234) +k1_in <- rlnorm(5, log(0.1), 0.3) +k2_in <- rlnorm(5, log(0.02), 0.3) +g_in <- plogis(rnorm(5, qlogis(0.5), 0.3)) +f_parent_to_m1_in <- plogis(rnorm(5, qlogis(0.3), 0.3)) +k_m1_in <- rlnorm(5, log(0.02), 0.3) + +pred_dfop_sfo <- function(k1, k2, g, f_parent_to_m1, k_m1) { + mkinpredict(dfop_sfo, + c(k1 = k1, k2 = k2, g = g, f_parent_to_m1 = f_parent_to_m1, k_m1 = k_m1), + c(parent = 100, m1 = 0), + sampling_times) +} + +ds_mean_dfop_sfo <- lapply(1:5, function(i) { + mkinpredict(dfop_sfo, + c(k1 = k1_in[i], k2 = k2_in[i], g = g_in[i], + f_parent_to_m1 = f_parent_to_m1_in[i], k_m1 = k_m1_in[i]), + c(parent = 100, m1 = 0), + sampling_times) +}) +names(ds_mean_dfop_sfo) <- paste("ds", 1:5) + +ds_syn_dfop_sfo <- lapply(ds_mean_dfop_sfo, function(ds) { + add_err(ds, + sdfunc = function(value) sqrt(1^2 + value^2 * 0.07^2), + n = 1)[[1]] +}) + +\dontrun{ +# Evaluate using mmkin and nlmixr +f_mmkin_dfop_sfo <- mmkin(list(dfop_sfo), ds_syn_dfop_sfo, + quiet = TRUE, error_model = "obs", cores = 5) +f_saemix_dfop_sfo <- mkin::saem(f_mmkin_dfop_sfo) +f_nlme_dfop_sfo <- mkin::nlme(f_mmkin_dfop_sfo) +f_nlmixr_dfop_sfo_saem <- nlmixr(f_mmkin_dfop_sfo, est = "saem") +#f_nlmixr_dfop_sfo_focei <- nlmixr(f_mmkin_dfop_sfo, est = "focei") +summary(f_nlmixr_dfop_sfo, data = TRUE) +} + +} +\author{ +Johannes Ranke for the mkin specific parts +nlmixr authors for the parts inherited from nlmixr. +} diff --git a/man/summary.saem.mmkin.Rd b/man/summary.saem.mmkin.Rd index 67cb3cbb..86938d31 100644 --- a/man/summary.saem.mmkin.Rd +++ b/man/summary.saem.mmkin.Rd @@ -1,30 +1,32 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/summary.saem.mmkin.R -\name{summary.saem.mmkin} -\alias{summary.saem.mmkin} +% Please edit documentation in R/summary.nlmixr.mmkin.R, R/summary.saem.mmkin.R +\name{print.summary.saem.mmkin} \alias{print.summary.saem.mmkin} +\alias{summary.saem.mmkin} \title{Summary method for class "saem.mmkin"} \usage{ +\method{print}{summary.saem.mmkin}(x, digits = max(3, getOption("digits") - 3), verbose = x$verbose, ...) + \method{summary}{saem.mmkin}(object, data = FALSE, verbose = FALSE, distimes = TRUE, ...) \method{print}{summary.saem.mmkin}(x, digits = max(3, getOption("digits") - 3), verbose = x$verbose, ...) } \arguments{ -\item{object}{an object of class \link{saem.mmkin}} +\item{x}{an object of class \link{summary.saem.mmkin}} -\item{data}{logical, indicating whether the full data should be included in -the summary.} +\item{digits}{Number of digits to use for printing} \item{verbose}{Should the summary be verbose?} -\item{distimes}{logical, indicating whether DT50 and DT90 values should be -included.} - \item{\dots}{optional arguments passed to methods like \code{print}.} -\item{x}{an object of class \link{summary.saem.mmkin}} +\item{object}{an object of class \link{saem.mmkin}} -\item{digits}{Number of digits to use for printing} +\item{data}{logical, indicating whether the full data should be included in +the summary.} + +\item{distimes}{logical, indicating whether DT50 and DT90 values should be +included.} } \value{ The summary function returns a list based on the \link[saemix:SaemixObject-class]{saemix::SaemixObject} -- cgit v1.2.3 From 8f015900156981ecc2f1f6a1d5a078277ef3f454 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 23 Jun 2021 17:02:20 +0200 Subject: Test log parameters by default when deriving saemix starting parameters --- R/saem.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'R/saem.R') diff --git a/R/saem.R b/R/saem.R index 5daf4be8..9db2c04a 100644 --- a/R/saem.R +++ b/R/saem.R @@ -115,7 +115,7 @@ saem <- function(object, ...) UseMethod("saem") saem.mmkin <- function(object, transformations = c("mkin", "saemix"), degparms_start = numeric(), - test_log_parms = FALSE, + test_log_parms = TRUE, conf.level = 0.6, solution_type = "auto", nbiter.saemix = c(300, 100), -- cgit v1.2.3 From d75378911cef79b3ed95daef71bf67db413d2ac8 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 17 Nov 2021 12:59:49 +0100 Subject: Update required saemix version, update tests --- R/saem.R | 5 +- test.log | 60 +- ...t-for-saem-object-with-mkin-transformations.svg | 2306 ++++++++++---------- ...for-saem-object-with-saemix-transformations.svg | 572 ++--- tests/testthat/print_sfo_saem_1.txt | 8 +- tests/testthat/setup_script.R | 2 +- tests/testthat/summary_saem_biphasic_s.txt | 38 +- vignettes/FOCUS_D.html | 8 +- vignettes/FOCUS_L.html | 52 +- 9 files changed, 1525 insertions(+), 1526 deletions(-) (limited to 'R/saem.R') diff --git a/R/saem.R b/R/saem.R index 9db2c04a..2c20f788 100644 --- a/R/saem.R +++ b/R/saem.R @@ -231,9 +231,8 @@ print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) { saemix_model <- function(object, solution_type = "auto", transformations = c("mkin", "saemix"), degparms_start = numeric(), test_log_parms = FALSE, verbose = FALSE, ...) { - if (packageVersion("saemix") < "3.1.9000") { - stop("To use the interface to saemix, you need to install a development version\n", - "preferably https://github.com/jranke/saemixextension@warp_combined") + if (packageVersion("saemix") < "3.0") { + stop("To use the interface to saemix, you need to install a version >= 3.0\n") } if (nrow(object) > 1) stop("Only row objects allowed") diff --git a/test.log b/test.log index f68ec45a..4d5012a0 100644 --- a/test.log +++ b/test.log @@ -1,42 +1,42 @@ ℹ Loading mkin Loading required package: parallel ℹ Testing mkin -✔ | OK F W S | Context -✔ | 5 | AIC calculation -✔ | 5 | Analytical solutions for coupled models [3.4 s] -✔ | 5 | Calculation of Akaike weights -✔ | 2 | Export dataset for reading into CAKE -✔ | 12 | Confidence intervals and p-values [1.0 s] -✔ | 14 | Error model fitting [4.6 s] -✔ | 5 | Time step normalisation -✔ | 4 | Calculation of FOCUS chi2 error levels [0.5 s] -✔ | 14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.8 s] -✔ | 4 | Test fitting the decline of metabolites from their maximum [0.3 s] -✔ | 1 | Fitting the logistic model [0.2 s] -✔ | 35 1 | Nonlinear mixed-effects models [26.3 s] +✔ | F W S OK | Context +✔ | 5 | AIC calculation +✔ | 5 | Analytical solutions for coupled models [3.5s] +✔ | 5 | Calculation of Akaike weights +✔ | 2 | Export dataset for reading into CAKE +✔ | 12 | Confidence intervals and p-values [1.0s] +✔ | 14 | Error model fitting [5.0s] +✔ | 5 | Time step normalisation +✔ | 4 | Calculation of FOCUS chi2 error levels [0.6s] +✔ | 14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.8s] +✔ | 4 | Test fitting the decline of metabolites from their maximum [0.3s] +✔ | 1 | Fitting the logistic model [0.2s] +✔ | 1 35 | Nonlinear mixed-effects models [26.8s] ──────────────────────────────────────────────────────────────────────────────── Skip (test_mixed.R:161:3): saem results are reproducible for biphasic fits Reason: Fitting with saemix takes around 10 minutes when using deSolve ──────────────────────────────────────────────────────────────────────────────── -✔ | 2 | Test dataset classes mkinds and mkindsg -✔ | 10 | Special cases of mkinfit calls [0.4 s] -✔ | 1 | mkinfit features [0.3 s] -✔ | 8 | mkinmod model generation and printing [0.2 s] -✔ | 3 | Model predictions with mkinpredict [0.3 s] -✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.4 s] -✔ | 9 | Nonlinear mixed-effects models with nlme [7.9 s] -✔ | 16 | Plotting [1.3 s] -✔ | 4 | Residuals extracted from mkinfit models -✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.5 s] -✔ | 7 | Fitting the SFORB model [3.7 s] -✔ | 1 | Summaries of old mkinfit objects -✔ | 4 | Summary [0.1 s] -✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.3 s] -✔ | 9 | Hypothesis tests [8.4 s] -✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.3 s] +✔ | 2 | Test dataset classes mkinds and mkindsg +✔ | 10 | Special cases of mkinfit calls [0.4s] +✔ | 1 | mkinfit features [0.4s] +✔ | 8 | mkinmod model generation and printing [0.2s] +✔ | 3 | Model predictions with mkinpredict [0.3s] +✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.5s] +✔ | 9 | Nonlinear mixed-effects models with nlme [8.2s] +✔ | 16 | Plotting [1.3s] +✔ | 4 | Residuals extracted from mkinfit models +✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.5s] +✔ | 7 | Fitting the SFORB model [3.8s] +✔ | 1 | Summaries of old mkinfit objects +✔ | 4 | Summary [0.1s] +✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.3s] +✔ | 9 | Hypothesis tests [8.6s] +✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 67.6 s +Duration: 69.5 s ── Skipped tests ────────────────────────────────────────────────────────────── • Fitting with saemix takes around 10 minutes when using deSolve (1) diff --git a/tests/testthat/_snaps/plot/mixed-model-fit-for-saem-object-with-mkin-transformations.svg b/tests/testthat/_snaps/plot/mixed-model-fit-for-saem-object-with-mkin-transformations.svg index d3ca239b..6346a383 100644 --- a/tests/testthat/_snaps/plot/mixed-model-fit-for-saem-object-with-mkin-transformations.svg +++ b/tests/testthat/_snaps/plot/mixed-model-fit-for-saem-object-with-mkin-transformations.svg @@ -96,7 +96,7 @@ - + @@ -157,7 +157,7 @@ - + @@ -176,7 +176,7 @@ - + @@ -213,7 +213,7 @@ - + @@ -250,7 +250,7 @@ - + @@ -269,7 +269,7 @@ - + @@ -288,7 +288,7 @@ - + @@ -343,7 +343,7 @@ - + @@ -416,7 +416,7 @@ - + @@ -471,7 +471,7 @@ - + @@ -526,7 +526,7 @@ - + @@ -563,7 +563,7 @@ - + @@ -618,7 +618,7 @@ - + @@ -673,7 +673,7 @@ - + @@ -710,7 +710,7 @@ - + @@ -729,7 +729,7 @@ - + @@ -739,30 +739,30 @@ - + - - - - - + + + + + 0 -20 -40 -60 -80 -100 - - - +20 +40 +60 +80 +100 + + + - - --4 --2 + + +-4 +-2 0 -2 -4 +2 +4 @@ -776,582 +776,582 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -1359,7 +1359,7 @@ - + @@ -1416,7 +1416,7 @@ - + @@ -1431,7 +1431,7 @@ - + @@ -1464,7 +1464,7 @@ - + @@ -1497,7 +1497,7 @@ - + @@ -1514,7 +1514,7 @@ - + @@ -1530,7 +1530,7 @@ - + @@ -1579,7 +1579,7 @@ - + @@ -1644,7 +1644,7 @@ - + @@ -1693,7 +1693,7 @@ - + @@ -1742,7 +1742,7 @@ - + @@ -1775,7 +1775,7 @@ - + @@ -1824,7 +1824,7 @@ - + @@ -1873,7 +1873,7 @@ - + @@ -1906,7 +1906,7 @@ - + @@ -1923,7 +1923,7 @@ - + @@ -1933,28 +1933,28 @@ - + - - - - + + + + 0 -10 -20 -30 -40 - - - +10 +20 +30 +40 + + + - - --4 --2 + + +-4 +-2 0 -2 -4 +2 +4 @@ -1968,515 +1968,515 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/testthat/_snaps/plot/mixed-model-fit-for-saem-object-with-saemix-transformations.svg b/tests/testthat/_snaps/plot/mixed-model-fit-for-saem-object-with-saemix-transformations.svg index 072154ee..13590b9b 100644 --- a/tests/testthat/_snaps/plot/mixed-model-fit-for-saem-object-with-saemix-transformations.svg +++ b/tests/testthat/_snaps/plot/mixed-model-fit-for-saem-object-with-saemix-transformations.svg @@ -51,7 +51,7 @@ - + @@ -108,7 +108,7 @@ - + @@ -127,7 +127,7 @@ - + @@ -160,7 +160,7 @@ - + @@ -201,7 +201,7 @@ - + @@ -218,7 +218,7 @@ - + @@ -228,30 +228,30 @@ - + - + - - + + 0 -20 +20 40 60 -80 -100 - - - +80 +100 + + + - - --4 --2 + + +-4 +-2 0 -2 -4 +2 +4 @@ -265,132 +265,132 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -398,7 +398,7 @@ - + @@ -453,7 +453,7 @@ - + @@ -470,7 +470,7 @@ - + @@ -499,7 +499,7 @@ - + @@ -536,7 +536,7 @@ - + @@ -551,7 +551,7 @@ - + @@ -561,30 +561,30 @@ - + - - - - - + + + + + 0 -5 -10 -15 -20 -25 - - - +5 +10 +15 +20 +25 + + + - - --4 --2 + + +-4 +-2 0 -2 -4 +2 +4 @@ -598,118 +598,118 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/testthat/print_sfo_saem_1.txt b/tests/testthat/print_sfo_saem_1.txt index 0c0e32ce..3fc9ca3b 100644 --- a/tests/testthat/print_sfo_saem_1.txt +++ b/tests/testthat/print_sfo_saem_1.txt @@ -7,15 +7,15 @@ Data: Likelihood computed by importance sampling AIC BIC logLik - 1310 1315 -649 + 1311 1315 -649 Fitted parameters: estimate lower upper -parent_0 1e+02 98.87 1e+02 +parent_0 1e+02 98.96 1e+02 k_parent 4e-02 0.03 4e-02 -Var.parent_0 1e+00 -1.72 5e+00 +Var.parent_0 8e-01 -1.94 3e+00 Var.k_parent 1e-01 0.03 2e-01 a.1 9e-01 0.75 1e+00 b.1 5e-02 0.04 5e-02 -SD.parent_0 1e+00 -0.12 3e+00 +SD.parent_0 9e-01 -0.67 2e+00 SD.k_parent 3e-01 0.20 4e-01 diff --git a/tests/testthat/setup_script.R b/tests/testthat/setup_script.R index 96e865d4..cb3713aa 100644 --- a/tests/testthat/setup_script.R +++ b/tests/testthat/setup_script.R @@ -185,7 +185,7 @@ ds_biphasic <- lapply(ds_biphasic_mean, function(ds) { # Mixed model fits saemix_available <- FALSE if (requireNamespace("saemix", quietly = TRUE)) { - if(packageVersion("saemix") >= "3.1.9000") saemix_available <- TRUE + if(packageVersion("saemix") >= "3.0") saemix_available <- TRUE } mmkin_sfo_1 <- mmkin("SFO", ds_sfo, quiet = TRUE, error_model = "tc", cores = n_cores) mmkin_dfop_1 <- mmkin("DFOP", ds_dfop, quiet = TRUE, cores = n_cores) diff --git a/tests/testthat/summary_saem_biphasic_s.txt b/tests/testthat/summary_saem_biphasic_s.txt index 8dfae367..bab4bf98 100644 --- a/tests/testthat/summary_saem_biphasic_s.txt +++ b/tests/testthat/summary_saem_biphasic_s.txt @@ -34,33 +34,33 @@ Results: Likelihood computed by importance sampling AIC BIC logLik - 2702 2711 -1338 + 2679 2689 -1327 Optimised parameters: est. lower upper parent_0 1.0e+02 1.0e+02 1.0e+02 -k_m1 4.7e-03 3.9e-03 5.6e-03 +k_m1 4.8e-03 4.1e-03 5.5e-03 f_parent_to_m1 4.8e-01 4.3e-01 5.2e-01 -k1 4.8e-02 3.1e-02 6.5e-02 -k2 1.3e-02 8.7e-03 1.7e-02 -g 5.0e-01 4.1e-01 5.8e-01 +k1 5.9e-02 4.6e-02 7.2e-02 +k2 1.1e-02 9.0e-03 1.3e-02 +g 4.9e-01 4.3e-01 5.4e-01 Correlation: prnt_0 k_m1 f_p__1 k1 k2 -k_m1 -0.152 -f_parent_to_m1 -0.143 0.366 -k1 0.097 -0.014 -0.021 -k2 0.022 0.083 0.023 0.101 -g -0.084 -0.144 -0.044 -0.303 -0.364 +k_m1 -0.168 +f_parent_to_m1 -0.141 0.379 +k1 0.139 -0.004 -0.024 +k2 0.055 0.154 0.033 0.246 +g -0.078 -0.206 -0.058 -0.435 -0.601 Random effects: - est. lower upper -SD.parent_0 1.22 0.316 2.12 -SD.k_m1 0.15 -0.079 0.38 -SD.f_parent_to_m1 0.32 0.191 0.44 -SD.k1 0.66 0.416 0.90 -SD.k2 0.59 0.368 0.80 -SD.g 0.16 -0.373 0.70 + est. lower upper +SD.parent_0 1.1986 0.28 2.12 +SD.k_m1 0.0034 -6.85 6.86 +SD.f_parent_to_m1 0.3369 0.21 0.46 +SD.k1 0.3790 0.24 0.52 +SD.k2 0.2666 0.16 0.37 +SD.g 0.0401 -0.67 0.75 Variance model: est. lower upper @@ -73,5 +73,5 @@ parent_sink 0.52 Estimated disappearance times: DT50 DT90 DT50back DT50_k1 DT50_k2 -parent 26 127 38 14 54 -m1 146 485 NA NA NA +parent 25 150 45 12 64 +m1 145 483 NA NA NA diff --git a/vignettes/FOCUS_D.html b/vignettes/FOCUS_D.html index ba514c18..d9127473 100644 --- a/vignettes/FOCUS_D.html +++ b/vignettes/FOCUS_D.html @@ -360,7 +360,7 @@ pre code {

    Example evaluation of FOCUS Example Dataset D

    Johannes Ranke

    -

    Last change 31 January 2019 (rebuilt 2021-09-16)

    +

    Last change 31 January 2019 (rebuilt 2021-11-17)

    @@ -435,9 +435,9 @@ print(FOCUS_2006_D)

    A comprehensive report of the results is obtained using the summary method for mkinfit objects.

    summary(fit)
    ## mkin version used for fitting:    1.1.0 
-## R version used for fitting:       4.1.1 
-## Date of fit:     Thu Sep 16 13:57:32 2021 
-## Date of summary: Thu Sep 16 13:57:33 2021 
+## R version used for fitting:       4.1.2 
+## Date of fit:     Wed Nov 17 12:15:48 2021 
+## Date of summary: Wed Nov 17 12:15:49 2021 
 ## 
 ## Equations:
 ## d_parent/dt = - k_parent * parent
diff --git a/vignettes/FOCUS_L.html b/vignettes/FOCUS_L.html
index b6ebb606..96a823cf 100644
--- a/vignettes/FOCUS_L.html
+++ b/vignettes/FOCUS_L.html
@@ -1513,7 +1513,7 @@ div.tocify {
 
 
    Example evaluation of FOCUS Laboratory Data L1 to L3
    
 
    Johannes Ranke
    
-
    Last change 17 November 2016 (rebuilt 2021-09-16)
    
+
    Last change 17 November 2016 (rebuilt 2021-11-17)
    @@ -1533,16 +1533,16 @@ FOCUS_2006_L1_mkin <- mkin_wide_to_long(FOCUS_2006_L1) 
    m.L1.SFO <- mkinfit("SFO", FOCUS_2006_L1_mkin, quiet = TRUE)
 summary(m.L1.SFO)
    ## mkin version used for fitting:    1.1.0 
-## R version used for fitting:       4.1.1 
-## Date of fit:     Thu Sep 16 13:57:35 2021 
-## Date of summary: Thu Sep 16 13:57:35 2021 
+## R version used for fitting:       4.1.2 
+## Date of fit:     Wed Nov 17 12:15:51 2021 
+## Date of summary: Wed Nov 17 12:15:51 2021 
 ## 
 ## Equations:
 ## d_parent/dt = - k_parent * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted using 133 model solutions performed in 0.031 s
+## Fitted using 133 model solutions performed in 0.032 s
 ## 
 ## Error model: Constant variance 
 ## 
@@ -1634,9 +1634,9 @@ summary(m.L1.SFO)
    ## Warning in cov2cor(ans$covar): diag(.) had 0 or NA entries; non-finite result is
 ## doubtful
    ## mkin version used for fitting:    1.1.0 
-## R version used for fitting:       4.1.1 
-## Date of fit:     Thu Sep 16 13:57:35 2021 
-## Date of summary: Thu Sep 16 13:57:35 2021 
+## R version used for fitting:       4.1.2 
+## Date of fit:     Wed Nov 17 12:15:51 2021 
+## Date of summary: Wed Nov 17 12:15:51 2021 
 ## 
 ## Equations:
 ## d_parent/dt = - (alpha/beta) * 1/((time/beta) + 1) * parent
@@ -1739,16 +1739,16 @@ plot(m.L2.FOMC, show_residuals = TRUE,
 
    
 
    summary(m.L2.FOMC, data = FALSE)
    
 
    ## mkin version used for fitting:    1.1.0 
-## R version used for fitting:       4.1.1 
-## Date of fit:     Thu Sep 16 13:57:35 2021 
-## Date of summary: Thu Sep 16 13:57:35 2021 
+## R version used for fitting:       4.1.2 
+## Date of fit:     Wed Nov 17 12:15:52 2021 
+## Date of summary: Wed Nov 17 12:15:52 2021 
 ## 
 ## Equations:
 ## d_parent/dt = - (alpha/beta) * 1/((time/beta) + 1) * parent
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted using 239 model solutions performed in 0.048 s
+## Fitted using 239 model solutions performed in 0.049 s
 ## 
 ## Error model: Constant variance 
 ## 
@@ -1817,9 +1817,9 @@ plot(m.L2.DFOP, show_residuals = TRUE, show_errmin = TRUE,
 
    
 
    summary(m.L2.DFOP, data = FALSE)
    
 
    ## mkin version used for fitting:    1.1.0 
-## R version used for fitting:       4.1.1 
-## Date of fit:     Thu Sep 16 13:57:36 2021 
-## Date of summary: Thu Sep 16 13:57:36 2021 
+## R version used for fitting:       4.1.2 
+## Date of fit:     Wed Nov 17 12:15:52 2021 
+## Date of summary: Wed Nov 17 12:15:52 2021 
 ## 
 ## Equations:
 ## d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -1828,7 +1828,7 @@ plot(m.L2.DFOP, show_residuals = TRUE, show_errmin = TRUE,
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted using 581 model solutions performed in 0.133 s
+## Fitted using 581 model solutions performed in 0.134 s
 ## 
 ## Error model: Constant variance 
 ## 
@@ -1917,9 +1917,9 @@ plot(mm.L3)
    
 
    We can extract the summary and plot for e.g. the DFOP fit, using square brackets for indexing which will result in the use of the summary and plot functions working on mkinfit objects.
    
 
    summary(mm.L3[["DFOP", 1]])
    
 
    ## mkin version used for fitting:    1.1.0 
-## R version used for fitting:       4.1.1 
-## Date of fit:     Thu Sep 16 13:57:36 2021 
-## Date of summary: Thu Sep 16 13:57:36 2021 
+## R version used for fitting:       4.1.2 
+## Date of fit:     Wed Nov 17 12:15:52 2021 
+## Date of summary: Wed Nov 17 12:15:52 2021 
 ## 
 ## Equations:
 ## d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 *
@@ -1928,7 +1928,7 @@ plot(mm.L3)
    
 ## 
 ## Model predictions using solution type analytical 
 ## 
-## Fitted using 376 model solutions performed in 0.079 s
+## Fitted using 376 model solutions performed in 0.08 s
 ## 
 ## Error model: Constant variance 
 ## 
@@ -2025,9 +2025,9 @@ plot(mm.L4)
    
 
    The χ2 error level of 3.3% as well as the plot suggest that the SFO model fits very well. The error level at which the χ2 test passes is slightly lower for the FOMC model. However, the difference appears negligible.
    
 
    summary(mm.L4[["SFO", 1]], data = FALSE)
    
 
    ## mkin version used for fitting:    1.1.0 
-## R version used for fitting:       4.1.1 
-## Date of fit:     Thu Sep 16 13:57:36 2021 
-## Date of summary: Thu Sep 16 13:57:36 2021 
+## R version used for fitting:       4.1.2 
+## Date of fit:     Wed Nov 17 12:15:53 2021 
+## Date of summary: Wed Nov 17 12:15:53 2021 
 ## 
 ## Equations:
 ## d_parent/dt = - k_parent * parent
@@ -2089,9 +2089,9 @@ plot(mm.L4)
    
 ## parent  106  352
    summary(mm.L4[["FOMC", 1]], data = FALSE)
    ## mkin version used for fitting:    1.1.0 
-## R version used for fitting:       4.1.1 
-## Date of fit:     Thu Sep 16 13:57:36 2021 
-## Date of summary: Thu Sep 16 13:57:36 2021 
+## R version used for fitting:       4.1.2 
+## Date of fit:     Wed Nov 17 12:15:53 2021 
+## Date of summary: Wed Nov 17 12:15:53 2021 
 ## 
 ## Equations:
 ## d_parent/dt = - (alpha/beta) * 1/((time/beta) + 1) * parent
-- 
cgit v1.2.3