cross-sex genetic correlation

Dear Jarrod and Paul,



Thank you for your reply. We used the suggested prior and model
specifications, but we also LOG transformed our weight data (L). The new
values are mostly negative. We ran:



prior <- list(R=list(V=diag(2), nu=0.02), G=list(G1=list(V=diag(2), nu=2,
alpha.mu=c(0,0),alpha.V=diag(2)*1000)))

model14 <- MCMCglmm(L~sex, random=~us(sex):animal, rcov=~idh(sex):units,
prior=prior, pedigree=Ped, data=Data1, nitt=100000, burnin=10000, thin=10)



The resulting summary is:



*Iterations = 10001:99991*

* Thinning interval  = 10*

* Sample size  = 9000 *



* DIC: -466.781 *



* G-structure:  ~us(sex):animal*



*                 post.mean   l-95% CI u-95% CI eff.samp*

*sex1:sex1.animal  0.003846  4.515e-10 0.009540     4761*

*sex2:sex1.animal  0.001122 -6.715e-04 0.003216     1436*

*sex1:sex2.animal  0.001122 -6.715e-04 0.003216     1436*

*sex2:sex2.animal  0.002096  1.310e-11 0.004439     5447*



* R-structure:  ~idh(sex):units*



*           post.mean l-95% CI u-95% CI eff.samp*

*sex1.units  0.019094 0.012842 0.025643     5700*

*sex2.units  0.007019 0.004553 0.009551     6510*



* Location effects: L ~ sex *



*            post.mean l-95% CI u-95% CI eff.samp  pMCMC    *

*(Intercept)   -0.9866  -1.0150  -0.9599     9521 <1e-04 ****

*sex2          -0.2536  -0.2843  -0.2227     9000 <1e-04 ****

*---*

*Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1*





With LOG values, the ?subscript out of bounds? problem in gone, herit()
analyses run and the resulting heritability and correlation values are
reasonable. However, the HPD interval is extremely wide.



*> **herit14<-model14$VCV[,'sex1:sex1.animal']/(model14$VCV[,'sex1:sex1.animal']+model14$VCV[,'sex1.units'])*

*> **herit15<-model14$VCV[,'sex2:sex2.animal']/(model14$VCV[,'sex2:sex2.animal']+model14$VCV[,'sex2.units'])*

*> **mean(herit14)*

*[1] 0.1643911*

*> **mean(herit15)*

*[1] 0.226494*

*> **corr.gen <- model14$VCV[,
'sex1:sex2.animal']/sqrt(model14$VCV[,'sex1:sex1.animal']*model14$VCV[,'sex2:sex2.animal'])*

*> **mean(corr.gen)*

*[1] 0.4729393*

*> **HPDinterval(herit14)*

*            lower     upper*

*var1 2.149316e-08 0.3883343*

*attr(,"Probability")*

*[1] 0.95*

*> **HPDinterval(herit15)*

*            lower     upper*

*var1 1.539724e-09 0.4509762*

*attr(,"Probability")*

*[1] 0.95*

*> **HPDinterval(corr.gen)*

*          lower    upper*

*var1 -0.1849416 0.999439*

*attr(,"Probability")*

*[1] 0.95*



We are starting to run out of ideas on why this is happening or where the
problem lies. We?d appreciate any further advice!



Eva and Simona

Hi,

The second way is a *much* better way of doing it but should give the same
answer. However, in both cases the residual covariance is not identifiable
(no individual is both male and female) and so you should use idh rather
than us.

The "subscript out of bounds" error is to do with your code that
post-processes the model output not an issue with MCMCglmm. Probably you
have used the wrong names for the (co)variance components.

Also, you haven't passed the prior to MCMCglmm, nor is the prior a valid
one for the problem as it specifies scalar variances rather than 2x2
covariance matrices. You could try

prior2 <- list(R=list(V=diag(2), nu=0.02), G=list(G1=list(V=diag(2), nu=2,
alpha.mu=c(0,0),alpha.V=diag(2)*1000)))

Cheers,

Jarrod




On 26/07/2017 13:33, Simona Kralj Fiser wrote:

model <- MCMCglmm(W~sex, random=~us(sex):animal, rcov=~us(sex):units,
prior=prior2, pedigree=Ped, data=Data1, nitt=100000, burnin=10000,
thin=10)

--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.

Doc. Dr. Simona Kralj-Fi?er

Institute of Biology
Scientific Research Centre of the Slovenian Academy of Sciences and Arts
Novi trg 2, P. O. Box 306, SI-1001 Ljubljana.
Phone: ++38614706333; fax: ++38614257797

http://ezlab.zrc-sazu.si/

http://bijh.zrc-sazu.si/sl/sodelavci/simona-kralj-fi%C5%A1er-sl#v

	[[alternative HTML version deleted]]