BayesFactor icon indicating copy to clipboard operation
BayesFactor copied to clipboard

`emmeans` integration

Open mattansb opened this issue 6 years ago • 9 comments
trafficstars

It would be great to be able to estimate contrasts / simple effects / simple slopes and their credible intervals from the posterior using emmeans (instead of doing this manually by summing up the correct combination of parameters)

This would require is adding (and exporting) two new functions: recover_data.BFBayesFactor and emm_basis.BFBayesFactor. Some more info can be found here.

Thanks!

mattansb avatar Jan 07 '19 06:01 mattansb

I'm not sure what I'd do for the dffun part of the emm_basis; it seems geared toward classical tests. Do you have any suggestions?

richarddmorey avatar Jan 08 '19 14:01 richarddmorey

Looking at the emm_basis.stanreg and emm_basis.brmsfit, dffun should be set to function(k, dfargs) Inf (e.g., here). It seems that when the output from emm_basis contains the argument post.beta, that represents sampling from the posterior distribution, emmeans changes its handling of the object from a Frequentist to a Bayesian framework.

mattansb avatar Jan 08 '19 17:01 mattansb

Since BayesFactor::posterior can be coerced into an mcmc object with BayesFactor:::as.mcmc.BFmcmc, it might be possible to pass the object to emmeans:::emm_basis.mcmc?

I've tried playing around with this, but was unsuccessful - but maybe you'll find this useful / inspiring?

mattansb avatar Jan 11 '19 09:01 mattansb

(sorry for spamming you - I seem to have found myself intensely working at this) This issue can be set as resolved (as can #133 ). I will add the code to my BFEffects package, unless you think it should be part of the main BayesFactor.

mattansb avatar Mar 11 '19 18:03 mattansb

Here is my attempt at this integration.

It's not perfect, but it's pretty close. Feel free to make any (or no) use of this going forward.

mattansb avatar May 12 '19 08:05 mattansb

Do you have some example code that uses this?

On Sun, May 12, 2019 at 9:09 AM Mattan S. Ben-Shachar < [email protected]> wrote:

Here is my attempt at this integration. https://gist.github.com/mattansb/f383af8c0dfe88922fb5c75e8572d03e

It's not perfect, but it's pretty close. Feel free to make any (or no) use of this going forward.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/richarddmorey/BayesFactor/issues/128#issuecomment-491574669, or mute the thread https://github.com/notifications/unsubscribe-auth/AAJZVWVJTH5KIGIRNQOPOOLPU7GB7ANCNFSM4GOL3ARA .

-- Richard D. Morey Reader, School of Psychology Cardiff University http://psych.cf.ac.uk/morey

richarddmorey avatar May 12 '19 16:05 richarddmorey

library(BayesFactor)
#> Loading required package: coda
#> Loading required package: Matrix
#> ************
#> Welcome to BayesFactor 0.9.12-4.2. If you have questions, please contact Richard Morey ([email protected]).
#> 
#> Type BFManual() to open the manual.
#> ************
library(emmeans)
library(magrittr)
library(purrr)
#> 
#> Attaching package: 'purrr'
#> The following object is masked from 'package:magrittr':
#> 
#>     set_names

source("https://gist.githubusercontent.com/mattansb/f383af8c0dfe88922fb5c75e8572d03e/raw/9616f4a2c959716a1e82bb4a20fd60a8fe8d0f3e/BFBayesFactor_2_emmeans.R")
options(contrasts=c('contr.sum', 'contr.poly'))

data(puzzles)
fit_freq <- aov(RT ~ shape*color + Error(ID/(shape*color)), data = puzzles)

fit_BF <- anovaBF(RT ~ shape*color + ID, data = puzzles, whichRandom = "ID")


emmip(fit_freq,color~shape, CIs = TRUE)
#> Registered S3 methods overwritten by 'ggplot2':
#>   method         from 
#>   [.quosures     rlang
#>   c.quosures     rlang
#>   print.quosures rlang

emmip(fit_BF[4],color~shape, CIs = TRUE)



emmeans(fit_freq,pairwise~shape|color)
#> $emmeans
#> color = color:
#>  shape  emmean    SE   df lower.CL upper.CL
#>  round      45 0.733 16.9     43.5     46.5
#>  square     44 0.733 16.9     42.5     45.5
#> 
#> color = monochromatic:
#>  shape  emmean    SE   df lower.CL upper.CL
#>  round      46 0.733 16.9     44.5     47.5
#>  square     45 0.733 16.9     43.5     46.5
#> 
#> Warning: EMMs are biased unless design is perfectly balanced 
#> Confidence level used: 0.95 
#> 
#> $contrasts
#> color = color:
#>  contrast       estimate    SE   df t.ratio p.value
#>  round - square        1 0.603 20.5 1.658   0.1125 
#> 
#> color = monochromatic:
#>  contrast       estimate    SE   df t.ratio p.value
#>  round - square        1 0.603 20.5 1.658   0.1125
emmeans(fit_BF[4],pairwise~shape|color)
#> $emmeans
#> color = color:
#>  shape  emmean lower.HPD upper.HPD
#>  round    45.0      43.6      46.5
#>  square   44.1      42.7      45.7
#> 
#> color = monochromatic:
#>  shape  emmean lower.HPD upper.HPD
#>  round    45.9      44.4      47.4
#>  square   45.0      43.6      46.6
#> 
#> Results are given on the : (not the response) scale. 
#> HPD interval probability: 0.95 
#> 
#> $contrasts
#> color = color:
#>  contrast       estimate lower.HPD upper.HPD
#>  round - square    0.866    -0.133      1.80
#> 
#> color = monochromatic:
#>  contrast       estimate lower.HPD upper.HPD
#>  round - square    0.847    -0.160      1.84
#> 
#> HPD interval probability: 0.95

Created on 2019-05-12 by the reprex package (v0.2.1)

As you can see, it basically works, but things would sometimes get tricky with categorical-by-continuous interactions, that I just couldn't figure out.

mattansb avatar May 12 '19 17:05 mattansb

Hi, nice to see the link to emmeans. Is the issue with categorical-by-continuous interactions still a thing? Cheers Jana

JanaJarecki avatar Dec 02 '19 15:12 JanaJarecki

No change on this end. It seems that the implementation in JASP has some method of posterior estimates, but I'm not sure.

mattansb avatar Dec 02 '19 18:12 mattansb