[R-meta] Assessing selection bias / multivariate meta-analysis

Dear James,

thank you so much for your detailed reply.

In fact, the negative effect size does not represent an improvement but a
deterioration. However, I would still need to change the code for fitting
the model as the problem of one-sided p-value and direction as described
in
your email remains. Unfortunately, I'm not sure which of your suggestions
(a
or b) I should go for. I have tried both (results below). Do you have
further recommendations?

Best,
Pia


a)
1.
# 3PSM-bootstrap
set.seed(20240916)

system.time(
  mod_3PSM_boot <- selection_model(
  data = DatMA,
  yi = ES_all_spv$yi*-1,
  sei = sei,
  cluster = id_database,
  selection_type = "step",
  steps = .025,
  CI_type = "percentile",
  bootstrap = "multinomial",
  R = 1999
  )
)
print(mod_3PSM_boot, transf_gamma = TRUE, transf_zeta = TRUE)

  param Est SE percentile_lower percentile_upper
  beta 0.702 0.159 0.3011 1.059
  tau2 0.325 0.186 0.0893 1.756
  lambda1 0.165 0.138 0.0244 0.772


2.
# 3PSM-bootstrap
set.seed(20240916)

system.time(
  mod_3PSM_boot <- selection_model(
  data = DatMA,
  yi = ES_LOR_spv$yi*-1,
  sei = sei,
  cluster = id_database,
  selection_type = "step",
  steps = .025,
  CI_type = "percentile",
  bootstrap = "multinomial",
  R = 1999
  )
)

print(mod_3PSM_boot, transf_gamma = TRUE, transf_zeta = TRUE)
  param Est SE percentile_lower percentile_upper
  beta 0.529 0.163 -2.11e-01 0.876
  tau2 0.210 0.119 3.07e-02 0.731
  lambda1 0.157 0.144 6.00e-17 0.697




b)
1.
system.time(
  mod_3PSM_boot <- selection_model(
  data = DatMA,
  yi = ES_all_spv$yi,
  sei = sei,
  cluster = id_database,
  selection_type = "step",
  steps = .975,
  CI_type = "percentile",
  bootstrap = "multinomial",
  R = 1999
  )
)

print(mod_3PSM_boot, transf_gamma = TRUE, transf_zeta = TRUE)
param Est SE percentile_lower percentile_upper
  beta -0.702 0.159 -1.0594 -0.301
  tau2 0.325 0.186 0.0893 1.756
  lambda1 6.073 5.092 1.2949 40.908



2.
# 3PSM-bootstrap
set.seed(20240916)
## Startpunkt des Zufallszahlengenerators = wird als Seed bezeichnet.

system.time(
  mod_3PSM_boot <- selection_model(
  data = DatMA,
  yi = ES_LOR_spv$yi,
  sei = sei,
  cluster = id_database,
  selection_type = "step",
  steps = .975,
  CI_type = "percentile",
  bootstrap = "multinomial",
  R = 1999
  )
)

print(mod_3PSM_boot, transf_gamma = TRUE, transf_zeta = TRUE)
param Est SE percentile_lower percentile_upper
  beta -0.529 0.163 -0.8756 2.11e-01
  tau2 0.210 0.119 0.0307 7.31e-01
  lambda1 6.363 5.848 1.4342 3.18e+15



On So, 1 Dez 2024 10:18:05 -0600
  James Pustejovsky <jepusto at gmail.com> wrote:

Hi Pia,

You can ignore the warning messages that you're getting.
(We haven't yet
worked out how to suppress them in the bootstrapping
code.)

