Chapter 15 Generalized Linear Mixed Effects Models

15.1 Exercise 1

Consider the RIKZdat data set from Chapter 2 of the book, which can be accessed using:

library(Data4Ecologists)
data(RIKZdat)

To help with model convergence, we will use a scaled and centered version of NAP. We will also create a new variable exposure.low that combines the lowest two exposure categories using:

library(dplyr)
RIKZdat <- RIKZdat %>% mutate(NAPc = as.numeric(scale(NAP)),
                              exposure.low = case_when(exposure ==  11 ~ 1, exposure < 11 ~ 0))
  1. Alter the code, below, used to fit a linear mixed effects model to the RIKZdat data set from Chapter 2 of the book so that it fits a model where, conditional on the random intercepts (and fixed effect covariates), the data are assumed to be Poisson distributed.
jags.lme.rc<-function(){
  
# Priors for the intercepts
 for (i in 1:n.groups){    
  # To allow for correlation between alpha[i] and beta[i], we need to model their joint (multivariate) distribution 
    alpha[i] <- B[i,1]  # Random intercepts
    beta1[i] <-B[i,2] # Rancom slopes for NAP
    B[i,1:2]~ dmnorm(B.hat[i,], Tau.B[,]) # distribution of the vector (alpha[i], beta[i])
    B.hat[i,1]<-mu.int # mean for intercepts
    B.hat[i,2]<-mu.slope # mean for slopes
 }

 # Hyperpriors for intercepts and slopes
  mu.int ~ dnorm(0, 0.001)       
  mu.slope ~dnorm(0,0.001)
 
 # Hyperpriors for Sigma =  var/cov matrix of the slope/intercept parameters
  sigma.int ~ dunif(0,100) #sd intercepts
  sigma.slope~dunif(0,100) # sd of slopes
  rho~dunif(-1,1) # correlation among intercepts and slopes
  Sigma.B[1,1]<-pow(sigma.int,2)
  Sigma.B[2,2]<-pow(sigma.slope,2)
  Sigma.B[1,2]<-rho*sigma.int*sigma.slope
  Sigma.B[2,1]<-Sigma.B[1,2]
 
 # Tau = inverse of Sigma (analogous to precision for univariate normal distribution)
  Tau.B[1:2,1:2]<-inverse(Sigma.B[,])
 
 # Priors for exposure regression parameter and within individual variance
  beta2 ~ dnorm(0, 0.001)       # Common slope exposure
  tau <- 1 / ( sigma * sigma)       # Residual precision
  sigma ~ dunif(0, 100)         # Residual standard deviation

# Likelihood
 for (i in 1:nobs) {
    mu[i] <- alpha[Beach[i]] + beta1[Beach[i]]*NAPc[i]+beta2*exposure.low[i] # Expectation
    Richness[i] ~ dnorm(mu[i], tau)     # The random variable
 }
}

  1. Describe the model using a set of equations.

  2. Fit the same model using the lmer function in the lme4 package and compare the estimates of the fixed effect and variance parameters.

15.2 Exercise 2

Consider the RIKZdat data set from Chapter 2 of the book, which can be accessed using:

library(Data4Ecologists)
data(RIKZdat)

To help with model convergence, we will use a scaled and centered version of NAP. We will also create a new variable exposure.low that combines the lowest two exposure categories using:

library(dplyr)
RIKZdat <- RIKZdat %>% mutate(NAPc = as.numeric(scale(NAP)),
                              exposure.low = case_when(exposure ==  11 ~ 1, exposure < 11 ~ 0))

  1. Fit a model relating species richness (Richness) to NAP to exposure.low. Include a random intercept for each beach. Assume, conditional on the random intercepts (and fixed effect covariates), that the data are Poisson distributed.

  2. Update the model from step 1 so that it also allows the effect of NAP to vary by beach using a random coefficient model. Why can we not allow the effect of exposure.low to also vary by beach?

  3. Use glmmTMB function in the glmmTMB package with family = nbinom2 to fit the model in step 2 but assuming that the distribution of species richness, conditional on the random effects (and fixed effect covariates) follows a negative binomial distribution.

  4. Compare the 3 models in steps 1-3 using AIC. Evaluate the reliability of the assumptions for the model with the lowest AIC using the check_model function in the performance package. Make sure to comment on the plots - particularly if you see potential assumption violations.

  5. Write down a set of equations to describe the chosen model in step 4.

