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 --- tests/testthat/FOCUS_2006_D.csf | 2 +- tests/testthat/setup_script.R | 10 +++++----- tests/testthat/test_confidence.R | 7 ++++--- tests/testthat/test_tests.R | 22 ++++++++++++++++++++-- 4 files changed, 30 insertions(+), 11 deletions(-) (limited to 'tests') diff --git a/tests/testthat/FOCUS_2006_D.csf b/tests/testthat/FOCUS_2006_D.csf index 528e2b61..09940aa3 100644 --- a/tests/testthat/FOCUS_2006_D.csf +++ b/tests/testthat/FOCUS_2006_D.csf @@ -5,7 +5,7 @@ Description: MeasurementUnits: % AR TimeUnits: days Comments: Created using mkin::CAKE_export -Date: 2019-11-05 +Date: 2019-11-09 Optimiser: IRLS [Data] diff --git a/tests/testthat/setup_script.R b/tests/testthat/setup_script.R index 9becdd2a..e33f4af7 100644 --- a/tests/testthat/setup_script.R +++ b/tests/testthat/setup_script.R @@ -32,9 +32,6 @@ f_1_mkin_trans <- mkinfit("SFO", FOCUS_2006_A, quiet = TRUE) f_1_mkin_notrans <- mkinfit("SFO", FOCUS_2006_A, quiet = TRUE, transform_rates = FALSE) -f_2_mkin <- mkinfit("DFOP", FOCUS_2006_C, quiet = TRUE) -f_2_nls <- nls(value ~ SSbiexp(time, A1, lrc1, A2, lrc2), data = FOCUS_2006_C) - # mmkin object of parent fits for tests models <- c("SFO", "FOMC", "DFOP", "HS") fits <- mmkin(models, @@ -62,11 +59,14 @@ f_sfo_sfo.ff <- mkinfit(SFO_SFO.ff, subset(FOCUS_2006_D, value != 0), quiet = TRUE) -# Two metabolites SFO_lin_a <- synthetic_data_for_UBA_2014[[1]]$data - DFOP_par_c <- synthetic_data_for_UBA_2014[[12]]$data +f_2_mkin <- mkinfit("DFOP", DFOP_par_c, quiet = TRUE) +f_2_nls <- nls(value ~ SSbiexp(time, A1, lrc1, A2, lrc2), data = subset(DFOP_par_c, name == "parent")) +f_2_anova <- lm(value ~ as.factor(time), data = subset(DFOP_par_c, name == "parent")) + +# Two metabolites m_synth_SFO_lin <- mkinmod( parent = mkinsub("SFO", "M1"), M1 = mkinsub("SFO", "M2"), diff --git a/tests/testthat/test_confidence.R b/tests/testthat/test_confidence.R index a2bf1401..e85fdb7a 100644 --- a/tests/testthat/test_confidence.R +++ b/tests/testthat/test_confidence.R @@ -54,11 +54,12 @@ test_that("Quadratic confidence intervals for rate constants are comparable to v # Another case: se_mkin_2 <- summary(f_2_mkin)$par[1:4, "Std. Error"] se_nls_2 <- summary(f_2_nls)$coefficients[, "Std. Error"] - # Here we the ratio of standard errors can be explained by the same + # Here the ratio of standard errors can be explained by the same # principle up to about 3% + nobs_DFOP_par_c_parent <- nrow(subset(DFOP_par_c, name == "parent")) expect_equivalent( se_nls_2[c("lrc1", "lrc2")] / se_mkin_2[c("log_k1", "log_k2")], - rep(sqrt(nrow(FOCUS_2006_C) / (nrow(FOCUS_2006_C) - 4)), 2), + rep(sqrt(nobs_DFOP_par_c_parent / (nobs_DFOP_par_c_parent - 4)), 2), tolerance = 0.03) }) @@ -73,7 +74,7 @@ test_that("Likelihood profile based confidence intervals work", { } f_mle <- stats4::mle(f_nll, start = as.list(parms(f)), nobs = nrow(FOCUS_2006_C)) - ci_mkin_1_p_0.95 <- confint(f, method = "profile", level = 0.95, + ci_mkin_1_p_0.95 <- confint(f, method = "profile", level = 0.95, cores = n_cores, quiet = TRUE) # Magically, we get very similar boundaries as stats4::mle diff --git a/tests/testthat/test_tests.R b/tests/testthat/test_tests.R index 5a522f8e..bdc72f08 100644 --- a/tests/testthat/test_tests.R +++ b/tests/testthat/test_tests.R @@ -1,5 +1,20 @@ context("Hypothesis tests") +test_that("The lack-of-fit test works and can be reproduced using nls", { + + expect_error(loftest(f_1_mkin_trans), "Not defined for fits to data without replicates") + + loftest_mkin <- loftest(f_2_mkin) + + # This code is inspired by Ritz and Streibig (2008) Nonlinear Regression using R, p. 64 + Q <- as.numeric(- 2 * (logLik(f_2_nls) - logLik(f_2_anova))) + df.Q <- df.residual(f_2_nls) - df.residual(f_2_anova) + p_nls <- 1 - pchisq(Q, df.Q) + + expect_equal(loftest_mkin[["2", "Pr(>Chisq)"]], p_nls, tolerance = 1e-5) + +}) + test_that("The likelihood ratio test works", { expect_error(lrtest(f_1_mkin_trans, f_2_mkin), "not been fitted to the same data") @@ -25,7 +40,7 @@ test_that("Updating fitted models works", { parent = mkinsub("DFOP", to = "A1"), A1 = mkinsub("SFO", to = "A2"), A2 = mkinsub("SFO"), - use_of_ff = "max" + use_of_ff = "max", quiet = TRUE ) f_soil_1_tc <- mkinfit(dfop_sfo_sfo, @@ -41,6 +56,9 @@ test_that("Updating fitted models works", { }) test_that("We can do a likelihood ratio test using an update specification", { + skip("This errors out if called by testthat while it works in a normal R session") test_2_mkin_k2 <- lrtest(f_2_mkin, fixed_parms = c(k2 = 0)) - expect_equivalent(test_2_mkin_k2[["2", "Pr(>Chisq)"]], 1.139e-6, tolerance = 1e-8) + expect_equivalent(test_2_mkin_k2[["2", "Pr(>Chisq)"]], 4.851e-8, tolerance = 1e-8) + test_2_mkin_tc <- lrtest(f_2_mkin, error_model = "tc") + expect_equivalent(test_2_mkin_tc[["2", "Pr(>Chisq)"]], 7.302e-5, tolerance = 1e-7) }) -- cgit v1.2.3 From f3757dfbaabe4e60ecac063437a75ea32999640a Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 13 Nov 2019 11:04:27 +0100 Subject: Remove leftover line --- tests/testthat/test_SFORB.R | 2 -- 1 file changed, 2 deletions(-) (limited to 'tests') diff --git a/tests/testthat/test_SFORB.R b/tests/testthat/test_SFORB.R index 49b3beed..bc9ab646 100644 --- a/tests/testthat/test_SFORB.R +++ b/tests/testthat/test_SFORB.R @@ -18,8 +18,6 @@ context("Fitting the SFORB model") -logistic <- mkinmod(parent = mkinsub("logistic")) - test_that("Fitting the SFORB model is equivalent to fitting DFOP", { f_sforb <- mkinfit("SFORB", FOCUS_2006_C, quiet = TRUE) f_dfop <- mkinfit("DFOP", FOCUS_2006_C, quiet = TRUE) -- cgit v1.2.3 From b0e529ff49dfa52568fe7c5e6e76439a83c62840 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Wed, 13 Nov 2019 16:22:17 +0100 Subject: Improve citation --- tests/testthat/test_tests.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) (limited to 'tests') diff --git a/tests/testthat/test_tests.R b/tests/testthat/test_tests.R index bdc72f08..ddf8e1a0 100644 --- a/tests/testthat/test_tests.R +++ b/tests/testthat/test_tests.R @@ -6,7 +6,8 @@ test_that("The lack-of-fit test works and can be reproduced using nls", { loftest_mkin <- loftest(f_2_mkin) - # This code is inspired by Ritz and Streibig (2008) Nonlinear Regression using R, p. 64 + # This code is a slightly modified version of that given in Ritz and Streibig + # (2008) Nonlinear Regression using R, p. 64 Q <- as.numeric(- 2 * (logLik(f_2_nls) - logLik(f_2_anova))) df.Q <- df.residual(f_2_nls) - df.residual(f_2_anova) p_nls <- 1 - pchisq(Q, df.Q) -- 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 'tests') 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
    +

  • ‘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’

  • ‘lrtest.mkinfit’: Improve naming of the compared fits in the case of fixed parameters

    • diff --git a/docs/reference/aw.html b/docs/reference/aw.html new file mode 100644 index 00000000..b6f3ce48 --- /dev/null +++ b/docs/reference/aw.html @@ -0,0 +1,214 @@ + + + + + + + + +Calculate Akaike weights for model averaging — aw • mkin + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
    +
    + + + + +
    + +
    +
    +
    +

    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).

    +
    + + 
    aw(object, ...)
