From f606838c5310f365eea1c0d6421f5c3636a4dc79 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Tue, 8 Dec 2020 22:08:38 +0100 Subject: mixed.mmkin and test coverage --- R/mixed.mmkin.R | 101 +++++++++++++++++++++++++++++++++++++++++++++++++ R/mkinfit.R | 4 +- R/plot.mixed.mmkin.R | 61 +++++++++++++++++++---------- R/saem.R | 4 +- R/summary.saem.mmkin.R | 75 +++++++++++++++++++----------------- 5 files changed, 188 insertions(+), 57 deletions(-) create mode 100644 R/mixed.mmkin.R (limited to 'R') diff --git a/R/mixed.mmkin.R b/R/mixed.mmkin.R new file mode 100644 index 00000000..6fe5130d --- /dev/null +++ b/R/mixed.mmkin.R @@ -0,0 +1,101 @@ +#' Create a mixed effects model from an mmkin row object +#' +#' @param object An [mmkin] row object +#' @param method The method to be used +#' @param \dots Currently not used +#' @examples +#' sampling_times = c(0, 1, 3, 7, 14, 28, 60, 90, 120) +#' n_biphasic <- 8 +#' err_1 = list(const = 1, prop = 0.07) +#' +#' DFOP_SFO <- mkinmod( +#' parent = mkinsub("DFOP", "m1"), +#' m1 = mkinsub("SFO"), +#' quiet = TRUE) +#' +#' set.seed(123456) +#' log_sd <- 0.3 +#' syn_biphasic_parms <- as.matrix(data.frame( +#' k1 = rlnorm(n_biphasic, log(0.05), log_sd), +#' k2 = rlnorm(n_biphasic, log(0.01), log_sd), +#' g = plogis(rnorm(n_biphasic, 0, log_sd)), +#' f_parent_to_m1 = plogis(rnorm(n_biphasic, 0, log_sd)), +#' k_m1 = rlnorm(n_biphasic, log(0.002), log_sd))) +#' +#' ds_biphasic_mean <- lapply(1:n_biphasic, +#' function(i) { +#' mkinpredict(DFOP_SFO, syn_biphasic_parms[i, ], +#' c(parent = 100, m1 = 0), sampling_times) +#' } +#' ) +#' +#' set.seed(123456L) +#' ds_biphasic <- lapply(ds_biphasic_mean, function(ds) { +#' add_err(ds, +#' sdfunc = function(value) sqrt(err_1$const^2 + value^2 * err_1$prop^2), +#' n = 1, secondary = "m1")[[1]] +#' }) +#' +#' \dontrun{ +#' f_mmkin <- mmkin(list("DFOP-SFO" = DFOP_SFO), ds_biphasic, error_model = "tc", quiet = TRUE) +#' +#' f_mixed <- mixed(f_mmkin) +#' print(f_mixed) +#' plot(f_mixed) +#' } +#' @export +mixed <- function(object, ...) { + UseMethod("mixed") +} + +#' @export +#' @rdname mixed +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) + + if (method[1] == "none") { + ds_list <- lapply(object, + function(x) x$data[c("variable", "time", "observed", "predicted", "residual")]) + + names(ds_list) <- ds_names + res$data <- purrr::map_dfr(ds_list, function(x) x, .id = "ds") + names(res$data)[1:4] <- c("ds", "name", "time", "value") + res$data$name <- as.character(res$data$name) + res$data$ds <- ordered(res$data$ds, levels = unique(res$data$ds)) + standardized <- unlist(lapply(object, residuals, standardized = TRUE)) + res$data$std <- res$data$residual / standardized + res$data$standardized <- standardized + + class(res) <- c("mixed.mmkin") + return(res) + } +} + +#' @export +#' @rdname mixed +#' @param x A mixed.mmkin object to print +#' @param digits Number of digits to use for printing. +print.mixed.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) { + cat("Kinetic model fitted by nonlinear regression to each dataset" ) + 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\n") + + print(x$mmkin) + + cat("\nMean fitted parameters:\n") + print(mean_degparms(x$mmkin)) + + invisible(x) +} diff --git a/R/mkinfit.R b/R/mkinfit.R index a6efc858..c659e446 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -885,8 +885,8 @@ mkinfit <- function(mkinmod, observed, fit$rss <- function(P) cost_function(P, OLS = TRUE, update_data = FALSE) # Log-likelihood with possibility to fix degparms or errparms - fit$ll <- function(P, fixed_degparms = FALSE, fixed_errparms = FALSE) { - - cost_function(P, trans = FALSE, fixed_degparms = fixed_degparms, + fit$ll <- function(P, fixed_degparms = FALSE, fixed_errparms = FALSE, trans = FALSE) { + - cost_function(P, trans = trans, fixed_degparms = fixed_degparms, fixed_errparms = fixed_errparms, OLS = FALSE, update_data = FALSE) } diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R index d8a1c2ac..db29376e 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 [saem.mmkin] or [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 @@ -63,31 +63,43 @@ plot.mixed.mmkin <- function(x, fit_1 <- x$mmkin[[1]] ds_names <- colnames(x$mmkin) + backtransform = TRUE + + if (identical(class(x), "mixed.mmkin")) { + degparms_pop <- mean_degparms(x$mmkin) + + degparms_tmp <- parms(x$mmkin, transformed = TRUE) + degparms_i <- as.data.frame(t(degparms_tmp[setdiff(rownames(degparms_tmp), names(fit_1$errparms)), ])) + residual_type = ifelse(standardized, "standardized", "residual") + residuals <- x$data[[residual_type]] + } + if (inherits(x, "nlme.mmkin")) { - degparms_optim <- coefficients(x) + degparms_i <- coefficients(x) degparms_pop <- nlme::fixef(x) residuals <- residuals(x, type = ifelse(standardized, "pearson", "response")) } if (inherits(x, "saem.mmkin")) { - degparms_optim <- saemix::psi(x$so) - rownames(degparms_optim) <- ds_names - degparms_optim_names <- setdiff(names(fit_1$par), names(fit_1$errparms)) - colnames(degparms_optim) <- degparms_optim_names - residual_type = ifelse(standardized, "residual", "standardized") + 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_optim_names + 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_optim), - matrix(rep(degparms_fixed, nrow(degparms_optim)), + degparms_all <- cbind(as.matrix(degparms_i), + matrix(rep(degparms_fixed, nrow(degparms_i)), ncol = length(degparms_fixed), - nrow = nrow(degparms_optim), byrow = TRUE)) - degparms_all_names <- c(names(degparms_optim), names(degparms_fixed)) + nrow = nrow(degparms_i), byrow = TRUE)) + degparms_all_names <- c(names(degparms_i), names(degparms_fixed)) colnames(degparms_all) <- degparms_all_names degparms_all_pop <- c(degparms_pop, degparms_fixed) @@ -106,10 +118,14 @@ plot.mixed.mmkin <- function(x, pred_ds <- purrr::map_dfr(i, function(ds_i) { odeparms_trans <- degparms_all[ds_i, odeparms_names] names(odeparms_trans) <- odeparms_names # needed if only one odeparm - odeparms <- backtransform_odeparms(odeparms_trans, - x$mkinmod, - transform_rates = fit_1$transform_rates, - transform_fractions = fit_1$transform_fractions) + if (backtransform) { + odeparms <- backtransform_odeparms(odeparms_trans, + x$mkinmod, + transform_rates = fit_1$transform_rates, + transform_fractions = fit_1$transform_fractions) + } else { + odeparms <- odeparms_trans + } odeini <- degparms_all[ds_i, odeini_names] names(odeini) <- gsub("_0", "", odeini_names) @@ -121,10 +137,15 @@ plot.mixed.mmkin <- function(x, }) odeparms_pop_trans <- degparms_all_pop[odeparms_names] - odeparms_pop <- backtransform_odeparms(odeparms_pop_trans, - x$mkinmod, - transform_rates = fit_1$transform_rates, - transform_fractions = fit_1$transform_fractions) + + if (backtransform) { + odeparms_pop <- backtransform_odeparms(odeparms_pop_trans, + x$mkinmod, + transform_rates = fit_1$transform_rates, + transform_fractions = fit_1$transform_fractions) + } else { + odeparms_pop <- odeparms_pop_trans + } odeini_pop <- degparms_all_pop[odeini_names] names(odeini_pop) <- gsub("_0", "", odeini_names) diff --git a/R/saem.R b/R/saem.R index e5634ac7..88b8e172 100644 --- a/R/saem.R +++ b/R/saem.R @@ -92,7 +92,7 @@ utils::globalVariables(c("predicted", "std")) #' plot(f_saem_fomc) #' } #' @export -saem <- function(object, control, ...) UseMethod("saem") +saem <- function(object, ...) UseMethod("saem") #' @rdname saem #' @export @@ -104,9 +104,11 @@ saem.mmkin <- function(object, cores = 1, verbose = FALSE, suppressPlot = TRUE, quiet = FALSE, ...) { + transformations <- match.arg(transformations) m_saemix <- saemix_model(object, cores = cores, verbose = verbose, solution_type = solution_type, transformations = transformations, ...) d_saemix <- saemix_data(object, verbose = verbose) + if (suppressPlot) { # We suppress the log-likelihood curve that saemix currently # produces at the end of the fit by plotting to a file diff --git a/R/summary.saem.mmkin.R b/R/summary.saem.mmkin.R index 97f9f2da..27c2ce6c 100644 --- a/R/summary.saem.mmkin.R +++ b/R/summary.saem.mmkin.R @@ -89,9 +89,42 @@ summary.saem.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes = confint_trans <- as.matrix(conf.int[pnames, c("estimate", "lower", "upper")]) colnames(confint_trans)[1] <- "est." - bp <- backtransform_odeparms(confint_trans[, "est."], object$mkinmod, - object$transform_rates, object$transform_fractions) - bpnames <- names(bp) + 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] @@ -111,34 +144,6 @@ summary.saem.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes = confint_errmod <- as.matrix(conf.int[enames, c("estimate", "lower", "upper")]) colnames(confint_errmod)[1] <- "est." - # Transform boundaries of CI for one parameter at a time, - # with the exception of sets of formation fractions (single fractions are OK). - f_names_skip <- character(0) - for (box in mod_vars) { # Figure out sets of fractions to skip - f_names <- grep(paste("^f", box, sep = "_"), pnames, value = TRUE) - n_paths <- length(f_names) - if (n_paths > 1) f_names_skip <- c(f_names_skip, f_names) - } - - confint_back <- matrix(NA, nrow = length(bp), ncol = 3, - dimnames = list(bpnames, colnames(confint_trans))) - confint_back[, "est."] <- bp - - for (pname in pnames) { - if (!pname %in% f_names_skip) { - par.lower <- confint_trans[pname, "lower"] - par.upper <- confint_trans[pname, "upper"] - names(par.lower) <- names(par.upper) <- pname - bpl <- backtransform_odeparms(par.lower, object$mkinmod, - object$transform_rates, - object$transform_fractions) - bpu <- backtransform_odeparms(par.upper, object$mkinmod, - object$transform_rates, - object$transform_fractions) - confint_back[names(bpl), "lower"] <- bpl - confint_back[names(bpu), "upper"] <- bpu - } - } object$confint_trans <- confint_trans object$confint_ranef <- confint_ranef @@ -213,7 +218,7 @@ print.summary.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3) print(data.frame(AIC = x$AIC, BIC = x$BIC, logLik = x$logLik, row.names = " "), digits = digits) - cat("\nOptimised, transformed parameters with symmetric confidence intervals:\n") + cat("\nOptimised parameters:\n") print(x$confint_trans, digits = digits) if (nrow(x$confint_trans) > 1) { @@ -228,8 +233,10 @@ print.summary.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3) cat("\nVariance model:\n") print(x$confint_errmod, digits = digits) - cat("\nBacktransformed parameters with asymmetric confidence intervals:\n") - print(x$confint_back, digits = digits) + if (x$transformations == "mkin") { + cat("\nBacktransformed parameters:\n") + print(x$confint_back, digits = digits) + } printSFORB <- !is.null(x$SFORB) if(printSFORB){ -- cgit v1.2.1