mcmcglmm_block_diagonal.R

library(MCMCglmm)
library(mvtnorm)

# separate covariances for each sex
Rf <- matrix(c(1, .4, .4, 1), nrow=2)
Rm <- matrix(c(1, 0, 0, 1), nrow=2)

# Simulate correlated phenotypes
females <- data.frame(rmvnorm(100, mean=c(0, 0), sigma=Rf))
males <- data.frame(rmvnorm(100, mean=c(0, 0), sigma=Rm))

names(females) <- c('x1', 'x2')
names(males) <- c('x1', 'x2')

females$sex <- 'f'
males$sex <- 'm'

# create single data frame
dat <- rbind(females, males)
dat$sex<-as.factor(dat$sex)

# prior R has two elements for each sex.
prior <- list(R=list(R1=list(V=diag(2), nu=3),
                     R2=list(V=diag(2), nu=3)))

# use at.level to match each row with the correct residual structure
m1 <- MCMCglmm(cbind(x1, x2) ~ trait-1 + trait:sex,
               rcov=~us(trait:at.level(sex, 'f')):units +
                     us(trait:at.level(sex, 'm')):units,
               prior=prior,
               family=rep('gaussian', times=2),
               data=dat, 
               verbose=FALSE)
summary(m1)
##  Iterations = 3001:12991
##  Thinning interval  = 10
##  Sample size  = 1000 
## 
##  DIC: 1150.074 
## 
##  R-structure:  ~us(trait:at.level(sex, "f")):units
## 
##                                                   post.mean l-95% CI u-95% CI eff.samp
## x1:at.level(sex, "f"):x1:at.level(sex, "f").units    0.8796   0.6646   1.1311     1000
## x2:at.level(sex, "f"):x1:at.level(sex, "f").units    0.5034   0.2926   0.7336     1068
## x1:at.level(sex, "f"):x2:at.level(sex, "f").units    0.5034   0.2926   0.7336     1068
## x2:at.level(sex, "f"):x2:at.level(sex, "f").units    1.2701   0.8957   1.6189     1390
## 
##                ~us(trait:at.level(sex, "m")):units
## 
##                                                   post.mean l-95% CI u-95% CI eff.samp
## x1:at.level(sex, "m"):x1:at.level(sex, "m").units    1.0809  0.81450   1.3813     1004
## x2:at.level(sex, "m"):x1:at.level(sex, "m").units    0.1664 -0.06324   0.3826     1000
## x1:at.level(sex, "m"):x2:at.level(sex, "m").units    0.1664 -0.06324   0.3826     1000
## x2:at.level(sex, "m"):x2:at.level(sex, "m").units    1.2258  0.91329   1.5702     1100
## 
##  Location effects: cbind(x1, x2) ~ trait - 1 + trait:sex 
## 
##              post.mean  l-95% CI  u-95% CI eff.samp pMCMC
## traitx1      -0.072394 -0.243784  0.109528     1000 0.422
## traitx2       0.006103 -0.197400  0.209340     1000 0.980
## traitx1:sexm  0.188854 -0.081251  0.463722     1095 0.210
## traitx2:sexm  0.173623 -0.172538  0.457807     1000 0.292