Using variance from a Jeffreys prior in a Prevalence meta-analysis

Question

I am doing a meta-analysis on a group of studies, where the observation is a rate (no. of successes / no of trials). Some of the studies have small sample sizes (n<10) and/or did not observe any successes. Because of this I did not do the typical normal approximation SE calculation on the proportion for inputting into the meta-analysis.

As an alternative I used Jeffreys Priors and calculated 95% Bayesian Credible intervals for each rate.

Then in stata using the metan function, I inputted the observed rate and two bayesian intervals. The function runs fine of course and I can generate estimates, but was curious if this is appropriate to do. I believe stata takes the confidence intervals and backwards calculates the SE/variance and uses this in weighting.

I was curious if anyone thought this is inappropriate to use bayesian intervals in the mixed effects Meta-analysis function?

Answer 1

Why invent new/untested methodology when there are appropriate methods for dealing with this situation? In particular, there is the binomial-normal model that does not invoke the usual normal approximation for the sampling distribution of the estimates and instead is directly based on the binomial distribution (the 'normal' part only comes in for the random effects used to model heterogeneity). See, for example:

Stijnen, T., Hamza, T. H., & Ozdemir, P. (2010). Random effects meta-analysis of event outcome in the framework of the generalized linear mixed model with applications in sparse data. Statistics in Medicine, 29, 3046–3067.

Here is an example using R:

library(metafor)
dat <- get(data(dat.pritz1997))
res <- rma.glmm(measure="PLO", xi=xi, ni=ni, data=dat)
summary(res)

This yields the output:

Random-Effects Model (k = 14; tau^2 estimator: ML)

  logLik  deviance       AIC       BIC      AICc  
-32.6054   30.1105   69.2109   70.3018   70.4890  

tau^2 (estimated amount of total heterogeneity): 0.5893
tau (square root of estimated tau^2 value):      0.7677
I^2 (total heterogeneity / total variability):   67.79%
H^2 (total variability / sampling variability):  3.10

Tests for Heterogeneity: 
Wld(df = 13) = 25.0979, p-val = 0.0224
LRT(df = 13) = 42.2324, p-val < .0001

Model Results:

estimate       se     zval     pval    ci.lb    ci.ub          
  1.3766   0.2884   4.7735   <.0001   0.8114   1.9418      *** 

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The analysis is done on the log-odds scale. For easier interpretation, one can back-transform the results:

predict(res, transf=transf.ilogit)

This yields:

   pred  ci.lb  ci.ub  cr.lb  cr.ub
 0.7984 0.6924 0.8746 0.4426 0.9518

So, an estimated improvement rate of almost 80% (CI: 69% to 87%). Due to considerable heterogeneity, the credible/prediction interval is much wider (44% to 95%).

If you want to stick to Stata, the following syntax should work:

xtmelogit xi || id:, binomial(ni)

(xi and ni are the number of cases and trials and id is a study id variable).

Here are the results from Stata using the same data:

Mixed-effects logistic regression               Number of obs      =        14
Binomial variable: ni
Group variable: id                              Number of groups   =        14

                                                Obs per group: min =         1
                                                               avg =       1.0
                                                               max =         1

Integration points =   7                        Wald chi2(0)       =         .
Log likelihood = -32.605433                     Prob > chi2        =         .

------------------------------------------------------------------------------
          xi |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _cons |     1.3766   .2883857     4.77   0.000     .8113745    1.941826
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
id: Identity                 |
                   sd(_cons) |   .7676753    .279763       .375816    1.568122
------------------------------------------------------------------------------
LR test vs. logistic regression: chibar2(01) =    12.12 Prob>=chibar2 = 0.0002

Except for the LRT (which is calculated in different ways in R and Stata), the results match up nicely.

Answer 2

Suppose you really want to go the route with Jeffreys prior. So, in a particular study, you have $y$ cases/successes out of $n$ trials. Assume that the random variable $y$ follows a Binomial distribution, so that: $$f(y,\pi) = {n \choose x} \pi^y (1-\pi)^{n-y},$$ where $\pi$ is the true probability of a success.

Now let's place Jeffreys' prior on $\pi$. It turns out that Jeffreys' prior in that case is the Beta distribution with parameters $a=1/2$ and $b=1/2$. And since the Beta distribution is conjugate in this case, the posterior is easy to obtain: It is Beta with parameters $\alpha = y + 1/2$ and $\beta = n - y + 1/2$.

Therefore, the posterior mean is $$E[\pi|y,a,b] = \frac{y + 1/2}{n + 1}$$ while the posterior variance is $$Var[\pi|y,a,b] = \frac{(y + 1/2)(n - y + 1/2)}{n(n+1)^2}.$$

So, you could compute these values for each study and feed them to the meta-analysis software of your choice (or take the square-root of the posterior variances if the software is asking you to supply the standard errors). However, this doesn't get your around the problem of normality. The usual inverse-variance method (or 'normal-normal' model) is based on the assumption that the sampling distributions are (at least approximately) normal. In this example, we know that the distributions are not normal, but in fact Beta. If the posterior mean is not too close to 0 or 1, approximating the posterior by a normal may be fine, but certainly not when you observe no or only very few successes.

For illustration, assume $n=10$ and let's examine what those posterior distributions look like for various values of $y$:

n <- 10
par(mfrow=c(3,3), mar=c(3,4,3,1))
for (y in c(0,1,2,3,5,7,8,9,10)) {
   hist(rbeta(100000, y+1/2, n-y+1/2), col="lightgray", xlab="", main=paste0("y = ", y), xlim=c(0,1), freq=FALSE)
   curve(dnorm(x, mean=(y+1/2)/(n+1), sd=sqrt((y+1/2)*(n-y+1/2)/(n*(n+1)^2))), from=0, to=1, add=TRUE, col="red", lwd=2)
}

So, unless $y$ is close to 5, the normal approximation doesn't work. So, using the binomial-normal model would probably be the better approach here.