[R-sig-ME] Help establishing mixed model equation for split, plot design
Paul Johnson
paul.johnson at glasgow.ac.uk
Wed Apr 4 11:14:34 CEST 2018
Hi,
> Sorry, lme only takes nested random effects
Other designs are possible, e.g. fitting crossed random effects in lme is possible but fiddly. Pinheiro & Bates (Mixed Effects Models in S and S-Plus, 2000, p163-167) show how to do this using pdBlocked.
Best wishes,
Paul
> On 4 Apr 2018, at 11:34, Samuel Knapp <samuel.knapp at tum.de> wrote:
>
> Hi Felipe,
>
>
> please also respond to the mailing list.
>
>
> Sorry, lme only takes nested random effects (see
> https://stats.stackexchange.com/questions/228800/crossed-vs-nested-random-effects-how-do-they-differ-and-how-are-they-specified).
> In order to use the rep:salinity effect, you could add a column for the
> main plots in your data set:
>
>
> sp.datos$rep_salinityF <- paste(sp.datos$rep,sp.datos$salinityF,paste="_")
>
> lmemodel <- lme(Plantsurvival ~ salinityF*immersionF*SpecTF, random=
> ~1|rep_salinityF, data=sp.datos)
>
>
> I hope, it works like this.
>
>
> I'm afraid, I don't fully understand your question regarding equal
> number of samples. Do you mean missing observations? In lme, you need to
> add na.action="na.omit", while lmer is automatically omitting missing
> observations.
>
> You're right, that the ANOVA is doing type III SS. If you're interested
> in other SS type, you could look at the Anova() function from the car
> package. Apparently it can also deal with mixed models. Furthermore, the
> new version of lmerTest can also handle other SS type. See
> https://github.com/runehaubo/lmerTestR/blob/master/pkg_notes/new_lmerTest.pdf
>
>
> Best,
>
> Samuel
>
>
>
> On 04/04/18 10:15, CALLEJA APESTEGUI, FELIPE FRANCISCO wrote:
>>
>> Hi Samuel,
>>
>>
>> first of all thanks a lot for the quick and clear response.
>>
>>
>> As you state, all models give the same results once they are
>> established following your examples. Now I have a clear view about
>> what is really doing the command and how to continue.
>>
>>
>> Just one quick comment. The line:
>>
>> *lmemodel <- lme(Plantsurvival ~ salinityF*immersionF*SpecTF, random=
>> ~**1|rep:salinityF, data=sp.datos)*
>>
>>
>> Throws the following error:
>>
>> *Error in getGroups.data.frame(dataMix, groups) : **invalid formula
>> for groups*
>>
>> I leave it here in case is interesting for you. I've decided to use
>> the lmer package because I have other variables that are not balanced,
>> so I prefer to use the same model in all cases, and that package is
>> working fine.
>>
>> Just one final question. Is there any consideration I should have when
>> working with variables without equal number of samples? I mean,
>> something specific in the command of the lmer I should be careful with
>> or investigate further on. Or does the function itself has all
>> considerations included. i see from the *anova(fit.okay3,
>> ddf="Kenward-Roger") *that the anova is made using the type III SS so
>> I wouldn't have to change that part. *
>> *
>>
>> Also, fyi,� the aov function in the case of not grouping by rep throws
>> the same result as the lmer.
>>
>> Again, I really appreciate your help, it has been a lifesaver.
>>
>> Best regards,
>>
>> *Felipe Calleja Ap�stegui*
>>
>> *Predoctoral researcher*
>>
>> *
>> *
>>
>> *Instituto de Hidr�ulica Ambiental "IH Cantabria"*
>>
>> C/ Isabel Torres, N� 15
>>
>> Parque Cient�fico y Tecnol�gico de Cantabria
>>
>> 39011 Santander (Espa�a)
>>
>> www.ihcantabria.es <http://www.ihcantabria.es/>
>>
>> Tel:� +34 942 20 16 16 Ext. 1153
>>
>> Fax: +34 942 26 63 61
>>
>> e-mail: f*elipe-francisco.calleja at alumnos.unican.es*
>>
>>
>> ------------------------------------------------------------------------
>> *De:* Samuel Knapp <samuel.knapp at tum.de>
>> *Enviado:* martes, 3 de abril de 2018 21:26:31
>> *Para:* r-sig-mixed-models at r-project.org; CALLEJA APESTEGUI, FELIPE
>> FRANCISCO
>> *Asunto:* Re: Help establishing mixed model equation for split, plot
>> design
>> Hi Felipe,
>>
>> if I understand the design correctly, it is a split-plot design with
>> salinityF as the whole-plot. Only, if you have arranged the boxes in a
>> way that replicate blocks are formed, you should include the rep effect
>> in a model. These replicate blocks could be, that you always have three
>> boxes with the three different salinity levels grouped together, e.g. on
>> one shelf. If you have simply replicated each box 5 times and they are
>> set up fully randomly, you should not include a repliate block.
>>
>> If you have no missing values, and you specify the right model, both
>> aov() and lme() should give you the same results. I would suggest the
>> following models:
>>
>>
>> # with replicate blocks
>>
>> aovmodel <- aov(Plantsurvival ~ salinityF*immersionF*SpecTF +
>> Error(rep/salinityF), data=sp.datos)
>>
>> lmemodel <- lme(Plantsurvival ~ salinityF*immersionF*SpecTF, random=
>> ~1|rep/salinityF, data=sp.datos)
>>
>> # or with lmer from lme4 package (library (lme4)
>>
>> lmermodel <- lmer(Plantsurvival ~ salinityF*immersionF*SpecTF
>> +(1|rep/salinityF), data=sp.datos)
>>
>> # without replicate blocks (I'm not completely sure, if aov will report
>> the exact same results here)
>>
>> aovmodel <- aov(Plantsurvival ~ salinityF*immersionF*SpecTF +
>> Error(rep:salinityF), data=sp.datos)
>>
>> lmemodel <- lme(Plantsurvival ~ salinityF*immersionF*SpecTF, random=
>> ~1|rep:salinityF, data=sp.datos)
>>
>> lmermodel <- lmer(Plantsurvival ~ salinityF*immersionF*SpecTF
>> +(1|rep:salinityF), data=sp.datos)
>>
>>
>> As rep/salinityF is simply a short form for rep+rep:salinityF the
>> difference will be that you have no rep block effect in the models
>> without replicate blocks.
>>
>>
>> If you have no missing values, the point about using a mixed model here
>> is less about variance estimation through a random effect, but
>> automatically getting the right F-Test of the main-plot effect in ANOVA
>> (which is the main part of doing a proper split-plot analysis!!!).
>>
>> My experience is that the anova() function from the lmerTest package
>> will give you the right denominator df when using the Kenward-Roger
>> method. You need to use lmer for this:
>>
>> library(lmerTest)
>>
>> anova(lmermodel,ddf="Kenward-Roger")
>>
>> For proper mean comparisons, I suggest to use the emmeans package, with
>> the Tukey test and letter display from the cld() function:
>>
>> library(emmeans)
>>
>> cld(emmeans(lmermodel,"salinityF") #replace salinityF by the factor you
>> are interested
>>
>>
>> Finally, if you have one or more missing values, aov() will return
>> strange anova tables, thus better use a mixed model!
>>
>>
>> Best regards,
>>
>> Samuel
>>
>> --
>> Samuel Knapp
>>
>> Lehrstuhl f�r Pflanzenern�hrung
>> Technische Universit�t M�nchen
>> (Chair of Plant Nutrition
>> Technical University of Munich)
>>
>> Emil-Ramann-Strasse 2
>> D-85354 Freising
>>
>> Tel. +49 8161 71-3578
>> samuel.knapp at tum.de
>> www.researchgate.net/profile/Samuel_Knapp
>> <http://www.researchgate.net/profile/Samuel_Knapp>
>>
>>
>>
>> On 03/04/18 19:16, r-sig-mixed-models-request at r-project.org wrote:
>>> Message: 1
>>> Date: Tue, 3 Apr 2018 07:53:15 +0000
>>> From: "CALLEJA APESTEGUI, FELIPE FRANCISCO"
>>> <felipe-francisco.calleja at alumnos.unican.es>
>>> To: "r-sig-mixed-models at r-project.org"
>>> ������� <r-sig-mixed-models at r-project.org>
>>> Subject: [R-sig-ME] Help establishing mixed model equation for split
>>> ������� plot design
>>> Message-ID:
>>>
>> <DB6P191MB0088D09AF911107B33C1DE5C87A50 at DB6P191MB0088.EURP191.PROD.OUTLOOK.COM>
>>>
>>> Content-Type: text/plain; charset="iso-8859-1"
>>>
>>> Hello,
>>>
>>>
>>> I'm looking for some help establishing a mixed model ANOVA using R,
>> for a split plot design I've made for an experiment of saltmarsh
>> germination. I'll explain as clear as possible the experimental design
>> and afterwards what I've done and my doubts. Hope you can help me. I
>> haven't worked very much with R so many of my doubts are about what
>> I'm "telling" it to do with one or other command.
>>>
>>>
>>> The experiment consists in meassuring the germination percentage of
>> one species in variable conditions of salinity, immersion time and
>> presence of other species. I sow seeds in soil core's, that are inside
>> plastic boxes that are filled periodically with saline water. There
>> are three factors: salinity (3 levels: 0, 5 and 18), immersion time (3
>> levels: 0, 20, 40%), species treatment (2 levels: Baccharis, Baccharis
>> + Juncus). Inside each box there are 6 cores that combine in a
>> complete random design the factors of immersion and species treatment.
>> The box is filled with water at one of the levels of salinity. Thus,
>> I'm using a split plot design with salinity as the whole plot factor,
>> and immersion and species treatment as the subplot factors. All
>> factors are considered fixed. Each box is repeated 5 times. Thus,
>> there are 15 boxes and 90 soil cores. The dependent variable is the
>> percentage of germination of the species Baccharis in each core.
>>>
>>>
>>> I have several doubts about how to analyze the array.
>>>
>>>
>>> 1 - As far as I understand, although I'm treating all my factors as
>> fixed, this is a mixed model because of the interaction between
>> subjects (the boxes I believe), and the "split plot nature" of the
>> array, right? In that sense, which function would be better to analyze
>> this, the aov of the stats package, the lme of the nlme package, the
>> lmer of lme4?
>>>
>>>
>>>
>>> �� 2 - I've had trouble calculating the degrees of freedom for the
>> residuals. The only reference I have is that the error of the whole
>> plot part should have 12 df's, and the within error should have 60
>> df's. With that reference I've established two possible R commands:
>>>
>>>
>>> fit.aov2 <- aov(Plantsurvival ~ salinityF*immersionF*SpecTF +
>> Error(rep:salinityF/immersionF:SpecTF), data=sp.datos)
>>>
>>>
>>> With rep being the number of repetition. This last one gives the 12
>> and 60 df's for the error terms.
>>>
>>>
>>> The other option is:
>>>
>>> fit.okay <- lme(Plantsurvival ~ salinityF*immersionF*SpecTF, random=
>> ~1|rep/salinityF, data=sp.datos)
>>>
>>> But in this last case, the df's are 8 and 60, which makes me suspect
>> maybe there is something wrong. But as I said, I haven't cleared my
>> head on which should be the correct df's.
>>>
>>>
>>> Questions: Is the aov line solving a mixed model adequate for my
>> design?,
>>>
>>>
>>>
>>> ���������������������� Is the lme line considering salinityF as a
>> random factor? If it is, how can I tell it to consider all factors as
>> fixed, but put salinity at the "higher level" of the whole plot and
>> the other ones in the "lower level" of the subplot?
>>>
>>>
>>> I hope the questions and the experimental array are clear. If there
>> is any doubt or need more information please let me know. I attach a
>> csv file with the data in case you want to see it.
>>>
>>>
>>> Finally, if is not too much to ask, I'm fairly new to the splitplot
>> anova's and R, so I would really appreciate if you could answer be
>> with as much detail as possible, to fully understand what's going on
>> and where to continue.
>>>
>>>
>>> Thanks a lot,
>>>
>>>
>>> Felipe Calleja Ap�stegui
>>>
>>> Predoctoral researcher
>>>
>>>
>>> Instituto de Hidr�ulica Ambiental "IH Cantabria"
>>>
>>> C/ Isabel Torres, N� 15
>>>
>>> Parque Cient�fico y Tecnol�gico de Cantabria
>>>
>>> 39011 Santander (Espa�a)
>>>
>>> www.ihcantabria.es<http://www.ihcantabria.es/>
>>>
>>> Tel:� +34 942 20 16 16 Ext. 1153
>>>
>>> Fax: +34 942 26 63 61
>>>
>>> e-mail: felipe-francisco.calleja at alumnos.unican.es
>>>
>>>
>>>
>>
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
More information about the R-sig-mixed-models
mailing list