How to deal with repeated measures in glmer with binary response?

Question

I would like to analyse the prevalence of a pesticide in nectar collected from fields in different years. The fields (ID) were sampled multiple times, therefore I used a GLMM with field ID as a random factor. I was surprised to see that the estimated prevalence was sometimes clearly lower than the prevalence calculated from (a) positive samples vs all samples (est_raw_sample) (b) positive fields vs all fields (est_raw_field).

Why is that? Did I make a mistake and what does the estimate of the glmer mean?

mydata = structure(list(ID = c("589", "10454", "9769", "10169", "13319", 
                               "10986", "8325", "4437", "768", "13015", "5922", "2443", "12901", 
                               "8325", "5585", "4403", "9801", "7391", "6855", "11679", "9329", 
                               "9643", "1257", "6622", "5596", "795", "1565", "12774", "7069", 
                               "4578", "7687", "1320", "4783", "6457", "11471", "6998", "11254", 
                               "10568", "5752", "7713", "7069", "7502", "2700", "10634", "8731", 
                               "12901", "5356", "6998", "6201", "6756", "1504", "9874", "16319", 
                               "2994", "16414", "4722", "2443", "7765", "12860", "289", "1242", 
                               "4722", "11535", "5910", "8325", "10536", "7168", "1497", "10435", 
                               "8076", "795", "6084", "5585", "497", "16414", "8423", "7765", 
                               "10568", "1565", "977", "4770", "6084", "8718", "5248", "8143", 
                               "13253", "4168", "2677", "3130", "14174", "6998", "6104", "768", 
                               "9477", "11934", "15370", "5844", "1320", "7873", "8423", "795", 
                               "10169", "1320", "4293", "4147", "16319", "796", "5889", "5372", 
                               "6683", "5059", "9095", "10657", "13253", "13251", "14132", "6934", 
                               "1895", "9643", "9477", "9641", "8738", "5667", "785", "12919", 
                               "9145", "6860", "9641", "6998", "14122", "4437", "10163", "6388", 
                               "4502", "8974", "12917", "4437", "1095", "920", "9848", "9630", 
                               "5519", "12924", "919", "656", "9433", "16319", "9482", "5498", 
                               "5566"), year = structure(c(4L, 3L, 2L, 3L, 3L, 1L, 3L, 3L, 3L, 
                                                           1L, 2L, 4L, 3L, 3L, 3L, 1L, 2L, 1L, 1L, 2L, 3L, 4L, 1L, 1L, 2L, 
                                                           3L, 3L, 4L, 4L, 2L, 4L, 3L, 2L, 3L, 1L, 3L, 4L, 3L, 1L, 1L, 4L, 
                                                           4L, 2L, 4L, 1L, 3L, 1L, 3L, 1L, 1L, 2L, 2L, 3L, 4L, 3L, 3L, 4L, 
                                                           3L, 4L, 4L, 2L, 3L, 3L, 1L, 3L, 4L, 2L, 3L, 2L, 2L, 3L, 3L, 3L, 
                                                           2L, 3L, 4L, 3L, 3L, 3L, 3L, 3L, 3L, 1L, 2L, 1L, 3L, 3L, 2L, 2L, 
                                                           3L, 3L, 3L, 3L, 3L, 2L, 3L, 1L, 3L, 2L, 4L, 3L, 3L, 3L, 2L, 2L, 
                                                           3L, 1L, 2L, 3L, 2L, 2L, 4L, 3L, 3L, 4L, 3L, 3L, 2L, 4L, 3L, 4L, 
                                                           2L, 1L, 1L, 3L, 1L, 2L, 4L, 3L, 3L, 3L, 3L, 1L, 2L, 1L, 4L, 3L, 
                                                           3L, 2L, 3L, 1L, 4L, 4L, 2L, 2L, 2L, 3L, 2L, 1L, 2L), .Label = c("2014", 
                                                                                                                           "2015", "2016", "2017"), class = "factor"), pesticide_found = c(1L, 
                                                                                                                                                                                           0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 
                                                                                                                                                                                           0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 
                                                                                                                                                                                           0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 
                                                                                                                                                                                           0L, 0L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 
                                                                                                                                                                                           0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 1L, 
                                                                                                                                                                                           0L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 
                                                                                                                                                                                           1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 
                                                                                                                                                                                           0L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 
                                                                                                                                                                                           1L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 
                                                                                                                                                                                           0L, 1L, 0L, 1L, 0L)), .Names = c("ID", "year", "pesticide_found"
                                                                                                                                                                                           ), row.names = c(391L, 315L, 151L, 310L, 347L, 74L, 290L, 227L, 
                                                                                                                                                                                                            174L, 81L, 124L, 415L, 338L, 288L, 240L, 15L, 152L, 50L, 46L, 
                                                                                                                                                                                                            162L, 297L, 394L, 5L, 42L, 119L, 179L, 206L, 432L, 419L, 107L, 
                                                                                                                                                                                                            430L, 196L, 109L, 256L, 77L, 268L, 425L, 319L, 25L, 52L, 474L, 
                                                                                                                                                                                                            435L, 98L, 485L, 59L, 337L, 17L, 265L, 31L, 45L, 93L, 153L, 379L, 
                                                                                                                                                                                                            426L, 382L, 230L, 402L, 279L, 433L, 408L, 91L, 231L, 331L, 29L, 
                                                                                                                                                                                                            289L, 383L, 133L, 203L, 156L, 141L, 180L, 246L, 242L, 85L, 380L, 
                                                                                                                                                                                                            392L, 278L, 317L, 211L, 188L, 235L, 245L, 58L, 116L, 55L, 344L, 
                                                                                                                                                                                                            224L, 97L, 100L, 361L, 269L, 251L, 175L, 301L, 165L, 365L, 27L, 
                                                                                                                                                                                                            200L, 139L, 436L, 181L, 309L, 199L, 104L, 103L, 378L, 4L, 122L, 
                                                                                                                                                                                                            237L, 129L, 114L, 403L, 322L, 345L, 441L, 357L, 264L, 95L, 428L, 
                                                                                                                                                                                                            300L, 431L, 146L, 22L, 3L, 340L, 61L, 132L, 461L, 266L, 353L, 
                                                                                                                                                                                                            226L, 306L, 38L, 106L, 60L, 437L, 229L, 191L, 90L, 303L, 65L, 
                                                                                                                                                                                                            396L, 390L, 89L, 86L, 148L, 377L, 150L, 20L, 118L), class = "data.frame")

