easystats / bayestestR

:ghost: Utilities for analyzing Bayesian models and posterior distributions
https://easystats.github.io/bayestestR/
GNU General Public License v3.0
576 stars 55 forks source link

Is a directional p_rope possible? #635

Closed qdread closed 11 months ago

qdread commented 11 months ago

I know that bayesfactor_rope() has a direction argument to specify a left-sided vs. right-sided test. I do not see a similar argument for p_rope(). Is there any way to do a one-sided test, for instance if you have an a priori assumption the effect must be positive, using p_rope()?

mattansb commented 11 months ago

bayesfactor_rope() has a direction to allow for comparing different hypotheses in the posterior(/prior)space. p_rope() does not compare different hypotheses, it "simply" computes the posterior probability of the parameter being in the ROPE.

If you have a priori assumption the effect must be positive, this can be tested to be true with p_direction(), and the ROPE can be defined with a lower bound of -Inf.

library(rstanarm)
library(bayestestR)

mod <- stan_glm(mpg ~ qsec, data = mtcars)

p_direction(mod)
#> Probability of Direction 
#> 
#> Parameter   |     pd
#> --------------------
#> (Intercept) | 67.88%
#> qsec        | 99.20%

p_rope(mod, range = c(-Inf, 0.5))
#> Proportion of samples inside the ROPE [-0.10, 0.10]
#> 
#> 
#> Parameter   | p (ROPE)
#> ----------------------
#> (Intercept) |    0.696
#> qsec        |    0.062

Or together in describe_posterior():

describe_posterior(mod, rope_ci = 1, rope_range = c(-Inf, 0.5))
#> Summary of Posterior Distribution 
#> 
#> Parameter   | Median |          95% CI |     pd |         ROPE | % in ROPE |  Rhat |     ESS
#> --------------------------------------------------------------------------------------------
#> (Intercept) |  -4.91 | [-25.37, 15.93] | 67.88% | [-Inf, 0.50] |    69.58% | 1.001 | 3213.00
#> qsec        |   1.39 | [  0.25,  2.55] | 99.20% | [-Inf, 0.50] |     6.15% | 1.001 | 3158.00

This would be reported as:

As hypothesized, the slope of effect of qsec had a probability of 0.992 of being positive, and a probability of 0.062 of being less than 0.5.


Alternatively, this a-priori information can be incorporated into the model in one of few ways. Here is an example, but see that link for more info and details.

library(posterior)

posteriors <- as_draws_rvars(mod)

iS <- posteriors$qsec > 0

posteriors_S <- lapply(posteriors, function(p) p[iS]) |> do.call(what = "c")

describe_posterior(posteriors_S, rope_ci = 1, rope_range = c(0, 0.5))
#> Summary of Posterior Distribution
#> 
#> Parameter      | Median |          95% CI |     pd |         ROPE | % in ROPE
#> -----------------------------------------------------------------------------
#> x[(Intercept)] |  -5.01 | [-25.39, 14.61] | 68.42% | [0.00, 0.50] |     1.71%
#> x[qsec]        |   1.40 | [  0.32,  2.55] |   100% | [0.00, 0.50] |     5.39%
#> x[sigma]       |   5.65 | [  4.44,  7.49] |   100% | [0.00, 0.50] |        0%

(Note that the probability of direction is trivially 100% now). This would be reported as:

When applying an order constraint allowing only for positive (>0) slopes of qsec, the slope had a probability 0.054 of being less than 0.5.