diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/confint.mkinfit.R | 86 | ||||
| -rw-r--r-- | R/mmkin.R | 5 | ||||
| -rw-r--r-- | R/nlme.R | 2 | ||||
| -rw-r--r-- | R/parms.mkinfit.R | 4 | ||||
| -rw-r--r-- | R/saemix.R | 31 |
5 files changed, 78 insertions, 50 deletions
diff --git a/R/confint.mkinfit.R b/R/confint.mkinfit.R index 78dda95d..53eb45ee 100644 --- a/R/confint.mkinfit.R +++ b/R/confint.mkinfit.R @@ -27,7 +27,10 @@ #' @param backtransform If we approximate the likelihood in terms of the #' transformed parameters, should we backtransform the parameters with #' their confidence intervals? -#' @param cores The number of cores to be used for multicore processing. +#' @param rel_tol If the method is 'profile', what should be the accuracy +#' of the lower and upper bounds, relative to the estimate obtained from +#' the quadratic method? +#' @param cores The number of cores to be used for multicore processing. #' On Windows machines, cores > 1 is currently not supported. #' @param quiet Should we suppress the message "Profiling the likelihood" #' @return A matrix with columns giving lower and upper confidence limits for @@ -121,7 +124,7 @@ confint.mkinfit <- function(object, parm, level = 0.95, alpha = 1 - level, cutoff, method = c("quadratic", "profile"), transformed = TRUE, backtransform = TRUE, - cores = round(detectCores()/2), quiet = FALSE, ...) + cores = parallel::detectCores(), rel_tol = 0.01, quiet = FALSE, ...) { tparms <- parms(object, transformed = TRUE) bparms <- parms(object, transformed = FALSE) @@ -140,50 +143,50 @@ confint.mkinfit <- function(object, parm, a <- c(alpha / 2, 1 - (alpha / 2)) - if (method == "quadratic") { + quantiles <- qt(a, object$df.residual) - quantiles <- qt(a, object$df.residual) - - covar_pnames <- if (missing(parm)) { - if (transformed) tpnames else bpnames - } else { - parm - } + covar_pnames <- if (missing(parm)) { + if (transformed) tpnames else bpnames + } else { + parm + } - return_parms <- if (backtransform) bparms[return_pnames] - else tparms[return_pnames] + return_parms <- if (backtransform) bparms[return_pnames] + else tparms[return_pnames] - covar_parms <- if (transformed) tparms[covar_pnames] - else bparms[covar_pnames] + covar_parms <- if (transformed) tparms[covar_pnames] + else bparms[covar_pnames] - if (transformed) { - covar <- try(solve(object$hessian), silent = TRUE) - } else { - covar <- try(solve(object$hessian_notrans), silent = TRUE) - } + if (transformed) { + covar <- try(solve(object$hessian), silent = TRUE) + } else { + covar <- try(solve(object$hessian_notrans), silent = TRUE) + } - # If inverting the covariance matrix failed or produced NA values - if (!is.numeric(covar) | is.na(covar[1])) { - ses <- lci <- uci <- rep(NA, p) - } else { - ses <- sqrt(diag(covar))[covar_pnames] - lci <- covar_parms + quantiles[1] * ses - uci <- covar_parms + quantiles[2] * ses - if (transformed & backtransform) { - lci_back <- backtransform_odeparms(lci, - object$mkinmod, object$transform_rates, object$transform_fractions) - uci_back <- backtransform_odeparms(uci, - object$mkinmod, object$transform_rates, object$transform_fractions) - - return_errparm_names <- intersect(names(object$errparms), return_pnames) - lci <- c(lci_back, lci[return_errparm_names]) - uci <- c(uci_back, uci[return_errparm_names]) - } + # If inverting the covariance matrix failed or produced NA values + if (!is.numeric(covar) | is.na(covar[1])) { + ses <- lci <- uci <- rep(NA, p) + } else { + ses <- sqrt(diag(covar))[covar_pnames] + lci <- covar_parms + quantiles[1] * ses + uci <- covar_parms + quantiles[2] * ses + if (transformed & backtransform) { + lci_back <- backtransform_odeparms(lci, + object$mkinmod, object$transform_rates, object$transform_fractions) + uci_back <- backtransform_odeparms(uci, + object$mkinmod, object$transform_rates, object$transform_fractions) + + return_errparm_names <- intersect(names(object$errparms), return_pnames) + lci <- c(lci_back, lci[return_errparm_names]) + uci <- c(uci_back, uci[return_errparm_names]) } - ci <- cbind(lower = lci, upper = uci) } + ci <- cbind(lower = lci, upper = uci) if (method == "profile") { + + ci_quadratic <- ci + if (!quiet) message("Profiling the likelihood") lci <- uci <- rep(NA, p) @@ -215,9 +218,14 @@ confint.mkinfit <- function(object, parm, (cutoff - (object$logLik - profile_ll(x)))^2 } - lci_pname <- optimize(cost, lower = 0, upper = all_parms[pname])$minimum + lower_quadratic <- ci_quadratic["lower"][pname] + upper_quadratic <- ci_quadratic["upper"][pname] + ltol <- if (!is.na(lower_quadratic)) rel_tol * lower_quadratic else .Machine$double.eps^0.25 + utol <- if (!is.na(upper_quadratic)) rel_tol * upper_quadratic else .Machine$double.eps^0.25 + lci_pname <- optimize(cost, lower = 0, upper = all_parms[pname], tol = ltol)$minimum uci_pname <- optimize(cost, lower = all_parms[pname], - upper = ifelse(grepl("^f_|^g$", pname), 1, 15 * all_parms[pname]))$minimum + upper = ifelse(grepl("^f_|^g$", pname), 1, 15 * all_parms[pname]), + tol = utol)$minimum return(c(lci_pname, uci_pname)) } ci <- t(parallel::mcmapply(get_ci, profile_pnames, mc.cores = cores)) @@ -12,7 +12,8 @@ #' @param cores The number of cores to be used for multicore processing. This #' is only used when the \code{cluster} argument is \code{NULL}. On Windows #' machines, cores > 1 is not supported, you need to use the \code{cluster} -#' argument to use multiple logical processors. +#' argument to use multiple logical processors. Per default, all cores +#' detected by [parallel::detectCores()] are used. #' @param cluster A cluster as returned by \code{\link{makeCluster}} to be used #' for parallel execution. #' @param \dots Further arguments that will be passed to \code{\link{mkinfit}}. @@ -62,7 +63,7 @@ #' #' @export mmkin mmkin <- function(models = c("SFO", "FOMC", "DFOP"), datasets, - cores = round(detectCores()/2), cluster = NULL, ...) + cores = detectCores(), cluster = NULL, ...) { parent_models_available = c("SFO", "FOMC", "DFOP", "HS", "SFORB", "IORE", "logistic") n.m <- length(models) @@ -125,7 +125,7 @@ nlme_function <- function(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? +#' 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) { diff --git a/R/parms.mkinfit.R b/R/parms.mkinfit.R index aae6fa52..a1f2d209 100644 --- a/R/parms.mkinfit.R +++ b/R/parms.mkinfit.R @@ -21,11 +21,13 @@ #' ds <- lapply(experimental_data_for_UBA_2019[6:10], #' function(x) subset(x$data[c("name", "time", "value")])) #' names(ds) <- paste("Dataset", 6:10) -#' fits <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE) +#' \dontrun{ +#' fits <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE, cores = 1) #' parms(fits["SFO", ]) #' parms(fits[, 2]) #' parms(fits) #' parms(fits, transformed = TRUE) +#' } #' @export parms <- function(object, ...) { @@ -5,25 +5,37 @@ #' list of mkinfit objects that have been obtained by fitting the same model to #' a list of datasets. #' +#' Starting values for the fixed effects (population mean parameters, argument psi0 of +#' [saemix::saemixModel()] are the mean values of the parameters found using +#' mmkin. Starting variances of the random effects (argument omega.init) are the +#' variances of the deviations of the parameters from these mean values. +#' #' @param object An mmkin row object containing several fits of the same model to different datasets +#' @param cores The number of cores to be used for multicore processing. +#' On Windows machines, cores > 1 is currently not supported. #' @rdname saemix #' @importFrom saemix saemixData saemixModel +#' @importFrom stats var #' @examples #' ds <- lapply(experimental_data_for_UBA_2019[6:10], #' function(x) subset(x$data[c("name", "time", "value")])) #' names(ds) <- paste("Dataset", 6:10) #' sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"), #' A1 = mkinsub("SFO")) -#' f_mmkin <- mmkin(list("SFO-SFO" = sfo_sfo), ds, quiet = TRUE, cores = 5) +#' \dontrun{ +#' f_mmkin <- mmkin(list("SFO-SFO" = sfo_sfo), ds, quiet = TRUE) +#' library(saemix) #' m_saemix <- saemix_model(f_mmkin) #' d_saemix <- saemix_data(f_mmkin) -#' saemix_options <- list(seed = 123456, save = FALSE, save.graphs = FALSE) -#' \dontrun{ -#' saemix(m_saemix, d_saemix, saemix_options) +#' saemix_options <- list(seed = 123456, +#' save = FALSE, save.graphs = FALSE, displayProgress = FALSE, +#' nbiter.saemix = c(200, 80)) +#' f_saemix <- saemix(m_saemix, d_saemix, saemix_options) +#' plot(f_saemix, plot.type = "convergence") #' } #' @return An [saemix::SaemixModel] object. #' @export -saemix_model <- function(object) { +saemix_model <- function(object, cores = parallel::detectCores()) { if (nrow(object) > 1) stop("Only row objects allowed") mkin_model <- object[[1]]$mkinmod @@ -81,14 +93,19 @@ saemix_model <- function(object) { out_values <- out_wide[out_index] } return(out_values) - }, mc.cores = 15) + }, mc.cores = cores) res <- unlist(res_list) return(res) } + raneff_0 <- mean_degparms(object, random = TRUE)$random$ds + var_raneff_0 <- apply(raneff_0, 2, var) + res <- saemixModel(model_function, psi0, "Mixed model generated from mmkin object", - transform.par = rep(0, length(degparms_optim))) + transform.par = rep(0, length(degparms_optim)), + omega.init = diag(var_raneff_0) + ) return(res) } |