Example from 5.1 of the MCMCglmm course notes

spending.R

library(MCMCglmm)

id<-gl(200,4)                  # 200 people recorded four times
av_wealth<-rlnorm(200, 0, 1)
ac_wealth<-av_wealth[id]+rlnorm(800, 0, 1)
# expected disposable incomes + some year to year variation

av_ratio<-rbeta(200,10,10)
ac_ratio<-rbeta(800, 2*(av_ratio[id]), 2*(1-av_ratio[id]))
# expected proportion spent on car + some year to year variation

y.car<-(ac_wealth*ac_ratio)^0.25 # disposable income * proportion spent on car
y.hol<-(ac_wealth*(1-ac_ratio))^0.25 # disposable income * proportion spent on holiday
Spending<-data.frame(y.hol=y.hol, y.car=y.car, id=id)

m5a.2 <- MCMCglmm(cbind(y.hol, y.car) ~ trait - 1, random = ~us(trait):id,
                  rcov = ~us(trait):units, data = Spending, family = c("gaussian",
                  "gaussian"), verbose = FALSE)

VCV-array.R

nitt <- dim(m5a.2$VCV)[1] # number of iterations
VCV <- array(as.vector(m5a.2$VCV), dim=c(nitt, 2, 2, 2),
             dimnames=list(itt=1:nitt, trait=c('y.hol', 'y.car'), trait=c('y.hol', 'y.car'), random=c('id', 'units')))

VCV.R

> head(m5a.2$VCV, 2)
Markov Chain Monte Carlo (MCMC) output:
Start = 3001 
End = 3021 
Thinning interval = 10 
     y.hol:y.hol.id y.car:y.hol.id y.hol:y.car.id y.car:y.car.id
[1,]     0.01625509    0.008000464    0.008000464     0.02351099
[2,]     0.01247109    0.008376878    0.008376878     0.02543445
[3,]     0.01807759    0.002596848    0.002596848     0.03515039
     y.hol:y.hol.units y.car:y.hol.units y.hol:y.car.units y.car:y.car.units
[1,]        0.07366457      -0.018532214      -0.018532214        0.08075113
[2,]        0.07256900      -0.016906268      -0.016906268        0.07619431
[3,]        0.07013422      -0.005917229      -0.005917229        0.07609141