search

danheck / metaBMA

Bayesian Model Averaging for Random and Fixed Effects Meta-Analysis
https://danheck.github.io/metaBMA/
28 stars 2 forks source link

changing `ci` value doesn't change posterior estimates #7

Closed IndrajeetPatil closed 3 years ago

IndrajeetPatil commented 3 years ago

The documentation mentions-

image

But changing the value of ci, doesn't seem to change estimates in the output. Is this by design?

library(metaBMA)
#> Loading required package: Rcpp
data(towels)

# default
set.seed(123)
mr <-
  meta_random(
    logOR,
    SE,
    study,
    data = towels,
    ci = 0.95
  )
#> Warning: There were 4 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

mr$estimates
#>          mean         sd        2.5%       50%     97.5% hpd95_lower
#> d   0.1820016 0.10363018 -0.03902097 0.1877645 0.3670940 -0.02095039
#> tau 0.1355887 0.09893691  0.03348296 0.1081037 0.4037003  0.02027280
#>     hpd95_upper  n_eff Rhat
#> d     0.3811258 3974.8    1
#> tau   0.3296088 3145.4    1

# change ci value
set.seed(123)
mr <-
  meta_random(
  logOR,
  SE,
  study,
  data = towels,
  ci = 0.99
)
#> Warning: There were 4 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

mr$estimates
#>          mean         sd        2.5%       50%     97.5% hpd95_lower
#> d   0.1820016 0.10363018 -0.03902097 0.1877645 0.3670940 -0.02095039
#> tau 0.1355887 0.09893691  0.03348296 0.1081037 0.4037003  0.02027280
#>     hpd95_upper  n_eff Rhat
#> d     0.3811258 3974.8    1
#> tau   0.3296088 3145.4    1
danheck commented 3 years ago

Thanks for noticing this issue! This is clearly a bug, the argument ci is not passed to the internal summary functions.

https://github.com/danheck/metaBMA/commit/dfa24c6c274752ee662b524b4323f981ec95960f provides a solution, but I have not checked this thoroughly yet. I will update the CRAN version as soon as possible.

IndrajeetPatil commented 3 years ago

Would the column names for the intervals still contain 95? Or will it change according to the specified confidence level?

I know this would be a breaking change, but i feel like the names for these columns should be more generic, like conf.low, ci.lower, etc.

What do you think?

danheck commented 3 years ago

Yes, the names will change depending on the value of ci.

On the one hand, generic names would be beneficial for accessing the estimates (e.g., in your package or within JASP). On the other hand, having the names in the columns will tell the user what level of ci was used. And changing this will require more time and checks.

IndrajeetPatil commented 3 years ago

Another option to have a generic column names and then a separate column called ci.width or conf.level that contains this information, if the user is curious.

IndrajeetPatil commented 3 years ago

By the way, do you plan to submit the patched version to CRAN any time soon?

danheck commented 3 years ago

Good idea, thanks for the suggestion.

Since the semester is starting again, I will require 1-2 weeks for this. Is this fast enough?

IndrajeetPatil commented 3 years ago

Unfortunately, I have a deadline of 5th of November to submit revised version of my package because the tests are currently failing.

That said, if you can already tell us (cc @strengejacke) what naming schema you will adopt, we can put in contingencies in our packages to make it work with both the CRAN and GitHub versions.

Sounds reasonable?

danheck commented 3 years ago

Okay, I will try to fix this asap and check whether renaming is a bigger issue or not.

EDIT: I will let you know until tomorrow about the naming scheme.

IndrajeetPatil commented 3 years ago

Thanks, Daniel! :)

strengejacke commented 3 years ago

I revised model_parameters(), so it should now work with flexible column names that include the ci-level (e.g. hdp95_upper or hdp89_upper...).

danheck commented 3 years ago

Great, this was fast. Just for your information, the column names when using ci=.99 are currently: mean sd 0.5% 50% 99.5% hpd99_lower hpd99_upper n_eff Rhat

Since there are both quantiles and the HPD interval, a consistent renaming of the table column names would require not only renaming but also reordering the columns, for instance, as follows: mean sd median ci bci_lower bci_upper hpd_lower hpd_upper n_eff Rhat

Alternatively, I could merely remove the information about the ci in the HPD labels. The ci level would then be implicit in the reported quantiles. This would be a compromise: mean sd 0.5% 50% 99.5% hpd_lower hpd_upper n_eff Rhat

IndrajeetPatil commented 3 years ago

Thanks @strengejacke!

I like this one-

mean sd median ci bci_lower bci_upper hpd_lower hpd_upper n_eff Rhat
strengejacke commented 3 years ago

Great, this was fast. Just for your information, the column names when using ci=.99 are currently: mean sd 0.5% 50% 99.5% hpd99_lower hpd99_upper n_eff Rhat

Since there are both quantiles and the HPD interval, a consistent renaming of the table column names would require not only renaming but also reordering the columns, for instance, as follows: mean sd median ci bci_lower bci_upper hpd_lower hpd_upper n_eff Rhat

