[R-sig-ME] Heterogeneous (co)variances in a bivariate model with MCMCglmm?
Vincent Careau
vcareau at ucr.edu
Thu Sep 6 19:25:29 CEST 2012
Thanks Dave,
I looked at your coding and modified it to the fictive data set (see below)
to fit a model with different variances and covariances in 2008 and 2009.
The following coding seems to work fine:
pr <- 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)))
BI.model<-MCMCglmm(cbind(x1, x2) ~ -1 + trait,
data = DATA,
random = ~ us(at.level(YEAR,1):trait):ID +
us(at.level(YEAR,2):trait):ID,
rcov = ~ us(trait):units,
family = c("gaussian","gaussian"),
nitt = 25000, burnin = 5000, thin = 20,
prior = pr)
As we can see, the covariance between x1 and x2 is different in 2008 and
2009. However, this code specifies common residual variances and covariances
in both years. I tried to modify the "rcov=" to allow heterogenous
residuals, but could not figure out how to make it work in a bivariate
context. I tried the same structure than I used in the random statement:
rcov = ~ us(at.level(YEAR,1):trait):units+
us(at.level(YEAR,2):trait):units,
but I get the following error:
"Error in MCMCglmm(cbind(x1, x2) ~ -1 + trait, data = DATA, random =
~us(at.level(YEAR, :
R-structure does not define unique residual for each data point"
The following seems to work fine:
rcov = ~ idh(at.level(YEAR,1):trait):units+
idh(at.level(YEAR,2):trait):units,
but when I call the model with the summary function I get this message:
"Error in rep(rep(1:length(object$Random$nrt), object$Random$nrt),
object$Random$nfl^2) :
invalid 'times' argument"
I would be surprised this is an overfitting problem because each individual
is measured twoce per year (changing it to 4 tiomes per year yields that
same error). Could it be related to prior specification?
Cheers,
Vincent
-----Original Message-----
From: r-sig-mixed-models-bounces at r-project.org
[mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf Of David Atkins
Sent: Thursday, September 06, 2012 3:26 AM
To: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] Heterogeneous (co)variances in a bivariate model
with MCMCglmm?
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]]
_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
More information about the R-sig-mixed-models
mailing list