diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/aw.R | 36 | ||||
| -rw-r--r-- | R/multistart.R | 59 | ||||
| -rw-r--r-- | R/saem.R | 8 |
3 files changed, 90 insertions, 13 deletions
@@ -30,6 +30,14 @@ #' @export aw <- function(object, ...) UseMethod("aw") +.aw <- function(all_objects) { + AIC_all <- sapply(all_objects, AIC) + delta_i <- AIC_all - min(AIC_all) + denom <- sum(exp(-delta_i/2)) + w_i <- exp(-delta_i/2) / denom + return(w_i) +} + #' @export #' @rdname aw aw.mkinfit <- function(object, ...) { @@ -43,11 +51,7 @@ aw.mkinfit <- function(object, ...) { } } all_objects <- list(object, ...) - AIC_all <- sapply(all_objects, AIC) - delta_i <- AIC_all - min(AIC_all) - denom <- sum(exp(-delta_i/2)) - w_i <- exp(-delta_i/2) / denom - return(w_i) + .aw(all_objects) } #' @export @@ -56,3 +60,25 @@ aw.mmkin <- function(object, ...) { if (ncol(object) > 1) stop("Please supply an mmkin column object") do.call(aw, object) } + +#' @export +#' @rdname aw +aw.mixed.mmkin <- function(object, ...) { + oo <- list(...) + data_object <- object$data[c("ds", "name", "time", "value")] + for (i in seq_along(oo)) { + if (!inherits(oo[[i]], "mixed.mmkin")) stop("Please supply objects inheriting from mixed.mmkin") + data_other_object <- oo[[i]]$data[c("ds", "name", "time", "value")] + if (!identical(data_object, data_other_object)) { + stop("It seems that the mixed.mmkin objects have not all been fitted to the same data") + } + } + all_objects <- list(object, ...) + .aw(all_objects) +} + +#' @export +#' @rdname aw +aw.multistart <- function(object, ...) { + do.call(aw, object) +} diff --git a/R/multistart.R b/R/multistart.R index a3afa08b..94292e82 100644 --- a/R/multistart.R +++ b/R/multistart.R @@ -14,7 +14,12 @@ #' @param n How many different combinations of starting parameters should be #' used? #' @param cores How many fits should be run in parallel? -#' @param \dots Passed to the update function. +#' @param \dots Passed to the update function, or to the basic plotting +#' function in the case of the graphical function. +#' @param x The multistart object to print +#' @param breaks Passed to [hist] +#' @param main title of the plot +#' @param lpos Positioning of the legend. #' @return A list of [saem.mmkin] objects, with class attributes #' 'multistart.saem.mmkin' and 'multistart'. #' @@ -39,6 +44,7 @@ multistart.saem.mmkin <- function(object, n = 50, cores = 1, ...) { res <- parallel::mclapply(1:n, function(x) { update(object, degparms_start = start_parms[x, ], ...) }, mc.cores = cores) + attr(res, "orig") <- object attr(res, "start_parms") <- start_parms class(res) <- c("multistart.saem.mmkin", "multistart") return(res) @@ -53,8 +59,53 @@ print.multistart <- function(x, ...) { #' @rdname multistart #' @export -summary.multistart.saem.mmkin <- function(object) { +parms.multistart <- function(object, ...) { + t(sapply(object, parms)) +} + +#' @rdname multistart +#' @importFrom stats median +#' @export +parhist <- function(object, lpos = "topleft", ...) { + orig <- attr(object, "orig") + orig_parms <- parms(orig) + start_parms <- orig$mean_dp_start + all_parms <- parms(object) + median_parms <- apply(all_parms, 2, median) + all_scaled_parms <- t(apply(all_parms, 1, function(x) x / median_parms)) + orig_scaled_parms <- orig_parms / median_parms + start_scaled_parms <- rep(NA_real_, length(orig_parms)) + names(start_scaled_parms) <- names(orig_parms) + start_scaled_parms[names(start_parms)] <- + start_parms / median_parms[names(start_parms)] + + boxplot(all_scaled_parms, log = "y", ...) + points(orig_scaled_parms, col = 2, cex = 2) + points(start_scaled_parms, col = 3, cex = 3) + legend(lpos, inset = c(0.05, 0.05), bty = "n", + pch = 1, col = 3:1, lty = c(NA, NA, 1), + legend = c( + "Starting parameters", + "Converged parameters", + "Multistart runs")) +} + +#' @rdname multistart +#' @importFrom KernSmooth bkde +#' @export +llhist <- function(object, breaks = "Sturges", main = "", lpos = "topleft", ...) { + ll <- sapply(object, logLik) + kde <- KernSmooth::bkde(ll) + h <- hist(ll, freq = TRUE, + xlim = range(kde$x), + xlab = "", main = main, + ylab = "Frequency of log likelihoods", breaks = breaks, ...) - parm_matrix <- sapply(object, parms) - parm_matrix + freq_factor <- h$counts[1] / h$density[1] + lines(kde$x, freq_factor * kde$y) + abline(v = logLik(attr(object, "orig")), col = 2) + legend(lpos, inset = c(0.05, 0.05), bty = "n", + lty = 1, col = c(2, 1), + legend = c("original log likelihood", + "kernel density estimate")) } @@ -568,7 +568,7 @@ saemix_data <- function(object, verbose = FALSE, ...) { 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, @@ -617,9 +617,9 @@ update.saem.mmkin <- function(object, ..., evaluate = TRUE) { #' @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(x, ci = FALSE, ...) { - conf.int <- x$so@results@conf.int[c("estimate", "lower", "upper")] - rownames(conf.int) <- x$so@results@conf.int[["name"]] +parms.saem.mmkin <- function(object, ci = FALSE, ...) { + 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, ] estimate <- conf.int[, "estimate"] |