mixed model for repeat obs (CL Pressland)
It gives me a sadistic pleasure that other people have similar headaches as me. A few comments: You can get estimates of the Year, Week and Site effects by extracting the random effects. random.effects(model1) This gives you an idea of the signifcances of random effects. Note that there are assumptions attached to your random effects in this kind of analysis. I noted also that you did not seem to specify a correlation structure of your random effects which will usually make the estimation of your fixed effects less reliable since your algorithm will assume no correlation between random effects (i.e. the number of butterflies in week 2 is completely independent of the number of butterflies in week 1 etc....) Considering that you have a yearly pattern the previous comment that it would be good to find a time series structure for your data seems appropiate. (I would, ideally, identify the actual correlation structure in the time series analysis programs and then introduce them into the lme). I had often more bad luck than good luck with the approach. There is another possibility to make the time series analysis for each series and compare parameters between sites. This is less elegant but might be easier to get to work. (I applied that approach in a paper by me in Tree Physiology in 2004). Note that time series analysis does not support missing or unequally spaced observations and usually at least 30 observations. Also, time series requires that your time series are stationary. Frank
r-sig-ecology-request at r-project.org wrote:
1. Re: mixed model for repeat obs (CL Pressland) 2. Re: mixed model for repeat obs (Peter Solymos) ---------------------------------------------------------------------- Message: 1 Date: Tue, 24 Mar 2009 15:35:27 +0000 From: CL Pressland <Kate.Pressland at bristol.ac.uk> Subject: Re: [R-sig-eco] mixed model for repeat obs To: r-sig-ecology at r-project.org Message-ID: <43704A29163BCC0D7C169BB5 at bio-mammal03.bio.bris.ac.uk> Content-Type: text/plain; charset=us-ascii; format=flowed I have a data set that is unbalanced and consists of: 67 SITEs measured over several YEARs every WEEK for butterflies (LEPS per m). I'm interested in the MANagement code assigned to each site, but I
have
also data on TEMPerature, average SUN and WIND. My guess is that a linear
mixed model would be most appropriate and have constructed this code:
model1<-lme(LEPS~MAN,random=~YEAR/WEEK|SITE)
The output gives me:
--------------------------------------------------------------------
Linear mixed-effects model fit by REML
Data: NULL
AIC BIC logLik
-37631.24 -37566.48 18824.62
Random effects:
Formula: ~YEAR/WEEK| SITE
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 5.875102e-03 (Intr) YEAR
YEAR 1.392439e-06 -0.164
YEAR:WEEK 5.068196e-07 0.531 0.301
Residual 3.532589e-02
Fixed effects: LEPS ~ MAN
Value Std.Error DF t-value p-value
(Intercept) 0.009866718 0.001428957 9793 6.904841 0.00
MAN 0.000028304 0.001127429 65 0.025105 0.98
Correlation:
(Intr)
MAN -0.685
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.70566579 -0.40089121 -0.18073723 0.05900735 19.16411466
Number of Observations: 9860
Number of Groups: 67
--------------------------------------------------------------------
This output confuses me greatly! I figure that this clearly means that
management has no effect on butterflies but how can I figure out what
effect SITE, YEAR and WEEK have on the data? Would I have to also include
them in the fixed effects side of the formula (I'm unsure if this is
allowed)? Also, how could I include my weather variables? Would they just
be placed on the fixed effect side of the formula? If they are correlated
(I expect the weather variables are) would I have to place them in an
interaction rather than separately?
Any help that would be given would be gratefully received!
Kate
------------------------------
Message: 2
Date: Tue, 24 Mar 2009 11:09:18 -0600
From: Peter Solymos <solymos at ualberta.ca>
Subject: Re: [R-sig-eco] mixed model for repeat obs
To: CL Pressland <Kate.Pressland at bristol.ac.uk>
Cc: r-sig-ecology at r-project.org
Message-ID:
<34896f090903241009n77b40848hd6d700a1eda598a at mail.gmail.com>
Content-Type: text/plain; charset=UTF-8
Hi Kate,
You can use time series analysis (ar, arima functions at first)
instead, because YEAR and WEEK clearly has structure (i.e.
observations are conditional on previous observations with some lag).
To control for SITE, you can use polynomials of the geographical
coordinates (or write a hierarchical space-time series model say in
WinBUGS). Anyway, you should decide what are your parameters of
interest, and what you consider as nuisance parameters, and formulate
the model accordingly.
Cheers,
P?ter
P?ter S?lymos, PhD
Postdoctoral Fellow
Department of Mathematical and Statistical Sciences
University of Alberta
Edmonton, Alberta, T6G 2G1
Canada
email <- paste("solymos", "ualberta.ca", sep = "@")
On Tue, Mar 24, 2009 at 9:35 AM, CL Pressland
<Kate.Pressland at bristol.ac.uk> wrote:
I have a data set that is unbalanced and consists of: 67 SITEs measured over several YEARs every WEEK for butterflies (LEPS per m). I'm interested in the MANagement code assigned to each site, but I
have
also data on TEMPerature, average SUN and WIND. My guess is that a linear mixed model would be most appropriate and have constructed this code: model1<-lme(LEPS~MAN,random=~YEAR/WEEK|SITE) The output gives me: -------------------------------------------------------------------- Linear mixed-effects model fit by REML Data: NULL ? ? ? AIC ? ? ? BIC ? logLik ?-37631.24 -37566.48 18824.62 Random effects: Formula: ~YEAR/WEEK| SITE Structure: General positive-definite, Log-Cholesky parametrization ? ? ? ? ? ? ? ? StdDev ? ? ? Corr (Intercept) ? ? 5.875102e-03 (Intr) YEAR YEAR ? ? ? ? ? ?1.392439e-06 -0.164 YEAR:WEEK ? ? ? ? ? ? ? 5.068196e-07 ?0.531 ?0.301 Residual ? ? ? ?3.532589e-02 Fixed effects: LEPS ~ MAN ? ? ? ? ? ? ? ? Value ? Std.Error ? DF ?t-value p-value (Intercept) 0.009866718 0.001428957 9793 6.904841 ? ?0.00 MAN ? ? ? ? ? ? 0.000028304 0.001127429 ? 65 0.025105 ? ?0.98 Correlation: ? ? ? (Intr) MAN -0.685 Standardized Within-Group Residuals: ? ? ? Min ? ? ? ? ?Q1 ? ? ? ? Med ? ? ? ? ?Q3 ? ? ? ? Max -2.70566579 -0.40089121 -0.18073723 ?0.05900735 19.16411466 Number of Observations: 9860 Number of Groups: 67 -------------------------------------------------------------------- This output confuses me greatly! I figure that this clearly means that management has no effect on butterflies but how can I figure out what
effect
SITE, YEAR and WEEK have on the data? Would I have to also include them in the fixed effects side of the formula (I'm unsure if this is allowed)?
Also,
how could I include my weather variables? Would they just be placed on the fixed effect side of the formula? If they are correlated (I expect the weather variables are) would I have to place them in an interaction rather than separately? Any help that would be given would be gratefully received! Kate
_______________________________________________ R-sig-ecology mailing list R-sig-ecology at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
------------------------------
_______________________________________________ R-sig-ecology mailing list R-sig-ecology at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-ecology End of R-sig-ecology Digest, Vol 12, Issue 8 ********************************************