?family9 GL(M)Ms
NoteIn this chapter
In this chapter, you will learn
- how a GLM connects the linear predictor to the response via a link function and a distribution, and which variant to use for which type of response (Section 9.2);
- how to check GLMs with simulated residuals and fix the typical problems: overdispersion and zero-inflation (Section 9.3);
- why effects and interactions depend on the scale (link vs. response scale), and how to report ratios, predictions and data-averaged effects on the response scale (Section 9.4), and what changes when you add random effects, including how to test them (Section 9.5, Section 9.5.2).
Workflow steps (Chapter 12): step 3 (choice of the distribution), step 5 (residual checks), step 6 (inference on the response scale).
Prerequisites: the linear model (Chapter 2), random effects (Chapter 4) for the GLMM parts.
9.1 Introduction to GLMs
Generalized linear models (GLMs) extend the linear model (LM) to other (i.e. non-normal) distributions. The idea is the following:
We want to keep the regression formula \(y \sim f(x)\) of the lm with all its syntax, including splines (Section 2.5.3), random effects (Chapter 4) and all that. In the GLM world, this is called the “linear predictor”.
However, in the GLM, we now allow other residual distributions than normal, e.g. binomial, Poisson, Gamma, etc.
Some families require a particular range of the predictions (e.g. binomial requires predictions between zero and one). To achieve this, we use a so-called link function to bring the results of the linear predictor to the desired range.
In R, the distribution is specified via the family() function, which distinguishes the glm from the lm function. If you look at the help of family, you will see that the link function is an argument of family(). If no link is provided, the default link for the respective family is chosen.
A full GLM structure is thus:
\[y \sim family[link^{-1}(f(x) + RE)]\]
Note
As you see, I noted the link function as \(link^{-1}\) in this equation. The reason is that traditionally, the link function is applied to the (expected value of the) left hand side of the equation, i.e.
\[link(E[y]) = f(x)\]
This is important to keep in mind when interpreting the names of the link function - the log link, for example, means that \(log(E[y]) = f(x)\), which actually means that we assume \(E[y] = exp(f(x))\), i.e. the result of the linear predictor enters the exponential function, which assures that we have strictly positive predictions.
The function \(f(x)\) itself can have all components that we discussed before, in particular
You can add random effects as in Section 4.2.3 (using functions
lme4::glmerorglmmTMB::glmmTMB)You can also use splines using
mgcv::gam(see Section 2.5.3)
9.2 Important GLM variants
The most important are
- Bernoulli / Binomial family = logistic regression with logit link
- Poisson family = Poisson regression with log link
- Beta regression for continuous proportions in (0, 1)
Of course, there are many additional distributions that you could consider for your response. Here an overview of the common choices:
Screenshot taken from Wikipedia. Content licensed under the Creative Commons Attribution-ShareAlike License 3.0.
9.2.1 Count data - Poisson regression
The standard model for count data (0, 1, 2, 3, …) is the Poisson regression, which implements
A log-link function -> \(y = exp(ax + b)\)
A Poisson distribution
As an example, we use the data set birdfeeding from the EcoData package - the dataset consists of observations of food given by parents to nestlings as a function of the attractiveness of the nestling.
library(EcoData)
#str(birdfeeding)
plot(feeding ~ attractiveness, data = birdfeeding)
To fit a Poisson regression to this data, we use the glm() function, specifying family = "poisson". Note again that the log link function is the default (see ?family), so it does not have to be specified.
fit = glm(feeding ~ attractiveness, data = birdfeeding, family = "poisson")
summary(fit)
Call:
glm(formula = feeding ~ attractiveness, family = "poisson", data = birdfeeding)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.47459 0.19443 7.584 3.34e-14 ***
attractiveness 0.14794 0.05437 2.721 0.00651 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 25.829 on 24 degrees of freedom
Residual deviance: 18.320 on 23 degrees of freedom
AIC: 115.42
Number of Fisher Scoring iterations: 4The output is very similar to the lm(), however, as the residuals are no longer assumed to scatter normally, all statistics based on R2 have been replaced by the deviance (defined in Note 9.1). So, we have
Deviance residuals on top
Instead of R2, we get null vs. residual deviance, and AIC. Based on the log-likelihood or the deviance, we can calculate a pseudo R2, e.g. McFadden, which is 1-[LogL(M)/LogL(M0)] (see the example box below)
Note 9.1: Different definitions of deviance
The term deviance is used with two related, but different definitions, which is a common source of confusion:
- Deviance of a GLM (general definition): \(D = 2\,[\log L(\text{saturated}) - \log L(\text{fitted})]\). The saturated model has one parameter per observation and fits the data (almost) perfectly, so \(D\) measures how far the fitted model is from the best possible fit. This is what
deviance()returns forglmobjects, and what is shown as null deviance (D of the intercept-only model) and residual deviance (D of the fitted model) insummary()and in theanova()tables of GLMs. For a normal linear model, \(D\) is the residual sum of squares (RSS; divided by \(\sigma^2\) for the scaled deviance), which is why the deviance plays the role of the sum of squares in GLMs (see Section 3.1.1 for the sums of squares in ANOVA). For 0/1 binomial data, \(\log L(\text{saturated}) = 0\), and thus \(D = -2 \log L\). - Deviance as \(-2 \log L\): in many mixed-model contexts, \(-2 \log L(\text{fitted})\) is called the deviance (older versions of
lme4did this in theiranova()tables, anddeviance()forlmerfits still returns it; newerlme4versions andglmmTMBprint it as-2*log(L)). It is the same as \(D\) without the term \(2 \log L(\text{saturated})\). AIC is built on this quantity, \(\mathrm{AIC} = -2 \log L + 2k\).
The two definitions differ by \(2 \log L(\text{saturated})\), which is the same for all models fitted to the same data. Therefore, differences between models agree for both definitions as long as the dispersion is fixed (e.g. Poisson, binomial), and likelihood ratio tests and AIC comparisons give the same results. The absolute values are different, and numbers calculated with the two definitions must not be compared. For normal linear models, differences in RSS and differences in \(-2 \log L\) are not identical, because the variance is estimated: at the maximum likelihood estimate, \(-2 \log L = n \log(\mathrm{RSS}/n) + n (1 + \log 2\pi)\).
Note
Example: we can calculate the deviance by hand, using a saturated model with one parameter per observation.
Deviance of M1
fit1 = glm(feeding ~ attractiveness, data = birdfeeding, family = "poisson")
summary(fit1)$deviance[1] 18.32001fitSat = glm(feeding ~ as.factor(1:nrow(birdfeeding)), data = birdfeeding, family = "poisson")
2*(logLik(fitSat) - logLik(fit1))'log Lik.' 18.32001 (df=25)Deviance of the null model (null model is for example needed to calculate the pseudo-R2 ):
summary(fit1)$null.deviance[1] 25.82928fitNull = glm(feeding ~ 1, data = birdfeeding, family = "poisson")
2*(logLik(fitSat) - logLik(fitNull))'log Lik.' 25.82928 (df=25)Calculation of Pseudo-R2 with deviance or logL
Common pseudo-R2 such as McFadden or Nagelkerke use the logL instead of the deviance, whereas the deviance-based pseudo-R2 (1 - Dres/Dnull, the “proportion of deviance explained”) uses the deviance. However, R2 calculated based on deviance or logL differs considerably, as the former is calculated in reference to the maximal achievable fit (saturated model) while the latter just compares the improvement of fit compared to the null model:
# Based on logL
McFadden = 1-(logLik(fit1))/(logLik(fitNull))
print(McFadden)'log Lik.' 0.06314239 (df=2)# Based on deviance
R2_dev = 1- (logLik(fitSat) - logLik(fit1)) / (logLik(fitSat) - logLik(fitNull))
print(R2_dev)'log Lik.' 0.2907272 (df=25)If we want to calculate model predictions, we have to transform to the response scale. Here we have a log link, i.e. we have to transform with exp(linear response).
exp(coef(fit)[1] + 3 * coef(fit)[2])(Intercept)
6.810034 Alternatively (and preferably), you can use the predict() function with type = "response"
dat = data.frame(attractiveness = 3)
predict(fit, newdata = dat) # linear predictor 1
1.918397 predict(fit, newdata = dat, type = "response") # response scale 1
6.810034 Because of the log link, effects are multiplicative: exp(coef(fit)) shows that each unit of attractiveness multiplies the expected number of feedings by about 1.16 (a rate ratio, see Table 12.1 for the interpretation of coefficients under the common link functions).
Effect plots work as before (see Section 2.3). Note that the effects package always transforms the y axis according to the link, so we have log scaling on the y axis, and the effect lines remain straight
library(effects)
plot(allEffects(fit))
9.2.1.1 Notes on the Poisson regression
Poisson vs. log transformed count data: For count data, even if the distribution is switched, the log link is nearly always the appropriate link function. Before GLMs were widely available, it was common to fit an lm to log-transformed counts (log(y + 1) if there are zeros). This is similar to, but not the same as, a Poisson GLM with log link: the lm assumes a constant variance on the log scale, the Poisson GLM a variance equal to the mean (see also Section 9.2.3).
Log offset: If there is a variable that controls the number of counts (e.g. time, area), this variable is usually added as an offset in the following form
fit = glm(y ~ x + offset(log(area)), family = "poisson")To see why this is useful, consider that with the log link
\[y = exp(a_0 + a_1 \cdot x + log(area) )\]
and then rewriting
\[y = exp(a_0 + a_1 \cdot x) \dot exp(log(area)) \]
i.e. this makes the expected counts proportional to area, or whatever variable is added as a log offset.
Interactions: As for all GLMs with nonlinear link functions, interpretation of the interactions is more complicated. See the section on interactions in GLMs (Section 9.4) for notes on this.
9.2.2 0/1 or k/n data - logistic regression
The standard model to fit binomial (0/1 or k/n) data is the logistic regression, which combines the binomial distribution with a logit link function. To get to know this model, let’s have a look at the titanic data set in EcoData:
library(EcoData)
#str(titanic)
#mosaicplot( ~ survived + sex + pclass, data = titanic)
titanic$pclass = as.factor(titanic$pclass)We want to analyze how survival in the titanic accident depended on other predictors. We could fit an lm, but the residual checks make it very evident that the data with a 0/1 response don’t fit to the assumption of an lm:
fit = lm(survived ~ sex * age, data = titanic)
summary(fit)
Call:
lm(formula = survived ~ sex * age, data = titanic)
Residuals:
Min 1Q Median 3Q Max
-0.8901 -0.2291 -0.1564 0.2612 0.9744
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.637645 0.046165 13.812 < 2e-16 ***
sexmale -0.321308 0.059757 -5.377 9.35e-08 ***
age 0.004006 0.001435 2.792 0.00534 **
sexmale:age -0.007641 0.001823 -4.192 3.01e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4115 on 1042 degrees of freedom
(263 observations deleted due to missingness)
Multiple R-squared: 0.3017, Adjusted R-squared: 0.2997
F-statistic: 150 on 3 and 1042 DF, p-value: < 2.2e-16par(mfrow = c(2, 2))
plot(fit)
Thus, let’s move to the logistic regression, which assumes a 0/1 response + logit link. In principle, this distribution is called Bernoulli, but in R both 0/1 and k/n are called “binomial”, as Bernoulli is the special case of binomial where n = 1.
m1 = glm(survived ~ sex*age, family = "binomial", data = titanic)
summary(m1)
Call:
glm(formula = survived ~ sex * age, family = "binomial", data = titanic)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.493381 0.254188 1.941 0.052257 .
sexmale -1.154139 0.339337 -3.401 0.000671 ***
age 0.022516 0.008535 2.638 0.008342 **
sexmale:age -0.046276 0.011216 -4.126 3.69e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1414.6 on 1045 degrees of freedom
Residual deviance: 1083.4 on 1042 degrees of freedom
(263 observations deleted due to missingness)
AIC: 1091.4
Number of Fisher Scoring iterations: 4
Note
The syntax here is for 0/1 data. If you have k/n data, you can either specify the response as cbind(k, n-k), or you can fit the glm with I(k/n) ~ x, weights = n (and family = binomial)
The coefficients are now on the logit scale (log odds, where odds = p/(1-p)); exp(coef(m1)) gives odds ratios (see Table 12.1). To transform to predictions, we have to use the inverse logit, which is in R:
cf = coef(m1)
plogis(cf[1] + cf[3] * 20) # Women, age 20.(Intercept)
0.7198448 plogis(cf[1] + cf[2] + 20 * (cf[3] + cf[4])) # Men, age 20(Intercept)
0.2430633 Alternatively, we can again use the predict function
newDat = data.frame(sex = as.factor(c("female", "male")), age = c(20,20))
predict(m1, newdata = newDat) # Linear predictor. 1 2
0.9436919 -1.1359580 predict(m1, newdata = newDat, type = "response") # Response scale. 1 2
0.7198448 0.2430633 Finally, the effect plots - note again the scaling of the y axis, which is now logit.
library(effects)
plot(allEffects(m1))
If you do an ANOVA on a glm, you should take care that you perform a Chisq and not an F-test. You notice that you have the right one if the ANOVA doesn’t supply SumSq statistics, but deviance. In the anova function, we have to set this by hand
anova(m1, test = "Chisq")Analysis of Deviance Table
Model: binomial, link: logit
Response: survived
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(>Chi)
NULL 1045 1414.6
sex 1 312.612 1044 1102.0 < 2.2e-16 ***
age 1 0.669 1043 1101.3 0.4133
sex:age 1 17.903 1042 1083.4 2.324e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova() on a single model gives sequential (type I) tests, and only for some model classes (e.g. lm, glm, gls; for glmer it gives no p-values, for glmmTMB no such method exists). The function car::Anova works for most model classes and uses a ChiSq test automatically (for glm this is a likelihood-ratio test; for glmer/glmmTMB, car::Anova computes Wald chi-squared tests, which are approximate). It gives type II tests by default (type III with type = "III", which needs sum-to-zero contrasts, Caution 9.1); as you probably want one of these anyway, prefer it.
car::Anova(m1)Analysis of Deviance Table (Type II tests)
Response: survived
LR Chisq Df Pr(>Chisq)
sex 310.044 1 < 2.2e-16 ***
age 0.669 1 0.4133
sex:age 17.903 1 2.324e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Caution 9.1: Before testing a GLM
- Check dispersion first. For Poisson and k/n binomial models, overdispersion (see the section on dispersion problems, Section 9.3.2) makes all p-values of
car::Anova()andsummary()too small. It is not unusual that a seemingly highly significant effect disappears once the dispersion is accounted for. Refit with a negative binomial or beta-binomial family (Section 9.3.2); quasi-families are discouraged (no residual checks, no AIC), and if you use one anyway, test withcar::Anova(m, test.statistic = "F"). - Coding and centering. The rules for type II vs. III, sum-to-zero contrasts and centering from the chapter on ANOVA (Section 3.2.4) also apply to GLMs; averages are taken on the link scale.
- Reporting on the response scale.
emmeans(m, ~ x, type = "response")back-transforms the means, andpairs()then gives ratios (log link) or odds ratios (logit link); for differences and data-averaged effects on the response scale, see Section 9.4.1.
For GLMMs, car::Anova() gives Wald tests that treat the estimated random-effect variances as exact and ignore the df problem of mixed models (Section 4.5.1), so with few groups the p-values are too small; for alternatives, see Table 9.1.
9.2.2.1 Notes on the logistic regression
Offset: for 0/1 data with different exposure times t, as in survival analysis, a binomial GLM with cloglog link and offset(log(t)) models a constant hazard exactly (P(event) = 1 - exp(-lambda t)). The logit link has no such interpretation.
Interactions: As for all GLMs with nonlinear link functions, interpretation of the interactions is more complicated. See the section on interactions in GLMs (Section 9.4) for notes on this.
Overdispersion: 0/1 binomial responses cannot be overdispersed, but k/n responses can be. However, 0/1 responses can show overdispersion if grouped to k/n. Note the section on residual checks (Section 9.3), in particular the dispersion problems (Section 9.3.2), as well as comments on testing binomial GLMs in the
DHARMa vignette.
More than two categories: if the response is one of several unordered categories (e.g. land-use type, type of food chosen), the logistic regression generalizes to the multinomial logistic regression (nnet::multinom; with random effects, brms with family categorical). It estimates one set of coefficients for each category, as log odds ratios relative to a reference category. The response can be one category per observation or the counts per category for each unit (e.g. number of food items of each type). If the categories are ordered (e.g. damage classes, Likert scales), use a cumulative link model (ordinal::clm, or clmm with random effects), which assumes that each predictor has the same effect on all thresholds between categories (check with ordinal::nominal_test()). If several categories can apply at the same time, the response is not multinomial, but several binary responses (see the multivariate GLMs, Section 11.6). Variants exist for choices with alternative-specific predictors (mlogit) and for matched case-control designs (survival::clogit).
TipSolution
A. Predictive analysis
library(MASS)
fit <- glm(presence ~ dist_roads + dem + ruggedness, data = elk, family = "binomial")
predictive_model = MASS::stepAIC(fit, direction = "both")Start: AIC=5109.03
presence ~ dist_roads + dem + ruggedness
Df Deviance AIC
- dist_roads 1 5101.9 5107.9
<none> 5101.0 5109.0
- ruggedness 1 5171.4 5177.4
- dem 1 5241.3 5247.3
Step: AIC=5107.94
presence ~ dem + ruggedness
Df Deviance AIC
<none> 5101.9 5107.9
+ dist_roads 1 5101.0 5109.0
- ruggedness 1 5172.0 5176.0
- dem 1 5324.8 5328.8summary(predictive_model)
Call:
glm(formula = presence ~ dem + ruggedness, family = "binomial",
data = elk)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -6.7890658 0.4917287 -13.807 <2e-16 ***
dem 0.0042951 0.0002994 14.343 <2e-16 ***
ruggedness -0.0289100 0.0035076 -8.242 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 5334.5 on 3847 degrees of freedom
Residual deviance: 5101.9 on 3845 degrees of freedom
AIC: 5107.9
Number of Fisher Scoring iterations: 4B. Causal analysis
The predictive model has actually dropped the variable of interest (distance to roads) which shows the risks of tools that select for the best predictive model such as AIC selection (see Section 8.3.1): Collinear variables that we need to adjust our effects, are often dropped.
For the causal model, we really need to think about the causal relationships between the variables:
We are interested in the effect of dist_roads on presence:
summary(glm(presence ~ dist_roads, data = elk, family = "binomial"))
Call:
glm(formula = presence ~ dist_roads, family = "binomial", data = elk)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.101e-01 6.026e-02 -6.806 1.0e-11 ***
dist_roads 3.204e-04 3.977e-05 8.056 7.9e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 5334.5 on 3847 degrees of freedom
Residual deviance: 5268.1 on 3846 degrees of freedom
AIC: 5272.1
Number of Fisher Scoring iterations: 4The effect of dist_roads is positive: the probability of finding elk increases with the distance to roads, i.e. there are fewer elk close to roads. This would fit the idea that elk avoid roads, but we can’t interpret this causally yet, because roads are not placed randomly in the landscape.
Altitude (dem) and the ruggedness probably affect both variables, presence and dist_roads (roads are mostly built in low and flat terrain), and thus they should be considered as confounders:
fit = glm(presence ~ dist_roads+ dem + ruggedness, data = elk, family = "binomial")
summary(fit)
Call:
glm(formula = presence ~ dist_roads + dem + ruggedness, family = "binomial",
data = elk)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -7.140e+00 6.153e-01 -11.604 <2e-16 ***
dist_roads -5.198e-05 5.446e-05 -0.954 0.34
dem 4.538e-03 3.941e-04 11.515 <2e-16 ***
ruggedness -2.901e-02 3.512e-03 -8.261 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 5334.5 on 3847 degrees of freedom
Residual deviance: 5101.0 on 3844 degrees of freedom
AIC: 5109
Number of Fisher Scoring iterations: 4After adjusting for altitude and ruggedness, the effect of dist_roads is no longer positive, but close to zero and n.s. Thus, the crude effect was mainly caused by the confounders: elk occur in terrain where there are fewer roads.
Let’s check the residuals:
library(DHARMa)
res <- simulateResiduals(fit)
plot(res, quantreg = TRUE) # quantreg adds the quantile lines (off by default for n > 2000)
plotResiduals(res, form = elk$dem, quantreg = TRUE)
plotResiduals(res, form = elk$ruggedness, quantreg = TRUE)
The functional forms of our confounders are not perfect.
Since we are not really interested in them, we model the confounders flexibly with a GAM, as in the section on statistical control (Section 7.2.2). Our main predictor dist_roads is still modeled as a linear effect.
library(mgcv)
fit2 <- gam(presence ~ dist_roads + s(dem) + s(ruggedness), data = elk, family = "binomial")
summary(fit2)
Family: binomial
Link function: logit
Formula:
presence ~ dist_roads + s(dem) + s(ruggedness)
Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.783e-01 8.229e-02 2.167 0.03025 *
dist_roads -1.798e-04 5.771e-05 -3.115 0.00184 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df Chi.sq p-value
s(dem) 8.283 8.845 220.3 <2e-16 ***
s(ruggedness) 8.510 8.918 128.3 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.114 Deviance explained = 9.27%
UBRE = 0.26754 Scale est. = 1 n = 3848Let’s take another look at the residual plots, in particular for the confounders.
res <- simulateResiduals(fit2)
plot(res, quantreg = TRUE) # with quantile lines
plotResiduals(res, form = elk$dem, quantreg = TRUE)
plotResiduals(res, form = elk$ruggedness, quantreg = TRUE)
Now, the quantile plots show no obvious misfit. With the confounders modeled flexibly, the effect of dist_roads is negative and significant: at the same altitude and ruggedness, elk are found slightly closer to roads. The crude effect has thus even changed its sign once the confounders were properly accounted for. Note that we did not choose the final model because the sign of the effect matched our expectations (it doesn’t), but because it includes the confounders and passes the residual checks. Judging a model by whether its estimates look plausible is a common, but dangerous, way to fool yourself.
9.2.3 Continuous positive response - Gamma regression
The Gamma regression is the natural choice for a strictly positive, continuous, right-skewed response whose variance grows with the mean (e.g. concentrations, waiting times, biomass, body sizes). A common alternative is to log-transform and fit an lm (as we do first below), but the two are not equivalent: the log-lm assumes constant variance on the log scale and, if you back-transform predictions, incurs a retransformation bias, whereas the Gamma GLM models the mean on the original scale with a variance that scales as \(\text{Var}(y) \propto \mu^2\). With a log link (which we set explicitly below, because R’s default for Gamma() is the inverse link, see the box below), coefficients are interpreted multiplicatively - exp(coef) is the factor by which the mean changes per unit of the predictor, as in the Poisson regression.
library(faraway)fit <- lm(log(resist) ~ x1 + x2 + x3 + x4, data = wafer)
summary(fit)
Call:
lm(formula = log(resist) ~ x1 + x2 + x3 + x4, data = wafer)
Residuals:
Min 1Q Median 3Q Max
-0.17572 -0.06222 0.01749 0.08765 0.10841
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.44048 0.05982 90.948 < 2e-16 ***
x1+ 0.12277 0.05350 2.295 0.042432 *
x2+ -0.29986 0.05350 -5.604 0.000159 ***
x3+ 0.17844 0.05350 3.335 0.006652 **
x4+ -0.05615 0.05350 -1.049 0.316515
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.107 on 11 degrees of freedom
Multiple R-squared: 0.8164, Adjusted R-squared: 0.7496
F-statistic: 12.22 on 4 and 11 DF, p-value: 0.0004915library(DHARMa)
simulateResiduals(fit, plot = T)
Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
Scaled residual values: 0.052 0.272 0.408 0.676 0.82 0.768 0.84 0.172 0.524 0.72 0.592 0.808 0.8 0.204 0.096 0.372Alternative: Gamma regression
fit <- glm(formula = resist ~ x1 + x2 + x3 + x4,
family = Gamma(link = "log"),
data = wafer)
summary(fit)
Call:
glm(formula = resist ~ x1 + x2 + x3 + x4, family = Gamma(link = "log"),
data = wafer)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.44552 0.05856 92.983 < 2e-16 ***
x1+ 0.12115 0.05238 2.313 0.041090 *
x2+ -0.30049 0.05238 -5.736 0.000131 ***
x3+ 0.17979 0.05238 3.432 0.005601 **
x4+ -0.05757 0.05238 -1.099 0.295248
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 0.01097542)
Null deviance: 0.69784 on 15 degrees of freedom
Residual deviance: 0.12418 on 11 degrees of freedom
AIC: 152.91
Number of Fisher Scoring iterations: 4
Note
For Gamma regression, different link functions are used. The canonical link, which is also the default in R, is the inverse link 1/x. However, also log and identity link are commonly used.
9.2.4 Continuous proportions: beta regression and alternatives
We already covered discrete proportions (k successes out of n trials), which can be modeled with a logistic regression (Section 9.2.2). Many proportions, however, are continuous: a fraction of an area, a concentration, a share of time spent in an activity. For these, there is no number of trials, so the binomial does not apply.
Before choosing a model, check how the proportion was measured:
- Counts out of a known total (e.g. germinated seeds out of seeds sown) are k/n data: use a binomial GLM, or a beta-binomial (
family = betabinomialin glmmTMB) if they are overdispersed (Douma and Weedon, 2019). - Ordered classes (e.g. 0-5%, 5-25%, 25-50%, …) are ordinal data (see the notes on the logistic regression, Section 9.2.2.1). Analyzing class midpoints as if they were measured values ignores how coarse the classes are.
- Continuous values are the topic of this section.
For continuous proportions, there are several reasonable options. The most common are presented in the following table:
| Model | R | Can handle 0/1 | Use when |
|---|---|---|---|
transformation + lm |
qlogis(y) |
only after adding a constant | quick analysis with few 0/1 (Warton and Hui, 2011); check sensitivity to the constant |
| beta regression | beta_family() |
only after squeezing into (0, 1); with many 0/1, ordered or hurdle beta are better | all values strictly between 0 and 1 |
| ordered beta | ordbeta() |
yes | 0 and 1 are the ends of the same scale (e.g. values below a detection limit) |
| hurdle beta | beta_family() + ziformula; with 1s: brms |
0 only (ziformula); 0 and 1: brms (zero_one_inflated_beta) |
the zeros arise from a different process than the other values |
The logit transform is a quick fix. The standard distribution for continuous proportions, however, is the beta distribution. It is bounded between 0 and 1 and has two parameters: the mean and the precision \(\phi\). A large \(\phi\) means that values scatter little around the mean. glmmTMB prints \(\phi\) as “Dispersion parameter”, which is confusing: for the beta, a large value means a narrow distribution. Like the binomial, the beta is naturally heteroskedastic: its variance is largest for a mean of 0.5 and shrinks towards 0 and 1. As in the logistic regression, the mean is modeled with a logit link, so coefficients are on the logit scale.
TipExample: nitrogen in leaves
We use the elemental data from EcoData (?elemental), which contain the nitrogen content of leaves as a proportion (N_dec; N_arc is its arcsine square root transformation, included in the data). There are no zeros or ones, so a standard beta regression applies.
The arcsine transformation was long the standard approach for proportions, but Warton and Hui (2011) showed that it should be avoided; we fit it here only for comparison:
abba <- elemental[elemental$Species == "ABBA", ]
m1 <- lm(N_arc ~ Year + Site, data = abba)
summary(m1)
Call:
lm(formula = N_arc ~ Year + Site, data = abba)
Residuals:
Min 1Q Median 3Q Max
-0.016553 -0.006151 -0.001069 0.004322 0.028149
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.719e+01 4.323e+00 -3.975 0.000127 ***
Year 8.570e-03 2.145e-03 3.996 0.000117 ***
SiteDP 3.880e-03 5.972e-03 0.650 0.517275
SiteTN 1.043e-02 6.191e-03 1.685 0.094856 .
SiteUNI 4.646e-04 5.971e-03 0.078 0.938124
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.00826 on 109 degrees of freedom
Multiple R-squared: 0.2968, Adjusted R-squared: 0.271
F-statistic: 11.5 on 4 and 109 DF, p-value: 7.936e-08Now a beta regression:
library(glmmTMB)
m2 <- glmmTMB(N_dec ~ Year + Site, family = beta_family, data = abba)
summary(m2) Family: beta ( logit )
Formula: N_dec ~ Year + Site
Data: abba
AIC BIC logLik -2*log(L) df.resid
-1148.0 -1131.6 580.0 -1160.0 108
Dispersion parameter for beta family (): 3.9e+03
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -361.33063 82.66891 -4.371 1.24e-05 ***
Year 0.17684 0.04101 4.313 1.61e-05 ***
SiteDP 0.08905 0.12851 0.693 0.4883
SiteTN 0.22115 0.13190 1.677 0.0936 .
SiteUNI 0.01559 0.12865 0.121 0.9035
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(simulateResiduals(m2))
Both models find an increase of leaf nitrogen over the years. The beta regression models the bounded response directly, and its residuals show no problems.
Warton, D.I. and Hui, F.K.C. (2011) “The arcsine is asinine: The analysis of proportions in ecology,” Ecology, 92(1), pp. 3–10.
Kubinec, R. (2023) “Ordered beta regression: A parsimonious, well-fitting model for continuous data with lower and upper bounds,” Political Analysis, 31(4), pp. 519–536.
Smithson, M. and Verkuilen, J. (2006) “A better lemon squeezer? Maximum-likelihood regression with beta-distributed dependent variables,” Psychological Methods, 11(1), pp. 54–71.
Dealing with 0/1 observations: The beta distribution is defined only strictly between 0 and 1, so exact zeros or ones need an extension. The hurdle beta (in glmmTMB: zero-inflated beta) works like zero-inflation for counts (Section 9.3.3), but because a beta never produces an exact 0, the zero part is a separate logistic model for the zeros. A predictor can therefore affect whether a value is zero and how large it is otherwise in different ways (Douma and Weedon, 2019). The ordered beta (Kubinec, 2023) uses a single linear predictor plus two thresholds on the logit scale that decide when a value is exactly 0 or 1. It assumes that a predictor that raises the mean of the continuous part also lowers the probability of a 0, and it needs fewer parameters when there are several predictors. The ordered beta was introduced only in 2022, so applied studies more often use zero/one-inflated or hurdle models. Other options are to squeeze the data into (0, 1) before a beta regression (Smithson and Verkuilen, 2006), a quasibinomial GLM (no likelihood, so no AIC and no DHARMa residuals), or quantile regression (Section 10.3.2). The beta models in glmmTMB (beta, ordered beta, hurdle beta) can include random effects, and their precision can depend on predictors (dispformula, Section 10.1.2). An example with zeros follows in the case study on plant cover (Section 9.7.2).
9.3 Residuals and their solutions in GLMs
First of all: everything we said about residual checks for LMs (Section 2.4) and about model selection (Chapter 8) also applies for GLMs, with only very few additions, so you should check your model in principle as before. However, there are a few tweaks you have to be aware of.
Let’s look again at the titanic example
m1 = glm(survived ~ sex*age, family = "binomial", data = titanic)How can we check the residuals of this model? The standard residual plots still run (Nerds: check class(m1)), but they look horrible, because they check for normally distributed residuals, whereas ours should be binomial:
par(mfrow = c(2, 2))
plot(m1)
9.3.1 DHARMa residual plots for GL(M)Ms
The DHARMa residuals solve this problem: they are simulated from the fitted model, so they work for any family and are read as for the linear model (Section 2.4.2).
library(DHARMa)
res = simulateResiduals(m1)Standard plot:
plot(res)
The Q-Q plot (left) follows the 1-1 line and the residuals vs. predicted plot (right) shows no pattern, so these two plots show no problem. The residuals against the predictors also show no particular problem:
par(mfrow = c(1, 2))
plotResiduals(res, form = model.frame(m1)$age)
plotResiduals(res, form = model.frame(m1)$sex)
However, residuals against the missing predictor pclass show a clear problem:
dataUsed = as.numeric(rownames(model.frame(m1)))
plotResiduals(res, form = titanic$pclass[dataUsed])
Thus, I should add passenger class to the model
m2 = glm(survived ~ sex*age + pclass, family = "binomial", data = titanic)
summary(m2)
Call:
glm(formula = survived ~ sex * age + pclass, family = "binomial",
data = titanic)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.790839 0.362822 7.692 1.45e-14 ***
sexmale -1.029755 0.358593 -2.872 0.00408 **
age -0.004084 0.009461 -0.432 0.66598
pclass2 -1.424582 0.241513 -5.899 3.67e-09 ***
pclass3 -2.388178 0.236380 -10.103 < 2e-16 ***
sexmale:age -0.052891 0.012025 -4.398 1.09e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1414.62 on 1045 degrees of freedom
Residual deviance: 961.92 on 1040 degrees of freedom
(263 observations deleted due to missingness)
AIC: 973.92
Number of Fisher Scoring iterations: 5res2 = simulateResiduals(m2)
plotResiduals(res2, form = model.frame(m2)$pclass)
Now, residuals look fine. Of course, if your model gets more complicated, you may want to do additional checks, for example for the distribution of random effects etc.
9.3.2 Dispersion problems in GLMs
One thing that is different between GLMs and LM is that GLMs can display overall dispersion problems. The most common GLMs to show overdispersion are the Poisson and the binomial with k/n data (0/1 data cannot be overdispersed).
The reason is that simple GLM distributions such as the Poisson or the Binomial (for k/n data) do not have a parameter for adjusting the spread of the observed data around the regression line (dispersion), but their variance is a fixed function of the mean.
There are good reasons for why this is the case (Poisson and Binomial describe particular processes, e.g. coin flip, for which the variance is a fixed function of the mean), but the fact is that when applying these GLMs on real data, we often find overdispersion (more dispersion than expected), and more rarely, underdispersion (less dispersion than expected).
To remove the assumptions of a fixed dispersion, there are three options, of which I usually prefer the third one:
- Quasi-distributions, which are available in glm. Those add a term to the likelihood that corrects the p-values for the dispersion, but they are not distributions, so there are no residual checks and no AIC. Therefore discouraged.
- Observation-level random effect (OLRE) - Add a separate random effect per observation. This effectively creates a normal random variate at the level of the linear predictor, increases variance on the responses.
- A GLM distribution with variable dispersion, for Poisson usually the negative binomial, for k/n binomial data the beta-binomial (
family = betabinomialin glmmTMB).
The reason why I prefer the 3rd option is that it allows better residual checks and to model the dispersion as a function of the predictors (see Section 10.2). OLRE is also a valid solution.
Note
Overdispersion is often created by model misfit. Thus, before moving to a variable dispersion GLM, you should check for / correct model misfit.
9.3.2.1 Recognizing overdispersion
To understand how to recognize overdispersion, let’s look at an example. We’ll use the Salamanders dataset from the package glmmTMB, staring with a simple Poisson glm:
library(glmmTMB)
library(lme4)
library(DHARMa)
m1 = glm(count ~ spp + mined, family = poisson, data = Salamanders)Overdispersion will be automatically highlighted in the standard DHARMa plots
res = simulateResiduals(m1, plot = T)
You see the dispersion problem by:
Dispersion test in the left plot significant
QQ plot S-shaped
Quantile lines in the right plots outside their expected quantiles
You can get a more detailed output with the testDispersion function, which also displays the direction of the dispersion problem (over or underdispersion)
testDispersion(res)
DHARMa nonparametric dispersion test via sd of residuals fitted vs.
simulated
data: simulationOutput
dispersion = 3.9152, p-value < 2.2e-16
alternative hypothesis: two.sidedOK, often the dispersion problem is caused by structural problems. Let’s add a random effect for site, which makes sense. We can do so using the function lme4::glmer, which adds the same extensions to glm as lmer adds to lm. This means we can use the usual random effect syntax from Section 4.2.1.
m2 = glmer(count ~ spp + mined + (1|site),
family = poisson, data = Salamanders)Testing the residuals
res = simulateResiduals(m2, plot = T)
shows that there is still some overdispersion. Actually, what the plots also show is heteroskedasticity, and we should probably deal with that as well, but we will only learn how in the next chapter on dispersion models (Chapter 10). For now, let’s switch to a negative binomial model. This could be fit with lme4, but it is more convenient to use the package glmmTMB, which has the same syntax as lme4, but more advanced options.
m3 = glmmTMB(count ~ spp + mined + (1|site), family = nbinom2, data = Salamanders)
res = simulateResiduals(m3, plot = T)
Now, the test is fine. The negative binomial is a Poisson distribution with a gamma-distributed rate. glmmTMB offers two parameterizations: nbinom2 (variance \(\mu + \mu^2/\theta\), quadratic in the mean, the usual starting point) and nbinom1 (variance \(\mu (1 + \alpha)\), linear in the mean); compare both by AIC if in doubt.
9.3.3 Zero-inflation
Another common problem in count data (Poisson / negative binomial), but also other GLMs (e.g. binomial, beta) is that the observed data has more zeros than expected by the fitted distribution. For the beta, 1-inflation, and for the k/n binomial, n-inflation is also common, and tested for / addressed in the same way.
To deal with this zero-inflation, one usually adds an additional model component that controls how many zeros are produced. The default way to do this is assuming two separate processes which act after one another:
- First, we have the normal GLM, predicting what values we would expect
- On top of that, we have a logistic regression, which decides whether the GLM prediction or a zero should be observed
Note that the result of 1 can also be zero, so there are two explanations for a zero in the data. Zero-inflated GLMMs can, for example, be fit with glmmTMB, using the ziformula argument.
9.3.3.1 Recognizing zero-inflation
Caution
The fact that you have a lot of zeros in your data does not indicate zero-inflation. Zero-inflation is with respect to the fitted model. You can only check for zero-inflation after fitting a model.
Let’s look at our last model - DHARMa has a special function to check for zero-inflation
testZeroInflation(res)
DHARMa zero-inflation test via comparison to expected zeros with
simulation under H0 = fitted model
data: simulationOutput
ratioObsSim = 1.0301, p-value = 0.28
alternative hypothesis: two.sidedThis shows no sign of zero-inflation. There are, however, two problems with this test:
- For mixed models, the test is less powerful with unconditional residuals, because the re-simulated random effects add variation (Tip 4.4). This happens if you set
simulateREs = "unconditional", or if DHARMa cannot simulate conditionally for the model’s covariance structure and falls back to unconditional simulations (it then says so in a message). - When there is really zero-inflation, variable dispersion models such as the negative Binomial often simply increase the dispersion to account for the zeros, leading to no apparent zero-inflation in the residuals, but rather underdispersion.
Thus, for zero-inflation, model selection, or simply fitting a ZIP model is often more reliable than residual checks. Compare the models with and without ziformula via AIC or likelihood ratio test. Do not use the p-value of the zero-inflation intercept for this, as it tests a zero-inflation probability of 0.5, not the absence of zero-inflation.
m4 = glmmTMB(count ~ spp + mined + (1|site), family = nbinom2, ziformula = ~1, data = Salamanders)
summary(m4) Family: nbinom2 ( log )
Formula: count ~ spp + mined + (1 | site)
Zero inflation: ~1
Data: Salamanders
AIC BIC logLik -2*log(L) df.resid
1674.4 1723.5 -826.2 1652.4 633
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
site (Intercept) 0.2944 0.5426
Number of obs: 644, groups: site, 23
Dispersion parameter for nbinom2 family (): 0.942
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.6832 0.2742 -6.140 8.28e-10 ***
sppPR -1.3197 0.2875 -4.591 4.42e-06 ***
sppDM 0.3686 0.2235 1.649 0.099047 .
sppEC-A -0.7098 0.2530 -2.806 0.005016 **
sppEC-L 0.5714 0.2191 2.608 0.009105 **
sppDES-L 0.7929 0.2166 3.660 0.000252 ***
sppDF 0.3120 0.2329 1.340 0.180329
minedno 2.2633 0.2838 7.975 1.53e-15 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -16.41 4039.11 -0.004 0.997In this case, we have no evidence for zero-inflation (AIC 1674.4 with vs. 1672.4 without zero-inflation, i.e. the simpler model is preferred; the estimated zero-inflation intercept of about -16 on the logit scale, i.e. a probability of about 0, with a huge standard error is a boundary estimate). To see an example where you can find zero-inflation, see the owls exercise (Caution 9.2).
9.4 Interpreting effects and interactions in GLMs
A significant problem with interpreting GLMs is the interpretation of effects, in particular of interactions. To see the problem, let’s simulate survival data in which the treatment and the concentration of a substance act additively on the probability scale: the treatment raises survival from 30% to 80%, and the concentration adds up to 15% in both groups. We might think that there is no interaction, but a logistic regression finds a clearly positive interaction coefficient (z ≈ 5):
library(effects)
set.seed(1)
n = 20000
trt = as.factor(sample(c("ctrl", "trt"), n, replace = TRUE))
concentration = runif(n)
response = 0.3 + 0.5 * (as.numeric(trt) - 1) + 0.15 * concentration
survival = rbinom(n, 1, prob = response)
dat = data.frame(trt = trt, concentration = concentration, survival = survival)
m2 = glm(survival ~ trt * concentration, data = dat, family = "binomial")
summary(m2)
Call:
glm(formula = survival ~ trt * concentration, family = "binomial",
data = dat)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.78369 0.04207 -18.626 < 2e-16 ***
trttrt 2.11054 0.06892 30.624 < 2e-16 ***
concentration 0.57043 0.07198 7.925 2.28e-15 ***
trttrt:concentration 0.69485 0.12958 5.362 8.22e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 26521 on 19999 degrees of freedom
Residual deviance: 20813 on 19996 degrees of freedom
AIC: 20821
Number of Fisher Scoring iterations: 4plot(allEffects(m2))
The reason is the link function: a coefficient is a constant change on the logit scale, and the same change on the logit scale translates into different changes of the probability, depending on where we start. If we plot the plogis() function, this becomes obvious: at a baseline near 0.4 (control), one unit on the logit scale changes the probability by a lot, at a baseline near 0.9 (treatment) by much less.
Turned around: for the same effect of concentration on the probability scale in both groups, the model needs an interaction on the logit scale. The effect plots show this as well, but note that they scale the y axis according to the link, so the lines look straight and the slopes look different. Whether a “real” interaction is one on the link scale or one on the response scale is a biological decision: one could argue that a change of 15% in survival means something different at 30% than at 80% survival. Decide this from your question before looking at the results; a link-scale interaction can appear or disappear depending on the scale, and choosing the scale after seeing which one gives the result you like is a forking path (Section 8.3).
Note
If effect directions change in sign, they will do so under any link function (as they are always monotone), so changes in effect direction are robust to this problem.
9.4.1 Ways to report an effect
The effect of a predictor on the response scale depends on all other predictors. As in the linear model (Section 2.11.2; effects plots to look at the model, emmeans for the numbers), you can report the curve or predictions at chosen values, slopes or differences at chosen values, and the data-averaged effect; in a GLM, a fourth option appears, the ratio, which is constant on the link scale. In a linear model without interactions or curvature, the slope at the mean equals the data-averaged slope; in a GLM it does not.
| What | Scale | Averaged over | Call |
|---|---|---|---|
| Curve or predictions at chosen values | response | nothing | effects plots; emmeans(m, ~ trt, at = list(...), type = "response") |
| Slopes or differences at chosen values | response | nothing | emtrends(m, ~ trt, var = "x", regrid = "response"); pairs(regrid(emmeans(...))) |
| Data-averaged effect: the difference (factor) or slope (numeric predictor) on the response scale, averaged over all observations | response | the observed data | factors: emmeans(m, ~ trt, counterfactuals = "trt", type = "response"); slopes: finite differences |
| Ratio: rate ratio (log link), odds ratio (logit link) | link scale, back-transformed | nothing (constant without interactions; with interactions it depends on the other predictors) | exp(coef(m)); pairs(emmeans(m, ~ trt, type = "response")) |
Ratios. For our model, the odds ratio of the treatment is
library(emmeans)Welcome to emmeans.
Caution: You lose important information if you filter this package's results.
See '? untidy'pairs(emmeans(m2, ~ trt, type = "response"), reverse = TRUE)NOTE: Results may be misleading due to involvement in interactions contrast odds.ratio SE df null z.ratio p.value
trt / ctrl 11.7 0.437 Inf 1 65.583 <0.0001
Tests are performed on the log odds ratio scale i.e. the odds of survival are about 12 times higher under treatment, evaluated at the mean concentration (emmeans may warn that the result can be misleading because of the interaction: the odds ratio changes with the concentration). Ratios arise because a logit-link model is additive in log-odds, and adding on a log scale means multiplying: a coefficient becomes a factor on the odds (with a log link, a factor on the expected count, the rate ratio). Ratios are compact and constant on the link scale, which is why they are the standard way to report GLM coefficients (Table 12.1), but they are hard to read: an odds ratio is not a risk ratio, and an odds ratio of 12 does not tell you whether survival went from 1% to 11% or from 30% to 84%.
Predictions on the response scale at chosen values of the other predictors answer this directly:
emmeans(m2, ~ trt | concentration, at = list(concentration = c(0.2, 0.8)), type = "response")concentration = 0.2:
trt prob SE df asymp.LCL asymp.UCL
ctrl 0.339 0.00680 Inf 0.325 0.352
trt 0.829 0.00546 Inf 0.818 0.840
concentration = 0.8:
trt prob SE df asymp.LCL asymp.UCL
ctrl 0.419 0.00714 Inf 0.405 0.433
trt 0.912 0.00407 Inf 0.904 0.920
Confidence level used: 0.95
Intervals are back-transformed from the logit scale The confidence intervals are back-transformed from the logit scale and therefore asymmetric and bounded by 0 and 1.
The data-averaged effect answers the question “by how much did the treatment change survival on average in this sample”. With counterfactuals, emmeans predicts every observation once as control and once as treated (with its own concentration), averages the predictions on the probability scale and compares them:
emmC = emmeans(m2, ~ trt, counterfactuals = "trt", type = "response")NOTE: Results may be misleading due to involvement in interactionsemmC trt prob SE df asymp.LCL asymp.UCL
ctrl 0.378 0.00481 Inf 0.369 0.388
trt 0.871 0.00335 Inf 0.864 0.877
Results are averaged over the levels of: actual_trt
Confidence level used: 0.95 confint(contrast(emmC, "revpairwise")) contrast estimate SE df asymp.LCL asymp.UCL
trt - ctrl 0.492 0.00586 Inf 0.481 0.504
Results are averaged over the levels of: actual_trt
Confidence level used: 0.95 The treatment raised survival by 0.49 on average (the true value is 0.50). This quantity is the G-computation (also called standardization) estimate of an average causal effect (Chapter 7), it refers to the covariate distribution of this sample, its standard error assumes that the model is right (check the residuals first, Section 9.3), and the predictor must be a factor (a 0/1 predictor is treated as numeric). The output label actual_trt refers to the observed treatment, over which the predictions are averaged.
For a numeric predictor, the data-averaged slope on the probability scale follows from finite differences as in Section 2.11.2, here per treatment group:
h = 0.01
d = (predict(m2, transform(dat, concentration = concentration + h), type = "response") -
predict(m2, transform(dat, concentration = concentration - h), type = "response")) / (2 * h)
tapply(d, dat$trt, mean) ctrl trt
0.1333440 0.1402435 emtrends() with regrid = "response" gives instead the slope of the probability at the mean of the other predictors (shown in the next subsection); in our example, the two differ by less than 0.01. The marginaleffects package (Arel-Bundock, Greifer and Heiss, 2024) computes such average effects and slopes directly for many model classes; in this book, we stay with emmeans.
Arel-Bundock, V., Greifer, N. and Heiss, A. (2024) “How to interpret statistical models using marginaleffects for R and Python,” Journal of Statistical Software, 111(9), pp. 1–32.
Which one to report? Always show the curve or the predictions, as for linear models; add one summary number. In a designed experiment with few covariates, a ratio or the difference at reference values is fine; for observational data, interactions or several covariates, the data-averaged effect is usually the most useful single number. For GLMMs, counterfactuals gives the data-averaged effect for a typical group, not a marginal one (Section 9.5.1). In all cases, say which scale and which averaging you used, e.g. “treatment increased survival by 0.49 (95% CI 0.48 to 0.50), averaged over the observed concentrations (data-averaged effect)”, or “the odds of survival were 12 times higher under treatment (95% CI 11 to 13) at the mean concentration”.
9.4.2 Interactions on the response scale
With these tools, we can now judge the interaction of our example on the scale we decided on. On the logit scale, the slopes of concentration clearly differ between the groups:
emtrends(m2, ~ trt, var = "concentration") trt concentration.trend SE df asymp.LCL asymp.UCL
ctrl 0.57 0.072 Inf 0.429 0.712
trt 1.27 0.108 Inf 1.054 1.476
Confidence level used: 0.95 On the probability scale, they do not:
emtrends(m2, ~ trt, var = "concentration", regrid = "response") trt concentration.trend SE df asymp.LCL asymp.UCL
ctrl 0.134 0.0169 Inf 0.101 0.167
trt 0.137 0.0113 Inf 0.115 0.159
Confidence level used: 0.95 pairs(emtrends(m2, ~ trt, var = "concentration", regrid = "response")) contrast estimate SE df z.ratio p.value
ctrl - trt -0.00304 0.0203 Inf -0.150 0.8810So the rule is: a significant interaction on the link scale says nothing about whether the effects are additive on the response scale, and vice versa. Judge an interaction by differences of slopes or of effects on the scale that you decided on; for two factors, contrast(regrid(emmeans(m, ~ A * B)), interaction = "pairwise") gives the difference of differences on the response scale. Note also that main effects and interactions can change if you change the link function, and that default links are not always natural from a biological viewpoint.
9.5 New considerations for GLMMs
As we saw, in principle, random effects can be added to GLMs very much in the same way as before (Chapter 4). Note, however, that there is a conceptual difference:
In an LMM, we had y = f(x) + RE + residual, where both RE and residual error were normal distributions
In a GLMM, we have \(y \sim dist[link^{-1}(f(x) + RE)]\), so the normally distributed RE goes through the inverse link into the response distribution
This has a number of consequences that may be unexpected if you don’t think about it.
9.5.1 Unconditional vs. marginal predictions
A consequence of the inverse link is that, with a nonlinear link function, unconditional predictions (fixed effects only) differ from marginal predictions (averaged over the random effects) and thus from the mean of your data; the three kinds of predictions are defined in Tip 4.3. To understand this, imagine we have a logistic regression, and our fixed effects predict a mean of 1. Then the unconditional prediction at the response scale
plogis(1)[1] 0.7310586Now, let’s try to see what the average (marginal) prediction over all random effects is. We simulate random intercepts from a normal distribution with sd = 1 (1e5 draws). In this case, there will be additional variation around the random intercept.
mean(plogis(1 + rnorm(1e5, sd = 1)))[1] 0.6966768The value (about 0.70) is lower than the unconditional prediction (0.73), and this value is also what you would approximately get if you take the mean of your data.
Whether this is a problem or not is a matter of perspective. From statistical viewpoint, assuming that your model assumptions correspond to the data-generating model, there is nothing wrong with the unconditional (fixed effect) prediction not corresponding to the mean of the data. From a practical side, however, many people are unsatisfied with this, because they want to show predictions against data. A simple and transparent way to create marginal predictions is the simulation above: create unconditional predictions on the link scale, add random effects simulated with the estimated variances, push them through the corresponding link function and average. emmeans does approximately the same with bias.adj = TRUE, given the standard deviation of the random effects (its confidence intervals ignore the uncertainty of this standard deviation), here for the cbpp data of the lme4 package:
library(lme4)
m = glmer(cbind(incidence, size - incidence) ~ period + (1 | herd), data = cbpp, family = binomial)
sdHerd = attr(VarCorr(m)$herd, "stddev")
emmeans(m, ~ period, type = "response") # unconditional period prob SE df asymp.LCL asymp.UCL
1 0.1981 0.0367 Inf 0.1357 0.280
2 0.0839 0.0236 Inf 0.0478 0.143
3 0.0740 0.0224 Inf 0.0404 0.132
4 0.0484 0.0196 Inf 0.0216 0.105
Confidence level used: 0.95
Intervals are back-transformed from the logit scale emmeans(m, ~ period, type = "response", bias.adj = TRUE, sigma = sdHerd) # marginal period prob SE df asymp.LCL asymp.UCL
1 0.2179 0.0371 Inf 0.1533 0.298
2 0.0971 0.0263 Inf 0.0562 0.161
3 0.0861 0.0251 Inf 0.0477 0.149
4 0.0570 0.0225 Inf 0.0258 0.120
Confidence level used: 0.95
Intervals are back-transformed from the logit scale
Bias adjustment applied based on sigma = 0.64207 Here, “marginal” means averaged over the random effects (for the other meanings of the word, see Tip 4.3); the data-averaged effect (Section 9.4.1) averages over the observations instead, with the random effects at zero: emmeans(m, ~ period, counterfactuals = "period", type = "response") gives the data-averaged effect for a typical herd (here identical to the unconditional means, because the model has no covariates), and the bias adjustment is not applied to counterfactual averages. A data-averaged effect that is also marginal over the random effects can be obtained by the simulation approach above, applied to each observation.
9.5.2 ANOVA and tests for GLMMs
The following table summarizes our recommendations for ANOVA-type questions in GLMMs, analogous to the table for linear mixed models (Table 4.2). As there, the table only lists which ANOVA types each function supports; which type to use is discussed in Section 3.2.3.2.
glmer, glmmTMB).
| Question | Recommended function | Notes |
|---|---|---|
| Test fixed-effect terms (ANOVA table) | car::Anova(m) |
Wald \(\chi^2\) tests, for glmer and glmmTMB models. Types II and III (type = "II", "III"; default II). For type I, there is no single-model test with p-values (anova(m) of lme4 gives none, and glmmTMB has no such method); compare a sequence of nested models, each adding one term, with anova(m1, m2, m3) (likelihood ratio tests). Check dispersion first (Caution 9.1). |
| … with fewer than about 10-20 groups, or a predictor that varies only between groups | afex::mixed(formula, data, family = ..., method = "PB") |
The Wald p-values are then too small. afex::mixed() computes the whole ANOVA table with p-values from a parametric bootstrap (slow; for glmer models; types II and III, default III; with pbkrtest 0.5.5, this fails for glmer models with an error, so check that it runs and otherwise use simulateLRT()). Alternative for single terms, also for glmmTMB: DHARMa::simulateLRT(m0, m1), with m0 = the model without the term. |
| How much variation do fixed and random effects explain? (description, not a test) | lme4::VarCorr(m), MuMIn::r.squaredGLMM(m) |
Random effects can be compared among each other by the ratios of their variances from VarCorr() (on the link scale). A GLMM has no residual variance on the link scale, so VPCs as for linear mixed models (Section 4.5.5) and R2 values are only approximate (they use a distribution-specific residual variance, e.g. \(\pi^2/3\) for the logit link); r.squaredGLMM() returns one row per approximation of this residual variance: “delta” (all families), “lognormal” and “trigamma” (log links; the package help recommends “trigamma” whenever available) and “theoretical” (binomial). Report “trigamma” for log links; otherwise state which row you report. |
| Test a random effect | DHARMa::simulateLRT(m0, m1) |
m0: the same fixed effects without the random effect (glm() if no random effect is left) (Section 4.5.3.3). If p = 0, report p < 1/(n+1) (Section 3.2.2). |
9.6 Take-home messages
Tip
- Choose the distribution from the type of response data (counts, 0/1, k/n, positive continuous, proportions), and then check the residuals of the fitted model.
- The estimates are on the scale of the link function. Report effects on an interpretable scale: rate ratios for log links, odds ratios for logit links, back-transformed predictions with CIs (
emmeans(..., type = "response")), or the data-averaged effect (counterfactuals); say which scale and averaging you used (Section 9.4.1). - Check dispersion before you interpret p-values: with overdispersed Poisson or k/n binomial data, p-values are too small. Overdispersion is often caused by model misfit, so check the model structure first; then use a negative binomial or beta-binomial family (if needed with a dispersion model that depends on predictors), or an observation-level random effect.
- Many zeros do not mean zero-inflation; zero-inflation is relative to the fitted model. Compare models with and without
ziformulainstead of relying on the residual test alone. - With a nonlinear link, interactions and effects depend on the scale: decide the scale from the biology before looking at the results, and judge interactions by slopes or differences on that scale. In a GLMM, predictions from the fixed effects alone (unconditional) need not match the average prediction over the random effects (marginal) or the mean of your data.
9.7 Case studies
The following case studies apply the models of this chapter to real data, including residual checks and fixes for dispersion and zero-inflation problems (Section 9.3). Try to solve the exercises yourself before looking at the solutions.
9.7.1 Count data - owls
Caution 9.2: Exercise — Owls
Look at the Owls data set in the glmmTMB package (?Owls). We want to know whether the begging of owl nestlings (SiblingNegotiation, the number of calls) depends on the food treatment and on the sex of the parent. A first model is
library(glmmTMB) # for the Owls data
m1 = glm(SiblingNegotiation ~ FoodTreatment*SexParent + offset(log(BroodSize)),
data = Owls , family = poisson)
res = simulateResiduals(m1)
plot(res)
The offset log(BroodSize) makes the expected number of calls proportional to the number of nestlings (see the log offset in Section 9.2.1.1).
Try to improve the model with everything we have discussed so far (offset, dispersion, zero-inflation, random effects). Hint: besides nbinom2, glmmTMB offers nbinom1, and the dispersion can depend on predictors (dispformula, see Section 10.1.2).
TipSolution
m1 = glmmTMB::glmmTMB(SiblingNegotiation ~ FoodTreatment * SexParent
+ (1|Nest) + offset(log(BroodSize)), data = Owls , family = nbinom1,
dispformula = ~ FoodTreatment + SexParent,
ziformula = ~ FoodTreatment + SexParent)
summary(m1) Family: nbinom1 ( log )
Formula:
SiblingNegotiation ~ FoodTreatment * SexParent + (1 | Nest) +
offset(log(BroodSize))
Zero inflation: ~FoodTreatment + SexParent
Dispersion: ~FoodTreatment + SexParent
Data: Owls
AIC BIC logLik -2*log(L) df.resid
3354.6 3402.9 -1666.3 3332.6 588
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
Nest (Intercept) 0.0876 0.296
Number of obs: 599, groups: Nest, 27
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.80028 0.09736 8.220 < 2e-16 ***
FoodTreatmentSatiated -0.46893 0.16760 -2.798 0.00514 **
SexParentMale -0.09127 0.09247 -0.987 0.32363
FoodTreatmentSatiated:SexParentMale 0.13087 0.19028 0.688 0.49158
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.9132 0.3269 -5.853 4.84e-09 ***
FoodTreatmentSatiated 1.0564 0.4072 2.594 0.00948 **
SexParentMale -0.4688 0.3659 -1.281 0.20012
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Dispersion model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.2122 0.2214 5.475 4.37e-08 ***
FoodTreatmentSatiated 0.7978 0.2732 2.920 0.0035 **
SexParentMale -0.1540 0.2399 -0.642 0.5209
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1res = simulateResiduals(m1, plot = T)
testDispersion(m1)
DHARMa nonparametric dispersion test via sd of residuals fitted vs.
simulated
data: simulationOutput
dispersion = 0.87257, p-value = 0.192
alternative hypothesis: two.sidedtestZeroInflation(m1)
DHARMa zero-inflation test via comparison to expected zeros with
simulation under H0 = fitted model
data: simulationOutput
ratioObsSim = 1.0594, p-value = 0.352
alternative hypothesis: two.sidedUsing brms (optional, requires Stan and takes several minutes)
The same model can also be fit with brms (here without dispersion and zero-inflation terms, just to show the syntax). A step-by-step Bayesian GLMM with posterior predictive checks, overdispersion and zero-inflation (in JAGS) is developed in the chapter on GLMMs of Hartig (2024). The brms code:
library(brms)
m2 = brms::brm(SiblingNegotiation ~ FoodTreatment * SexParent
+ (1|Nest) + offset(log(BroodSize)),
data = Owls ,
family = negbinomial)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.000135 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.35 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.84 seconds (Warm-up)
Chain 1: 0.679 seconds (Sampling)
Chain 1: 1.519 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 5.5e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.55 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 31.439 seconds (Warm-up)
Chain 2: 40.276 seconds (Sampling)
Chain 2: 71.715 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 6.8e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.68 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 1.195 seconds (Warm-up)
Chain 3: 0.812 seconds (Sampling)
Chain 3: 2.007 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 6.9e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.69 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.846 seconds (Warm-up)
Chain 4: 0.652 seconds (Sampling)
Chain 4: 1.498 seconds (Total)
Chain 4: summary(m2) Family: negbinomial
Links: mu = log
Formula: SiblingNegotiation ~ FoodTreatment * SexParent + (1 | Nest) + offset(log(BroodSize))
Data: Owls (Number of observations: 599)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Multilevel Hyperparameters:
~Nest (Number of levels: 27)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.30 0.19 0.00 0.59 1.59 7 11
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept -7.91 14.95 -33.79 0.99 1.57
FoodTreatmentSatiated 0.44 2.12 -1.09 4.11 1.59
SexParentMale 0.25 0.50 -0.30 1.08 1.59
FoodTreatmentSatiated:SexParentMale 1.12 1.68 -0.21 4.01 1.53
Bulk_ESS Tail_ESS
Intercept 7 12
FoodTreatmentSatiated 7 12
SexParentMale 7 11
FoodTreatmentSatiated:SexParentMale 7 45
Further Distributional Parameters:
Estimate Est.Error
shape 259394536299280214785340135702528.00 449340687172916610793870660206592.00
l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 0.72 1037578145197123597329933984071680.00 1.53 7 NA
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).plot(m2, ask = FALSE)
Hartig, F. (2024) “Introduction to Bayesian statistics with R.” Available at: https://florianhartig.github.io/LearningBayes/.
9.7.2 Continuous proportions - plant cover
Plant cover is one of the most common continuous proportions in ecology, and it comes with some special issues (Damgaard and Irvine, 2019):
Damgaard, C.F. and Irvine, K.M. (2019) “Using the beta distribution to analyse plant cover data,” Journal of Ecology, 107(6), pp. 2747–2759.
- How was cover measured? Pin-point (point-intercept) data are hits out of N pins, i.e. k/n data; as plants grow in clumps, they are often overdispersed (beta-binomial). Cover classes (Braun-Blanquet, Londo, …) are ordinal. Visual estimates in % are continuous, but often rounded to 1, 5 or 10%.
- Zeros: the choice between ordered and hurdle beta depends on what a zero means (Section 9.2.4): “too little to be recorded” suits the ordered beta, “the species cannot grow there” the hurdle beta.
- Totals: because plants overlap and there is bare ground, the cover of all species in a plot does not sum to 100%. Only relative cover (share of total vegetation) sums to 1.
- Several species: the cover of species in the same plot is usually correlated (Section 9.2.4.1). For joint models of many species, Korhonen et al. (2024) found the hurdle beta most accurate in simulations, with no noticeable difference to ordered beta for real data.
Korhonen, P. et al. (2024) “A comparison of joint species distribution models for percent cover data,” Methods in Ecology and Evolution, 15(12), pp. 2359–2372.
The data set varespec (package vegan) contains visually estimated cover (in %) of 44 species in 24 reindeer-grazed lichen pastures in Finland; varechem contains soil chemistry for the same plots. We ask whether the cover of heather (Calluna vulgaris, Callvulg) depends on soil aluminum. Heather is absent in 9 of the 24 plots:
library(vegan)
data(varespec)
data(varechem)
# Al is scaled (per SD); [, 1] drops the attributes, so that predict() works
# with new values given in SD units, e.g. Al = c(-1, 0, 1)
d <- data.frame(cover = varespec$Callvulg / 100, Al = scale(varechem$Al)[, 1])
table(d$cover == 0)
FALSE TRUE
15 9 A beta regression fails because of the zeros:
glmmTMB(cover ~ Al, family = beta_family, data = d)Error:
! y values must be 0 < y < 1A logit transform needs a constant to deal with the zeros. The result depends on this arbitrary choice:
sapply(c(0.001, 0.01, 0.05), function(eps) {
coef(lm(qlogis((d$cover + eps) / (1 + 2 * eps)) ~ Al, data = d))["Al"]
}) Al Al Al
0.4871682 0.2695578 0.1447610 The ordered beta handles the zeros directly:
m_ord <- glmmTMB(cover ~ Al, family = ordbeta(), data = d)
summary(m_ord) Family: ordbeta ( logit )
Formula: cover ~ Al
Data: d
AIC BIC logLik -2*log(L) df.resid
-46.8 -40.9 28.4 -56.8 19
Dispersion parameter for ordbeta family (): 12.7
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.5472 0.4216 -8.413 <2e-16 ***
Al 0.2571 0.2370 1.085 0.278
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(simulateResiduals(m_ord))
The hurdle beta separates presence and cover. Note that the zero part (ziformula) models the probability of absence:
m_hurdle <- glmmTMB(cover ~ Al, ziformula = ~ Al, family = beta_family, data = d)
summary(m_hurdle) Family: beta ( logit )
Formula: cover ~ Al
Zero inflation: ~Al
Data: d
AIC BIC logLik -2*log(L) df.resid
-49.4 -43.5 29.7 -59.4 19
Dispersion parameter for beta family (): 12.6
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.45345 0.41401 -8.341 <2e-16 ***
Al 0.02052 0.30000 0.068 0.945
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.6527 0.4861 -1.343 0.1794
Al -0.9854 0.5917 -1.666 0.0958 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The data hint that more aluminum raises the probability that heather is present (negative coefficient on absence, p = 0.10), with no sign of an effect on its cover where it is present.
For the expected cover at low, mean and high aluminum (-1, 0, 1 SD), we can use predict() for the hurdle model:
newd <- data.frame(Al = c(-1, 0, 1))
predict(m_hurdle, newdata = newd, type = "response")[1] 0.01255334 0.02016658 0.02619169For the ordered beta, predict(type = "response") gives only the mean of the continuous part, ignoring the zeros; the expected cover has to be computed from the thresholds with a small helper function (see the box):
TipExpected values for the ordered beta
The ordered beta has a linear predictor \(\eta\) and two thresholds \(\psi_1 < \psi_2\) (stored as psi). The probability of an exact 0 is \(1 - \text{logit}^{-1}(\eta - \psi_1)\), of an exact 1 it is \(\text{logit}^{-1}(\eta - \psi_2)\), and in between the values follow a beta with mean \(\text{logit}^{-1}(\eta)\). The expected value combines the three parts:
ordbeta_mean <- function(m, newdata) {
eta <- predict(m, newdata = newdata, type = "link", re.form = NA)
psi <- m$fit$par[names(m$fit$par) == "psi"] # the two thresholds
p0 <- 1 - plogis(eta - psi[1]) # P(y = 0)
p1 <- plogis(eta - psi[2]) # P(y = 1)
p1 + (1 - p0 - p1) * plogis(eta)
}re.form = NA sets random effects to zero, so the function also works for mixed models. If the data contain no 1s, the upper threshold cannot be estimated well and gets a very large value (here about 19), which does no harm.
ordbeta_mean(m_ord, newd)[1] 0.01231587 0.01755623 0.02459877Both models give similar predictions; the hurdle beta additionally shows where the effect acts.
TipSolution
d2 <- data.frame(cover = varespec$Vaccmyrt / 100,
Hum = scale(varechem$Humdepth)[, 1]) # scaled, as Al above
table(d2$cover == 0)
FALSE TRUE
10 14 m_ord2 <- glmmTMB(cover ~ Hum, family = ordbeta(), data = d2)
summary(m_ord2) Family: ordbeta ( logit )
Formula: cover ~ Hum
Data: d2
AIC BIC logLik -2*log(L) df.resid
-7.1 -1.2 8.5 -17.1 19
Dispersion parameter for ordbeta family (): 11.7
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.4648 0.4744 -7.303 2.81e-13 ***
Hum 0.6604 0.2528 2.612 0.009 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(simulateResiduals(m_ord2))
m_hur2 <- glmmTMB(cover ~ Hum, ziformula = ~ Hum, family = beta_family, data = d2)
summary(m_hur2) Family: beta ( logit )
Formula: cover ~ Hum
Zero inflation: ~Hum
Data: d2
AIC BIC logLik -2*log(L) df.resid
-10.8 -4.9 10.4 -20.8 19
Dispersion parameter for beta family (): 11.5
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.2088 0.4896 -6.554 5.61e-11 ***
Hum 0.3675 0.3281 1.120 0.263
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Zero-inflation model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.4904 0.5173 0.948 0.343
Hum -1.7533 0.7710 -2.274 0.023 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(simulateResiduals(m_hur2))
The residuals of both models show no problems. The ordered beta finds a clear positive effect of humus depth. The hurdle beta shows where it acts: deeper humus lowers the probability that bilberry is absent (zero part, p = 0.02), but has no clear effect on its cover where it is present. The predicted probability of presence and the expected cover at shallow, mean and deep humus (-1, 0, 1 SD):
newd2 <- data.frame(Hum = c(-1, 0, 1))
1 - predict(m_hur2, newdata = newd2, type = "zprob") # probability of presence[1] 0.09589356 0.37979132 0.77951464predict(m_hur2, newdata = newd2, type = "response") # expected cover, hurdle[1] 0.002609822 0.014748874 0.042975591ordbeta_mean(m_ord2, newd2) # expected cover, ordered beta[1] 0.004196112 0.012423085 0.032719285Both models predict that the expected cover increases with humus depth (about tenfold, but from very low values; it stays below 5%), but only the hurdle beta tells us that this is mainly because bilberry occurs more often on deep humus.
9.7.3 Binomial GLM - snails
The following dataset contains information about the occurrences of 3 freshwater snail species and their infection rate with schistosomiasis (freshwater snails are intermediate hosts). The data and the analysis is from Rabone et al. (2019).
Rabone, M. et al. (2019) “Freshwater snails of biomedical importance in the niger river valley: Evidence of temporal and spatial patterns in abundance, distribution and infection with schistosoma spp.” Parasites & vectors, 12(1), pp. 1–20.
library(EcoData)
?EcoData::snailsThe dataset has data on three snail species: Biomphalaria pfeifferi (BP), Bulinus forskalii (BF), and Bulinus truncatus (BT). For this task, we will focus only on Bulinus truncatus (BT).
We want to answer two questions:
- What are the variables that explain the occurrence of the Bulinus truncatus species (=Species distribution model)
- What are the variables that influence the infection of the snails with Schistosomiasis
TipHints
Prepare dataset before the analysis:
data = EcoData::snails
data$sTemp_Water = scale(data$Temp_Water)
data$spH = scale(data$pH)
data$swater_speed_ms = scale(data$water_speed_ms)
data$swater_depth = scale(data$water_depth)
data$sCond = scale(data$Cond)
data$swmo_prec = scale(data$wmo_prec)
data$syear = scale(data$year)
data$sLat = scale(data$Latitude)
data$sLon = scale(data$Longitude)
data$sTemp_Air = scale(data$Temp_Air)
# Remove NAs
vars = c("Temp_Water", "pH", "Latitude", "Longitude", "Cond", "seas_wmo", "wmo_prec", "water_speed_ms", "water_depth", "Temp_Air", "year", "duration", "locality", "site_irn", "coll_date")
data = data[complete.cases(data[, vars]), ]Factors which drive the occurrence of BT
Minimal model:
model1 = glm(bt_pres~ sTemp_Water, data = data, family = binomial)Add environmental variables such as Temperature, pH, water metrics (e.g. water speed etc). Which variables could be potential confounders (year? locality? Season?….)
Factors which affect the infection rate
We use a k/n binomial model:
model2 = glm(cbind(BT_pos_tot, BT_tot - BT_pos_tot )~ sTemp_Water , data = data[data$BT_tot > 0, ], family = binomial)Add environmental variables such as Temperature, pH, water metrics (e.g. water speed etc). Which variables could be potential confounders (year? locality? Season?….)
TipSolution
Species distribution model
Environmental variables (part of our hypothesis): Temp_water, Temp_Air, pH, Cond, swmo_prec, water_speed_ms, and water_depth.
Potential confounders: site_type, year, seas_wmo (season) (, and locality)
Our model:
model1 = glm(bt_pres~ site_type + sTemp_Water + spH +
sCond + swmo_prec + swater_speed_ms + duration +
sTemp_Air + seas_wmo + syear + swater_depth,
data = data, family = binomial)
summary(model1)
Call:
glm(formula = bt_pres ~ site_type + sTemp_Water + spH + sCond +
swmo_prec + swater_speed_ms + duration + sTemp_Air + seas_wmo +
syear + swater_depth, family = binomial, data = data)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.815224 0.194147 4.199 2.68e-05 ***
site_typecanal.3 -0.332009 0.121838 -2.725 0.00643 **
site_typepond -1.256373 0.198187 -6.339 2.31e-10 ***
site_typerice.p -2.549160 0.226259 -11.267 < 2e-16 ***
site_typeriver -0.937624 0.178915 -5.241 1.60e-07 ***
site_typerivulet -0.655550 0.214947 -3.050 0.00229 **
site_typespillway -2.948800 0.410426 -7.185 6.73e-13 ***
sTemp_Water -0.147118 0.089085 -1.651 0.09865 .
spH -0.058152 0.050195 -1.159 0.24665
sCond 0.226405 0.057109 3.964 7.36e-05 ***
swmo_prec -0.156915 0.059027 -2.658 0.00785 **
swater_speed_ms -0.129785 0.054930 -2.363 0.01814 *
duration -0.007267 0.006853 -1.060 0.28899
sTemp_Air -0.207746 0.075345 -2.757 0.00583 **
seas_wmowet -0.070866 0.134537 -0.527 0.59838
syear -0.487076 0.063083 -7.721 1.15e-14 ***
swater_depth -0.050282 0.051463 -0.977 0.32854
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 2862.7 on 2071 degrees of freedom
Residual deviance: 2469.6 on 2055 degrees of freedom
AIC: 2503.6
Number of Fisher Scoring iterations: 4Intercept corresponds to seas_wmo == dry and site_type == can.2
Coefficients of confounders (here year, site type and season) are not interpreted (Table 2 fallacy, Section 7.4.1). We only interpret the variables of our hypothesis: conductivity has a clear positive effect, precipitation, water speed and air temperature clear negative effects. Let’s look at the effects:
plot(allEffects(model1))
it appears that conductivity has the largest positive effect on occurrence rates. It is reasonable to assume that there are interactions between freshwater (site_type) and factors such as conductivity:
model1 = glm(bt_pres~ site_type * sCond -sCond + swater_speed_ms+spH + sTemp_Water + swmo_prec + duration +
sTemp_Air + seas_wmo + syear + swater_depth,
data = data, family = binomial)
summary(model1)
Call:
glm(formula = bt_pres ~ site_type * sCond - sCond + swater_speed_ms +
spH + sTemp_Water + swmo_prec + duration + sTemp_Air + seas_wmo +
syear + swater_depth, family = binomial, data = data)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.814583 0.196286 4.150 3.33e-05 ***
site_typecanal.3 -0.353643 0.122684 -2.883 0.003945 **
site_typepond -1.441902 0.216914 -6.647 2.98e-11 ***
site_typerice.p -2.729871 0.251280 -10.864 < 2e-16 ***
site_typeriver -1.129095 0.238930 -4.726 2.29e-06 ***
site_typerivulet -0.744506 0.290013 -2.567 0.010254 *
site_typespillway -2.691676 0.419567 -6.415 1.40e-10 ***
swater_speed_ms -0.136708 0.055190 -2.477 0.013248 *
spH -0.070772 0.050585 -1.399 0.161788
sTemp_Water -0.156144 0.089931 -1.736 0.082517 .
swmo_prec -0.156697 0.059477 -2.635 0.008423 **
duration -0.007841 0.006913 -1.134 0.256646
sTemp_Air -0.205400 0.075886 -2.707 0.006796 **
seas_wmowet -0.020764 0.137386 -0.151 0.879868
syear -0.472567 0.063643 -7.425 1.13e-13 ***
swater_depth -0.055160 0.052236 -1.056 0.290978
site_typecan.2:sCond 0.283216 0.115382 2.455 0.014105 *
site_typecanal.3:sCond -0.044265 0.116647 -0.379 0.704333
site_typepond:sCond 0.478142 0.132729 3.602 0.000315 ***
site_typerice.p:sCond 0.717542 0.205410 3.493 0.000477 ***
site_typeriver:sCond -0.195577 0.347940 -0.562 0.574049
site_typerivulet:sCond 0.039191 0.417782 0.094 0.925263
site_typespillway:sCond 0.025064 0.148041 0.169 0.865554
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 2862.7 on 2071 degrees of freedom
Residual deviance: 2449.3 on 2049 degrees of freedom
AIC: 2495.3
Number of Fisher Scoring iterations: 4“-sCond” allows us to test the slopes for the different site_types against 0 (no reference level)
Residual checks:
library(DHARMa)
sim = simulateResiduals(model1, plot = TRUE)
Check misfits by plotting predictors against residuals:
plotResiduals(sim, data$sCond)
plotResiduals(sim, data$sTemp_Air)
They look okay!
Bonus - Spatial autocorrelation:
A common problem with spatial data is spatial autocorrelation. We can use DHARMa to test if our residuals are spatially autocorrelated
res2 = recalculateResiduals(sim, group = c(data$site_irn))
groupLocations = aggregate(cbind(data$sLat, data$sLon ), list( data$site_irn), mean)
testSpatialAutocorrelation(res2, x = groupLocations$V1, y = groupLocations$V2)
DHARMa Moran's I test for distance-based autocorrelation
data: res2
observed = 0.224538, expected = -0.011364, sd = 0.067266, p-value =
0.0004532
alternative hypothesis: Distance-based autocorrelationThe test is significant, i.e. we should assume that the data is spatially autocorrelated. Addressing this issue, however, is covered in chapter 4C
TipSolution
Infection rate
Environmental variables (part of our hypothesis): Temp_water, Temp_Air, pH, Cond, swmo_prec, water_speed_ms, and water_depth.
Potential confounders: site_type, year, seas_wmo (season), duration, ( and locality)
Our model:
data_inf = data[data$BT_tot > 0 , ] # only observations
model2 = glm(cbind(BT_pos_tot, BT_tot - BT_pos_tot )~ site_type + sTemp_Water + spH +
sCond + log(swmo_prec+1) + swater_speed_ms + log(duration) +
sTemp_Air + seas_wmo + syear + swater_depth,
data = data_inf, family = binomial)
summary(model2)
Call:
glm(formula = cbind(BT_pos_tot, BT_tot - BT_pos_tot) ~ site_type +
sTemp_Water + spH + sCond + log(swmo_prec + 1) + swater_speed_ms +
log(duration) + sTemp_Air + seas_wmo + syear + swater_depth,
family = binomial, data = data_inf)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.84336 0.52079 -5.460 4.77e-08 ***
site_typecanal.3 0.06733 0.14076 0.478 0.632399
site_typepond 2.42003 0.20151 12.010 < 2e-16 ***
site_typerice.p 1.06829 0.39734 2.689 0.007175 **
site_typeriver 0.11099 0.46215 0.240 0.810207
site_typerivulet 1.16396 0.33047 3.522 0.000428 ***
site_typespillway 1.20405 1.02591 1.174 0.240542
sTemp_Water 0.23679 0.10178 2.327 0.019991 *
spH 0.03519 0.06160 0.571 0.567794
sCond -0.10822 0.05731 -1.888 0.059000 .
log(swmo_prec + 1) -1.14756 0.33413 -3.434 0.000594 ***
swater_speed_ms 0.09253 0.04145 2.232 0.025595 *
log(duration) -0.90117 0.17268 -5.219 1.80e-07 ***
sTemp_Air -0.29302 0.09338 -3.138 0.001702 **
seas_wmowet 0.17373 0.16907 1.028 0.304160
syear -0.11184 0.06612 -1.691 0.090746 .
swater_depth -0.16272 0.06577 -2.474 0.013353 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1532.4 on 1106 degrees of freedom
Residual deviance: 1264.9 on 1090 degrees of freedom
AIC: 1597.1
Number of Fisher Scoring iterations: 6Intercept corresponds to seas_wmo == dry and site_type == can.2
Again, coefficients of confounders (year, site type, season) are not interpreted (Table 2 fallacy, Section 7.4.1). Let’s look at the estimates of our hypothesis variables:
plot(allEffects(model2))
With partial residuals:
plot(allEffects(model2, partial.residuals = TRUE))
Reminder: The partial residuals don’t tell us anything because the variance of the residuals (remember binomial) depends on the probability of the binomial process AND on k (k/n binomial model). We need DHARMa to check the residuals:
library(DHARMa)
sim = simulateResiduals(model2, plot = TRUE)
Residuals look okayish for a k/n model.
Check misfits by plotting predictors against residuals:
plotResiduals(sim, log(data_inf$swmo_prec+1))
plotResiduals(sim, data_inf$sTemp_Water)
swmo_prec doesn’t look good but this is caused by the many zeros in the variable
9.7.4 Making spatial predictions with a logistic regression - elephant SDM
To demonstrate how to fit a species distribution model and extrapolate it in space, let’s have a look at the elephant dataset, which is contained in EcoData
library(EcoData)
?elephantThe object elephant contains two subdatasets
elephant$occurenceData contains presence / absence data as well as bioclim variables (environmental predictors) for the African elephant
elephant$predictionData data with environmental predictors for spatial predictions
The environmental data consists of 19 environmental variables, called bio1 through bio19, which are public and globally available bioclimatic variables (see https://www.worldclim.org/data/bioclim.html for a description of the variables). For example, bio1 is the mean annual temperature. No understanding of these variables is required for the task, the only difficulty is that many of them are highly correlated because they encode similar information (e.g. there are several temperature variables).
The goal of this exercise is to fit a logistic regression model based on the observed presence / absences, and then make new predictions of habitat suitability in space across Africa based on the fitted model. Thus, our workflow consists of two steps:
- building and optimizing the predictive model, and
- using the predictive model to make predictions on new data and visualizing the results.
Here an example of how you could do this
Build predictive model:
model = glm(Presence~bio1, data = elephant$occurenceData, family = binomial()) # minimum modelTo check the predictive power of the model, you can either use
AIC(model) [1] 3511.176library(pROC)
auc(elephant$occurenceData$Presence, predict(model, type = "response"))Area under the curve: 0.7053The AUC is a common measure of how well a model separates presences from absences: the probability that a randomly chosen presence gets a higher predicted value than a randomly chosen absence (0.5 = random, 1 = perfect). However, it doesn’t penalize for complexity, as the AIC. Better to look at it under cross-validation. As explained in the section on cross-validation (Section 8.4), we collect the out-of-sample predictions of all folds and calculate the AUC once. Because the loop uses update(), it works for any model you fit to these data:
dat = elephant$occurenceData
set.seed(123)
fold = sample(rep(1:5, length.out = nrow(dat))) # random assignment to 5 folds
predCV = numeric(nrow(dat))
for (k in 1:5) {
fit = update(model, data = dat[fold != k, ])
predCV[fold == k] = predict(fit, newdata = dat[fold == k, ], type = "response")
}
auc(dat$Presence, predCV)Setting levels: control = 0, case = 1Setting direction: controls < casesArea under the curve: 0.7048Currently, fit and validation AUC are nearly identical (0.71), because we have a very simple model. However, for more complex models, this could change.
Drop some of the highly correlated variables (don’t use all of them).
Use quadratic effects
Use
stepAICor thedredgefunction (Tip 8.2) to optimize the predictive model
Make new predictions
The data for making spatial predictions is in elephant$predictionData. This new dataset is not a data.frame but a raster object, which is a special data class for spatial data. You can plot one of the predictors in the following way.
library(sp)
library(raster)
plot(elephant$predictionData$bio1)
As our new_data object is not a typical data.frame, we are not using the standard predict function for a glm, which is ?predict.glm, but the predict function from the raster object (which internally transforms the new_data into a classical data.frame, pass then the data.frame to our model, and then transforms the output back to a raster object). Therefore, the syntax is slightly different to how we previously used predict().
predictions = predict(elephant$predictionData, model = model, type = "response")
head(as.data.frame(predictions)) layer
1 0.1195466
2 0.1148339
3 0.1102836
4 0.1080683
5 0.1037547
6 0.1016556The advantage of the raster object is that we can directly use it to create a map (the raster object has coordinates for each observation):
spplot(predictions, colorkey = list(space = "left") )