R-sig-mixed-models Digest, Vol 33, Issue 23
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 **************************************************