Your code looks good to me except for one subtle but
potentially
consequential issue. Based on the multivariate summary
meta-analyses, it
looks like you have a strongly negative average effect
size. If the effects
are coded so that negative values represent improvement,
then this needs to
be taken into account when fitting the selection models.
The models
implemented in metaselection are based on one-side
p-values, where the null
is mu <= 0 and the alternative is mu > 0 (i.e., positive
effects are
improvements). We have not yet implemented an option to
change the
direction of the null and alternative hypotheses
(although it's high on my
to-do list). In the mean time, there are a few
alternative ways that you
could modify the existing code to fit a more appropriate
model. Either:
a) Recode effect sizes so that positive values
correspond to improvement
(i.e., take yi = -yi) and interpret the beta estimate
accordingly.
or
b) Change the threshold of the step function to steps =
.975 and interpret
the lambda (selection parameter) estimate as the
relative probability that
a statistically significant (and negative) effect size
is reported compared
to the probability that a non-significant or
counter-therapeutic effect
size is reported. For instance, if lambda = 4 then this
means that a
statistically significant, therapeutic effect is four
times more likely to
be reported than a non-significant or non-therapeutic
effect.
Again, you only need do (a) or (b) but not both. It's
possible that making
this change will shift the results in a meaningful way
because the revised
model will be capturing a more plausible form of
selective reporting.

James

On Thu, Nov 28, 2024 at 11:24?AM Pia-Magdalena Schmidt <
pia-magdalena.schmidt at uni-bonn.de> wrote:

Dear James & Wolfgang, dear all,


I have tried the metaselection package for some of my
meta-analyses and
would be very grateful if you could have a look at the
results to see if
it
looks plausible to you?


James, you implemented several models, did I understand
correctly that you
recommend using the modified 3PSM with bootstrapping?

Furthermore, I received warnings after fitting the
mod_3PSM_boot (see
below). May I ignore them?


Results from fitted models with rma.mv, cluster robust &
mod_3PSM
bootstrapping for two analyses:


1.
1.1 Results rma.mv
res_spv <- rma.mv(yi=ES_all_spv$yi, V, random = ~ 1 |
id_database/effect_id,
data = dat)


Multivariate Meta-Analysis Model (k = 45; method: REML)
  logLik Deviance AIC BIC AICc
-52.7333 105.4667 111.4667 116.8192 112.0667

Variance Components:
  estim sqrt nlvls fixed factor
sigma^2.1 0.0619 0.2488 30 no id_database
sigma^2.2 0.2018 0.4493 45 no id_database/effect_id

Test for Heterogeneity:
Q(df = 44) = 187.0060, p-val < .0001

Model Results:
estimate se zval pval ci.lb ci.ub
  -1.0496 0.1073 -9.7845 <.0001 -1.2599 -0.8394 ***


1.2 Results rma.mv cluster robust
res_spv.robust <- robust(res_spv, cluster = id_database,
clubSandwich =
TRUE)


Multivariate Meta-Analysis Model (k = 45; method: REML)
  logLik Deviance AIC BIC AICc
-52.7333 105.4667 111.4667 116.8192 112.0667

Variance Components:
  estim sqrt nlvls fixed factor
sigma^2.1 0.0619 0.2488 30 no id_database
sigma^2.2 0.2018 0.4493 45 no id_database/effect_id

Test for Heterogeneity:
Q(df = 44) = 187.0060, p-val < .0001

Number of estimates: 45
Number of clusters: 30
Estimates per cluster: 1-3 (mean: 1.50, median: 1)

Model Results:
estimate se? tval? df? pval? ci.lb? ci.ub?
  -1.0496 0.1071 -9.7997 24.28 <.0001 -1.2706 -0.8287
***
1) results based on cluster-robust inference (var-cov
estimator: CR2,
  approx t-test and confidence interval, df:
Satterthwaite approx)


1.3 Results mod_3PSM bootstrapping
set.seed(20240916)
system.time(
  mod_3PSM_boot <- selection_model(
  data = dat,
  yi = ES_all_spv$yi,
  sei = sei,
  cluster = id_database,
  selection_type = "step",
  steps = .025,
  CI_type = "percentile",
  bootstrap = "multinomial",
  R = 1999
  )
)
print(mod_3PSM_boot, transf_gamma = TRUE, transf_zeta =
TRUE)


user system elapsed
269.861 0.799 270.760

