including prior information in model_parameters output for metaBMA objects

[IndrajeetPatil]

library(metaBMA)
#> Loading required package: Rcpp

data(towels)
set.seed(123)
mr <- meta_random(logOR, SE, study, data = towels,
                  d = prior("norm", c(mean=0, sd=.3), lower = 0),
                  tau = prior("invgamma", c(shape = 1, scale = 0.15)),
                  rel.tol=.Machine$double.eps^.15, iter=1000)

insight::get_priors(mr)
#>     Parameter Distribution Location Scale
#> 1 (Intercept)         norm        0  0.30
#> 2         tau     invgamma        1  0.15

parameters::model_parameters(mr)
#> # Studies
#> 
#> Parameter                                           | Coefficient |   SE |        95% CI | Weight
#> -------------------------------------------------------------------------------------------------
#> Goldstein, Cialdini, & Griskevicius (2008), Exp. 1  |        0.38 | 0.20 | [-0.01, 0.77] |  25.59
#> Goldstein, Cialdini, & Griskevicius  (2008), Exp. 2 |        0.30 | 0.14 | [ 0.04, 0.57] |  53.97
#> Schultz, Khazian, & Zaleski (2008), Exp. 2          |        0.21 | 0.19 | [-0.17, 0.58] |  27.24
#> Schultz, Khazian, & Zaleski (2008), Exp. 3          |        0.25 | 0.17 | [-0.08, 0.58] |  34.57
#> Mair & Bergin-Seers (2010), Exp. 1                  |        0.29 | 0.82 | [-1.33, 1.90] |   1.47
#> Bohner & Schluter (2014), Exp. 1                    |       -0.12 | 0.25 | [-0.61, 0.36] |  16.25
#> Bohner & Schluter (2014), Exp. 2                    |       -1.46 | 0.76 | [-2.95, 0.03] |   1.73
#> 
#> # Meta-Parameters
#> 
#> Parameter | Coefficient |   SE |       95% CI |   BF |  Rhat |     ESS
#> ----------------------------------------------------------------------
#>           |        0.19 | 0.09 | [0.02, 0.35] | 3.81 | 1.001 | 1077.40
#> tau       |        0.13 | 0.09 | [0.02, 0.32] |      | 1.002 |  928.20
#> Using highest density intervals as credible intervals.

^{Created on 2020-11-24 by the reprex package (v0.3.0.9001)}

[strengejacke]

library(parameters)
library(metaBMA)
#> Loading required package: Rcpp

data(towels)
set.seed(123)
mr <- meta_random(logOR, SE, study, data = towels,
                  d = prior("norm", c(mean=0, sd=.3), lower = 0),
                  tau = prior("invgamma", c(shape = 1, scale = 0.15)),
                  rel.tol=.Machine$double.eps^.15, iter=1000)

mr2 <- meta_fixed(logOR, SE, study, data = towels)

model_parameters(mr)
#> # Studies
#> 
#> Parameter                                           | Coefficient |   SE |        95% CI | Weight
#> -------------------------------------------------------------------------------------------------
#> Goldstein, Cialdini, & Griskevicius (2008), Exp. 1  |        0.38 | 0.20 | [-0.01, 0.77] |  25.59
#> Goldstein, Cialdini, & Griskevicius  (2008), Exp. 2 |        0.30 | 0.14 | [ 0.04, 0.57] |  53.97
#> Schultz, Khazian, & Zaleski (2008), Exp. 2          |        0.21 | 0.19 | [-0.17, 0.58] |  27.24
#> Schultz, Khazian, & Zaleski (2008), Exp. 3          |        0.25 | 0.17 | [-0.08, 0.58] |  34.57
#> Mair & Bergin-Seers (2010), Exp. 1                  |        0.29 | 0.82 | [-1.33, 1.90] |   1.47
#> Bohner & Schluter (2014), Exp. 1                    |       -0.12 | 0.25 | [-0.61, 0.36] |  16.25
#> Bohner & Schluter (2014), Exp. 2                    |       -1.46 | 0.76 | [-2.95, 0.03] |   1.73
#> 
#> # Meta-Parameters
#> 
#> Parameter | Coefficient |   SE |       95% CI |   BF |  Rhat |     ESS |                Prior
#> ---------------------------------------------------------------------------------------------
#> Overall   |        0.19 | 0.09 | [0.02, 0.35] | 3.81 | 1.001 | 1077.40 |     Norm (0 +- 0.30)
#> tau       |        0.13 | 0.09 | [0.02, 0.32] |      | 1.002 |  928.20 | Invgamma (1 +- 0.15)
#> Using highest density intervals as credible intervals.