Code:

# Check which fields were positive
library(dplyr)
mydata_field = mydata %>% group_by(ID, year) %>% summarise(
    pesticide_found_in_field = sum(pesticide_found)) %>% transform(
        pesticide_found_in_field = ifelse(pesticide_found_in_field > 0, 1, 0)
    )
summary(mydata_field)

mydata$pesticide_found = factor(mydata$pesticide_found)
mydata_field$pesticide_found_in_field = factor(mydata_field$pesticide_found_in_field)

pesticide_tab = with(mydata, table(year, pesticide_found)) 
pesticide_tab_field = with(mydata_field, table(year, pesticide_found_in_field)) #

library(lme4)
library(lmerTest)
pesticide_glmer = glmer(pesticide_found ~ year + (1|ID), family = "binomial", data = mydata)
summary(pesticide_glmer) # Strange that it lists each year separately

# Function to convert from odds ratio to probability of success/contamination
unlogit = function(y){
    exp(y)/(1+ exp(y))
}

pesticide_ctab = data.frame(year = 2014:2017, est = unlogit(fixef(pesticide_glmer)))

# Add estimates based on the raw data per sample and per field
positive_sample = pesticide_tab[,2]
n_sample = pesticide_tab[,1]+pesticide_tab[,2]
pesticide_ctab$est_raw_sample = positive_sample/n_sample

