[R-sig-ME] GLMM for Combined experiments and overdispersed data
Thierry Onkelinx
thierry.onkelinx at inbo.be
Mon Apr 24 16:46:13 CEST 2017
Dear Peter,
Both models will yield identical results in case tree_id uses unique codes
over the blocks.
Best regards,
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium
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
2017-04-24 15:55 GMT+02:00 Peter Claussen <dakotajudo op mac.com>:
> Juan,
>
> I would model this as
>
> m3 = glmer(resp ~ trt * farm + (1| bk/tree), family = binomial, data=df)
> or
> m3 = glmer(resp ~ trt * farm + (1| bk) + (1| tree_id), family =
> binomial, data=df)
> (I can’t say off the top of my head if what the difference would be if
> you’re dealing with over-dispersion).
>
> 1. I’m assuming that block is a somewhat uniform grouping of trees, so
> that including block in the model gives you an estimate of spatial
> variability in the response, and if that is important relative to
> tree-to-tree variation.
>
> 2. You will most certainly want to include trt*farm to test for
> treatment-by-environment interaction. If interaction is not significant,
> you may choose to exclude interaction from the model. If there is
> interaction, then you will want to examine each farm to determine if
> cross-over interaction present.
>
> If your experiment is to determine the “best” fungicide spraying system,
> and cross-over interaction is present, then you have no “best” system. You
> might have cross-over arising because, say, system 1 ranks “best” on farm
> 1, but system 2 ranks “best” on farm 2.
>
> There is extensive literature on the topic, mostly from the plant breeding
> genotype-by-environment interaction side. Some of the associated statistics
> implemented in the agricolae package, i.e. AMMI.
>
> Peter
>
> > On Apr 24, 2017, at 6:56 AM, Juan Pablo Edwards Molina <
> edwardsmolina op gmail.com> wrote:
> >
> > I´m sorry... I´m new in the list, and when I figured out that the
> question
> > would suit best in the mixed model list I had already post it in general
> > R-help. I don´t know if there´s a way to "cancel a question"... I will
> take
> > care of it from now on.
> >
> > Dear Thierry, thanks for your answer.
> > Yes, I am not interested in the effect of a specific farm, they simply
> > represent the total of farms from the region where I want to suggest the
> > best treatments.
> >
> > I Followed your suggestions, but still have a couple of doubts,
> >
> > 1- May "farm" be include as a simple fixed effect or interacting with the
> > treatment?
> >
> > m3 = glmer(resp ~ trt * farm + (1|tree_id), family = binomial, data=df)
> > m4 = glmer(resp ~ trt + farm + (1|tree_id), family = binomial, data=df)
> >
> > 2 -
> > In case of significant
> > [ trt * farm ], should I report the results for each farm?
> >
> > Thanks again Thierry,
> >
> > Juan Edwards
> >
> >
> > *Juan*
> >
> > On Mon, Apr 24, 2017 at 4:29 AM, Thierry Onkelinx <
> thierry.onkelinx op inbo.be>
> > wrote:
> >
> >> Dear Juan,
> >>
> >> Use unique id's for random effects variables. So each bk should only be
> >> present in one farm. And each tree_id should be present in only one bk.
> In
> >> case each block has different treatments then each tree_id should be
> unique
> >> to one combination of bk and trt.
> >>
> >> Farm has too few levels to be a random effects. So either model is as a
> >> fixed effect or drop it. In case you drop it, the information will be
> >> picked up by bk. Note that trt + (1|farm) is less complex than trt *
> farm.
> >>
> >> Assuming that you are not interested in the effect of a specific farm,
> you
> >> could use sum, polynomial or helmert contrasts for the farms. Unlike the
> >> default treatment contrast, these type of contrasts sum to zero. Thus
> the
> >> effect of trt will be that for the average farm instead of the reference
> >> farm.
> >>
> >> Best regards,
> >>
> >> ir. Thierry Onkelinx
> >> Instituut voor natuur- en bosonderzoek / Research Institute for Nature
> and
> >> Forest
> >> team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
> >> Kliniekstraat 25
> >> 1070 Anderlecht
> >> Belgium
> >>
> >> 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
> >>
> >> 2017-04-21 22:32 GMT+02:00 Juan Pablo Edwards Molina <
> >> edwardsmolina op gmail.com>:
> >>
> >>> I am analyzing data from 3 field experiments (farms=3) for a citrus
> flower
> >>> disease: response variable is binomial because the flower can only be
> >>> diseased or healthy.
> >>>
> >>> I have particular interest in comparing 5 fungicide spraying systems
> >>> (trt=5).
> >>>
> >>> Each farm had 4 blocks (bk=4) including 2 trees as subsamples (tree=2)
> in
> >>> which I assessed 100 flowers each one. This is a quick look of the
> data:
> >>>
> >>> farm trt bk tree dis tot <fctr> <fctr> <fctr>
> >>> <fctr> <int> <int>
> >>> iaras cal 1 1 0 100
> >>> iaras cal 1 2 1 100
> >>> iaras cal 2 1 1 100
> >>> iaras cal 2 2 3 100
> >>> iaras cal 3 1 0 100
> >>> iaras cal 3 2 5 100...
> >>>
> >>> The model I considered was:
> >>>
> >>> resp <- with(df, cbind(dis, tot-dis))
> >>>
> >>> m1 = glmer(resp ~ trt + (1|farm/bk) , family = binomial, data=df)
> >>>
> >>> I tested the overdispersion with the overdisp_fun() from GLMM page
> >>> <http://glmm.wikidot.com/faq>
> >>>
> >>> chisq ratio p logp
> >>> 4.191645e+02 3.742540e+00 4.804126e-37 -8.362617e+01
> >>>
> >>> As ratio (residual dev/residual df) > 1, and the p-value < 0.05, I
> >>> considered to add the observation level random effect (link
> >>> <http://r.789695.n4.nabble.com/Question-on-
> overdispersion-td3049898.html
> >>>> )
> >>> to deal with the overdispersion.
> >>>
> >>> farm trt bk tree dis tot tree_id <fctr> <fctr>
> >>> <fctr> <fctr> <int> <int> <fctr>
> >>> iaras cal 1 1 0 100 1
> >>> iaras cal 1 2 1 100 2
> >>> iaras cal 2 1 1 100 3...
> >>>
> >>> so now was added a random effect for each row (tree_id) to the model,
> but
> >>> I
> >>> am not sure of how to include it. This is my approach:
> >>>
> >>> m2 = glmer(resp ~ trt + (1|farm/bk) + (1|tree_id), family = binomial,
> >>> data=df)
> >>>
> >>> I also wonder if farm should be a fixed effect, since it has only 3
> >>> levels...
> >>>
> >>> m3 = glmer(resp ~ trt * farm + (1|farm:bk) + (1|tree_id), family =
> >>> binomial, data=df)
> >>>
> >>> I really appreciate your suggestions about my model specifications...
> >>>
> >>>
> >>>
> >>> *Juan Edwards- - - - - - - - - - - - - - - - - - - - - - - -# PhD
> student
> >>> - ESALQ-USP/Brazil*
> >>>
> >>> [[alternative HTML version deleted]]
> >>>
> >>> _______________________________________________
> >>> R-sig-mixed-models op r-project.org mailing list
> >>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >>
> >>
> >>
> >
> > [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-mixed-models op r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
> _______________________________________________
> R-sig-mixed-models op r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list