Skip to content

[R-meta] sd of blups vs tau in RE model

15 messages · Reza Norouzian, James Pustejovsky, Wolfgang Viechtbauer +1 more

Original
Yefeng Yang
# 
Dear experts,


I kindly request your guidance in resolving my matter.

To simplify my query, I would like to focus solely on the random-effects model for now. My question pertains to the conceptual equivalence between the tau (standard deviation of the mean/overall effect) and the standard deviation of the study-specific effects within the studies included in a meta-analysis. Some researchers refer to these study-specific effects as BLUPs.

To explore this further, I conducted an analysis using a specific dataset in metafor. However, my findings seem to suggest otherwise. I present a reproducible example below for your reference:

# calculate es and var
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)

# RE model
res <- rma(yi, vi, data=dat)

# calculate blups
blup_value <- blup(res)

# test whether the overall effect from res is equal to the simple mean of blup_value
res$beta[1] # -0.7145323
mean(blup_value$pred) # -0.7145323

We see the two values are equal.

# test whether tau from res is equal to the var of blup_value
sqrt(res$tau2) # 0.5596815
sd(blup_value$pred) # 0.4816293

We see the two values are not equal and seem to have a big gap.
Any comments and insights?

Best,
Yefeng
James Pustejovsky
# 
Yes, a known property of empirical Bayes (BLUP) estimates is that their SD
is not equal to the (estimated) SD of the random effects. See the following
articles and references therein:

Louis, T. A. (1984). Estimating a population of parameter values using
Bayes and Empirical Bayes methods. Journal of the American Statistical
Association, 79, 393?398.

Howard S. Bloom, Stephen W. Raudenbush, Michael J. Weiss & Kristin Porter
(2016): Using Multisite Experiments to Study Cross-Site Variation in
Treatment Effects: A Hybrid Approach With Fixed Intercepts and a Random
Treatment Coefficient, Journal of Research on Educational Effectiveness.
https://doi.org/10.1080/19345747.2016.1264518

Personally, I don't think this is a big problem because the SD of the BLUPs
is not the same thing as the SD of the population distribution. If it's a
concern, you could try putting a kernel density smooth over the BLUPs. The
resulting density estimate will have larger variance than the BLUPs. (I've
wondered for a while whether there's a way to choose the smoothing
bandwidth to make it match the random effects SD, but haven't had time to
look into it.)

James

On Thu, Jun 29, 2023 at 7:10?AM Yefeng Yang via R-sig-meta-analysis <
r-sig-meta-analysis at r-project.org> wrote:

            

      
      
      
    

        

      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Reza Norouzian
    

    

      
      #
    

  


  

      
-- Yefeng,

Along the same lines, estimates of variance and covariance by the
model are obtained from the unconditional/marginal distribution of the
random effects independent of the observed effect sizes.

BLUPs, on the other hand, represent the mean of random effects
conditional on the observed effect sizes. Thus, if you hand-calculate
a statistic (sd etc.) from the BLUPs, it won't likely match that
estimated by the model.

-- Just curious, James,

In my field, sometimes there is an overlience on hand-calculating
certain statistics from BLUPs when those estimates can be more
accurately estimated by the model. I sometimes think it might be a bit
problematic.

For example, in the following paper (doi:10.1017/S027226312300027X),
authors calculated a split-half reliability by fitting two separate
models; each for one set of items and then correlated the relevant
BLUPs to arrive at their reliability estimate (~.84; relatively good
reliability).

If instead, they fit a single multivariate model, their estimate would
have been ~.47 (relatively low reliability).

Possibly, that could be an issue.


Reza

On Thu, Jun 29, 2023 at 9:04?AM James Pustejovsky via

        
R-sig-meta-analysis <r-sig-meta-analysis at r-project.org> wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Yefeng Yang
    

    

      
      #
    

  


  

      
Hi James,

Very much thank you for providing your insights, two papers, and a potential solution. Before carefully looking at the papers you provided, I still have two related questions hope you can help address a bit:


  1.   Would you like to explain a bit "the SD of the BLUPs is not the same thing as the SD of the population distribution". Conceptually, I thought they are the same - both representing the marginal distribution after accounting for sampling error.
  2.  If I use some dispersion-relevant stats derived from the distribution of BLUPs to represent the heterogeneity, is it reasonable?

Best,
Yefeng

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      James Pustejovsky
    

    

      
      #
    

  


  

          

      
      
      
    

        
population consists of effect sizes for all the people on the listserv. We
estimate a meta-analysis using effect size estimates from you, me,
Wolfgang, and Gerta. The BLUPs are then estimates of the true effects for
you, me, Wolfgang, and Gerta. They're not estimating a feature of the
population distribution, they're estimating a feature of the effects for
the observed meta-analysts. Thus, there's no reason that the SD of the
effects from you, me, Wolfgang, and Gerta should be exactly the same as the
SD of the population of all the meta-analysts on the listserv.

            

      
      
      
    

        
of what I was getting at with my comment about the kernel density estimator.

  
  


  


      

    

    

      
      

        


  

        

      Yefeng Yang
    

    

      
      #
    

  


  

      
Hi James,

Thanks for your clarification. Your explanations are very clear. Actually, the SD of BLUPs and tau will converge when the within-study replicates are getting large.

A quick update for you:  I investigated KDE a little bit using density() in R. It seems to work well with BLUPs - it can show the empirical distribution without relying on normal distribution assumption. There is an augment in density() that can be used to tune smoothing bandwidth, although this parameter can not be tuned directly in density().


Also appreciate your comments, Reza.

Best,
Yefeng

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Reza Norouzian
    

    

      
      #
    

  


  

      