Alternatively, I could merely remove the information about the ci in the HPD labels. The ci level would then be implicit in the reported quantiles. This would be a compromise: mean sd 0.5% 50% 99.5% hpd_lower hpd_upper n_eff Rhat

Yes, the pattern to detect the appropriate CI columns in parameters::model_parameters() should be OK:

ci_method <- "hdi"
ci <- .95

switch(
  toupper(ci_method),
  "HDI" = sprintf(c("hpd%i_lower", "hpd%i_upper"), 100 * ci),
  c(sprintf("%g%%", (100 * (1 - ci)) / 2), sprintf("%g%%", 100 - (100 * (1 - ci)) / 2))
)
#> [1] "hpd95_lower" "hpd95_upper"

ci_method <- "hdi"
ci <- .99

switch(
  toupper(ci_method),
  "HDI" = sprintf(c("hpd%i_lower", "hpd%i_upper"), 100 * ci),
  c(sprintf("%g%%", (100 * (1 - ci)) / 2), sprintf("%g%%", 100 - (100 * (1 - ci)) / 2))
)
#> [1] "hpd99_lower" "hpd99_upper"

ci_method <- "eti"
ci <- .95

switch(
  toupper(ci_method),
  "HDI" = sprintf(c("hpd%i_lower", "hpd%i_upper"), 100 * ci),
  c(sprintf("%g%%", (100 * (1 - ci)) / 2), sprintf("%g%%", 100 - (100 * (1 - ci)) / 2))
)
#> [1] "2.5%"  "97.5%"

ci_method <- "eti"
ci <- .99

switch(
  toupper(ci_method),
  "HDI" = sprintf(c("hpd%i_lower", "hpd%i_upper"), 100 * ci),
  c(sprintf("%g%%", (100 * (1 - ci)) / 2), sprintf("%g%%", 100 - (100 * (1 - ci)) / 2))
)
#> [1] "0.5%"  "99.5%"

Created on 2020-10-26 by the reprex package (v0.3.0)

strengejacke commented 3 years ago

I'm ok with both approaches. It would be good to know your decision the next days, so I can already take care of this in my package update, which should be submitted the next days.

If your current implementation is mean sd 0.5% 50% 99.5% hpd99_lower hpd99_upper n_eff Rhat, it will already work.

danheck commented 3 years ago

Okay, I think I will just leave the names as they are. I am a bit hesitant due to other dependencies (e.g., JASP).

Keeping the current naming scheme also has the benefit that it is not necessary that metaBMA is submitted before your package. You could disable unit tests with ci values other than 95% with:
if(packageVersion("metaBMA") < "0.6.4") ...

Anyways, I will of course still check everything and submit the revision asap. Thanks again for pointing out the issue.

strengejacke commented 3 years ago

This is how it looks (and works) with current dev of metaBMA.

library(metaBMA)
library(parameters)
data(towels)

set.seed(123)
m1 <- meta_bma(logOR, SE, study, data = towels, ci = 0.95)
model_parameters(m1)
#> Parameter | Coefficient |   SE |        95% CI |    BF |  Rhat |     ESS
#> ------------------------------------------------------------------------
#> averaged  |        0.20 | 0.09 | [ 0.03, 0.37] |       |       |        
#> fixed     |        0.21 | 0.07 | [ 0.07, 0.36] | 11.97 | 1.000 | 3712.10
#> random    |        0.18 | 0.11 | [-0.03, 0.37] |  5.37 | 1.002 | 3670.80
#> 
#> Using highest density intervals as credible intervals.

set.seed(123)
m2 <- meta_bma(logOR, SE, study, data = towels, ci = 0.99)
model_parameters(m2)
#> Parameter | Coefficient |   SE |        99% CI |    BF |  Rhat |     ESS
#> ------------------------------------------------------------------------
#> averaged  |        0.20 | 0.09 | [-0.06, 0.43] |       |       |        
#> fixed     |        0.21 | 0.07 | [ 0.02, 0.39] | 11.97 | 1.000 | 3712.10
#> random    |        0.18 | 0.11 | [-0.16, 0.44] |  5.37 | 1.002 | 3670.80
#> 
#> Using highest density intervals as credible intervals.

set.seed(123)
m3 <- meta_random(logOR, SE, study, data = towels, ci = 0.95)
model_parameters(m3, ci_method = "eti")
#> Parameter | Coefficient |   SE |        95% CI |   BF |  Rhat |     ESS
#> -----------------------------------------------------------------------
#> Overall   |        0.18 | 0.10 | [-0.04, 0.37] | 2.00 | 1.000 | 3974.80
#> tau       |        0.14 | 0.10 | [ 0.03, 0.40] |      | 1.000 | 3145.40
#> 
#> Using equal-tailed intervals as credible intervals.

