[R-sig-ME] correct specification and interpretation of multivariate mixed-model in lme
Thierry Onkelinx
thierry.onkelinx at inbo.be
Mon Jun 22 16:47:48 CEST 2015
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 op 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]]
>
> _______________________________________________
> R-sig-mixed-models op r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
More information about the R-sig-mixed-models
mailing list