How to summarise binomial GLM output?

Question

I have a binomial GLM:

all.fit <-  glm(Presence/Total~Season*ToD*Site, family = binomial, weights = Total, data = all.dt)
Call:  glm(formula = Presence/Total ~ Season * ToD * Site, family = binomial, 
    data = all.dt, weights = Total)
Coefficients:
                         (Intercept)                          SeasonSpring

                            -4.66074                               0.60478

                        SeasonSummer                          SeasonWinter

                            -0.82682                               1.26868

                              ToDDay                               ToDDusk

                             0.09865                               2.48420

                            ToDNight                             SiteKawau

                             2.65294                               3.09495

                          SiteNoises                        SiteTawharanui

                             3.48048                               2.94694

                        SiteTiritiri                   SeasonSpring:ToDDay

                             3.25976                               0.90540

                 SeasonSummer:ToDDay                   SeasonWinter:ToDDay

                             2.16685                              -0.12465

                SeasonSpring:ToDDusk                  SeasonSummer:ToDDusk

                            -2.78750                              -1.52114

                SeasonWinter:ToDDusk                 SeasonSpring:ToDNight

                            -0.90286                              -1.66162

               SeasonSummer:ToDNight                 SeasonWinter:ToDNight

                            -1.85826                              -0.21758

              SeasonSpring:SiteKawau                SeasonSummer:SiteKawau

                            -0.67217                               1.24760

              SeasonWinter:SiteKawau               SeasonSpring:SiteNoises

                            -1.68237                              -0.79023

             SeasonSummer:SiteNoises               SeasonWinter:SiteNoises

                             0.42570                              -0.73055

         SeasonSpring:SiteTawharanui           SeasonSummer:SiteTawharanui

                            -0.25773                               1.32456

         SeasonWinter:SiteTawharanui             SeasonSpring:SiteTiritiri

                            -0.51055                              -0.51381

           SeasonSummer:SiteTiritiri             SeasonWinter:SiteTiritiri

                             0.95804                              -1.29898

                    ToDDay:SiteKawau                     ToDDusk:SiteKawau

                             0.70933                              -2.35320

                  ToDNight:SiteKawau                     ToDDay:SiteNoises

                            -5.33351                               0.27852

                  ToDDusk:SiteNoises                   ToDNight:SiteNoises

                            -2.53814                              -4.04432

               ToDDay:SiteTawharanui                ToDDusk:SiteTawharanui

                            -0.16455                              -2.76975

             ToDNight:SiteTawharanui                   ToDDay:SiteTiritiri

                            -2.62701                               0.81969

                ToDDusk:SiteTiritiri                 ToDNight:SiteTiritiri

                            -2.24361                              -4.14062

       SeasonSpring:ToDDay:SiteKawau         SeasonSummer:ToDDay:SiteKawau

                            -0.19831                              -1.61316

       SeasonWinter:ToDDay:SiteKawau        SeasonSpring:ToDDusk:SiteKawau

                             0.57193                               2.48229

      SeasonSummer:ToDDusk:SiteKawau        SeasonWinter:ToDDusk:SiteKawau

                             1.37786                               0.81702

     SeasonSpring:ToDNight:SiteKawau       SeasonSummer:ToDNight:SiteKawau

                             1.41405                               2.06227

     SeasonWinter:ToDNight:SiteKawau        SeasonSpring:ToDDay:SiteNoises

                             0.08689                              -0.28007

      SeasonSummer:ToDDay:SiteNoises        SeasonWinter:ToDDay:SiteNoises

                            -1.48112                              -0.06044

     SeasonSpring:ToDDusk:SiteNoises       SeasonSummer:ToDDusk:SiteNoises

                             1.93847                               0.63232

     SeasonWinter:ToDDusk:SiteNoises      SeasonSpring:ToDNight:SiteNoises

                             0.83278                               0.69621

    SeasonSummer:ToDNight:SiteNoises      SeasonWinter:ToDNight:SiteNoises

                             0.06911                               0.48690

  SeasonSpring:ToDDay:SiteTawharanui    SeasonSummer:ToDDay:SiteTawharanui

                            -0.23150                              -1.96624

  SeasonWinter:ToDDay:SiteTawharanui   SeasonSpring:ToDDusk:SiteTawharanui

                             0.43487                               2.04614

 SeasonSummer:ToDDusk:SiteTawharanui   SeasonWinter:ToDDusk:SiteTawharanui

                             0.59597                               0.73280

SeasonSpring:ToDNight:SiteTawharanui  SeasonSummer:ToDNight:SiteTawharanui

                             1.02467                               1.19957

SeasonWinter:ToDNight:SiteTawharanui      SeasonSpring:ToDDay:SiteTiritiri

                            -0.16571                              -0.48960

    SeasonSummer:ToDDay:SiteTiritiri      SeasonWinter:ToDDay:SiteTiritiri

                            -1.63190                               0.35354

   SeasonSpring:ToDDusk:SiteTiritiri     SeasonSummer:ToDDusk:SiteTiritiri

                             1.77115                               0.79347

   SeasonWinter:ToDDusk:SiteTiritiri    SeasonSpring:ToDNight:SiteTiritiri

                             0.48411                               1.51686

  SeasonSummer:ToDNight:SiteTiritiri    SeasonWinter:ToDNight:SiteTiritiri

                             1.27621                               0.77385
