Closed IndrajeetPatil closed 3 years ago
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)
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)
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)
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?
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.
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)
.
Created on 2020-11-24 by the reprex package (v0.3.0.9001)