15.3 Exercise 3

Consider again the data set from exercise 12.4 containing observations of bird species richness in habitat patches classified according to different landscape types (landscape = Agriculture, Bauxite, Forest, or Urban). The data can be accessed using:

library(Data4Ecologists)
data(birds, package = "Data4Ecologists")

It may also help to drop observations with missing species richness data (e.g., when specifying the model using JAGS):

library(dplyr)
birds <- birds %>% filter(is.na(S)!=TRUE)

In Exercise 12.4, we considered the following model:

\[\begin{gather} S_i \sim Poisson(\lambda_i) \\ \log(\lambda_i) = \beta_0 + \beta_1 I(landscape = Bauxite)_i + \beta_2 I(landscape = Forest)_i + \\ \beta_3 I(landscape = Urban)_i +\beta_4 \log(Area)_i + \beta_5 \log(Area)^2_i \end{gather}\]

However, we noted that each habitat patch was monitored in 3 separate years. In this exercise, we will attempt to address the fact that each patch was observed repeatedly.

  1. Modify the model to include a random intercept for each habitat patch and fit the model using the glmer function in the lme4 package. Compare this model to the Poisson model without a random intercept. Comment on your findings.

  2. Fit the model using JAGS. Note, if you want to compare coefficients in the two models, make sure to use the same parameterization (i.e., either effects or means coding). Also, you will want to make sure to use raw (or orthogonal) polynomials in both models!

  3. Describe the model in step 2 using a set of equations. Make sure to also describe the assumed prior distributions.

15.4 Exercise 4

For this exercise, you will explore differences between subject-specific (conditional) and population-averaged (marginal) means estimated using generalized linear mixed effects models and generalized estimating equations.

Thompson et al. (2016)9 counted birds in a set of waterfowl production areas in Minnesota to evaluate the effects of tree removal treatments. For this exercise, we will ignore the treatments and just explore how total counts of woodland birds changed over the course of the study.

The data for this exercise can be accessed using:

library(Data4Ecologists)
data(d.woodall)
  1. Determine the mean number of birds counted in each year across the different waterfowl production areas. This can be accomplished using the code below (you will want to save this output to an R object so you can plot the counts over time in step 6):
d.woodall %>% group_by(year) %>%
  summarize(meancount = mean(total))
## # A tibble: 7 × 2
##    year meancount
##   <dbl>     <dbl>
## 1     0     1.26 
## 2     1     0.896
## 3     2     0.818
## 4     3     0.775
## 5     4     0.686
## 6     5     0.794
## 7     6     0.761
  1. Fit the following generalized linear mixed effect model to the counts of woodland birds:

glmer(total ~ as.factor(year) + (1 | wpa), family = poisson(), data = d.woodall)

Estimate the mean count for a “typical waterfowl production area” in each year (e.g., using the predict function).

  1. Describe the model from step 2 using a set of equations.

  2. Fit two generalized estimating equation models to these same data, first using an independence working correlation structure and then an exchangeable working correlation structure.

You can use the geeglm function in the geepack package:

geeglm(total ~ as.factor(year), id = wpa, data = d.woodall, family = poisson(), corstr = "independence") geeglm(total ~ as.factor(year), id = wpa, data = d.woodall, family = poisson(), corstr = "exchangeable")

Estimate the population mean count (i.e., across the population of waterfowl production areas) in each year (e.g., using the predict function) for both of these models.

  1. Use the model from step 2 and the results from Section 19.7.2 of the book to estimate the (marginal) population mean count in each year.

  2. Compare the observed means from step 1 to the estimated means in steps 2, 4 and 5. Provide an explanation for any observed differences.

