Dear all,

I'm having trouble with pedigreemm and simulated binary data. From the 
simulation, I know the heritability should be 0.5 and MCMCglmm give me 
for this dataset 0.53. Yet, pedigreemm yields a h² of 0.22. I'm not sure 
of how much bias I should expect from a Gauss-Hermite quadrature 
algorithm, but this still annoys me somehow...

Thus I would like to check with people on this list if every step of my 
R script is correct.

Here is what my data looks like:*
 > head(data)
   animal phen
1      1    0
2      2    1
3      3    1
4      4    0
5      5    1
6      6    0*

They come along with a pedigree:
*> head(ped) > tail(ped)
   animal sire dam                     animal sire dam
1      1   NA  NA                    995     995  773 870
2      2   NA  NA                    996     996  758 819
3      3   NA  NA                    997     997  766 821
4      4   NA  NA                    998     998  846 854
5      5   NA  NA                    999     999  822 867
6      6   NA  NA                    1000   1000  766 821*

Now, in order to use pedigreemm, I have to transform my data.frame "ped" 
into a pedigree class variable "ped2":
*ped2<-pedigree(ped$sire,ped$dam,ped$animal)*

My model is a binomial using a probit link with only the additive 
genetic random effect (Va variance component):
*m1<-pedigreemm(phen~(1|animal),data=data,family="binomial"(link="probit"),nAGQ=10,REML=TRUE,pedigree=list(animal=ped2),verbose=TRUE)*

I get a warning: /Number of levels of a grouping factor for the random 
effects is *equal* to n, the number of observations/
Which makes sense since I have an effect for each individual, but this 
is normal with an animal model.

A summary of the model:
*Generalized linear mixed model fit by the adaptive Gaussian Hermite 
approximation
Formula: phen ~ (1 | animal)
    Data: data
   AIC  BIC logLik deviance
  1617 1627 -806.5     1613
Random effects:
  Groups Name        Variance Std.Dev.
  animal (Intercept) 0.28499  0.53385
Number of obs: 1000, groups: animal, 1000

Fixed effects:
             Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.23902    0.08041  -2.973  0.00295 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 *

Now, I can extract the additive variance Va:
*va<-VarCorr(m1)$animal[1]*
And, because the model of pedigreemm assume no residual variance and I 
have a probit link, I can calculate the heritability with:
*h<-va/(va+1)*

And h equals 0.22 (not only for this dataset, but for all of them I get

Thanks a lot for any help you could provide !!

Regards,
Pierre de Villemereuil.








	[[alternative HTML version deleted]]

