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