Skip to content

Meta-analysis for heritability using MCMCglmm?

2 messages · Jackie Wood, Jarrod Hadfield

Original
Jackie Wood
# 
Hi Jarrod,

We were finally able to dig into the statistical analysis of our
heritability data (as a reminder, we are conducting a meta-analysis
investigating heritability in relation to population size) and of course a
few questions have come up.

You had mentioned examining the residuals for the model. This may seem like
a "beginner" question, how does one extract the residuals from MCMCglmm?
The residuals.mcmcglmm function does not work. Summary plots of the MCMC
parameter estimates appear to be roughly normally distributed, and the
traces seem fine. Is this what you were referring to? We ran an unweighted
analysis (heritability estimates without SEs) in both MCMCglmm and lmer to
see if they give similar results, and they are basically concordant as
well. However, the residual distribution for the models run in lmer
(treating heritability as gaussian) is slightly skewed. These are from an
unweighted analysis, though, as I have also read that lme4 is unsuitable
for conducting formal weighted meta-analyses.

We also had a few questions that we thought would be worth discussing about
some methodological issues relating to incorporating common estimates of
heritability in the literature. Bayesian methodologies have become
increasingly popular to use when estimating trait heritabilities, but
bayesian estimates do not provide typical standard error or variance
estimates, as parent-offspring/ANOVA/REML methods do. Published bayesian
heritability estimates typically only include asymmetric confidence
intervals, and we unsure whether these can be translated into variance
estimates that can be used to weight our meta-analysis. For now, we plan on
performing a weighted meta-analysis using heritability estimates that
provide S.E.s, and an additional unweighted analysis that will include the
bayesian point-estimates we have collected from the literature. We were
wondering if you resolved this issue in your own heritability meta-analysis
and knew of a way to incorporate bayesian estimates (which form a
considerable proportion of the suitable heritability estimates available,
at least in recent history) into a formal weighted meta-analysis.

Additionally, we were wondering about the suitability of DIC to conduct
model selection for our analysis of heritability. I recall reading on this
SIG list that you had mentioned that there were potential issues using DIC
for hierarchical models as well as non-gaussian data. Since we're treating
heritability as gaussian, would it still be appropriate?

Any advise would be much appreciated!
On Thu, Jan 15, 2015 at 3:41 PM, Jackie Wood <jackiewood7 at gmail.com> wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


    

    
12 days later

    

      
      

        


  

        

      Jarrod Hadfield
    

    

      
      #
    

  


  

      
Hi Jackie,

1/ For the residuals you can do:

my_data$y-predict(my_model)

In the code above the predictions are marginal with respect to the  
random effects. If you want the predictions to include the random  
effects then use:

my_data$y-predict(my_model, marginal=NULL)


2/ I don't have a good solution for the weighted analysis. In fact we  
didn't do a weighted analysis in the end because only about 50% of  
studies provided standard errors. We figured that we lost more power  
by discarding the 50% than we would gain by incorporating information  
about the relative precision of different estimates. REML standard  
errors are based on the sampling distribution of estimates tending to  
the normal as n goes to infinity. In reality the sampling distribution  
will be skewed (particularly when h2 and/or sample size is low) like  
the posterior distribution. Its far from perfect but you could use the  
standard error given by:

0.5*((u95-pm)/qnorm(0.975)+(l95-pm)/qnorm(0.025))

where l95 and u95 are the lower and upper 95% credible intervals and  
pm is the posterior mean/mode.  The motivation behind this is *if* the  
posterior was normal and the asymmetry in the 95% credible intervals  
was just due to Monte Carlo error then this would be the best estimate  
of the posterior standard deviation (sort of like a standard error)  
given the information at hand. Not great justification!


3/ DIC (as focussed in MCMCglmm) is not a reliable model selection  
criterion for scientific inference. For Gaussian resopnses you could  
refocus DIC at the correct level post-analysis and then it would be  
more reliable. I should probably write a function to do this ....

Cheers,

Jarrod


Quoting Jackie Wood <jackiewood7 at gmail.com> on Sun, 1 Feb 2015 12:10:51 -0500: