[R-sig-ME] correct specification and interpretation of multivariate mixed-model in lme

[David Villegas Ríos]

Dear list,
I have a dataset of 8 response variables measured every month during a
period of 3 years for 290 individuals. Not all individuals have the same
number of replicates (range=3-18, mean=8). I'm trying to specify a
multivariate mixed-model using lme because *I'm interested in getting the
among-individual variace-covariance matrix.* I have made some attempts in
MCMCglmm using the following model:

mcmc1=MCMCglmm(cbind(t1,t2,t3,t4,t5,t6,t7,t8)~trait-1,random=~us(trait):id,rcov=~us(trait):units,data=wide2,

 family=c("gaussian","gaussian","gaussian","gaussian","gaussian","gaussian","gaussian","gaussian"),verbose=FALSE,prior=prior1,pr=TRUE)

However, there is strong temporal autocorrelation inside each individual
time series, and this is reflected in the residuals of the MCMCglmm model.
So my only alternative is to use lme (I guess). The univariate lme models
for each trait separately included an autocorrelation term (corAR1 in most
cases) which completely accounted for the temporal autocorrelation. So I
wanted to use the same approach for the multivariate model. Nevertheless I
haven't found much information on how to fit a multivariate lme. I have
checked the following sources:

   -
   http://stats.stackexchange.com/questions/108779/multivariate-model-in-lme-with-independent-random-effect-similar-to-mcmcglmm
   -
   http://scs.math.yorku.ca/index.php/Mixed_Models_with_R/Multivariate_Mixed_Models.Rmd
   -
   http://rstudio-pubs-static.s3.amazonaws.com/3336_03636030d93d47de9131e625b72f58c6.html
   - Alain Zuur book on "Mixed effects models and extensions in Ecology
   with R"


However I'm still not sure about how to specify the model I need.

This is the head() and str() of my data (which is already in long format):

> head(long)
           key trait      value   id month_year month year
1 2011-05 5491   t1   0.5039553 5491    2011-05     5 2011
2 2011-05 5492   t1   1.1132514 5492    2011-05     5 2011
3 2011-05 5493   t1   1.7398036 5493    2011-05     5 2011
4 2011-05 5494   t1  -0.1328723 5494    2011-05     5 2011
5 2011-05 5495   t1   0.5433441 5495    2011-05     5 2011
6 2011-05 5496   t1   0.8299277 5496    2011-05     5 2011

> str(long)
'data.frame': 17992 obs. of  7 variables:
 $ key       : Factor w/ 2249 levels "2011-05 5491",..: 1 2 3 4 5 6 7 8 9
10 ...
 $ trait     : Factor w/ 8 levels "t1","t2","t3",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ value     : num  0.504 1.113 1.74 -0.133 0.543 ...
 $ id        : Factor w/ 274 levels "5488","5489",..: 3 4 5 6 7 8 23 24 25
26 ...
 $ month_year: chr  "2011-05" "2011-05" "2011-05" "2011-05" ...
 $ month     : num  5 5 5 5 5 5 5 5 5 5 ...
 $ year      : Factor w/ 4 levels "2011","2012",..: 1 1 1 1 1 1 1 1 1 1 ...


Ideally, I want to fit exactly the same model I fit with MCMCglmm but
including the correlation term. This is my best attempt so far:

model=lme(value~trait-1, random=~trait-1|id,
weights=varIdent(form=~1|trait),
correlation=corAR1(form=trait~1|id),data=long)

My questions are:

   1. Correlation term: I need to account for temporal autocorrelation
   within each response variable and for each individual. Is it correctly
   specified? Or should I create a "timer" variable that specifies the
   temporal order of the data, and include it in the correlation term? Because
   right now I'm not telling the correlation term that the replicates follow a
   temporal order I think. How could I specify that? In the univariate models
   my correlation term was corAR1(form=~timer), being timer a dummy variable
   reflecting the "time sequence".
   2. My aim including the weights term is to account for different
   variances of the different response variables. I guess this is strictly
   neccesary, correct?
   3. There are no fixed effects in the model yet, so the fixed part is
   reduced to trait-1. Correct? In case I want to account for month variation
   within each reponse variable, how should I include it in the fixed part of
   the model: (trait*month-1)?

Of course I'm aware that the model might not converge with 8 response
variables (actually it doesn't...but maybe because it is not correctly
specified), but once I know how to specifiy what I want properly, I can try
with less variables (perhaps 3-4). However, one last very general question:

      4. In general, if a MCMCglmm model with 8 responses converge, should
I expect convergence as well in lme?

I'd trully appreciate your advise. Please let me know if providing the
dataset helps.

Many thanks,

David

	[[alternative HTML version deleted]]