Unexpected behavior when using `quantities = “quantile”` at `powerscale_plot_quantities`

[lcgodoy]

Hello everyone,

Thanks for making such an amazing package available.

I use priorsense to run a sensitivity analysis for a consulting project. For this specific project, we were interested in checking the changes in the confidence intervals. However, I get an error when I try to run powerscale_plot_quantities with quantities = “quantile”. See the reprex below.

library(priorsense)

normal_model <- example_powerscale_model("univariate_normal")

fit <- rstan::stan(
  model_code = normal_model$model_code,
  data = normal_model$data,
  refresh = FALSE,
  seed = 1234
  )

pss <- powerscale_sequence(fit)
#> Loading required namespace: testthat

powerscale_plot_quantities(
  pss,
  quantities = c("quantile"),
  probs = c(.05, .95),
  div_measure = "cjs_dist",
  variables = c("mu", "sigma")
)
#> Error in rbind(deparse.level, ...): numbers of columns of arguments do not match

I did some debugging, and there are (at least) two issues here. Firstly, powerscale_plot_quantities ignores probs. This leads to the first error, the quantile_weighted function (used for the power scaled priors) has probs = c(.05, .95) as its default, while the stats::quantile function uses probs = seq(0, 1, 0.25). Next, if we rewrite the quantile_weight setting probs = seq(0, 1, 0.25) as the default we still get an error. The error is because the function's output is renamed differently than the output of the stats::quantile function.

In summary, if I write

quantile_weighted <- function(x, weights,
                              probs = seq(0, 1, 0.25),
                              type = "7", ...) {
  if (type == "7") {
    cdf_fun <- function(n, p) {
      return(function(cdf_probs) {
        h <- p * (n - 1) + 1
        u <- pmax((h - 1) / n, pmin(h / n, cdf_probs))
        u * n - h + 1
      })
    }
  }
  else if (type == "hd") {
    cdf_fun <- function(n, p) {
      return(function(cdf_probs) {
        stats::pbeta(cdf_probs, (n + 1) * p, (n + 1) *
                                             (1 - p))
      })
    }
  }
  quants <- priorsense:::.quantile_weighted(x = x, weights = weights,
                                            probs = probs,
                                            cdf_fun = cdf_fun)
  names(quants) <- paste0(probs * 100, "%")
  return(quants)
}

The function call below works almost with no issues. However, it displays the results for probs = seq(0, 1, 0.25), ignoring the input probs = c(.05, .95).

powerscale_plot_quantities(
  pss,
  quantities = c("quantile"),
  probs = c(.05, .95),
  div_measure = "cjs_dist",
  variables = c("mu", "sigma")
)

A possible solution would be adding a quantity_args parameter to deal with arguments from summary functions in the same spirit as what you have done with measure_args.

I’d be happy to help with a pull request. Is there a contribution guide?

Also, I’m looking forward to the release of the separate_scaling features.

Best, Lucas.

Session info

