Chapter 22 Extending use cases of GLM
In the previous chapter, we used count data (bounded at 0 and integer-valued) to understand how generalised linear models work. This chapter extends our understanding by looking at another data type: binomial data.
22.1 Binomial response data
There are three common formats for representing binomial data, all of which involve the fundamental concepts of “success” and “failure.” It is up to you, the modeller, to define what constitutes a success or failure. These can be any data that can be grouped into two discrete categories, such as pass/fail, survived/died, presence/absence, yes/no, male/female, red/green.
The three formats provide flexibility for modelling binomial data in various contexts:
Binary Data (0/1): The data can be a two-level factor, often coded as zeros and ones. For example, this could represent whether an individual died (0) or survived (1) in a study, or two sexes (male or female) in a study predicting sex based on other variables (e.g., size, morphology).
Proportions: The data can be a numeric vector of proportions bounded between 0 and 1. For instance, the proportion of individuals in a group that survive some event can be represented this way. In such cases, the total number of trials contributing to the proportion can be included as
weightsin the regression model, ensuring that points representing larger sample sizes have more influence.Success/Failure Counts (Two-Column Matrix): The data can be expressed as a two-column matrix, where the first column gives the number of “successes” and the second gives the number of “failures.” For example, 25 successes and 15 failures in a study can be represented as
cbind(successes, failures). This format is commonly used when raw counts are available.
The aim of a generalized linear model (GLM) for binomial data is usually to estimate the probability of “success” (e.g., survival, passing, scoring a goal, presence). The predicted response (fit of the model) is a value between 0 and 1 on the natural scale, which can be interpreted as the probability of success.
To achieve this, binomial GLMs commonly use the logit link function, which transforms the expected probability (\(p\)) from the range (0, 1) to the range \((- \infty, \infty)\). On this “logit” (or log-odds) scale, the S-shaped relationship between a predictor and the probability becomes a straight line, so the model can be written as logit(p) = intercept + slope × x. As with the Poisson GLM, it is the expected value that is transformed, not the data: the 0/1 data themselves are never transformed (which is just as well, because the logit of 0 or 1 is minus or plus infinity). This type of regression is known as logistic regression.
The logit transformation is defined as:
\[ y = \log\left(\frac{p}{1-p}\right), \]
where \(p\) is the probability of success. Its inverse, the anti-logit, maps linear predictions (\(y\)) back to probabilities and is defined as:
\[ p = \frac{\exp(y)}{1 + \exp(y)}. \]
This transformation enables logistic regression to handle binary or proportional data effectively, with predictions that are interpretable as probabilities (e.g., of survival, success, or presence).
Why not simply transform the data and use lm? For binary (0/1) data this is impossible, because the logit of 0 or 1 is infinite. For proportions we could try it (as long as there were no 0s or 1s), but (like with Poisson regression) the other assumption of the ordinary linear model, homoscedasticity, would cause problems. With S-shaped binomial relationships, the expected variance is small at the two ends of the range and large in the middle. This contrasts strongly with the constant variance assumption of ordinary linear models. Therefore it is wise to account for this using a Generalised Linear Model that explicitly models this variance structure.
Let’s try a couple of examples.
In this chapter we use these packages (and, as always, work inside your course RStudio project):
22.2 Example: NFL field goals
In this example you will be dealing with binary data (0/1, failure/success) from the American National Football League (NFL). The data are a record of field goals, which are a relatively rare method of scoring where someone kicks the ball over the crossbar and through the uprights of the goal during play.
Our aim is to estimate how the probability of success changes with distance from the goal. We already have a good idea that success probability will decline with increasing distance! But at what rate does the probability decline? And at what distance is there a probability of 0.5 (i.e. 50% chance of success)?
First, import the data and convert the distance from yards to metres.
Next, plot the data with ggplot (geom_jitter would be a good option, but you might like to adjust the height and alpha arguments, e.g. geom_jitter(height = 0.05,alpha = 0.5)). You can see that the data are distributed in two bands with 1 representing success and 0 representing failure.
Now we can fit a GLM using an appropriate binomial error structure.
As is standard practice, we begin by examining the model diagnostics. In this case, they look quite problematic. But be careful: the main reason for the odd patterns in these plots is a limitation of the plots themselves. They were designed for linear models with normally distributed errors, and they give misleading results for logistic regression with binary (0/1) data, even when the model is perfectly appropriate. (Poor diagnostics can also be a sign that important variables are missing from a model, but as we will see below, that is not the explanation here.)
Given these limitations, it is often better to avoid relying heavily on traditional diagnostic plots for logistic regression. Instead, it is more prudent to evaluate your model choice based on theoretical principles, particularly those addressing the bounded nature of the response variable.
However, a newer diagnostic method specifically suited to models like logistic regression exists: the DHARMa package offers an alternative simulation-based approach that provides clearer and more reliable insights.
22.2.1 DHARMa
The DHARMa package is an alternative approach to model diagnostics, particularly where traditional tools like ggfortify::autoplot fall short, especially for models involving non-Gaussian data types such as binomial or Poisson distributions. These types of data often challenge standard diagnostic methods, making interpretation difficult.
At the heart of DHARMa’s approach is the simulateResiduals() function. The function simulates standardised residuals based on the model’s assumptions, including the specified error distribution and link function. DHARMa then constructs diagnostic plots that can reveal model deficiencies such as mismatched link functions, unmodelled predictors, and overdispersion.
The two plots are the Q-Q plot and the residuals versus predicted values plot. The Q-Q plot compares the observed residuals against an expected distribution generated from the simulations. A well-fitted model is reflected by points aligning along the diagonal. The residuals versus predicted values plot allows users to check for patterns in residuals, with a random scatter indicating a well-fitting model. Deviations from these expectations signal issues like incorrect link functions, overdispersion, or omitted variables.
Here is how to implement this approach:
The text on the Q-Q plot reports the results of three tests: the KS test, Dispersion test, and Outlier test, each targeting specific aspects of model fit. The KS (Kolmogorov-Smirnov) test checks whether the distribution of residuals aligns with the expected distribution, flagging deviations that suggest model misspecifications. The Dispersion test evaluates whether the variance in the residuals matches the model’s assumptions, identifying overdispersion or underdispersion. Finally, the Outlier test highlights residuals that fall outside expected ranges, pointing to potential anomalies or influential data points. Together, these tests provide a comprehensive diagnostic of model performance and in this case indicate that there are no problems.
As you can see, the two plots, and the in-built tests, indicate that the model structure is appropriate: The observed and expected points fall on the QQ diagonal line and the residuals versus predicted values plot looks like a random star-field without any patterns.
This example illustrates a general principle. DHARMa is a reliable, general-purpose diagnostic tool for GLMs. Because it works by simulation, it produces trustworthy residuals whatever the error family (binomial, Poisson, and so on) and whatever the structure of the data. The traditional autoplot diagnostics, by contrast, become less and less reliable as the response data become more discrete. They are perfectly readable for ordinary linear models and remain acceptable for GLMs with aggregated or high-count responses, but they can be very misleading for binary (0/1) data like the NFL field goals. The reason is that, for binary data, the residual at any given fitted value can only take one of two possible values, so the Q-Q plot looks alarming even when the model is perfectly appropriate, which is exactly what we saw with autoplot(nflMod) earlier.
Which diagnostic when?
-
autoplot(deviance residuals): fine for ordinary linear models, and acceptable for GLMs with aggregated or high-count data (e.g. grouped binomialcbinddata, or Poisson counts that are not dominated by zeros). -
DHARMa(simulated residuals): reliable for all GLMs. Essential for binary (0/1) data, whereautoplotcannot be trusted.
When in doubt, use DHARMa. If autoplot and DHARMa agree,
you can be confident; if they disagree, believe DHARMa.
22.2.2 Continuing the analysis
So, based on theoretical principles about the response variable’s bounds or the results of this DHARMa analysis, we can be pretty sure that the binomial family is the most appropriate family to use. Let’s continue.
We proceed in the normal way by obtaining the ANOVA table for the model. We need to specify that we want to calculate p-values using a “Chi” squared test. This shows us that indeed distance is an important variable in determining probability of success (not so surprising!)
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Success
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 947 955.38
## Distance 1 137.8 946 817.58 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We get more useful information from the coefficient summary of the model. This gives the intercept and slope of the model on the scale of the linear predictor (see the figure above).
##
## Call:
## glm(formula = Success ~ Distance, family = binomial, data = NFL)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.68958 0.45021 12.64 <2e-16 ***
## Distance -0.13118 0.01263 -10.38 <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: 955.38 on 947 degrees of freedom
## Residual deviance: 817.58 on 946 degrees of freedom
## AIC: 821.58
##
## Number of Fisher Scoring iterations: 5The slope (-0.131) is on the logit (log-odds) scale. Back-transforming it with exp() gives an odds ratio: the factor by which the odds of success are multiplied for each extra metre.
## [1] 0.8770628## 2.5 % 97.5 %
## 0.8556137 0.8990496We could report this something like this:
The probability of scoring a field goal declined with distance: each additional metre multiplied the odds of success by 0.877 (95% CI 0.856 to 0.899), i.e. reduced the odds by about 12% per metre (binomial GLM: \(\chi^2\) = 137.8, d.f. = 1, p < 0.001). The equation of the fitted line was logit(probability of success) = 5.690 -0.131 \(\times\) distance (m), and the probability of success fell to 50% at about 43.4 m (Figure 22.1).
The equation of the model if you want to express it on the natural scale works out to be:
\(y =\frac{1}{1+exp(-(\beta_0 + \beta_1 x))}\), or \(probability =\frac{1}{1+exp(-(5.690 - 0.131\times distance))}\).
Let’s make sure that works by creating a set of predicted data from this equation and plotting it onto the graph:
d1 <- data.frame(x = 10:50) %>%
mutate(y = 1 / (1 + exp(-(5.690 - 0.131 * x))))
A + geom_line(data = d1, aes(x, y), colour = "red")
Rather than using the equation, there is an easier way by using the R’s predict function. This was covered in the previous chapter, but for convenience, here are the steps:
create a data frame for the model what to predict
from.use the
predictfunction to predict fitted values and standard errors (SE).put these predicted values into a
newDatadata frame, multiplying the SE values these by 1.96 to get the 95% confidence intervals of the fitted values.back-transform the fitted values to the natural scale using the appropriate “inverse” of the link function.
So, this how to do it in this case. First we obtain the predictions and CI on the scale of the linear predictor (logit scale). We calculate 95% CI from the standard errors simply by multiplying by 1.96.
# Dataset to predict FROM
newData <- data.frame(Distance = 14:52)
# Get predictions from the model
predVals <- predict(nflMod, newdata = newData, se.fit = TRUE)
# Add those predictions to newDat
newData <- newData %>%
mutate(fit_LP = predVals$fit) %>%
mutate(lowerCI_LP = predVals$fit - 1.96 * predVals$se.fit) %>%
mutate(upperCI_LP = predVals$fit + 1.96 * predVals$se.fit)Now we can first obtain the inverse link from the model object family(nflMod)$linkinv, and use that to backtransform the data onto the natural scale ready for plotting.
# Get the inverse link function
inverseFunction <- family(nflMod)$linkinv
# transform predicted data to the natural scale
newData <- newData %>%
mutate(
fit = inverseFunction(fit_LP),
lowerCI = inverseFunction(lowerCI_LP),
upperCI = inverseFunction(upperCI_LP)
)Now we can finally plot the model predictions and the 95% confidence intervals for them.
# The plot and ribbon
(A <- ggplot(NFL, aes(x = Distance, y = Success)) +
geom_ribbon(
data = newData, aes(x = Distance, ymin = lowerCI, ymax = upperCI),
fill = "grey75", alpha = 0.5, inherit.aes = FALSE
) +
geom_line(data = newData, aes(x = Distance, y = fit)) +
geom_jitter(height = 0.05, alpha = 0.5))
Figure 22.1: Probability of scoring an NFL field goal in relation to distance, with the fitted binomial GLM and its 95% confidence interval.
So, at what distance does the probability of success fall to just 50%? You could read this directly from the plot as “approximately 44 m”. Alternatively, you could obtain the value from the newData dataset you created above with predictions from the model, by filtering it. This confirms the probability of success reaches 0.5 somewhere between 43 and 44 m.
## Distance fit
## 1 42 0.5449184
## 2 43 0.5122432
## 3 44 0.4794630The most precise way is to use the model equation. The probability is 0.5 when the odds are 1, i.e. when the logit is 0. Setting intercept + slope × distance = 0 gives distance = -intercept / slope:
## [1] 43.3734122.3 Example: Sex ratio in turtles
In this example we will look at sex ratio of hawksbill turtles (Eretmochelys imbricata)16. The data are counts of males and females in clutches of eggs incubated at different temperatures.
The sex ratio in the species varies with temperature during incubation. We are interested in the “tipping point” temperature, at which the sex ratio is 50:50. This is known as the pivotal temperature.
This is a small dataset, you can look at the whole thing:
## Temperature Total Nmale Nfemale
## 1 27.9 12 12 0
## 2 28.4 29 29 0
## 3 28.9 31 26 5
## 4 29.3 20 18 2
## 5 29.5 22 14 8
## 6 29.8 10 2 8
## 7 30.4 10 0 10We are interested in sex ratio which we can calculate as the proportion of the population that is female (i.e. number of females divided by total number).
Let’s plot that data. We can use the trick of telling R to plot the points different sizes depending on the sample size (using the size = argument):
You can see that the proportion of females increases with temperature. Let’s fit a model to these data to better understand them. We could use the propFemale data as the response variable, but there is a big problem with that: we would be giving equal weight to each of the data points, even though the sample size for each one ranges from 10 to 31. This is not good because we would have much more faith in the very large sample sizes than the small ones.
Instead we can bind the data into a two column matrix using cbind. The first column is our “success” and the second column is our “failure”. If we put females in the first column the model will predict “probability of being female”, which is what we want. This two-column approach provides the model with the sample size which it can use to weight the regression appropriately.
Now lets fit the model.
As ever, we should take a quick look at the model diagnostic plots. Note the important difference from the NFL example: these turtle data are aggregated binomial data (each row summarises many individual turtles via cbind), rather than binary 0/1 data. As explained in the “Which diagnostic when?” box above, this is exactly the situation where the traditional autoplot deviance-residual plots are trustworthy, and here they look fine.
To confirm this, we can also check the model with DHARMa.
The three formal tests all pass: the KS test, the dispersion test, and the outlier test are each non-significant, so DHARMa agrees with autoplot that the binomial model is appropriate. The non-significant dispersion test is especially reassuring: it suggests that the data are not seriously overdispersed. (The rule-of-thumb index, residual deviance divided by its degrees of freedom, is 1.7, which falls in the “caution” band. But with only 5 residual degrees of freedom that index is very unreliable, which is a good reason to rely on the simulation-based test here.) You will, however, notice that the right-hand panel highlights a “quantile deviation” in red. Don’t panic: the combined adjusted quantile test (reported at the top of that panel) is itself non-significant, and with such a tiny dataset (only 7 rows) DHARMa’s quantile lines are very easily wiggled by one or two points, so an isolated red flag like this should not be over-interpreted. The contrast with the NFL example is the instructive part: there, autoplot and DHARMa genuinely disagreed and DHARMa was the one to trust; here both tools agree the model is fine, and the lone red flag is a small-sample artefact rather than a real problem.
Now we can look at the anova table. This will tell us what we already know - there is a strong effect of temperature on sex ratio.
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: y
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 6 70.353
## Temperature 1 61.869 5 8.484 3.671e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Now for the coefficients:
##
## Call:
## glm(formula = y ~ Temperature, family = binomial, data = hawksbill)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -111.7200 22.4260 -4.982 6.30e-07 ***
## Temperature 3.7754 0.7625 4.951 7.37e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 70.3528 on 6 degrees of freedom
## Residual deviance: 8.4841 on 5 degrees of freedom
## AIC: 24.184
##
## Number of Fisher Scoring iterations: 5This table gives us the coefficients for the formula of our relationship. We could use these to produce a formula of the form \(y =\frac{1}{1+exp(-(\beta_0 + \beta_1 x))}\).
It is perhaps more useful to plot the model fit onto the plot. First we need to create a data frame to predict from:
Now we can predict the values (and 95% CI) on the scale of the linear predictor (logit).
pv <- predict(modA, newDat, se.fit = TRUE)
newDat <- newDat %>%
mutate(
propFemale_LP = pv$fit,
lowerCI_LP = pv$fit - 1.96 * pv$se.fit,
upperCI_LP = pv$fit + 1.96 * pv$se.fit
)Now we can use the inverse link function to backtransform to the natural probability scale.
# Get the inverse link function
inverseFunction <- family(modA)$linkinv
# transform predicted data to the natural scale
newDat <- newDat %>%
mutate(
propFemale = inverseFunction(propFemale_LP),
ymin = inverseFunction(lowerCI_LP),
ymax = inverseFunction(upperCI_LP)
)We can add these values to the plot like this:
# The plot and ribbon
(A <- ggplot(hawksbill, aes(
x = Temperature, y = propFemale,
size = Total
)) +
geom_ribbon(
data = newDat, aes(x = Temperature, ymin = ymin, ymax = ymax),
fill = "grey75", alpha = 0.5, inherit.aes = FALSE
) +
geom_smooth(
data = newDat, aes(x = Temperature, y = propFemale),
stat = "identity", inherit.aes = FALSE
) +
geom_point())
Can you read off the graph the approximate temperature at which the sex ratio is 50:50 (the pivotal temperature)? As with the NFL example, you can calculate it exactly from the model as -intercept / slope:
## [1] 29.59139
Try refitting the model using simply the proportion female
(propFemale) data rather than the two-column
(cbind) approach. Then try writing up the result in the
same way as shown for the NFL field goals example.
22.4 Example: Smoking
As I mentioned above, binomial regression can be applied to anything where the data can be classified into two groups. I’ll illustrate that now with an example about smoking.
The data set is very small and looks like this:
| Student smokes | Student does not smoke | |
|---|---|---|
| Parent(s) smoke | 816 | 3203 |
| No parents smoke | 188 | 1168 |
The dataset is the number of students smoking and not smoking grouped according to whether their parents smoke. These data are binomial/proportion data because the values in the cells of the table are constrained by the overall total and students must fall into one of the categories. We can use these data to calculate the probability of the student being a smoker.
Before fitting a GLM lets just work out these probabilities by hand. What is the probability that a child of smoking parents is a smoker themselves? This is simply 816/(816+3203) = 0.2030. (i.e. it is the number of smokers divided by the total number). Similarly, the probability that the child of non-smokers smokes is 188/(188+1168) = 0.1386.
But is this a significant difference? That is what we are trying to find out using a GLM.
To do this, we can turn these data into a two column matrix of success (yes - smoker) and failure (no - non-smoker).
## 1_yes 0_no
## [1,] 816 3203
## [2,] 188 1168So we can see “success” on the left and “failure” on the right”.
We now create a (tiny) data.frame for the parental status (smoker = “1_yes”, non-smoker = “0_no”).
Now we can fit the model. Pause now and think about what the NULL hypothesis is here. It is that parental smoking does not have any influence on whether the child smokes, and that the probability that the student smokes is unrelated to parental smoking.
Diagnostic plots will not tell us anything useful for this model. It has two parameters (an intercept and the effect of parental smoking) estimated from only two rows of data, so it fits the data perfectly: there are no residual degrees of freedom left.
As usual, we first ask for the (Analysis of Deviance table using anova.
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: y
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 1 29.121
## parentSmoke 1 29.121 0 0.000 6.801e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1This gives clear evidence that the probability of a student smoking differs according to whether their parents smoke.
We can find out what this pattern is by examining the summary table.
##
## Call:
## glm(formula = y ~ parentSmoke, family = binomial, data = smoke)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.82661 0.07858 -23.244 < 2e-16 ***
## parentSmoke1_yes 0.45918 0.08782 5.228 1.71e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2.9121e+01 on 1 degrees of freedom
## Residual deviance: 1.7186e-13 on 0 degrees of freedom
## AIC: 19.242
##
## Number of Fisher Scoring iterations: 2This shows us the estimates on the logit scale. We can use predict to get a sense for these predictions on the more intuitive probability scale.
First we calculate the fitted values and 95% confidence intervals on the scale of the linear predictor (logit):
pv <- predict(smokeMod, smoke, se.fit = TRUE)
smoke <- smoke %>%
mutate(
prob_LP = pv$fit,
lowerCI_LP = pv$fit - 1.96 * pv$se.fit,
upperCI_LP = pv$fit + 1.96 * pv$se.fit
)Then we can backtransform these values to the probability scale using the inverse link function:
# Get the inverse link function
inverseFunction <- family(smokeMod)$linkinv
# transform predicted data to the natural scale
smoke <- smoke %>%
mutate(
prob = inverseFunction(prob_LP),
ymin = inverseFunction(lowerCI_LP),
ymax = inverseFunction(upperCI_LP)
)
smoke## parentSmoke prob_LP lowerCI_LP upperCI_LP prob ymin ymax
## 1 1_yes -1.367429 -1.444287 -1.290570 0.2030356 0.1908823 0.2157563
## 2 0_no -1.826606 -1.980629 -1.672583 0.1386431 0.1212518 0.1580801This table shows us that the probability of the students smoking is 0.2030 (95% CI = 0.191-0.216) for children of smokers, and 0.139 (95% CI = 0.121-0.158) for children of non-smokers. Another common way of expressing the difference is as an odds ratio, exp() of the parentSmoke1_yes coefficient: the odds of smoking are 1.58 times higher for children of smokers. We can plot the probabilities using ggplot like this.
ggplot(smoke, aes(
x = parentSmoke, y = prob, ymin = ymin,
ymax = ymax
)) +
geom_point() +
geom_segment(aes(xend = parentSmoke, y = ymin, yend = ymax))
22.5 Overdispersion in binomial models
Just as count data modelled with Poisson errors can suffer from overdispersion, so can binomial data. Recall from the previous chapter that overdispersion occurs when the response varies more than the error family assumes, and that it is often caused by an important variable being left out of the model, or by non-independent (clustered) observations. Exactly the same problem, with the same causes, can affect a binomial GLM.
Two points are worth emphasising for the binomial case:
- It only applies to aggregated binomial data, that is, proportion data or a two-column
cbind(successes, failures)response, like the turtle and smoking examples. Overdispersion cannot be diagnosed from binary (0/1) data such as the NFL field goals: with only a single trial per row the residual deviance carries no information about dispersion, so the usual checks are meaningless. Overdispersion is only a concern when each row summarises several trials. - How to detect it. The methods are the ones you already know. You can use the DHARMa dispersion test (as we did for the turtles, where it was comfortably non-significant), or the “rule of thumb” of dividing the residual deviance by its degrees of freedom and checking whether the result is much greater than 1 (values above about 1.5 call for caution, and values above about 2 are a clear warning sign).
What to do when it is detected. The fix mirrors the Poisson case. Refit the model exactly as before, but replace family = binomial with family = quasibinomial. The quasibinomial family estimates a dispersion parameter and widens the standard errors (and therefore the p-values) accordingly, while leaving the coefficient estimates themselves unchanged. Because a dispersion parameter is now being estimated, you should assess significance with an F test rather than the chi-squared test (i.e. anova(model, test = "F")). The negative-binomial and zero-inflation options mentioned in the previous chapter are specific to count data, but the quasi-likelihood idea carries over directly to binomial data through the quasibinomial family. The net effect, as with quasi-Poisson, is a more honest and more cautious picture of how certain you can be about the relationship.
Try it: GLM Family Lab: binomial models (open the lab)
Choose the binomial family.
- With 1 trial per observation (0/1 data), why does the lab give no verdict on the dispersion ratio?
- Use the odds ratio and the 50% point to describe the fitted curve in words.
- Increase the trials to 20 and add extra variation. What happens to the dispersion ratio and to the binomial confidence interval for the slope?
22.6 Exercise: Forensic footprints
A crime has been committed. At the scene, investigators recover a single clear footprint, and from it a forensic technician estimates that the person who left it had a foot length of 271 mm. Was the culprit male or female?
To answer this, you have access to a reference database, morphometry.csv, containing measurements from 155 people of known sex. For each person the data record their Sex (Male or Female) and three body measurements in millimetres: Height, HandLength and FootLength. Your job is to build a “forensic tool” - a binomial GLM (logistic regression) that predicts sex from foot length - evaluate how accurate it is, and then use it to make a prediction about the person who left the print.
This is the kind of binary (Male vs. Female) single-predictor problem you saw in the NFL field goals example, so use that worked example as your template. Unlike the examples above, this one is not worked through for you - figure out the steps yourself, and ask for help if you get stuck.
- Import the reference data. Your response variable needs to be binary (0/1), so create a new column that is
1when the person is male and0when they are female (ifelseis one way to do this). You will therefore be modelling the probability of being male. - Plot the data (
geom_jitterworks well here - remember to adjust theheightandalphaarguments as in the NFL example). - Fit an appropriate binomial GLM predicting your 0/1 variable from
FootLength. - Check the model using
DHARMa(simulateResiduals). Does the binomial family look appropriate? - Get the analysis of deviance table (
anovawithtest = "Chi") and thesummarytable. Is foot length a significant predictor of sex? (Remember the coefficients are on the logit scale.) - How good is your forensic tool? Use the model to predict the probability of being male for everyone in the reference database (
predictwithtype = "response"), classify each person as male if that probability is greater than 0.5, and compare your classification with their true sex. What proportion of people does your tool classify correctly? - The crime scene. Use
predictto estimate the probability (and its 95% confidence interval) that a person with a foot length of 271 mm is male. What is your conclusion, and how confident are you? Would this be strong enough evidence to present in court? - A second, fainter print at the scene suggests a foot length of about 249 mm. Predict the probability of being male for this print too. Why is it much harder to draw a conclusion here? (Think about the foot length at which the probability of being male is 0.5.)
data from: Godfrey et al. (1999) Can. J. Zool. 77: 1465–1473↩︎