15.5 Exercise 5

For this exercise, we will consider data from the following paper (see abstract below):

Dakin and Montgomeric (2014). Deceptive copulation calls attract female visitors to Peacock leks. The American Naturalist 558-564.

Abstract: Theory holds that dishonest signaling can be stable if it is rare. We report here that some peacocks perform specialized copulation calls (hoots) when females are not present and the peacocks are clearly not attempting to copulate. Because these solo hoots are almost always given out of view of females, they may be dishonest signals of male mating attempts. These dishonest calls are surprisingly common, making up about a third of all hoot calls in our study populations. Females are more likely to visit males after they give a solo hoot call, and we confirm using a playback experiment that females are attracted to the sound of the hoot. Our findings suggest that both sexes use the hoot call tactically: females to locate potential mates and males to attract female visitors. We suggest that the solo hoot may be a deceptive signal that is acquired and maintained through reward-based learning.

The data can be accessed using:

library(Data4Ecologists)
data(peacock)

We will use a centered and scaled version of Julian date to make it easier to fit the model and also change some of the character variables to factors. It will also help with interpretation if we change the reference level for the PrePost variable to be Pre using the following code:

# Character to factor variables
peacock$Place <- as.factor(peacock$Place) 
peacock$PrePost <- as.factor(peacock$PrePost )
peacock$DofY <- scale(peacock$DayofYear) # centers and scales this variable
peacock$PrePost <- relevel(peacock$PrePost, ref = "Pre") # ensures the reference level is "pre hoot"
  1. Use ggplot to explore how the probability of a female being located within 5m of a focal male depends on Julian date, site (i.e., place), and if the observation took place before or after the male produced a hoot call.
library(ggplot2)
library(ggthemes) #colorblind pallete
ggplot(peacock, aes(x = DayofYear, y = FemaleVisit.5m, colour = PrePost)) +
     geom_point(position = position_jitter(w = 2, h = 0.05), size=3) +
     geom_smooth(method="glm",  method.args = list(family = "binomial"))  + 
     facet_wrap(~Place, scales="free") +  scale_colour_colorblind()
## `geom_smooth()` using formula = 'y ~ x'

  1. Use the glmer function in the lme4 library to fit a generalized linear mixed effects model in which the log odds of a male attracting a female is modeled as a linear function of Julian date (DofY), site (i.e., Place), and whether the observation period occurred before or after a hoot call (PrePost). Include a random intercept for each male (MaleID).

  2. Interpret the fixed effects parameters and the variance parameter reported by the summary function.

  3. Use the confint function to get profile-likelhood based and Wald-based (i.e., \(\pm 1.96SE\)) confidence intervals. See ?confint.merMod for the correct syntax. (Note: bootstrap and profile-likelihood based intervals should be “better” than Wald intervals, but there are tradeoffs in terms of time required to get an answer. I tried to include bootstrap confidence intervals when producing the answer key, but it took so long I eventually stopped the computation!).

  4. Fit the same model using JAGS.[Optional: assess goodness-of-fit using Pearson residuals.]. The MCMC samples will be highly autocorrelated so it is important to run the model for longer than most of the examples we have looked at to date. Run the sampler for 80,000 iterations, with a burn-in of 10,000 iterations. Again, use 3 chains but this time with thin = 3 (keep every 3rd observation). Assess convergence using traceplots, density plots, and Gelman-Rubin statistic (R).

  5. Use the geeglm function in the geepack library to fit a generalized estimating equations model in which logit(p) is modeled as a linear function of date, site, and pre- versus post-hoot. Try both independence and exchangeable working correlation structures. In both cases, add an argument scale.fix = TRUE so that R does not estimate an overdispersion parameter.

  6. Describe/interpret the parameters in the GEE model (fit with geepack). Note: interpretation is the same for both working correlation structures, so you only need to do this once!

  7. Compare 95% confidence intervals from both GEE models (with working independence and ex- changeable correlation structures) to 95% credible intervals for the glmm (fit using JAGS) and 95% confidence intervals (profile-likelhood and Wald-based) you calculated using confint and glmer. Comment on your findings. HINT: you can access the standard errors for the coefficients in the gee model using coef(summary([name of gee model]))$Std.err. You can then get Wald-based intervals by taking \(\pm 1.96SE\).

