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

Additional bayesfactor_savagedickey methods #103

Closed mattansb closed 5 years ago

mattansb commented 5 years ago

Add

mattansb commented 5 years ago

Currently working on .brmsfit.. It might only work for fixed effects.

mattansb commented 5 years ago

@DominiqueMakowski @strengejacke I've added the 2 methods (brmsfit and data.frame). Seem to have an issue I'm having trouble with in .brmsfit where I can only get priors for fixed effects parameters. If neither of you can fix this, should the .stanreg method also only return fixed effects? (I'm inclined to say yes).

DominiqueMakowski commented 5 years ago

If I am not mistaken currently the other functions (for rstanarm) return by default only the fixed effects, but this behaviour is adjustable via a series of arguments

mattansb commented 5 years ago

And for brmsfit?

DominiqueMakowski commented 5 years ago

@strengejacke is the expert of brmsfit 😬

strengejacke commented 5 years ago

I try to look at it the next days...

DominiqueMakowski commented 5 years ago

Would it make sense to add a method for BFBayesFactor objects? Maybe it could simply extract the BF of the BF of the object?

(I am thinking about this because currently describe_posterior(BFBayesFactor) fails if test="all" (because it tries to run a method for BFBayesFactor))

mattansb commented 5 years ago

Currently bayesfactor_models does just that - it simply extracts the BF for each model in the object:

library(BayesFactor)
library(bayestestR)

data(puzzles)
result <- anovaBF(RT ~ shape*color + ID, data = puzzles, whichRandom = "ID",
                  whichModels = 'top', progress=FALSE)
result
#> Bayes factor top-down analysis
#> --------------
#> When effect is omitted from shape + color + shape:color + ID , BF is...
#> [1] Omit color:shape : 2.878847  ±2.54%
#> [2] Omit color       : 0.2457026 ±2.66%
#> [3] Omit shape       : 0.2462095 ±2.36%
#> 
#> Against denominator:
#>   RT ~ shape + color + shape:color + ID 
#> ---
#> Bayes factor type: BFlinearModel, JZS

bayesfactor_models(result)
#> Bayes factor analysis
#> --------------                                    
#> [2] shape + color + ID          2.88
#> [3] shape + shape:color + ID    0.25
#> [4] color + shape:color + ID    0.25
#> 
#> Against denominator:
#>   [1] shape + color + shape:color + ID   
#> ---
#> Bayes factor type:  JZS (BayesFactor)

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

Generally, there's a problem computing sdBF for BFBayesFactor, as there is no way to extract / sample from priors for the parameters...

mattansb commented 5 years ago

So it seems that brmsfit cannot sample from random-effects' priors, which means that sdBF cannot be computed for random effects (in brms::hypothesis the sdBF is called Evid.Ratio):

library(brms)
library(bayestestR)

get_prior(extra ~ group + (1|ID), sleep)

priors <- 
  set_prior("student_t(3, 0, 1)",class = "b") + 
  set_prior("student_t(3, 0, 1)",class = "sd", group = "ID")
priors

fit <- brm(extra ~ group + (1|ID), sleep, sample_prior = TRUE, prior = priors)

hypothesis(fit, "Intercept=0", group = "ID", scope = "ranef")
#> Hypothesis Tests for class :
#>    Group      Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
#> 1      1 (Intercept) = 0    -0.19      0.80    -1.79     1.41         NA        NA     
#> 2      2 (Intercept) = 0    -1.45      0.91    -3.19     0.31         NA        NA     
#> 3      3 (Intercept) = 0    -0.82      0.84    -2.51     0.83         NA        NA     
#> 4      4 (Intercept) = 0    -1.54      0.92    -3.40     0.23         NA        NA     
#> 5      5 (Intercept) = 0    -1.23      0.89    -3.02     0.47         NA        NA     
#> 6      6 (Intercept) = 0     1.76      0.96    -0.10     3.64         NA        NA     
#> 7      7 (Intercept) = 0     2.28      1.04     0.02     4.24         NA        NA    *
#> 8      8 (Intercept) = 0    -0.25      0.81    -1.87     1.38         NA        NA     
#> 9      9 (Intercept) = 0     0.56      0.86    -1.13     2.34         NA        NA     
#> 10    10 (Intercept) = 0     0.86      0.84    -0.77     2.54         NA        NA     
#> ---
#> '*': The expected value under the hypothesis lies outside the 95%-CI.
#> Posterior probabilities of point hypotheses assume equal prior probabilities.

Thus my solution is to give a warning, and return NA:

bayesfactor_savagedickey(fit, effects = "random")
#> # Bayes Factor (Savage-Dickey density ratio)
#> 
#>             Parameter Bayes Factor
#> 1   r_ID.1.Intercept.           NA
#> 2   r_ID.2.Intercept.           NA
#> 3   r_ID.3.Intercept.           NA
#> 4   r_ID.4.Intercept.           NA
#> 5   r_ID.5.Intercept.           NA
#> 6   r_ID.6.Intercept.           NA
#> 7   r_ID.7.Intercept.           NA
#> 8   r_ID.8.Intercept.           NA
#> 9   r_ID.9.Intercept.           NA
#> 10 r_ID.10.Intercept.           NA
#> ---
#> Test Value: 0
#> Warning message:
#> In bayesfactor_savagedickey.brmsfit(fit, effects = "random") :
#>   Cannot compute Savage-Dickey BFs for random effect from 'brmsfit' objects.

Other than that, for the fixed effect, bayesfactor_savagedickey.brmsfit returns equivalent results to brms::hypothesis (not that Evid.Ratio = posterior/prior whereas our sdBF = prior/posterior):

hypothesis(fit, c("b_Intercept=0","b_group2=0"), class = NULL)
#> Hypothesis Tests for class :
#>          Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
#> 1 (b_Intercept) = 0     0.89      0.64    -0.38     2.16       0.59      0.37     
#> 2    (b_group2) = 0     1.31      0.51     0.23     2.23       0.12      0.11    *
#> ---
#> '*': The expected value under the hypothesis lies outside the 95%-CI.
#> Posterior probabilities of point hypotheses assume equal prior probabilities.

bayesfactor_savagedickey(fit, effects = "fixed")
#> # Bayes Factor (Savage-Dickey density ratio)
#> 
#>     Parameter Bayes Factor
#> 1 b_Intercept          1.7
#> 2    b_group2          8.7
#> ---
#> Test Value: 0

1/0.59
#> [1] 1.694915

1/0.12
#> [1] 8.333333
DominiqueMakowski commented 5 years ago

Currently bayesfactor_models does just that - it simply extracts the BF for each model in the object:

Does it make sense for you, to run (in describe_posterior) bayesfactor_models() for BFBayesFactor objects, when BFsds are run for all other types of objects? What would you proceed?

mattansb commented 5 years ago

I suggest not supporting sdBF for BFBayesFactor in describe_posterior at this time.

Maybe in the future we can create our own function to extract / sample priors for BFBayesFactor, at which point we would be able to add a sdBF method for BFBayesFactor.

mattansb commented 5 years ago

Also, it seems that describe_posterior currently doesn't really support BFBayesFactor with linear models, only ttests / correlations.

How about I add a non exported sdBF function for these cases as a temporary* solution?

(* In Hebrew we have a saying: What is temporary, is permanent...)

strengejacke commented 5 years ago

So it seems that brmsfit cannot sample from random-effects' priors,

I think this is because the prior for random effects is set to sd (set_prior("student_t(3, 0, 1)",class = "sd", group = "ID")), that's why you can only test the variance of the random effects, not their "intercepts".

library(brms)
priors <- 
  set_prior("student_t(3, 0, 1)",class = "b") + 
  set_prior("student_t(3, 0, 1)",class = "sd", group = "ID")

fit <- brm(extra ~ group + (1|ID), sleep, sample_prior = TRUE, prior = priors)

hypothesis(fit, "sd_ID__Intercept = 0", class = NULL)
#> Hypothesis Tests for class :
#>               Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
#> 1 (sd_ID__Intercept) = 0     1.49      0.55     0.26     2.66       0.15
#>   Post.Prob Star
#> 1      0.13    *
#> ---
#> '*': The expected value under the hypothesis lies outside the 95%-CI.
#> Posterior probabilities of point hypotheses assume equal prior probabilities.

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

That said, I think it should be possible, because there's an example in ?hypothesis, which, however, does not work (or at least, it also returns only NA):

library(brms)
#> Loading required package: Rcpp
prior <- c(set_prior("normal(0,2)", class = "b"),
           set_prior("student_t(10,0,1)", class = "sigma"),
           set_prior("student_t(10,0,1)", class = "sd"))

## fit a linear mixed effects models
fit <- brm(time ~ age + sex + disease + (1 + age|patient),
           data = kidney, family = lognormal(),
           prior = prior, sample_prior = TRUE, 
           control = list(adapt_delta = 0.95), cores = 4)

hypothesis(fit, "age = 0", scope = "coef", group = "patient")
#> Hypothesis Tests for class :
#>    Group Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
#> 1      1  (age) = 0    -0.01      0.02    -0.04     0.03         NA
#> 2      2  (age) = 0    -0.01      0.02    -0.04     0.02         NA
#> 3      3  (age) = 0    -0.01      0.02    -0.04     0.03         NA
#> 4      4  (age) = 0     0.00      0.02    -0.04     0.03         NA
#> 5      5  (age) = 0    -0.01      0.02    -0.04     0.03         NA
#> 6      6  (age) = 0    -0.01      0.02    -0.04     0.03         NA
#> 7      7  (age) = 0    -0.01      0.02    -0.04     0.03         NA
#> 8      8  (age) = 0     0.00      0.02    -0.03     0.03         NA
#> 9      9  (age) = 0     0.00      0.02    -0.03     0.03         NA
#> 10    10  (age) = 0     0.00      0.02    -0.04     0.03         NA
#> 11    11  (age) = 0    -0.01      0.02    -0.04     0.03         NA
#> 12    12  (age) = 0    -0.01      0.02    -0.04     0.03         NA
#> 13    13  (age) = 0    -0.01      0.02    -0.04     0.03         NA
#> 14    14  (age) = 0     0.00      0.02    -0.04     0.03         NA
#> 15    15  (age) = 0    -0.01      0.02    -0.04     0.03         NA
#> 16    16  (age) = 0    -0.01      0.02    -0.04     0.02         NA
#> 17    17  (age) = 0     0.00      0.02    -0.03     0.03         NA
#> 18    18  (age) = 0     0.00      0.02    -0.03     0.03         NA
#> 19    19  (age) = 0    -0.01      0.02    -0.04     0.03         NA
#> 20    20  (age) = 0    -0.01      0.02    -0.04     0.02         NA
#> 21    21  (age) = 0     0.00      0.02    -0.03     0.04         NA
#> 22    22  (age) = 0    -0.01      0.02    -0.04     0.03         NA
#> 23    23  (age) = 0    -0.01      0.02    -0.04     0.02         NA
#> 24    24  (age) = 0    -0.01      0.02    -0.04     0.03         NA
#> 25    25  (age) = 0    -0.01      0.02    -0.04     0.03         NA
#> 26    26  (age) = 0     0.00      0.02    -0.03     0.03         NA
#> 27    27  (age) = 0    -0.01      0.02    -0.04     0.03         NA
#> 28    28  (age) = 0    -0.01      0.02    -0.04     0.03         NA
#> 29    29  (age) = 0    -0.01      0.02    -0.04     0.02         NA
#> 30    30  (age) = 0    -0.01      0.02    -0.04     0.03         NA
#> 31    31  (age) = 0    -0.01      0.02    -0.04     0.03         NA
#> 32    32  (age) = 0    -0.01      0.02    -0.05     0.02         NA
#> 33    33  (age) = 0    -0.01      0.02    -0.04     0.02         NA
#> 34    34  (age) = 0    -0.01      0.02    -0.04     0.03         NA
#> 35    35  (age) = 0    -0.01      0.02    -0.04     0.03         NA
#> 36    36  (age) = 0    -0.01      0.02    -0.04     0.03         NA
#> 37    37  (age) = 0    -0.01      0.02    -0.05     0.02         NA
#> 38    38  (age) = 0    -0.01      0.02    -0.04     0.02         NA
#> ---
#> '*': The expected value under the hypothesis lies outside the 95%-CI.
#> Posterior probabilities of point hypotheses assume equal prior probabilities.

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

Tagging @paul-buerkner - Paul, is it expected that the example shown in the docs of brms::hypothesis() (see above) only returns NA as evidence ratio?

paul-buerkner commented 5 years ago

Yes, it is expected since brms does not allow to sample from hierarchical priors if sample_prior = "yes".

strengejacke commented 5 years ago

Ok, thanks for clarification!

mattansb commented 5 years ago

Thanks @paul-buerkner !

I think I solved this (using a update(..., sample_prior = "only") to get all priors!):

library(brms)
library(bayestestR)

priors <- 
  set_prior("student_t(3, 0, 1)",class = "b") + 
  set_prior("student_t(3, 0, 1)",class = "sd", group = "ID")

fit <- brm(extra ~ group + (1|ID), sleep, prior = priors)

bayesfactor_savagedickey(fit, effects = "all")
#> Start sampling
#> # Bayes Factor (Savage-Dickey density ratio)
#> 
#>             Parameter Bayes Factor
#> 1         b_Intercept         0.16
#> 2            b_group2         8.66
#> 3   r_ID.1.Intercept.         2.25
#> 4   r_ID.2.Intercept.         8.52
#> 5   r_ID.3.Intercept.         3.74
#> 6   r_ID.4.Intercept.         9.98
#> 7   r_ID.5.Intercept.         5.34
#> 8   r_ID.6.Intercept.         8.63
#> 9   r_ID.7.Intercept.        16.54
#> 10  r_ID.8.Intercept.         2.37
#> 11  r_ID.9.Intercept.         3.11
#> 12 r_ID.10.Intercept.         4.88
#> ---
#> Evidence Against Test Value: 0
DominiqueMakowski commented 5 years ago

Is this issue solved?

mattansb commented 5 years ago

Yes