Heterogeneous (co)variances in a bivariate model with MCMCglmm?

Hi Vincent--

Hmm, we may need to let Jarrod sound off on this; you are fitting  
(what I would call) a stratified model for the random-effects, and I  
have never tried to do something similar for the residual variances.

Have you tried either of:

rcov =   ~ us(Year*trait):units  # possibly '-1 + Year*trait'

or

rcov =   ~ idh(Year*trait):units  # possibly '-1 + Year*trait'

and, I would think you would want to alter your R prior to have:

V = diag(4), nu = 3.022

perhaps as a starting place.

I *think* those should fit separate residual error terms for both  
outcomes by each year.  The former might be under-identified given  
your data (as it allows covariances among all residuals).

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

On 9/6/12 7:25 PM, Vincent Careau wrote:

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