• Minor issues formatting related to bat sticks here in details

Now that the initial release rush is over, I've checked your plotting methods, which have a very ingenious way of retrieving the original data.

However, I feel like this would need a more systematic implementation (so that the data retrieval attempt works for all functions, p_direction, hdi, etc.).

Do you think adding the name of the object as an attribute would work in most cases? (I believe it won't work in cases such as hdi(rnorm(1000)). Maybe it would be simpler (not from an efficiency point of view obviously) to store the original data in the objects (especially when these are non-models)?

Would be nice to unify the printing output up to the standard of ROPE/equivalence test.

• p_direction: dataframe-like output needed
• p_map: everything needed

Do we need insight listed under "remotes" when it's also listed under "imports"?

Maybe we manage it to make functions like equivalence_test() also work for multivariate response models (like stanmvreg), however, this seems more like a long-term goal to me.

Once I'm done with the remaining stuff for insight and it's submitted to CRAN, I would try to write some print / plot methods for some of the bayestestR functions.

I think especially equivalence_test() would benefit from a nice print-method. What would you say?

equivalence_test() returns results for the complete posterior distribiution. However, if you fit your model with brms and use sample_prior = TRUE, then you have parameters for the priors as well, which are most likely of no interest for functions like equivalence_test().

Hence, remove all paramaters starting with prior_.

I can work on this issue...

Hum, I just realised that the bayes_indices.csv dataset, used in the vignettes, is loaded with the package. Which is unintended and unwanted from a weight perspective. Moreover, the vignettes work by downloading the dataset from github, so no need for it to be stored in the package.

From here:

• If you want to store parsed data, but not make it available to the user, put it in R/sysdata.rda. This is the best place to put data that your functions need.
• If you want to store raw data, put it in inst/extdata.

I would expect the p-value to have a range from 0-1, not 0 to 100.

it would be good if we managed to improve testing coverage to > 90% before submitting to JOSS or somewhere else, which might necessitate to try making testing work for brms models.

See email.

I ran a Bayesian multilevel model in brms and tried to plot the results using the equivalence_test() function:

library(bayestestR)
x <- equivalence_test(model)
plot(x)

and got the following error:
Error in FUN(X[[i]], ...) : object 'grp' not found

I then used the function equi_test from the sjstats package and it worked:

library(sjstats)
equi_test(x, rope = c(-0.05, 0.05), plot = TRUE, out = 'plot')

which did not give me the error. I wonder if I am doing something wrong or if this might be a bug?

the individual y values of the density function are not commonly used (as their meaning is in relationship with portions of the function from which to estimate probabilities).
Nevertheless, it seems like the density value of the MAP is sensitive to the general precision (I.e., range and spread) of the density function. The higher it is, and the more "precise" is the posterior. Hence, I wonder how this index would behave in relationship to other markers such as rope percentage, pd and such. Could be interesting to include it in a revised simulation study.

• I'll check the new hdi() for some models tonight, and if it works, we may extend the arguments effects and component to rope() as well.
• The print()-method for equivalence_test() for models should be revised, so it does not print redundant information, but includes useful information (HDI for each parameter) instead.
• The print()-method for equivalence_test() should work for multiple CI-values.
• Then, the ci()-method can be targeted, and this PR is "complete".

Originally posted by @strengejacke in #29 (comment)

https://easystats.github.io/bayestestR/reference/rope_range.html

Could be a bit more text, I think.

@strengejacke I am having issues for ROPE support (see bottom of README for example). It seems the rope object based on models does contain all the information regarding HDI_area. It has (in attributes) only one value for it, instead of several for parameters and levels of ci. Is it possible to address it here, or modify the thing in bayestestR?

It seems that the extension of ci and equivalence_test in parameters works, although documenting their methods separately is not straightforward. Users will probably be shown bayestestR's documentation instead of the one specific to frequentist models.

Thus, I think we will have to add a link to the parameters' website "for frequentist models, you can find the documentation here" and/or directly a small paragraph in details documenting the method for frequentist models.

Altough for now, there is only the "kruschke2018" method, I might be interesting to prepare the terrain for an easy switch between methods of default values selection, in a similar way to the interpret_* functions in report. As with every rules of thumb, I can foresee a debate surrounding the appropriate values in the literature, so we might as well support them here :)

Related to #5 #6

Currently, all functions in the package are "model agnostic", i.e., they don't work on models but directly on the posterior vectors.

These are fundamental operations that would be used in functions dealing with Bayesian models.

