[R-sig-ME] R-sig-mixed-models Digest, Vol 33, Issue 23
Iasonas Lamprianou
lamprianou at yahoo.com
Tue Sep 22 21:55:25 CEST 2009
Dear friends,
I am wondering if anyone has a reference, or could give us his/her own opinion, on whether the multilevel models (MLM) and Generalizability Theory (GT) should give the same results. In essence, the GT uses ANOVA techniques to compute the variance components. Therefore, I would expect the models to give slightly different results when our data design is heavily unbalanced. Anyone else with a more theoretical opinion on the issue? Also, can the GT considered to be a special case of a general MLM?
There is a relevant paper but is not very illuminative:
Estimating reliability of school-level scores using multilevel and generalizability theory models, Asia Pacific Education Review, Volume 10, Number 2 / June, 2009, by Min-Jeong Jeon, Guemin Lee , Jeong-Won Hwang and Sang-Jin Kang
Thank you for your response
P.S. When I say MLM, I mean a generalized mixed effects model such the ones we run with lmer
Dr. Iasonas Lamprianou
Assistant Professor (Educational Research and Evaluation)
Department of Education Sciences
European University-Cyprus
P.O. Box 22006
1516 Nicosia
Cyprus
Tel.: +357-22-713178
Fax: +357-22-590539
Honorary Research Fellow
Department of Education
The University of Manchester
Oxford Road, Manchester M13 9PL, UK
Tel. 0044 161 275 3485
iasonas.lamprianou at manchester.ac.uk
--- On Tue, 22/9/09, r-sig-mixed-models-request at r-project.org <r-sig-mixed-models-request at r-project.org> wrote:
> From: r-sig-mixed-models-request at r-project.org <r-sig-mixed-models-request at r-project.org>
> Subject: R-sig-mixed-models Digest, Vol 33, Issue 23
> To: r-sig-mixed-models at r-project.org
> Date: Tuesday, 22 September, 2009, 3:29 PM
> Send R-sig-mixed-models mailing list
> submissions to
> r-sig-mixed-models at r-project.org
>
> To subscribe or unsubscribe via the World Wide Web, visit
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> or, via email, send a message with subject or body 'help'
> to
> r-sig-mixed-models-request at r-project.org
>
> You can reach the person managing the list at
> r-sig-mixed-models-owner at r-project.org
>
> When replying, please edit your Subject line so it is more
> specific
> than "Re: Contents of R-sig-mixed-models digest..."
>
>
> Today's Topics:
>
> 1. ranef, plot and dotplot in lmer and
> lme (Dunbar, Michael)
> 2. Re: [R-sig-eco] Problem with
> correlation structure in lme
> (Jonathan Myers)
> 3. Re: [R-sig-eco] Problem with
> correlation structure in lme
> (ONKELINX, Thierry)
>
>
> ----------------------------------------------------------------------
>
> Message: 1
> Date: Tue, 22 Sep 2009 11:31:08 +0100
> From: "Dunbar, Michael" <mdu at ceh.ac.uk>
> Subject: [R-sig-ME] ranef, plot and dotplot in lmer and
> lme
> To: "r-sig-mixed-models at r-project.org"
> <r-sig-mixed-models at r-project.org>
> Message-ID:
> <9C29AFE27E3FBA4DB1742CC8E8F74ECD04682E346B at nerckwmb1.ad.nerc.ac.uk>
> Content-Type: text/plain
>
>
> I have a two level-model with a random intercept and two
> random slope components. I can fit this with lme and then do
> plot(ranef(model), layout=c(3,1,1)), and somewhere in the
> code, it recognises that the x-axis scales are different and
> does the plot accordingly, and as expected.
>
> With lmer, things seem more complicated. I can do:
>
> rr1 <- ranef(model, postVar=TRUE)
> dotplot(rr1,scales = list(x = list(relation = 'free'),
> layout=c(3,1,1)))[["groupingfactor"]]
>
> However it seems like the explicit specification of
> relation=free then causes layout to be ignored, generally
> leading to a 2x2 panel plot, which causes too much bunching
> in the groups to be able to read them, and generally not
> using the space on the page in an optimal manner.
>
> It is possible to go for
> dotplot(rr1 , layout=c(3,1,1)))[["groupingfactor"]]
> and get the layout right but the x-axis scales are then
> equal for each panel, which will force the random slopes to
> be plotted on the same scale as the random intercepts :-(
>
> Has anyone else also had these problems, and perhaps found
> a solution?
>
> regards
>
> Mike Dunbar
>
>
>
>
> --
> This message (and any attachments) is for the recipient
> ...{{dropped:10}}
>
>
>
> ------------------------------
>
> Message: 2
> Date: Tue, 22 Sep 2009 08:55:59 -0500
> From: Jonathan Myers <jmyer19 at lsu.edu>
> Subject: Re: [R-sig-ME] [R-sig-eco] Problem with
> correlation structure
> in lme
> To: "ONKELINX, Thierry" <Thierry.ONKELINX at inbo.be>
> Cc: r-sig-mixed-models at r-project.org
> Message-ID:
> <ac10b4f30909220655i4907ca78u3bf96ac64a2a22dd at mail.gmail.com>
> Content-Type: text/plain
>
> 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 at 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 at 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 at 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 at gmail.com
> [mailto:jonamyers at gmail.com]
> *Namens *Jonathan
> > Myers
> > *Verzonden:* maandag 21 september 2009 18:56
> > *Aan:* ONKELINX, Thierry
> > *CC:* r-sig-mixed-models at 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 at 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 at 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 at 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 at gmail.com
> [mailto:jonamyers at gmail.com]
> *Namens *Jonathan
> >> Myers
> >> *Verzonden:* vrijdag 18 september 2009 21:55
> >> *Aan:* Philip Dixon
> >> *CC:* r-sig-ecology at r-project.org;
> ONKELINX, Thierry; mdu at 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 at 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 at 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 at 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 at 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.
> >
> >
>
> [[alternative HTML version deleted]]
>
>
>
> ------------------------------
>
> Message: 3
> Date: Tue, 22 Sep 2009 16:29:22 +0200
> From: "ONKELINX, Thierry" <Thierry.ONKELINX at inbo.be>
> Subject: Re: [R-sig-ME] [R-sig-eco] Problem with
> correlation structure
> in lme
> To: <jmyer19 at lsu.edu>
> Cc: r-sig-mixed-models at r-project.org
> Message-ID:
> <2E9C414912813E4EB981326983E0A10406B95CD6 at inboexch.inbo.be>
> Content-Type: text/plain
>
> 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 at 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 at gmail.com
> [mailto:jonamyers at gmail.com]
> Namens Jonathan
> Myers
> Verzonden: dinsdag 22 september 2009 15:56
> Aan: ONKELINX, Thierry
> CC: r-sig-mixed-models at 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 at 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 at 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 at 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 at gmail.com
> [mailto:jonamyers at gmail.com]
> Namens
> Jonathan Myers
>
> Verzonden: maandag 21 september 2009
> 18:56
> Aan: ONKELINX, Thierry
> CC: r-sig-mixed-models at 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 at 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 at 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 at 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 at gmail.com
> [mailto:jonamyers at gmail.com]
> Namens Jonathan Myers
> Verzonden: vrijdag 18
> september 2009 21:55
> Aan: Philip Dixon
> CC: r-sig-ecology at r-project.org;
> ONKELINX, Thierry;
> mdu at 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 at 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 at 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 at 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 at 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.
>
> [[alternative HTML version deleted]]
>
>
>
> ------------------------------
>
> _______________________________________________
> R-sig-mixed-models mailing list
> R-sig-mixed-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
>
> End of R-sig-mixed-models Digest, Vol 33, Issue 23
> **************************************************
>
More information about the R-sig-mixed-models
mailing list