Generalized linear (mixed) model, binomial - help!

Question

I work in biology and I´ve done an experiment exposing an invertebrate to a pesticide at different temperatures. One of my endpoints is hatching success of their eggs. The animals lay clutches of eggs, usually 10-30 eggs. My impression is that I shouldn’t model the data using proportions as datapoints (e.g., 0.7 proportion of the eggs in one clutch hatched), but rather use raw data (1/0 for hatched/not hatched, per egg). Then I need to add “clutch-ID” as a random factor. I, therefore, plan to use a generalized linear mixed model (glmer() from the lme4 package) with binomial family link = "logit".

My model:

glmer(hatchingsuccess ~ Temperature + Pesticide + (1 | clutch.id), 
      data = eggdata, family = binomial(link = "logit"))

The sample size is approx. 8000 eggs (in about 500 clutches). I also first included an interaction term with Temperature * Pesticide, but then nothing was significant, so I eliminated that.

Summary:

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial ( logit )
Formula: hatchingsuccess ~ Temperature + Pesticide + (1 | clutch.id)
   Data: eggdata

     AIC      BIC   logLik deviance df.resid 
  6879.3   6914.3  -3434.6   6869.3     8139 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.2772 -0.2007  0.2111  0.4408  3.7051 

Random effects:
 Groups   Name        Variance Std.Dev.
 Batch.id (Intercept) 5.121    2.263   
Number of obs: 8144, groups:  Batch.id, 560

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)   
(Intercept)       0.62213    0.43920   1.416  0.15663   
Temperature       0.08385    0.02709   3.096  0.00196 **
lowPesticide -0.29740    0.25770  -1.154  0.24847   
highPesticide -0.50853    0.26549  -1.915  0.05543 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
optimizer (Nelder_Mead) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.0023972 (tol = 0.002, component 1)

The proportions of eggs hatched look like this (the different colours represent pesticide treatments, and the three plots are for three different temperatures.

These are diagnostic plots using the package DHARMa:

1. Does this R output from the model and diagnostics look ok? I got an answer about the converge problem and will look into that.

2. Is it ok to present these proportion plots and refer to model statistics, or do I need to plot single datapoints on the y-axis and a line for the model? I think these proportion plots more clearly show the data/results, but I need to have some stats to refer to.

Answer 1

Regarding your second question, Is it ok to present these proportion plots, I would say yes, this is ok, but your plot is missing error bars. As you may know, most journal editors now ask for a measure of precision in plots, such as 95% confidence intervals.

In your case, it is straightforward to have a plot with 95% confidence interval if you use analysis of proportions library of mine ANOPA.

To illustrate, I will suppose that you have the data in a data.frame() called dta with 0s and 1s and also that your factors are called pesticide (three levels, say A, B, and C) and temperature (with three levels as well, say cold, warm and hot). The beginning of the data frame could be, e.g.,

head(dta)
#    id pesticide temperature s
# 1   1         A        cold 1
# 2   2         A        cold 1
# 3   3         A        cold 1
# 4   4         A        cold 0
# 5   5         A        cold 1
# 6   6         A        cold 1

From there, you can do:

library(ANOPA)
w <- anopa( s ~ pesticide + temperature, dta)
anopaPlot(w)

from which you get

of if you prefer bars:

anopaPlot(w, plotStyle="bar"

This plot is a ggplot so it can be followed with any ggplot attributes, e.g.,

library(ggplot2)
anopaPlot(w, plotStyle="bar" ) + theme_bw()

---

BTW, if you want to know if there are effects in the simulated data that were illustrated above, simply run

summary(w)

In the present example, the sample sizes of the data that I simulated were moderately large (100 per group, for a total of 900 simulated success or failures, 0s or 1s). Here, as was hinted by the plot, the effect of pesticide is significant ($F(2, \infty) = 7.08$, $p = .0008$), as seen in the ANOPA table:

                            MS  df        F   pvalue correction    Fcorr pvalcorr
pesticide             0.017643   2 7.092461 0.000831   1.002222 7.076735 0.000845
temperature           0.005343   2 2.147827 0.116738   1.002222 2.143065 0.117295
pesticide:temperature 0.002202   4 0.885204 0.471699   1.033333 0.856649 0.489127
Error(between)        0.002488 Inf                      

You can request a simple effect plot of the factor pesticide alone with

 anopaPlot(w, ~ pesticide) + theme_bw()

which shows