15.6 Exercise 6

For this exercise, you will consider a data set from Howard et al. (2022).

Howard, S. R., Greentree, J., Avarguès-Weber, A., Garcia, J. E., Greentree, A. D., & Dyer, A. G. (2022). Numerosity Categorization by Parity in an Insect and Simple Neural Network. Frontiers in Ecology and Evolution, 252.

The authors of this study evaluated whether: 1) free-flying honeybees could be trained to differentiate between odd and even numbers of geometric shapes, and 2) whether they could transfer this learning to new situations. Bees were initially trained to prefer odd or even numbers of shapes by associating stimuli (cards with 1-10 shapes) with either sugar or quinine. After training, each bee participated in 20 learning-test trials in which they made a choice between stimuli with an even or odd number of shapes. These trials again used stimuli with between 1 and 10 shapes. Lastly, each bee also participated in 20 transfer-test trials where they chose between stimuli consisting of either 11 or 12 shapes. The full experimental design is illustrated in the figure below and more details can be found in the paper if you are interested.

Figure 1 from Howard et al. (2022). Sequence of the different experimental phases and examples of stimuli which could be presented to a bee during the preference test, learning phase, learning test, and transfer test.

Figure 15.1: Figure 1 from Howard et al. (2022). Sequence of the different experimental phases and examples of stimuli which could be presented to a bee during the preference test, learning phase, learning test, and transfer test.

Data from the experiments can be accessed using:

library(Data4Ecologists)
data(beedat)

Key variables in the data set include:

  • Bee: a unique identifier for each bee
  • Choice: equal to 1 if the bee made a correct choice and 0 otherwise. For example, a bee trained to prefer even numbers of shapes would have a correct choice if it choose a stimuli with 2, 4, …, 12 shapes.
  • Test: identifies the type of trial based on whether bees were trained to select even or odd numbers and whether the trial was a learning or transfer test. Test = 1 for learning tests involving even trained bees. Test = 2 for learning tests involving odd trained bees. Test = 3 for transfer tests involving even trained bees. Latly Test = 4 for transfer tests involving odd trained bees.
  • Ieven: an indicator variable equal to 1 for even trained bees and 0 for odd trained bees
  • Itransfer: an indicator variable equal to 1 if the observation was from a transfer test and 0 for a learning test

Note: the Test variable should be treated as a factor using:

beedat$Test <- as.factor(beedat$Test)

  1. Use a generalized linear mixed effects model to evaluate: a) if it is easier for bees to learn to find even patterns relative to their ability to discern odd patterns; b) if bees performed better during the learning phase relative to the transfer phase; and c) if bees do better than random at picking the right stimuli. Provide statistical evidence to support your conclusions (e.g., confidence intervals or p-values for appropriate hypothesis tests or comparison of models using AIC). (11 pts)

  2. Using a set of equations, describe the model(s) you fit in part 1. Identify the estimated parameters in R output when using the summary function. (4 pts)

  3. Use your fitted model from 4 to estimate the probability that a “typical bee” trained to select stimuli with an even number of shapes chooses correctly on a transfer test. (3 pts)

  4. Explain why the question in 5 refers to a “typical bee”. (2 pts)

If you are curious, I’ve included a video, below, showing an example transfer trial!

  1. Thompson, S. J., Arnold, T. W., Fieberg, J., Granfors, D. A., Vacek, S., and Palaia, N. (2016). Grassland birds demonstrate delayed response to large‐scale tree removal in central North America. Journal of Applied Ecology, 53(1), 284-294.↩︎