positive_field = pesticide_tab_field[,2]
n_field = pesticide_tab_field[,1]+pesticide_tab_field[,2]
pesticide_ctab$est_raw_field = positive_field/n_field

Answer 1

The problem is that you're overlooking the contrasts that are imposed by R on the factor year when it fits the model. By default, R uses so-called "treatment" contrasts: the coefficients from the GLMM are the differences in expected value (on a logit scale) between each level of the factor and an arbitrarily selected reference level. Thus, the coefficient (Intercept) is the expected value of year 2014 (on a logit scale), the coefficient year2015 is the difference in expected value between year 2015 and year 2014, and so on. This is one of many ways to parameterize a categorical covariate in a linear model.

Just to be a tad more explicit: let $\mathbb{E}[y_{\mathrm{year}}]$ be the expectation of the binomial probability in a given year (ie. probability of a "success" event in that year, which I guess in your case is the presence of some nasty neonic in floral nectar). The default contrast parameterization looks like: $$\text{logit }\mathbb{E}[y_{2014}] = \alpha$$ $$\text{logit }\mathbb{E}[y_{2015}] = \alpha + \beta_{2015}$$ $$\text{logit }\mathbb{E}[y_{2016}] = \alpha + \beta_{2016}$$ and so on. The coefficient $\alpha$ is what's called the (Intercept) in your output, the coefficient $\beta_{2015}$ is what's called year2015 ... and so on. A more concise way to write this for a the $i$th datapoint is: $$\text{logit }\mathbb{E}[y_i] = \alpha + d_i^{2015} \beta_{2015} + d_i^{2016} \beta_{2016} + ...$$ Where the $d_i^{\mathrm{year}}$ are a set of dummy variables associated with the $i$th observation, that will equal 1 if the observation is from the indicated particular year and 0 otherwise.

Because of this contrast parameterization, it's completely meaningless for you to back-transform these coefficients to a $[0,1]$ scale. The appropriate approach, which will work regardless of the chosen type of contrast, is to first calculate the predicted values for each year (on the logit scale), and then back-transform. Also I'll note that the back-transform is just the CDF of a logistic distribution, no need to define a new function.

pred_vals       <- data.frame(year=factor(2014:2017))
pred_vals$glmer <- plogis(predict(pesticide_glmer, pred_vals, re.form=NA))
pred_vals$raw   <- pesticide_ctab$est_raw_sample

Note that the model predictions, while close to those calculated from the raw data, do differ a bit. This is because of the random effect, which accounts for dependence between observations with the same level of ID (effectively discounting dependent observations when estimating the fixed effects). Were you to fit a GLM instead, the model predictions would be identical to those from the raw data.

pesticide_glm   <- glm(pesticide_found~year,family="binomial",data=mydata)
pred_vals$glm <- plogis(predict(pesticide_glm, pred_vals))
pred_vals
#  year      glmer       raw       glm
#1 2014 0.78909629 0.7692308 0.7692308
#2 2015 0.09918552 0.1142857 0.1142857
#3 2016 0.66740090 0.6562500 0.6562500
#4 2017 0.13853975 0.1600000 0.1600000

Why is the "treatment" contrast useful enough to be a default? Well, in many cases, we are actually interested in the difference between various manipulated treatment levels and a control group. Thus the estimated coefficients are actually quantities we care about, and the standard errors give a useful summary of our uncertainty w.r.t these estimates. In your case, you might find it more informative to not use contrasts, and instead have the coefficient for each year represent the logit expected value for that year. In R, this can accomplished by adding a 0 or -1 to the model formula:

pesticide_glmer_nocon = glmer(pesticide_found ~ 0 + year + (1|ID), family = "binomial", data = mydata)
plogis(fixef(pesticide_glmer_nocon))
#  year2014  year2015  year2016  year2017 
# 0.7890976 0.0991851 0.6674002 0.1385394