diff options
Diffstat (limited to 'R/saem.R')
| -rw-r--r-- | R/saem.R | 402 |
1 files changed, 318 insertions, 84 deletions
@@ -15,35 +15,46 @@ utils::globalVariables(c("predicted", "std")) #' #' @importFrom utils packageVersion #' @param object An [mmkin] row object containing several fits of the same -#' [mkinmod] model to different datasets +#' [mkinmod] model to different datasets #' @param verbose Should we print information about created objects of -#' type [saemix::SaemixModel] and [saemix::SaemixData]? +#' 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. +#' in mkin. If this argument is set to 'saemix', parameter transformations +#' are done in 'saemix' for the supported cases, i.e. (as of version 1.1.2) +#' SFO, FOMC, DFOP and HS without fixing `parent_0`, and SFO or DFOP with +#' one SFO metabolite. #' @param degparms_start Parameter values given as a named numeric vector will -#' be used to override the starting values obtained from the 'mmkin' object. +#' 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]. +#' 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'. +#' 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 +#' automatic choice is not desired +#' @param no_random_effect Character vector of degradation parameters for +#' which there should be no variability over the groups. Only used +#' if the covariance model is not explicitly specified. +#' @param covariance.model Will be passed to [saemix::SaemixModel()]. Per +#' default, uncorrelated random effects are specified for all degradation +#' parameters. +#' @param covariates A data frame with covariate data for use in +#' 'covariate_models', with dataset names as row names. +#' @param covariate_models A list containing linear model formulas with one explanatory +#' variable, i.e. of the type 'parameter ~ covariate'. Covariates must be available +#' in the 'covariates' data frame. #' @param fail_with_errors Should a failure to compute standard errors -#' from the inverse of the Fisher Information Matrix be a failure? +#' 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? +#' and the end of the optimisation process? #' @param nbiter.saemix Convenience option to increase the number of -#' iterations +#' 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'. +#' [saemix::SaemixObject] as a list component named 'so'. The +#' object also inherits from 'mixed.mmkin'. #' @seealso [summary.saem.mmkin] [plot.mixed.mmkin] #' @examples #' \dontrun{ @@ -58,7 +69,13 @@ utils::globalVariables(c("predicted", "std")) #' f_saem_sfo <- saem(f_mmkin_parent["SFO", ]) #' f_saem_fomc <- saem(f_mmkin_parent["FOMC", ]) #' f_saem_dfop <- saem(f_mmkin_parent["DFOP", ]) +#' anova(f_saem_sfo, f_saem_fomc, f_saem_dfop) +#' anova(f_saem_sfo, f_saem_dfop, test = TRUE) +#' illparms(f_saem_dfop) +#' f_saem_dfop_red <- update(f_saem_dfop, no_random_effect = "g_qlogis") +#' anova(f_saem_dfop, f_saem_dfop_red, test = TRUE) #' +#' anova(f_saem_sfo, f_saem_fomc, f_saem_dfop) #' # The returned saem.mmkin object contains an SaemixObject, therefore we can use #' # functions from saemix #' library(saemix) @@ -70,7 +87,7 @@ utils::globalVariables(c("predicted", "std")) #' #' 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) +#' anova(f_saem_fomc, f_saem_fomc_tc, test = TRUE) #' #' sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"), #' A1 = mkinsub("SFO")) @@ -118,6 +135,10 @@ saem.mmkin <- function(object, test_log_parms = TRUE, conf.level = 0.6, solution_type = "auto", + covariance.model = "auto", + covariates = NULL, + covariate_models = NULL, + no_random_effect = NULL, nbiter.saemix = c(300, 100), control = list(displayProgress = FALSE, print = FALSE, nbiter.saemix = nbiter.saemix, @@ -125,56 +146,68 @@ saem.mmkin <- function(object, fail_with_errors = TRUE, verbose = FALSE, quiet = FALSE, ...) { + call <- match.call() 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) + transformations = transformations, + covariance.model = covariance.model, + covariates = covariates, + covariate_models = covariate_models, + no_random_effect = no_random_effect, + ...) + d_saemix <- saemix_data(object, covariates = covariates, verbose = verbose) + fit_failed <- FALSE + FIM_failed <- NULL fit_time <- system.time({ - utils::capture.output(f_saemix <- saemix::saemix(m_saemix, d_saemix, control), split = !quiet) - FIM_failed <- NULL + utils::capture.output(f_saemix <- try(saemix::saemix(m_saemix, d_saemix, control)), split = !quiet) + if (inherits(f_saemix, "try-error")) fit_failed <- TRUE + }) + + return_data <- nlme_data(object) + + if (!fit_failed) { 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") + FIM_failed <- c(FIM_failed, "random effects and error model 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 + 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 - } + 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) - saemix_data_ds <- f_saemix@data@data$ds - mkin_ds_order <- as.character(unique(return_data$ds)) - saemix_ds_order <- unique(saemix_data_ds) - - psi <- saemix::psi(f_saemix) - rownames(psi) <- saemix_ds_order - return_data$predicted <- f_saemix@model@model( - psi = psi[mkin_ds_order, ], - 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) + saemix_data_ds <- f_saemix@data@data$ds + mkin_ds_order <- as.character(unique(return_data$ds)) + saemix_ds_order <- unique(saemix_data_ds) + + psi <- saemix::psi(f_saemix) + rownames(psi) <- saemix_ds_order + return_data$predicted <- f_saemix@model@model( + psi = psi[mkin_ds_order, ], + 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, @@ -183,14 +216,14 @@ saem.mmkin <- function(object, transformations = transformations, transform_rates = object[[1]]$transform_rates, transform_fractions = object[[1]]$transform_fractions, + sm = m_saemix, so = f_saemix, + call = call, time = fit_time, + FIM_failed = FIM_failed, mean_dp_start = attr(m_saemix, "mean_dp_start"), - bparms.optim = bparms_optim, bparms.fixed = object[[1]]$bparms.fixed, data = return_data, - mkin_ds_order = mkin_ds_order, - saemix_ds_order = saemix_ds_order, err_mod = object[[1]]$err_mod, date.fit = date(), saemixversion = as.character(utils::packageVersion("saemix")), @@ -198,6 +231,12 @@ saem.mmkin <- function(object, Rversion = paste(R.version$major, R.version$minor, sep=".") ) + if (!fit_failed) { + result$mkin_ds_order <- mkin_ds_order + result$saemix_ds_order <- saemix_ds_order + result$bparms.optim <- bparms_optim + } + class(result) <- c("saem.mmkin", "mixed.mmkin") return(result) } @@ -217,18 +256,20 @@ print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) { 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"]] - conf.int.var <- grepl("^Var\\.", rownames(conf.int)) - print(conf.int[!conf.int.var, ], digits = digits) + if (inherits(x$so, "try-error")) { + cat("\nFit did not terminate successfully\n") + } else { + 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 <- parms(x, ci = TRUE) + print(conf.int, digits = digits) + } invisible(x) } @@ -236,14 +277,23 @@ print.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3), ...) { #' @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, conf.level = 0.6, verbose = FALSE, ...) +saemix_model <- function(object, solution_type = "auto", + transformations = c("mkin", "saemix"), degparms_start = numeric(), + covariance.model = "auto", no_random_effect = NULL, + covariates = NULL, covariate_models = NULL, + test_log_parms = FALSE, conf.level = 0.6, verbose = FALSE, ...) { 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) + na_degparms <- names(which(is.na(degparms_optim))) + if (length(na_degparms) > 0) { + message("Did not find valid starting values for ", paste(na_degparms, collapse = ", "), "\n", + "Now trying with test_log_parms = FALSE") + degparms_optim <- mean_degparms(object, test_log_parms = FALSE) + } if (transformations == "saemix") { degparms_optim <- backtransform_odeparms(degparms_optim, object[[1]]$mkinmod, @@ -252,7 +302,8 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk } degparms_fixed <- object[[1]]$bparms.fixed - # Transformations are done in the degradation function + # Transformations are done in the degradation function by default + # (transformations = "mkin") transform.par = rep(0, length(degparms_optim)) odeini_optim_parm_names <- grep('_0$', names(degparms_optim), value = TRUE) @@ -275,7 +326,10 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk # Parent only if (length(mkin_model$spec) == 1) { parent_type <- mkin_model$spec[[1]]$type - if (length(odeini_fixed) == 1) { + if (length(odeini_fixed) == 1 && !grepl("_bound$", names(odeini_fixed))) { + if (transformations == "saemix") { + stop("saemix transformations are not supported for parent fits with fixed initial parent value") + } if (parent_type == "SFO") { stop("saemix needs at least two parameters to work on.") } @@ -303,6 +357,9 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk } } } else { + if (length(odeini_fixed) == 2) { + stop("SFORB with fixed initial parent value is not supported") + } if (parent_type == "SFO") { if (transformations == "mkin") { model_function <- function(psi, id, xidep) { @@ -316,8 +373,15 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk } } 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 (transformations == "mkin") { + model_function <- function(psi, id, xidep) { + psi[id, 1] / (xidep[, "time"]/exp(psi[id, 3]) + 1)^exp(psi[id, 2]) + } + } else { + model_function <- function(psi, id, xidep) { + psi[id, 1] / (xidep[, "time"]/psi[id, 3] + 1)^psi[id, 2] + } + transform.par = c(0, 1, 1) } } if (parent_type == "DFOP") { @@ -338,14 +402,57 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk transform.par = c(0, 1, 1, 3) } } + if (parent_type == "SFORB") { + if (transformations == "mkin") { + model_function <- function(psi, id, xidep) { + k_12 <- exp(psi[id, 3]) + k_21 <- exp(psi[id, 4]) + k_1output <- exp(psi[id, 2]) + t <- xidep[, "time"] + + sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21) + b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp + b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp + + psi[id, 1] * (((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) + + ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t)) + } + } else { + model_function <- function(psi, id, xidep) { + k_12 <- psi[id, 3] + k_21 <- psi[id, 4] + k_1output <- psi[id, 2] + t <- xidep[, "time"] + + sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21) + b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp + b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp + + psi[id, 1] * (((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) + + ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t)) + } + transform.par = c(0, 1, 1, 1) + } + } 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))) + if (transformations == "mkin") { + 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))) + } + } else { + model_function <- function(psi, id, xidep) { + tb <- psi[id, 4] + t <- xidep[, "time"] + psi[id, 1] * ifelse(t <= tb, + exp(- psi[id, 2] * t), + exp(- psi[id, 2] * tb) * exp(- psi[id, 3] * (t - tb))) + } + transform.par = c(0, 1, 1, 1) } } } @@ -518,8 +625,52 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk degparms_psi0 <- degparms_optim degparms_psi0[names(degparms_start)] <- degparms_start - psi0_matrix <- matrix(degparms_psi0, nrow = 1) - colnames(psi0_matrix) <- names(degparms_psi0) + psi0_matrix <- matrix(degparms_psi0, nrow = 1, + dimnames = list("(Intercept)", names(degparms_psi0))) + + if (covariance.model[1] == "auto") { + covariance_diagonal <- rep(1, length(degparms_optim)) + if (!is.null(no_random_effect)) { + degparms_no_random <- which(names(degparms_psi0) %in% no_random_effect) + covariance_diagonal[degparms_no_random] <- 0 + } + covariance.model = diag(covariance_diagonal) + } + + if (is.null(covariate_models)) { + covariate.model <- matrix(nrow = 0, ncol = 0) # default in saemixModel() + } else { + degparms_dependent <- sapply(covariate_models, function(x) as.character(x[[2]])) + covariates_in_models = unique(unlist(lapply( + covariate_models, function(x) + colnames(attr(terms(x), "factors")) + ))) + covariates_not_available <- setdiff(covariates_in_models, names(covariates)) + if (length(covariates_not_available) > 0) { + stop("Covariate(s) ", paste(covariates_not_available, collapse = ", "), + " used in the covariate models not available in the covariate data") + } + psi0_matrix <- rbind(psi0_matrix, + matrix(0, nrow = length(covariates), ncol = ncol(psi0_matrix), + dimnames = list(names(covariates), colnames(psi0_matrix)))) + covariate.model <- matrix(0, nrow = length(covariates), + ncol = ncol(psi0_matrix), + dimnames = list( + covariates = names(covariates), + degparms = colnames(psi0_matrix))) + if (transformations == "saemix") { + stop("Covariate models with saemix transformations currently not supported") + } + parms_trans <- as.data.frame(t(sapply(object, parms, transformed = TRUE))) + for (covariate_model in covariate_models) { + covariate_name <- as.character(covariate_model[[2]]) + model_data <- cbind(parms_trans, covariates) + ini_model <- lm(covariate_model, data = model_data) + ini_coef <- coef(ini_model) + psi0_matrix[names(ini_coef), covariate_name] <- ini_coef + covariate.model[names(ini_coef)[-1], covariate_name] <- as.numeric(as.logical(ini_coef[-1])) + } + } res <- saemix::saemixModel(model_function, psi0 = psi0_matrix, @@ -527,6 +678,8 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk transform.par = transform.par, error.model = error.model, verbose = verbose, + covariance.model = covariance.model, + covariate.model = covariate.model, ... ) attr(res, "mean_dp_start") <- degparms_optim @@ -534,26 +687,107 @@ saemix_model <- function(object, solution_type = "auto", transformations = c("mk } #' @rdname saem +#' @importFrom rlang !!! #' @return An [saemix::SaemixData] object. #' @export -saemix_data <- function(object, verbose = FALSE, ...) { +saemix_data <- function(object, covariates = NULL, 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_all <- vctrs::vec_rbind(!!!ds_list, .names_to = "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) + if (!is.null(covariates)) { + name.covariates <- names(covariates) + covariates$ds <- rownames(covariates) + ds_saemix <- merge(ds_saemix, covariates, sort = FALSE) + } else { + name.covariates <- character(0) + } res <- saemix::saemixData(ds_saemix, name.group = "ds", name.predictors = c("time", "name"), name.response = "value", + name.covariates = name.covariates, verbose = verbose, ...) return(res) } + +#' logLik method for saem.mmkin objects +#' +#' @param object The fitted [saem.mmkin] object +#' @param \dots Passed to [saemix::logLik.SaemixObject] +#' @param method Passed to [saemix::logLik.SaemixObject] +#' @export +logLik.saem.mmkin <- function(object, ..., method = c("is", "lin", "gq")) { + method <- match.arg(method) + return(logLik(object$so, method = method)) +} + +#' @export +update.saem.mmkin <- function(object, ..., evaluate = TRUE) { + call <- object$call + # For some reason we get saem.mmkin in the call when using mhmkin + # so we need to fix this so we do not have to export saem.mmkin in + # addition to the S3 method + call[[1]] <- saem + + # We also need to provide the mmkin object in the call, so it + # will also be found when called by testthat or pkgdown + call[[2]] <- object$mmkin + + update_arguments <- match.call(expand.dots = FALSE)$... + + if (length(update_arguments) > 0) { + update_arguments_in_call <- !is.na(match(names(update_arguments), names(call))) + } + + for (a in names(update_arguments)[update_arguments_in_call]) { + call[[a]] <- update_arguments[[a]] + } + + update_arguments_not_in_call <- !update_arguments_in_call + if(any(update_arguments_not_in_call)) { + call <- c(as.list(call), update_arguments[update_arguments_not_in_call]) + call <- as.call(call) + } + if(evaluate) eval(call, parent.frame()) + else call +} + +#' @export +#' @rdname saem +#' @param ci Should a matrix with estimates and confidence interval boundaries +#' be returned? If FALSE (default), a vector of estimates is returned. +parms.saem.mmkin <- function(object, ci = FALSE, ...) { + cov.mod <- object$sm@covariance.model + n_cov_mod_parms <- sum(cov.mod[upper.tri(cov.mod, diag = TRUE)]) + n_parms <- length(object$sm@name.modpar) + + n_cov_mod_parms + + length(object$sm@name.sigma) + + if (inherits(object$so, "try-error")) { + conf.int <- matrix(rep(NA, 3 * n_parms), ncol = 3) + colnames(conf.int) <- c("estimate", "lower", "upper") + } else { + conf.int <- object$so@results@conf.int[c("estimate", "lower", "upper")] + rownames(conf.int) <- object$so@results@conf.int[["name"]] + conf.int.var <- grepl("^Var\\.", rownames(conf.int)) + conf.int <- conf.int[!conf.int.var, ] + conf.int.cov <- grepl("^Cov\\.", rownames(conf.int)) + conf.int <- conf.int[!conf.int.cov, ] + } + estimate <- conf.int[, "estimate"] + + names(estimate) <- rownames(conf.int) + + if (ci) return(conf.int) + else return(estimate) +} |