diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/endpoints.R | 8 | ||||
| -rw-r--r-- | R/plot.mixed.mmkin.R | 23 | ||||
| -rw-r--r-- | R/saem.R | 512 | ||||
| -rw-r--r-- | R/summary.saem.mmkin.R | 268 |
4 files changed, 806 insertions, 5 deletions
diff --git a/R/endpoints.R b/R/endpoints.R index b5872e68..f1f47581 100644 --- a/R/endpoints.R +++ b/R/endpoints.R @@ -10,8 +10,8 @@ #' 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] or [nlme.mmkin] -#' or another object that has list components +#' @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. @@ -20,8 +20,8 @@ #' 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 [summary.mkinfit] -#' and [summary.nlme.mmkin] +#' @note The function is used internally by [summary.mkinfit], +#' [summary.nlme.mmkin] and [summary.saem.mmkin]. #' @author Johannes Ranke #' @examples #' diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R index 4c1f1531..21399496 100644 --- a/R/plot.mixed.mmkin.R +++ b/R/plot.mixed.mmkin.R @@ -2,7 +2,7 @@ utils::globalVariables("ds") #' Plot predictions from a fitted nonlinear mixed model obtained via an mmkin row object #' -#' @param x An object of class [mixed.mmkin], [nlme.mmkin] +#' @param x An object of class [mixed.mmkin], [saem.mmkin] or [nlme.mmkin] #' @param i A numeric index to select datasets for which to plot the individual predictions, #' in case plots get too large #' @inheritParams plot.mkinfit @@ -39,6 +39,15 @@ utils::globalVariables("ds") #' f_nlme <- nlme(f, control = list(pnlsMaxIter = 120, tolerance = 1e-3)) #' plot(f_nlme) #' +#' f_saem <- saem(f, transformations = "saemix") +#' plot(f_saem) +#' +#' # We can overlay the two variants if we generate predictions +#' pred_nlme <- mkinpredict(dfop_sfo, +#' f_nlme$bparms.optim[-1], +#' c(parent = f_nlme$bparms.optim[[1]], A1 = 0), +#' seq(0, 180, by = 0.2)) +#' plot(f_saem, pred_over = list(nlme = pred_nlme)) #' } #' @export plot.mixed.mmkin <- function(x, @@ -82,6 +91,18 @@ plot.mixed.mmkin <- function(x, type = ifelse(standardized, "pearson", "response")) } + if (inherits(x, "saem.mmkin")) { + if (x$transformations == "saemix") backtransform = FALSE + degparms_i <- saemix::psi(x$so) + rownames(degparms_i) <- ds_names + degparms_i_names <- setdiff(x$so@results@name.fixed, names(fit_1$errparms)) + colnames(degparms_i) <- degparms_i_names + residual_type = ifelse(standardized, "standardized", "residual") + residuals <- x$data[[residual_type]] + degparms_pop <- x$so@results@fixed.effects + names(degparms_pop) <- degparms_i_names + } + 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 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) +} diff --git a/R/summary.saem.mmkin.R b/R/summary.saem.mmkin.R new file mode 100644 index 00000000..e92c561c --- /dev/null +++ b/R/summary.saem.mmkin.R @@ -0,0 +1,268 @@ +#' 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{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 +#' saemix authors for the parts inherited from saemix. +#' @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 saem +#' f_mmkin_dfop_sfo <- mmkin(list(dfop_sfo), ds_syn_dfop_sfo, +#' quiet = TRUE, error_model = "tc", cores = 5) +#' f_saem_dfop_sfo <- saem(f_mmkin_dfop_sfo) +#' summary(f_saem_dfop_sfo, 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." + + # In case objects were produced by earlier versions of saem + if (is.null(object$transformations)) object$transformations <- "mkin" + + if (object$transformations == "mkin") { + 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 + } + } + } else { + confint_back <- confint_trans + } + + # 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 <- if (object$err_mod == "const") "a.1" else c("a.1", "b.1") + confint_errmod <- as.matrix(conf.int[enames, c("estimate", "lower", "upper")]) + colnames(confint_errmod)[1] <- "est." + + + 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 # boolean: Should we print the data? + 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) +} |