In fact, I first considered dissociating the "model-related" part with the "posterior-related" part, but writing about this, I am not sure it makes much sense. First, the scope (and name) of the package suggests its ability to deal with models. And second, it would nice to have the same function names for posteriors and models (hdi(rnorm(1000)) and hdi(stanmodel)).

The other possibility would be to move the model-targetted functions to the "parameters" (such as sjstats::tidy_stan) and "model performance" (e..g. for sjstats::msce) packages.

What is your feeling about this?

We already have a ci() method currently defined in parameters. However, parameters does depend on bayestestR. Thus, we can try defining the method in bayestestR and expanding it to frequentist models in parameters. We just have to find a way to put the links of the online documentation of each method in the docstring of bayestestR's definition.

See strengejacke/sjstats#72

I fit an errors-in-variables model in brms (i.e., with the term me(age, age_se)) and sjstats parameter extraction misfired (it shows a bunch of parameters in tab_model and equi_test that are not shown when calling summary.brmsfit. It would be cool, if you permitted a pars argument like many brms and bayesplot functions to allow customisation.

This would not just be useful for misfiring like above, but also to trim the tested parameters in equi_test if not all are of substantial interest.

Would you consider adding a function for a Savage-Dickey density ratio Bayes factor?

Something like this, perhaps?

savage_dickey_BF <- function(xpost,xprior,null = 0,direction = 0){
  if (!require(logspline)) {
    stop('Need logspline to run')
  }
  f_post <- suppressWarnings(logspline(xpost))
  f_prior <- suppressWarnings(logspline(xprior))
  
  d_post <- dlogspline(null,f_post)
  d_prior <- dlogspline(null,f_prior)
  
  norm_post <- norm_prior <- 1
  if (direction < 0) {
    norm_post <- plogspline(null,f_post)
    norm_prior <- plogspline(null,f_prior)
  } else if (direction > 0) {
    norm_post <- 1-plogspline(null,f_post)
    norm_prior <- 1-plogspline(null,f_prior)
  }
  
  (d_post/norm_post)/(d_prior/norm_prior)
}

Move here from strengejacke/sjstats#68

@IndrajeetPatil due to the re-organization of some packages, all Bayes-related stuff will now be implemented in the bayestestR-package. So I'm moving your issue over here.

I can't check the package because vignettes can't be built. The vignette index is missing. Any ideas? Or is it because I'm testing the dev branch?

https://easystats.github.io/bayestestR/reference/rope.html

References to Kruschke need to be added.

My impression is that bayestestR is not far away from its initial CRAN submission. Would you agree?
Let's put together a to-do list to check what needs to be done before submission:

• implement ci() (related to #30) (DL)
• maybe implement plot() for equivalence_test()? I really would like to have this feature, because I think that (visual) equivalence-testing is rising more and more awareness, and I've been already contacted two or three times regarding this issue (related to #15 and #34) (DL)
• check test-coverage (DL / DM)
• compile vignettes (to see if it works and they look good?) (related to #32) (DM)
• remove Remotes field from DESCRIPTION (DM)
• related to the previous point, do we need to fix the test-files, i.e. to make them not run on CRAN? (also related to #19) (DL / DM)
• check docs / helpfiles (related to #33 ) (DL / DM)
• allow specification of pars (related to #37) (DL)
• final check / website

What else? I also made some suggestions who takes over which tasks, feel free to add / remove yourself from the list. ;-)

The covered percentage in ROPE differs between those two functions. We need to check which of these both are probanly wrong, or if the computational basis just differs.

Thanks,

we see:
Size of tarball: 5048980 bytes

Not more than 5 MB for a CRAN package, please.

If there are any references describing the methods in your package,
please add these in the Description field of your DESCRIPTION file in
the form
authors (year) doi:...
authors (year) arXiv:...
authors (year, ISBN:...)
with no space after 'doi:', 'arXiv:' and angle brackets for auto-linking.

Please fix and resubmit.

model <- rstanarm::stan_glm(mpg ~ wt + cyl, data = mtcars)
attributes(rope(model))$HDI_area
#> $`(Intercept)`
#>     CI  HDI_low HDI_high
#> 1 0.90 36.60966 42.66055
#> 2 0.95 36.12163 43.45867
#> 
#> $wt
#>     CI   HDI_low  HDI_high
#> 1 0.90 -4.432082 -1.776227
#> 2 0.95 -4.668768 -1.481225
#> 
#> $cyl
#>     CI   HDI_low   HDI_high
#> 1 0.90 -2.254774 -0.8010249
#> 2 0.95 -2.396010 -0.6578996

Since we return one or multiple dataframes for models (which I like), could we also return such dataframe for the cases where a vector is passed to ROPE?

library(bayestestR)
attributes(rope(x = rnorm(1000, 1, 1)))$HDI_area
#> [[1]]
#> [1] -0.5071979  2.7443272
attributes(rope(x = rnorm(1000, 1, 1), ci = c(.90, .95)))$HDI_area
#> [[1]]
#> [1] -0.8034466  2.5204936
#> 
#> [[2]]
#> [1] -0.8451218  2.9656530

So the output would like:

attributes(rope(x = rnorm(1000, 1, 1), ci = c(.90, .95)))$HDI_area
#>     CI  HDI_low HDI_high
#> 1 0.90 36.60966 42.66055
#> 2 0.95 36.12163 43.45867

Sorry to bother with this, I would gladly change it myself but I am afraid of breaking everything since the attribute code is still not perfectly clear to me (I am studying it tho :)

Originally posted by @DominiqueMakowski in #69 (comment)

I found a strange inconsistency for hdi() with brms-models, but I'm not sure if this is only happening for brms-models or not? I'm not sure in how far your implementation of HDI differs from mine (which is adapted from Kruscke 2015).

stanreg = identical
brms = differences

m1 <- insight::download_model("stanreg_lm_1")

as.data.frame(bayestestR::hdi(m1, ci = .95))
#>     Parameter CI    CI_low   CI_high
#> 1 (Intercept) 95 33.641431 41.069859
#> 2          wt 95 -6.418331 -4.236012

as.data.frame(sjstats::hdi(m1, prob = .95))
#>          term   hdi.low  hdi.high
#> 1 (Intercept) 33.641431 41.069859
#> 2          wt -6.418331 -4.236012

set.seed(123)
library(brms)
b4 <- brms::brm(mpg ~ wt + cyl + gear, data = mtcars, iter = 500, chains = 1)

as.data.frame(bayestestR::hdi(b4, ci = .95))
#>     Parameter CI    CI_low    CI_high
#> 1 b_Intercept 95 34.996964 50.9281495
#> 2        b_wt 95 -5.339557 -1.9348616
#> 3       b_cyl 95 -2.380895 -0.6869611
#> 4      b_gear 95 -1.790944  1.0685669

as.data.frame(sjstats::hdi(b4, prob = .95))
#>          term   hdi.low   hdi.high
#> 1 b_Intercept 34.549999 50.8355624
#> 2        b_wt -5.373982 -1.9348616
#> 3       b_cyl -2.295808 -0.5923283
#> 4      b_gear -1.845044  1.0685669

^{Created on 2019-03-14 by the reprex package (v0.2.1)}

I think it would be nice to have methods for stanfit or brms-objects as well, so the functions not only expect a single vector, but may also cope with model objects.

@strengejacke I am still unsure about the name of the equivalence test, currently called rope_test() (as it is directly related to the rope() function. I know that you named it equi_test in sjstats (although to be used on models instead of vectors of posteriors). What is your opinion?

@DominiqueMakowski This commit requires some attention! Since map_estimate() returns different things, depending on the input and the arguments, I tried to add the section for the return-value to the docs. However, I'm not sure what density does. so can you check this please, and add a short description to the @return-part of the docs?

I do have this error when testing hdi() for brms and stanreg models:

> bayestestR::hdi(x)
 Hide Traceback
 
 Rerun with Debug
 Error in data[[variable]] : subscript out of bounds 
5.
which(data[[variable]] == value) at utils.R#22
4.
.select_rows(dat, "Group", "fixed") at utils_hdi_ci.R#79
3.
.compute_interval_stanreg(x, ci, effects, parameters, verbose, 
    fun = "hdi") at hdi.R#72
2.
hdi.stanreg(x) at hdi.R#54
1.
bayestestR::hdi(x) 

I think I have the other packages correctly updated (insight)
do you have it too?

I think we should also add a function that returns the credible intervals (not HDI) as well, for convenience reasons, where the output meets our needs.

In case of transformation, eg having odds ratios, you can also exp the CI, but not the HDI.

See https://discourse.mc-stan.org/t/exponentiation-or-transformation-of-point-estimates/7848

Especially for models, the function has some redundant output (CI and rope bounds are identical for each parameter). Maybe we can also add the HDI? So users see the HDI and the rope, and can see which part lies inside/outside the rope? We should probably think about a print()-method, to have some user-friendly output.

library(bayestestR)
library(insight)
m <- insight::download_model("stanreg_glm_5")

class(hdi(m))
#> [1] "hdi"        "hdi"        "hdi"        "data.frame"

class(ci(m, ci = c(.5, .8)))
#> [1] "ci"         "data.frame"

class(rope(m))
#> [1] "rope"       "data.frame"

^{Created on 2019-04-18 by the reprex package (v0.2.1)}

Have found this only for hdi() yet. We should make sure that class-attributes simply are unique().

Would you say that it makes sense to "copy" this vignette from sjstats:
https://strengejacke.github.io/sjstats/articles/bayesian-statistics.html

Or maybe revise it slightly?

Turns out vignettes files can't start with numbers, which was added for dealing with their order in the website. Must specify it in pkgdown parameters:

• example

It would be interesting to include some methods to check the parameters sensitivity to different priors specification. I am not sure if there is any "standard default" method for doing so.

The ROPE-based p-value represents the maximum Credible Interval (HDI) that does not contain (positive values) or is entirely contained (negative values) in the negligible values space defined by the ROPE. A ROPE-based p of 97

1. it's not only positive, but also only negative values.
2. the description ends abrupt.

@strengejacke I would be grateful if you could give me your input and opinion regarding this package and its functions 😺

Could you please tell me how to plot the figures that showed in the example, these are very nice, thanks.

Related to the fact that ROPE bounds are selected automatically for models with the 0.1 * SD(response) formula. How to adapt this to binomial models?

Log odds can be converted to standardized differences (see formula), so maybe through something involving it?

Is the absence of variable expected?

What do you think about preparing (and submitting) a paper for JOSS?

https://joss.theoj.org/

It would be a very short paper, mostly similar to a readme-file, but would have the advantage of a "nicer" citation for the package. The journal has an ISSN, so it's being considered as "citable".

I have revised the docs slightly in this commit: 856cd25

1. I changed Value or vector of HDI probability to alue or vector of probability of the interval as we inherits the docs from hdi() also in ci(), so this is more general.

2. In ci(). I removed For Bayesian Credible Intervals, CIs are often computed by the hdi method. I'm not sure if HDI's are really more often computed than equal-tailed credible intervals, so unless we don't have a proof for this statement, we may better remove it.

This issue to plan the documentation, to avoid redundancy and making the documentation useful, clear and easily findable

• README
  • All main functions with their docstrings?
• Vignettes
  • Get started
    • Why use the Bayesian framework?
    • A very simple linear model
  • Example: Bayesian GLMs
  • Example: Bayesian mixed models
  • In depth: Point estimates comparison
  • In depth: Indices of effect existence comparison
  • Reporting guidelines

@strengejacke How would do you see it?

Altough I feel like analyses examples should be in this repo, the home for Bayesian-related stuff, at the same time it would make then sense to add here the other aspects related to Bayesian models present in easystats (such as model's performance, contrasts estimates, report (and visreport?) etc.

Ressources (for inspiration):

• https://strengejacke.github.io/sjstats/articles/bayesian-statistics.html

I would have thought that p_rope() is similar or even equivalent to sjstats::equi_test(). In this model, 28.6% of the values of the posterior from the parameter am are inside the rope of (-.1, .1). p_rope() returns a p-value of 8.5. How do these values relate to each other?

library(sjstats)
library(bayestestR)
library(rstanarm)
set.seed(123)
m <- stan_glm(mpg ~ wt + am, data = mtcars, chains = 1)

equi_test(m)
#> 
#> # Test for Practical Equivalence of Model Predictors
#> 
#>   Effect Size: 0.10
#>          ROPE: [-0.60 0.60]
#>       Samples: 1000
#> 
#>                         H0 %inROPE      HDI(95%)
#>  (Intercept) (*)    reject    0.00 [30.54 42.76]
#>  wt (*)             reject    0.00 [-6.96 -3.79]
#>  am (*)          undecided   28.60 [-2.92  2.98]
#>  sigma (*)          reject    0.00 [ 2.47  4.08]
#> 
#> (*) the number of effective samples may be insufficient for some parameters

p_rope(as.data.frame(m)$am)
#> [1] 8.5

Enable the ROPE parameter to "default", which will be -0.1 0.1 for rope.numeric and 0.1 * SD(response) for models... But what in the case of logistic models? What's the best way to return the ROPE bounds?