Dear Thierry,
Thanks again for the help and feedback. I think I've finally settled on a
model I am happy with!

lme(x ~ (water+fire+seed+census)^3 + block, random = ~1 | plot/subplot, data
= spprich5, correlation = corAR1(form = ~year), weights = varIdent(form = ~1
| fire))

#note: census and year refer to the same thing (year), but census is coded
as a factor (for the fixed effects component of the model) and year is coded
as an integer (for the correlation structure)

All the best,

Jonathan

On Tue, Sep 22, 2009 at 9:29 AM, ONKELINX, Thierry <Thierry.ONKELINX@inbo.be
> wrote:

>  Dear Jonathan,
>
> Note that you can use year as a factor in your fixed effect. But it needs
> to be a number for the corAR1() structure. I would even recommend to include
> is as a factor in the fixed effects, unless you have sound evidence for a
> linear trend along year.
>
> You could something like this:
>
> lme(x ~ (water+fire+seed+factor(year))^3 + block, random = ~1 |
> plot/subplot, data = spprich5, correlation = corAR1(form = ~year), weights =
> varIdent(form = ~1 | fire))
>
> Or you could add year twice to the dataset. Once as a number and once as a
> factor.
>
> Note that I moved block to the fixed effects because you have only two
> levels. You won't get good variance estimates with only two levels.
>
> HTH,
>
> Thierry
>
>
> ----------------------------------------------------------------------------
> ir. Thierry Onkelinx
> Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
> Forest
> Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
> methodology and quality assurance
> Gaverstraat 4
> 9500 Geraardsbergen
> Belgium
> tel. + 32 54/436 185
> Thierry.Onkelinx@inbo.be
> www.inbo.be
>
> 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
>
>
>
>  ------------------------------
> *Van:* jonamyers@gmail.com [mailto:jonamyers@gmail.com] *Namens *Jonathan
> Myers
> *Verzonden:* dinsdag 22 september 2009 15:56
>
> *Aan:* ONKELINX, Thierry
> *CC:* r-sig-mixed-models@r-project.org
> *Onderwerp:* Re: [R-sig-eco] Problem with correlation structure in lme
>
> Dear Thierry,
> Thanks very much for your reply. I changed year from a factor to an integer
> and that change seems to have fixed the problem, although it's not clear to
> me why lme accepts an integer as a fixed effect (perhaps lme treats numeric
> variables as factors when included as fixed effects?). I also added
> heterogeneous variance structure for one of the fixed effects (fire
> treatment):
>
> model <- lme(x ~ (water+fire+seed+year)^3, random = ~1 |
> block/plot/subplot, data = spprich5, correlation = corAR1(form = ~year),
> weights = varIdent(form = ~1 | fire))
>
> The output from this model looks appropriate: I now have 53 denDF for
> whole-plot factors and interactions, 54 denDF for split-factors and
> interactions, 230 denDF for year, and 230 denDF for interactions including
> year; there are a total of 360 observations (120 subplots sampled per year x
> 3 years) in the data set. The correlation structure for the repeated
> measurements was:
>
> Correlation Structure: AR(1)
>  Formula: ~year | block/plot/subplot
>  Parameter estimate(s):
>       Phi
> 0.5232107
>
> Thanks again for the help!
>
> All the best,
>
> Jonathan
>
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> Jonathan A. Myers
> Department of Biological Sciences
> Division of Systematics, Ecology, and Evolution
> Louisiana State University
> Baton Rouge, LA 70803 USA
>
> E-mail: jmyer19@lsu.edu
> Telephone: 225-578-7567
> Fax: 225-578-2597
>
> Website: http://www.biology.lsu.edu/labpages/harmslab/jmyers/index.html
> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>
> On Tue, Sep 22, 2009 at 3:00 AM, ONKELINX, Thierry <
> Thierry.ONKELINX@inbo.be> wrote:
>
>>  Dear Jonathan,
>>
>> How is year coded? As a number or as a factor? I woudl expect is to be a
>> factor, since that would make sense in your fixed effects. corAR1() will
>> probably require a number.
>>
>> Hopefully that might do the trick.
>>
>> HTH,
>>
>> Thierry
>>
>>
>> ----------------------------------------------------------------------------
>> ir. Thierry Onkelinx
>> Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
>> Forest
>> Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
>> methodology and quality assurance
>> Gaverstraat 4
>> 9500 Geraardsbergen
>> Belgium
>> tel. + 32 54/436 185
>> Thierry.Onkelinx@inbo.be
>> www.inbo.be
>>
>> 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
>>
>>
>>
>>  ------------------------------
>>  *Van:* jonamyers@gmail.com [mailto:jonamyers@gmail.com] *Namens *Jonathan
>> Myers
>> *Verzonden:* maandag 21 september 2009 18:56
>> *Aan:* ONKELINX, Thierry
>> *CC:* r-sig-mixed-models@r-project.org
>> *Onderwerp:* Re: [R-sig-eco] Problem with correlation structure in lme
>>
>>   Dear Thierry and List Members,
>> Thank you very much for the help. I think I'm very close to having an
>> appropriate model, but I cannot figure out why I get the following error
>> message:
>>
>> MODEL:
>>  lme(species.richness ~ (water+fuel+seed+year)^3, random = ~1 |
>> block/plot/subplot, correlation = corAR1(form = ~year))
>>
>>  ERROR MESSAGE:
>> Error in attributes(.Data) <- c(attributes(.Data), attrib) :
>>   'names' attribute [240] must be the same length as the vector [0]
>>
>> I would greatly appreciate advice on how to specify the correlation
>> structure properly. The model works fine when I remove the optional
>> correlation argument.
>>
>> My overall goal is to have a model that accounts for repeated measurements
>> of the subplots in each of the three years in the data set (i.e., to avoid
>> pseudoreplication). I've included a description of the experiment below.
>>
>> Thanks very much,
>>
>> Jonathan
>>
>> My experiment consists of three categorical treatments (fire, water, seed)
>> arranged in a split-plot design. The fire treatment (2 levels) and water
>> treatment (3 levels) were assigned to plots, and the seed treatment (2
>> levels) was assigned to two subplots within each plot. There are 60 total
>> plots (120 total subplots), divided equally among two large blocks (30 plots
>> per block). I measured species richness in each subplot in three separate
>> years. I would like to test for main effects of the three treatments, a main
>> effect of year, and all 3-way interactions. To avoid pseudoreplication, I
>> would also like to account for repeated measurements of the subplots in each
>> of the three years.
>>
>> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>> Jonathan A. Myers
>> Department of Biological Sciences
>> Division of Systematics, Ecology, and Evolution
>> Louisiana State University
>> Baton Rouge, LA 70803 USA
>>
>> E-mail: jmyer19@lsu.edu
>> Telephone: 225-578-7567
>> Fax: 225-578-2597
>>
>> Website: http://www.biology.lsu.edu/labpages/harmslab/jmyers/index.html
>> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>>
>>
>>
>>
>> On Mon, Sep 21, 2009 at 4:50 AM, ONKELINX, Thierry <
>> Thierry.ONKELINX@inbo.be> wrote:
>>
>>>  Dear Jonathan,
>>>
>>> Since you have multiple measurements (from several years) in your
>>> subplots, it does makes sense to include that in the random effect.
>>>
>>> There are several things you can do with the information about the years.
>>> First you could assume that there is a different trend (random slope along
>>> year) in each plot and subpplot. Here is an example with some tay data.
>>>
>>>
>>> years *<-* rnorm(
>>> 3)
>>>
>>> sdPlot *<-* 2
>>>
>>> sdSubplot *<-*
>>> 1
>>>
>>> sdYear *<-*
>>> 0.5
>>>
>>> sdError *<-*
>>> 0.1
>>>
>>> dataset *<-* expand.grid(year =
>>> 1:3, plot = factor(LETTERS[1:26]), subplot = factor(letters[1:26]))
>>>
>>> dataset
>>> $Y *<-* with(dataset, years[year] + rnorm(length(unique(plot)), sd =
>>> sdPlot)[plot] + year * rnorm(length(unique(plot)), sd = sdYear)[plot] +rnorm(length(unique(subplot)), sd = sdSubplot)
>>> [subplot] + rnorm(nrow(dataset), sd = sdError))
>>>
>>> library(nlme)
>>>
>>> lme(Y
>>> ~ factor(year), random = ~year|plot/subplot, data = dataset)
>>>
>>> Or you can assume that the random slope is only at the plot level. Since
>>> you have only a few subplots per plot, that will be more reasonable.
>>>
>>> lme(Y
>>> ~ factor(year), random = list(plot = pdSymm(form = ~year), subplot =
>>> pdSymm(form = ~1)), data = dataset)
>>>
>>> lme(Y
>>> ~ factor(year), random = list(plot = pdDiag(form = ~year), subplot =
>>> pdSymm(form = ~1)), data = dataset)
>>>
>>> Another option is that you assume that the residuals are correlated in
>>> time (e.g. AR1 correlation).
>>>
>>> lme(Y
>>> ~ factor(year), random = ~1|plot/subplot, data = dataset, correlation =
>>> corAR1(form = ~year))
>>>
>>> And you can even combine both a random slope and a correlation structure.
>>>
>>> lme(Y
>>> ~ factor(year), random = list(plot = pdDiag(form = ~year), subplot =
>>> pdSymm(form = ~1)), data = dataset, correlation = corAR1(form = ~year))
>>>
>>>
>>>
>>> Note that in the SAS the correlation structure defines the correlation
>>> between the random effects. This is in nlme equivalent with the pdClasses
>>> (see ?pdClasses for more info). The correlation structures in in nlme work
>>> on the residuals. See ?corClasses
>>>
>>> HTH,
>>>
>>> Thierry
>>>
>>> PS. R-sig-mixed models is a better list for this kind of questions.
>>>
>>>
>>> ----------------------------------------------------------------------------
>>> ir. Thierry Onkelinx
>>> Instituut voor natuur- en bosonderzoek / Research Institute for Nature
>>> and Forest
>>> Cel biometrie, methodologie en kwaliteitszorg / Section biometrics,
>>> methodology and quality assurance
>>> Gaverstraat 4
>>> 9500 Geraardsbergen
>>> Belgium
>>> tel. + 32 54/436 185
>>> Thierry.Onkelinx@inbo.be
>>> www.inbo.be
>>>
>>> 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
>>>
>>>
>>>
>>>  ------------------------------
>>> *Van:* jonamyers@gmail.com [mailto:jonamyers@gmail.com] *Namens *Jonathan
>>> Myers
>>> *Verzonden:* vrijdag 18 september 2009 21:55
>>> *Aan:* Philip Dixon
>>> *CC:* r-sig-ecology@r-project.org; ONKELINX, Thierry; mdu@ceh.ac.uk
>>> *Onderwerp:* Re: [R-sig-eco] mixed effects model in lme
>>>
>>>   Dear List Members,
>>> Thank you very much for your helpful replies and advice.
>>>
>>> I would greatly appreciate any additional advice on how to rewrite my lme
>>> model to account for repeated measurements of the subplots in each of the
>>> three years. I imagine this would involve modifying the random-effects
>>> component of the model and/or using the "CorStruct" function to specify the
>>> correlation structure (e.g., compound symmetry, autoregressive-moving
>>> average [AR1], etc.) for the repeated measures. My current model is:
>>>
>>> model = lme(species.richness ~ (water+fuel+seed+year)^3, random = ~1 |
>>> block/plot)
>>>
>>> Note that the model now includes only 3-way interactions for the fixed
>>> effects (I removed the 4-way interaction). According to Crawley (2007, The R
>>> Book, pg. 632), it is not necessary to specify the smallest nested spatial
>>> scale (subplot) in the random-effects component of the model; i.e., "random
>>> = ~1 | block/plot" will produce results identical to "random = ~1 |
>>> block/plot/subplot."
>>>
>>> Can anyone provide advice for how to rewrite this model to account for
>>> repeated measurements of the subplots in the three years? Either code for R
>>> or code for SAS would be helpful at this point, as I may ultimately have to
>>> perform the analysis in SAS to figure out how to do it using lme in R.
>>>
>>> I've pasted my original message below that includes the description of
>>> the experiment.
>>>
>>> Thanks very much!
>>>
>>> Jonathan
>>>
>>>
>>> Original message:
>>>
>>>  Dear List Members,
>>>
>>> I am using a mixed-effects model in lme and would like to know whether I
>>> am using the proper structure for the random-effects component of the model.
>>> My experiment consists of three categorical treatments (fire, water, seed)
>>> arranged in a split-plot design. The fire treatment (2 levels) and water
>>> treatment (3 levels) were assigned to plots, and the seed treatment (2
>>> levels) was assigned to two subplots within each plot. There are 60 total
>>> plots (120 total subplots), divided equally among two large blocks (30 plots
>>> per block). I measured species richness in each subplot in three separate
>>> years. My goal is to test for main effects of the three treatments, a main
>>> effect of year, and all interactions. My current model consists of four
>>> factorial fixed effects (fire, water, seed, year) and 1 random effect
>>> (block), with plots nested within blocks (to account for the split-plot
>>> structure of the experiment):
>>>
>>> model = lme(species.richness ~ water*fuel*seed*year, random = ~1 |
>>> block/plot)
>>>
>>> The ANOVA output includes two denominator degrees of freedom (denDF): 53
>>> denDF for plot factors (fire, water, fire x water interaction) and 270 denDF
>>> for split-plot factors (everything else).
>>>
>>> I would greatly appreciate feedback as to whether the random-effects
>>> component of the model looks appropriate, and if not, how it should be
>>> modified.
>>>
>>> Thanks very much!
>>>
>>> Cheers,
>>>
>>> Jonathan
>>>
>>> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>>> Jonathan A. Myers
>>> Department of Biological Sciences
>>> Division of Systematics, Ecology, and Evolution
>>> Louisiana State University
>>> Baton Rouge, LA 70803 USA
>>> Telephone: 225-578-7567
>>> Fax: 225-578-2597
>>> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>> On Fri, Sep 18, 2009 at 9:25 AM, Philip Dixon <pdixon@iastate.edu>wrote:
>>>
>>>> I believe you have three levels of variation:
>>>>
>>>> between plots assigned to the same fire/water treatment
>>>> between subplots assigned to the same fire/water/seed treatment
>>>> and between measurements (fire/water/seed/year)
>>>>
>>>> Your models (and the previous replies) appear to be ignoring the
>>>> repeated
>>>> measurements on the same subplot.  You will probably need to explore
>>>> various
>>>> choices of correlation structure among years (e.g. Compound Symmetry =
>>>> split/split plot, ar(1), ...)
>>>>
>>>> One key is that tests of split plot factors have 270 error DF (denDF).
>>>>  You
>>>> only have 120 subplots.
>>>>
>>>> Also, I suggest you retain at least all the fire/water/seed interactions
>>>> in
>>>> the model.   Each combination of fire/water/seed represents a specific
>>>> "thing"
>>>> (some call it treatment, but treatment has many different meanings) done
>>>> to a
>>>> subplot.   The full model with all these interactions corresponds to a
>>>> cell
>>>> means model.  Marginal means (i.e. main effects) correspond to averages
>>>> of
>>>> cell means.  If you drop those interactions, you are assuming that the
>>>> true
>>>> interaction is zero, so those dropped interactions only inform you about
>>>> the
>>>> error variation.  I believe this is somewhat dangerous because tests of
>>>> interactions have low power (i.e. lower than tests of main effects).
>>>> I recognize that there are many other opinions here.
>>>>
>>>> The above argument does not apply to interactions with year because year
>>>> is
>>>> not randomly assigned to a subplot or plot.  Many different opinions.
>>>>
>>>> I prefer to fit models like this in SAS, because it is much easier to
>>>> get
>>>> meaningful estimates (i.e. not just tests) in SAS than in lme.
>>>>
>>>> Best wishes,
>>>> Philip Dixon
>>>>
>>>> _______________________________________________
>>>> R-sig-ecology mailing list
>>>> R-sig-ecology@r-project.org
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>>>>
>>>
>>>
>>>
>>> --
>>> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>>> Jonathan A. Myers
>>> Department of Biological Sciences
>>> Division of Systematics, Ecology, and Evolution
>>> Louisiana State University
>>> Baton Rouge, LA 70803 USA
>>>
>>> E-mail: jmyer19@lsu.edu
>>> Telephone: 225-578-7567
>>> Fax: 225-578-2597
>>>
>>> Website: http://www.biology.lsu.edu/labpages/harmslab/jmyers/index.html
>>> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>>>
>>> Druk dit bericht a.u.b. niet onnodig af.
>>> Please do not print this message unnecessarily.
>>>
>>> Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer
>>> en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is
>>> door een geldig ondertekend document. The views expressed in  this message
>>> and any annex are purely those of the writer and may not be regarded as stating
>>> an official position of INBO, as long as the message is not confirmed by a duly
>>> signed document.
>>>
>>>
>>
>>
>>
>> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>> Jonathan A. Myers
>> Department of Biological Sciences
>> Division of Systematics, Ecology, and Evolution
>> Louisiana State University
>> Baton Rouge, LA 70803 USA
>>
>> E-mail: jmyer19@lsu.edu
>> Telephone: 225-578-7567
>> Fax: 225-578-2597
>>
>> Website: http://www.biology.lsu.edu/labpages/harmslab/jmyers/index.html
>> ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
>>
>> Druk dit bericht a.u.b. niet onnodig af.
>> Please do not print this message unnecessarily.
>>
>> Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer
>> en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is
>> door een geldig ondertekend document. The views expressed in  this message
>> and any annex are purely those of the writer and may not be regarded as stating
>> an official position of INBO, as long as the message is not confirmed by a duly
>> signed document.
>>
>>
>
>
>
>
> Druk dit bericht a.u.b. niet onnodig af.
> Please do not print this message unnecessarily.
>
> Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer
> en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is
> door een geldig ondertekend document. The views expressed in  this message
> and any annex are purely those of the writer and may not be regarded as stating
> an official position of INBO, as long as the message is not confirmed by a duly
> signed document.
>
>


-- 
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Jonathan A. Myers
Department of Biological Sciences
Division of Systematics, Ecology, and Evolution
Louisiana State University
Baton Rouge, LA 70803 USA

E-mail: jmyer19@lsu.edu
Telephone: 225-578-7567
Fax: 225-578-2597

Website: http://www.biology.lsu.edu/labpages/harmslab/jmyers/index.html
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	[[alternative HTML version deleted]]


