diff --git a/.Rbuildignore b/.Rbuildignore index 596581d9..fd782bab 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -33,4 +33,5 @@ ^vignettes/used_by*$ ^revdep$ ^data-raw$ -^paper$ \ No newline at end of file +^paper$ +^dev$ \ No newline at end of file diff --git a/.github/CONTRIBUTING.md b/.github/CONTRIBUTING.md new file mode 100644 index 00000000..62d32faf --- /dev/null +++ b/.github/CONTRIBUTING.md @@ -0,0 +1,58 @@ +# Contributing to priorsense + +This outlines how to propose a change to priorsense and is based on similar +instructions for tidyverse packages, including the contributing guidelines +generated by `usethis::use_tidy_contributing()`. + +## Fixing typos + +You can fix typos, spelling mistakes, or grammatical errors in the documentation +directly using the GitHub web interface, as long as the changes are made in the +_source_ file. This generally means you'll need to edit +[roxygen2 comments](https://roxygen2.r-lib.org/articles/roxygen2.html) in an `.R`, +not a `.Rd` file. You can find the `.R` file that generates the `.Rd` by reading +the comment in the first line. + +## Bigger changes + +If you want to make a bigger change, it's a good idea to first file an issue and +make sure someone from the team agrees that it’s needed. If you’ve found a bug, +please file an issue that illustrates the bug with a minimal reproducible +example (see e.g. the [tidyverse reprex instructions](https://www.tidyverse.org/help/#reprex)). +The tidyverse guide on [how to create a great issue](https://code-review.tidyverse.org/issues/) +has more advice. + +### Pull request process + +If you are new to creating pull requests here are some tips. Using the functions +from the `usethis` package is not required but can be helpful if this process is +new to you. + +* Fork the package and clone onto your computer. If you haven't done this before, we recommend using `usethis::create_from_github("stan-dev/priorsense", fork = TRUE)`. + +* Install all development dependencies with `devtools::install_dev_deps()`, and then make sure the package passes R CMD check by running `devtools::check()`. + If R CMD check doesn't pass cleanly, it's a good idea to ask for help before continuing. +* Create a Git branch for your pull request (PR). We recommend using `usethis::pr_init("brief-description-of-change")`. + +* Make your changes, commit to git, and then create a PR by running `usethis::pr_push()`, and following the prompts in your browser. + The title of your PR should briefly describe the change. + The body of your PR should contain `Fixes #issue-number`. + +* For user-facing changes, add a bullet to the top of `NEWS.md` (i.e. just below the first header). Follow the style already used in `NEWS.md`. + +### Code style + +* New code should attempt to follow the style used in the package. When in doubt follow the tidyverse [style guide](https://style.tidyverse.org). + +* We use [roxygen2](https://cran.r-project.org/package=roxygen2), with [Markdown syntax](https://cran.r-project.org/web/packages/roxygen2/vignettes/rd-formatting.html), for documentation. + +* We use [testthat](https://cran.r-project.org/package=testthat) for unit tests. + Contributions with test cases included are easier to accept. + +## Code of Conduct + +Please note that the priorsense project follows the Stan project's +[Code of Conduct](https://discourse.mc-stan.org/t/announcing-our-new-stan-code-of-conduct/23764). +By contributing to this project you agree to abide by its terms. + +All contributions must follow the [Stan AI Contribution Policy](https://github.com/stan-dev/stan/wiki/AI-Contribution-Policy). diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index e62ca74f..1afc2400 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -1,4 +1,4 @@ -# Workflow derived from https://github.com/r-lib/actions/tree/master/examples +# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples # Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help on: push: @@ -6,7 +6,9 @@ on: pull_request: branches: [main, development] -name: R-CMD-check +name: R-CMD-check.yaml + +permissions: read-all jobs: R-CMD-check: @@ -18,12 +20,9 @@ jobs: fail-fast: false matrix: config: - #- {os: macOS-latest, r: 'devel', rtools: ''} - #- {os: macOS-latest, r: 'release', rtools: ''} - #- {os: windows-latest, r: 'devel', rtools: '45'} - #- {os: windows-latest, r: 'release', rtools: '45'} - - {os: ubuntu-22.04, r: 'devel', rtools: ''} - - {os: ubuntu-22.04, r: 'release', rtools: ''} + - {os: ubuntu-latest, r: 'devel', http-user-agent: 'release'} + - {os: ubuntu-latest, r: 'release'} + - {os: ubuntu-latest, r: 'oldrel-1'} env: R_REMOTES_NO_ERRORS_FROM_WARNINGS: true @@ -40,12 +39,13 @@ jobs: with: token: ${{ secrets.GITHUB_TOKEN }} - - uses: actions/checkout@v4 + - uses: actions/checkout@v7 - uses: r-lib/actions/setup-r@v2 with: r-version: ${{ matrix.config.r }} rtools-version: ${{ matrix.config.rtools }} + http-user-agent: ${{ matrix.config.http-user-agent }} use-public-rspm: true - uses: r-lib/actions/setup-pandoc@v2 with: @@ -85,6 +85,6 @@ jobs: - uses: r-lib/actions/check-r-package@v2 with: - args: 'c("--no-manual", "--as-cran")' - error-on: '"warning"' - check-dir: '"check"' + upload-snapshots: true + args: 'c("--no-manual")' + build_args: 'c("--no-manual")' diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml index 1057d11c..deea88f7 100644 --- a/.github/workflows/pkgdown.yaml +++ b/.github/workflows/pkgdown.yaml @@ -29,50 +29,40 @@ jobs: - uses: r-lib/actions/setup-pandoc@v2 - uses: quarto-dev/quarto-actions/setup@v2 - - - uses: r-lib/actions/setup-r@v2 + + - name: Set up R + uses: r-lib/actions/setup-r@v2 with: use-public-rspm: true - - - uses: r-lib/actions/setup-r-dependencies@v2 + extra-repositories: | + https://stan-dev.r-universe.dev + https://topipa.r-universe.dev + + - name: Install pkgdown dependencies + uses: r-lib/actions/setup-r-dependencies@v2 with: - extra-packages: | + cache-version: 4 + packages: | + deps::. + any::sessioninfo any::pkgdown - local::. - XML - stan-dev/cmdstanr - paul-buerkner/brms - topipa/iwmm - rcmdcheck - rstan - checkmate - jsonlite - nimble - posterior - processx - R6 + needs: website + dependencies: '"all"' + extra-packages: | BH - R2jags RcppEigen StanHeaders RcppParallel - withr - testthat - quarto - any::XML + stan-dev/cmdstanr + topipa/iwmm + stan-dev/pkgdown-config any::textshaping - needs: website - - - - name: Build site + + - name: Build pkgdown site run: | - withr::with_envvar( - c("NOT_CRAN" = "true"), # this should already be set by setup-r@v2? keeping because vignettes don't build otherwise - pkgdown::build_site_github_pages( - lazy = FALSE, # change to TRUE if runner times out. - run_dont_run = TRUE, - new_process = TRUE - ) + pkgdown::build_site_github_pages( + install = TRUE, + new_process = FALSE ) shell: Rscript {0} diff --git a/.gitignore b/.gitignore index a5d43d49..ad6361b4 100644 --- a/.gitignore +++ b/.gitignore @@ -8,6 +8,7 @@ priorsense.Rcheck .DS_Store docs /doc/ +/dev/ /Meta/ /pkgdown/ /..Rcheck/ diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md deleted file mode 100644 index 7e1fd129..00000000 --- a/CONTRIBUTING.md +++ /dev/null @@ -1,6 +0,0 @@ -# Contributing - -Contributions are welcome! If you find an bug or have an idea for a feature, open an issue. If you are able to fix an issue, fork the repository and make a pull request to the `development` branch. - -## Lifecycle -priorsense is in a stable state of development, with some degree of active subsequent development as envisioned by the primary authors and in response to user feedback. diff --git a/DESCRIPTION b/DESCRIPTION index eb160572..4d9d7b34 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: priorsense Title: Prior Diagnostics and Sensitivity Analysis -Version: 1.3.1 -Authors@R: c(person("Noa", "Kallioinen", email = "noa.kallioinen@helsinki.fi", role = c("aut", "cre", "cph"), comment = c(ORCID = "0000-0003-1586-8382")), +Version: 1.4.0 +Authors@R: c(person("Noa", "Kallioinen", email = "noa.kallioinen@helsinki.fi", role = c("aut", "cre"), comment = c(ORCID = "0000-0003-1586-8382")), person("Topi", "Paananen", role = c("aut"), comment = c(ORCID = "0000-0002-6542-407X")), person("Paul-Christian", "Bürkner", role = c("aut"), comment = c(ORCID = "0000-0001-5765-8995")), person("Aki", "Vehtari", role = c("aut"), comment = c(ORCID = "0000-0003-2164-9469")), @@ -14,9 +14,9 @@ License: GPL (>= 3) Encoding: UTF-8 LazyData: true Roxygen: list(markdown = TRUE, roclets = c ("namespace", "rd", "srr::srr_stats_roclet")) -RoxygenNote: 8.0.0 Imports: checkmate (>= 2.3.4), + cli, ggdist (>= 3.3.3), ggh4x (>= 0.3.1), ggplot2 (>= 4.0.3), @@ -29,18 +29,18 @@ Imports: tibble (>= 3.3.1), utils Suggests: - bayesplot (>= 1.15.0), + bayesplot (>= 1.16.0), brms (>= 2.23.0), cmdstanr (>= 0.8.1), iwmm (>= 0.0.1), - nimble (>= 1.4.2), + nimble (>= 1.4.3), philentropy (>= 0.10.0), - quarto (>= 1.4.4), + quarto (>= 1.5.1), R2jags (>= 0.8), - rstan (>= 2.32.6), - testthat (>= 3.0.0), + rstan (>= 2.32.7), + testthat (>= 3.3.2), transport (>= 0.15), - vdiffr (>= 1.0.8) + vdiffr (>= 1.0.9) Config/testthat/edition: 3 Depends: R (>= 4.1.0) @@ -48,6 +48,7 @@ VignetteBuilder: quarto Additional_repositories: https://topipa.r-universe.dev, https://stan-dev.r-universe.dev -URL: https://n-kall.github.io/priorsense/, https://github.com/n-kall/priorsense -BugReports: https://github.com/n-kall/priorsense/issues +URL: https://mc-stan.org/priorsense/, https://mc-stan.org/priorsense +BugReports: https://github.com/stan-dev/priorsense/issues Config/Needs/website: quarto +Config/roxygen2/version: 8.1.0 diff --git a/NEWS.md b/NEWS.md index 3701c435..b5c91b94 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,9 @@ +priorsense 1.3.1.9000 +--- ++ `priorsense` is now an official Stan package, and links to pages and repositories have been updated ++ Plot help text now prints to the console rather than on the plot ++ `powerscale_plot_quantities()` now highlights the line segment closer to 1 when k-hat is high + priorsense 1.3.1 --- + Add new JOSS paper as primary citation for the package. diff --git a/R/log_prior_draws.R b/R/log_prior_draws.R index d3dede04..74479564 100644 --- a/R/log_prior_draws.R +++ b/R/log_prior_draws.R @@ -20,93 +20,93 @@ ##' ##' @export log_prior_draws <- function(x, ...) { - UseMethod("log_prior_draws") + UseMethod("log_prior_draws") } ##' @rdname log_prior_draws ##' @export log_prior_draws.stanfit <- function( - x, - joint = FALSE, - log_prior_name = "lprior", - ... + x, + joint = FALSE, + log_prior_name = "lprior", + ... ) { - if (!inherits(x, "stanfit")) { - stop("Not a stanfit object.", call. = FALSE) - } - if (x@mode != 0) { - stop("Stan model does not contain posterior draws.", call. = FALSE) - } - if (!requireNamespace("rstan", quietly = TRUE)) { - stop("Please load the 'rstan' package.", call. = FALSE) - } + if (!inherits(x, "stanfit")) { + stop("Not a stanfit object.", call. = FALSE) + } + if (x@mode != 0) { + stop("Stan model does not contain posterior draws.", call. = FALSE) + } + if (!requireNamespace("rstan", quietly = TRUE)) { + stop("Please load the 'rstan' package.", call. = FALSE) + } - checkmate::assert_logical(joint, len = 1) - checkmate::assert_character(log_prior_name, len = 1) + checkmate::assert_logical(joint, len = 1) + checkmate::assert_character(log_prior_name) - log_prior <- posterior::subset_draws( - posterior::as_draws_array(x), - variable = paste0("^", log_prior_name), - regex = TRUE - ) + log_prior <- posterior::subset_draws( + posterior::as_draws_array(x), + variable = paste0("^", log_prior_name), + regex = TRUE + ) - if (joint) { - log_prior <- rowsums_draws(log_prior) - posterior::variables(log_prior) <- log_prior_name - } + if (joint) { + log_prior <- rowsums_draws(log_prior) + posterior::variables(log_prior) <- log_prior_name + } - return(log_prior) + return(log_prior) } ##' @rdname log_prior_draws ##' @export log_prior_draws.CmdStanFit <- function( - x, - joint = FALSE, - log_prior_name = "lprior", - ... + x, + joint = FALSE, + log_prior_name = "lprior", + ... ) { - checkmate::assert_logical(joint, len = 1) - checkmate::assert_character(log_prior_name, len = 1) + checkmate::assert_logical(joint, len = 1) + checkmate::assert_character(log_prior_name, len = 1) - all_draws <- x$draws() + all_draws <- x$draws() - log_prior <- posterior::subset_draws( - all_draws, - variable = paste0("^", log_prior_name), - regex = TRUE - ) + log_prior <- posterior::subset_draws( + all_draws, + variable = paste0("^", log_prior_name), + regex = TRUE + ) - if (joint) { - log_prior <- rowsums_draws(log_prior) - posterior::variables(log_prior) <- log_prior_name - } + if (joint) { + log_prior <- rowsums_draws(log_prior) + posterior::variables(log_prior) <- log_prior_name + } - return(log_prior) + return(log_prior) } ##' @rdname log_prior_draws ##' @export log_prior_draws.draws <- function( - x, - joint = FALSE, - log_prior_name = "lprior", - ... + x, + joint = FALSE, + log_prior_name = "lprior", + ... ) { - checkmate::assert_logical(joint, len = 1) - checkmate::assert_character(log_prior_name, len = 1) + checkmate::assert_logical(joint, len = 1) + checkmate::assert_character(log_prior_name) - log_prior <- posterior::subset_draws( - x, - variable = paste0("^", log_prior_name), - regex = TRUE - ) + log_prior <- posterior::subset_draws( + x, + variable = paste0("^", log_prior_name), + regex = TRUE + ) - if (joint) { - log_prior <- rowsums_draws(log_prior) - posterior::variables(log_prior) <- log_prior_name - } + if (joint) { + log_prior <- rowsums_draws(log_prior) + posterior::variables(log_prior) <- log_prior_name + } - return(log_prior) + return(log_prior) } diff --git a/R/plots.R b/R/plots.R index 1049983f..e4c89d61 100644 --- a/R/plots.R +++ b/R/plots.R @@ -42,112 +42,114 @@ NULL ##' @keywords internal ##' @noRd prepare_plot_data <- function(x, variable, resample, ...) { - base_draws <- posterior::merge_chains(x$base_draws) - - if (!resample && !(x$resampled)) { - base_draws <- posterior::weight_draws( - x = base_draws, - weights = rep( - 1 / posterior::ndraws(base_draws), - posterior::ndraws(base_draws) - ) - ) - } + base_draws <- posterior::merge_chains(x$base_draws) + + if (!resample && !(x$resampled)) { + base_draws <- posterior::weight_draws( + x = base_draws, + weights = rep( + 1 / posterior::ndraws(base_draws), + posterior::ndraws(base_draws) + ) + ) + } - likelihood_draws <- list() - prior_draws <- list() - base_draws_prior <- data.frame() - base_draws_lik <- data.frame() + likelihood_draws <- list() + prior_draws <- list() + base_draws_prior <- data.frame() + base_draws_lik <- data.frame() - if (!(is.null(x$prior_scaled))) { - prior_scaled <- x$prior_scaled$draws_sequence + if (!(is.null(x$prior_scaled))) { + prior_scaled <- x$prior_scaled$draws_sequence - for (i in seq_along(prior_scaled)) { - prior_draws[[i]] <- posterior::merge_chains(prior_scaled[[i]]) + for (i in seq_along(prior_scaled)) { + prior_draws[[i]] <- posterior::merge_chains(prior_scaled[[i]]) - if (resample && !x$resampled) { - prior_draws[[i]] <- posterior::resample_draws(prior_draws[[i]]) - } + if (resample && !x$resampled) { + prior_draws[[i]] <- posterior::resample_draws(prior_draws[[i]]) + } - prior_ps_details <- get_powerscaling_details(prior_scaled[[i]]) + prior_ps_details <- get_powerscaling_details(prior_scaled[[i]]) - prior_draws[[i]][[".powerscale_alpha"]] <- prior_ps_details$alpha - prior_draws[[i]]$component <- "prior" - prior_draws[[i]]$pareto_k <- prior_ps_details$diagnostics$khat - prior_draws[[i]]$pareto_k_threshold <- - prior_ps_details$diagnostics$khat_threshold + prior_draws[[i]][[".powerscale_alpha"]] <- prior_ps_details$alpha + prior_draws[[i]]$component <- "prior" + prior_draws[[i]]$pareto_k <- prior_ps_details$diagnostics$khat + prior_draws[[i]]$pareto_k_threshold <- + prior_ps_details$diagnostics$khat_threshold + } + + base_draws_prior <- base_draws + base_draws_prior[[".powerscale_alpha"]] <- 1 + base_draws_prior$component <- "prior" + base_draws_prior$pareto_k <- -Inf + base_draws_prior$pareto_k_threshold <- Inf } - base_draws_prior <- base_draws - base_draws_prior[[".powerscale_alpha"]] <- 1 - base_draws_prior$component <- "prior" - base_draws_prior$pareto_k <- -Inf - base_draws_prior$pareto_k_threshold <- Inf - } + if (!(is.null(x$likelihood_scaled))) { + likelihood_scaled <- x$likelihood_scaled$draws_sequence - if (!(is.null(x$likelihood_scaled))) { - likelihood_scaled <- x$likelihood_scaled$draws_sequence + for (i in seq_along(likelihood_scaled)) { + likelihood_draws[[i]] <- posterior::merge_chains( + likelihood_scaled[[i]] + ) - for (i in seq_along(likelihood_scaled)) { - likelihood_draws[[i]] <- posterior::merge_chains( - likelihood_scaled[[i]] - ) + if (resample && !x$resampled) { + likelihood_draws[[i]] <- posterior::resample_draws( + likelihood_draws[[i]] + ) + } - if (resample && !x$resampled) { - likelihood_draws[[i]] <- posterior::resample_draws( - likelihood_draws[[i]] - ) - } + likelihood_ps_details <- get_powerscaling_details( + likelihood_scaled[[i]] + ) - likelihood_ps_details <- get_powerscaling_details( - likelihood_scaled[[i]] - ) + likelihood_draws[[i]][[".powerscale_alpha"]] <- + likelihood_ps_details$alpha - likelihood_draws[[i]][[".powerscale_alpha"]] <- - likelihood_ps_details$alpha + likelihood_draws[[i]]$component <- "likelihood" - likelihood_draws[[i]]$component <- "likelihood" + likelihood_draws[[ + i + ]]$pareto_k <- likelihood_ps_details$diagnostics$khat - likelihood_draws[[i]]$pareto_k <- likelihood_ps_details$diagnostics$khat + likelihood_draws[[i]]$pareto_k_threshold <- + likelihood_ps_details$diagnostics$khat_threshold + } - likelihood_draws[[i]]$pareto_k_threshold <- - likelihood_ps_details$diagnostics$khat_threshold + base_draws_lik <- base_draws + base_draws_lik[[".powerscale_alpha"]] <- 1 + base_draws_lik$component <- "likelihood" + base_draws_lik$pareto_k <- -Inf + base_draws_lik$pareto_k_threshold <- Inf } - base_draws_lik <- base_draws - base_draws_lik[[".powerscale_alpha"]] <- 1 - base_draws_lik$component <- "likelihood" - base_draws_lik$pareto_k <- -Inf - base_draws_lik$pareto_k_threshold <- Inf - } - - d <- rbind( - do.call("rbind", prior_draws), - do.call("rbind", likelihood_draws), - base_draws_lik, - base_draws_prior - ) - - d$pareto_k_value <- ifelse(d$pareto_k > d$pareto_k_threshold, "High", "OK") - - d$pareto_k_value <- factor( - d$pareto_k_value, - levels = c("OK", "High") - ) - - d$component <- factor(d$component, levels = c("prior", "likelihood")) - - # prepare for plotting - d <- stats::reshape( - data = as.data.frame(d), - varying = variable, - direction = "long", - times = variable, - v.names = "value", - timevar = "variable" - ) - - return(d) + d <- rbind( + do.call("rbind", prior_draws), + do.call("rbind", likelihood_draws), + base_draws_lik, + base_draws_prior + ) + + d$pareto_k_value <- ifelse(d$pareto_k > d$pareto_k_threshold, "High", "OK") + + d$pareto_k_value <- factor( + d$pareto_k_value, + levels = c("OK", "High") + ) + + d$component <- factor(d$component, levels = c("prior", "likelihood")) + + # prepare for plotting + d <- stats::reshape( + data = as.data.frame(d), + varying = variable, + direction = "long", + times = variable, + v.names = "value", + timevar = "variable" + ) + + return(d) } ##' prepare plot ##' @param d data frame of data from plotting @@ -159,816 +161,901 @@ prepare_plot_data <- function(x, variable, resample, ...) { ##' @keywords internal ##' @noRd prepare_plot <- function(d, resample, variable, colors, ...) { - if (resample) { - p <- ggplot2::ggplot( - data = d, - ggplot2::aes( - x = .data$value, - group = .data[[".powerscale_alpha"]], - color = .data[[".powerscale_alpha"]], - linetype = .data$pareto_k_value - ) - ) - } else { - p <- ggplot2::ggplot( - data = d, - ggplot2::aes( - x = .data$value, - weight = exp(.data$.log_weight), - group = .data[[".powerscale_alpha"]], - color = .data[[".powerscale_alpha"]], - linetype = .data$pareto_k_value - ) - ) - } - - p <- p + - ggplot2::scale_linetype_manual( - values = c("solid", "dashed"), - drop = TRUE, - name = "Pareto k" - ) + - ggplot2::scale_color_gradientn( - name = "Power-scaling alpha", - colours = colors[1:3], - trans = "log", - limits = c( - min(d[[".powerscale_alpha"]]) - 0.01, - max(d[[".powerscale_alpha"]]) + 0.01 - ), - breaks = c( - min(d[[".powerscale_alpha"]]), - 1, - max(d[[".powerscale_alpha"]]) - ), - labels = c( - round(min(d[[".powerscale_alpha"]]), digits = 3), - "1", - round(max(d[[".powerscale_alpha"]]), digits = 3) - ) - ) + - ggplot2::scale_fill_gradientn( - colours = c(colors[1:3]), - trans = "log", - limits = c( - min(d[[".powerscale_alpha"]]) - 0.01, - max(d[[".powerscale_alpha"]]) + 0.01 - ), - breaks = c( - min(d[[".powerscale_alpha"]]), - 1, - max(d[[".powerscale_alpha"]]) - ), - labels = c( - round(min(d[[".powerscale_alpha"]]), digits = 3), - "1", - round(max(d[[".powerscale_alpha"]]), digits = 3) - ) - ) + if (resample) { + p <- ggplot2::ggplot( + data = d, + ggplot2::aes( + x = .data$value, + group = .data[[".powerscale_alpha"]], + color = .data[[".powerscale_alpha"]], + linetype = .data$pareto_k_value + ) + ) + } else { + p <- ggplot2::ggplot( + data = d, + ggplot2::aes( + x = .data$value, + weight = exp(.data$.log_weight), + group = .data[[".powerscale_alpha"]], + color = .data[[".powerscale_alpha"]], + linetype = .data$pareto_k_value + ) + ) + } - if (length(unique(d[[".powerscale_alpha"]])) == 3) { p <- p + - ggplot2::guides( - color = ggplot2::guide_legend( - override.aes = ggplot2::aes(linetype = "solid") + ggplot2::scale_linetype_manual( + values = c("solid", "dashed"), + drop = TRUE + ) + + ggplot2::scale_color_gradientn( + name = "Power-scaling alpha", + colours = colors[1:3], + trans = "log", + limits = c( + min(d[[".powerscale_alpha"]]) - 0.01, + max(d[[".powerscale_alpha"]]) + 0.01 + ), + breaks = c( + min(d[[".powerscale_alpha"]]), + 1, + max(d[[".powerscale_alpha"]]) + ), + labels = c( + round(min(d[[".powerscale_alpha"]]), digits = 3), + "1", + round(max(d[[".powerscale_alpha"]]), digits = 3) + ) + ) + + ggplot2::scale_fill_gradientn( + name = "Power-scaling alpha", + colours = c(colors[1:3]), + trans = "log", + limits = c( + min(d[[".powerscale_alpha"]]) - 0.01, + max(d[[".powerscale_alpha"]]) + 0.01 + ), + breaks = c( + min(d[[".powerscale_alpha"]]), + 1, + max(d[[".powerscale_alpha"]]) + ), + labels = c( + round(min(d[[".powerscale_alpha"]]), digits = 3), + "1", + round(max(d[[".powerscale_alpha"]]), digits = 3) + ) ) - ) - } - if (!(any(d$pareto_k_value == "High"))) { - p <- p + - ggplot2::guides( - linetype = "none" - ) - } + if (length(unique(d[[".powerscale_alpha"]])) == 3) { + p <- p + + ggplot2::guides( + color = ggplot2::guide_legend( + title = "Power-scaling alpha", + override.aes = ggplot2::aes(linetype = "solid"), + order = 5 + ) + ) + } else { + p <- p + + ggplot2::guides( + fill = ggplot2::guide_legend( + title = "Power-scaling alpha", + order = 5 + ) + ) + } - return(p) + if (any(d$pareto_k_value == "High")) { + p <- p + + ggplot2::guides( + linetype = ggplot2::guide_legend( + title = "Pareto k", + order = 10 + ) + ) + } else { + p <- p + + ggplot2::guides(linetype = "none") + } + return(p) } ##' @rdname powerscale-plots ##' @export powerscale_plot_dens <- function(x, ...) { - UseMethod("powerscale_plot_dens") + UseMethod("powerscale_plot_dens") } ##' @export powerscale_plot_dens.default <- - function( - x, - variable = NULL, - variables = NULL, - length = 3, - resample = FALSE, - intervals = c(0.5, 0.8, 0.95), - trim = NULL, - facet_rows = "component", - help_text = getOption("priorsense.plot_help_text", TRUE), - colors = NULL, - colours = NULL, - variables_per_page = 6, - ... - ) { - ps <- powerscale_sequence(x, length = length, ...) - powerscale_plot_dens( - ps, - variable = variable, - variables = variables, - length = length, - resample = resample, - intervals = intervals, - trim = trim, - facet_rows = facet_rows, - help_text = help_text, - colors = colors, - colours = colours, - variables_per_page = variables_per_page - ) - } - -draw_key_path2 <- function(data, params, size) { - grid::segmentsGrob( - x0 = 0.1, - x1 = 0.9, - y0 = 0.5, - y1 = 0.5, - gp = grid::gpar(col = data$colour) - ) -} + function( + x, + variable = NULL, + variables = NULL, + length = 3, + resample = FALSE, + intervals = c(0.5, 0.8, 0.95), + trim = NULL, + facet_rows = "component", + help_text = getOption("priorsense.plot_help_text", TRUE), + colors = NULL, + colours = NULL, + variables_per_page = 6, + ... + ) { + ps <- powerscale_sequence(x, length = length, ...) + powerscale_plot_dens( + ps, + variable = variable, + variables = variables, + length = length, + resample = resample, + intervals = intervals, + trim = trim, + facet_rows = facet_rows, + help_text = help_text, + colors = colors, + colours = colours, + variables_per_page = variables_per_page + ) + } ##' @export powerscale_plot_dens.powerscaled_sequence <- - function( - x, - variable = NULL, - variables = NULL, - resample = FALSE, - intervals = c(0.5, 0.8, 0.95), - trim = NULL, - facet_rows = "component", - help_text = getOption("priorsense.plot_help_text", TRUE), - colors = NULL, - colours = NULL, - variables_per_page = getOption( - "priorsense.plot_variables_per_page", - 6 - ), - ... - ) { - # input checks - if (!is.null(variable) && !is.null(variables)) { - checkmate::assert( - if (identical(variable, variables)) { - TRUE - } else { - "must be identical if both provided" - }, - .var.name = "`variable` and `variables`" - ) - } - if (is.null(variable)) { - variable <- variables - } - - if (!is.null(colors) && !is.null(colours)) { - checkmate::assert( - if (identical(colors, colours)) { - TRUE - } else { - "must be identical if both provided" - }, - .var.name = "`colors` and `colours`" - ) - } - if (is.null(colors)) { - colors <- colours - } - - checkmate::assert_character(variable, null.ok = TRUE) - checkmate::assert_logical(resample, len = 1) - checkmate::assert_logical(help_text, len = 1) - checkmate::assert_character(colors, len = 3, null.ok = TRUE) - checkmate::assert_numeric(intervals, null.ok = TRUE) - checkmate::assert_choice(facet_rows, c("component", "variable")) - checkmate::assert_number(variables_per_page, lower = 1, null.ok = TRUE) - - if (is.null(colors)) { - colors <- default_priorsense_colors()[1:3] - } - - if (is.null(variable)) { - variable <- posterior::variables(x$base_draws) - } else { - variable <- posterior::variables( - posterior::subset_draws(x$base_draws, variable = variable) - ) - } - - nvars <- length(variable) + function( + x, + variable = NULL, + variables = NULL, + resample = FALSE, + intervals = c(0.5, 0.8, 0.95), + trim = NULL, + facet_rows = "component", + help_text = getOption("priorsense.plot_help_text", TRUE), + colors = NULL, + colours = NULL, + variables_per_page = getOption( + "priorsense.plot_variables_per_page", + 6 + ), + ... + ) { + # input checks + if (!is.null(variable) && !is.null(variables)) { + checkmate::assert( + if (identical(variable, variables)) { + TRUE + } else { + "must be identical if both provided" + }, + .var.name = "`variable` and `variables`" + ) + } + if (is.null(variable)) { + variable <- variables + } - if (is.null(variables_per_page) || is.infinite(variables_per_page)) { - variables_per_page <- nvars - } + if (!is.null(colors) && !is.null(colours)) { + checkmate::assert( + if (identical(colors, colours)) { + TRUE + } else { + "must be identical if both provided" + }, + .var.name = "`colors` and `colours`" + ) + } + if (is.null(colors)) { + colors <- colours + } - variables_per_page <- floor(variables_per_page) + checkmate::assert_character(variable, null.ok = TRUE) + checkmate::assert_logical(resample, len = 1) + checkmate::assert_logical(help_text, len = 1) + checkmate::assert_character(colors, len = 3, null.ok = TRUE) + checkmate::assert_numeric(intervals, null.ok = TRUE) + checkmate::assert_choice(facet_rows, c("component", "variable")) + checkmate::assert_number(variables_per_page, lower = 1, null.ok = TRUE) - n_plots <- ceiling(nvars / variables_per_page) - plots <- vector(mode = "list", length = n_plots) + if (is.null(colors)) { + colors <- default_priorsense_colors()[1:3] + } - d <- prepare_plot_data(x, variable = variable, resample = resample, ...) + if (is.null(variable)) { + variable <- posterior::variables(x$base_draws) + } else { + variable <- posterior::variables( + posterior::subset_draws(x$base_draws, variable = variable) + ) + } - interval_positions <- data.frame( - ".powerscale_alpha" = unique(d[[".powerscale_alpha"]]), - interval_y = as.numeric(as.factor(unique(d[[".powerscale_alpha"]]))) - ) + nvars <- length(variable) - interval_positions$interval_y <- (0.5 - (interval_positions$interval_y)) / - (6 * nrow(interval_positions)) + if (is.null(variables_per_page) || is.infinite(variables_per_page)) { + variables_per_page <- nvars + } - d <- merge(d, interval_positions) + variables_per_page <- floor(variables_per_page) - n_components <- length(unique(d$component)) + n_plots <- ceiling(nvars / variables_per_page) + plots <- vector(mode = "list", length = n_plots) - for (i in seq_len(n_plots)) { - sub <- ((i - 1) * - variables_per_page + - 1):min(i * variables_per_page, nvars) - sub_variable <- variable[sub] - - if (resample || x$resample) { - resample <- TRUE - } - - dsub <- d[d$variable %in% sub_variable, ] - - plot <- prepare_plot(dsub, resample = resample, colors = colors, ...) + - ggplot2::ylab("Density") - - # here we have to draw 2 stat slabs (one with alpha 0 and black fill, - # one with fill as NA) to get the legend correct see - # https://github.com/mjskay/ggdist/issues/134 - plot <- plot + - ggdist::stat_slab( - fill = NA, - linewidth = 0.5, - trim = FALSE, - normalize = "xy", - key_glyph = draw_key_path2 - ) + - ggplot2::xlab(NULL) + - ggplot2::ylab(NULL) - - if (!is.null(intervals)) { - plot <- plot + - ggdist::stat_pointinterval( - ggplot2::aes(y = .data$interval_y), - .width = intervals, - fill = NA, - trim = FALSE, - normalize = "xy", - show.legend = FALSE - ) - } - - if (facet_rows == "component") { - plot <- plot + - ggh4x::facet_grid2( - rows = ggplot2::vars(.data$component), - cols = ggplot2::vars(.data$variable), - labeller = ggplot2::labeller( - component = c( - likelihood = "Likelihood\npower-scaling", - prior = "Prior\npower-scaling" - ) - ), - independent = "all", - scales = "free", - switch = "y" - ) - } else { - plot <- plot + - ggh4x::facet_grid2( - rows = ggplot2::vars(.data$variable), - cols = ggplot2::vars(.data$component), - labeller = ggplot2::labeller( - component = c( - likelihood = "Likelihood\npower-scaling", - prior = "Prior\npower-scaling" - ) - ), - independent = "all", - scales = "free", - switch = "y" - ) - } - - if (help_text) { - plot <- plot + - ggplot2::ggtitle( - label = "Power-scaling sensitivity", - subtitle = paste0( - "Posterior density estimates depending on ", - "amount of power-scaling (alpha).\n", - "Overlapping lines indicate low sensitivity.\n", - "Wider gaps between lines indicate greater sensitivity.\n", - "Estimates with high Pareto k (dashed lines) may be inaccurate." - ) - ) - } - - # additional theming - plot <- plot + - ggplot2::theme( - axis.line.y = ggplot2::element_blank(), - axis.text.y = ggplot2::element_blank(), - axis.ticks.y = ggplot2::element_blank() - ) + d <- prepare_plot_data(x, variable = variable, resample = resample, ...) - if (facet_rows == "component") { - plot <- plot + - ggplot2::theme(legend.position = "bottom") - } - - if (!is.null(trim)) { - position_scales <- lapply( - variable, - FUN = function(.x, prob) { - limits <- posterior::quantile2( - x$base_draws[[.x]], - probs = c((1 - prob) / 2, prob + (1 - prob) / 2) - ) - return(ggplot2::scale_x_continuous(limits = limits)) - }, - prob = trim + interval_positions <- data.frame( + ".powerscale_alpha" = unique(d[[".powerscale_alpha"]]), + interval_y = as.numeric(as.factor(unique(d[[".powerscale_alpha"]]))) ) - if (facet_rows == "component") { - plot <- plot + - ggh4x::facetted_pos_scales( - x = rep( - position_scales, - times = 2 - ) - ) - } else { - plot <- plot + - ggh4x::facetted_pos_scales( - x = rep( - position_scales, - each = 2 - ) - ) + interval_positions$interval_y <- (0.5 - + (interval_positions$interval_y)) / + (6 * nrow(interval_positions)) + + d <- merge(d, interval_positions) + + n_components <- length(unique(d$component)) + + for (i in seq_len(n_plots)) { + sub <- ((i - 1) * + variables_per_page + + 1):min(i * variables_per_page, nvars) + sub_variable <- variable[sub] + + if (resample || x$resample) { + resample <- TRUE + } + + dsub <- d[d$variable %in% sub_variable, ] + + plot <- prepare_plot( + dsub, + resample = resample, + colors = colors, + ... + ) + + ggplot2::ylab("Density") + + # here we have to draw 2 stat slabs (one with alpha 0 and black fill, + # one with fill as NA) to get the legend correct see + # https://github.com/mjskay/ggdist/issues/134 + plot <- plot + + ggdist::stat_slab( + fill = NA, + linewidth = 0.5, + trim = FALSE, + normalize = "xy", + key_glyph = draw_key_path2 + ) + + ggplot2::xlab(NULL) + + ggplot2::ylab(NULL) + + if (!is.null(intervals)) { + plot <- plot + + ggdist::stat_pointinterval( + ggplot2::aes(y = .data$interval_y), + .width = intervals, + fill = NA, + trim = FALSE, + normalize = "xy", + show.legend = FALSE + ) + } + + if (facet_rows == "component") { + plot <- plot + + ggh4x::facet_grid2( + rows = ggplot2::vars(.data$component), + cols = ggplot2::vars(.data$variable), + labeller = ggplot2::labeller( + component = c( + likelihood = "Likelihood\npower-scaling", + prior = "Prior\npower-scaling" + ) + ), + independent = "all", + scales = "free", + switch = "y" + ) + } else { + plot <- plot + + ggh4x::facet_grid2( + rows = ggplot2::vars(.data$variable), + cols = ggplot2::vars(.data$component), + labeller = ggplot2::labeller( + component = c( + likelihood = "Likelihood\npower-scaling", + prior = "Prior\npower-scaling" + ) + ), + independent = "all", + scales = "free", + switch = "y" + ) + } + + if (help_text) { + plot_help_text(type = "density") + } + + # additional theming + plot <- plot + + ggplot2::theme( + axis.line.y = ggplot2::element_blank(), + axis.text.y = ggplot2::element_blank(), + axis.ticks.y = ggplot2::element_blank() + ) + + if (facet_rows == "component") { + plot <- plot + + ggplot2::theme(legend.position = "bottom") + } + + if (!is.null(trim)) { + position_scales <- lapply( + variable, + FUN = function(.x, prob) { + limits <- posterior::quantile2( + x$base_draws[[.x]], + probs = c((1 - prob) / 2, prob + (1 - prob) / 2) + ) + return(ggplot2::scale_x_continuous(limits = limits)) + }, + prob = trim + ) + + if (facet_rows == "component") { + plot <- plot + + ggh4x::facetted_pos_scales( + x = rep( + position_scales, + times = 2 + ) + ) + } else { + plot <- plot + + ggh4x::facetted_pos_scales( + x = rep( + position_scales, + each = 2 + ) + ) + } + } + + plots[[i]] <- plot } - } - plots[[i]] <- plot - } + class(plots) <- c("priorsense_plot", class(plots)) - class(plots) <- c("priorsense_plot", class(plots)) + if (length(plots) == 1) { + plots <- plots[[1]] + } - if (length(plots) == 1) { - plots <- plots[[1]] + return(plots) } - return(plots) - } +draw_key_path2 <- function(data, params, size) { + grid::segmentsGrob( + x0 = 0.1, + x1 = 0.9, + y0 = 0.5, + y1 = 0.5, + gp = grid::gpar( + col = data$colour, + lty = data$linetype, + lwd = 2 * data$linewidth + ) + ) +} ##' @rdname powerscale-plots ##' @export powerscale_plot_ecdf <- function(x, ...) { - UseMethod("powerscale_plot_ecdf") + UseMethod("powerscale_plot_ecdf") } ##' @export powerscale_plot_ecdf.default <- - function( - x, - variable = NULL, - variables = NULL, - length = 3, - resample = FALSE, - facet_rows = "component", - help_text = getOption("priorsense.plot_help_text", TRUE), - colors = NULL, - colours = NULL, - variables_per_page = getOption( - "priorsense.plot_variables_per_page", - 6 - ), - ... - ) { - ps <- powerscale_sequence(x, length = length, ...) - powerscale_plot_ecdf( - ps, - variable = variable, - variables = variables, - resample = resample, - facet_rows = facet_rows, - help_text = help_text, - colors = colors, - colours = colours, - variables_per_page = variables_per_page - ) - } + function( + x, + variable = NULL, + variables = NULL, + length = 3, + resample = FALSE, + facet_rows = "component", + help_text = getOption("priorsense.plot_help_text", TRUE), + colors = NULL, + colours = NULL, + variables_per_page = getOption( + "priorsense.plot_variables_per_page", + 6 + ), + ... + ) { + ps <- powerscale_sequence(x, length = length, ...) + powerscale_plot_ecdf( + ps, + variable = variable, + variables = variables, + resample = resample, + facet_rows = facet_rows, + help_text = help_text, + colors = colors, + colours = colours, + variables_per_page = variables_per_page + ) + } ##' @rdname powerscale-plots ##' @export powerscale_plot_ecdf.powerscaled_sequence <- - function( - x, - variable = NULL, - variables = NULL, - resample = FALSE, - length = 3, - facet_rows = "component", - help_text = getOption("priorsense.plot_help_text", TRUE), - colors = NULL, - colours = NULL, - variables_per_page = getOption( - "priorsense.plot_variables_per_page", - 6 - ), - ... - ) { - # input checks - if (!is.null(variable) && !is.null(variables)) { - checkmate::assert( - if (identical(variable, variables)) { - TRUE - } else { - "must be identical if both provided" - }, - .var.name = "`variable` and `variables`" - ) - } - if (is.null(variable)) { - variable <- variables - } + function( + x, + variable = NULL, + variables = NULL, + resample = FALSE, + length = 3, + facet_rows = "component", + help_text = getOption("priorsense.plot_help_text", TRUE), + colors = NULL, + colours = NULL, + variables_per_page = getOption( + "priorsense.plot_variables_per_page", + 6 + ), + ... + ) { + # input checks + if (!is.null(variable) && !is.null(variables)) { + checkmate::assert( + if (identical(variable, variables)) { + TRUE + } else { + "must be identical if both provided" + }, + .var.name = "`variable` and `variables`" + ) + } + if (is.null(variable)) { + variable <- variables + } - if (!is.null(colors) && !is.null(colours)) { - checkmate::assert( - if (identical(colors, colours)) { - TRUE + if (!is.null(colors) && !is.null(colours)) { + checkmate::assert( + if (identical(colors, colours)) { + TRUE + } else { + "must be identical if both provided" + }, + .var.name = "`colors` and `colours`" + ) + } + if (is.null(colors)) { + colors <- colours + } + + checkmate::assert_character(variable, null.ok = TRUE) + checkmate::assert_logical(resample, len = 1) + checkmate::assert_logical(help_text, len = 1) + checkmate::assertCharacter(colors, len = 3, null.ok = TRUE) + checkmate::assert_choice(facet_rows, c("component", "variable")) + checkmate::assert_number(variables_per_page, lower = 1, null.ok = TRUE) + + if (is.null(colors)) { + colors <- default_priorsense_colors()[1:3] + } + + if (is.null(variable)) { + variable <- posterior::variables(x$base_draws) } else { - "must be identical if both provided" - }, - .var.name = "`colors` and `colours`" - ) - } - if (is.null(colors)) { - colors <- colours - } + variable <- posterior::variables( + posterior::subset_draws(x$base_draws, variable = variable) + ) + } - checkmate::assert_character(variable, null.ok = TRUE) - checkmate::assert_logical(resample, len = 1) - checkmate::assert_logical(help_text, len = 1) - checkmate::assertCharacter(colors, len = 3, null.ok = TRUE) - checkmate::assert_choice(facet_rows, c("component", "variable")) - checkmate::assert_number(variables_per_page, lower = 1, null.ok = TRUE) + d <- prepare_plot_data(x, variable = variable, resample = resample, ...) - if (is.null(colors)) { - colors <- default_priorsense_colors()[1:3] - } + n_components <- length(unique(d$component)) - if (is.null(variable)) { - variable <- posterior::variables(x$base_draws) - } else { - variable <- posterior::variables( - posterior::subset_draws(x$base_draws, variable = variable) - ) - } + if (resample || x$resample) { + resample <- TRUE + } - d <- prepare_plot_data(x, variable = variable, resample = resample, ...) + nvars <- length(variable) - n_components <- length(unique(d$component)) + if (is.null(variables_per_page) || is.infinite(variables_per_page)) { + variables_per_page <- nvars + } - if (resample || x$resample) { - resample <- TRUE - } + variables_per_page <- floor(variables_per_page) - nvars <- length(variable) + n_plots <- ceiling(nvars / variables_per_page) + plots <- vector(mode = "list", length = n_plots) - if (is.null(variables_per_page) || is.infinite(variables_per_page)) { - variables_per_page <- nvars - } + for (i in seq_len(n_plots)) { + sub <- ((i - 1) * variables_per_page + 1):min( + i * variables_per_page, + nvars + ) + sub_variable <- variable[sub] + + dsub <- d[d$variable %in% sub_variable, ] + + p <- prepare_plot(dsub, resample = resample, colors = colors, ...) + + ggplot2::ylab("ECDF") + + ggplot2::xlab(NULL) + + p <- p + + ggplot2::stat_ecdf(ggplot2::aes( + color = .data[[".powerscale_alpha"]] + )) + + if (facet_rows == "component") { + p <- p + + ggh4x::facet_grid2( + rows = ggplot2::vars(.data$component), + cols = ggplot2::vars(.data$variable), + labeller = ggplot2::labeller( + component = c( + likelihood = "Likelihood\npower-scaling", + prior = "Prior\npower-scaling" + ) + ), + scales = "free", + independent = "all", + switch = "y" + ) + } else { + p <- p + + ggh4x::facet_grid2( + rows = ggplot2::vars(.data$variable), + cols = ggplot2::vars(.data$component), + labeller = ggplot2::labeller( + component = c( + likelihood = "Likelihood\npower-scaling", + prior = "Prior\npower-scaling" + ) + ), + scales = "free", + independent = "all", + switch = "y" + ) + } + + p <- p + + ggplot2::guides( + colour = ggplot2::guide_legend( + title = "Power-scaling alpha", + order = 5 + ) + ) + + if (any(d$pareto_k_value == "High")) { + p <- p + + ggplot2::guides( + linetype = ggplot2::guide_legend( + title = "Pareto k", + order = 10 + ) + ) + } else { + p <- p + + ggplot2::guides(linetype = "none") + } + + if (help_text) { + plot_help_text("ECDF") + } + if (facet_rows == "component") { + p <- p + + ggplot2::theme(legend.position = "bottom") + } + + plots[[i]] <- p + } - variables_per_page <- floor(variables_per_page) + class(plots) <- c("priorsense_plot", class(plots)) - n_plots <- ceiling(nvars / variables_per_page) - plots <- vector(mode = "list", length = n_plots) + if (length(plots) == 1) { + plots <- plots[[1]] + } - for (i in seq_len(n_plots)) { - sub <- ((i - 1) * variables_per_page + 1):min( - i * variables_per_page, - nvars - ) - sub_variable <- variable[sub] - - dsub <- d[d$variable %in% sub_variable, ] - - p <- prepare_plot(dsub, resample = resample, colors = colors, ...) + - ggplot2::guides( - linetype = ggplot2::guide_legend( - title = "Pareto k" - ) - ) + - ggplot2::ylab("ECDF") + - ggplot2::xlab(NULL) + return(plots) + } - p <- p + - ggplot2::stat_ecdf(ggplot2::aes(color = .data[[".powerscale_alpha"]])) - if (facet_rows == "component") { - p <- p + - ggh4x::facet_grid2( - rows = ggplot2::vars(.data$component), - cols = ggplot2::vars(.data$variable), - labeller = ggplot2::labeller( - component = c( - likelihood = "Likelihood\npower-scaling", - prior = "Prior\npower-scaling" - ) - ), - scales = "free", - independent = "all", - switch = "y" - ) - } else { - p <- p + - ggh4x::facet_grid2( - rows = ggplot2::vars(.data$variable), - cols = ggplot2::vars(.data$component), - labeller = ggplot2::labeller( - component = c( - likelihood = "Likelihood\npower-scaling", - prior = "Prior\npower-scaling" - ) - ), - scales = "free", - independent = "all", - switch = "y" - ) - } +powerscale_quantity_segments <- function(summaries) { + group_vars <- interaction( + summaries$variable, + summaries$quantity, + summaries$component, + drop = TRUE, + lex.order = TRUE + ) - if (!(any(d$pareto_k_value == "High"))) { - p <- p + - ggplot2::guides(linetype = "none") - } + grouped_summaries <- split(summaries, group_vars) - if (help_text) { - p <- p + - ggplot2::ggtitle( - label = "Power-scaling sensitivity", - subtitle = paste0( - "Posterior ECDF depending on amount of power-scaling (alpha).\n", - "Overlapping lines indicate low sensitivity.\n", - "Wider gaps between lines indicate greater sensitivity.\n", - "Estimates with high Pareto k (dashed lines) may be inaccurate." + segments <- lapply(grouped_summaries, function(group_data) { + group_data <- group_data[ + order(group_data[[".powerscale_alpha"]]), + , + drop = FALSE + ] + + if (nrow(group_data) < 2) { + return(NULL) + } + + segment_data <- group_data[-nrow(group_data), , drop = FALSE] + + segment_data$x <- group_data[[".powerscale_alpha"]][-nrow(group_data)] + segment_data$xend <- group_data[[".powerscale_alpha"]][-1] + segment_data$y <- group_data$value[-nrow(group_data)] + segment_data$yend <- group_data$value[-1] + + lower_distance <- abs(log(segment_data$x)) + upper_distance <- abs(log(segment_data$xend)) + + lower_status <- group_data$pareto_k_value[-nrow(group_data)] + upper_status <- group_data$pareto_k_value[-1] + + segment_status <- ifelse( + lower_distance > upper_distance, + as.character(lower_status), + ifelse( + upper_distance > lower_distance, + as.character(upper_status), + ifelse( + lower_status == "High" | upper_status == "High", + "High", + "OK" + ) ) - ) - } - if (facet_rows == "component") { - p <- p + - ggplot2::theme(legend.position = "bottom") - } + ) - plots[[i]] <- p - } + segment_data$pareto_k_value <- factor( + segment_status, + levels = c("OK", "High") + ) - class(plots) <- c("priorsense_plot", class(plots)) + segment_data + }) - if (length(plots) == 1) { - plots <- plots[[1]] - } + segments <- Filter(Negate(is.null), segments) - return(plots) - } + if (length(segments) == 0) { + return(summaries[FALSE, , drop = FALSE]) + } + do.call(rbind, segments) +} ##' @rdname powerscale-plots ##' @export powerscale_plot_quantities <- function(x, ...) { - UseMethod("powerscale_plot_quantities") + UseMethod("powerscale_plot_quantities") } ##' @export powerscale_plot_quantities.default <- - function( - x, - variable = NULL, - variables = NULL, - quantity = c("mean", "sd"), - div_measure = "cjs_dist", - length = 11, - resample = FALSE, - measure_args = NULL, - mcse = TRUE, - quantity_args = NULL, - help_text = getOption("priorsense.plot_help_text", TRUE), - colors = NULL, - colours = NULL, - variables_per_page = getOption( - "priorsense.plot_variables_per_page", - 6 - ), - ... - ) { - ps <- powerscale_sequence(x, length = length, ...) - - powerscale_plot_quantities( - ps, - variable = variable, - variables = variables, - quantity = quantity, - div_measure = div_measure, - resample = resample, - measure_args = measure_args, - mcse = mcse, - quantity_args = quantity_args, - help_text = help_text, - colors = colors, - colours = colours, - variables_per_page = variables_per_page - ) - } + function( + x, + variable = NULL, + variables = NULL, + quantity = c("mean", "sd"), + div_measure = "cjs_dist", + length = 11, + resample = FALSE, + measure_args = NULL, + mcse = TRUE, + quantity_args = NULL, + help_text = getOption("priorsense.plot_help_text", TRUE), + colors = NULL, + colours = NULL, + variables_per_page = getOption( + "priorsense.plot_variables_per_page", + 6 + ), + ... + ) { + ps <- powerscale_sequence(x, length = length, ...) + + powerscale_plot_quantities( + ps, + variable = variable, + variables = variables, + quantity = quantity, + div_measure = div_measure, + resample = resample, + measure_args = measure_args, + mcse = mcse, + quantity_args = quantity_args, + help_text = help_text, + colors = colors, + colours = colours, + variables_per_page = variables_per_page + ) + } ##' @rdname powerscale-plots ##' @export powerscale_plot_quantities.powerscaled_sequence <- - function( - x, - variable = NULL, - variables = NULL, - quantity = c("mean", "sd"), - div_measure = "cjs_dist", - resample = FALSE, - measure_args = NULL, - mcse = TRUE, - quantity_args = NULL, - help_text = getOption("priorsense.plot_help_text", TRUE), - colors = NULL, - colours = NULL, - variables_per_page = getOption( - "priorsense.plot_variables_per_page", - 6 - ), - ... - ) { - if (!is.null(variable) && !is.null(variables)) { - checkmate::assert( - if (identical(variable, variables)) { - TRUE - } else { - "must be identical if both provided" - }, - .var.name = "`variable` and `variables`" - ) - } - if (is.null(variable)) { - variable <- variables - } + function( + x, + variable = NULL, + variables = NULL, + quantity = c("mean", "sd"), + div_measure = "cjs_dist", + resample = FALSE, + measure_args = NULL, + mcse = TRUE, + quantity_args = NULL, + help_text = getOption("priorsense.plot_help_text", TRUE), + colors = NULL, + colours = NULL, + variables_per_page = getOption( + "priorsense.plot_variables_per_page", + 6 + ), + ... + ) { + if (!is.null(variable) && !is.null(variables)) { + checkmate::assert( + if (identical(variable, variables)) { + TRUE + } else { + "must be identical if both provided" + }, + .var.name = "`variable` and `variables`" + ) + } + if (is.null(variable)) { + variable <- variables + } - if (!is.null(colors) && !is.null(colours)) { - checkmate::assert( - if (identical(colors, colours)) { - TRUE - } else { - "must be identical if both provided" - }, - .var.name = "`colors` and `colours`" - ) - } - if (is.null(colors)) { - colors <- colours - } + if (!is.null(colors) && !is.null(colours)) { + checkmate::assert( + if (identical(colors, colours)) { + TRUE + } else { + "must be identical if both provided" + }, + .var.name = "`colors` and `colours`" + ) + } + if (is.null(colors)) { + colors <- colours + } - checkmate::assertCharacter(variable, null.ok = TRUE) - checkmate::assertCharacter(quantity) - checkmate::assertCharacter(div_measure, null.ok = TRUE) - checkmate::assertLogical(resample, len = 1) - checkmate::assertList(measure_args, null.ok = TRUE) - checkmate::assertLogical(mcse, len = 1) - checkmate::assertList(quantity_args, null.ok = TRUE) - checkmate::assertLogical(help_text, len = 1) - checkmate::assertCharacter(colors, len = 2, null.ok = TRUE) - checkmate::assert_number(variables_per_page, lower = 1, null.ok = TRUE) - - if (is.null(colors)) { - colors <- default_priorsense_colors()[4:5] - } + checkmate::assertCharacter(variable, null.ok = TRUE) + checkmate::assertCharacter(quantity) + checkmate::assertCharacter(div_measure, null.ok = TRUE) + checkmate::assertLogical(resample, len = 1) + checkmate::assertList(measure_args, null.ok = TRUE) + checkmate::assertLogical(mcse, len = 1) + checkmate::assertList(quantity_args, null.ok = TRUE) + checkmate::assertLogical(help_text, len = 1) + checkmate::assertCharacter(colors, len = 2, null.ok = TRUE) + checkmate::assert_number(variables_per_page, lower = 1, null.ok = TRUE) + + if (is.null(colors)) { + colors <- default_priorsense_colors()[4:5] + } - names(quantity) <- quantity + names(quantity) <- quantity - summ <- summarise_draws( - x, - quantity, - .args = quantity_args, - resample = resample, - div_measures = div_measure, - measure_args = measure_args - ) + summ <- summarise_draws( + x, + quantity, + .args = quantity_args, + resample = resample, + div_measures = div_measure, + measure_args = measure_args + ) - if (is.null(variable)) { - variable <- posterior::variables(x$base_draws) - } else { - # Efficient variable expansion avoiding heavy subset_draws operations - base_vars <- posterior::variables(x$base_draws) + if (is.null(variable)) { + variable <- posterior::variables(x$base_draws) + } else { + # Efficient variable expansion avoiding heavy subset_draws operations + base_vars <- posterior::variables(x$base_draws) - expanded_vars <- lapply(variable, function(v) { - # Keep variable if it matches exactly - if (v %in% base_vars) { - return(v) - } + expanded_vars <- lapply(variable, function(v) { + # Keep variable if it matches exactly + if (v %in% base_vars) { + return(v) + } + + # Check for indexed variables (e.g., "a" -> "a[1]", "a[2]") + prefix <- paste0(v, "[") + indexed_matches <- base_vars[startsWith(base_vars, prefix)] - # Check for indexed variables (e.g., "a" -> "a[1]", "a[2]") - prefix <- paste0(v, "[") - indexed_matches <- base_vars[startsWith(base_vars, prefix)] + if (length(indexed_matches) > 0) { + return(indexed_matches) + } - if (length(indexed_matches) > 0) { - return(indexed_matches) + # Return original string if no match (downstream handles errors) + return(v) + }) + + variable <- unique(unlist(expanded_vars)) } - # Return original string if no match (downstream handles errors) - return(v) - }) + if (mcse) { + quants <- setdiff( + colnames(summ[[1]]), + c( + "variable", + ".powerscale_alpha", + "component", + "pareto_k", + "pareto_kf", + "pareto_k_threshold", + "n_eff", + div_measure + ) + ) - variable <- unique(unlist(expanded_vars)) - } + mcse_functions <- paste0("mcse_", quantity) - if (mcse) { - quants <- setdiff( - colnames(summ[[1]]), - c( - "variable", - ".powerscale_alpha", - "component", - "pareto_k", - "pareto_kf", - "pareto_k_threshold", - "n_eff", - div_measure - ) - ) + base_quantities <- summ[[1]][ + which(summ[[1]][[".powerscale_alpha"]] == 1), + ] - mcse_functions <- paste0("mcse_", quantity) + base_quantities <- unique(base_quantities[c("variable", quants)]) - base_quantities <- summ[[1]][ - which(summ[[1]][[".powerscale_alpha"]] == 1), - ] + base_q <- as.data.frame(base_quantities) - base_quantities <- unique(base_quantities[c("variable", quants)]) + base_q <- as.data.frame(base_quantities) - base_q <- as.data.frame(base_quantities) + base_q <- stats::reshape( + data = base_q, + varying = quants, + direction = "long", + times = quants, + v.names = "value", + timevar = "quantity", + idvar = "variable" + ) - base_q <- as.data.frame(base_quantities) + base_mcse <- posterior::summarise_draws( + x$base_draws, + mcse_functions, + .args = quantity_args + ) + base_mcse <- base_mcse[which(base_mcse$variable %in% variable), ] + base_mcse <- as.data.frame(base_mcse) + + mcse_names <- colnames(base_mcse)[-1] + + base_mcse <- stats::reshape( + data = base_mcse, + varying = mcse_names, + direction = "long", + times = quants, + v.names = "mcse", + timevar = "quantity", + idvar = "variable" + ) - base_q <- stats::reshape( - data = base_q, - varying = quants, - direction = "long", - times = quants, - v.names = "value", - timevar = "quantity", - idvar = "variable" - ) - - base_mcse <- posterior::summarise_draws( - x$base_draws, - mcse_functions, - .args = quantity_args - ) - base_mcse <- base_mcse[which(base_mcse$variable %in% variable), ] - base_mcse <- as.data.frame(base_mcse) - - mcse_names <- colnames(base_mcse)[-1] - - base_mcse <- stats::reshape( - data = base_mcse, - varying = mcse_names, - direction = "long", - times = quants, - v.names = "mcse", - timevar = "quantity", - idvar = "variable" - ) - - base_mcse <- merge(base_q, base_mcse) - base_mcse$mcse_min <- base_mcse$value - 2 * base_mcse$mcse - base_mcse$mcse_max <- base_mcse$value + 2 * base_mcse$mcse - } else { - base_mcse <- NULL - } + base_mcse <- merge(base_q, base_mcse) + base_mcse$mcse_min <- base_mcse$value - 2 * base_mcse$mcse + base_mcse$mcse_max <- base_mcse$value + 2 * base_mcse$mcse + } else { + base_mcse <- NULL + } - powerscale_summary_plot( - summ, - variable = variable, - base_mcse = base_mcse, - help_text = help_text, - colors = colors, - variables_per_page = variables_per_page, - ... - ) - } + powerscale_summary_plot( + summ, + variable = variable, + base_mcse = base_mcse, + help_text = help_text, + colors = colors, + variables_per_page = variables_per_page, + ... + ) + } ##' power-scale summary plot ##' ##' internal function for powerscale_plot_quantities @@ -985,212 +1072,223 @@ powerscale_plot_quantities.powerscaled_sequence <- ##' @keywords internal ##' @noRd powerscale_summary_plot <- function( - x, - variable, - base_mcse = NULL, - help_text, - colors, - variables_per_page, - ... + x, + variable, + base_mcse = NULL, + help_text, + colors, + variables_per_page, + ... ) { - nvars <- length(variable) - - if (is.null(variables_per_page) || is.infinite(variables_per_page)) { - variables_per_page <- nvars - } - - variables_per_page <- floor(variables_per_page) - - n_plots <- ceiling(nvars / variables_per_page) - plots <- vector(mode = "list", length = n_plots) - - pareto_k_colours <- colors - - # get default quantities - quantities <- setdiff( - colnames(x[[1]]), - c( - "variable", - ".powerscale_alpha", - "component", - "pareto_k", - "pareto_kf", - "n_eff", - "pareto_k_threshold" - ) - ) - - for (i in seq_len(n_plots)) { - sub <- ((i - 1) * variables_per_page + 1):min(i * variables_per_page, nvars) - sub_variable <- variable[sub] - # select only specified variables - xsub <- x[[1]][x[[1]][["variable"]] %in% sub_variable, ] - - sub_mcse <- base_mcse[base_mcse$variable %in% sub_variable, ] - - # reshape quantities for plotting - summaries <- stats::reshape( - data = xsub, - varying = quantities, - direction = "long", - times = quantities, - v.names = "value", - timevar = "quantity" - ) + nvars <- length(variable) - summaries$pareto_k_value <- ifelse( - summaries$pareto_k > summaries$pareto_k_threshold, - "High", - "OK" - ) + if (is.null(variables_per_page) || is.infinite(variables_per_page)) { + variables_per_page <- nvars + } - summaries$pareto_k_value <- factor( - summaries$pareto_k_value, - levels = c("OK", "High") + variables_per_page <- floor(variables_per_page) + + n_plots <- ceiling(nvars / variables_per_page) + plots <- vector(mode = "list", length = n_plots) + + pareto_k_colours <- colors + + # get default quantities + quantities <- setdiff( + colnames(x[[1]]), + c( + "variable", + ".powerscale_alpha", + "component", + "pareto_k", + "pareto_kf", + "n_eff", + "pareto_k_threshold" + ) ) - # subset for plotting points at ends of lines - points <- summaries[ - summaries[[".powerscale_alpha"]] == - min(summaries[[".powerscale_alpha"]]) | - summaries[[".powerscale_alpha"]] == - max(summaries[[".powerscale_alpha"]]), - ] - - p <- ggplot2::ggplot( - data = summaries, - mapping = ggplot2::aes(x = .data[[".powerscale_alpha"]], y = .data$value) - ) + - ggplot2::geom_line(ggplot2::aes( - color = .data$pareto_k_value, - group = .data$component - )) + - ggh4x::facet_grid2( - rows = ggplot2::vars(factor( - .data$variable, - levels = unique(.data$variable) - )), - cols = ggplot2::vars(factor( - .data$quantity, - levels = unique(.data$quantity) - )), - scales = "free", - switch = "y", - independent = "all" - ) + - ggplot2::geom_point( - ggplot2::aes( - x = .data[[".powerscale_alpha"]], - y = .data$value, - shape = .data$component, - colour = .data$pareto_k_value - ), - fill = "white", - size = 3, - data = points - ) + - ggplot2::scale_shape_manual( - values = c("likelihood" = 22, "prior" = 15) - ) + - ggplot2::scale_color_manual(values = pareto_k_colours) + - ggplot2::guides( - color = ggplot2::guide_legend( - title = "Pareto k", - override.aes = list(shape = 15) + for (i in seq_len(n_plots)) { + sub <- ((i - 1) * variables_per_page + 1):min( + i * variables_per_page, + nvars + ) + sub_variable <- variable[sub] + # select only specified variables + xsub <- x[[1]][x[[1]][["variable"]] %in% sub_variable, ] + + sub_mcse <- base_mcse[base_mcse$variable %in% sub_variable, ] + + # reshape quantities for plotting + summaries <- stats::reshape( + data = xsub, + varying = quantities, + direction = "long", + times = quantities, + v.names = "value", + timevar = "quantity" ) - ) + - ggplot2::ylab(NULL) + - ggplot2::scale_x_continuous( - trans = "log2", - limits = c( - 0.95 * min(summaries[[".powerscale_alpha"]]), - 1.05 * max(summaries[[".powerscale_alpha"]]) - ), - breaks = c( - min(summaries[[".powerscale_alpha"]]), - 1, - max(summaries[[".powerscale_alpha"]]) - ), - labels = round( - c( - min(summaries[[".powerscale_alpha"]]), - 1, - max(summaries[[".powerscale_alpha"]]) - ), - digits = 3 - ), - name = "Power-scaling alpha" - ) - if (!(any(summaries$pareto_k_value == "High"))) { - p <- p + - ggplot2::guides( - colour = "none" + summaries$pareto_k_value <- ifelse( + summaries$pareto_k > summaries$pareto_k_threshold, + "High", + "OK" ) - } - if (help_text) { - p <- p + - ggplot2::ggtitle( - label = "Power-scaling sensitivity", - subtitle = paste0( - "Posterior quantities depending on amount of power-scaling (alpha).\n", - "Horizontal lines indicate low sensitivity.\n", - "Steeper lines indicate greater sensitivity.\n", - "Estimates with high Pareto k (highlighted) may be inaccurate." - ) + summaries$pareto_k_value <- factor( + summaries$pareto_k_value, + levels = c("OK", "High") ) - } - if (!is.null(sub_mcse)) { - p <- p + - ggplot2::scale_linetype_manual(values = "dashed", name = NULL) + - ggplot2::geom_hline( - ggplot2::aes( - yintercept = .data$mcse_min, - linetype = "+/-2MCSE" - ), - data = sub_mcse, - color = "black" + # subset for plotting points at ends of lines + points <- summaries[ + summaries[[".powerscale_alpha"]] == + min(summaries[[".powerscale_alpha"]]) | + summaries[[".powerscale_alpha"]] == + max(summaries[[".powerscale_alpha"]]), + ] + + line_segments <- powerscale_quantity_segments(summaries) + + p <- ggplot2::ggplot( + data = summaries, + mapping = ggplot2::aes( + x = .data[[".powerscale_alpha"]], + y = .data$value + ) ) + - ggplot2::geom_hline( - ggplot2::aes( - yintercept = .data$mcse_max, - linetype = "+/-2MCSE" - ), - data = sub_mcse, - color = "black" - ) - } + ggplot2::geom_segment( + data = line_segments, + mapping = ggplot2::aes( + x = .data$x, + xend = .data$xend, + y = .data$y, + yend = .data$yend, + color = .data$pareto_k_value + ) + ) + + ggh4x::facet_grid2( + rows = ggplot2::vars(factor( + .data$variable, + levels = unique(.data$variable) + )), + cols = ggplot2::vars(factor( + .data$quantity, + levels = unique(.data$quantity) + )), + scales = "free", + switch = "y", + independent = "all" + ) + + ggplot2::geom_point( + ggplot2::aes( + x = .data[[".powerscale_alpha"]], + y = .data$value, + shape = .data$component, + colour = .data$pareto_k_value + ), + fill = "white", + size = 3, + data = points + ) + + ggplot2::scale_shape_manual( + values = c("likelihood" = 22, "prior" = 15) + ) + + ggplot2::scale_color_manual(values = pareto_k_colours) + + ggplot2::guides( + color = ggplot2::guide_legend( + title = "Pareto k", + override.aes = list(shape = 15), + order = 2 + ), + shape = ggplot2::guide_legend( + title = "Component", + order = 1 + ) + ) + + ggplot2::ylab(NULL) + + ggplot2::scale_x_continuous( + trans = "log2", + limits = c( + 0.95 * min(summaries[[".powerscale_alpha"]]), + 1.05 * max(summaries[[".powerscale_alpha"]]) + ), + breaks = c( + min(summaries[[".powerscale_alpha"]]), + 1, + max(summaries[[".powerscale_alpha"]]) + ), + labels = round( + c( + min(summaries[[".powerscale_alpha"]]), + 1, + max(summaries[[".powerscale_alpha"]]) + ), + digits = 3 + ), + name = "Power-scaling alpha" + ) - plots[[i]] <- p - } + if (!(any(summaries$pareto_k_value == "High", na.rm = TRUE))) { + p <- p + + ggplot2::guides( + colour = "none" + ) + } - class(plots) <- c("priorsense_plot", class(plots)) + if (help_text) { + plot_help_text("quantities") + } + + if (!is.null(sub_mcse)) { + p <- p + + ggplot2::scale_linetype_manual(values = "dashed", name = NULL) + + ggplot2::geom_hline( + ggplot2::aes( + yintercept = .data$mcse_min, + linetype = "+/-2MCSE" + ), + data = sub_mcse, + color = "black" + ) + + ggplot2::geom_hline( + ggplot2::aes( + yintercept = .data$mcse_max, + linetype = "+/-2MCSE" + ), + data = sub_mcse, + color = "black" + ) + + ggplot2::guides(linetype = ggplot2::guide_legend(order = 3)) + } + + plots[[i]] <- p + } + + class(plots) <- c("priorsense_plot", class(plots)) - if (length(plots) == 1) { - plots <- plots[[1]] - } + if (length(plots) == 1) { + plots <- plots[[1]] + } - return(plots) + return(plots) } ##' @exportS3Method plot.priorsense_plot <- function( - x, - ask = getOption("priorsense.plot_ask", TRUE), - ... + x, + ask = getOption("priorsense.plot_ask", TRUE), + ... ) { - grDevices::devAskNewPage(ask = FALSE) - on.exit(grDevices::devAskNewPage(ask = FALSE)) + grDevices::devAskNewPage(ask = FALSE) + on.exit(grDevices::devAskNewPage(ask = FALSE)) - for (i in seq_along(x)) { - plot(x[[i]], newpage = TRUE) - if (i == 1) { - grDevices::devAskNewPage(ask = ask) + for (i in seq_along(x)) { + plot(x[[i]], newpage = TRUE, help_text = FALSE) + if (i == 1) { + grDevices::devAskNewPage(ask = ask) + } } - } - invisible(x) + invisible(x) } @@ -1201,10 +1299,34 @@ print.priorsense_plot <- plot.priorsense_plot ##' @exportS3Method ##' @rdname powerscale-plots plot.powerscaled_sequence <- function( - x, - type = c("dens", "ecdf", "quantities"), - ... + x, + type = c("dens", "ecdf", "quantities"), + ... ) { - type <- match.arg(type) - do.call(paste0("powerscale_plot_", type), args = list(x = x, ...)) + type <- match.arg(type) + do.call(paste0("powerscale_plot_", type), args = list(x = x, ...)) +} + + +plot_help_text <- function(type) { + cli::cli_h2("Power-scaling sensitivity {type} plot:") + if (type %in% c("density", "ECDF")) { + cli::cli_text( + "The plot shows posterior {type} depending on ", + "the degree of power-scaling (alpha).\n", + "Overlapping lines indicate low sensitivity.\n", + "Wider gaps between lines indicate greater sensitivity.\n", + "Estimates with high Pareto k (dashed lines) may be inaccurate.\n" + ) + } else if (type == "quantities") { + cli::cli_text( + "Posterior quantities depending on amount of power-scaling (alpha).\n", + "Horizontal lines indicate low sensitivity.\n", + "Steeper lines indicate greater sensitivity.\n", + "Estimates with high Pareto k (highlighted) may be inaccurate.\n" + ) + } + cli::cli_text( + "Disable this help text with `help_text = FALSE` or `options(priorsense.plot_help_text = FALSE)`" + ) } diff --git a/R/powerscale_sensitivity.R b/R/powerscale_sensitivity.R index 9eb07605..abc2b629 100644 --- a/R/powerscale_sensitivity.R +++ b/R/powerscale_sensitivity.R @@ -8,6 +8,7 @@ ##' @param x Model fit object or priorsense_data object. ##' @param ... Further arguments passed to functions. ##' @param variable Character vector of variables to check. +##' @param variables alias of `variable`. ##' @param lower_alpha Lower alpha value for gradient calculation. ##' @param upper_alpha Upper alpha value for gradient calculation. ##' @param component Character vector specifying component(s) to scale @@ -32,195 +33,213 @@ ##' powerscale_sensitivity(ex$draws) ##' @export powerscale_sensitivity <- function(x, ...) { - UseMethod("powerscale_sensitivity") + UseMethod("powerscale_sensitivity") } ##' @rdname powerscale-sensitivity ##' @export powerscale_sensitivity.default <- function( - x, - variable = NULL, - lower_alpha = 0.99, - upper_alpha = 1.01, - div_measure = "cjs_dist", - measure_args = list(), - component = c( - "prior", - "likelihood" - ), - sensitivity_threshold = 0.05, - moment_match = FALSE, - k_threshold = 0.5, - resample = FALSE, - transform = NULL, - prediction = NULL, - prior_selection = NULL, - likelihood_selection = NULL, - log_prior_name = "lprior", - log_lik_name = "log_lik", - separator = "_", - num_args = NULL, - ... -) { - psd <- create_priorsense_data( - x = x, - log_prior_name = log_prior_name, - log_lik_name = log_lik_name, - ... - ) - - powerscale_sensitivity.priorsense_data( - psd, - variable = variable, - lower_alpha = lower_alpha, - upper_alpha = upper_alpha, - div_measure = div_measure, - measure_args = measure_args, - component = component, - sensitivity_threshold = sensitivity_threshold, - moment_match = moment_match, - k_threshold = k_threshold, - resample = resample, - transform = transform, - prediction = prediction, - prior_selection = prior_selection, - likelihood_selection = likelihood_selection, - separator = separator, - num_args = num_args, + x, + variable = NULL, + variables = NULL, + lower_alpha = 0.99, + upper_alpha = 1.01, + div_measure = "cjs_dist", + measure_args = list(), + component = c( + "prior", + "likelihood" + ), + sensitivity_threshold = 0.05, + moment_match = FALSE, + k_threshold = 0.5, + resample = FALSE, + transform = NULL, + prediction = NULL, + prior_selection = NULL, + likelihood_selection = NULL, + log_prior_name = "lprior", + log_lik_name = "log_lik", + separator = "_", + num_args = NULL, ... - ) +) { + psd <- create_priorsense_data( + x = x, + log_prior_name = log_prior_name, + log_lik_name = log_lik_name, + ... + ) + + powerscale_sensitivity.priorsense_data( + psd, + variable = variable, + variables = variables, + lower_alpha = lower_alpha, + upper_alpha = upper_alpha, + div_measure = div_measure, + measure_args = measure_args, + component = component, + sensitivity_threshold = sensitivity_threshold, + moment_match = moment_match, + k_threshold = k_threshold, + resample = resample, + transform = transform, + prediction = prediction, + prior_selection = prior_selection, + likelihood_selection = likelihood_selection, + separator = separator, + num_args = num_args, + ... + ) } ##' @rdname powerscale-sensitivity ##' @export powerscale_sensitivity.priorsense_data <- function( - x, - variable = NULL, - lower_alpha = 0.99, - upper_alpha = 1.01, - div_measure = "cjs_dist", - measure_args = list(), - component = c( - "prior", - "likelihood" - ), - sensitivity_threshold = 0.05, - moment_match = FALSE, - k_threshold = 0.5, - resample = FALSE, - transform = NULL, - prediction = NULL, - prior_selection = NULL, - likelihood_selection = NULL, - separator = "_", - num_args = NULL, - ... -) { - component <- tolower(component) - - # input checks - checkmate::assertCharacter(variable, null.ok = TRUE) - checkmate::assertNumber(lower_alpha, lower = 0, upper = 1) - checkmate::assertNumber(upper_alpha, lower = 1) - checkmate::assertCharacter(div_measure, len = 1) - checkmate::assertList(measure_args) - checkmate::assertLogical(moment_match, len = 1) - checkmate::assertSubset(component, c("prior", "likelihood")) - checkmate::assertNumber(sensitivity_threshold, lower = 0) - checkmate::assertNumber(k_threshold, null.ok = TRUE) - checkmate::assertLogical(resample, len = 1) - checkmate::assertCharacter(transform, null.ok = TRUE, len = 1) - checkmate::assertFunction(prediction, null.ok = TRUE) - checkmate::assertCharacter(separator) - - gradients <- powerscale_gradients( - x = x, - variable = variable, - component = component, - type = "divergence", - lower_alpha = lower_alpha, - upper_alpha = upper_alpha, - moment_match = moment_match, - div_measure = div_measure, - measure_args = measure_args, - transform = transform, - resample = resample, - prediction = prediction, - prior_selection = prior_selection, - likelihood_selection = likelihood_selection, - separator = separator, + x, + variable = NULL, + variables = NULL, + lower_alpha = 0.99, + upper_alpha = 1.01, + div_measure = "cjs_dist", + measure_args = list(), + component = c( + "prior", + "likelihood" + ), + sensitivity_threshold = 0.05, + moment_match = FALSE, + k_threshold = 0.5, + resample = FALSE, + transform = NULL, + prediction = NULL, + prior_selection = NULL, + likelihood_selection = NULL, + separator = "_", + num_args = NULL, ... - ) - - prior_sense <- gradients$divergence$prior[[2]] - lik_sense <- gradients$divergence$likelihood[[2]] - - if (is.null(lik_sense)) { - lik_sense <- NA - } - - if (is.null(prior_sense)) { - prior_sense <- NA - } - - varnames <- unique(c( - as.character(gradients$divergence$prior$variable), - as.character(gradients$divergence$likelihood$variable) - )) - - sense <- data.frame( - variable = varnames, - prior = prior_sense, - likelihood = lik_sense - ) - - # categorise variables has prior-data conflict or uninformative - # likelihood - - sense$diagnosis <- ifelse( - sense$prior >= sensitivity_threshold & - sense$likelihood >= sensitivity_threshold, - "potential prior-likelihood conflict", - ifelse( - sense$prior > sensitivity_threshold & - sense$likelihood < sensitivity_threshold, - "potential strong prior / weak likelihood", - "-" +) { + component <- tolower(component) + + # input checks + checkmate::assertCharacter(variable, null.ok = TRUE) + checkmate::assertCharacter(variables, null.ok = TRUE) + checkmate::assertNumber(lower_alpha, lower = 0, upper = 1) + checkmate::assertNumber(upper_alpha, lower = 1) + checkmate::assertCharacter(div_measure, len = 1) + checkmate::assertList(measure_args) + checkmate::assertLogical(moment_match, len = 1) + checkmate::assertSubset(component, c("prior", "likelihood")) + checkmate::assertNumber(sensitivity_threshold, lower = 0) + checkmate::assertNumber(k_threshold, null.ok = TRUE) + checkmate::assertLogical(resample, len = 1) + checkmate::assertCharacter(transform, null.ok = TRUE, len = 1) + checkmate::assertFunction(prediction, null.ok = TRUE) + checkmate::assertCharacter(separator) + + if (!is.null(variable) && !is.null(variables)) { + checkmate::assert( + if (identical(variable, variables)) { + TRUE + } else { + "must be identical if both provided" + }, + .var.name = "`variable` and `variables`" + ) + } + if (is.null(variable)) { + variable <- variables + } + + gradients <- powerscale_gradients( + x = x, + variable = variable, + component = component, + type = "divergence", + lower_alpha = lower_alpha, + upper_alpha = upper_alpha, + moment_match = moment_match, + div_measure = div_measure, + measure_args = measure_args, + transform = transform, + resample = resample, + prediction = prediction, + prior_selection = prior_selection, + likelihood_selection = likelihood_selection, + separator = separator, + ... + ) + + prior_sense <- gradients$divergence$prior[[2]] + lik_sense <- gradients$divergence$likelihood[[2]] + + if (is.null(lik_sense)) { + lik_sense <- NA + } + + if (is.null(prior_sense)) { + prior_sense <- NA + } + + varnames <- unique(c( + as.character(gradients$divergence$prior$variable), + as.character(gradients$divergence$likelihood$variable) + )) + + sense <- data.frame( + variable = varnames, + prior = prior_sense, + likelihood = lik_sense + ) + + # categorise variables has prior-data conflict or uninformative + # likelihood + + sense$diagnosis <- ifelse( + sense$prior >= sensitivity_threshold & + sense$likelihood >= sensitivity_threshold, + "potential prior-likelihood conflict", + ifelse( + sense$prior > sensitivity_threshold & + sense$likelihood < sensitivity_threshold, + "potential strong prior / weak likelihood", + "-" + ) ) - ) - out <- sense + out <- sense - class(out) <- c("powerscaled_sensitivity_summary", class(out)) + class(out) <- c("powerscaled_sensitivity_summary", class(out)) - attr(out, "num_args") <- num_args - attr(out, "div_measure") <- div_measure - attr(out, "loadings") <- gradients$loadings - attr(out, "prior_selection") <- prior_selection - attr(out, "likelihood_selection") <- likelihood_selection + attr(out, "num_args") <- num_args + attr(out, "div_measure") <- div_measure + attr(out, "loadings") <- gradients$loadings + attr(out, "prior_selection") <- prior_selection + attr(out, "likelihood_selection") <- likelihood_selection - return(out) + return(out) } ##' @rdname powerscale-sensitivity ##' @export powerscale_sensitivity.CmdStanFit <- function(x, ...) { - psd <- create_priorsense_data.CmdStanFit(x) + psd <- create_priorsense_data.CmdStanFit(x) - powerscale_sensitivity.priorsense_data( - psd, - ... - ) + powerscale_sensitivity.priorsense_data( + psd, + ... + ) } ##' @rdname powerscale-sensitivity ##' @export powerscale_sensitivity.stanfit <- function(x, ...) { - psd <- create_priorsense_data.stanfit(x, ...) + psd <- create_priorsense_data.stanfit(x, ...) - powerscale_sensitivity.priorsense_data( - psd, - ... - ) + powerscale_sensitivity.priorsense_data( + psd, + ... + ) } diff --git a/README.md b/README.md index 4a52b48d..76cd803f 100644 --- a/README.md +++ b/README.md @@ -9,10 +9,10 @@ state and is being actively developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/#active) [![priorsense status -badge](https://n-kall.r-universe.dev/badges/priorsense)](https://n-kall.r-universe.dev) +badge](https://stan-dev.r-universe.dev/badges/priorsense)](https://stan-dev.r-universe.dev) [![priorsense CRAN badge](https://www.r-pkg.org/badges/version/priorsense)](https://cran.r-project.org/package=priorsense) -[![R-CMD-check](https://github.com/n-kall/priorsense/workflows/R-CMD-check/badge.svg)](https://github.com/n-kall/priorsense/actions) +[![R-CMD-check](https://github.com/stan-dev/priorsense/workflows/R-CMD-check/badge.svg)](https://github.com/stan-dev/priorsense/actions) [![Status at rOpenSci Software Peer Review](https://badges.ropensci.org/704_status.svg)](https://github.com/ropensci/software-review/issues/704) [![DOI](https://joss.theoj.org/papers/10.21105/joss.11036/status.svg)](https://doi.org/10.21105/joss.11036) @@ -42,11 +42,20 @@ al. (2023)](https://doi.org/10.1007/s11222-023-10366-5). ### Resources - Check the [getting started - vignette](https://n-kall.github.io/priorsense/articles/getting_started.html) + vignette](https://mc-stan.org/priorsense/articles/getting_started.html) for a simple example - For a more detailed modelling example see - [here](https://n-kall.github.io/priorsense/articles/airquality.html) + [here](https://mc-stan.org/priorsense/articles/airquality.html) + +### Resources + +- [mc-stan.org/priorsense](https://mc-stan.org/priorsense) (online + documentation, vignettes) +- [Ask a question](https://discourse.mc-stan.org) (Stan Forums on + Discourse) +- [Open an issue](https://github.com/stan-dev/priorsense/issues) (GitHub + issues for bug reports, feature requests) ### Installation @@ -61,14 +70,16 @@ with: ``` r # install.packages("pak") -pak::pkg_install("n-kall/priorsense@development") +pak::pkg_install("stan-dev/priorsense@development") ``` ### Contributing Contributions are welcome! If you find a bug or have an idea for a feature, open an issue. If you are able to fix an issue, fork the -repository and make a pull request to the `development` branch. +repository and make a pull request to the `development` branch. Read +[CONTRIBUTING.md](https://github.com/stan-dev/priorsense/blob/main/.github/CONTRIBUTING.md) +for more details. ### References diff --git a/README.qmd b/README.qmd index 935067ac..5439c458 100644 --- a/README.qmd +++ b/README.qmd @@ -13,9 +13,9 @@ ggplot2::theme_set(bayesplot::theme_default(base_family = "sans")) [![Project Status: Active – The project has reached a stable, usable state and is being actively developed.](https://www.repostatus.org/badges/latest/active.svg)](https://www.repostatus.org/#active) -[![priorsense status badge](https://n-kall.r-universe.dev/badges/priorsense)](https://n-kall.r-universe.dev) +[![priorsense status badge](https://stan-dev.r-universe.dev/badges/priorsense)](https://stan-dev.r-universe.dev) [![priorsense CRAN badge](https://www.r-pkg.org/badges/version/priorsense)](https://cran.r-project.org/package=priorsense) -[![R-CMD-check](https://github.com/n-kall/priorsense/workflows/R-CMD-check/badge.svg)](https://github.com/n-kall/priorsense/actions) +[![R-CMD-check](https://github.com/stan-dev/priorsense/workflows/R-CMD-check/badge.svg)](https://github.com/stan-dev/priorsense/actions) [![Status at rOpenSci Software Peer Review](https://badges.ropensci.org/704_status.svg)](https://github.com/ropensci/software-review/issues/704) [![DOI](https://joss.theoj.org/papers/10.21105/joss.11036/status.svg)](https://doi.org/10.21105/joss.11036) @@ -41,11 +41,17 @@ from other software. Power-scaling sensitivity analysis checks are described in ### Resources - Check the [getting started - vignette](https://n-kall.github.io/priorsense/articles/getting_started.html) + vignette](https://mc-stan.org/priorsense/articles/getting_started.html) for a simple example - For a more detailed modelling example see - [here](https://n-kall.github.io/priorsense/articles/airquality.html) + [here](https://mc-stan.org/priorsense/articles/airquality.html) + +### Resources + +* [mc-stan.org/priorsense](https://mc-stan.org/priorsense) (online documentation, vignettes) +* [Ask a question](https://discourse.mc-stan.org) (Stan Forums on Discourse) +* [Open an issue](https://github.com/stan-dev/priorsense/issues) (GitHub issues for bug reports, feature requests) ### Installation @@ -61,14 +67,14 @@ Download the development version from [GitHub](https://github.com/) with: ```{r} #| eval: false # install.packages("pak") -pak::pkg_install("n-kall/priorsense@development") +pak::pkg_install("stan-dev/priorsense@development") ``` ### Contributing -Contributions are welcome! If you find a bug or have an idea for a feature, -open an issue. If you are able to fix an issue, fork the repository and make a -pull request to the `development` branch. +Contributions are welcome! If you find a bug or have an idea for a feature, open +an issue. If you are able to fix an issue, fork the repository and make a pull +request to the `development` branch. Read [CONTRIBUTING.md](https://github.com/stan-dev/priorsense/blob/main/.github/CONTRIBUTING.md) for more details. ### References diff --git a/_pkgdown.yml b/_pkgdown.yml index 4f6227aa..81e9b0bf 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -1,4 +1,97 @@ -url: https://n-kall.github.io/priorsense/ +url: https://mc-stan.org/priorsense/ + +destination: "." + +development: + mode: auto + template: - bootstrap: 5 + package: pkgdownconfig + +navbar: + title: "priorsense" + + structure: + left: [home, articles, functions, news, pkgs, stan] + right: [search, bluesky, forum, github, lightswitch] + + components: + pkgs: + text: Other Packages + menu: + - text: bayesplot + href: https://mc-stan.org/bayesplot + - text: cmdstanr + href: https://mc-stan.org/cmdstanr + - text: loo + href: https://mc-stan.org/loo + - text: posterior + href: https://mc-stan.org/posterior + - text: projpred + href: https://mc-stan.org/projpred + - text: rstan + href: https://mc-stan.org/rstan + - text: rstanarm + href: https://mc-stan.org/rstanarm + - text: rstantools + href: https://mc-stan.org/rstantools + - text: shinystan + href: https://mc-stan.org/shinystan + +articles: + - title: Getting started + desc: | + These pages demonstrate how to use the **priorsense** package to check prior and likelihood sensitivity for Bayesian models fit using MCMC. + contents: + - getting_started + - priorsense_with_stan + - priorsense_with_brms + - priorsense_with_jags + - priorsense_with_nimble + - title: Additional examples + desc: | + These pages demonstrate how to use the **priorsense** package in more detail. + contents: + - airquality + - title: Further details + contents: + - sensitivity_diagnostic + - quantity_of_interest + - papers + +reference: + - title: Package description + contents: + - priorsense-package + - title: Numerical power-scaling sensitivity checks + desc: | + Numerical sensitivity checks for prior and likelihood + contents: + - powerscale_sensitivity + - powerscale_derivative + - title: Graphical checks + desc: | + Plots for power-scaling sensitivity checks + contents: + - powerscale_plot_dens + - powerscale_plot_ecdf + - powerscale_plot_quantities + - title: Example models + desc: | + Provides example model code and data for different PPLs compatible with priorsense + contents: + - example_powerscale_model + - title: Other functions for exploring sensitivity + contents: + - powerscale + - powerscale_sequence + - powerscale_gradients + - create_priorsense_data + - title: Helper functions + contents: + - create_priorsense_data + - log_lik_draws + - log_prior_draws + - predictions_as_draws + - cjs_dist diff --git a/man/figures/logo.svg b/man/figures/logo.svg new file mode 100644 index 00000000..496f0402 --- /dev/null +++ b/man/figures/logo.svg @@ -0,0 +1 @@ + \ No newline at end of file diff --git a/man/powerscale-sensitivity.Rd b/man/powerscale-sensitivity.Rd index ff071f65..26647b72 100644 --- a/man/powerscale-sensitivity.Rd +++ b/man/powerscale-sensitivity.Rd @@ -14,6 +14,7 @@ powerscale_sensitivity(x, ...) \method{powerscale_sensitivity}{default}( x, variable = NULL, + variables = NULL, lower_alpha = 0.99, upper_alpha = 1.01, div_measure = "cjs_dist", @@ -37,6 +38,7 @@ powerscale_sensitivity(x, ...) \method{powerscale_sensitivity}{priorsense_data}( x, variable = NULL, + variables = NULL, lower_alpha = 0.99, upper_alpha = 1.01, div_measure = "cjs_dist", @@ -66,6 +68,8 @@ powerscale_sensitivity(x, ...) \item{variable}{Character vector of variables to check.} +\item{variables}{alias of \code{variable}.} + \item{lower_alpha}{Lower alpha value for gradient calculation.} \item{upper_alpha}{Upper alpha value for gradient calculation.} diff --git a/man/priorsense-package.Rd b/man/priorsense-package.Rd index f1faec00..b26db283 100644 --- a/man/priorsense-package.Rd +++ b/man/priorsense-package.Rd @@ -59,11 +59,11 @@ Computing}. 31(16). \code{doi:10.1007/s11222-020-09982-2} \link{powerscale-plots} } \author{ -\strong{Maintainer}: Noa Kallioinen \email{noa.kallioinen@helsinki.fi} (\href{https://orcid.org/0000-0003-1586-8382}{ORCID}) [copyright holder] +\strong{Maintainer}: Noa Kallioinen \email{noa.kallioinen@helsinki.fi} (\href{https://orcid.org/0000-0003-1586-8382}{ORCID}) Authors: \itemize{ - \item Noa Kallioinen \email{noa.kallioinen@helsinki.fi} (\href{https://orcid.org/0000-0003-1586-8382}{ORCID}) [copyright holder] + \item Noa Kallioinen \email{noa.kallioinen@helsinki.fi} (\href{https://orcid.org/0000-0003-1586-8382}{ORCID}) \item Topi Paananen (\href{https://orcid.org/0000-0002-6542-407X}{ORCID}) \item Paul-Christian Bürkner (\href{https://orcid.org/0000-0001-5765-8995}{ORCID}) \item Aki Vehtari (\href{https://orcid.org/0000-0003-2164-9469}{ORCID}) diff --git a/pkgdown/favicon/favicon.svg b/pkgdown/favicon/favicon.svg new file mode 100644 index 00000000..c9fa0f5c --- /dev/null +++ b/pkgdown/favicon/favicon.svg @@ -0,0 +1,3 @@ + \ No newline at end of file diff --git a/tests/testthat/_snaps/plots/normal-model-density-plot.svg b/tests/testthat/_snaps/plots/normal-model-density-plot.svg index c699a612..87d8e911 100644 --- a/tests/testthat/_snaps/plots/normal-model-density-plot.svg +++ b/tests/testthat/_snaps/plots/normal-model-density-plot.svg @@ -21,246 +21,246 @@ - - + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - + + - - -mu + + +mu - - + + - - -sigma + + +sigma - - + + - - -Prior -power-scaling + + +Prior +power-scaling - - + + - - -Likelihood -power-scaling + + +Likelihood +power-scaling @@ -281,24 +281,24 @@ 1.0 1.5 2.0 - - - - - -8.0 -8.5 -9.0 -9.5 -10.0 - - - - -0.5 -1.0 -1.5 -2.0 + + + + + +8.0 +8.5 +9.0 +9.5 +10.0 + + + + +0.5 +1.0 +1.5 +2.0 Power-scaling alpha @@ -310,10 +310,6 @@ 0.8 1 1.25 -Posterior density estimates depending on amount of power-scaling (alpha). -Overlapping lines indicate low sensitivity. -Wider gaps between lines indicate greater sensitivity. -Estimates with high Pareto k (dashed lines) may be inaccurate. -Power-scaling sensitivity +Normal model density plot diff --git a/tests/testthat/_snaps/plots/normal-model-ecdf-plot.svg b/tests/testthat/_snaps/plots/normal-model-ecdf-plot.svg index af711148..53bb07d8 100644 --- a/tests/testthat/_snaps/plots/normal-model-ecdf-plot.svg +++ b/tests/testthat/_snaps/plots/normal-model-ecdf-plot.svg @@ -21,166 +21,166 @@ - - + + - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + - - + + - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + - - + + - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + - - + + - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + - - + + - - -mu + + +mu - - + + - - -sigma + + +sigma - - + + - - -Prior -power-scaling + + +Prior +power-scaling - - + + - - -Likelihood -power-scaling + + +Likelihood +power-scaling @@ -197,61 +197,61 @@ 0.5 1.0 1.5 - - - - -8.5 -9.0 -9.5 -10.0 - - - -0.5 -1.0 -1.5 -0.00 -0.25 -0.50 -0.75 -1.00 - - - - - -0.00 -0.25 -0.50 -0.75 -1.00 - - - - - -0.00 -0.25 -0.50 -0.75 -1.00 - - - - - -0.00 -0.25 -0.50 -0.75 -1.00 - - - - - -ECDF + + + + +8.5 +9.0 +9.5 +10.0 + + + +0.5 +1.0 +1.5 +0.00 +0.25 +0.50 +0.75 +1.00 + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 + + + + + +0.00 +0.25 +0.50 +0.75 +1.00 + + + + + +ECDF Power-scaling alpha @@ -263,10 +263,6 @@ 0.8 1 1.25 -Posterior ECDF depending on amount of power-scaling (alpha). -Overlapping lines indicate low sensitivity. -Wider gaps between lines indicate greater sensitivity. -Estimates with high Pareto k (dashed lines) may be inaccurate. -Power-scaling sensitivity +Normal model ecdf plot diff --git a/tests/testthat/_snaps/plots/normal-model-quantities-plot.svg b/tests/testthat/_snaps/plots/normal-model-quantities-plot.svg index 954c72a2..a9503f83 100644 --- a/tests/testthat/_snaps/plots/normal-model-quantities-plot.svg +++ b/tests/testthat/_snaps/plots/normal-model-quantities-plot.svg @@ -21,167 +21,179 @@ - - + + - - - - - - - - - - - + + + + + + + + + + + + + - - + + - - - - - - - - - - - + + + + + + + + + + + + + - - + + - - - - - - - - - + + + + + + + + + + + - - + + - - - - - - - - - - - + + + + + + + + + + + + + - - + + - - - - - - - - - - - + + + + + + + + + + + + + - - + + - - - - - - - - - + + + + + + + + + + + - - + + - - -mean + + +mean - - + + - - -sd + + +sd - - + + - - -cjs_dist + + +cjs_dist - - + + - - -mu + + +mu - - + + - - -sigma + + +sigma @@ -202,94 +214,90 @@ 0.8 1 1.25 - - - -0.8 -1 -1.25 - - - -0.8 -1 -1.25 - - - -0.8 -1 -1.25 -0.00 -0.05 -0.10 -0.15 -0.20 - - - - - -0.00 -0.05 -0.10 -0.15 -0.20 - - - - - -0.20 -0.25 -0.30 -0.35 - - - - -0.15 -0.20 -0.25 - - - -9.35 -9.40 -9.45 -9.50 -9.55 -9.60 - - - - - - -0.85 -0.90 -0.95 - - - + + + +0.8 +1 +1.25 + + + +0.8 +1 +1.25 + + + +0.8 +1 +1.25 +0.00 +0.05 +0.10 +0.15 +0.20 + + + + + +0.00 +0.05 +0.10 +0.15 +0.20 + + + + + +0.20 +0.25 +0.30 +0.35 + + + + +0.15 +0.20 +0.25 + + + +9.35 +9.40 +9.45 +9.50 +9.55 +9.60 + + + + + + +0.85 +0.90 +0.95 + + + Power-scaling alpha - - - - -+/-2MCSE - -component - - - - -prior -likelihood -Posterior quantities depending on amount of power-scaling (alpha). -Horizontal lines indicate low sensitivity. -Steeper lines indicate greater sensitivity. -Estimates with high Pareto k (highlighted) may be inaccurate. -Power-scaling sensitivity + +Component + + + + +prior +likelihood + + + + ++/-2MCSE +Normal model quantities plot diff --git a/tests/testthat/test_component_name.R b/tests/testthat/test_component_name.R index c21a5e49..745bd894 100644 --- a/tests/testthat/test_component_name.R +++ b/tests/testthat/test_component_name.R @@ -1,26 +1,73 @@ ex <- example_powerscale_model()$draws ex_renamed <- posterior::rename_variables( - ex, - log_prior = lprior, - log_prior_sigma = lprior_sigma, - log_prior_mu = lprior_mu + ex, + ll = log_lik, + log_prior = lprior, + log_prior_sigma = lprior_sigma, + log_prior_mu = lprior_mu ) psd <- create_priorsense_data(ex) -psd_r <- create_priorsense_data(ex_renamed, log_prior_name = "log_prior") +psd_r <- create_priorsense_data( + ex_renamed, + log_prior_name = "log_prior", + log_lik_name = "ll" +) testthat::expect_error(powerscale_sensitivity(ex, log_lik_name = "ll")) testthat::expect_error(powerscale_sensitivity(ex_renamed)) testthat::expect_equal( - powerscale_sensitivity(psd_r), - powerscale_sensitivity(psd) + powerscale_sensitivity(psd_r), + powerscale_sensitivity(psd) +) + +testthat::expect_equal( + powerscale_sensitivity(ex), + powerscale_sensitivity( + ex_renamed, + log_prior_name = "log_prior", + log_lik_name = "ll" + ) +) + + +ex_new_var <- posterior::mutate_variables( + ex, + log_prior = lprior, + log_lik1 = `log_lik[1]`, + log_lik2 = `log_lik[1]`, +) + +testthat::expect_equal( + powerscale_sensitivity( + ex_new_var, + log_lik_name = "log_lik1", + variable = "mu" + ), + powerscale_sensitivity( + ex_new_var, + log_lik_name = "log_lik2", + variable = "mu" + ) ) testthat::expect_equal( - powerscale_sensitivity(ex), - powerscale_sensitivity(ex_renamed, log_prior_name = "log_prior") + powerscale_sensitivity( + ex_new_var, + log_lik_name = "log_lik", + likelihood_selection = "1", + separator = "", + variable = "mu" + )$likelihood, + powerscale_sensitivity( + ex_new_var, + log_lik_name = "log_lik", + likelihood_selection = "2", + separator = "", + variable = "mu" + )$likelihood ) diff --git a/tests/testthat/test_plots.R b/tests/testthat/test_plots.R index 98c56cdd..03ce1132 100644 --- a/tests/testthat/test_plots.R +++ b/tests/testthat/test_plots.R @@ -3,121 +3,82 @@ eight_schools_example <- example_powerscale_model("eight_schools") ps <- powerscale_sequence(eight_schools_example$draws, length = 3) test_that("diagnostic plots give no errors", { - expect_error( - powerscale_plot_ecdf( - ps, - variable = c("mu", "tau") - ), - NA - ) - expect_error( - powerscale_plot_dens( - x = ps, - variable = c("mu", "tau") - ), - NA - ) - expect_error( - powerscale_plot_quantities( - ps, - variable = c("mu", "tau") - ), - NA - ) + expect_error( + powerscale_plot_ecdf( + ps, + variable = c("mu", "tau") + ), + NA + ) + expect_error( + powerscale_plot_dens( + x = ps, + variable = c("mu", "tau") + ), + NA + ) + expect_error( + powerscale_plot_quantities( + ps, + variable = c("mu", "tau") + ), + NA + ) }) test_that("plots contain expected data", { - psq <- powerscale_plot_quantities( - ps, - variable = c("mu"), - quantity = c("quantile", "mean"), - quantity_args = list(probs = c(0.1, 0.9)) - ) - expect_equal( - colnames(psq$data), - c( - "variable", - ".powerscale_alpha", - "pareto_k_threshold", - "pareto_k", - "component", - "quantity", - "value", - "id", - "pareto_k_value" + psq <- powerscale_plot_quantities( + ps, + variable = c("mu"), + quantity = c("quantile", "mean"), + quantity_args = list(probs = c(0.1, 0.9)) + ) + expect_equal( + colnames(psq$data), + c( + "variable", + ".powerscale_alpha", + "pareto_k_threshold", + "pareto_k", + "component", + "quantity", + "value", + "id", + "pareto_k_value" + ) ) - ) - expect_equal( - unique(psq$data$quantity), - c("q10", "q90", "mean", "cjs_dist") - ) + expect_equal( + unique(psq$data$quantity), + c("q10", "q90", "mean", "cjs_dist") + ) }) -test_that("help_text behaves as expected in plots", { - psq_title <- powerscale_plot_quantities( - ps, - variable = c("mu"), - help_text = TRUE - ) - - psq_notitle <- powerscale_plot_quantities( - ps, - variable = c("mu"), - help_text = FALSE - ) - - expect_false(is.null(psq_title$labels$title)) - expect_false(is.null(psq_title$labels$subtitle)) - - expect_null(psq_notitle$labels$title) - expect_null(psq_notitle$labels$subtitle) - - psecdf_title <- powerscale_plot_ecdf(ps, variable = "mu") - - psecdf_notitle <- powerscale_plot_ecdf(ps, variable = "mu", help_text = FALSE) - - expect_false(is.null(psecdf_title$labels$title)) - expect_false(is.null(psecdf_title$labels$subtitle)) - - expect_null(psecdf_notitle$labels$title) - expect_null(psecdf_notitle$labels$subtitle) - - psdens_title <- powerscale_plot_dens(ps, variable = "mu") - psdens_notitle <- powerscale_plot_dens(ps, variable = "mu", help_text = FALSE) - - expect_false(is.null(psdens_title$labels$title)) - expect_false(is.null(psdens_title$labels$subtitle)) - - expect_null(psdens_notitle$labels$title) - expect_null(psdens_notitle$labels$subtitle) -}) - test_that("pagination of plots works as expected", { - expect_length( - powerscale_plot_quantities( - ps, - variables_per_page = 1 - ), - 18 - ) + expect_length( + powerscale_plot_quantities( + ps, + variables_per_page = 1 + ), + 18 + ) - expect_length( - powerscale_plot_quantities( - ps, - variables_per_page = 2 - ), - 9 - ) + expect_length( + powerscale_plot_quantities( + ps, + variables_per_page = 2 + ), + 9 + ) - expect_length( - powerscale_plot_quantities( - ps, - variables_per_page = Inf - ), - 1 - ) + expect_length( + powerscale_plot_quantities( + ps, + variables_per_page = Inf + ), + 1 + ) }) @@ -125,17 +86,17 @@ test_that("pagination of plots works as expected", { ps_normal <- powerscale_sequence(example_powerscale_model()$draws) vdiffr::expect_doppelganger( - "Normal model density plot", - powerscale_plot_dens(ps_normal) + "Normal model density plot", + powerscale_plot_dens(ps_normal) ) vdiffr::expect_doppelganger( - "Normal model ecdf plot", - powerscale_plot_ecdf(ps_normal) + "Normal model ecdf plot", + powerscale_plot_ecdf(ps_normal) ) vdiffr::expect_doppelganger( - "Normal model quantities plot", - powerscale_plot_quantities(ps_normal) + "Normal model quantities plot", + powerscale_plot_quantities(ps_normal) ) diff --git a/vignettes/used_by.bib b/vignettes/papers.bib similarity index 100% rename from vignettes/used_by.bib rename to vignettes/papers.bib diff --git a/vignettes/used_by.csl b/vignettes/papers.csl similarity index 100% rename from vignettes/used_by.csl rename to vignettes/papers.csl diff --git a/vignettes/used_by.qmd b/vignettes/papers.qmd similarity index 57% rename from vignettes/used_by.qmd rename to vignettes/papers.qmd index c98101fd..c0630d0e 100644 --- a/vignettes/used_by.qmd +++ b/vignettes/papers.qmd @@ -1,7 +1,11 @@ --- title: "Published papers using priorsense" -bibliography: used_by.bib -csl: used_by.csl +vignette: > + %\VignetteIndexEntry{Published papers using priorsense} + %\VignetteEngine{quarto::html} + %\VignetteEncoding{UTF-8} +bibliography: papers.bib +csl: papers.csl nocite: | @* ---