diff options
| -rw-r--r-- | NAMESPACE | 3 | ||||
| -rw-r--r-- | NEWS.md | 2 | ||||
| -rw-r--r-- | R/endpoints.R | 87 | ||||
| -rw-r--r-- | R/nlme.mmkin.R | 15 | ||||
| -rw-r--r-- | R/saemix.R | 37 | ||||
| -rw-r--r-- | R/summary.nlme.mmkin.R | 22 | ||||
| -rw-r--r-- | R/summary.saem.mmkin.R | 250 | ||||
| -rw-r--r-- | _pkgdown.yml | 1 | ||||
| -rw-r--r-- | man/endpoints.Rd | 15 | ||||
| -rw-r--r-- | man/saem.Rd | 12 | ||||
| -rw-r--r-- | man/summary.saem.mmkin.Rd | 89 | ||||
| -rw-r--r-- | test.log | 6 | ||||
| -rw-r--r-- | tests/testthat/FOCUS_2006_D.csf | 2 |
13 files changed, 459 insertions, 82 deletions
@@ -27,10 +27,12 @@ S3method(print,nafta) S3method(print,nlme.mmkin) S3method(print,summary.mkinfit) S3method(print,summary.nlme.mmkin) +S3method(print,summary.saem.mmkin) S3method(residuals,mkinfit) S3method(saem,mmkin) S3method(summary,mkinfit) S3method(summary,nlme.mmkin) +S3method(summary,saem.mmkin) S3method(update,mkinfit) S3method(update,mmkin) S3method(update,nlme.mmkin) @@ -123,5 +125,6 @@ importFrom(stats,residuals) importFrom(stats,rnorm) importFrom(stats,shapiro.test) importFrom(stats,update) +importFrom(stats,vcov) importFrom(utils,getFromNamespace) importFrom(utils,write.table) @@ -1,5 +1,7 @@ # mkin 0.9.50.4 (unreleased) +- 'saem' generic function to fit saemix models, with a method 'saem.mmkin' and a corresponding 'summary.saem.mmkin' + - 'transform_odeparms', 'backtransform_odeparms': Use logit transformation for solitary fractions like the g parameter of the DFOP model, or formation fractions for a pathway to only one target variable - 'update' method for 'mmkin' objects diff --git a/R/endpoints.R b/R/endpoints.R index e4813db9..f1f47581 100644 --- a/R/endpoints.R +++ b/R/endpoints.R @@ -7,16 +7,22 @@ #' are equivalent to the rate constants of the DFOP model, but with the #' advantage that the SFORB model can also be used for metabolites. #' -#' @param fit An object of class \code{\link{mkinfit}} or -#' \code{\link{nlme.mmkin}} +#' Additional DT50 values are calculated from the FOMC DT90 and k1 and k2 from +#' HS and DFOP, as well as from Eigenvalues b1 and b2 of any SFORB models +#' +#' @param fit An object of class [mkinfit], [nlme.mmkin] or +#' [saem.mmkin]. Or another object that has list components +#' mkinmod containing an [mkinmod] degradation model, and two numeric vectors, +#' bparms.optim and bparms.fixed, that contain parameter values +#' for that model. #' @importFrom stats optimize #' @return A list with a matrix of dissipation times named distimes, #' and, if applicable, a vector of formation fractions named ff #' and, if the SFORB model was in use, a vector of eigenvalues #' of these SFORB models, equivalent to DFOP rate constants -#' @note The function is used internally by \code{\link{summary.mkinfit}}. +#' @note The function is used internally by [summary.mkinfit], +#' [summary.nlme.mmkin] and [summary.saem.mmkin]. #' @author Johannes Ranke -#' @keywords manip #' @examples #' #' fit <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE) @@ -30,26 +36,9 @@ #' #' @export endpoints <- function(fit) { - # Calculate dissipation times DT50 and DT90 and formation - # fractions as well as SFORB eigenvalues from optimised parameters - # Additional DT50 values are calculated from the FOMC DT90 and k1 and k2 from - # HS and DFOP, as well as from Eigenvalues b1 and b2 of any SFORB models ep <- list() - if (inherits(fit, "mkinfit")) { - mkinmod <- fit$mkinmod - parms.all <- c(fit$bparms.optim, fit$bparms.fixed) - } else { - if (inherits(fit, "nlme.mmkin")) { - mkinmod <- fit$mmkin_orig[[1]]$mkinmod - bparms.optim <- backtransform_odeparms(fit$coefficients$fixed, - mkinmod, - transform_rates = fit$mmkin_orig[[1]]$transform_rates, - transform_fractions = fit$mmkin_orig[[1]]$transform_fractions) - parms.all <- c(bparms.optim, fit$bparms.fixed) - } else { - stop("Only implemented for mkinfit and nlme.mmkin objects") - } - } + mkinmod <- fit$mkinmod + degparms <- c(fit$bparms.optim, fit$bparms.fixed) obs_vars <- names(mkinmod$spec) ep$ff <- vector() ep$SFORB <- vector() @@ -61,9 +50,9 @@ endpoints <- function(fit) { type = names(mkinmod$map[[obs_var]])[1] # Get formation fractions if directly fitted, and calculate remaining fraction to sink - f_names = grep(paste("^f", obs_var, sep = "_"), names(parms.all), value=TRUE) + f_names = grep(paste("^f", obs_var, sep = "_"), names(degparms), value=TRUE) if (length(f_names) > 0) { - f_values = parms.all[f_names] + f_values = degparms[f_names] f_to_sink = 1 - sum(f_values) names(f_to_sink) = ifelse(type == "SFORB", paste(obs_var, "free", "sink", sep = "_"), @@ -76,34 +65,34 @@ endpoints <- function(fit) { # Get the rest if (type == "SFO") { - k_names = grep(paste("^k", obs_var, sep="_"), names(parms.all), value=TRUE) - k_tot = sum(parms.all[k_names]) + k_names = grep(paste("^k", obs_var, sep="_"), names(degparms), value=TRUE) + k_tot = sum(degparms[k_names]) DT50 = log(2)/k_tot DT90 = log(10)/k_tot if (mkinmod$use_of_ff == "min" && length(obs_vars) > 1) { for (k_name in k_names) { - ep$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / k_tot + ep$ff[[sub("k_", "", k_name)]] = degparms[[k_name]] / k_tot } } } if (type == "FOMC") { - alpha = parms.all["alpha"] - beta = parms.all["beta"] + alpha = degparms["alpha"] + beta = degparms["beta"] DT50 = beta * (2^(1/alpha) - 1) DT90 = beta * (10^(1/alpha) - 1) DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011 ep$distimes[obs_var, c("DT50back")] = DT50_back } if (type == "IORE") { - k_names = grep(paste("^k__iore", obs_var, sep="_"), names(parms.all), value=TRUE) - k_tot = sum(parms.all[k_names]) + k_names = grep(paste("^k__iore", obs_var, sep="_"), names(degparms), value=TRUE) + k_tot = sum(degparms[k_names]) # From the NAFTA kinetics guidance, p. 5 - n = parms.all[paste("N", obs_var, sep = "_")] + n = degparms[paste("N", obs_var, sep = "_")] k = k_tot # Use the initial concentration of the parent compound source_name = mkinmod$map[[1]][[1]] - c0 = parms.all[paste(source_name, "0", sep = "_")] + c0 = degparms[paste(source_name, "0", sep = "_")] alpha = 1 / (n - 1) beta = (c0^(1 - n))/(k * (n - 1)) DT50 = beta * (2^(1/alpha) - 1) @@ -113,14 +102,14 @@ endpoints <- function(fit) { if (mkinmod$use_of_ff == "min") { for (k_name in k_names) { - ep$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / k_tot + ep$ff[[sub("k_", "", k_name)]] = degparms[[k_name]] / k_tot } } } if (type == "DFOP") { - k1 = parms.all["k1"] - k2 = parms.all["k2"] - g = parms.all["g"] + k1 = degparms["k1"] + k2 = degparms["k2"] + g = degparms["g"] f <- function(log_t, x) { t <- exp(log_t) fraction <- g * exp( - k1 * t) + (1 - g) * exp( - k2 * t) @@ -144,9 +133,9 @@ endpoints <- function(fit) { ep$distimes[obs_var, c("DT50_k2")] = DT50_k2 } if (type == "HS") { - k1 = parms.all["k1"] - k2 = parms.all["k2"] - tb = parms.all["tb"] + k1 = degparms["k1"] + k2 = degparms["k2"] + tb = degparms["tb"] DTx <- function(x) { DTx.a <- (log(100/(100 - x)))/k1 DTx.b <- tb + (log(100/(100 - x)) - k1 * tb)/k2 @@ -165,11 +154,11 @@ endpoints <- function(fit) { } if (type == "SFORB") { # FOCUS kinetics (2006), p. 60 f - k_out_names = grep(paste("^k", obs_var, "free", sep="_"), names(parms.all), value=TRUE) + k_out_names = grep(paste("^k", obs_var, "free", sep="_"), names(degparms), value=TRUE) k_out_names = setdiff(k_out_names, paste("k", obs_var, "free", "bound", sep="_")) - k_1output = sum(parms.all[k_out_names]) - k_12 = parms.all[paste("k", obs_var, "free", "bound", sep="_")] - k_21 = parms.all[paste("k", obs_var, "bound", "free", sep="_")] + k_1output = sum(degparms[k_out_names]) + k_12 = degparms[paste("k", obs_var, "free", "bound", sep="_")] + k_21 = degparms[paste("k", obs_var, "bound", "free", sep="_")] sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21) b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp @@ -201,7 +190,7 @@ endpoints <- function(fit) { for (k_out_name in k_out_names) { - ep$ff[[sub("k_", "", k_out_name)]] = parms.all[[k_out_name]] / k_1output + ep$ff[[sub("k_", "", k_out_name)]] = degparms[[k_out_name]] / k_1output } # Return the eigenvalues for comparison with DFOP rate constants @@ -214,9 +203,9 @@ endpoints <- function(fit) { } if (type == "logistic") { # FOCUS kinetics (2014) p. 67 - kmax = parms.all["kmax"] - k0 = parms.all["k0"] - r = parms.all["r"] + kmax = degparms["kmax"] + k0 = degparms["k0"] + r = degparms["r"] DT50 = (1/r) * log(1 - ((kmax/k0) * (1 - 2^(r/kmax)))) DT90 = (1/r) * log(1 - ((kmax/k0) * (1 - 10^(r/kmax)))) diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R index 526cb10b..af92e8a1 100644 --- a/R/nlme.mmkin.R +++ b/R/nlme.mmkin.R @@ -181,14 +181,21 @@ nlme.mmkin <- function(model, data = sys.frame(sys.parent()), fit_time <- system.time(val <- do.call("nlme.formula", thisCall)) val$time <- fit_time - val$mean_dp_start <- mean_dp_start - val$mmkin_orig <- model - val$data <- thisCall[["data"]] val$mkinmod <- model[[1]]$mkinmod - val$err_mode <- error_model + val$data <- thisCall[["data"]] + val$mmkin_orig <- model + val$mean_dp_start <- mean_dp_start val$transform_rates <- model[[1]]$transform_rates val$transform_fractions <- model[[1]]$transform_fractions val$solution_type <- model[[1]]$solution_type + val$err_mode <- error_model + + val$bparms.optim <- backtransform_odeparms(val$coefficients$fixed, + val$mkinmod, + transform_rates = val$transform_rates, + transform_fractions = val$transform_fractions) + + val$bparms.fixed <- model[[1]]$bparms.fixed val$date.fit <- date() val$nlmeversion <- as.character(utils::packageVersion("nlme")) val$mkinversion <- as.character(utils::packageVersion("mkin")) @@ -1,8 +1,8 @@ #' 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 -#' to the expectation maximisation algorithm (SAEM). +#' 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]. @@ -23,7 +23,9 @@ #' @param control Passed to [saemix::saemix] #' @param \dots Further parameters passed to [saemix::saemixData] #' and [saemix::saemixModel]. -#' @return An [saemix::SaemixObject]. +#' @return An S3 object of class 'saem.mmkin', containing the fitted +#' [saemix::SaemixObject] as a list component named 'so'. +#' @seealso [summary.saem.mmkin] #' @examples #' \dontrun{ #' ds <- lapply(experimental_data_for_UBA_2019[6:10], @@ -57,7 +59,7 @@ #' # Using a single core, the following takes about 6 minutes, using 10 cores #' # it is slower instead of faster #' f_saem_des <- saem(f_mmkin_des, cores = 1) -#' compare.saemix(list(f_saemix$so, f_saemix_des$so)) +#' compare.saemix(list(f_saem$so, f_saem_des$so)) #' } #' @export saem <- function(object, control, ...) UseMethod("saem") @@ -79,18 +81,40 @@ saem.mmkin <- function(object, tmp <- tempfile() grDevices::png(tmp) } - f_saemix <- saemix::saemix(m_saemix, d_saemix, control) + fit_time <- system.time({ + f_saemix <- saemix::saemix(m_saemix, d_saemix, control) + f_saemix <- saemix::saemix.predict(f_saemix) + }) if (suppressPlot) { grDevices::dev.off() unlink(tmp) } + transparms_optim = f_saemix@results@fixed.effects + names(transparms_optim) = f_saemix@results@name.fixed + bparms_optim <- backtransform_odeparms(transparms_optim, + object[[1]]$mkinmod, + object[[1]]$transform_rates, + object[[1]]$transform_fractions) + result <- list( mkinmod = object[[1]]$mkinmod, mmkin = object, solution_type = object[[1]]$solution_type, transform_rates = object[[1]]$transform_rates, transform_fractions = object[[1]]$transform_fractions, - so = f_saemix) + 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 = nlme_data(object), + 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) <- "saem.mmkin" return(result) } @@ -256,6 +280,7 @@ saemix_model <- function(object, cores = 1, verbose = FALSE, ...) { error.init = error.init, verbose = verbose ) + attr(res, "mean_dp_start") <- degparms_optim return(res) } diff --git a/R/summary.nlme.mmkin.R b/R/summary.nlme.mmkin.R index 9fdd3f73..7e404e00 100644 --- a/R/summary.nlme.mmkin.R +++ b/R/summary.nlme.mmkin.R @@ -75,7 +75,7 @@ summary.nlme.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes = confint_trans <- intervals(object, which = "fixed", level = 1 - alpha)$fixed attr(confint_trans, "label") <- NULL pnames <- rownames(confint_trans) - confint_trans[, "est."] + bp <- backtransform_odeparms(confint_trans[, "est."], object$mkinmod, object$transform_rates, object$transform_fractions) bpnames <- names(bp) @@ -127,14 +127,12 @@ summary.nlme.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes = object$diffs <- object$mkinmod$diffs object$print_data <- data - if (data) { - object$data[["observed"]] <- object$data[["value"]] - object$data[["value"]] <- NULL - object$data[["predicted"]] <- predict(object) - object$data[["residual"]] <- residuals(object, type = "response") - object$data[["std"]] <- object$sigma <- 1/attr(object$modelStruct$varStruct, "weights") - object$data[["standardized"]] <- residuals(object, type = "pearson") - } + object$data[["observed"]] <- object$data[["value"]] + object$data[["value"]] <- NULL + object$data[["predicted"]] <- predict(object) + object$data[["residual"]] <- residuals(object, type = "response") + object$data[["std"]] <- object$sigma <- 1/attr(object$modelStruct$varStruct, "weights") + object$data[["standardized"]] <- residuals(object, type = "pearson") object$verbose <- verbose object$fixed <- object$mmkin_orig[[1]]$fixed @@ -200,11 +198,13 @@ print.summary.nlme.mmkin <- function(x, digits = max(3, getOption("digits") - 3) print(corr, title = "\nCorrelation:", ...) } + cat("\n") # Random effects + print(summary(x$modelStruct), sigma = x$sigma, + reEstimates = x$coef$random, verbose = verbose, ...) + cat("\nBacktransformed parameters with asymmetric confidence intervals:\n") print(x$confint_back) - print(summary(x$modelStruct), sigma = x$sigma, - reEstimates = x$coef$random, verbose = verbose, ...) printSFORB <- !is.null(x$SFORB) if(printSFORB){ diff --git a/R/summary.saem.mmkin.R b/R/summary.saem.mmkin.R new file mode 100644 index 00000000..f7110dd0 --- /dev/null +++ b/R/summary.saem.mmkin.R @@ -0,0 +1,250 @@ +#' Summary method for class "saem.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 [saem.mmkin] +#' @param x an object of class [summary.saem.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 based on the [saemix::SaemixObject] +#' obtained in the fit, with at least the following additional components +#' \item{saemixversion, mkinversion, Rversion}{The saemix, 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{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 +#' saemix authors for the parts inherited from saemix. +#' @examples +#' # Generate five datasets following SFO kinetics +#' sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120) +#' dt50_sfo_in_pop <- 50 +#' k_in_pop <- log(2) / dt50_sfo_in_pop +#' set.seed(1234) +#' k_in <- rlnorm(5, log(k_in_pop), 0.5) +#' SFO <- mkinmod(parent = mkinsub("SFO")) +#' +#' pred_sfo <- function(k) { +#' mkinpredict(SFO, +#' c(k_parent = k), +#' c(parent = 100), +#' sampling_times) +#' } +#' +#' ds_sfo_mean <- lapply(k_in, pred_sfo) +#' names(ds_sfo_mean) <- paste("ds", 1:5) +#' +#' ds_sfo_syn <- lapply(ds_sfo_mean, function(ds) { +#' add_err(ds, +#' sdfunc = function(value) sqrt(1^2 + value^2 * 0.07^2), +#' n = 1)[[1]] +#' }) +#' +#' \dontrun{ +#' # Evaluate using mmkin and saem +#' f_mmkin <- mmkin("SFO", ds_sfo_syn, quiet = TRUE, error_model = "tc", cores = 1) +#' f_saem <- saem(f_mmkin) +#' summary(f_saem, data = TRUE) +#' } +#' +#' @export +summary.saem.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 <- object$so@results@conf.int + rownames(conf.int) <- conf.int$name + confint_trans <- as.matrix(conf.int[pnames, c("estimate", "lower", "upper")]) + colnames(confint_trans)[1] <- "est." + + bp <- backtransform_odeparms(confint_trans[, "est."], object$mkinmod, + object$transform_rates, object$transform_fractions) + bpnames <- names(bp) + + # Correlation of fixed effects (inspired by summary.nlme) + varFix <- vcov(object$so)[1:np, 1:np] + stdFix <- sqrt(diag(varFix)) + object$corFixed <- array( + t(varFix/stdFix)/stdFix, + dim(varFix), + list(pnames, pnames)) + + # Random effects + rnames <- paste0("SD.", pnames) + confint_ranef <- as.matrix(conf.int[rnames, c("estimate", "lower", "upper")]) + colnames(confint_ranef)[1] <- "est." + + # Error model + enames <- object$so@results@name.sigma + confint_errmod <- as.matrix(conf.int[enames, c("estimate", "lower", "upper")]) + colnames(confint_errmod)[1] <- "est." + + # 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 + } + } + + object$confint_trans <- confint_trans + object$confint_ranef <- confint_ranef + object$confint_errmod <- confint_errmod + object$confint_back <- confint_back + + object$date.summary = date() + object$use_of_ff = object$mkinmod$use_of_ff + object$error_model_algorithm = object$mmkin_orig[[1]]$error_model_algorithm + err_mod = object$mmkin_orig[[1]]$err_mod + + object$diffs <- object$mkinmod$diffs + object$print_data <- data + so_pred <- object$so@results@predictions + + object$data[["observed"]] <- object$data[["value"]] + object$data[["value"]] <- NULL + object$data[["predicted"]] <- so_pred$ipred + object$data[["residual"]] <- so_pred$ires + object$data[["standardized"]] <- so_pred$iwres + 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) + + cat("\nFixed degradation parameter values:\n") + if(length(x$fixed$value) == 0) cat("None\n") + else print(x$fixed) + + 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 = " ")) + + cat("\nOptimised, transformed parameters with symmetric confidence intervals:\n") + print(x$confint_trans) + + if (nrow(x$confint_trans) > 1) { + corr <- x$corFixed + class(corr) <- "correlation" + print(corr, title = "\nCorrelation:", ...) + } + + cat("\nRandom effects:\n") + print(x$confint_ranef) + + cat("\nVariance model:\n") + print(x$confint_errmod) + + cat("\nBacktransformed parameters with asymmetric confidence intervals:\n") + print(x$confint_back) + + 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/_pkgdown.yml b/_pkgdown.yml index 3668750f..b036ed59 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -45,6 +45,7 @@ reference: - plot.nlme.mmkin - summary.nlme.mmkin - saem.mmkin + - summary.saem.mmkin - nlme_function - get_deg_func - saemix_model diff --git a/man/endpoints.Rd b/man/endpoints.Rd index 9c354ae9..72487717 100644 --- a/man/endpoints.Rd +++ b/man/endpoints.Rd @@ -8,8 +8,11 @@ with mkinfit} endpoints(fit) } \arguments{ -\item{fit}{An object of class \code{\link{mkinfit}} or -\code{\link{nlme.mmkin}}} +\item{fit}{An object of class \link{mkinfit}, \link{nlme.mmkin} or +\link{saem.mmkin}. Or another object that has list components +mkinmod containing an \link{mkinmod} degradation model, and two numeric vectors, +bparms.optim and bparms.fixed, that contain parameter values +for that model.} } \value{ A list with a matrix of dissipation times named distimes, @@ -24,8 +27,13 @@ for one of the parents or metabolites, the Eigenvalues are returned. These are equivalent to the rate constants of the DFOP model, but with the advantage that the SFORB model can also be used for metabolites. } +\details{ +Additional DT50 values are calculated from the FOMC DT90 and k1 and k2 from +HS and DFOP, as well as from Eigenvalues b1 and b2 of any SFORB models +} \note{ -The function is used internally by \code{\link{summary.mkinfit}}. +The function is used internally by \link{summary.mkinfit}, +\link{summary.nlme.mmkin} and \link{summary.saem.mmkin}. } \examples{ @@ -42,4 +50,3 @@ The function is used internally by \code{\link{summary.mkinfit}}. \author{ Johannes Ranke } -\keyword{manip} diff --git a/man/saem.Rd b/man/saem.Rd index 39b66448..96b8b55a 100644 --- a/man/saem.Rd +++ b/man/saem.Rd @@ -43,7 +43,8 @@ used.} by the saemix function?} } \value{ -An \link[saemix:SaemixObject-class]{saemix::SaemixObject}. +An S3 object of class 'saem.mmkin', containing the fitted +\link[saemix:SaemixObject-class]{saemix::SaemixObject} as a list component named 'so'. An \link[saemix:SaemixModel-class]{saemix::SaemixModel} object. @@ -51,8 +52,8 @@ An \link[saemix:SaemixData-class]{saemix::SaemixData} object. } \description{ This function uses \code{\link[saemix:saemix]{saemix::saemix()}} as a backend for fitting nonlinear mixed -effects models created from \link{mmkin} row objects using the stochastic approximation -to the expectation maximisation algorithm (SAEM). +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 @@ -95,6 +96,9 @@ f_mmkin_des <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, solution_typ # Using a single core, the following takes about 6 minutes, using 10 cores # it is slower instead of faster f_saem_des <- saem(f_mmkin_des, cores = 1) -compare.saemix(list(f_saemix$so, f_saemix_des$so)) +compare.saemix(list(f_saem$so, f_saem_des$so)) } } +\seealso{ +\link{summary.saem.mmkin} +} diff --git a/man/summary.saem.mmkin.Rd b/man/summary.saem.mmkin.Rd new file mode 100644 index 00000000..0f0d8264 --- /dev/null +++ b/man/summary.saem.mmkin.Rd @@ -0,0 +1,89 @@ +% 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} +\alias{print.summary.saem.mmkin} +\title{Summary method for class "saem.mmkin"} +\usage{ +\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{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.saem.mmkin}} + +\item{digits}{Number of digits to use for printing} +} +\value{ +The summary function returns a list based on the \link[saemix:SaemixObject-class]{saemix::SaemixObject} +obtained in the fit, with at least the following additional components +\item{saemixversion, mkinversion, Rversion}{The saemix, 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{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 SFO kinetics +sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120) +dt50_sfo_in_pop <- 50 +k_in_pop <- log(2) / dt50_sfo_in_pop +set.seed(1234) +k_in <- rlnorm(5, log(k_in_pop), 0.5) +SFO <- mkinmod(parent = mkinsub("SFO")) + +pred_sfo <- function(k) { + mkinpredict(SFO, + c(k_parent = k), + c(parent = 100), + sampling_times) +} + +ds_sfo_mean <- lapply(k_in, pred_sfo) +names(ds_sfo_mean) <- paste("ds", 1:5) + +ds_sfo_syn <- lapply(ds_sfo_mean, function(ds) { + add_err(ds, + sdfunc = function(value) sqrt(1^2 + value^2 * 0.07^2), + n = 1)[[1]] +}) + +\dontrun{ +# Evaluate using mmkin and saem +f_mmkin <- mmkin("SFO", ds_sfo_syn, quiet = TRUE, error_model = "tc", cores = 1) +f_saem <- saem(f_mmkin) +summary(f_saem, data = TRUE) +} + +} +\author{ +Johannes Ranke for the mkin specific parts +saemix authors for the parts inherited from saemix. +} @@ -5,7 +5,7 @@ Testing mkin ✔ | 2 | Export dataset for reading into CAKE ✔ | 14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.9 s] ✔ | 4 | Calculation of FOCUS chi2 error levels [0.4 s] -✔ | 7 | Fitting the SFORB model [3.3 s] +✔ | 7 | Fitting the SFORB model [3.4 s] ✔ | 5 | Analytical solutions for coupled models [3.0 s] ✔ | 5 | Calculation of Akaike weights ✔ | 10 | Confidence intervals and p-values [1.0 s] @@ -25,7 +25,7 @@ Reason: getRversion() < "4.1.0" is TRUE test_nafta.R:53: skip: Test data from Appendix D are correctly evaluated Reason: getRversion() < "4.1.0" is TRUE ──────────────────────────────────────────────────────────────────────────────── -✔ | 9 | Nonlinear mixed-effects models [7.9 s] +✔ | 9 | Nonlinear mixed-effects models [7.8 s] ✔ | 0 1 | Plotting [0.7 s] ──────────────────────────────────────────────────────────────────────────────── test_plot.R:24: skip: Plotting mkinfit and mmkin objects is reproducible @@ -40,7 +40,7 @@ Reason: getRversion() < "4.1.0" is TRUE ✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.5 s] ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 37.7 s +Duration: 37.8 s OK: 146 Failed: 0 diff --git a/tests/testthat/FOCUS_2006_D.csf b/tests/testthat/FOCUS_2006_D.csf index dc369ee1..a0e63657 100644 --- a/tests/testthat/FOCUS_2006_D.csf +++ b/tests/testthat/FOCUS_2006_D.csf @@ -5,7 +5,7 @@ Description: MeasurementUnits: % AR TimeUnits: days Comments: Created using mkin::CAKE_export -Date: 2020-11-07 +Date: 2020-11-08 Optimiser: IRLS [Data] |