Yefeng,

Another question/comment that I hear when I distinguish between
statistics (sd, correlation etc.) hand-calculated from BLUPs and and
their counterparts estimated directly by the model (which is not
conditioned on observed data) is that given that the BLUPs are
conditioned on observed data and don't necessarily reflect the
population, then would it be possible or make sense to, for example,
construct a confidence interval for a statistic (sd, correlation etc.)
hand-calculate from BLUPs.


On Thu, Jun 29, 2023 at 8:48?PM Yefeng Yang via R-sig-meta-analysis

        
<r-sig-meta-analysis at r-project.org> wrote:

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      James Pustejovsky
    

    

      
      #
    

  


  

          

      
      
      
    

        
Can you say more about this? Is this claim based on simulations or
something? I see the intuition, but it also seems like this property might
depend not only on the within-study replicates all being large, but also on
their _relative_ sizes.

  
  


  


      

    

    

      
      

        


  

        

      Yefeng Yang
    

    

      
      #
    

  


  

      
Hi both,

I happen to come across a paper, which can answer both of your comments.

Eq. 1 and the following Eqs. show the derivation of the equivalence mentioned by my earlier email.

Wang C C, Lee W C. A simple method to estimate prediction intervals and predictive distributions: summarizing meta$B!>(Banalyses beyond means and confidence intervals[J]. Research Synthesis Methods, 2019, 10(2): 255-266.


Best,
Yefeng

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      James Pustejovsky
    

    

      
      #
    

  


  

      
The approximations there are predicated on k (the number of studies) being
large enough that the estimated heterogeneity (tau-hat) converges to the
true heterogeneity parameter.

On Thu, Jun 29, 2023 at 11:29?PM Yefeng Yang <yefeng.yang1 at unsw.edu.au>
wrote:

            

      
      
      
    

        

  
  


  


      

    

    

      
      

        


  

        

      Yefeng Yang
    

    

      
      #
    

  


  

      
Exactly. I was meant to the number of studies. I have no idea why I typed within-study replicates. Sorry for the confusion.

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Wolfgang Viechtbauer
    

    

      
      #
    

  


  

      
Dear Yefeng,

This was actually recently discussed on this mailing list:

https://stat.ethz.ch/pipermail/r-sig-meta-analysis/2023-May/004627.html

It is NOT true that the variance of the BLUPs will be equal to or approximate tau^2 when k is large. As I explain in that post, one can decompose tau^2 into two parts by the law of total variance. The variance of the BLUPs is only one part of this. To demonstrate:

library(metafor)

tau2 <- .02
k <- 2500
vi <- runif(k, .002, .05)
yi <- rnorm(k, 0, sqrt(vi + tau2))

res <- rma(yi, vi)
res$tau2
blups <- ranef(res)

# variance of the BLUPs (way too small as this is essentially only var(E(u_i|y_i)))
var(blups$pred)

# by adding what is essentially E(var(u_i|y_i)), we get (approximately) tau^2
var(blups$pred) + mean(blups$se^2)

The larger tau^2 is relative to the sampling variances, the less relevant mean(blups$se^2) will be (try running the code above with tau2 <- .2). But still, the variance of the BLUPs will underestimate tau^2.

Best,
Wolfgang

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Yefeng Yang
    

    

      
      #
    

  


  

      
Dear Wolfgang,

Excellent!  Sorry that I did not realize this post was on the mailing list (the title of this post really hard let me relate it to my question - anyway sorry for repeating the question)

I am arguing a bit. As pointed out by Wolfgang, the total variance of the population can be decomposed into two parts:
tau^2 = var(u_i) = E(var(u_i|y_i)) + var(E(u_i|y_i))

If k is larger enough, the expectation of the variance of the random effects should be zero or trivial, assuming the random effects are normally distributed. So, the first part of the above formula E(var(u_i|y_i)) is decreasing to 0 as K -> infinite. This is my understanding. I might be silly.

A relevant question is do you think it is meaningful to use the distribution of BLUPs to calculate something like "the proportion of the true effects above a certain value" to represent the heterogeneity?

Let's use your numerical example to show my point:
Say we want to know the proportion of the true effects (which are denoted by BLUPs) above 0,  we get a proportion of 0.5046695. Can we say the data is moderately heterogeneous?
BLUP = 0
Z = (BLUP - res$beta[1]) / sqrt(res$tau2)
1 - pnorm(Z)

Best,
Yefeng

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Wolfgang Viechtbauer
    

    

      
      #
    

  


  

      
Please see below for my responses.

Best,
Wolfgang

            

      
      
      
    

        
But this is not true as demonstrated with the code I provided.

            

      
      
      
    

        
Please try out what happens when you change k to different values. It has no systematic effect on mean(blups$se^2).

            

      
      
      
    

        
I am a bit confused by your question, as you are not using the BLUPs to compute this; you are using the estimated coefficient and value of tau^2. Also, looking at the proportion larger than 0 isn't the best example because this will be ~50% eithre way. So let's compare:

pnorm(0.1, coef(res), sqrt(res$tau2), lower.tail=FALSE)
mean(blups$pred > 0.1)

As you will see, these are not the same. The latter yields a smaller proportion since the distribution of the BLUPs is not as wide as the estimated distribution of true effects based on the model estimates. We can seee this with:

hist(blups$pred, freq=FALSE, breaks=50)
curve(dnorm(x, coef(res), sqrt(res$tau2)), add=TRUE, lwd=3)

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Yefeng Yang
    

    

      
      #
    

  


  

      
Dear Wolfgang,

Good points. I really appreciate your clarification and example.

Best,
Yefeng