[R-sig-ME] Heterogeneous (co)variances in a bivariate model with MCMCglmm?

[David Atkins]

Vincent--

I've been doing a bit of work on multivariate mixed models using 
MCMCglmm with some colleagues (actually, for a tutorial paper on that 
topic...).

The following will fit random intercepts and slopes for two outcomes but 
*not* allowing correlations between outcomes (note that "j" is the 
grouping variable in our data).  Data is for simple two treatment over 
time study (common in psychiatry / psychology) with two outcomes.

# random intercept and slopes, but not correlated across outcomes
prior2 <- list(R = list(V = diag(2), nu = 1.002),
                G = list(G1 = list(V = diag(2), nu = 1.002),
                         G2 = list(V = diag(2), nu = 1.002)))

mult.mcmc3 <- MCMCglmm(cbind(y, y2) ~ -1 + trait*(tx*time),
                        data = data.wide,
                        random = ~ us(at.level(trait, 1) +
                                   at.level(trait, 1):time):j +
                                   us(at.level(trait, 2) +
                                   at.level(trait, 2):time):j,
                        rcov = ~ us(trait):units,
                        family = c("gaussian","gaussian"),
                        verbose = TRUE,
                        nitt = 25000, burnin = 5000, thin = 20,
                        prior = prior2)


Fitting correlated effects is actually a bit less "wordy" in terms of 
syntax:

# unstructured matrix between outcomes
prior3 <- list(R = list(V = .2*diag(2), nu = 3),
                G = list(G1 = list(V = diag(c(.2,.2,.1,.1)), nu = 5)))

mult.mcmc3.1 <- MCMCglmm(cbind(y, y2) ~ -1 + trait*(tx*time),
                        data = data.wide,
                        random = ~ us(-1 + trait*time):j,
                        rcov = ~ us(trait):units,
                        family = c("gaussian","gaussian"),
                        verbose = TRUE,
                        nitt = 25000, burnin = 5000, thin = 20,
                        prior = prior3)

Note that we're using simulated data and have been futzing a bit with 
priors, so I would definitely *not* suggest "prior3" for off-the-shelf 
usage.

Hope that helps.

cheers, Dave


-- 
Dave Atkins, PhD
University of Washington
datkins at u.washington.edu
http://depts.washington.edu/cshrb/

August 1 - October 30:

Universitat Zurich
Psychologisches Institut
Klinische Psychologie
Binzmuhlestrasse 14/23
CH-8050 Zurich

+41 44 635 71 75

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