set.seed(123)
m4 <- meta_random(logOR, SE, study, data = towels, ci = 0.99)
model_parameters(m4, ci_method = "eti")
#> Parameter | Coefficient |   SE |        99% CI |   BF |  Rhat |     ESS
#> -----------------------------------------------------------------------
#> Overall   |        0.18 | 0.10 | [-0.16, 0.43] | 2.00 | 1.000 | 3974.80
#> tau       |        0.14 | 0.10 | [ 0.02, 0.60] |      | 1.000 | 3145.40
#> 
#> Using equal-tailed intervals as credible intervals.

set.seed(123)
m5 <- meta_fixed(logOR, SE, study, data = towels, ci = 0.95)
model_parameters(m5)
#> Parameter | Coefficient |   SE |       95% CI |    BF
#> -----------------------------------------------------
#> Overall   |        0.21 | 0.08 | [0.06, 0.36] | 11.97
#> 
#> Using highest density intervals as credible intervals.

set.seed(123)
m6 <- meta_fixed(logOR, SE, study, data = towels, ci = 0.99)
model_parameters(m6)
#> Parameter | Coefficient |   SE |       99% CI |    BF
#> -----------------------------------------------------
#> Overall   |        0.21 | 0.08 | [0.01, 0.41] | 11.97
#> 
#> Using highest density intervals as credible intervals.

Created on 2020-10-26 by the reprex package (v0.3.0)

strengejacke commented 3 years ago

This is how it looks (and works) with current dev of metaBMA.

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

set.seed(123)
m1 <- suppressWarnings(meta_bma(logOR, SE, study, data = towels, ci = 0.95))
model_parameters(m1)
#> # 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
#> ------------------------------------------------------------------------
#> averaged  |        0.20 | 0.09 | [ 0.03, 0.37] |       |       |        
#> fixed     |        0.21 | 0.07 | [ 0.07, 0.36] | 11.97 | 1.000 | 3712.10
#> random    |        0.18 | 0.11 | [-0.03, 0.37] |  5.37 | 1.002 | 3670.80
#> 
#> Using highest density intervals as credible intervals.

set.seed(123)
m2 <- suppressWarnings(meta_bma(logOR, SE, study, data = towels, ci = 0.99))
model_parameters(m2)
#> # Studies 
#> 
#> Parameter                                           | Coefficient |   SE |        99% 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 |        99% CI |    BF |  Rhat |     ESS
#> ------------------------------------------------------------------------
#> averaged  |        0.20 | 0.09 | [-0.06, 0.43] |       |       |        
#> fixed     |        0.21 | 0.07 | [ 0.02, 0.39] | 11.97 | 1.000 | 3712.10
#> random    |        0.18 | 0.11 | [-0.16, 0.44] |  5.37 | 1.002 | 3670.80
#> 
#> Using highest density intervals as credible intervals.

set.seed(123)
m3 <- suppressWarnings(meta_random(logOR, SE, study, data = towels, ci = 0.95))
model_parameters(m3, ci_method = "eti")
#> # 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
#> -----------------------------------------------------------------------
#> Overall   |        0.18 | 0.10 | [-0.04, 0.37] | 2.00 | 1.000 | 3974.80
#> tau       |        0.14 | 0.10 | [ 0.03, 0.40] |      | 1.000 | 3145.40
#> 
#> Using equal-tailed intervals as credible intervals.

set.seed(123)
m4 <- suppressWarnings(meta_random(logOR, SE, study, data = towels, ci = 0.99))
model_parameters(m4, ci_method = "eti")
#> # Studies 
#> 
#> Parameter                                           | Coefficient |   SE |        99% 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 |        99% CI |   BF |  Rhat |     ESS
#> -----------------------------------------------------------------------
#> Overall   |        0.18 | 0.10 | [-0.16, 0.43] | 2.00 | 1.000 | 3974.80
#> tau       |        0.14 | 0.10 | [ 0.02, 0.60] |      | 1.000 | 3145.40
#> 
#> Using equal-tailed intervals as credible intervals.

set.seed(123)
m5 <- suppressWarnings(meta_fixed(logOR, SE, study, data = towels, ci = 0.95))
model_parameters(m5)
#> # 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
#> -----------------------------------------------------
#> Overall   |        0.21 | 0.08 | [0.06, 0.36] | 11.97
#> 
#> Using highest density intervals as credible intervals.

set.seed(123)
m6 <- suppressWarnings(meta_fixed(logOR, SE, study, data = towels, ci = 0.99))
model_parameters(m6)
#> # Studies 
#> 
#> Parameter                                           | Coefficient |   SE |        99% 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 |       99% CI |    BF
#> -----------------------------------------------------
#> Overall   |        0.21 | 0.08 | [0.01, 0.41] | 11.97
#> 
#> Using highest density intervals as credible intervals.

Created on 2020-10-26 by the reprex package (v0.3.0)

IndrajeetPatil commented 3 years ago

closed in https://github.com/danheck/metaBMA/commit/dfa24c6c274752ee662b524b4323f981ec95960f