[R-sig-ME] glmer vs. MCMCglmm
Jarrod Hadfield
j.hadfield at ed.ac.uk
Sat May 2 17:53:16 CEST 2009
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.
More information about the R-sig-mixed-models
mailing list