aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
Diffstat (limited to 'R')
-rw-r--r--R/aw.R36
-rw-r--r--R/multistart.R59
-rw-r--r--R/saem.R8
3 files changed, 90 insertions, 13 deletions
diff --git a/R/aw.R b/R/aw.R
index f46b20ec..b3992f94 100644
--- a/R/aw.R
+++ b/R/aw.R
@@ -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"))
}
diff --git a/R/saem.R b/R/saem.R
index 875456c3..0d0d9b8a 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -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"]

Contact - Imprint