[R-sig-ME] residual covariance in trivariate MCMCglmm model with different error families
claudia.kasper at iee.unibe.ch
claudia.kasper at iee.unibe.ch
Fri May 20 10:49:41 CEST 2016
Dear list members,
I have a question concerning a way to obtain the residual covariance from a trivariate model with different error families. I am using MCMCglmm to see if 3 traits (one count, two binary) correlate on a genetic, maternal, and individual level. The model contains six random effects and an inverse Wishart prior that distributes the variation across the variance components and the residual variation had produced the best results (in terms of convergence, chain length, etc). To run the model I fixed the residual variance in a way that the variances of the binary traits as well as their covariances are fixed, and only the unit variance of the count trait and its covariances are estimated (see prior below).
Now I wonder how I calculate the phenotypic correlations? My guess would be summing up all (co)variance matrices including the residual (co)variance matrix. I followed Nakagawa's & Schielzeth's (2010) guide for the residual variances for binary and count data that have to be calculated separately after running the model, and I assume the same applies to covariances between two binary or two count traits (I add (pi^2)/3 to the binomial and log(1/exp(intercept)+1) to the Poisson residual (co)variance)). However, I am not sure what I should do for covariances of count traits with binary traits? Can I use sqrt((pi^2)/3) * sqrt(log(1/exp(intercept)+1))?
Probably there is a more standard way to get the phenotypic correlations that I am not aware of. In this case I would be very grateful if somebody could point me towards web resources or published papers.
Below are the prior specifications and the model formula I used:
priorIW <- list(R = list(V = diag(3)/7, n = 4, fix=2), G = list(G1 = list(V = diag(3)/7, n = 4), G2 = list(V = diag(3)/7, n = 4), G3 = list(V = diag(3)/7, n = 4), G4 = list(V = diag(3)/7, n = 4), G5 = list(V = diag(3)/7, n = 4), G6 = list(V = diag(3)/7, n = 4)))
m.GenCorr.1 <- MCMCglmm(cbind(Df,Cl,Dg) ~ trait-1, random=~us(trait):animal + us(trait):DamID + us(trait):TIDec + us(trait):Group + us(trait):DIDec + us(trait):Pop, rcov=~us(trait):units, data=dataDefense, pedigree=ped, family=c("poisson","categorical","categorical"), verbose = F, prior=priorIW, nitt=10000000, burnin = 100000, thin = 5000)
Kind regards,
Claudia Kasper
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
C¦A¦T¦T¦C¦G¦A¦T¦T¦C¦T¦A¦A¦G¦G¦C¦A¦T¦¦T¦T¦T¦C¦G¦A¦T
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
~~ Claudia Kasper
~~ PhD candidate
~~ Division of Behavioural Ecology
~~ Institute of Ecology and Evolution
~~ University of Bern, Switzerland
><((°> <°))>< ><((°> <°))>< ><((°> <°))>< ><((°> <°))><
More information about the R-sig-mixed-models
mailing list