search

easystats / bayestestR

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

no more plotting for equivalence_test on difference between posterior distributions of model parameters #263

Closed aschetti closed 4 years ago

aschetti commented 4 years ago

Describe the bug Hi, first off, thanks for the package: I have been using it for a couple of projects and it is very user-friendly! I don't know if what I am about to describe can be considered a bug, but I recently found out that the plot function does not work anymore for equivalence_test. Below you can find a reproducible example. After running a model in brms, I compute the difference between posterior distributions of the parameters of interest and run equivalence_test on the resulting posterior. I try to plot the results but get the following error:

Error in find_parameters.data.frame(x, effects = "all", component = "all", : A data frame is no valid object for this function.

A note in the help says: "[...] there is a plot()-method to visualize the results from the equivalence-test (for models only).". However, the code below was working until package version 0.3 (if I remember correctly).

To Reproduce

Here is the reprex output (using simulated data):


library(tidyverse)
library(brms)
library(bayestestR)
library(see)

seed_reprex <- 1
set.seed(seed_reprex)

# create fake data
data <- tibble(
  participant = as_factor(1:100),
  condition = as_factor(rep(1:2, each = 50)),
  variable = c(
    rnorm(50, mean = 0, sd = 1),
    rnorm(50, mean = 1, sd = 1)
  )
)

# fit model
model <- brm(variable ~ 0 + condition,
  data = data,
  seed = seed_reprex
)

# difference between posterior distributions of conditions
model_posterior_samples <- posterior_samples(model) %>%
  as_tibble() %>%
  mutate(
    "cond1-cond2" = b_condition1 - b_condition2
  ) %>%
  dplyr::select("cond1-cond2")

# equivalence test
model_ropeHDI <- equivalence_test(model_posterior_samples,
  range = c(-0.05, 0.05), # ROPE
  ci = 1
)

# decision on H0
plot(model_ropeHDI)
#> Error in find_parameters.data.frame(x, effects = "all", component = "all", : A data frame is no valid object for this function.

Created on 2019-12-06 by the reprex package (v0.3.0)

Session info ``` r devtools::session_info() #> ─ Session info ─────────────────────────────────────────────────────────────── #> setting value #> version R version 3.6.1 (2019-07-05) #> os Ubuntu 18.04.3 LTS #> system x86_64, linux-gnu #> ui X11 #> language (EN) #> collate en_US.UTF-8 #> ctype en_US.UTF-8 #> tz Europe/Amsterdam #> date 2019-12-06 #> #> ─ Packages ─────────────────────────────────────────────────────────────────── #> package * version date lib source #> abind 1.4-5 2016-07-21 [1] CRAN (R 3.6.0) #> assertthat 0.2.1 2019-03-21 [1] CRAN (R 3.6.0) #> backports 1.1.5 2019-10-02 [1] CRAN (R 3.6.1) #> base64enc 0.1-3 2015-07-28 [1] CRAN (R 3.6.0) #> bayesplot 1.7.1 2019-12-01 [1] CRAN (R 3.6.1) #> bayestestR * 0.4.0 2019-10-20 [1] CRAN (R 3.6.1) #> bridgesampling 0.7-2 2019-07-21 [1] CRAN (R 3.6.1) #> brms * 2.10.0 2019-08-29 [1] CRAN (R 3.6.1) #> Brobdingnag 1.2-6 2018-08-13 [1] CRAN (R 3.6.0) #> broom 0.5.2 2019-04-07 [1] CRAN (R 3.6.0) #> callr 3.3.2 2019-09-22 [1] CRAN (R 3.6.1) #> cellranger 1.1.0 2016-07-27 [1] CRAN (R 3.6.0) #> cli 1.1.0 2019-03-19 [1] CRAN (R 3.6.0) #> coda 0.19-3 2019-07-05 [1] CRAN (R 3.6.0) #> codetools 0.2-16 2018-12-24 [1] CRAN (R 3.6.0) #> colorspace 1.4-1 2019-03-18 [1] CRAN (R 3.6.0) #> colourpicker 1.0 2017-09-27 [1] CRAN (R 3.6.0) #> crayon 1.3.4 2017-09-16 [1] CRAN (R 3.6.0) #> crosstalk 1.0.0 2016-12-21 [1] CRAN (R 3.6.0) #> DBI 1.0.0 2018-05-02 [1] CRAN (R 3.6.0) #> dbplyr 1.4.2 2019-06-17 [1] CRAN (R 3.6.0) #> desc 1.2.0 2018-05-01 [1] CRAN (R 3.6.0) #> devtools 2.2.1 2019-09-24 [1] CRAN (R 3.6.1) #> digest 0.6.23 2019-11-23 [1] CRAN (R 3.6.1) #> dplyr * 0.8.3 2019-07-04 [1] CRAN (R 3.6.0) #> DT 0.10 2019-11-12 [1] CRAN (R 3.6.1) #> dygraphs 1.1.1.6 2018-07-11 [1] CRAN (R 3.6.0) #> effectsize 0.0.1 2019-11-15 [1] CRAN (R 3.6.1) #> ellipsis 0.3.0 2019-09-20 [1] CRAN (R 3.6.1) #> evaluate 0.14 2019-05-28 [1] CRAN (R 3.6.0) #> fastmap 1.0.1 2019-10-08 [1] CRAN (R 3.6.1) #> forcats * 0.4.0 2019-02-17 [1] CRAN (R 3.6.0) #> fs 1.3.1 2019-05-06 [1] CRAN (R 3.6.0) #> generics 0.0.2 2018-11-29 [1] CRAN (R 3.6.0) #> ggplot2 * 3.2.1 2019-08-10 [1] CRAN (R 3.6.1) #> ggridges 0.5.1 2018-09-27 [1] CRAN (R 3.6.0) #> glue 1.3.1.9000 2019-11-17 [1] Github (tidyverse/glue@c02d7d4) #> gridExtra 2.3 2017-09-09 [1] CRAN (R 3.6.0) #> gtable 0.3.0 2019-03-25 [1] CRAN (R 3.6.0) #> gtools 3.8.1 2018-06-26 [1] CRAN (R 3.6.0) #> haven 2.2.0 2019-11-08 [1] CRAN (R 3.6.1) #> highr 0.8 2019-03-20 [1] CRAN (R 3.6.0) #> hms 0.5.2 2019-10-30 [1] CRAN (R 3.6.1) #> htmltools 0.4.0 2019-10-04 [1] CRAN (R 3.6.1) #> htmlwidgets 1.5.1 2019-10-08 [1] CRAN (R 3.6.1) #> httpuv 1.5.2 2019-09-11 [1] CRAN (R 3.6.1) #> httr 1.4.1 2019-08-05 [1] CRAN (R 3.6.1) #> igraph 1.2.4.2 2019-11-27 [1] CRAN (R 3.6.1) #> inline 0.3.15 2018-05-18 [1] CRAN (R 3.6.0) #> insight 0.7.1 2019-11-28 [1] CRAN (R 3.6.1) #> jsonlite 1.6 2018-12-07 [1] CRAN (R 3.6.0) #> knitr 1.26 2019-11-12 [1] CRAN (R 3.6.1) #> later 1.0.0 2019-10-04 [1] CRAN (R 3.6.1) #> lattice 0.20-38 2018-11-04 [1] CRAN (R 3.6.0) #> lazyeval 0.2.2 2019-03-15 [1] CRAN (R 3.6.0) #> lifecycle 0.1.0 2019-08-01 [1] CRAN (R 3.6.1) #> loo 2.1.0 2019-03-13 [1] CRAN (R 3.6.0) #> lubridate 1.7.4 2018-04-11 [1] CRAN (R 3.6.0) #> magrittr 1.5 2014-11-22 [1] CRAN (R 3.6.0) #> markdown 1.1 2019-08-07 [1] CRAN (R 3.6.1) #> Matrix 1.2-18 2019-11-27 [1] CRAN (R 3.6.1) #> matrixStats 0.55.0 2019-09-07 [1] CRAN (R 3.6.1) #> memoise 1.1.0 2017-04-21 [1] CRAN (R 3.6.0) #> mime 0.7 2019-06-11 [1] CRAN (R 3.6.0) #> miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 3.6.0) #> modelr 0.1.5 2019-08-08 [1] CRAN (R 3.6.1) #> munsell 0.5.0 2018-06-12 [1] CRAN (R 3.6.0) #> nlme 3.1-142 2019-11-07 [1] CRAN (R 3.6.1) #> parameters 0.3.0 2019-11-20 [1] CRAN (R 3.6.1) #> pillar 1.4.2 2019-06-29 [1] CRAN (R 3.6.0) #> pkgbuild 1.0.6 2019-10-09 [1] CRAN (R 3.6.1) #> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 3.6.1) #> pkgload 1.0.2 2018-10-29 [1] CRAN (R 3.6.0) #> plyr 1.8.4 2016-06-08 [1] CRAN (R 3.6.0) #> prettyunits 1.0.2 2015-07-13 [1] CRAN (R 3.6.0) #> processx 3.4.1 2019-07-18 [1] CRAN (R 3.6.1) #> promises 1.1.0 2019-10-04 [1] CRAN (R 3.6.1) #> ps 1.3.0 2018-12-21 [1] CRAN (R 3.6.0) #> purrr * 0.3.3 2019-10-18 [1] CRAN (R 3.6.1) #> R6 2.4.1 2019-11-12 [1] CRAN (R 3.6.1) #> Rcpp * 1.0.3 2019-11-08 [1] CRAN (R 3.6.1) #> readr * 1.3.1 2018-12-21 [1] CRAN (R 3.6.0) #> readxl 1.3.1 2019-03-13 [1] CRAN (R 3.6.0) #> remotes 2.1.0 2019-06-24 [1] CRAN (R 3.6.0) #> reprex 0.3.0 2019-05-16 [1] CRAN (R 3.6.0) #> reshape2 1.4.3 2017-12-11 [1] CRAN (R 3.6.0) #> rlang 0.4.2 2019-11-23 [1] CRAN (R 3.6.1) #> rmarkdown 1.18 2019-11-27 [1] CRAN (R 3.6.1) #> rprojroot 1.3-2 2018-01-03 [1] CRAN (R 3.6.0) #> rsconnect 0.8.15 2019-07-22 [1] CRAN (R 3.6.1) #> rstan 2.19.2 2019-07-09 [1] CRAN (R 3.6.1) #> rstantools 2.0.0 2019-09-15 [1] CRAN (R 3.6.1) #> rvest 0.3.5 2019-11-08 [1] CRAN (R 3.6.1) #> scales 1.1.0 2019-11-18 [1] CRAN (R 3.6.1) #> see * 0.3.0 2019-11-19 [1] CRAN (R 3.6.1) #> sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.6.0) #> shiny 1.4.0 2019-10-10 [1] CRAN (R 3.6.1) #> shinyjs 1.0 2018-01-08 [1] CRAN (R 3.6.0) #> shinystan 2.5.0 2018-05-01 [1] CRAN (R 3.6.0) #> shinythemes 1.1.2 2018-11-06 [1] CRAN (R 3.6.0) #> StanHeaders 2.19.0 2019-09-07 [1] CRAN (R 3.6.1) #> stringi 1.4.3 2019-03-12 [1] CRAN (R 3.6.0) #> stringr * 1.4.0 2019-02-10 [1] CRAN (R 3.6.0) #> testthat 2.3.1 2019-12-01 [1] CRAN (R 3.6.1) #> threejs 0.3.1 2017-08-13 [1] CRAN (R 3.6.0) #> tibble * 2.1.3 2019-06-06 [1] CRAN (R 3.6.0) #> tidyr * 1.0.0 2019-09-11 [1] CRAN (R 3.6.1) #> tidyselect 0.2.5 2018-10-11 [1] CRAN (R 3.6.0) #> tidyverse * 1.3.0 2019-11-21 [1] CRAN (R 3.6.1) #> usethis 1.5.1 2019-07-04 [1] CRAN (R 3.6.0) #> vctrs 0.2.0 2019-07-05 [1] CRAN (R 3.6.0) #> withr 2.1.2 2018-03-15 [1] CRAN (R 3.6.0) #> xfun 0.11 2019-11-12 [1] CRAN (R 3.6.1) #> xml2 1.2.2 2019-08-09 [1] CRAN (R 3.6.1) #> xtable 1.8-4 2019-04-21 [1] CRAN (R 3.6.0) #> xts 0.11-2 2018-11-05 [1] CRAN (R 3.6.0) #> yaml 2.2.0 2018-07-25 [1] CRAN (R 3.6.0) #> zeallot 0.1.0 2018-01-28 [1] CRAN (R 3.6.0) #> zoo 1.8-6 2019-05-28 [1] CRAN (R 3.6.0) #> #> [1] /home/aschetti/R/x86_64-pc-linux-gnu-library/3.6 #> [2] /usr/local/lib/R/site-library #> [3] /usr/lib/R/site-library #> [4] /usr/lib/R/library ```

Expected behaviour A plot similar to what is displayed in your vignette, not only for posterior distributions of model parameters but also for posterior distributions of their difference.

Specifications (please complete the following information): bayestestR 0.4.0. For additional information, see sessionInfo() in reprex output.

strengejacke commented 4 years ago

I'm not sure if we really had a plot()-method for equivalence_test() with data frames?

It works with rope(), but I think we never had a plot-method for equivalence_test with data frames as input.

plot(rope(model_posterior_samples, range = c(-0.05, 0.05), ci = 1))
aschetti commented 4 years ago

Hi @strengejacke, thanks for the quick reply! As I said, I seem to recall that there was a plot()-method for equivalence_test() with data frames, but I trust your memory more than mine :grin: The code you posted works perfectly, I'm just going to use that. Again, thanks to you and your collaborators for such a great set of packages!

strengejacke commented 4 years ago

but I trust your memory more than mine

I wouldn't bet on my memory ;-) but good to know if the alternative code with rope() works for you.

aschetti commented 4 years ago

Hi @strengejacke, my apologies, but I think I was overly optimistic when I closed this issue. I rechecked my old code and found out that I could produce the plot below with

fit_lm_ropeHDI <- equivalence_test(fit_lm_posterior_samples,
                 range = c(-0.05, 0.05), # ROPE
                  ci = 1)

plot(fit_lm_ropeHDI) +
  theme_custom

plot_equivalence_test

That is what I would like to reproduce, because it also visually helps on the decision on H0. So, it seems that plot() could work with data frames. It would actually be very helpful also when using the function p_direction, whose output can be plotted if the input is a (brms) model but not a data frame with samples from the posteriors.

pd <- p_direction(fit_lm_posterior_samples)
plot(pd)

#> Error: Failed at retrieving data :( Please provide original model or data through the `data` argument

Thanks and sorry for reopening the issue!