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

[David Atkins]

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