From 20b9c584e7c43ecbb708459e531c24a1a4751e17 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Sat, 9 Nov 2019 01:05:51 +0100 Subject: Add a lack-of-fit test - Switch an example dataset in the test setup to a dataset with replicates, adapt tests - Skip the test for lrtest with an update specification as it does not only fail when pkgdown generates static help pages, but also in testthat --- test.log | 26 +++++++++++++++----------- 1 file changed, 15 insertions(+), 11 deletions(-) (limited to 'test.log') diff --git a/test.log b/test.log index 806c72b3..bc6d26ae 100644 --- a/test.log +++ b/test.log @@ -2,32 +2,36 @@ Loading mkin Testing mkin ✔ | OK F W S | Context ✔ | 2 | Export dataset for reading into CAKE -✔ | 10 | Confidence intervals and p-values [9.7 s] -✔ | 14 | Error model fitting [36.5 s] +✔ | 10 | Confidence intervals and p-values [10.1 s] +✔ | 14 | Error model fitting [40.5 s] ✔ | 4 | Calculation of FOCUS chi2 error levels [2.2 s] -✔ | 13 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [3.3 s] +✔ | 13 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [3.4 s] ✔ | 6 | Test fitting the decline of metabolites from their maximum [0.7 s] ✔ | 1 | Fitting the logistic model [0.9 s] ✔ | 1 | Test dataset class mkinds used in gmkin ✔ | 12 | Special cases of mkinfit calls [2.4 s] ✔ | 9 | mkinmod model generation and printing [0.2 s] ✔ | 3 | Model predictions with mkinpredict [0.3 s] -✔ | 16 | Evaluations according to 2015 NAFTA guidance [4.0 s] +✔ | 16 | Evaluations according to 2015 NAFTA guidance [4.1 s] ✔ | 4 | Calculation of maximum time weighted average concentrations (TWAs) [2.3 s] ✔ | 3 | Summary ✔ | 11 | Plotting [0.6 s] ✔ | 4 | AIC calculation ✔ | 2 | Residuals extracted from mkinfit models -✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [5.3 s] -✔ | 4 | Fitting the SFORB model [1.7 s] +✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [5.6 s] +✔ | 4 | Fitting the SFORB model [1.8 s] ✔ | 1 | Summaries of old mkinfit objects -✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [7.1 s] -✔ | 6 | Hypothesis tests [31.2 s] +✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [7.5 s] +✔ | 7 1 | Hypothesis tests [34.1 s] +──────────────────────────────────────────────────────────────────────────────── +test_tests.R:59: skip: We can do a likelihood ratio test using an update specification +Reason: This errors out if called by testthat while it works in a normal R session +──────────────────────────────────────────────────────────────────────────────── ══ Results ═════════════════════════════════════════════════════════════════════ -Duration: 108.3 s +Duration: 116.9 s -OK: 132 +OK: 133 Failed: 0 Warnings: 0 -Skipped: 0 +Skipped: 1 -- cgit v1.2.3 From c3700bec3a704660d3ade7a54c56b7084beb02b4 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 13 Nov 2019 21:15:35 +0100 Subject: Calculate Akaike weights --- NAMESPACE | 3 + NEWS.md | 2 + R/aw.R | 60 +++++++++++ _pkgdown.yml | 1 + docs/news/index.html | 1 + docs/reference/aw.html | 214 ++++++++++++++++++++++++++++++++++++++++ docs/reference/index.html | 6 ++ docs/sitemap.xml | 3 + man/aw.Rd | 47 +++++++++ test.log | 21 ++-- tests/testthat/FOCUS_2006_D.csf | 2 +- tests/testthat/test_aw.R | 12 +++ 12 files changed, 361 insertions(+), 11 deletions(-) create mode 100644 R/aw.R create mode 100644 docs/reference/aw.html create mode 100644 man/aw.Rd create mode 100644 tests/testthat/test_aw.R (limited to 'test.log') diff --git a/NAMESPACE b/NAMESPACE index e561621b..26995055 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,6 +3,8 @@ S3method("[",mmkin) S3method(AIC,mmkin) S3method(BIC,mmkin) +S3method(aw,mkinfit) +S3method(aw,mmkin) S3method(confint,mkinfit) S3method(loftest,mkinfit) S3method(logLik,mkinfit) @@ -30,6 +32,7 @@ export(IORE.solution) export(SFO.solution) export(SFORB.solution) export(add_err) +export(aw) export(backtransform_odeparms) export(endpoints) export(ilr) diff --git a/NEWS.md b/NEWS.md index 965105f4..28cf76ad 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,7 @@ # mkin 0.9.49.8 (unreleased) +- 'aw': Generic function for calculating Akaike weights, methods for mkinfit objects and mmkin columns + - 'loftest': Add a lack-of-fit test - 'plot_res', 'plot_sep' and 'mkinerrplot': Add the possibility to show standardized residuals and make it the default for fits with error models other than 'const' diff --git a/R/aw.R b/R/aw.R new file mode 100644 index 00000000..24078f38 --- /dev/null +++ b/R/aw.R @@ -0,0 +1,60 @@ +#' Calculate Akaike weights for model averaging +#' +#' Akaike weights are calculated based on the relative +#' expected Kullback-Leibler information as specified +#' by Burnham and Anderson (2004). +#' +#' @param object An mmkin column object, containing two or more +#' \code{\link{mkinfit}} models that have been fitted to the same data, +#' or an mkinfit object. In the latter case, further mkinfit +#' objects fitted to the same data should be specified +#' as dots arguments. +#' @param \dots Not used in the method for mmkin column objects, +#' further mkinfit objects in the method for mkinfit objects. +#' @references Burnham KP and Anderson DR (2004) Multimodel +#' Inference: Understanding AIC and BIC in Model Selection +#' Sociological Methods & Research 33(2) 261-304 +#' @examples +#' \dontrun{ +#' f_sfo <- mkinfit("SFO", FOCUS_2006_D, quiet = TRUE) +#' f_dfop <- mkinfit("DFOP", FOCUS_2006_D, quiet = TRUE) +#' aw_sfo_dfop <- aw(f_sfo, f_dfop) +#' sum(aw_sfo_dfop) +#' aw_sfo_dfop # SFO gets more weight as it has less parameters and a similar fit +#' f <- mmkin(c("SFO", "FOMC", "DFOP"), list("FOCUS D" = FOCUS_2006_D), cores = 1, quiet = TRUE) +#' aw(f) +#' sum(aw(f)) +#' aw(f[c("SFO", "DFOP")]) +#' } +#' @export +aw <- function(object, ...) UseMethod("aw") + +#' @export +#' @rdname aw +aw.mkinfit <- function(object, ...) { + oo <- list(...) + data_object <- object$data[c("time", "variable", "observed")] + for (i in seq_along(oo)) { + if (!inherits(oo[[i]], "mkinfit")) stop("Please supply only mkinfit objects") + data_other_object <- oo[[i]]$data[c("time", "variable", "observed")] + if (!identical(data_object, data_other_object)) { + stop("It seems that the mkinfit objects have not all been fitted to the same data") + } + } + 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) +} + +#' @export +#' @rdname aw +aw.mmkin <- function(object, ...) { + if (ncol(object) > 1) stop("Please supply an mmkin column object") + do.call(aw, object) +} + + + diff --git a/_pkgdown.yml b/_pkgdown.yml index c298256f..6bb05b3e 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -25,6 +25,7 @@ reference: - loftest - mkinerrmin - endpoints + - aw - CAKE_export - title: Work with mmkin objects desc: Functions working with aggregated results diff --git a/docs/news/index.html b/docs/news/index.html index 9aa2e18b..6b0b89fa 100644 --- a/docs/news/index.html +++ b/docs/news/index.html @@ -134,6 +134,7 @@ mkin 0.9.49.8 (unreleased) Unreleased