[R-sig-ME] correct specification and interpretation of multivariate mixed-model in lme
David Villegas Ríos
chirleu at gmail.com
Tue Jun 23 09:00:18 CEST 2015
Dear Thierry, thanks for the comments (I'll try what you suggest) and the
link to INLA, which I didn't know - but seems really useful for some of my
data.
David
2015-06-22 16:47 GMT+02:00 Thierry Onkelinx <thierry.onkelinx at inbo.be>:
> 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]]
> >
> > _______________________________________________
> > R-sig-mixed-models at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list