correct specification and interpretation of multivariate mixed-model in lme

Dear David,

1. It's best to use a "timer" variable. You want either
corAR1(form=~timer|id) or corAR1(form=~timer|id:trait). Note that the
correlation only works on the residuals within the same grouping
level. Residuals among levels are assumed to be be independent;
2. You allow for a difference residual variance among trait. It's up
to you to decide whether this makes sense.
3. trait-1 is equivalent with a different intercept for each trait. It
makes sense to always include the interaction between trait and the
other covariates.
4. You'll have to try it.

Note that the INLA packages allows for multivariate mixed models too.
It allows for correlated random effects (rather than correlated
residuals). You can find INLA at R-inla.org

Best regards,

ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature
and Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium

To call in the statistician after the experiment is done may be no
more than asking him to perform a post-mortem examination: he may be
able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does
not ensure that a reasonable answer can be extracted from a given body
of data. ~ John Tukey


2015-06-22 11:16 GMT+02:00 David Villegas R?os <chirleu at gmail.com>:

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