[R-sig-ME] Heterogeneous (co)variances in a bivariate model with MCMCglmm?
Jarrod Hadfield
j.hadfield at ed.ac.uk
Tue Sep 11 18:29:02 CEST 2012
Hi,
The rcov function does not allow the specification of direct products,
unlike the random structure which does fit what you want.
Being able to specify:
rcov=~ us(at.level(YEAR,1):trait):units + us(at.level(YEAR,2):trait):units
would be useful, but I think I have been asked this before (where YEAR
was GENDER) and when I looked into the code I realised that it would
actually take quite a bit of work to make possible. Sorry.
Cheers,
Jarrod
Quoting David Atkins <datkins at u.washington.edu> on Fri, 07 Sep 2012
08:50:02 +0200:
>
> 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
>>
>>
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
More information about the R-sig-mixed-models
mailing list