[R-sig-ME] MCMCglmm R-structure problem in heteroskedastic model

[Jarrod Hadfield]

Hi Paul,

When you fit an animal model, MCMCglmm introduces "phantom" parents  
with missing data. These phantom parents allow the inverse of A (the  
relationship matrix) to be a) calculated easily b) very sparse and c)  
permuted for solving the mixed model equations more quickly.  When an  
interaction is fitted between animal and a factor like asthma.family  
phantom records are set up for all asthma.family/animal  combinations,  
hence  6730 missing records in models 2 and 3. If the pedigree links  
do not exist between individuals phenotyped in the two groups, then  
there would be more efficient ways of setting up the model, but they  
are not implemented. When the number of missing records is large it  
can even be more efficient (in terms of effective sample size per  
second) to do away with phantom parents and use a denser inverse A  
(specify nodes="TIPS" in the call to MCMCglmm).
If some individuals in mfs.fam have not been phenotyped you can use  
prunePed to remove useless individuals (make sure make.base=TRUE).

In model3 the phantom parents have not been assigned an asthma.family  
because you have a simple random term ~animal.  When it comes to build  
the R-structure it does not "know" what to do with the individuals.  
These phantom parents could be assigned any asthma.family and it would  
make no difference to the posterior, but unfortunatley I did not think  
that far head when writing MCMCglmm. I will try and fix this in later  
versions. The work around is to add the phantom parents to mfs.fam,  
assigning them anything for the fixed effects but having NA for the  
response. It should then run as intended.

Cheers,

Jarrod







On 2 Feb 2011, at 01:11, Paul Johnson wrote:

> Hello,
>
> I'm using MCMCglmm to estimate heritability of lung function in a  
> data set with 1346 nuclear families, each family having 2 parents  
> and at least 1 offspring (total n=4777). I'm having trouble with a  
> particular model:
>
> ########################
>>    prior.model3<- 
>> list(R=list(nu=0.02,V=diag(2)),G=list(G1=list(nu=0.02,V=1)))
>>    model3<-
> +       MCMCglmm(fixed.form,
> +         random=~animal,
> +         pedigree=pedigree,
> +         rcov=~idh(asthma.family):units,
> +         prior=prior.model3,
> +         data=mfs.fam,nitt=nitt,thin=thin,burnin=burnin)
> Warning in MCMCglmm(fixed.form, random = ~animal, pedigree =  
> pedigree, rcov = ~idh(asthma.family):units,  :
>  some combinations in animal do not exist and 2692 missing records  
> have been generated
> Error in MCMCglmm(fixed.form, random = ~animal, pedigree = pedigree,  
> rcov = ~idh(asthma.family):units,  :
>  R-structure does not define unique residual for each data point
> ##########################
>
> What I'm trying to do is estimate two different residual variances,  
> one in families with at least 1 asthmatic (n=722 individuals) and  
> one in the rest (n=4055 individuals), but as you can see I get an  
> error. Can anyone see what the problem is?
>
> Similar models work fine, including
> 1. a model with two variances [random=~animal, rcov=~units]
> 2. a model with 4 variances [random=~idh(asthma.family):animal,  
> rcov=~idh(asthma.family):units].
> 3. a 3-variance model with a single residual variance and two animal  
> variances [random=~idh(asthma.family):animal, rcov=~units].
>
> (An odd thing about 2 and 3 above is that MCMCglmm warns me that  
> "6730 missing records have been generated" which is exactly 2.5  
> times the 2692 that I'd expect).
>
> Thanks for any help,
> Paul Johnson
>
> PS other details:
>
>> R.Version()$version.string
> [1] "R version 2.12.0 (2010-10-15)"
>
>> Sys.info()
>     sysname      release      version     nodename       
> machine        login         user
>   "Windows"      "7 x64" "build 7600"    "PAUL-PC"         
> "x86"       "Paul"       "Paul"
>
>> cbind(installed.packages()["MCMCglmm",])
>          [,1]
> Package   "MCMCglmm"
> LibPath   "C:\\Users\\Paul\\Documents/R/win-library/2.12"
> Version   "2.10"
> Priority  NA
> Depends   "tensorA, Matrix, coda, ape, corpcor"
> Imports   NA
> LinkingTo NA
> Suggests  "rgl, combinat, mvtnorm, orthopolynom"
> Enhances  NA
> OS_type   NA
> License   "GPL (>= 2)"
> Archs     "i386, x64"
> Built     "2.12.1"
>
>
> The University of Glasgow, charity number SC004401
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20110202/2435007d/attachment-0004.pl>