Wouldn’t it make sense to just pool the adjusted prevalance estimates since the uncertainty should take the population size into account? It is a bit tricky, but I would suggest a beta regression for pooling. So, you need to figure out the Beta distribution parameters for each adjusted prevalence case, in the parameterization that brms uses.
This is the parameterization that brms uses: \mu=\frac{a}{a+b} (between 0 and 1), and \phi=a+b (>0). Or in reverse a=\mu \phi, b=\phi(1-\mu). \phi is interpretable as the population size, in your case an adjusted population size.
So, we need to turn the point estimate and the standard error/confidence interval into a mean (which is the point estimate) and into dispersion (\phi). If the interval is symmetrical and relies on the normal approximation we can use your formula for \sigma (standard deviation) i.e. the standard error.
Some other useful identities:
\sigma = \frac{\sqrt{\frac{\alpha \beta}{\alpha+\beta+1}}}{\alpha+\beta}=\sqrt{1-\mu} \sqrt{\frac{\mu}{\phi+1}}
From this follows:
\phi = \frac{\mu}{\sigma^2}-\left(\frac{\mu}{\sigma}\right)^2-1
Then the model would look something like this:
brm(formula = bf(mean ~ 1 + (1|study), phi ~ 0 + offset(derived_phi)), data = dataset, family = Beta(link_phi = "identity")) -> m1
EDIT:
An additional step that need to be done is to derive the mean for the logit-normal for each draw, since the random intercept is assumed to be normally distributed on the logit scale, and that distribution doesn’t have an analytical solution. We can do that with something like this:
# Function to calculate the mean of the logit-normal distribution
mean_logit_normal <- function(mu, sigma) {
# Compute mean using numerical integration
# This is the by the way the law of the unconscious statistician (LOTUS)
# plogis is expit, i.e. inverse logit
return(integrate(function(x) plogis(x) * dnorm(x, mu, sigma), lower = -Inf, upper = Inf)$value)
}
# Calculate the mean of the logit-normal distribution for each pair of mu and sigma
tibble(mu = m0 %>% as_draws_rvars() %>% pluck("Intercept") %>% draws_of() %>% as.vector(),
sigma = m0 %>% as_draws_rvars() %>% pluck("sd_study__Intercept") %>% draws_of() %>% as.vector()) %>%
mutate(mean_Y = map2_dbl(mu, sigma, mean_logit_normal, .progress = T)) ->
parameters
parameters$mean_Y %>% rvar()
parameters$mean_Y then contains the posterior draws for the pooled estimate.