diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/dimethenamid_2018.R | 50 | ||||
| -rw-r--r-- | R/endpoints.R | 8 | ||||
| -rw-r--r-- | R/intervals.R | 179 | ||||
| -rw-r--r-- | R/mean_degparms.R | 62 | ||||
| -rw-r--r-- | R/mixed.mmkin.R | 3 | ||||
| -rw-r--r-- | R/mkinsub.R | 5 | ||||
| -rw-r--r-- | R/nlme.R | 26 | ||||
| -rw-r--r-- | R/nlme.mmkin.R | 11 | ||||
| -rw-r--r-- | R/nlmixr.R | 557 | ||||
| -rw-r--r-- | R/plot.mixed.mmkin.R | 48 | ||||
| -rw-r--r-- | R/saem.R | 542 | ||||
| -rw-r--r-- | R/summary.nlmixr.mmkin.R | 250 | ||||
| -rw-r--r-- | R/summary.saem.mmkin.R | 268 | ||||
| -rw-r--r-- | R/tffm0.R | 48 | ||||
| -rw-r--r-- | R/transform_odeparms.R | 13 |
15 files changed, 2022 insertions, 48 deletions
diff --git a/R/dimethenamid_2018.R b/R/dimethenamid_2018.R index 0d8f681f..2bf5fb77 100644 --- a/R/dimethenamid_2018.R +++ b/R/dimethenamid_2018.R @@ -15,7 +15,55 @@ #' @source Rapporteur Member State Germany, Co-Rapporteur Member State Bulgaria (2018) #' Renewal Assessment Report Dimethenamid-P Volume 3 - B.8 Environmental fate and behaviour #' Rev. 2 - November 2017 -#' \href{https://open.efsa.europa.eu}{https://open.efsa.europa.eu/study-inventory/EFSA-Q-2014-00716} +#' \url{https://open.efsa.europa.eu/study-inventory/EFSA-Q-2014-00716} #' @examples #' print(dimethenamid_2018) +#' dmta_ds <- lapply(1:7, function(i) { +#' ds_i <- dimethenamid_2018$ds[[i]]$data +#' ds_i[ds_i$name == "DMTAP", "name"] <- "DMTA" +#' ds_i$time <- ds_i$time * dimethenamid_2018$f_time_norm[i] +#' ds_i +#' }) +#' names(dmta_ds) <- sapply(dimethenamid_2018$ds, function(ds) ds$title) +#' dmta_ds[["Elliot"]] <- rbind(dmta_ds[["Elliot 1"]], dmta_ds[["Elliot 2"]]) +#' dmta_ds[["Elliot 1"]] <- NULL +#' dmta_ds[["Elliot 2"]] <- NULL +#' \dontrun{ +#' dfop_sfo3_plus <- mkinmod( +#' DMTA = mkinsub("DFOP", c("M23", "M27", "M31")), +#' M23 = mkinsub("SFO"), +#' M27 = mkinsub("SFO"), +#' M31 = mkinsub("SFO", "M27", sink = FALSE), +#' quiet = TRUE +#' ) +#' f_dmta_mkin_tc <- mmkin( +#' list("DFOP-SFO3+" = dfop_sfo3_plus), +#' dmta_ds, quiet = TRUE, error_model = "tc") +#' nlmixr_model(f_dmta_mkin_tc) +#' # The focei fit takes about four minutes on my system +#' system.time( +#' f_dmta_nlmixr_focei <- nlmixr(f_dmta_mkin_tc, est = "focei", +#' control = nlmixr::foceiControl(print = 500)) +#' ) +#' summary(f_dmta_nlmixr_focei) +#' plot(f_dmta_nlmixr_focei) +#' # Using saemix takes about 18 minutes +#' system.time( +#' f_dmta_saemix <- saem(f_dmta_mkin_tc, test_log_parms = TRUE) +#' ) +#' +#' # nlmixr with est = "saem" is pretty fast with default iteration numbers, most +#' # of the time (about 2.5 minutes) is spent for calculating the log likelihood at the end +#' # The likelihood calculated for the nlmixr fit is much lower than that found by saemix +#' # Also, the trace plot and the plot of the individual predictions is not +#' # convincing for the parent. It seems we are fitting an overparameterised +#' # model, so the result we get strongly depends on starting parameters and control settings. +#' system.time( +#' f_dmta_nlmixr_saem <- nlmixr(f_dmta_mkin_tc, est = "saem", +#' control = nlmixr::saemControl(print = 500, logLik = TRUE, nmc = 9)) +#' ) +#' traceplot(f_dmta_nlmixr_saem$nm) +#' summary(f_dmta_nlmixr_saem) +#' plot(f_dmta_nlmixr_saem) +#' } "dimethenamid_2018" diff --git a/R/endpoints.R b/R/endpoints.R index b5872e68..6bf52f07 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], [saem.mmkin] or +#' [nlmixr.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/intervals.R b/R/intervals.R new file mode 100644 index 00000000..e2d342f0 --- /dev/null +++ b/R/intervals.R @@ -0,0 +1,179 @@ +#' @importFrom nlme intervals +#' @export +nlme::intervals + +#' Confidence intervals for parameters in saem.mmkin objects +#' +#' @param object The fitted saem.mmkin object +#' @param level The confidence level. Must be the default of 0.95 as this is what +#' is available in the saemix object +#' @param backtransform In case the model was fitted with mkin transformations, +#' should we backtransform the parameters where a one to one correlation +#' between transformed and backtransformed parameters exists? +#' @return An object with 'intervals.saem.mmkin' and 'intervals.lme' in the +#' class attribute +#' @export +intervals.saem.mmkin <- function(object, level = 0.95, backtransform = TRUE, ...) +{ + if (!identical(level, 0.95)) { + stop("Confidence intervals are only available for a level of 95%") + } + + mod_vars <- names(object$mkinmod$diffs) + + pnames <- names(object$mean_dp_start) + + # Confidence intervals are available in the SaemixObject, so + # we just need to extract them and put them into a list modelled + # after the result of nlme::intervals.lme + + conf.int <- object$so@results@conf.int + rownames(conf.int) <- conf.int$name + colnames(conf.int)[2] <- "est." + confint_trans <- as.matrix(conf.int[pnames, c("lower", "est.", "upper")]) + + # Fixed effects + # In case objects were produced by earlier versions of saem + if (is.null(object$transformations)) object$transformations <- "mkin" + + if (object$transformations == "mkin" & backtransform) { + 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 + } + } + confint_ret <- confint_back + } else { + confint_ret <- confint_trans + } + attr(confint_ret, "label") <- "Fixed effects:" + + # Random effects + ranef_ret <- as.matrix(conf.int[paste0("SD.", pnames), c("lower", "est.", "upper")]) + rownames(ranef_ret) <- paste0(gsub("SD\\.", "sd(", rownames(ranef_ret)), ")") + attr(ranef_ret, "label") <- "Random effects:" + + + # Error model + enames <- if (object$err_mod == "const") "a.1" else c("a.1", "b.1") + err_ret <- as.matrix(conf.int[enames, c("lower", "est.", "upper")]) + + res <- list( + fixed = confint_ret, + random = ranef_ret, + errmod = err_ret + ) + class(res) <- c("intervals.saemix.mmkin", "intervals.lme") + attr(res, "level") <- level + return(res) +} + +#' Confidence intervals for parameters in nlmixr.mmkin objects +#' +#' @param object The fitted saem.mmkin object +#' @param level The confidence level. +#' @param backtransform Should we backtransform the parameters where a one to +#' one correlation between transformed and backtransformed parameters exists? +#' @importFrom nlme intervals +#' @return An object with 'intervals.saem.mmkin' and 'intervals.lme' in the +#' class attribute +#' @export +intervals.nlmixr.mmkin <- function(object, level = 0.95, backtransform = TRUE, ...) +{ + + # Fixed effects + mod_vars <- names(object$mkinmod$diffs) + + conf.int <- confint(object$nm) + dpnames <- setdiff(rownames(conf.int), names(object$mean_ep_start)) + ndp <- length(dpnames) + + confint_trans <- as.matrix(conf.int[dpnames, c(3, 1, 4)]) + colnames(confint_trans) <- c("lower", "est.", "upper") + + if (backtransform) { + 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 = "_"), dpnames, 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 dpnames) { + 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 + } + } + confint_ret <- confint_back + } else { + confint_ret <- confint_trans + } + attr(confint_ret, "label") <- "Fixed effects:" + + # Random effects + ranef_ret <- as.matrix(data.frame(lower = NA, + "est." = sqrt(diag(object$nm$omega)), upper = NA)) + rownames(ranef_ret) <- paste0(gsub("eta\\.", "sd(", rownames(ranef_ret)), ")") + attr(ranef_ret, "label") <- "Random effects:" + + # Error model + enames <- names(object$nm$sigma) + err_ret <- as.matrix(conf.int[enames, c(3, 1, 4)]) + colnames(err_ret) <- c("lower", "est.", "upper") + + res <- list( + fixed = confint_ret, + random = ranef_ret, + errmod = err_ret + ) + class(res) <- c("intervals.nlmixr.mmkin", "intervals.lme") + attr(res, "level") <- level + return(res) +} diff --git a/R/mean_degparms.R b/R/mean_degparms.R new file mode 100644 index 00000000..ec20c068 --- /dev/null +++ b/R/mean_degparms.R @@ -0,0 +1,62 @@ +#' Calculate mean degradation parameters for an mmkin row object +#' +#' @return If random is FALSE (default), a named vector containing mean values +#' of the fitted degradation model parameters. If random is TRUE, a list with +#' fixed and random effects, in the format required by the start argument of +#' nlme for the case of a single grouping variable ds. +#' @param object An mmkin row object containing several fits of the same model to different datasets +#' @param random Should a list with fixed and random effects be returned? +#' @param test_log_parms If TRUE, log parameters are only considered in +#' the mean calculations if their untransformed counterparts (most likely +#' rate constants) pass the t-test for significant difference from zero. +#' @param conf.level Possibility to adjust the required confidence level +#' for parameter that are tested if requested by 'test_log_parms'. +#' @export +mean_degparms <- function(object, random = FALSE, test_log_parms = FALSE, conf.level = 0.6) +{ + if (nrow(object) > 1) stop("Only row objects allowed") + parm_mat_trans <- sapply(object, parms, transformed = TRUE) + + if (test_log_parms) { + parm_mat_dim <- dim(parm_mat_trans) + parm_mat_dimnames <- dimnames(parm_mat_trans) + + log_parm_trans_names <- grep("^log_", rownames(parm_mat_trans), value = TRUE) + log_parm_names <- gsub("^log_", "", log_parm_trans_names) + + t_test_back_OK <- matrix( + sapply(object, function(o) { + suppressWarnings(summary(o)$bpar[log_parm_names, "Pr(>t)"] < (1 - conf.level)) + }), nrow = length(log_parm_names)) + rownames(t_test_back_OK) <- log_parm_trans_names + + parm_mat_trans_OK <- parm_mat_trans + for (trans_parm in log_parm_trans_names) { + parm_mat_trans_OK[trans_parm, ] <- ifelse(t_test_back_OK[trans_parm, ], + parm_mat_trans[trans_parm, ], NA) + } + } else { + parm_mat_trans_OK <- parm_mat_trans + } + + mean_degparm_names <- setdiff(rownames(parm_mat_trans), names(object[[1]]$errparms)) + degparm_mat_trans <- parm_mat_trans[mean_degparm_names, , drop = FALSE] + degparm_mat_trans_OK <- parm_mat_trans_OK[mean_degparm_names, , drop = FALSE] + + fixed <- apply(degparm_mat_trans_OK, 1, mean, na.rm = TRUE) + if (random) { + random <- t(apply(degparm_mat_trans[mean_degparm_names, , drop = FALSE], 2, function(column) column - fixed)) + # If we only have one parameter, apply returns a vector so we get a single row + if (nrow(degparm_mat_trans) == 1) random <- t(random) + rownames(random) <- levels(nlme_data(object)$ds) + + # For nlmixr we can specify starting values for standard deviations eta, and + # we ignore uncertain parameters if test_log_parms is FALSE + eta <- apply(degparm_mat_trans_OK, 1, stats::sd, na.rm = TRUE) + + return(list(fixed = fixed, random = list(ds = random), eta = eta)) + } else { + return(fixed) + } +} + diff --git a/R/mixed.mmkin.R b/R/mixed.mmkin.R index 7aa5edd5..682a9a34 100644 --- a/R/mixed.mmkin.R +++ b/R/mixed.mmkin.R @@ -3,6 +3,8 @@ #' @param object An [mmkin] row object #' @param method The method to be used #' @param \dots Currently not used +#' @return An object of class 'mixed.mmkin' which has the observed data in a +#' single dataframe which is convenient for plotting #' @examples #' sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120) #' n_biphasic <- 8 @@ -54,7 +56,6 @@ mixed.mmkin <- function(object, method = c("none"), ...) { if (nrow(object) > 1) stop("Only row objects allowed") method <- match.arg(method) - if (method == "default") method = c("naive", "nlme") ds_names <- colnames(object) res <- list(mmkin = object, mkinmod = object[[1]]$mkinmod) diff --git a/R/mkinsub.R b/R/mkinsub.R index 886f712c..93af3f16 100644 --- a/R/mkinsub.R +++ b/R/mkinsub.R @@ -1,8 +1,3 @@ -#' Function to set up a kinetic submodel for one state variable -#' -#' This is a convenience function to set up the lists used as arguments for -#' \code{\link{mkinmod}}. -#' #' @rdname mkinmod #' @param submodel Character vector of length one to specify the submodel type. #' See \code{\link{mkinmod}} for the list of allowed submodel names. @@ -36,7 +36,7 @@ #' nlme_f <- nlme_function(f) #' # These assignments are necessary for these objects to be #' # visible to nlme and augPred when evaluation is done by -#' # pkgdown to generated the html docs. +#' # pkgdown to generate the html docs. #' assign("nlme_f", nlme_f, globalenv()) #' assign("grouped_data", grouped_data, globalenv()) #' @@ -125,30 +125,6 @@ nlme_function <- function(object) { } #' @rdname nlme -#' @return If random is FALSE (default), a named vector containing mean values -#' of the fitted degradation model parameters. If random is TRUE, a list with -#' fixed and random effects, in the format required by the start argument of -#' nlme for the case of a single grouping variable ds. -#' @param random Should a list with fixed and random effects be returned? -#' @export -mean_degparms <- function(object, random = FALSE) { - if (nrow(object) > 1) stop("Only row objects allowed") - parm_mat_trans <- sapply(object, parms, transformed = TRUE) - mean_degparm_names <- setdiff(rownames(parm_mat_trans), names(object[[1]]$errparms)) - degparm_mat_trans <- parm_mat_trans[mean_degparm_names, , drop = FALSE] - fixed <- apply(degparm_mat_trans, 1, mean) - if (random) { - random <- t(apply(degparm_mat_trans[mean_degparm_names, , drop = FALSE], 2, function(column) column - fixed)) - # If we only have one parameter, apply returns a vector so we get a single row - if (nrow(degparm_mat_trans) == 1) random <- t(random) - rownames(random) <- levels(nlme_data(object)$ds) - return(list(fixed = fixed, random = list(ds = random))) - } else { - return(fixed) - } -} - -#' @rdname nlme #' @importFrom purrr map_dfr #' @return A \code{\link{groupedData}} object #' @export diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R index ff1f2fff..7049a9a1 100644 --- a/R/nlme.mmkin.R +++ b/R/nlme.mmkin.R @@ -24,7 +24,7 @@ get_deg_func <- function() { #' This functions sets up a nonlinear mixed effects model for an mmkin row #' object. An mmkin row object is essentially a list of mkinfit objects that #' have been obtained by fitting the same model to a list of datasets. -#' +#' #' Note that the convergence of the nlme algorithms depends on the quality #' of the data. In degradation kinetics, we often only have few datasets #' (e.g. data for few soils) and complicated degradation models, which may @@ -34,10 +34,9 @@ get_deg_func <- function() { #' @param data Ignored, data are taken from the mmkin model #' @param fixed Ignored, all degradation parameters fitted in the #' mmkin model are used as fixed parameters -#' @param random If not specified, correlated random effects are set up -#' for all optimised degradation model parameters using the log-Cholesky -#' parameterization [nlme::pdLogChol] that is also the default of -#' the generic [nlme] method. +#' @param random If not specified, no correlations between random effects are +#' set up for the optimised degradation model parameters. This is +#' achieved by using the [nlme::pdDiag] method. #' @param groups See the documentation of nlme #' @param start If not specified, mean values of the fitted degradation #' parameters taken from the mmkin object are used @@ -135,7 +134,7 @@ nlme.mmkin <- function(model, data = "auto", function(el) eval(parse(text = paste(el, 1, sep = "~")))), random = pdDiag(fixed), groups, - start = mean_degparms(model, random = TRUE), + start = mean_degparms(model, random = TRUE, test_log_parms = TRUE), correlation = NULL, weights = NULL, subset, method = c("ML", "REML"), na.action = na.fail, naPattern, diff --git a/R/nlmixr.R b/R/nlmixr.R new file mode 100644 index 00000000..fd12f555 --- /dev/null +++ b/R/nlmixr.R @@ -0,0 +1,557 @@ +utils::globalVariables(c("predicted", "std", "ID", "TIME", "CMT", "DV", "IPRED", "IRES", "IWRES")) + +#' @export +nlmixr::nlmixr + +#' Fit nonlinear mixed models using nlmixr +#' +#' This function uses [nlmixr::nlmixr()] as a backend for fitting nonlinear mixed +#' effects models created from [mmkin] row objects using the Stochastic Approximation +#' Expectation Maximisation algorithm (SAEM) or First Order Conditional +#' Estimation with Interaction (FOCEI). +#' +#' An mmkin row object is essentially a list of mkinfit objects that have been +#' obtained by fitting the same model to a list of datasets using [mkinfit]. +#' +#' @importFrom nlmixr nlmixr tableControl +#' @importFrom dplyr %>% +#' @param object An [mmkin] row object containing several fits of the same +#' [mkinmod] model to different datasets +#' @param data Not used, the data are extracted from the mmkin row object +#' @param est Estimation method passed to [nlmixr::nlmixr] +#' @param degparms_start Parameter values given as a named numeric vector will +#' be used to override the starting values obtained from the 'mmkin' object. +#' @param eta_start Standard deviations on the transformed scale given as a +#' named numeric vector will be used to override the starting values obtained +#' from the 'mmkin' object. +#' @param test_log_parms If TRUE, an attempt is made to use more robust starting +#' values for population parameters fitted as log parameters in mkin (like +#' rate constants) by only considering rate constants that pass the t-test +#' when calculating mean degradation parameters using [mean_degparms]. +#' @param conf.level Possibility to adjust the required confidence level +#' for parameter that are tested if requested by 'test_log_parms'. +#' @param data Not used, as the data are extracted from the mmkin row object +#' @param table Passed to [nlmixr::nlmixr] +#' @param error_model Optional argument to override the error model which is +#' being set based on the error model used in the mmkin row object. +#' @param control Passed to [nlmixr::nlmixr] +#' @param \dots Passed to [nlmixr_model] +#' @param save Passed to [nlmixr::nlmixr] +#' @param envir Passed to [nlmixr::nlmixr] +#' @return An S3 object of class 'nlmixr.mmkin', containing the fitted +#' [nlmixr::nlmixr] object as a list component named 'nm'. The +#' object also inherits from 'mixed.mmkin'. +#' @seealso [summary.nlmixr.mmkin] [plot.mixed.mmkin] +#' @examples +#' \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 <- mmkin(c("SFO", "FOMC", "DFOP", "HS"), ds, quiet = TRUE, cores = 1) +#' f_mmkin_parent_tc <- mmkin(c("SFO", "FOMC", "DFOP"), ds, error_model = "tc", +#' cores = 1, quiet = TRUE) +#' +#' library(nlmixr) +#' f_nlmixr_sfo_saem <- nlmixr(f_mmkin_parent["SFO", ], est = "saem", +#' control = saemControl(print = 0)) +#' f_nlmixr_sfo_focei <- nlmixr(f_mmkin_parent["SFO", ], est = "focei", +#' control = foceiControl(print = 0)) +#' +#' f_nlmixr_fomc_saem <- nlmixr(f_mmkin_parent["FOMC", ], est = "saem", +#' control = saemControl(print = 0)) +#' f_nlmixr_fomc_focei <- nlmixr(f_mmkin_parent["FOMC", ], est = "focei", +#' control = foceiControl(print = 0)) +#' +#' f_nlmixr_dfop_saem <- nlmixr(f_mmkin_parent["DFOP", ], est = "saem", +#' control = saemControl(print = 0)) +#' f_nlmixr_dfop_focei <- nlmixr(f_mmkin_parent["DFOP", ], est = "focei", +#' control = foceiControl(print = 0)) +#' +#' f_nlmixr_hs_saem <- nlmixr(f_mmkin_parent["HS", ], est = "saem", +#' control = saemControl(print = 0)) +#' f_nlmixr_hs_focei <- nlmixr(f_mmkin_parent["HS", ], est = "focei", +#' control = foceiControl(print = 0)) +#' +#' f_nlmixr_fomc_saem_tc <- nlmixr(f_mmkin_parent_tc["FOMC", ], est = "saem", +#' control = saemControl(print = 0)) +#' f_nlmixr_fomc_focei_tc <- nlmixr(f_mmkin_parent_tc["FOMC", ], est = "focei", +#' control = foceiControl(print = 0)) +#' +#' AIC( +#' f_nlmixr_sfo_saem$nm, f_nlmixr_sfo_focei$nm, +#' f_nlmixr_fomc_saem$nm, f_nlmixr_fomc_focei$nm, +#' f_nlmixr_dfop_saem$nm, f_nlmixr_dfop_focei$nm, +#' f_nlmixr_hs_saem$nm, f_nlmixr_hs_focei$nm, +#' f_nlmixr_fomc_saem_tc$nm, f_nlmixr_fomc_focei_tc$nm) +#' +#' AIC(nlme(f_mmkin_parent["FOMC", ])) +#' AIC(nlme(f_mmkin_parent["HS", ])) +#' +#' # The FOCEI fit of FOMC with constant variance or the tc error model is best +#' plot(f_nlmixr_fomc_saem_tc) +#' +#' sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"), +#' A1 = mkinsub("SFO"), quiet = TRUE) +#' fomc_sfo <- mkinmod(parent = mkinsub("FOMC", "A1"), +#' A1 = mkinsub("SFO"), quiet = TRUE) +#' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "A1"), +#' A1 = mkinsub("SFO"), quiet = TRUE) +#' +#' f_mmkin_const <- mmkin(list( +#' "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo), +#' ds, quiet = TRUE, error_model = "const") +#' f_mmkin_obs <- mmkin(list( +#' "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo), +#' ds, quiet = TRUE, error_model = "obs") +#' f_mmkin_tc <- mmkin(list( +#' "SFO-SFO" = sfo_sfo, "FOMC-SFO" = fomc_sfo, "DFOP-SFO" = dfop_sfo), +#' ds, quiet = TRUE, error_model = "tc") +#' +#' nlmixr_model(f_mmkin_const["SFO-SFO", ]) +#' +#' # A single constant variance is currently only possible with est = 'focei' in nlmixr +#' f_nlmixr_sfo_sfo_focei_const <- nlmixr(f_mmkin_const["SFO-SFO", ], est = "focei") +#' f_nlmixr_fomc_sfo_focei_const <- nlmixr(f_mmkin_const["FOMC-SFO", ], est = "focei") +#' f_nlmixr_dfop_sfo_focei_const <- nlmixr(f_mmkin_const["DFOP-SFO", ], est = "focei") +#' +#' # Variance by variable is supported by 'saem' and 'focei' +#' f_nlmixr_fomc_sfo_saem_obs <- nlmixr(f_mmkin_obs["FOMC-SFO", ], est = "saem") +#' f_nlmixr_fomc_sfo_focei_obs <- nlmixr(f_mmkin_obs["FOMC-SFO", ], est = "focei") +#' f_nlmixr_dfop_sfo_saem_obs <- nlmixr(f_mmkin_obs["DFOP-SFO", ], est = "saem") +#' f_nlmixr_dfop_sfo_focei_obs <- nlmixr(f_mmkin_obs["DFOP-SFO", ], est = "focei") +#' +#' # Identical two-component error for all variables is only possible with +#' # est = 'focei' in nlmixr +#' f_nlmixr_fomc_sfo_focei_tc <- nlmixr(f_mmkin_tc["FOMC-SFO", ], est = "focei") +#' f_nlmixr_dfop_sfo_focei_tc <- nlmixr(f_mmkin_tc["DFOP-SFO", ], est = "focei") +#' +#' # Two-component error by variable is possible with both estimation methods +#' # Variance by variable is supported by 'saem' and 'focei' +#' f_nlmixr_fomc_sfo_saem_obs_tc <- nlmixr(f_mmkin_tc["FOMC-SFO", ], est = "saem", +#' error_model = "obs_tc") +#' f_nlmixr_fomc_sfo_focei_obs_tc <- nlmixr(f_mmkin_tc["FOMC-SFO", ], est = "focei", +#' error_model = "obs_tc") +#' f_nlmixr_dfop_sfo_saem_obs_tc <- nlmixr(f_mmkin_tc["DFOP-SFO", ], est = "saem", +#' error_model = "obs_tc") +#' f_nlmixr_dfop_sfo_focei_obs_tc <- nlmixr(f_mmkin_tc["DFOP-SFO", ], est = "focei", +#' error_model = "obs_tc") +#' +#' AIC( +#' f_nlmixr_sfo_sfo_focei_const$nm, +#' f_nlmixr_fomc_sfo_focei_const$nm, +#' f_nlmixr_dfop_sfo_focei_const$nm, +#' f_nlmixr_fomc_sfo_saem_obs$nm, +#' f_nlmixr_fomc_sfo_focei_obs$nm, +#' f_nlmixr_dfop_sfo_saem_obs$nm, +#' f_nlmixr_dfop_sfo_focei_obs$nm, +#' f_nlmixr_fomc_sfo_focei_tc$nm, +#' f_nlmixr_dfop_sfo_focei_tc$nm, +#' f_nlmixr_fomc_sfo_saem_obs_tc$nm, +#' f_nlmixr_fomc_sfo_focei_obs_tc$nm, +#' f_nlmixr_dfop_sfo_saem_obs_tc$nm, +#' f_nlmixr_dfop_sfo_focei_obs_tc$nm +#' ) +#' # Currently, FOMC-SFO with two-component error by variable fitted by focei gives the +#' # lowest AIC +#' plot(f_nlmixr_fomc_sfo_focei_obs_tc) +#' summary(f_nlmixr_fomc_sfo_focei_obs_tc) +#' } +#' @export +nlmixr.mmkin <- function(object, data = NULL, + est = NULL, control = list(), + table = tableControl(), + error_model = object[[1]]$err_mod, + test_log_parms = TRUE, + conf.level = 0.6, + degparms_start = "auto", + eta_start = "auto", + ..., + save = NULL, + envir = parent.frame() +) +{ + m_nlmixr <- nlmixr_model(object, est = est, + error_model = error_model, add_attributes = TRUE, + test_log_parms = test_log_parms, conf.level = conf.level, + degparms_start = degparms_start, eta_start = eta_start + ) + d_nlmixr <- nlmixr_data(object) + + mean_dp_start <- attr(m_nlmixr, "mean_dp_start") + mean_ep_start <- attr(m_nlmixr, "mean_ep_start") + + attributes(m_nlmixr) <- NULL + + fit_time <- system.time({ + f_nlmixr <- nlmixr(m_nlmixr, d_nlmixr, est = est, control = control) + }) + + if (is.null(f_nlmixr$CMT)) { + nlmixr_df <- as.data.frame(f_nlmixr[c("ID", "TIME", "DV", "IPRED", "IRES", "IWRES")]) + nlmixr_df$CMT <- as.character(object[[1]]$data$variable[1]) + } else { + nlmixr_df <- as.data.frame(f_nlmixr[c("ID", "TIME", "DV", "CMT", "IPRED", "IRES", "IWRES")]) + } + + return_data <- nlmixr_df %>% + dplyr::transmute(ds = ID, name = CMT, time = TIME, value = DV, + predicted = IPRED, residual = IRES, + std = IRES/IWRES, standardized = IWRES) %>% + dplyr::arrange(ds, name, time) + + bparms_optim <- backtransform_odeparms(f_nlmixr$theta, + object[[1]]$mkinmod, + object[[1]]$transform_rates, + object[[1]]$transform_fractions) + + result <- list( + mkinmod = object[[1]]$mkinmod, + mmkin = object, + transform_rates = object[[1]]$transform_rates, + transform_fractions = object[[1]]$transform_fractions, + nm = f_nlmixr, + est = est, + time = fit_time, + mean_dp_start = mean_dp_start, + mean_ep_start = mean_ep_start, + bparms.optim = bparms_optim, + bparms.fixed = object[[1]]$bparms.fixed, + data = return_data, + err_mod = error_model, + date.fit = date(), + nlmixrversion = as.character(utils::packageVersion("nlmixr")), + mkinversion = as.character(utils::packageVersion("mkin")), + Rversion = paste(R.version$major, R.version$minor, sep=".") + ) + + class(result) <- c("nlmixr.mmkin", "mixed.mmkin") + return(result) +} + +#' @export +#' @rdname nlmixr.mmkin +#' @param x An nlmixr.mmkin object to print +#' @param digits Number of digits to use for printing +print.nlmixr.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) { + cat("Kinetic nonlinear mixed-effects model fit by", x$est, "using nlmixr") + cat("\nStructural model:\n") + diffs <- x$mmkin[[1]]$mkinmod$diffs + nice_diffs <- gsub("^(d.*) =", "\\1/dt =", diffs) + writeLines(strwrap(nice_diffs, exdent = 11)) + cat("\nData:\n") + cat(nrow(x$data), "observations of", + length(unique(x$data$name)), "variable(s) grouped in", + length(unique(x$data$ds)), "datasets\n") + + cat("\nLikelihood:\n") + print(data.frame( + AIC = AIC(x$nm), + BIC = BIC(x$nm), + logLik = logLik(x$nm), + row.names = " "), digits = digits) + + cat("\nFitted parameters:\n") + print(x$nm$parFixed, digits = digits) + + invisible(x) +} + +#' @rdname nlmixr.mmkin +#' @param add_attributes Should the starting values used for degradation model +#' parameters and their distribution and for the error model parameters +#' be returned as attributes? +#' @return An function defining a model suitable for fitting with [nlmixr::nlmixr]. +#' @export +nlmixr_model <- function(object, + est = c("saem", "focei"), + degparms_start = "auto", + eta_start = "auto", + test_log_parms = TRUE, conf.level = 0.6, + error_model = object[[1]]$err_mod, add_attributes = FALSE) +{ + if (nrow(object) > 1) stop("Only row objects allowed") + est = match.arg(est) + + mkin_model <- object[[1]]$mkinmod + obs_vars <- names(mkin_model$spec) + + if (error_model == object[[1]]$err_mod) { + + if (length(object[[1]]$mkinmod$spec) > 1 & est == "saem") { + if (error_model == "const") { + message( + "Constant variance for more than one variable is not supported for est = 'saem'\n", + "Changing the error model to 'obs' (variance by observed variable)") + error_model <- "obs" + } + if (error_model =="tc") { + message( + "With est = 'saem', a different error model is required for each observed variable", + "Changing the error model to 'obs_tc' (Two-component error for each observed variable)") + error_model <- "obs_tc" + } + } + } + + degparms_mmkin <- mean_degparms(object, + test_log_parms = test_log_parms, + conf.level = conf.level, random = TRUE) + + degparms_optim <- degparms_mmkin$fixed + + degparms_optim_ilr_names <- grep("^f_.*_ilr", names(degparms_optim), value = TRUE) + obs_vars_ilr <- unique(gsub("f_(.*)_ilr.*$", "\\1", degparms_optim_ilr_names)) + degparms_optim_noilr <- degparms_optim[setdiff(names(degparms_optim), + degparms_optim_ilr_names)] + + degparms_optim_back <- backtransform_odeparms(degparms_optim, + object[[1]]$mkinmod, + object[[1]]$transform_rates, + object[[1]]$transform_fractions) + + if (degparms_start[1] == "auto") { + degparms_start <- degparms_optim_noilr + for (obs_var_ilr in obs_vars_ilr) { + ff_names <- grep(paste0("^f_", obs_var_ilr, "_"), + names(degparms_optim_back), value = TRUE) + f_tffm0 <- tffm0(degparms_optim_back[ff_names]) + f_tffm0_qlogis <- qlogis(f_tffm0) + names(f_tffm0_qlogis) <- paste0("f_", obs_var_ilr, + "_tffm0_", 1:length(f_tffm0), "_qlogis") + degparms_start <- c(degparms_start, f_tffm0_qlogis) + } + } + + if (eta_start[1] == "auto") { + eta_start <- degparms_mmkin$eta[setdiff(names(degparms_optim), + degparms_optim_ilr_names)] + for (obs_var_ilr in obs_vars_ilr) { + ff_n <- length(grep(paste0("^f_", obs_var_ilr, "_"), + names(degparms_optim_back), value = TRUE)) + eta_start_ff <- rep(0.3, ff_n) + names(eta_start_ff) <- paste0("f_", obs_var_ilr, + "_tffm0_", 1:ff_n, "_qlogis") + eta_start <- c(eta_start, eta_start_ff) + } + } + + + degparms_fixed <- object[[1]]$bparms.fixed + + odeini_optim_parm_names <- grep('_0$', names(degparms_optim), value = TRUE) + odeini_fixed_parm_names <- grep('_0$', names(degparms_fixed), value = TRUE) + + odeparms_fixed_names <- setdiff(names(degparms_fixed), odeini_fixed_parm_names) + odeparms_fixed <- degparms_fixed[odeparms_fixed_names] + + odeini_fixed <- degparms_fixed[odeini_fixed_parm_names] + names(odeini_fixed) <- gsub('_0$', '', odeini_fixed_parm_names) + + # Definition of the model function + f <- function(){} + + ini_block <- "ini({" + + # Initial values for all degradation parameters + for (parm_name in names(degparms_start)) { + # As initials for state variables are not transformed, + # we need to modify the name here as we want to + # use the original name in the model block + ini_block <- paste0( + ini_block, + parm_name, " = ", + as.character(signif(degparms_start[parm_name], 2)), + "\n", + "eta.", parm_name, " ~ ", + as.character(signif(eta_start[parm_name], 2)), + "\n" + ) + } + + # Error model parameters + error_model_mkin <- object[[1]]$err_mod + + errparm_names_mkin <- names(object[[1]]$errparms) + errparms_mkin <- sapply(errparm_names_mkin, function(parm_name) { + mean(sapply(object, function(x) x$errparms[parm_name])) + }) + + sigma_tc_mkin <- errparms_ini <- errparms_mkin[1] + + mean(unlist(sapply(object, function(x) x$data$observed)), na.rm = TRUE) * + errparms_mkin[2] + + if (error_model == "const") { + if (error_model_mkin == "tc") { + errparms_ini <- sigma_tc_mkin + } else { + errparms_ini <- mean(errparms_mkin) + } + names(errparms_ini) <- "sigma" + } + + if (error_model == "obs") { + errparms_ini <- switch(error_model_mkin, + const = rep(errparms_mkin["sigma"], length(obs_vars)), + obs = errparms_mkin, + tc = sigma_tc_mkin) + names(errparms_ini) <- paste0("sigma_", obs_vars) + } + + if (error_model == "tc") { + if (error_model_mkin != "tc") { + stop("Not supported") + } else { + errparms_ini <- errparms_mkin + } + } + + if (error_model == "obs_tc") { + if (error_model_mkin != "tc") { + stop("Not supported") + } else { + errparms_ini <- rep(errparms_mkin, length(obs_vars)) + names(errparms_ini) <- paste0( + rep(names(errparms_mkin), length(obs_vars)), + "_", + rep(obs_vars, each = 2)) + } + } + + for (parm_name in names(errparms_ini)) { + ini_block <- paste0( + ini_block, + parm_name, " = ", + as.character(signif(errparms_ini[parm_name], 2)), + "\n" + ) + } + + ini_block <- paste0(ini_block, "})") + + body(f)[2] <- parse(text = ini_block) + + model_block <- "model({" + + # Population initial values for the ODE state variables + for (parm_name in odeini_optim_parm_names) { + model_block <- paste0( + model_block, + parm_name, "_model = ", + parm_name, " + eta.", parm_name, "\n", + gsub("(.*)_0", "\\1(0)", parm_name), " = ", parm_name, "_model\n") + } + + # Population initial values for log rate constants + for (parm_name in grep("^log_", names(degparms_start), value = TRUE)) { + model_block <- paste0( + model_block, + gsub("^log_", "", parm_name), " = ", + "exp(", parm_name, " + eta.", parm_name, ")\n") + } + + # Population initial values for logit transformed parameters + for (parm_name in grep("_qlogis$", names(degparms_start), value = TRUE)) { + parm_name_new <- names( + backtransform_odeparms(degparms_start[parm_name], + object[[1]]$mkinmod, + object[[1]]$transform_rates, + object[[1]]$transform_fractions)) + model_block <- paste0( + model_block, + parm_name_new, " = ", + "expit(", parm_name, " + eta.", parm_name, ")\n") + } + + # Calculate formation fractions from tffm0 transformed values + for (obs_var_ilr in obs_vars_ilr) { + ff_names <- grep(paste0("^f_", obs_var_ilr, "_"), + names(degparms_optim_back), value = TRUE) + pattern <- paste0("^f_", obs_var_ilr, "_to_(.*)$") + target_vars <- gsub(pattern, "\\1", + grep(paste0("^f_", obs_var_ilr, "_to_"), names(degparms_optim_back), value = TRUE)) + for (i in 1:length(target_vars)) { + ff_name <- ff_names[i] + ff_line <- paste0(ff_name, " = f_", obs_var_ilr, "_tffm0_", i) + if (i > 1) { + for (j in (i - 1):1) { + ff_line <- paste0(ff_line, " * (1 - f_", obs_var_ilr, "_tffm0_", j , ")") + } + } + model_block <- paste0( + model_block, + ff_line, + "\n" + ) + } + } + + # Differential equations + model_block <- paste0( + model_block, + paste( + gsub("d_(.*) =", "d/dt(\\1) =", mkin_model$diffs), + collapse = "\n"), + "\n" + ) + + # Error model + if (error_model == "const") { + model_block <- paste0(model_block, + paste(paste0(obs_vars, " ~ add(sigma)"), collapse = "\n")) + } + if (error_model == "obs") { + model_block <- paste0(model_block, + paste(paste0(obs_vars, " ~ add(sigma_", obs_vars, ")"), collapse = "\n"), + "\n") + } + if (error_model == "tc") { + model_block <- paste0(model_block, + paste(paste0(obs_vars, " ~ add(sigma_low) + prop(rsd_high)"), collapse = "\n"), + "\n") + } + if (error_model == "obs_tc") { + model_block <- paste0(model_block, + paste( + paste0(obs_vars, " ~ add(sigma_low_", obs_vars, ") + ", + "prop(rsd_high_", obs_vars, ")"), collapse = "\n"), + "\n") + } + + model_block <- paste0(model_block, "})") + + body(f)[3] <- parse(text = model_block) + + if (add_attributes) { + attr(f, "mean_dp_start") <- degparms_optim + attr(f, "eta_start") <- degparms_mmkin$eta + attr(f, "mean_ep_start") <- errparms_ini + } + + return(f) +} + +#' @rdname nlmixr.mmkin +#' @return An dataframe suitable for use with [nlmixr::nlmixr] +#' @export +nlmixr_data <- function(object, ...) { + if (nrow(object) > 1) stop("Only row objects allowed") + d <- lapply(object, function(x) x$data) + compartment_map <- 1:length(object[[1]]$mkinmod$spec) + names(compartment_map) <- names(object[[1]]$mkinmod$spec) + ds_names <- colnames(object) + + ds_list <- lapply(object, function(x) x$data[c("time", "variable", "observed")]) + names(ds_list) <- ds_names + ds_nlmixr <- purrr::map_dfr(ds_list, function(x) x, .id = "ds") + ds_nlmixr$variable <- as.character(ds_nlmixr$variable) + ds_nlmixr_renamed <- data.frame( + ID = ds_nlmixr$ds, + TIME = ds_nlmixr$time, + AMT = 0, EVID = 0, + CMT = ds_nlmixr$variable, + DV = ds_nlmixr$observed, + stringsAsFactors = FALSE) + + return(ds_nlmixr_renamed) +} diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R index 4c1f1531..1ac62b07 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 @@ -10,6 +10,10 @@ utils::globalVariables("ds") #' `resplot = "time"`. #' @param pred_over Named list of alternative predictions as obtained #' from [mkinpredict] with a compatible [mkinmod]. +#' @param test_log_parms Passed to [mean_degparms] in the case of an +#' [mixed.mmkin] object +#' @param conf.level Passed to [mean_degparms] in the case of an +#' [mixed.mmkin] object #' @param rel.height.legend The relative height of the legend shown on top #' @param rel.height.bottom The relative height of the bottom plot row #' @param ymax Vector of maximum y axis values @@ -36,9 +40,23 @@ utils::globalVariables("ds") #' #' # For this fit we need to increase pnlsMaxiter, and we increase the #' # tolerance in order to speed up the fit for this example evaluation +#' # It still takes 20 seconds to run #' f_nlme <- nlme(f, control = list(pnlsMaxIter = 120, tolerance = 1e-3)) #' plot(f_nlme) #' +#' f_saem <- saem(f, transformations = "saemix") +#' plot(f_saem) +#' +#' f_obs <- mmkin(list("DFOP-SFO" = dfop_sfo), ds, quiet = TRUE, error_model = "obs") +#' f_nlmix <- nlmix(f_obs) +#' plot(f_nlmix) +#' +#' # We can overlay the two variants if we generate predictions +#' pred_nlme <- mkinpredict(dfop_sfo, +#' f_nlme$bparms.optim[-1], +#' 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, @@ -49,6 +67,8 @@ plot.mixed.mmkin <- function(x, xlim = range(x$data$time), resplot = c("predicted", "time"), pred_over = NULL, + test_log_parms = FALSE, + conf.level = 0.6, ymax = "auto", maxabs = "auto", ncol.legend = ifelse(length(i) <= 3, length(i) + 1, ifelse(length(i) <= 8, 3, 4)), nrow.legend = ceiling((length(i) + 1) / ncol.legend), @@ -67,7 +87,7 @@ plot.mixed.mmkin <- function(x, backtransform = TRUE if (identical(class(x), "mixed.mmkin")) { - degparms_pop <- mean_degparms(x$mmkin) + degparms_pop <- mean_degparms(x$mmkin, test_log_parms = test_log_parms, conf.level = conf.level) degparms_tmp <- parms(x$mmkin, transformed = TRUE) degparms_i <- as.data.frame(t(degparms_tmp[setdiff(rownames(degparms_tmp), names(fit_1$errparms)), ])) @@ -82,6 +102,30 @@ 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 + } + + if (inherits(x, "nlmixr.mmkin")) { + eta_i <- random.effects(x$nm)[-1] + names(eta_i) <- gsub("^eta.", "", names(eta_i)) + degparms_i <- eta_i + degparms_pop <- x$nm$theta + for (parm_name in names(degparms_i)) { + degparms_i[parm_name] <- eta_i[parm_name] + degparms_pop[parm_name] + } + residual_type = ifelse(standardized, "standardized", "residual") + residuals <- x$data[[residual_type]] + } + degparms_fixed <- fit_1$fixed$value names(degparms_fixed) <- rownames(fit_1$fixed) degparms_all <- cbind(as.matrix(degparms_i), diff --git a/R/saem.R b/R/saem.R new file mode 100644 index 00000000..2c20f788 --- /dev/null +++ b/R/saem.R @@ -0,0 +1,542 @@ +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]. +#' +#' @importFrom utils packageVersion +#' @param object An [mmkin] row object containing several fits of the same +#' [mkinmod] model to different datasets +#' @param verbose Should we print information about created objects of +#' 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 test_log_parms If TRUE, an attempt is made to use more robust starting +#' values for population parameters fitted as log parameters in mkin (like +#' rate constants) by only considering rate constants that pass the t-test +#' when calculating mean degradation parameters using [mean_degparms]. +#' @param conf.level Possibility to adjust the required confidence level +#' for parameter that are tested if requested by 'test_log_parms'. +#' @param solution_type Possibility to specify the solution type in case the +#' automatic choice is not desired +#' @param fail_with_errors Should a failure to compute standard errors +#' from the inverse of the Fisher Information Matrix be a failure? +#' @param quiet Should we suppress the messages saemix prints at the beginning +#' and the end of the optimisation process? +#' @param nbiter.saemix Convenience option to increase the number of +#' iterations +#' @param control Passed to [saemix::saemix]. +#' @param \dots Further parameters passed to [saemix::saemixModel]. +#' @return An S3 object of class 'saem.mmkin', containing the fitted +#' [saemix::SaemixObject] as a list component named 'so'. The +#' 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(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(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(), + test_log_parms = TRUE, + conf.level = 0.6, + solution_type = "auto", + nbiter.saemix = c(300, 100), + control = list(displayProgress = FALSE, print = FALSE, + nbiter.saemix = nbiter.saemix, + save = FALSE, save.graphs = FALSE), + fail_with_errors = TRUE, + verbose = FALSE, quiet = FALSE, ...) +{ + transformations <- match.arg(transformations) + m_saemix <- saemix_model(object, verbose = verbose, + degparms_start = degparms_start, + test_log_parms = test_log_parms, conf.level = conf.level, + solution_type = solution_type, + transformations = transformations, ...) + d_saemix <- saemix_data(object, verbose = verbose) + + fit_time <- system.time({ + utils::capture.output(f_saemix <- saemix::saemix(m_saemix, d_saemix, control), split = !quiet) + FIM_failed <- NULL + if (any(is.na(f_saemix@results@se.fixed))) FIM_failed <- c(FIM_failed, "fixed effects") + if (any(is.na(c(f_saemix@results@se.omega, f_saemix@results@se.respar)))) { + FIM_failed <- c(FIM_failed, "random effects and residual error parameters") + } + if (!is.null(FIM_failed) & fail_with_errors) { + stop("Could not invert FIM for ", paste(FIM_failed, collapse = " and ")) + } + }) + + transparms_optim <- f_saemix@results@fixed.effects + 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 = value - predicted, + 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(), test_log_parms = FALSE, verbose = FALSE, ...) +{ + if (packageVersion("saemix") < "3.0") { + stop("To use the interface to saemix, you need to install a version >= 3.0\n") + } + + if (nrow(object) > 1) stop("Only row objects allowed") + + mkin_model <- object[[1]]$mkinmod + + degparms_optim <- mean_degparms(object, test_log_parms = test_log_parms) + 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.nlmixr.mmkin.R b/R/summary.nlmixr.mmkin.R new file mode 100644 index 00000000..a023f319 --- /dev/null +++ b/R/summary.nlmixr.mmkin.R @@ -0,0 +1,250 @@ +#' Summary method for class "nlmixr.mmkin" +#' +#' Lists model equations, initial parameter values, optimised parameters +#' for fixed effects (population), random effects (deviations from the +#' population mean) and residual error model, as well as the resulting +#' endpoints such as formation fractions and DT50 values. Optionally +#' (default is FALSE), the data are listed in full. +#' +#' @importFrom stats confint sd +#' @param object an object of class [nlmixr.mmkin] +#' @param x an object of class [summary.nlmixr.mmkin] +#' @param data logical, indicating whether the full data should be included in +#' the summary. +#' @param verbose Should the summary be verbose? +#' @param distimes logical, indicating whether DT50 and DT90 values should be +#' included. +#' @param digits Number of digits to use for printing +#' @param \dots optional arguments passed to methods like \code{print}. +#' @return The summary function returns a list obtained in the fit, with at +#' least the following additional components +#' \item{nlmixrversion, mkinversion, Rversion}{The nlmixr, mkin and R versions used} +#' \item{date.fit, date.summary}{The dates where the fit and the summary were +#' produced} +#' \item{diffs}{The differential equations used in the degradation model} +#' \item{use_of_ff}{Was maximum or minimum use made of formation fractions} +#' \item{data}{The data} +#' \item{confint_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 +#' nlmixr authors for the parts inherited from nlmixr. +#' @examples +#' # Generate five datasets following DFOP-SFO kinetics +#' sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120) +#' dfop_sfo <- mkinmod(parent = mkinsub("DFOP", "m1"), +#' m1 = mkinsub("SFO"), quiet = TRUE) +#' set.seed(1234) +#' k1_in <- rlnorm(5, log(0.1), 0.3) +#' k2_in <- rlnorm(5, log(0.02), 0.3) +#' g_in <- plogis(rnorm(5, qlogis(0.5), 0.3)) +#' f_parent_to_m1_in <- plogis(rnorm(5, qlogis(0.3), 0.3)) +#' k_m1_in <- rlnorm(5, log(0.02), 0.3) +#' +#' pred_dfop_sfo <- function(k1, k2, g, f_parent_to_m1, k_m1) { +#' mkinpredict(dfop_sfo, +#' c(k1 = k1, k2 = k2, g = g, f_parent_to_m1 = f_parent_to_m1, k_m1 = k_m1), +#' c(parent = 100, m1 = 0), +#' sampling_times) +#' } +#' +#' ds_mean_dfop_sfo <- lapply(1:5, function(i) { +#' mkinpredict(dfop_sfo, +#' c(k1 = k1_in[i], k2 = k2_in[i], g = g_in[i], +#' f_parent_to_m1 = f_parent_to_m1_in[i], k_m1 = k_m1_in[i]), +#' c(parent = 100, m1 = 0), +#' sampling_times) +#' }) +#' names(ds_mean_dfop_sfo) <- paste("ds", 1:5) +#' +#' ds_syn_dfop_sfo <- lapply(ds_mean_dfop_sfo, function(ds) { +#' add_err(ds, +#' sdfunc = function(value) sqrt(1^2 + value^2 * 0.07^2), +#' n = 1)[[1]] +#' }) +#' +#' \dontrun{ +#' # Evaluate using mmkin and nlmixr +#' f_mmkin_dfop_sfo <- mmkin(list(dfop_sfo), ds_syn_dfop_sfo, +#' quiet = TRUE, error_model = "tc", cores = 5) +#' f_saemix_dfop_sfo <- mkin::saem(f_mmkin_dfop_sfo) +#' f_nlme_dfop_sfo <- mkin::nlme(f_mmkin_dfop_sfo) +#' f_nlmixr_dfop_sfo_saem <- nlmixr(f_mmkin_dfop_sfo, est = "saem") +#' # The following takes a very long time but gives +#' f_nlmixr_dfop_sfo_focei <- nlmixr(f_mmkin_dfop_sfo, est = "focei") +#' AIC(f_nlmixr_dfop_sfo_saem$nm, f_nlmixr_dfop_sfo_focei$nm) +#' summary(f_nlmixr_dfop_sfo_sfo, data = TRUE) +#' } +#' +#' @export +summary.nlmixr.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes = TRUE, ...) { + + mod_vars <- names(object$mkinmod$diffs) + + conf.int <- confint(object$nm) + dpnames <- setdiff(rownames(conf.int), names(object$mean_ep_start)) + ndp <- length(dpnames) + + confint_trans <- as.matrix(conf.int[dpnames, c(1, 3, 4)]) + colnames(confint_trans) <- c("est.", "lower", "upper") + + bp <- backtransform_odeparms(confint_trans[, "est."], object$mkinmod, + object$transform_rates, object$transform_fractions) + bpnames <- names(bp) + + # Transform boundaries of CI for one parameter at a time, + # with the exception of sets of formation fractions (single fractions are OK). + f_names_skip <- character(0) + for (box in mod_vars) { # Figure out sets of fractions to skip + f_names <- grep(paste("^f", box, sep = "_"), dpnames, 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 dpnames) { + if (!pname %in% f_names_skip) { + par.lower <- confint_trans[pname, "lower"] + par.upper <- confint_trans[pname, "upper"] + names(par.lower) <- names(par.upper) <- pname + bpl <- backtransform_odeparms(par.lower, object$mkinmod, + object$transform_rates, + object$transform_fractions) + bpu <- backtransform_odeparms(par.upper, object$mkinmod, + object$transform_rates, + object$transform_fractions) + confint_back[names(bpl), "lower"] <- bpl + confint_back[names(bpu), "upper"] <- bpu + } + } + + # Correlation of fixed effects (inspired by summary.nlme) + varFix <- vcov(object$nm) + stdFix <- sqrt(diag(varFix)) + object$corFixed <- array( + t(varFix/stdFix)/stdFix, + dim(varFix), + list(dpnames, dpnames)) + + object$confint_trans <- confint_trans + object$confint_back <- confint_back + + object$date.summary = date() + object$use_of_ff = object$mkinmod$use_of_ff + + object$diffs <- object$mkinmod$diffs + object$print_data <- data # boolean: Should we print the data? + + names(object$data)[4] <- "observed" # rename value to observed + + object$verbose <- verbose + + object$fixed <- object$mmkin_orig[[1]]$fixed + object$AIC = AIC(object$nm) + object$BIC = BIC(object$nm) + object$logLik = logLik(object$nm) + + 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.nlmixr.mmkin") + return(object) +} + +#' @rdname summary.nlmixr.mmkin +#' @export +print.summary.nlmixr.mmkin <- function(x, digits = max(3, getOption("digits") - 3), verbose = x$verbose, ...) { + cat("nlmixr version used for fitting: ", x$nlmixrversion, "\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("\nDegradation model predictions using RxODE\n") + + cat("\nFitted in", x$time[["elapsed"]], "s\n") + + cat("\nVariance model: ") + cat(switch(x$err_mod, + const = "Constant variance", + obs = "Variance unique to each observed variable", + tc = "Two-component variance function", + obs_tc = "Two-component variance unique to each observed variable"), "\n") + + cat("\nMean of starting values for individual parameters:\n") + print(x$mean_dp_start, digits = digits) + + cat("\nMean of starting values for error model parameters:\n") + print(x$mean_ep_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 calculated by", nlmixr::getOfvType(x$nm), " \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 (omega):\n") + print(x$nm$omega, digits = digits) + + cat("\nVariance model:\n") + print(x$nm$sigma, digits = digits) + + cat("\nBacktransformed parameters:\n") + print(x$confint_back, digits = digits) + + printSFORB <- !is.null(x$SFORB) + if(printSFORB){ + cat("\nEstimated Eigenvalues of SFORB model(s):\n") + print(x$SFORB, digits = digits,...) + } + + printff <- !is.null(x$ff) + if(printff){ + cat("\nResulting formation fractions:\n") + print(data.frame(ff = x$ff), digits = digits,...) + } + + printdistimes <- !is.null(x$distimes) + if(printdistimes){ + cat("\nEstimated disappearance times:\n") + print(x$distimes, digits = digits,...) + } + + if (x$print_data){ + cat("\nData:\n") + print(format(x$data, digits = digits, ...), row.names = FALSE) + } + + invisible(x) +} diff --git a/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) +} diff --git a/R/tffm0.R b/R/tffm0.R new file mode 100644 index 00000000..bb5f4cf5 --- /dev/null +++ b/R/tffm0.R @@ -0,0 +1,48 @@ +#' Transform formation fractions as in the first published mkin version +#' +#' The transformed fractions can be restricted between 0 and 1 in model +#' optimisations. Therefore this transformation was used originally in mkin. It +#' was later replaced by the [ilr] transformation because the ilr transformed +#' fractions can assumed to follow normal distribution. As the ilr +#' transformation is not available in [RxODE] and can therefore not be used in +#' the nlmixr modelling language, this transformation is currently used for +#' translating mkin models with formation fractions to more than one target +#' compartment for fitting with nlmixr in [nlmixr_model]. However, +#' this implementation cannot be used there, as it is not accessible +#' from RxODE. +#' +#' @param ff Vector of untransformed formation fractions. The sum +#' must be smaller or equal to one +#' @param ff_trans Vector of transformed formation fractions that can be +#' restricted to the interval from 0 to 1 +#' @return A vector of the transformed formation fractions +#' @export +#' @examples +#' ff_example <- c( +#' 0.10983681, 0.09035905, 0.08399383 +#' ) +#' ff_example_trans <- tffm0(ff_example) +#' invtffm0(ff_example_trans) +tffm0 <- function(ff) { + n <- length(ff) + res <- numeric(n) + f_remaining <- 1 + for (i in 1:n) { + res[i] <- ff[i]/f_remaining + f_remaining <- f_remaining - ff[i] + } + return(res) +} +#' @rdname tffm0 +#' @export +#' @return A vector of backtransformed formation fractions for natural use in degradation models +invtffm0 <- function(ff_trans) { + n <- length(ff_trans) + res <- numeric(n) + f_remaining <- 1 + for (i in 1:n) { + res[i] <- ff_trans[i] * f_remaining + f_remaining <- f_remaining - res[i] + } + return(res) +} diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R index 4fe4e5c2..174e7c2d 100644 --- a/R/transform_odeparms.R +++ b/R/transform_odeparms.R @@ -229,13 +229,18 @@ backtransform_odeparms <- function(transparms, mkinmod, if (length(trans_f) > 0) { if(transform_fractions) { if (any(grepl("qlogis", names(trans_f)))) { - parms[f_names] <- plogis(trans_f) + f_tmp <- plogis(trans_f) + if (any(grepl("_tffm0_.*_qlogis$", names(f_tmp)))) { + parms[f_names] <- invtffm0(f_tmp) + } else { + parms[f_names] <- f_tmp + } } else { - f <- invilr(trans_f) + f_tmp <- invilr(trans_f) if (spec[[box]]$sink) { - parms[f_names] <- f[1:length(f)-1] + parms[f_names] <- f_tmp[1:length(f_tmp)-1] } else { - parms[f_names] <- f + parms[f_names] <- f_tmp } } } else { |