Degrees of Freedom: 22998 Total (i.e. Null);  22919 Residual
Null Deviance:      105700 
Residual Deviance: 63170    AIC: 87340

And then I do an anova:

> anova(all.fit,test = "Chisq")
Analysis of Deviance Table
Model: binomial, link: logit
Response: Presence/Total
Terms added sequentially (first to last)
            Df Deviance Resid. Df Resid. Dev  Pr(&gt;Chi)    

NULL                            22998     105714

Season           3   1005.8     22995     104709 < 2.2e-16 ***
ToD              3  16404.9     22992      88304 < 2.2e-16 ***
Site             4   7664.5     22988      80639 < 2.2e-16 ***
Season:ToD       9   2793.5     22979      77846 < 2.2e-16 ***
Season:Site     12   3890.0     22967      73956 < 2.2e-16 ***
ToD:Site        12   9804.5     22955      64151 < 2.2e-16 ***
Season:ToD:Site 36    981.0     22919      63170 < 2.2e-16 ***

Signif. codes:  0 ‘*’ 0.001 ‘’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

In this example, all of the interactioms are significant. But, I don't want to report a huge list of significant combinations of factors in my report. Since I guess if the interaction is important, the single effects don't matter? How could I state the significance of my findings more simply?

Could I use the output of anova(all.fit)? If I did so, what would I then be reporting?

Your code doesn't make sense. lmer() doesn't take a family= argument. Did you mean to use glmer()? However, you don't indicate any random effects in the formula (eg, something like +(1:id)). Do you have such? Otherwise, it would be better to use glm(). In fact, for the purposes of this question, glm() would be much preferred, since it's simpler & nothing in your Q relates to the greater complexities associated w/ the other functions. — gung - Reinstate Monica, Jun 17 '21 at 14:42
@gung-ReinstateMonica Sorry, I have corrected my code above! I see that I can use
anova(all.fit,test="LRT")

but I am not sure what this is telling me. Or how I should report it. I am guessing I should include the df and the residual df for each factor, the significance, but should I also report the deviance? — Louise Wilson, Jun 17 '21 at 20:50
Thanks. It might help if you post example output. Are you just wondering about which p-values to list in the results section, or how to best convey the meaning of it all to an audience? (Something else? What?) Who is the audience? What is the context? Is this for a poster at a scientific conference, a presentation (to whom?), a submission to a scientific journal, something else? More information will help. — gung - Reinstate Monica, Jun 18 '21 at 11:05
@gung-ReinstateMonica Thanks! I have updated the question with example outputs. I want to report the significance of the interaction, and yes, to convey the meaning at a conference and within a report, a presentation to other scientists doing this same research. Aswell as stating the significance of the interaction, I want to be able to look at comparisons between the different levels of each variable, and so for example say, 'summer at one site had a higher proportion of 'presence' than winter'. Does that make sense? — Louise Wilson, Jun 20 '21 at 06:19

score 1 · Accepted Answer · answered Jun 20 '21 at 14:17

You have all categorical variables, and they are all nested within the top level, three-way interaction. Make sure you understand how categorical variables are coded by default and what the standard output means (see: here and here). In R, the output from summary(model) will differ from anova(model). The latter will output a sequential (Type I, see: here) analysis of deviance table (since you have a GLM). In general, p-values will differ between these two (see: here), although in your case the summary output doesn't actually include the test of any of the variables when taken as a whole (again, read up on how variables are coded and presented). In general, drop1(model, test="LRT") will better accord with people's intuitions, but since you only have categorical variables all nested within a three-way interaction, drop1 will only show the p-value for the interaction, which will be the same as listed in the output from anova. As it happens, that's really the only p-value you need. Any other p-values wouldn't really mean what people think they mean.

After that, you have 3x3x4=36 combinations of levels. I would simply make a table with the predicted means, or represent the means with a plot (like a dotplot or a barplot). You'll want to order the cells / dots / bars in such a way to allow you to conveniently point out the comparisons of interest.

Thank you! Turns out the p-value and deviance output from anova(all.fit,test = "Chisq") is the same as LRT output of drop1(model, test="LRT"). Should I report this as (LRT 36, 250944 = 976.67, p = < 0.005) ? — Louise Wilson, Jun 22 '21 at 03:21
@LouiseWilson, in R "LRT" & "Chisq" are synonyms in drop1() & anova(). I would report it as, ($\chi^2=981.0$, $df=36$, $p<0.0001$. — gung - Reinstate Monica, Jun 22 '21 at 11:25

score 0 · Answer 2 · answered Jun 17 '21 at 09:46

If you want to hold on to your current presentation, it is necessary to include at least the third interaction term. If you cannot reduce that, it will not make sense to reduce any of the lower level interaction terms. To test the significance of time, you can remove time completely from from the model, and then test from one model to the other. Similar for season and site.

If you variables time, site, season are all discreate, you can combine them in a new variables and use anova, to show that not all groups can be considered as one. You can then test for reduction to a model without time, without site, and without season, to consider them separately. I am aware that this might not be exactly the answer you are looking for, but I think this is the simplest solution for interpretation.

I hope it answers your question.

How to summarise binomial GLM output?

2 Answers2