diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/create_deg_func.R | 38 | ||||
| -rw-r--r-- | R/illparms.R | 10 | ||||
| -rw-r--r-- | R/mhmkin.R | 28 | ||||
| -rw-r--r-- | R/mkinfit.R | 6 | ||||
| -rw-r--r-- | R/mkinpredict.R | 4 | ||||
| -rw-r--r-- | R/multistart.R | 4 | ||||
| -rw-r--r-- | R/nlme.R | 2 | ||||
| -rw-r--r-- | R/parplot.R | 19 | ||||
| -rw-r--r-- | R/plot.mixed.mmkin.R | 2 | ||||
| -rw-r--r-- | R/status.R | 18 |
10 files changed, 88 insertions, 43 deletions
diff --git a/R/create_deg_func.R b/R/create_deg_func.R index 5794c65c..d3e11f78 100644 --- a/R/create_deg_func.R +++ b/R/create_deg_func.R @@ -61,16 +61,21 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) { ), ")") - if (length(obs_vars) >= 2) { - supported <- FALSE # except for the following cases + if (length(obs_vars) >= 2) supported <- FALSE + # Except for the following cases: + + if (length(obs_vars) == 2) { n1 <- obs_vars[1] n2 <- obs_vars[2] n10 <- paste0("odeini['", parent, "']") n20 <- paste0("odeini['", n2, "']") # sfo_sfo - if (all(spec[[1]]$sink == FALSE, length(obs_vars) == 2, - spec[[1]]$type == "SFO", spec[[2]]$type == "SFO")) { + if (all( + spec[[1]]$sink == FALSE, + spec[[1]]$type == "SFO", + spec[[2]]$type == "SFO", + is.null(spec[[2]]$to))) { supported <- TRUE k1 <- k_parent k2 <- paste0("k_", n2, if(min_ff) "_sink" else "") @@ -80,8 +85,12 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) { } # sfo_f12_sfo - if (all(use_of_ff == "max", spec[[1]]$sink == TRUE, length(obs_vars) == 2, - spec[[1]]$type == "SFO", spec[[2]]$type == "SFO")) { + if (all( + use_of_ff == "max", + spec[[1]]$sink == TRUE, + spec[[1]]$type == "SFO", + spec[[2]]$type == "SFO", + is.null(spec[[2]]$to))) { supported <- TRUE k1 <- k_parent k2 <- paste0("k_", n2) @@ -92,8 +101,12 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) { } # sfo_k120_sfo - if (all(use_of_ff == "min", spec[[1]]$sink == TRUE, length(obs_vars) == 2, - spec[[1]]$type == "SFO", spec[[2]]$type == "SFO")) { + if (all( + use_of_ff == "min", + spec[[1]]$sink == TRUE, + spec[[1]]$type == "SFO", + spec[[2]]$type == "SFO", + is.null(spec[[2]]$to))) { supported <- TRUE k12 <- paste0("k_", n1, "_", n2) k10 <- paste0("k_", n1, "_sink") @@ -104,8 +117,12 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) { } # dfop_f12_sfo - if (all(use_of_ff == "max", spec[[1]]$sink == TRUE, length(obs_vars) == 2, - spec[[1]]$type == "DFOP", spec[[2]]$type == "SFO")) { + if (all( + use_of_ff == "max", + spec[[1]]$sink == TRUE, + spec[[1]]$type == "DFOP", + spec[[2]]$type == "SFO", + is.null(spec[[2]]$to))) { supported <- TRUE f12 <- paste0("f_", n1, "_to_", n2) k2 <- paste0("k_", n2) @@ -119,7 +136,6 @@ create_deg_func <- function(spec, use_of_ff = c("min", "max")) { } - if (supported) { deg_func <- function(observed, odeini, odeparms) {} diff --git a/R/illparms.R b/R/illparms.R index 68a6bff6..b4b37fbb 100644 --- a/R/illparms.R +++ b/R/illparms.R @@ -102,12 +102,14 @@ illparms.saem.mmkin <- function(object, conf.level = 0.95, random = TRUE, errmod ints <- intervals(object, conf.level = conf.level) ill_parms <- character(0) if (random) { - ill_parms_random <- ints$random[, "lower"] < 0 - ill_parms <- c(ill_parms, names(which(ill_parms_random))) + ill_parms_random_i <- which(ints$random[, "lower"] < 0) + ill_parms_random <- rownames(ints$random)[ill_parms_random_i] + ill_parms <- c(ill_parms, ill_parms_random) } if (errmod) { - ill_parms_errmod <- ints$errmod[, "lower"] < 0 & ints$errmod[, "upper"] > 0 - ill_parms <- c(ill_parms, names(which(ill_parms_errmod))) + ill_parms_errmod_i <- which(ints$errmod[, "lower"] < 0 & ints$errmod[, "upper"] > 0) + ill_parms_errmod <- rownames(ints$errmod)[ill_parms_errmod_i] + ill_parms <- c(ill_parms, ill_parms_errmod) } if (slopes) { if (is.null(object$so)) stop("Slope testing is only implemented for the saemix backend") @@ -219,11 +219,22 @@ print.mhmkin <- function(x, ...) { print(status(x)) } +#' Check if fit within an mhmkin object failed +#' @param x The object to be checked +check_failed <- function(x) { + if (inherits(x, "try-error")) { + return(TRUE) + } else { + if (inherits(x$so, "try-error")) { + return(TRUE) + } else { + return(FALSE) + } + } +} + #' @export AIC.mhmkin <- function(object, ..., k = 2) { - if (inherits(object[[1]], "saem.mmkin")) { - check_failed <- function(x) if (inherits(x$so, "try-error")) TRUE else FALSE - } res <- sapply(object, function(x) { if (check_failed(x)) return(NA) else return(AIC(x$so, k = k)) @@ -235,9 +246,6 @@ AIC.mhmkin <- function(object, ..., k = 2) { #' @export BIC.mhmkin <- function(object, ...) { - if (inherits(object[[1]], "saem.mmkin")) { - check_failed <- function(x) if (inherits(x$so, "try-error")) TRUE else FALSE - } res <- sapply(object, function(x) { if (check_failed(x)) return(NA) else return(BIC(x$so)) @@ -280,7 +288,13 @@ anova.mhmkin <- function(object, ..., if (identical(model.names, "auto")) { model.names <- outer(rownames(object), colnames(object), paste) } - rlang::inject(anova(!!!(object), method = method, test = test, + failed_index <- which(sapply(object, check_failed), arr.ind = TRUE) + if (length(failed_index > 0)) { + rlang::inject(anova(!!!(object[-failed_index]), method = method, test = test, + model.names = model.names[-failed_index])) + } else { + rlang::inject(anova(!!!(object), method = method, test = test, model.names = model.names)) + } } diff --git a/R/mkinfit.R b/R/mkinfit.R index c851fddb..52053685 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -21,7 +21,8 @@ utils::globalVariables(c("name", "time", "value")) #' "FOMC", "DFOP", "HS", "SFORB", "IORE"). If a shorthand name is given, a #' parent only degradation model is generated for the variable with the #' highest value in \code{observed}. -#' @param observed A dataframe with the observed data. The first column called +#' @param observed A dataframe or an object coercible to a dataframe +#' (e.g. a \code{tibble}) with the observed data. The first column called #' "name" must contain the name of the observed variable for each data point. #' The second column must contain the times of observation, named "time". #' The third column must be named "value" and contain the observed values. @@ -292,6 +293,9 @@ mkinfit <- function(mkinmod, observed, # Get the names of observed variables obs_vars <- names(mkinmod$spec) + # Coerce observed data to a dataframe + observed <- as.data.frame(observed) + # Subset observed data with names of observed data in the model and remove NA values observed <- subset(observed, name %in% obs_vars) observed <- subset(observed, !is.na(value)) diff --git a/R/mkinpredict.R b/R/mkinpredict.R index 60456fb2..4c6d7862 100644 --- a/R/mkinpredict.R +++ b/R/mkinpredict.R @@ -22,7 +22,7 @@ #' variables. The third possibility "eigen" is fast in comparison to uncompiled #' ODE models, but not applicable to some models, e.g. using FOMC for the #' parent compound. -#' @param method.ode The solution method passed via [mkinpredict] to [ode]] in +#' @param method.ode The solution method passed via [mkinpredict] to `deSolve::ode()` in #' case the solution type is "deSolve" and we are not using compiled code. #' When using compiled code, only lsoda is supported. #' @param use_compiled If set to \code{FALSE}, no compiled version of the @@ -36,7 +36,7 @@ #' the observed variables (default) or for all state variables (if set to #' FALSE). Setting this to FALSE has no effect for analytical solutions, #' as these always return mapped output. -#' @param na_stop Should it be an error if [ode] returns NaN values +#' @param na_stop Should it be an error if `deSolve::ode()` returns NaN values #' @param \dots Further arguments passed to the ode solver in case such a #' solver is used. #' @return A matrix with the numeric solution in wide format diff --git a/R/multistart.R b/R/multistart.R index aeea2d81..4ba82a43 100644 --- a/R/multistart.R +++ b/R/multistart.R @@ -40,7 +40,7 @@ #' f_mmkin <- mmkin("DFOP", dmta_ds, error_model = "tc", cores = 7, quiet = TRUE) #' f_saem_full <- saem(f_mmkin) #' f_saem_full_multi <- multistart(f_saem_full, n = 16, cores = 16) -#' parplot(f_saem_full_multi, lpos = "topleft") +#' parplot(f_saem_full_multi, lpos = "topleft", las = 2) #' illparms(f_saem_full) #' #' f_saem_reduced <- update(f_saem_full, no_random_effect = "log_k2") @@ -50,7 +50,7 @@ #' library(parallel) #' cl <- makePSOCKcluster(12) #' f_saem_reduced_multi <- multistart(f_saem_reduced, n = 16, cluster = cl) -#' parplot(f_saem_reduced_multi, lpos = "topright", ylim = c(0.5, 2)) +#' parplot(f_saem_reduced_multi, lpos = "topright", ylim = c(0.5, 2), las = 2) #' stopCluster(cl) #' } multistart <- function(object, n = 50, @@ -126,7 +126,7 @@ nlme_function <- function(object) { #' @rdname nlme #' @importFrom rlang !!! -#' @return A \code{\link{groupedData}} object +#' @return A `nlme::groupedData` object #' @export nlme_data <- function(object) { if (nrow(object) > 1) stop("Only row objects allowed") diff --git a/R/parplot.R b/R/parplot.R index 3da4b51a..a33112a5 100644 --- a/R/parplot.R +++ b/R/parplot.R @@ -35,9 +35,6 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, llquant = NA, scale = c("best", "median"), lpos = "bottomleft", main = "", ...) { - oldpar <- par(no.readonly = TRUE) - on.exit(par(oldpar, no.readonly = TRUE)) - orig <- attr(object, "orig") orig_parms <- parms(orig) start_degparms <- orig$mean_dp_start @@ -59,11 +56,10 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, llquant = NA, selected <- which(ll > llmin) selected_parms <- all_parms[selected, ] - par(las = 1) if (orig$transformations == "mkin") { degparm_names_transformed <- names(start_degparms) degparm_index <- which(names(orig_parms) %in% degparm_names_transformed) - orig_parms[degparm_names_transformed] <- backtransform_odeparms( + orig_degparms <- backtransform_odeparms( orig_parms[degparm_names_transformed], orig$mmkin[[1]]$mkinmod, transform_rates = orig$mmkin[[1]]$transform_rates, @@ -74,14 +70,17 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, llquant = NA, transform_fractions = orig$mmkin[[1]]$transform_fractions) degparm_names <- names(start_degparms) - names(orig_parms) <- c(degparm_names, names(orig_parms[-degparm_index])) + orig_parms_back <- orig_parms + orig_parms_back[degparm_index] <- orig_degparms + names(orig_parms_back)[degparm_index] <- degparm_names + orig_parms <- orig_parms_back selected_parms[, degparm_names_transformed] <- t(apply(selected_parms[, degparm_names_transformed], 1, backtransform_odeparms, orig$mmkin[[1]]$mkinmod, transform_rates = orig$mmkin[[1]]$transform_rates, transform_fractions = orig$mmkin[[1]]$transform_fractions)) - colnames(selected_parms)[1:length(degparm_names)] <- degparm_names + colnames(selected_parms)[degparm_index] <- degparm_names } start_errparms <- orig$so@model@error.init @@ -99,6 +98,12 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, llquant = NA, # Boxplots of all scaled parameters selected_scaled_parms <- t(apply(selected_parms, 1, function(x) x / parm_scale)) + i_negative <- selected_scaled_parms <= 0 + parms_with_negative_scaled_values <- paste(names(which(apply(i_negative, 2, any))), collapse = ", ") + if (any(i_negative)) { + warning("There are negative values for ", parms_with_negative_scaled_values, " which are set to NA for plotting") + } + selected_scaled_parms[i_negative] <- NA boxplot(selected_scaled_parms, log = "y", main = main, , ylab = "Normalised parameters", ...) diff --git a/R/plot.mixed.mmkin.R b/R/plot.mixed.mmkin.R index d6c3d0de..f05f1110 100644 --- a/R/plot.mixed.mmkin.R +++ b/R/plot.mixed.mmkin.R @@ -93,7 +93,7 @@ plot.mixed.mmkin <- function(x, nrow.legend = ceiling((length(i) + 1) / ncol.legend), rel.height.legend = 0.02 + 0.07 * nrow.legend, rel.height.bottom = 1.1, - pch_ds = 1:length(i), + pch_ds = c(1:25, 33, 35:38, 40:41, 47:57, 60:90)[1:length(i)], col_ds = pch_ds + 1, lty_ds = col_ds, frame = TRUE, ... @@ -74,15 +74,19 @@ print.status.mmkin <- function(x, ...) { status.mhmkin <- function(object, ...) { if (inherits(object[[1]], "saem.mmkin")) { test_func <- function(fit) { - if (inherits(fit$so, "try-error")) { - return("E") + if (inherits(fit, "try-error")) { + return("E") } else { - if (!is.null(fit$FIM_failed)) { - return_values <- c("fixed effects" = "Fth", - "random effects and error model parameters" = "FO") - return(paste(return_values[fit$FIM_failed], collapse = ", ")) + if (inherits(fit$so, "try-error")) { + return("E") } else { - return("OK") + if (!is.null(fit$FIM_failed)) { + return_values <- c("fixed effects" = "Fth", + "random effects and error model parameters" = "FO") + return(paste(return_values[fit$FIM_failed], collapse = ", ")) + } else { + return("OK") + } } } } |