diff options
Diffstat (limited to 'R')
| -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 |
5 files changed, 341 insertions, 70 deletions
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) +} |