Warning messages:
1: In optimx.run(par, optcfg$ufn, optcfg$ugr,
optcfg$uhess, lower, :
  Hessian is reported non-symmetric with asymmetry ratio
1.24952908365386e-11
2: Hessian forced symmetric
3: In optimx.run(par, optcfg$ufn, optcfg$ugr,
optcfg$uhess, lower, :
  Hessian is reported non-symmetric with asymmetry ratio
5.66259792559018e-11
4: Hessian forced symmetric
5: In optimx.check(par, optcfg$ufn, optcfg$ugr,
optcfg$uhess, lower, :
  Parameters or bounds appear to have different
scalings.
  This can cause poor performance in optimization.
  It is important for derivative free methods like
BOBYQA, UOBYQA, NEWUOA.
6: In optimx.check(par, optcfg$ufn, optcfg$ugr,
optcfg$uhess, lower, :
  Parameters or bounds appear to have different
scalings.
  This can cause poor performance in optimization.
  It is important for derivative free methods like
BOBYQA, UOBYQA, NEWUOA.

print(mod_3PSM_boot, transf_gamma = TRUE, transf_zeta

= TRUE)
  param Est SE percentile_lower percentile_upper
  beta -1.06e+00 0.1703 -1.63e+00 -8.09e-01
  tau2 3.07e-01 0.1781 9.40e-02 2.20e+00
  lambda1 8.07e+13 0.0604 1.13e+12 4.87e+14


2.
2.1 Results rma.mv
res_LOR <- rma.mv(yi=ES_LOR_spv$yi, V, random = ~ 1 |
id_database/effect_id,
data = dat)


Multivariate Meta-Analysis Model (k = 19; method: REML)
  logLik Deviance AIC BIC AICc
-14.8934 29.7867 35.7867 38.4579 37.5010

Variance Components:
  estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 13 no id_database
sigma^2.2 0.1889 0.4347 19 no id_database/effect_id

Test for Heterogeneity:
Q(df = 18) = 88.4226, p-val < .0001

Model Results:
estimate se zval pval ci.lb ci.ub
  -0.8414 0.1257 -6.6954 <.0001 -1.0877 -0.5951 ***


2.2 Results rma.mv cluster robust
res_LOR.robust <- robust(res_LOR, cluster = id_database,
clubSandwich =
TRUE)


Multivariate Meta-Analysis Model (k = 19; method: REML)

  logLik Deviance AIC BIC AICc
-14.8934 29.7867 35.7867 38.4579 37.5010

Variance Components:
  estim sqrt nlvls fixed factor
sigma^2.1 0.0000 0.0000 13 no id_database
sigma^2.2 0.1889 0.4347 19 no id_database/effect_id

Test for Heterogeneity:
Q(df = 18) = 88.4226, p-val < .0001

Number of estimates: 19
Number of clusters: 13
Estimates per cluster: 1-3 (mean: 1.46, median: 1)

Model Results:
estimate se? tval? df? pval? ci.lb? ci.ub?
  -0.8414 0.1178 -7.1438 9.08 <.0001 -1.1075 -0.5754 ***
1) results based on cluster-robust inference (var-cov
estimator: CR2,
  approx t-test and confidence interval, df:
Satterthwaite approx)


2.3 Results mod_3PSM bootstrapping
set.seed(20240916)
system.time(
  mod_3PSM_boot <- selection_model(
  data = dat,
  yi = ES_LOR_spv$yi,
  sei = sei,
  cluster = id_database,
  selection_type = "step",
  steps = .025,
  CI_type = "percentile",
  bootstrap = "multinomial",
  R = 1999
  )
)
print(mod_3PSM_boot, transf_gamma = TRUE, transf_zeta =
TRUE)

user system elapsed
262.519 0.324 262.812
Warning messages:
1: In optimx.run(par, optcfg$ufn, optcfg$ugr,
optcfg$uhess, lower, :
  Hessian is reported non-symmetric with asymmetry ratio
9.13585399793946e-13
2: Hessian forced symmetric
3: In optimx.run(par, optcfg$ufn, optcfg$ugr,
optcfg$uhess, lower, :
  Hessian is reported non-symmetric with asymmetry ratio
1.06505934682683e-11
4: Hessian forced symmetric

print(mod_3PSM_boot, transf_gamma = TRUE, transf_zeta

= TRUE)
  param Est SE percentile_lower percentile_upper
  beta -8.13e-01 0.1227 -1.14e+00 -6.07e-01
  tau2 1.56e-01 0.0892 1.87e-02 3.67e-01
  lambda1 1.12e+13 0.0566 1.36e+11 2.31e+14


Best,
Pia

On Do, 21 Nov 2024 17:18:28 +0000
  Viechtbauer, Wolfgang (NP) via R-sig-meta-analysis
<r-sig-meta-analysis at r-project.org> wrote:

Thanks for providing additional info about
metaselection. I think this package will be a very
important tool in the meta-analytic toolbox!

I was inspired by your recent blog post:

https://jepusto.com/posts/beta-density-selection-models/

and added the truncated beta selection model also to
selmodel(). I took a very quick peak at your code and I
think you are using analytic gradients, which helps to
speed up / stabilize the model fitting. But it's nice

to

have two parallel implementations to cross-check the
results.

Best,
Wolfgang

-----Original Message-----
From: James Pustejovsky <jepusto at gmail.com>
Sent: Thursday, November 21, 2024 15:10
To: R Special Interest Group for Meta-Analysis
<r-sig-meta-analysis at r-
project.org>
Cc: Viechtbauer, Wolfgang (NP)
<wolfgang.viechtbauer at maastrichtuniversity.nl>
Subject: Re: [R-meta] Assessing selection bias /
multivariate meta-analysis

Was going to chime in about the metaselection package
(https://github.com/jepusto/metaselection)---it's

still

under development but
the core functionality and documentation is in place.
The package implements the
bootstrapped selection model as demonstrated in my

blog

post

(https://jepusto.com/posts/cluster-bootstrap-selection-model/),

but with a much
easier interface and faster calculation; it also
implements selection models
with cluster-robust standard errors, though these

seem

to be not as accurate as
bootstrapping. Folks are welcome to give the package

a

try and to reach out
with questions or potential bugs if you run into
anything. We are working on a
paper describing the methods implemented in the

package

and reporting pretty
extensive simulations about their performance.

My student Man Chen (now on the faculty at UT Austin)
has studied a whole bunch
of the available methods for selective reporting bias
correction, looking
specifically at how they perform in meta-analyses

with

dependent effect sizes,
and proposing adaptations of some of the methods to
better acknowledge
dependency. Our working paper (on this is
here: https://osf.io/preprints/metaarxiv/jq52s

Pia asked about a few other possible techniques:
- The Egger's test / PET-PEESE approach with
cluster-robust variance estimation
is reasonable but, as Wolfgang noted, it is not
specifically diagnostic about
missing studies vs. missing effects. If the effect

sizes

nested within a given
study tend to have similar standard errors, then it

will

mostly be picking up on
association between study sample size and study-level
average effect size. And
of course, it also has the limitation that this
small-study association can be
caused by things other than selective reporting.
- Mathur & Vanderweele's sensitivity analysis is

quite

useful, though it does
not provide an estimate of the severity of selective
reporting (instead, it
provides information about the degree of potential

bias

assuming a specific
level of selection).
- For 3PSM, the cluster-bootstrap technique

implemented

in the metaselection
package is a way to deal with dependent effects, so

it

is no longer necessary to
use ad hoc approaches like ignoring dependence,
aggregating to the study level,
or selecting a single effect per study.

James

On Thu, Nov 21, 2024 at 6:37?AM Viechtbauer, Wolfgang
(NP) via R-sig-meta-
analysis <mailto:r-sig-meta-analysis at r-project.org>
wrote:
And I just stumbled across this:

https://github.com/jepusto/metaselection

James, don't hide all your good work from us!

Best,
Wolfgang

-----Original Message-----
From: R-sig-meta-analysis

<mailto:r-sig-meta-analysis-bounces at r-project.org>
On Behalf

Of Viechtbauer, Wolfgang (NP) via

R-sig-meta-analysis

Sent: Thursday, November 21, 2024 13:21
To: R Special Interest Group for Meta-Analysis

<r-sig-meta-analysis at r-

http://project.org>
Cc: Viechtbauer, Wolfgang (NP)

<mailto:wolfgang.viechtbauer at maastrichtuniversity.nl>

Subject: Re: [R-meta] Assessing selection bias /

multivariate meta-analysis

Dear Pia,

Generally, I don't think there really is any method

that is going to be a
great

choice here. The 'Egger sandwich' (i.e., an Egger

type

regression model using

cluster-robust inference methods) is a decent

option,

since it logically

generalizes the standard Egger regression method to

this context, but it is

unclear what kind of bias/selection effect this may

pick up (missing studies,

missing estimates within studies, a combination

thereof).

Yes, for the 3PSM, you would have to either ignore

the

dependencies or select

one estimate per study (and maybe repeat the latter

a

large number of times
for

different subsets).

I assume you are familiar with these papers. If

not,

they are directly
relevant:

Rodgers, M. A., & Pustejovsky, J. E. (2021).

Evaluating meta-analytic methods
to

detect selective reporting in the presence of

dependent effect sizes.

Psychological Methods, 26(2), 141-160.

https://doi.org/10.1037/met0000300

Fern?ndez-Castilla, B., Declercq, L., Jamshidi, L.,

Beretvas, S. N., Onghena,

P., & Van den Noortgate, W. (2021). Detecting

selection bias in meta-analyses

with multiple outcomes: A simulation study. The

Journal of Experimental

Education, 89(1), 125-144.

https://doi.org/10.1080/00220973.2019.1582470

Nakagawa, S., Lagisz, M., Jennions, M. D.,

Koricheva,

J., Noble, D. W. A.,

Parker, T. H., S?nchez-T?jar, A., Yang, Y., &

O'Dea,

R. E. (2022). Methods for

testing publication bias in ecological and

evolutionary meta-analyses. Methods

in Ecology and Evolution, 13(1), 4-21.

https://doi.org/10.1111/2041-210X.13724

I think James is working on some methods related to

this topic:

https://jepusto.com/posts/cluster-bootstrap-selection-model/

Best,
Wolfgang

-----Original Message-----
From: R-sig-meta-analysis

<mailto:r-sig-meta-analysis-bounces at r-project.org>
On

Behalf

Of Pia-Magdalena Schmidt via R-sig-meta-analysis
Sent: Wednesday, November 20, 2024 21:58
To: mailto:r-sig-meta-analysis at r-project.org
Cc: Pia-Magdalena Schmidt

<mailto:pia-magdalena.schmidt at uni-bonn.de>

Subject: [R-meta] Assessing selection bias /

multivariate meta-analysis

Dear all,
Although this topic has been discussed several

times

and I read the archives

and referenced papers, I?m still not sure how to

assess and possibly correct

for selection bias in multivariate meta-analyses.

I used the metafor package and ran meta-analyses

with SMCC as effect size

(all studies used within-designs) and fitted

http://rma.mv models as several

studies report more than one effect size.

Furthermore, I used cluster-robust

methods to examine the robustness of the models.
For a subset of my data, I used meta-regressions

with one continuous

moderator.
All effect sizes are from published journal

articles. The range of included

studies is between 30 and 6 with a number of

effect

sizes between 45 and 10.

Since I want to take the dependencies into

account,

I would not use funnel

plots or trim and fill. I wonder if using Egger's

regression test adjusted

for http://rma.mv as well as PET-PEESE and

perhaps

the sensitivity analysis

suggested by Mathur & Vanderweele (2020) as well

as

3PSM would be a

reasonable way to go? Although the latter would

only

use one effect size per

study or an aggregated effect size, right?

I would be very grateful for any recommendations!
Best,
Pia

Below is an excerpt from my code:
ES_all <- escalc(measure="SMCC", m1i= m1i, sd1i=

sd1i, ni = ni, m2i= m2i,

sd2i= sd2i, pi= pi, ri = ri, data= dat)
V <- vcalc(vi=ES_all$vi, cluster=id_database, obs

=

effect_id, rho =0.605,

data=dat)
res <- http://rma.mv(yi=ES_all$yi, V, random = ~

1 |

id_database/effect_id,
data =

dat)
res.robust <- robust(res, cluster = id_database,

clubSandwich = TRUE)

# subset
res_LOR <- http://rma.mv(yi=ES_LOR$yi, V, random

= ~

1 |
id_database/effect_id,

mods = ~ dose, data = dat)

to: