diff options
| -rw-r--r-- | NAMESPACE | 13 | ||||
| -rw-r--r-- | NEWS.md | 8 | ||||
| -rw-r--r-- | R/illparms.R | 70 | ||||
| -rw-r--r-- | R/mhmkin.R | 184 | ||||
| -rw-r--r-- | R/saem.R | 34 | ||||
| -rw-r--r-- | log/test.log | 21 | ||||
| -rw-r--r-- | man/illparms.Rd | 35 | ||||
| -rw-r--r-- | man/mhmkin.Rd | 75 | ||||
| -rw-r--r-- | man/saem.Rd | 9 | ||||
| -rw-r--r-- | tests/testthat/convergence_fits_synth_const.txt | 0 | ||||
| -rw-r--r-- | tests/testthat/convergence_hfits_synth.txt | 0 | ||||
| -rw-r--r-- | tests/testthat/illparms_hfits_synth.txt | 0 | ||||
| -rw-r--r-- | tests/testthat/setup_script.R | 4 | ||||
| -rw-r--r-- | tests/testthat/test_mhmkin.R | 31 |
14 files changed, 460 insertions, 24 deletions
@@ -1,21 +1,29 @@ # Generated by roxygen2: do not edit by hand +S3method("[",mhmkin) S3method("[",mmkin) +S3method(AIC,mhmkin) S3method(AIC,mmkin) +S3method(BIC,mhmkin) S3method(BIC,mmkin) S3method(aw,mkinfit) S3method(aw,mmkin) S3method(confint,mkinfit) +S3method(convergence,mhmkin) S3method(convergence,mmkin) S3method(f_time_norm_focus,mkindsg) S3method(f_time_norm_focus,numeric) +S3method(illparms,mhmkin) S3method(illparms,mkinfit) S3method(illparms,mmkin) +S3method(illparms,saem.mmkin) S3method(intervals,saem.mmkin) S3method(loftest,mkinfit) S3method(logLik,mkinfit) +S3method(logLik,saem.mmkin) S3method(lrtest,mkinfit) S3method(lrtest,mmkin) +S3method(mhmkin,list) S3method(mixed,mmkin) S3method(mkinpredict,mkinfit) S3method(mkinpredict,mkinmod) @@ -27,8 +35,11 @@ S3method(plot,mixed.mmkin) S3method(plot,mkinfit) S3method(plot,mmkin) S3method(plot,nafta) +S3method(print,convergence.mhmkin) S3method(print,convergence.mmkin) +S3method(print,illparms.mhmkin) S3method(print,illparms.mmkin) +S3method(print,mhmkin) S3method(print,mixed.mmkin) S3method(print,mkinds) S3method(print,mkindsg) @@ -50,6 +61,7 @@ S3method(summary,saem.mmkin) S3method(update,mkinfit) S3method(update,mmkin) S3method(update,nlme.mmkin) +S3method(update,saem.mmkin) export(CAKE_export) export(DFOP.solution) export(FOMC.solution) @@ -78,6 +90,7 @@ export(max_twa_hs) export(max_twa_parent) export(max_twa_sfo) export(mean_degparms) +export(mhmkin) export(mixed) export(mkin_long_to_wide) export(mkin_wide_to_long) @@ -1,10 +1,14 @@ # mkin 1.1.2 +- 'R/mhmkin.R': New method for performing multiple hierarchical mkin fits in one function call, optionally in parallel. + - 'R/saem.R': Implement and test saemix transformations for FOMC and HS. Also, error out if saemix transformations are requested but not supported. -- 'R/convergence.R': New generic to show convergence information with a method for 'mmin' objects. +- 'R/saem.R': 'logLik' and 'update' methods for 'saem.mmkin' objects. + +- 'R/convergence.R': New generic to show convergence information with methods for 'mmkin' and 'mhmkin' objects. -- 'R/illparms.R': New generic to show ill-defined parameters with methods for 'mkinfit' and 'mmkin' objects. +- 'R/illparms.R': New generic to show ill-defined parameters with methods for 'mkinfit', 'mmkin', 'saem.mmkin' and 'mhmkin' objects. - 'R/summary.mmkin.R': Summary method for mmkin objects. diff --git a/R/illparms.R b/R/illparms.R index f23f1cae..e4b28c56 100644 --- a/R/illparms.R +++ b/R/illparms.R @@ -1,12 +1,29 @@ #' Method to get the names of ill-defined parameters #' +#' The method for generalised nonlinear regression fits as obtained +#' with [mkinfit] and [mmkin] checks if the degradation parameters +#' pass the Wald test (in degradation kinetics often simply called t-test) for +#' significant difference from zero. For this test, the parameterisation +#' without parameter transformations is used. +#' +#' The method for hierarchical model fits, also known as nonlinear +#' mixed-effects model fits as obtained with [saem] and [mhmkin] +#' checks if any of the confidence intervals for the random +#' effects expressed as standard deviations include zero, and if +#' the confidence intervals for the error model parameters include +#' zero. +#' #' @param object The object to investigate #' @param x The object to be printed #' @param conf.level The confidence level for checking p values #' @param \dots For potential future extensions -#' @return For [mkinfit] objects, a character vector of parameter names -#' For [mmkin] objects, an object of class 'illparms.mmkin' with a -#' suitable printing method. +#' @param random For hierarchical fits, should random effects be tested? +#' @param errmod For hierarchical fits, should error model parameters be +#' tested? +#' @return For [mkinfit] or [saem] objects, a character vector of parameter +#' names. For [mmkin] or [mhmkin] objects, a matrix like object of class +#' 'illparms.mmkin' or 'illparms.mhmkin'. The latter objects have a suitable +#' printing method. #' @export illparms <- function(object, ...) { @@ -60,3 +77,50 @@ print.illparms.mmkin <- function(x, ...) { class(x) <- NULL print(x, quote = FALSE) } + +#' @rdname illparms +#' @export +illparms.saem.mmkin <- function(object, conf.level = 0.95, random = TRUE, errmod = TRUE, ...) { + if (inherits(object, "try-error")) { + return(NA) + } else { + ints <- intervals(object, conf.level = conf.level) + failed <- NULL + if (random) { + failed_random <- ints$random[, "lower"] < 0 + failed <- c(failed, names(which(failed_random))) + } + if (errmod) { + failed_errmod <- ints$errmod[, "lower"] < 0 & ints$errmod[, "est."] > 0 + failed <- c(failed, names(which(failed_errmod))) + } + } + return(failed) +} + +#' @rdname illparms +#' @export +illparms.mhmkin <- function(object, conf.level = 0.95, random = TRUE, errmod = TRUE, ...) { + result <- lapply(object, + function(fit) { + if (inherits(fit, "try-error")) return("E") + ill <- illparms(fit, conf.level = conf.level, random = random, errmod = errmod) + if (length(ill) > 0) { + return(paste(ill, collapse = ", ")) + } else { + return("") + } + }) + result <- unlist(result) + dim(result) <- dim(object) + dimnames(result) <- dimnames(object) + class(result) <- "illparms.mhmkin" + return(result) +} + +#' @rdname illparms +#' @export +print.illparms.mhmkin <- function(x, ...) { + class(x) <- NULL + print(x, quote = FALSE) +} diff --git a/R/mhmkin.R b/R/mhmkin.R new file mode 100644 index 00000000..61d67204 --- /dev/null +++ b/R/mhmkin.R @@ -0,0 +1,184 @@ +#' Fit nonlinear mixed-effects models built from one or more kinetic +#' degradation models and one or more error models +#' +#' The name of the methods expresses that (**m**ultiple) **h**ierarchichal +#' (also known as multilevel) **m**ulticompartment **kin**etic models are +#' fitted. Our kinetic models are nonlinear, so we can use various nonlinear +#' mixed-effects model fitting functions. +#' +#' @param objects A list of [mmkin] objects containing fits of the same +#' degradation models to the same data, but using different error models. +#' @param backend The backend to be used for fitting. Currently, only saemix is +#' supported +#' @param algorithm The algorithm to be used for fitting (currently not used) +#' @param \dots Further arguments that will be passed to the nonlinear mixed-effects +#' model fitting function. +#' @param cores The number of cores to be used for multicore processing. This +#' is only used when the \code{cluster} argument is \code{NULL}. On Windows +#' machines, cores > 1 is not supported, you need to use the \code{cluster} +#' argument to use multiple logical processors. Per default, all cores detected +#' by [parallel::detectCores()] are used, except on Windows where the default +#' is 1. +#' @param cluster A cluster as returned by [makeCluster] to be used for +#' parallel execution. +#' @importFrom parallel mclapply parLapply detectCores +#' @return A two-dimensional [array] of fit objects and/or try-errors that can +#' be indexed using the degradation model names for the first index (row index) +#' and the error model names for the second index (column index), with class +#' attribute 'mhmkin'. +#' @author Johannes Ranke +#' @seealso \code{\link{[.mhmkin}} for subsetting [mhmkin] objects +#' @export +mhmkin <- function(objects, backend = "saemix", algorithm = "saem", ...) { + UseMethod("mhmkin") +} + +#' @export +#' @rdname mhmkin +mhmkin.list <- function(objects, backend = "saemix", + ..., + cores = if (Sys.info()["sysname"] == "Windows") 1 else parallel::detectCores(), cluster = NULL) +{ + call <- match.call() + dot_args <- list(...) + backend_function <- switch(backend, + saemix = "saem" + ) + + deg_models <- lapply(objects[[1]][, 1], function(x) x$mkinmod) + names(deg_models) <- dimnames(objects[[1]])$model + n.deg <- length(deg_models) + + ds <- lapply(objects[[1]][1, ], function(x) x$data) + + for (other in objects[-1]) { + # Check if the degradation models in all objects are the same + for (deg_model_name in names(deg_models)) { + if (!all.equal(other[[deg_model_name, 1]]$mkinmod$spec, + deg_models[[deg_model_name]]$spec)) + { + stop("The mmkin objects have to be based on the same degradation models") + } + } + # Check if they have been fitted to the same dataset + other_object_ds <- lapply(other[1, ], function(x) x$data) + for (i in 1:length(ds)) { + if (!all.equal(ds[[i]][c("time", "variable", "observed")], + other_object_ds[[i]][c("time", "variable", "observed")])) + { + stop("The mmkin objects have to be fitted to the same datasets") + } + } + } + + n.o <- length(objects) + + error_models = sapply(objects, function(x) x[[1]]$err_mod) + n.e <- length(error_models) + + n.fits <- n.deg * n.e + fit_indices <- matrix(1:n.fits, ncol = n.e) + dimnames(fit_indices) <- list(degradation = names(deg_models), + error = error_models) + + fit_function <- function(fit_index) { + w <- which(fit_indices == fit_index, arr.ind = TRUE) + deg_index <- w[1] + error_index <- w[2] + mmkin_row <- objects[[error_index]][deg_index, ] + res <- try(do.call(backend_function, args = c(list(mmkin_row), dot_args))) + return(res) + } + + fit_time <- system.time({ + if (is.null(cluster)) { + results <- parallel::mclapply(as.list(1:n.fits), fit_function, + mc.cores = cores, mc.preschedule = FALSE) + } else { + results <- parallel::parLapply(cluster, as.list(1:n.fits), fit_function) + } + }) + + attributes(results) <- attributes(fit_indices) + attr(results, "call") <- call + attr(results, "time") <- fit_time + class(results) <- "mhmkin" + return(results) +} + +#' Subsetting method for mhmkin objects +#' +#' @param x An [mhmkin] object. +#' @param i Row index selecting the fits for specific models +#' @param j Column index selecting the fits to specific datasets +#' @param drop If FALSE, the method always returns an mhmkin object, otherwise +#' either a list of fit objects or a single fit object. +#' @return An object of class \code{\link{mhmkin}}. +#' @rdname mhmkin +#' @export +`[.mhmkin` <- function(x, i, j, ..., drop = FALSE) { + class(x) <- NULL + x_sub <- x[i, j, drop = drop] + + if (!drop) { + class(x_sub) <- "mhmkin" + } + return(x_sub) +} + +#' Print method for mhmkin objects +#' +#' @rdname mhmkin +#' @export +print.mhmkin <- function(x, ...) { + cat("<mhmkin> object\n") + cat("Status of individual fits:\n\n") + print(convergence(x)) +} + +#' @export +convergence.mhmkin <- function(object, ...) { + all_summary_warnings <- character() + + result <- lapply(object, + function(fit) { + if (inherits(fit, "try-error")) return("E") + else { + return("OK") + } + }) + result <- unlist(result) + dim(result) <- dim(object) + dimnames(result) <- dimnames(object) + + class(result) <- "convergence.mhmkin" + return(result) +} + +#' @export +print.convergence.mhmkin <- function(x, ...) { + class(x) <- NULL + print(x, quote = FALSE) + cat("\n") + if (any(x == "OK")) cat("OK: Fit terminated successfully\n") + if (any(x == "E")) cat("E: Error\n") +} + +#' @export +AIC.mhmkin <- function(object, ..., k = 2) { + res <- sapply(object, function(x) { + if (inherits(x, "try-error")) return(NA) + else return(AIC(x$so, k = k)) + }) + dim(res) <- dim(object) + dimnames(res) <- dimnames(object) + return(res) +} + +#' @export +BIC.mhmkin <- function(object, ...) { + res <- sapply(object, function(x) BIC(x$so)) + dim(res) <- dim(object) + dimnames(res) <- dimnames(object) + return(res) +} @@ -58,6 +58,9 @@ utils::globalVariables(c("predicted", "std")) #' f_saem_sfo <- saem(f_mmkin_parent["SFO", ]) #' f_saem_fomc <- saem(f_mmkin_parent["FOMC", ]) #' f_saem_dfop <- saem(f_mmkin_parent["DFOP", ]) +#' illparms(f_saem_dfop) +#' update(f_saem_dfop, covariance.model = diag(c(1, 1, 1, 0))) +#' AIC(f_saem_dfop) #' #' # The returned saem.mmkin object contains an SaemixObject, therefore we can use #' # functions from saemix @@ -125,6 +128,7 @@ saem.mmkin <- function(object, fail_with_errors = TRUE, verbose = FALSE, quiet = FALSE, ...) { + call <- match.call() transformations <- match.arg(transformations) m_saemix <- saemix_model(object, verbose = verbose, degparms_start = degparms_start, @@ -184,6 +188,7 @@ saem.mmkin <- function(object, transform_rates = object[[1]]$transform_rates, transform_fractions = object[[1]]$transform_fractions, so = f_saemix, + call = call, time = fit_time, mean_dp_start = attr(m_saemix, "mean_dp_start"), bparms.optim = bparms_optim, @@ -579,3 +584,32 @@ saemix_data <- function(object, verbose = FALSE, ...) { ...) return(res) } + +#' @export +logLik.saem.mmkin <- function(object, ...) return(logLik(object$so)) + +#' @export +update.saem.mmkin <- function(object, ..., evaluate = TRUE) { + call <- object$call + # For some reason we get saem.mmkin in the call when using mhmkin + # so we need to fix this in order to avoid exporting saem.mmkin + # in addition to the S3 method + call[[1]] <- saem + update_arguments <- match.call(expand.dots = FALSE)$... + + if (length(update_arguments) > 0) { + update_arguments_in_call <- !is.na(match(names(update_arguments), names(call))) + } + + for (a in names(update_arguments)[update_arguments_in_call]) { + call[[a]] <- update_arguments[[a]] + } + + update_arguments_not_in_call <- !update_arguments_in_call + if(any(update_arguments_not_in_call)) { + call <- c(as.list(call), update_arguments[update_arguments_not_in_call]) + call <- as.call(call) + } + if(evaluate) eval(call, parent.frame()) + else call +} diff --git a/log/test.log b/log/test.log index 67ba8841..8be4a512 100644 --- a/log/test.log +++ b/log/test.log @@ -1,11 +1,11 @@ ℹ Testing mkin ✔ | F W S OK | Context ✔ | 5 | AIC calculation -✔ | 5 | Analytical solutions for coupled models [3.3s] +✔ | 5 | Analytical solutions for coupled models [3.4s] ✔ | 5 | Calculation of Akaike weights ✔ | 3 | Export dataset for reading into CAKE ✔ | 12 | Confidence intervals and p-values [1.0s] -✔ | 1 12 | Dimethenamid data from 2018 [32.9s] +✔ | 1 12 | Dimethenamid data from 2018 [32.4s] ──────────────────────────────────────────────────────────────────────────────── Skip (test_dmta.R:98:3): Different backends get consistent results for SFO-SFO3+, dimethenamid data Reason: Fitting this ODE model with saemix takes about 15 minutes on my system @@ -16,6 +16,7 @@ Reason: Fitting this ODE model with saemix takes about 15 minutes on my system ✔ | 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.4s] ✔ | 1 | Fitting the logistic model [0.2s] +✔ | 5 | Batch fitting and diagnosing hierarchical kinetic models [14.5s] ✔ | 1 12 | Nonlinear mixed-effects models [0.2s] ──────────────────────────────────────────────────────────────────────────────── Skip (test_mixed.R:68:3): saemix results are reproducible for biphasic fits @@ -27,23 +28,23 @@ Reason: Fitting with saemix takes around 10 minutes when using deSolve ✔ | 8 | mkinmod model generation and printing [0.2s] ✔ | 3 | Model predictions with mkinpredict [0.4s] ✔ | 16 | Evaluations according to 2015 NAFTA guidance [1.8s] -✔ | 9 | Nonlinear mixed-effects models with nlme [8.3s] -✔ | 16 | Plotting [10.6s] +✔ | 9 | Nonlinear mixed-effects models with nlme [8.7s] +✔ | 16 | Plotting [10.1s] ✔ | 4 | Residuals extracted from mkinfit models -✔ | 28 | saemix parent models [182.2s] +✔ | 28 | saemix parent models [181.2s] ✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4s] -✔ | 7 | Fitting the SFORB model [3.9s] +✔ | 7 | Fitting the SFORB model [3.7s] ✔ | 1 | Summaries of old mkinfit objects ✔ | 5 | Summary [0.2s] ✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.1s] -✔ | 9 | Hypothesis tests [7.7s] -✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s] +✔ | 9 | Hypothesis tests [7.8s] +✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.1s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 266.8 s +Duration: 279.6 s ── Skipped tests ────────────────────────────────────────────────────────────── • Fitting this ODE model with saemix takes about 15 minutes on my system (1) • Fitting with saemix takes around 10 minutes when using deSolve (1) -[ FAIL 0 | WARN 0 | SKIP 2 | PASS 228 ] +[ FAIL 0 | WARN 0 | SKIP 2 | PASS 233 ] diff --git a/man/illparms.Rd b/man/illparms.Rd index 7f70229d..90adf2bb 100644 --- a/man/illparms.Rd +++ b/man/illparms.Rd @@ -5,6 +5,9 @@ \alias{illparms.mkinfit} \alias{illparms.mmkin} \alias{print.illparms.mmkin} +\alias{illparms.saem.mmkin} +\alias{illparms.mhmkin} +\alias{print.illparms.mhmkin} \title{Method to get the names of ill-defined parameters} \usage{ illparms(object, ...) @@ -14,6 +17,12 @@ illparms(object, ...) \method{illparms}{mmkin}(object, conf.level = 0.95, ...) \method{print}{illparms.mmkin}(x, ...) + +\method{illparms}{saem.mmkin}(object, conf.level = 0.95, random = TRUE, errmod = TRUE, ...) + +\method{illparms}{mhmkin}(object, conf.level = 0.95, random = TRUE, errmod = TRUE, ...) + +\method{print}{illparms.mhmkin}(x, ...) } \arguments{ \item{object}{The object to investigate} @@ -23,14 +32,32 @@ illparms(object, ...) \item{conf.level}{The confidence level for checking p values} \item{x}{The object to be printed} + +\item{random}{For hierarchical fits, should random effects be tested?} + +\item{errmod}{For hierarchical fits, should error model parameters be +tested?} } \value{ -For \link{mkinfit} objects, a character vector of parameter names -For \link{mmkin} objects, an object of class 'illparms.mmkin' with a -suitable printing method. +For \link{mkinfit} or \link{saem} objects, a character vector of parameter +names. For \link{mmkin} or \link{mhmkin} objects, a matrix like object of class +'illparms.mmkin' or 'illparms.mhmkin'. The latter objects have a suitable +printing method. } \description{ -Method to get the names of ill-defined parameters +The method for generalised nonlinear regression fits as obtained +with \link{mkinfit} and \link{mmkin} checks if the degradation parameters +pass the Wald test (in degradation kinetics often simply called t-test) for +significant difference from zero. For this test, the parameterisation +without parameter transformations is used. +} +\details{ +The method for hierarchical model fits, also known as nonlinear +mixed-effects model fits as obtained with \link{saem} and \link{mhmkin} +checks if any of the confidence intervals for the random +effects expressed as standard deviations include zero, and if +the confidence intervals for the error model parameters include +zero. } \examples{ fit <- mkinfit("FOMC", FOCUS_2006_A, quiet = TRUE) diff --git a/man/mhmkin.Rd b/man/mhmkin.Rd new file mode 100644 index 00000000..0ef1599e --- /dev/null +++ b/man/mhmkin.Rd @@ -0,0 +1,75 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/mhmkin.R +\name{mhmkin} +\alias{mhmkin} +\alias{mhmkin.list} +\alias{[.mhmkin} +\alias{print.mhmkin} +\title{Fit nonlinear mixed-effects models built from one or more kinetic +degradation models and one or more error models} +\usage{ +mhmkin(objects, backend = "saemix", algorithm = "saem", ...) + +\method{mhmkin}{list}( + objects, + backend = "saemix", + ..., + cores = if (Sys.info()["sysname"] == "Windows") 1 else parallel::detectCores(), + cluster = NULL +) + +\method{[}{mhmkin}(x, i, j, ..., drop = FALSE) + +\method{print}{mhmkin}(x, ...) +} +\arguments{ +\item{objects}{A list of \link{mmkin} objects containing fits of the same +degradation models to the same data, but using different error models.} + +\item{backend}{The backend to be used for fitting. Currently, only saemix is +supported} + +\item{algorithm}{The algorithm to be used for fitting (currently not used)} + +\item{\dots}{Further arguments that will be passed to the nonlinear mixed-effects +model fitting function.} + +\item{cores}{The number of cores to be used for multicore processing. This +is only used when the \code{cluster} argument is \code{NULL}. On Windows +machines, cores > 1 is not supported, you need to use the \code{cluster} +argument to use multiple logical processors. Per default, all cores detected +by \code{\link[parallel:detectCores]{parallel::detectCores()}} are used, except on Windows where the default +is 1.} + +\item{cluster}{A cluster as returned by \link{makeCluster} to be used for +parallel execution.} + +\item{x}{An \link{mhmkin} object.} + +\item{i}{Row index selecting the fits for specific models} + +\item{j}{Column index selecting the fits to specific datasets} + +\item{drop}{If FALSE, the method always returns an mhmkin object, otherwise +either a list of fit objects or a single fit object.} +} +\value{ +A two-dimensional \link{array} of fit objects and/or try-errors that can +be indexed using the degradation model names for the first index (row index) +and the error model names for the second index (column index), with class +attribute 'mhmkin'. + +An object of class \code{\link{mhmkin}}. +} +\description{ +The name of the methods expresses that (\strong{m}ultiple) \strong{h}ierarchichal +(also known as multilevel) \strong{m}ulticompartment \strong{kin}etic models are +fitted. Our kinetic models are nonlinear, so we can use various nonlinear +mixed-effects model fitting functions. +} +\seealso{ +\code{\link{[.mhmkin}} for subsetting \link{mhmkin} objects +} +\author{ +Johannes Ranke +} diff --git a/man/saem.Rd b/man/saem.Rd index 0c066dd2..a202f52f 100644 --- a/man/saem.Rd +++ b/man/saem.Rd @@ -49,9 +49,9 @@ saemix_data(object, verbose = FALSE, ...) \item{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.} +are done in 'saemix' for the supported cases, i.e. (as of version 1.1.2) +SFO, FOMC, DFOP and HS without fixing \code{parent_0}, and SFO or DFOP with +one SFO metabolite.} \item{degparms_start}{Parameter values given as a named numeric vector will be used to override the starting values obtained from the 'mmkin' object.} @@ -120,6 +120,9 @@ 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", ]) +illparms(f_saem_dfop) +update(f_saem_dfop, covariance.model = diag(c(1, 1, 1, 0))) +AIC(f_saem_dfop) # The returned saem.mmkin object contains an SaemixObject, therefore we can use # functions from saemix diff --git a/tests/testthat/convergence_fits_synth_const.txt b/tests/testthat/convergence_fits_synth_const.txt new file mode 100644 index 00000000..e69de29b --- /dev/null +++ b/tests/testthat/convergence_fits_synth_const.txt diff --git a/tests/testthat/convergence_hfits_synth.txt b/tests/testthat/convergence_hfits_synth.txt new file mode 100644 index 00000000..e69de29b --- /dev/null +++ b/tests/testthat/convergence_hfits_synth.txt diff --git a/tests/testthat/illparms_hfits_synth.txt b/tests/testthat/illparms_hfits_synth.txt new file mode 100644 index 00000000..e69de29b --- /dev/null +++ b/tests/testthat/illparms_hfits_synth.txt diff --git a/tests/testthat/setup_script.R b/tests/testthat/setup_script.R index 62cdf78e..af059668 100644 --- a/tests/testthat/setup_script.R +++ b/tests/testthat/setup_script.R @@ -8,8 +8,8 @@ if (identical(Sys.getenv("NOT_CRAN"), "true")) { n_cores <- 1 } -# We are only allowed one core on travis, but they also set NOT_CRAN=true -if (Sys.getenv("TRAVIS") != "") n_cores = 1 +# Use the two available cores on travis +if (Sys.getenv("TRAVIS") != "") n_cores = 2 # On Windows we would need to make a cluster first if (Sys.info()["sysname"] == "Windows") n_cores = 1 diff --git a/tests/testthat/test_mhmkin.R b/tests/testthat/test_mhmkin.R new file mode 100644 index 00000000..5243f971 --- /dev/null +++ b/tests/testthat/test_mhmkin.R @@ -0,0 +1,31 @@ +context("Batch fitting and diagnosing hierarchical kinetic models") + +test_that("Multiple hierarchical kinetic models can be fitted and diagnosed", { + + skip_on_cran() + fits_synth_const <- suppressWarnings( + mmkin(c("SFO", "FOMC"), ds_sfo[1:6], cores = n_cores, quiet = TRUE)) + + fits_synth_tc <- suppressWarnings( + update(fits_synth_const, error_model = "tc")) + + expect_known_output( + convergence(fits_synth_const), + "convergence_fits_synth_const.txt") + + hfits <- mhmkin(list(fits_synth_const, fits_synth_tc)) + + expect_known_output( + convergence(hfits), + "convergence_hfits_synth.txt") + + expect_known_output( + illparms(hfits), + "illparms_hfits_synth.txt") + + expect_equal(which.min(AIC(hfits)), 3) + + hfit_sfo_tc <- update(hfits[["SFO", "tc"]], + covariance.model = diag(c(0, 1))) + expect_equal(illparms(hfit_sfo_tc), character(0)) +}) |