Testing moderator comparisons in bayesian meta-analysis

Interfaces brms
,
1

I am using brms to fit a bayesian multilevel meta-analysis, including a moderator model with a 3-level moderator variable (i.e., Object, Body, Face). ex_data.csv (6.4 KB)

str(df)
tibble [135 × 5] (S3: tbl_df/tbl/data.frame)
  study_index : int [1:135] 1 1 1 1 2 3 3 3 3 3 ...
  es_num      : int [1:135] 1 2 3 4 5 6 7 8 9 10 ...
  yi          : num [1:135] 1.849 0.395 0.845 1.2 1.337 ...
  sei         : num [1:135] 0.44 0.277 0.311 0.35 0.16 ...
  stimuli_type: Factor w/ 3 levels "object","body",..: 2 2 2 3 2 2 2 3 3 1 ...

To fit the moderator model, I used the following code. Note - I actually used different simulation parameters (i.e., 20,000 iterations, tree depth =15)

fit <- brm(yi | se(sei) ~ 0 + Intercept + stimuli_type + (1|study_index) + (1|es_num), 
           data = df, 
           prior = c(prior(normal(0, 1), class = b),prior(cauchy(0,.3), class = sd)), 
           sample_prior = "yes")

Which results in summary output that looks like this:

> summary(fit)
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: yi | se(sei) ~ 0 + Intercept + stimuli_type + (1 | study_index) + (1 | es_num) 
   Data: df (Number of observations: 135) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Group-Level Effects: 
~es_num (Number of levels: 135) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.47      0.05     0.38     0.58 1.00      548     1618

~study_index (Number of levels: 36) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.20      0.09     0.02     0.39 1.01      279      366

Population-Level Effects: 
                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept            0.26      0.14    -0.01     0.53 1.00      850     1550
stimuli_typebody     0.45      0.14     0.18     0.74 1.00      898     1424
stimuli_typeface     0.79      0.18     0.44     1.15 1.00     1163     1907

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.00      0.00     0.00     0.00 1.00     4000     4000

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

To visualize the posterior distributions of this model, I used the tidybayes package:

fit_draws <- spread_draws(fit, `b_.*` ,regex = TRUE) %>%
  mutate(Object = b_Intercept,
         Body = b_Intercept + b_stimuli_typebody,
         Face = b_Intercept + b_stimuli_typeface) %>%
  dplyr::select(-c(b_Intercept, b_stimuli_typeface, b_stimuli_typebody)) %>% 
  pivot_longer(cols = -c(.chain,.iteration,.draw)) %>% arrange(name,.draw)

fit_draws$name <- factor(fit_draws$name,levels = c("Object", "Body", "Face"))
ggplot(fit_draws, aes(fill = name, x = value, color = name)) +
  stat_slab(alpha = .5)

image
image547×527 28.5 KB

To test different hypotheses about these, I used the hypothesis function to first evaluate if they were all above zero:

hypothesis(fit, "Intercept > 0") # Object greater than 0?
hypothesis(fit, "Intercept + stimuli_typebody > 0") # Body greater than 0?
hypothesis(fit, "Intercept + stimuli_typeface > 0") # face greater than 0?

and then compared them to each other using:

hypothesis(fit, "Intercept + stimuli_typebody > Intercept") # Body greater than Object?
hypothesis(fit, "Intercept + stimuli_typeface > Intercept + stimuli_typebody") # Face greater than body?
hypothesis(fit, "Intercept + stimuli_typeface > Intercept") # face greater than Object?

Is this the appropriate way to compare levels of this moderator variable to each other? If not, how would I compare these to each other (i.e., pairwise comparisons)? The reason I ask is because there is some overlap in these distributions, but the output from the hypothesis functions show that the posterior probability is .99 for hypothesis(fit, "Intercept + stimuli_typeface > Intercept + stimuli_typebody"); however, according the visualization, it would not appear that 99 percent of the face (blue) distribution is greater than the body distribution (green).

Your help is much appreciated!

1 Like
2

Good question. I am a bit short on time, so just quick notes:

That’s a bit suspicious and makes me worry whether the model does not have some additional problem (if you can’t get good results without large treedepth it usually means something is not going great). But the Rhat and ESS are not completely terrible, so the model is at least somewhat trustworthy.

If you don’t have a bug in your code (which I don’t see), then this is IMHO sensible.

One reason the marginal distributions might be misleading is:
a) For overlap to happen you need e.g. “Object” in the upper tail and “Body” in the lower tail at the same time which has lower probability than the variables being in the tails individually
b) there might correlation between the posteriors of the variables. Notably, if you have strong positive correlation between two variables, the sign of their difference can be determined with high certainty, even if their marginal distributions overlap substantially.

Best of luck!