[R-sig-ME] glmer vs. MCMCglmm
David Duffy
David.Duffy at qimr.edu.au
Fri Apr 3 04:45:54 CEST 2009
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
More information about the R-sig-mixed-models
mailing list