``` r sessioninfo::session_info() #> ─ Session info ─────────────────────────────────────────────────────────────── #> setting value #> version R version 4.3.1 (2023-06-16) #> os macOS Sonoma 14.2.1 #> system aarch64, darwin20 #> ui X11 #> language (EN) #> collate en_US.UTF-8 #> ctype en_US.UTF-8 #> tz America/Costa_Rica #> date 2024-01-16 #> pandoc 3.1.6.1 @ /opt/homebrew/bin/ (via rmarkdown) #> #> ─ Packages ─────────────────────────────────────────────────────────────────── #> package * version date (UTC) lib source #> abind 1.4-5 2016-07-21 [1] CRAN (R 4.3.0) #> backports 1.4.1 2021-12-13 [1] CRAN (R 4.3.0) #> brio 1.1.3 2021-11-30 [1] CRAN (R 4.3.0) #> callr 3.7.3 2022-11-02 [1] CRAN (R 4.3.0) #> checkmate 2.3.1 2023-12-04 [1] CRAN (R 4.3.1) #> cli 3.6.2 2023-12-11 [1] CRAN (R 4.3.1) #> codetools 0.2-19 2023-02-01 [1] CRAN (R 4.3.1) #> colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.3.0) #> crayon 1.5.2 2022-09-29 [1] CRAN (R 4.3.0) #> curl 5.0.2 2023-08-14 [1] CRAN (R 4.3.0) #> digest 0.6.34 2024-01-11 [1] CRAN (R 4.3.1) #> distributional 0.3.2 2023-03-22 [1] CRAN (R 4.3.0) #> dplyr 1.1.4 2023-11-17 [1] CRAN (R 4.3.1) #> evaluate 0.21 2023-05-05 [1] CRAN (R 4.3.0) #> fansi 1.0.6 2023-12-08 [1] CRAN (R 4.3.1) #> farver 2.1.1 2022-07-06 [1] CRAN (R 4.3.0) #> fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.3.0) #> fs 1.6.3 2023-07-20 [1] CRAN (R 4.3.0) #> generics 0.1.3 2022-07-05 [1] CRAN (R 4.3.0) #> ggplot2 3.4.4 2023-10-12 [1] CRAN (R 4.3.1) #> glue 1.7.0 2024-01-09 [1] CRAN (R 4.3.1) #> gridExtra 2.3 2017-09-09 [1] CRAN (R 4.3.0) #> gtable 0.3.4 2023-08-21 [1] CRAN (R 4.3.0) #> htmltools 0.5.6 2023-08-10 [1] CRAN (R 4.3.0) #> inline 0.3.19 2021-05-31 [1] CRAN (R 4.3.0) #> jsonlite 1.8.7 2023-06-29 [1] CRAN (R 4.3.0) #> knitr 1.43 2023-05-25 [1] CRAN (R 4.3.0) #> lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.3.1) #> loo 2.6.0 2023-03-31 [1] CRAN (R 4.3.0) #> magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.3.0) #> matrixStats 1.2.0 2023-12-11 [1] CRAN (R 4.3.1) #> munsell 0.5.0 2018-06-12 [1] CRAN (R 4.3.0) #> pillar 1.9.0 2023-03-22 [1] CRAN (R 4.3.0) #> pkgbuild 1.4.2 2023-06-26 [1] CRAN (R 4.3.0) #> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.3.0) #> posterior 1.5.0 2023-10-31 [1] CRAN (R 4.3.1) #> prettyunits 1.2.0 2023-09-24 [1] CRAN (R 4.3.1) #> priorsense * 0.0.0.9000 2024-01-15 [1] Github (n-kall/priorsense@f6c6bca) #> processx 3.8.2 2023-06-30 [1] CRAN (R 4.3.0) #> ps 1.7.5 2023-04-18 [1] CRAN (R 4.3.0) #> QuickJSR 1.0.8 2023-11-24 [1] CRAN (R 4.3.1) #> R6 2.5.1 2021-08-19 [1] CRAN (R 4.3.0) #> Rcpp 1.0.11 2023-07-06 [1] CRAN (R 4.3.1) #> RcppParallel 5.1.7 2023-02-27 [1] CRAN (R 4.3.1) #> reprex 2.1.0 2024-01-11 [1] CRAN (R 4.3.1) #> rlang 1.1.3 2024-01-10 [1] CRAN (R 4.3.1) #> rmarkdown 2.24 2023-08-14 [1] CRAN (R 4.3.0) #> rstan 2.32.3 2023-10-15 [1] CRAN (R 4.3.1) #> scales 1.3.0 2023-11-28 [1] CRAN (R 4.3.1) #> sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.3.0) #> StanHeaders 2.26.28 2023-09-07 [1] CRAN (R 4.3.1) #> tensorA 0.36.2.1 2023-12-13 [1] CRAN (R 4.3.1) #> testthat 3.1.10 2023-07-06 [1] CRAN (R 4.3.0) #> tibble 3.2.1 2023-03-20 [1] CRAN (R 4.3.0) #> tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.3.0) #> utf8 1.2.4 2023-10-22 [1] CRAN (R 4.3.1) #> V8 4.3.3 2023-07-18 [1] CRAN (R 4.3.0) #> vctrs 0.6.5 2023-12-01 [1] CRAN (R 4.3.1) #> withr 2.5.2 2023-10-30 [1] CRAN (R 4.3.1) #> xfun 0.40 2023-08-09 [1] CRAN (R 4.3.0) #> yaml 2.3.7 2023-01-23 [1] CRAN (R 4.3.0) #> #> [1] /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library #> #> ────────────────────────────────────────────────────────────────────────────── ```

[n-kall]

Thanks for reporting this. I think the issue arose from quantile_weighted trying to match posterior::quantile2 rather than stats::quantile. So using "quantile2" would have worked (but still would only give q5 and q95). I have now adjusted things such that "quantile" should also work, and added a quantity_args as you suggested which can be used for e.g. probs = c(0.1, 0.9). Let me know if there are still issues

[lcgodoy]

It worked great, thanks a lot!