[R-sig-ME] glmer vs. MCMCglmm

[David Duffy]

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