Hi,

Is it possible to fit a bivariate model with heterogeneous variance
components in MCMCglmm? 

 

Let's say we measured two traits (x1 and x2) on 50 individuals, twice per
year over two years (2008 and 2009). CREATE THE DATA:

ID<-c(rbind(seq(1,50),seq(1,50),seq(1,50),seq(1,50)))

YEAR<-rep(2008:2009,100)

x1<-ID/10+rnorm(200, mean = 0, sd = 2)

tmp<-data.frame(ID, YEAR,x1)

tmp$ID<-as.factor(tmp$ID)

tmp$YEAR<-as.factor(tmp$YEAR)

plot(x1~ID,data=tmp)

tmp2008<-subset(tmp,YEAR=="2008")

tmp2009<-subset(tmp,YEAR=="2009")

tmp2008$x2<-rnorm(100, mean = 0, sd = 3)+tmp2008$x1

tmp2009$x2<-rnorm(100, mean = 0, sd = 3)-tmp2009$x1+2*(mean(tmp2009$x1))

DATA<-rbind(tmp2008,tmp2009)

plot(x1~x2,DATA, col=YEAR)

 

In section 3.4 of the MCMCglmm course notes ("Heterogenous Residual
Variance"), we can see how to estimate a separate variance for ID and
residuals according to YEAR in a univariate model:

pr = list(R = list(V = diag(2), nu = 1), G = list(G1 = list(V = diag(2),nu =
1)))

model<-MCMCglmm(x1~1,

random=~idh(YEAR):ID,

rcov=~idh(YEAR):units,

data=DATA, nitt = 130000, thin = 100, burnin = 30000,prior=pr)

 

We can also fit a bivariate model:

BI.model<-MCMCglmm(cbind(x1, x2)~1,random=~us(trait):ID,

                                      rcov=~us(trait):units,

                                      family=c("gaussian","gaussian"),

                                      prior=pr,data=DATA)

 

However, I cannot find how to combine the two models above, where we
estimate separate variance (and covariances) in a bivariate model. Is this
possible at all? For sake of comparison, the code for such a model in
ASRreml-R is:

 

BI.model<-asreml(cbind(x1,x2)~trait,
random=~at(YEAR):us(trait):ID,rcov=~at(YEAR):units:us(trait),data=DATA,maxit
er=500)

 

cheers,

Vincent

 


	[[alternative HTML version deleted]]

