Flexible ROPE values for describe_posterior

[strengejacke]

Fixes #643

[strengejacke]

I think this is ready to review, if you like, @mattansb . But everything should be ok, I also added some tests. (test-coverage is really bad in this package... 😬 )

[strengejacke]

Ok, let me fix errors first

[strengejacke]

Ok @mattansb, I think you can review. Here are some examples for relevant changes / new features added in this PR, addressing #643

I added this feature

• rope()
• describe_posterior()
• equivalence_test()
• p_significance()

library(bayestestR)
model <- rstanarm::stan_glm(
  mpg ~ wt + cyl,
  data = mtcars,
  chains = 2, refresh = 0
)

p_significance(model)
#> Practical Significance (threshold: 0.60)
#> 
#> Parameter   |   ps
#> ------------------
#> (Intercept) | 1.00
#> wt          | 1.00
#> cyl         | 0.98

# multiple thresholds - asymmetric, symmetric, default
p_significance(model, threshold = list(c(-10, 5), 0.2, "default"))
#> Practical Significance
#> 
#> Parameter   |   ps |           ROPE
#> -----------------------------------
#> (Intercept) | 1.00 | [-10.00, 5.00]
#> wt          | 1.00 | [ -0.20, 0.20]
#> cyl         | 0.98 | [ -0.60, 0.60]

# names thresholds
p_significance(model, threshold = list(wt = 0.2, `(Intercept)` = c(-10, 5)))
#> Practical Significance
#> 
#> Parameter   |   ps |           ROPE
#> -----------------------------------
#> (Intercept) | 1.00 | [-10.00, 5.00]
#> wt          | 1.00 | [ -0.20, 0.20]
#> cyl         | 0.98 | [ -0.60, 0.60]

equivalence_test(model)
#> Possible multicollinearity between cyl and wt (r = 0.8). This might lead
#>   to inappropriate results. See 'Details' in '?equivalence_test'.
#> # Test for Practical Equivalence
#> 
#>   ROPE: [-0.60 0.60]
#> 
#> Parameter   |       H0 | inside ROPE |        95% HDI
#> -----------------------------------------------------
#> (Intercept) | Rejected |      0.00 % | [36.19, 43.03]
#> wt          | Rejected |      0.00 % | [-4.72, -1.66]
#> cyl         | Rejected |      0.00 % | [-2.35, -0.64]

# multiple ROPE ranges - asymmetric, symmetric, default
equivalence_test(model, range = list(c(10, 40), c(-5, -4), "default"))
#> Possible multicollinearity between cyl and wt (r = 0.8). This might lead
#>   to inappropriate results. See 'Details' in '?equivalence_test'.
#> # Test for Practical Equivalence
#> 
#> Parameter   |        H0 | inside ROPE |        95% HDI |           ROPE
#> -----------------------------------------------------------------------
#> (Intercept) | Undecided |     58.05 % | [36.19, 43.03] | [10.00, 40.00]
#> wt          | Undecided |     14.26 % | [-4.72, -1.66] | [-5.00, -4.00]
#> cyl         |  Rejected |      0.00 % | [-2.35, -0.64] |  [-0.10, 0.10]

# named
equivalence_test(model, range = list(wt = c(-5, -4)))
#> Possible multicollinearity between cyl and wt (r = 0.8). This might lead
#>   to inappropriate results. See 'Details' in '?equivalence_test'.
#> # Test for Practical Equivalence
#> 
#> Parameter   |        H0 | inside ROPE |        95% HDI |           ROPE
#> -----------------------------------------------------------------------
#> (Intercept) |  Rejected |      0.00 % | [36.19, 43.03] |  [-0.10, 0.10]
#> wt          | Undecided |     14.26 % | [-4.72, -1.66] | [-5.00, -4.00]
#> cyl         |  Rejected |      0.00 % | [-2.35, -0.64] |  [-0.10, 0.10]

describe_posterior(model)
#> Summary of Posterior Distribution
#> 
#> Parameter   | Median |         95% CI |     pd |          ROPE | % in ROPE |  Rhat |     ESS
#> --------------------------------------------------------------------------------------------
#> (Intercept) |  39.68 | [36.19, 43.03] |   100% | [-0.60, 0.60] |        0% | 1.000 | 2234.00
#> wt          |  -3.17 | [-4.72, -1.66] |   100% | [-0.60, 0.60] |        0% | 1.000 | 1008.00
#> cyl         |  -1.51 | [-2.35, -0.64] | 99.95% | [-0.60, 0.60] |        0% | 0.999 | 1050.00

describe_posterior(model, rope_range = list(c(10, 40), c(-5, -4), "default"))
#> Summary of Posterior Distribution
#> 
#> Parameter   | Median |         95% CI |     pd |           ROPE | % in ROPE |  Rhat |     ESS
#> ---------------------------------------------------------------------------------------------
#> (Intercept) |  39.68 | [36.19, 43.03] |   100% | [10.00, 40.00] |    58.05% | 1.000 | 2234.00
#> wt          |  -3.17 | [-4.72, -1.66] |   100% | [-5.00, -4.00] |    14.26% | 1.000 | 1008.00
#> cyl         |  -1.51 | [-2.35, -0.64] | 99.95% | [-0.10,  0.10] |        0% | 0.999 | 1050.00

describe_posterior(model, rope_range = list(wt = c(-5, -4)))
#> Summary of Posterior Distribution
#> 
#> Parameter   | Median |         95% CI |     pd |           ROPE | % in ROPE |  Rhat |     ESS
#> ---------------------------------------------------------------------------------------------
#> (Intercept) |  39.68 | [36.19, 43.03] |   100% | [-0.10,  0.10] |        0% | 1.000 | 2234.00
#> wt          |  -3.17 | [-4.72, -1.66] |   100% | [-5.00, -4.00] |    14.26% | 1.000 | 1008.00
#> cyl         |  -1.51 | [-2.35, -0.64] | 99.95% | [-0.10,  0.10] |        0% | 0.999 | 1050.00

# errors - wrong length
describe_posterior(model, rope_range = list(c(1, 2)))
#> Error: Length of `range` (i.e. number of ROPE limits) should match the number
#>   of parameters.

# errors - wrong string
describe_posterior(model, rope_range = list(c(10, 40), c(-5, -4), "a"))
#> Error: `range` should be 'default' or a vector of 2 numeric values (e.g.,
#>   c(-0.1, 0.1)).

# errors - wrong name
describe_posterior(model, rope_range = list(wx = c(-5, -4)))
#> Error: Not all elements of `range` were found in the parameters. Please check
#>   following names: wx

^{Created on 2024-09-17 with reprex v2.1.1}