[R-sig-ME] Heterogeneous (co)variances in a bivariate model with MCMCglmm?
David Atkins
datkins at u.washington.edu
Thu Sep 6 12:25:36 CEST 2012
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]]
More information about the R-sig-mixed-models
mailing list