[R-sig-ME] cross-sex genetic correlation

[Simona Kralj Fiser]

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

On 26 July 2017 at 14:42, Jarrod Hadfield <j.hadfield at ed.ac.uk> wrote:

> 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]]