[R-sig-ME] glmer vs. MCMCglmm

[Jarrod Hadfield]

Hi,

Sorry for not replying to this earlier - I'm in the middle of a field  
season at the moment so don't have much spare time. Thankfully  
Shinichi Nakagawa sent me the relevant information for reconciling  
MCMCglmm's results with sib-pair/lmer etc in the context of binary data.

The intra class correlation should be calculated as

VA/(VA+VE+pi^2/3)

where VE is a user-defined (and usually fixed) residual variance. This  
proportion is invariant to the choice of VE and for the PlodiaRB  
example gives an intra-class correlation of 0.15 which is the same as  
other programs.

If VE is set to 0.0001 (as in David's second MCMCglmm example), the  
chain will mix very poorly and not give valid inferences. This should  
be clear from basic mcmc diagnostics including plotting.

See Browne (2005) Journal of the Royal Statistical Society A 168  
599-613 for more info.

Cheers,

Jarrod



Quoting David Duffy <David.Duffy at qimr.edu.au>:

> On Thu, 2 Apr 2009, Jarrod Hadfield wrote:
>
>> Hi,
>>
>> As David suggested, the problem with binary data is that any   
>> residual variance in the underlying probability cannot be observed,  
>>  unlike other types of data. For example, if data are generated   
>> according to the process
>>
>> y ~ binom(n=1, p=inv.logit(mu + e_i))
>>
>> then the data look the same irrespective of the variance of the e's:
>>
>> y1<-rbinom(1000, 1, inv.logit(rnorm(1000, 0, 1)))
>> y2<-rbinom(1000, 1, inv.logit(rnorm(1000, 0, 10)))
>>
>> Because there is no information on VE (residual variance) we have   
>> to make an assumption about its value.  In terms of the fit of the   
>> model to the data the choice of VE should make no difference   
>> because the parameter is redundant. I think (if some one could   
>> clarify this that would be great) that lmer fixes VE to zero, where  
>>  with MCMCglmm the user is free to fix VE at any value.  Other than  
>>  that I believe the models should be identical. However, as VE   
>> approaches zero in MCMCglmm mixing becomes  problem, and actually   
>> when VE=0 the chain no longer mixes at all.
>>
>
> I'm not sure if this is of general interest or not.  I have been
> playing with that plodiaRB example, which has the advantage of just
> being sibship/cluster data.
>
> random intercept model
>            intercept  V_RE    intraclass r   LRT V_RE=0  Wald test SD_RE=0
> lmer          -0.986  0.5602  0.145 (1)
> glmmML        -0.986  0.5602  0.145          35.80       33.87
> Sib-pair      -1.020  0.6085  0.155 (1)                  31.60
> MCMCglmm(2)   -1.173  0.8830  0.469                      34.57
> MCMCglmm(3)    0.023  0.3590
>
> ANOVA icc estimator           0.098
> Tarone test = 73.61
> Variance chi-square = 127.54 (df=48)
>
> (1) taken as V/(V+pi^2/3)
> (2) units set to 1
> (3) units set to 0.0001
>
> animal model
>            intercept  V_A     h2
> Sib-pair      -1.068  0.7956  0.195 (1)
> MCMCglmm(2)   -1.339  2.2174  0.679
> MCMCglmm(3)   -0.046  0.3744   . Solar MFT (4)   .     0.375   0.375
>
> (1) taken as VA/(VA+pi^2/3)
> (2) units set to 1
> (3) units set to 0.0001
> (4) Multifactorial threshold model fitted using the SOLAR package
>
> For the one traditional geneticists model, which uses the binary (ie Pearson)
> correlations and a linear model for probabilities, the heritability should
> be ~0.2.  Using the other, the MFT, which (essentially) fits to the
> tetrachoric correlations, it is 0.38 (the relative magnitudes are the
> usual pattern I see).  For nongeneticists, the heritability can be
> thought of as the predicted intraclass correlation for identical twins,
> here one on the observed scale, the other on the linear predictor or
> latent variable scale.
>
> Cheers, David Duffy.
>
> -- 
> | David Duffy (MBBS PhD)                                         ,-_|\
> | email: davidD at qimr.edu.au  ph: INT+61+7+3362-0217 fax: -0101  /     *
> | Epidemiology Unit, Queensland Institute of Medical Research   \_,-._/
> | 300 Herston Rd, Brisbane, Queensland 4029, Australia  GPG 4D0B994A v



-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.