model_parameters(mr2)
#> # Studies
#> 
#> Parameter                                           | Coefficient |   SE |        95% CI | Weight
#> -------------------------------------------------------------------------------------------------
#> Goldstein, Cialdini, & Griskevicius (2008), Exp. 1  |        0.38 | 0.20 | [-0.01, 0.77] |  25.59
#> Goldstein, Cialdini, & Griskevicius  (2008), Exp. 2 |        0.30 | 0.14 | [ 0.04, 0.57] |  53.97
#> Schultz, Khazian, & Zaleski (2008), Exp. 2          |        0.21 | 0.19 | [-0.17, 0.58] |  27.24
#> Schultz, Khazian, & Zaleski (2008), Exp. 3          |        0.25 | 0.17 | [-0.08, 0.58] |  34.57
#> Mair & Bergin-Seers (2010), Exp. 1                  |        0.29 | 0.82 | [-1.33, 1.90] |   1.47
#> Bohner & Schluter (2014), Exp. 1                    |       -0.12 | 0.25 | [-0.61, 0.36] |  16.25
#> Bohner & Schluter (2014), Exp. 2                    |       -1.46 | 0.76 | [-2.95, 0.03] |   1.73
#> 
#> # Meta-Parameters
#> 
#> Parameter | Coefficient |   SE |       95% CI |    BF |            Prior
#> ------------------------------------------------------------------------
#> Overall   |        0.21 | 0.08 | [0.06, 0.36] | 11.97 | Norm (0 +- 0.30)
#> Using highest density intervals as credible intervals.

^{Created on 2020-11-24 by the reprex package (v0.3.0)}

[IndrajeetPatil]

I think there is something wrong here. The prior information doesn't match with prior information in the model object.

I would have expected the prior here to be 0.707.

# setup
set.seed(123)
library(metaBMA)
#> Loading required package: Rcpp

# creating a dataframe
df1 <-
  structure(
    .Data = list(
      term = c("1", "2", "3", "4", "5"),
      estimate = c(
        0.382047603321706,
        0.780783111514665,
        0.425607573765058,
        0.558365541235078,
        0.956473848429961
      ),
      std.error = c(
        0.0465576338644502,
        0.0330218199731529,
        0.0362834986178494,
        0.0480571500648261,
        0.062215818388157
      )
    ),
    row.names = c(NA, -5L),
    class = c("tbl_df", "tbl", "data.frame")
  )

# model
set.seed(123)
mod <- 
  meta_random(
  y = df1$estimate,
    SE = df1$std.error,
  iter = 1000, 
  summarize = "integrate"
)
#> Warning: There were 1 divergent transitions after warmup. See
#> http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
#> to find out why this is a problem and how to eliminate them.
#> Warning: Examine the pairs() plot to diagnose sampling problems
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> http://mc-stan.org/misc/warnings.html#bulk-ess

mod$jzs
#> $rscale_contin
#> [1] 0.5
#> 
#> $rscale_discrete
#> [1] 0.7071068
#> 
#> $centering
#> [1] TRUE

insight::get_priors(mod)
#>     Parameter Distribution Location Scale
#> 1 (Intercept)         norm        0  0.30
#> 2         tau     invgamma        1  0.15

^{Created on 2020-11-24 by the reprex package (v0.3.0.9001)}

[strengejacke]

I think we had this discussion before. insight::get_priors() returns what you get from print() or summary(), thus being in line with the "official" output:

mod
#> ### Bayesian Random-Effects Meta-Analysis ### 
#>    Prior on d:      'norm' (mean=0, sd=0.3) with support on the interval [-Inf,Inf]. 
#>    Prior on tau:    'invgamma' (shape=1, scale=0.15) with support on the interval [0,Inf]. 
#> 
#> # Bayes factors:
#>            (denominator)
#> (numerator) random_H0 random_H1
#>   random_H0       1.0    0.0354
#>   random_H1      28.2    1.0000
#> 
#> # Posterior summary statistics of random-effects model:
#>      mean    sd  2.5%   50% 97.5% hpd95_lower hpd95_upper n_eff Rhat
#> d   0.518 0.140 0.180 0.536 0.743       0.219       0.766    NA   NA
#> tau 0.288 0.141 0.134 0.253 0.653       0.106       0.553    NA   NA

^{Created on 2020-11-24 by the reprex package (v0.3.0)}

[IndrajeetPatil]

Hmm, I see.

Okay, I will still be curious to hear from @danheck about what's the official take on this.
What should we be returning in the Prior_Scale column of the output?

[danheck]

Sorry for the confusion - the internal numbers rscale_contin and rscale_discrete refer to the scale of the JZS prior for moderator analysis. This is a feature not yet documented. The numbers are defaults ignored in the standard analyses without moderators documented above.

However, one could argue that the default prior normal(0, 0.30) is not a good choice. This was reasonable in one of the first publications in which we used this approach. But it seems that JASP (based on metaBMA) now uses the more common cauchy(0, 0.707) . It would be relatively easy to change the default i the package as well, and I think the consequences would be rather minor for users.

[IndrajeetPatil]

Thanks a lot, @danheck! This is really helpful and clears up the confusion on my side.

I would also vote for defaulting to cauchy(0, 0.707).