+
+# S3 method for mkinfit
+aw(object, ...)
+
+# S3 method for mmkin
+aw(object, ...)
    + +

    Arguments

    + + + + + + + + + + +
    object

    An mmkin column object, containing two or more +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.

    ...

    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)
    #> [1] 1
    aw_sfo_dfop # SFO gets more weight as it has less parameters and a similar fit
    #> [1] 0.5970258 0.4029742
    f <- mmkin(c("SFO", "FOMC", "DFOP"), list("FOCUS D" = FOCUS_2006_D), cores = 1, quiet = TRUE)
+aw(f)
    #> [1] 0.4808722 0.1945539 0.3245740
    sum(aw(f))
    #> [1] 1
    aw(f[c("SFO", "DFOP")])
    #> [1] 0.5970258 0.4029742
    # }
+
    +
    + +
    + + +
    +
    +

    Developed by Johannes Ranke.

    +
    + +
    +

    Site built with pkgdown 1.4.1.

    +
    + +
    +
    + + + + + + + + diff --git a/docs/reference/index.html b/docs/reference/index.html index 4398469b..3d417267 100644 --- a/docs/reference/index.html +++ b/docs/reference/index.html @@ -233,6 +233,12 @@ more datasets

    with mkinfit

    + +

    aw()

    + +

    Calculate Akaike weights for model averaging

    + +

    CAKE_export()

    diff --git a/docs/sitemap.xml b/docs/sitemap.xml index 66b776b2..a8d6dfa4 100644 --- a/docs/sitemap.xml +++ b/docs/sitemap.xml @@ -54,6 +54,9 @@ https://pkgdown.jrwb.de/mkin/reference/add_err.html + + https://pkgdown.jrwb.de/mkin/reference/aw.html + https://pkgdown.jrwb.de/mkin/reference/confint.mkinfit.html diff --git a/man/aw.Rd b/man/aw.Rd new file mode 100644 index 00000000..f0994b94 --- /dev/null +++ b/man/aw.Rd @@ -0,0 +1,47 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/aw.R +\name{aw} +\alias{aw} +\alias{aw.mkinfit} +\alias{aw.mmkin} +\title{Calculate Akaike weights for model averaging} +\usage{ +aw(object, ...) + +\method{aw}{mkinfit}(object, ...) + +\method{aw}{mmkin}(object, ...) +} +\arguments{ +\item{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.} + +\item{\dots}{Not used in the method for mmkin column objects, +further mkinfit objects in the method for mkinfit objects.} +} +\description{ +Akaike weights are calculated based on the relative +expected Kullback-Leibler information as specified +by Burnham and Anderson (2004). +} +\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")]) +} +} +\references{ +Burnham KP and Anderson DR (2004) Multimodel + Inference: Understanding AIC and BIC in Model Selection + Sociological Methods & Research 33(2) 261-304 +} diff --git a/test.log b/test.log index bc6d26ae..c51d06b8 100644 --- a/test.log +++ b/test.log @@ -1,11 +1,12 @@ Loading mkin Testing mkin ✔ | OK F W S | Context +✔ | 5 | Calculation of Akaike weights ✔ | 2 | Export dataset for reading into CAKE -✔ | 10 | Confidence intervals and p-values [10.1 s] -✔ | 14 | Error model fitting [40.5 s] +✔ | 10 | Confidence intervals and p-values [9.7 s] +✔ | 14 | Error model fitting [36.9 s] ✔ | 4 | Calculation of FOCUS chi2 error levels [2.2 s] -✔ | 13 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [3.4 s] +✔ | 13 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [3.3 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 @@ -18,20 +19,20 @@ Testing mkin ✔ | 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.6 s] -✔ | 4 | Fitting the SFORB model [1.8 s] +✔ | 2 | Complex test case from Schaefer et al. (2007) Piacenza paper [5.3 s] +✔ | 4 | Fitting the SFORB model [1.7 s] ✔ | 1 | Summaries of old mkinfit objects -✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [7.5 s] -✔ | 7 1 | Hypothesis tests [34.1 s] +✔ | 4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [7.2 s] +✔ | 7 1 | Hypothesis tests [32.3 s] ──────────────────────────────────────────────────────────────────────────────── -test_tests.R:59: skip: We can do a likelihood ratio test using an update specification +test_tests.R:60: 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: 116.9 s +Duration: 110.2 s -OK: 133 +OK: 138 Failed: 0 Warnings: 0 Skipped: 1 diff --git a/tests/testthat/FOCUS_2006_D.csf b/tests/testthat/FOCUS_2006_D.csf index 09940aa3..358b50e3 100644 --- a/tests/testthat/FOCUS_2006_D.csf +++ b/tests/testthat/FOCUS_2006_D.csf @@ -5,7 +5,7 @@ Description: MeasurementUnits: % AR TimeUnits: days Comments: Created using mkin::CAKE_export -Date: 2019-11-09 +Date: 2019-11-13 Optimiser: IRLS [Data] diff --git a/tests/testthat/test_aw.R b/tests/testthat/test_aw.R new file mode 100644 index 00000000..0a493893 --- /dev/null +++ b/tests/testthat/test_aw.R @@ -0,0 +1,12 @@ +context("Calculation of Akaike weights") + +test_that("Akaike weights sum to one", { + skip_on_cran() + aw_1 <- aw(fit_nw_1, fit_obs_1, fit_tc_1) + expect_error(aw(fit_nw_1, f_2_mkin), "same data") + expect_error(aw(fit_nw_1, 3), "mkinfit objects") + expect_equal(sum(aw_1), 1) + aw_2 <- aw(fits[c("SFO", "DFOP"), "FOCUS_D"]) + expect_equal(sum(aw_2), 1) + expect_error(aw(fits), "mmkin column object") +}) -- cgit v1.2.3