Chapter 18 Final model

The same methods we just saw can be used with a GLMM or LMM to choose between models with various random intercepts and/or random slopes and to choose fixed effects to keep in final model.

Recall the two approaches :

an information theoretic approach (e.g., AICc - Workshop 5)
a frequentist approach (where the significance of each term is evaluated using a likelihood ratio test; LRT, with the anova() function)

We first code potential models and compare them using AICc:

mpl1 <- glmer(total.fruits ~ nutrient * amd + rack + status +
    (1 | X) + (1 | popu) + (1 | gen), data = dat.tf, family = "poisson",
    control = glmerControl(optimizer = "bobyqa"))

mpl2 <- update(mpl1, . ~ . - rack)  # model without rack
mpl3 <- update(mpl1, . ~ . - status)  # model without status
mpl4 <- update(mpl1, . ~ . - amd:nutrient)  # without amd:nutrient interaction

aic_tab <- MuMIn::model.sel(mpl1, mpl2, mpl3, mpl4)
(round(aic_table <- aic_tab[, c("AICc", "delta", "df")], digits = 2))

##         AICc delta df
## mpl1 5015.73  0.00 10
## mpl4 5017.11  1.38  9
## mpl3 5017.22  1.49  8
## mpl2 5070.75 55.02  9

Note: We do not cover all possible models above, however, the interaction amd:nutrient can only be evaluated if both amd and nutrient (i.e., the main effects) are included in the model.

Alternatively, we can use drop1() and dfun() functions to evaluate our fixed effects (dfun() converts the AIC values returned by the drop1() into \(\Delta\)AIC values)

dd_LRT <- drop1(mpl1, test = "Chisq")
(dd_AIC <- dfun(drop1(mpl1)))

## Single term deletions
## 
## Model:
## total.fruits ~ nutrient * amd + rack + status + (1 | X) + (1 | 
##     popu) + (1 | gen)
##              npar   dAIC
## <none>             0.000
## rack            1 55.083
## status          2  1.612
## nutrient:amd    1  1.444

Strong rack effect (dAIC = 55.08 if we remove this variable)
Effects of status and interaction term are weak (dAIC < 2)

Let’s start by removing the non-significant interaction term to test main effects of nutrient and clipping

mpl2 <- update(mpl1, . ~ . - amd:nutrient)

mpl3 <- update(mpl2, . ~ . - rack)  # pas de rack ou d'interaction
mpl4 <- update(mpl2, . ~ . - status)  # pas de status ou d'interaction
mpl5 <- update(mpl2, . ~ . - nutrient)  # pas de nutrient ou d'interaction
mpl6 <- update(mpl2, . ~ . - amd)  # pas d'herbivorie ou d'interaction

Choose the method you want to select the best model:

Method with AICc

aic_tab2 <- MuMIn::model.sel(mpl2, mpl3, mpl4, mpl5, mpl6)
(round(aic_table2 <- aic_tab2[, c("AICc", "delta", "df")], digits = 2))

##         AICc  delta df
## mpl2 5017.11   0.00  9
## mpl4 5018.28   1.17  7
## mpl6 5027.27  10.16  8
## mpl3 5071.28  54.17  8
## mpl5 5152.74 135.63  8

Method with drop1()

dd_LRT2 <- drop1(mpl2, test = "Chisq")
dd_AIC2 <- dfun(drop1(mpl2))

What are our conclusions ?

Both the main effects of nutrient and clipping are strong (large change in AIC of \(135.6\) (mpl5) and \(10.2\) (mpl6) if either nutrient or clipping are dropped, respectively).

Our final model includes :

Fixed effects * nutrients * clipping * rack

Random effects * observation-level random effect (1|X) * populations (1|popu) * genotypes (1|gen)

18.1 Challenge 9

Use the inverts dataset (larval development times (PLD) of 74 marine invertebrate and vertebrate species reared at different temperatures and time), answer the following questions:

1. What is the effect of feeding type and climate (fixed effects) on PLD?
1. Does this relationship vary among taxa (random effects)?
1. What is the best distribution family for this count data?
1. Finally, once you determined the best distribution family, re-evaluate your random and fixed effects.

Challenge 9 Solution:

# inverts <- read.csv('data/inverts.csv', header = TRUE)
# head(inverts) table(inverts$temp, inverts$feeding.type)

# mod.glm <- glm(PLD ~ temp + feeding.type, family =
# poisson(), data = inverts) summary(mod.glm)
# drop1(mod.glm, test = 'Chisq')

# boxplot(PLD ~ temp, data = inverts) boxplot(PLD ~
# feeding.type, data = inverts)

# boxplot(predict(mod.glm, type = 'response')~inverts$temp)

# plot()

# modglm <- glm(PLD ~ temp + feeding.type, family =
# poisson(), data = inverts)

# r2 <- ranef(mpl1, condVar = TRUE) d2 <- dotplot(r2,
# par.settings = pp)

# plot(aggregate(PLD ~ taxon, FUN = mean, data =
# inverts)[,2], aggregate(PLD ~ taxon, FUN = var, data =
# inverts)[,2], pch = 19) abline(a = 0, b = 1, lty = 2)

# mod.glmer <- glmer.nb(PLD ~ temp + feeding.type +
# (1|taxon), data = inverts) mod.glm <- glm.nb(PLD ~ temp +
# feeding.type, family = poisson(), data = inverts)

# plot(aggregate(PLD ~ taxon, FUN = var, data =
# inverts)[,2], aggregate(PLD ~ taxon, FUN = mean, data =
# inverts)[,2]) abline(a = 0, b = 1, lty = 2 )