[R-sig-eco] mixed model for repeat obs (CL Pressland)

Mike Dunbar mdu at ceh.ac.uk
Thu Mar 26 15:11:23 CET 2009


>From a quick look only, I'm not convinced that the structure of the model as suggested by Kate, is correct.

There is only one random effect in the model, site. Year and Week are specified as random slopes. But they are not in the fixed part of the model, hence I don't think this makes sense. Also if butterflies are measured every week, then what does week represent. It can't be a specific random effect as there is no replication within week. 

The question is how to treat time in the model. If there is autocorrelation in time for the residuals then it would make sense to model  this. But lme can only estimate a single common AR parameter across all sites. If there are seasonal patterns then again it would make sense to include this, either using some kind of sine or cosine functions or simply a factor representing the seasons. 

It would clearly make sense to include the other relevant covariates (temp etc) and this might reduce the autocorrelation in the residuals. 

Mike


>>> Frank Berninger <frankberninger at gmail.com> 26/03/2009 13:53 >>>
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
> ********************************************
>

_______________________________________________
R-sig-ecology mailing list
R-sig-ecology at r-project.org 
https://stat.ethz.ch/mailman/listinfo/r-sig-ecology

-- 
This message (and any attachments) is for the recipient ...{{dropped:6}}



More information about the R-sig-ecology mailing list