[R-sig-ME] MCMCglmm R-structure problem in heteroskedastic model
Jarrod Hadfield
j.hadfield at ed.ac.uk
Wed Feb 2 12:21:35 CET 2011
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>
More information about the R-sig-mixed-models
mailing list