Workshop 10: Bernoulli GLMs

BI3010 Statistics for Biologists

Author

Roos & Pinard

Published

9 August 2026

Glossary / 词汇表

Probabilities and odds

Bernoulli distribution 伯努利分布

A probability distribution for binary outcomes: observations can only take the value 0 (failure) or 1 (success). The distribution has a single parameter p, the probability of success. It is a special case of the Binomial distribution with one trial. In R, a Bernoulli GLM is fitted using family = binomial(link = "logit"), which handles 0/1 responses automatically.

用于二元结果的概率分布：观测值只能取0（失败）或1（成功）。该分布有一个参数p，即成功的概率。它是试验次数为1的二项分布的特例。在R中，伯努利GLM通过family = binomial(link = "logit")拟合，可自动处理0/1响应变量。

Odds 几率

The ratio of the probability of success to the probability of failure: p / (1 − p). If the probability of success is 0.75, the odds are 0.75 / 0.25 = 3, meaning success is three times as likely as failure. Odds range from 0 (impossible) to +∞ (certain), unlike probabilities which range from 0 to 1.

成功概率与失败概率之比：p / (1 − p)。若成功概率为0.75，则几率为0.75 / 0.25 = 3，表示成功是失败可能性的三倍。与概率（范围0到1）不同，几率的范围从0（不可能）到+∞（必然）。

Log-odds 对数几率

The natural logarithm of the odds: log(p / (1 − p)). Also called the logit of p. Log-odds can range from −∞ (probability of 0) to +∞ (probability of 1), making them suitable as the linear predictor in a Bernoulli GLM. A log-odds of 0 corresponds to a probability of 0.5; positive values correspond to probabilities above 0.5; negative values to probabilities below 0.5.

几率的自然对数：log(p / (1 − p))，也称为p的logit。对数几率的范围从−∞（概率为0）到+∞（概率为1），因此适合作为伯努利GLM的线性预测子。对数几率为0对应概率0.5；正值对应高于0.5的概率；负值对应低于0.5的概率。

The logit link

Logit link logit连接

The link function used in a Bernoulli GLM. It transforms the probability p to the log-odds scale: logit(p) = log(p / (1 − p)). This allows the linear predictor to range from −∞ to +∞ while keeping predicted probabilities between 0 and 1. The logit link transforms the linear predictor, not the observed 0/1 data.

伯努利GLM中使用的连接函数。它将概率p变换到对数几率尺度：logit(p) = log(p / (1 − p))。这使得线性预测子可以在−∞到+∞范围内取值，同时将预测概率保持在0和1之间。logit连接变换的是线性预测子，而非观测到的0/1数据。

Back-transformation (logit) 反变换（logit）

Converting a prediction from the log-odds scale back to a probability. The inverse of the logit function is the logistic function, applied in R with plogis(). For a log-odds value of x, the corresponding probability is plogis(x) = exp(x) / (1 + exp(x)). Confidence interval limits must also be back-transformed after being calculated on the log-odds scale.

将预测值从对数几率尺度转换回概率。logit函数的逆函数是logistic函数，在R中通过plogis()应用。对于对数几率值x，对应的概率为plogis(x) = exp(x) / (1 + exp(x))。置信区间上下限也必须在对数几率尺度上计算后再进行反变换。

Diagnostics

Dispersion (Bernoulli) 离散度（伯努利）

Unlike the Poisson distribution, the Bernoulli distribution has no free dispersion parameter. The variance is entirely determined by the mean: var = p(1 − p). There is therefore no overdispersion to check for in a Bernoulli GLM: the residual deviance / degrees of freedom check used in Workshop 9 does not apply here.

与泊松分布不同，伯努利分布没有自由的离散参数。方差完全由均值决定：var = p(1 − p)。因此，伯努利GLM中无需检查过度离散，第9次研讨课中使用的残差偏差/自由度检验在此不适用。

plogis() plogis()函数

The R function that converts log-odds to a probability. plogis(0) = 0.5, plogis(2) ≈ 0.88, plogis(-2) ≈ 0.12. It is the inverse of qlogis() (which converts a probability to log-odds). Use it whenever you need to back-transform predictions or confidence interval limits from a Bernoulli GLM.

将对数几率转换为概率的R函数。plogis(0) = 0.5，plogis(2) ≈ 0.88，plogis(-2) ≈ 0.12。它是qlogis()（将概率转换为对数几率）的逆函数。每当需要将伯努利GLM的预测值或置信区间上下限进行反变换时，使用此函数。

check_model() / performance check_model() / performance 包

A function from the performance package (first met in Workshop 4) that produces model-appropriate diagnostic plots. For a Bernoulli GLM it automatically switches to diagnostics suited to binary data, including a Posterior Predictive Check and a Binned Residuals plot, and prints guidance text on each panel describing what a passing plot looks like.

来自 performance 包（第4次研讨课首次介绍）的函数，可生成与模型类型相匹配的诊断图。对于伯努利GLM，它会自动切换到适合二元数据的诊断图，包括后验预测检验和分箱残差图，并在每张图上打印说明文字，描述通过的图应有的样子。

Binned residuals 分箱残差

A residual plot designed for binary responses. Instead of plotting all the individual 0/1 residuals (which fall into two uninformative bands), observations are sorted by their fitted probability, grouped into bins, and the average residual of each bin is plotted against the average fitted value. Averaging cancels the 0/1 noise, so a systematic trend becomes visible. Most points should fall within the shaded error bounds and scatter around zero with no slope. This is the panel to focus on when checking a Bernoulli GLM.

专为二元响应设计的残差图。它不绘制所有单个的0/1残差（那样会形成两条毫无信息的带状），而是将观测值按其拟合概率排序、分组到各个区间（箱），并将每个箱的平均残差对平均拟合值作图。取平均消除了0/1噪声，从而使系统性趋势变得可见。大多数点应落在阴影误差范围内，围绕零散布且没有斜率。这是检查伯努利GLM时应重点关注的面板。

Interactions

Interaction 交互作用

When the effect of one predictor depends on the value of another. Fitted in R with * (e.g. genotype * stratification), which adds the two main effects plus their interaction term. On the log-odds scale the interaction coefficient is the difference between the two groups in how much the second predictor changes the log-odds. Here, cold stratification raises germination far more in wild-type seeds than in the dormancy mutant: the effect of cold depends on genotype.

当一个预测变量的效应取决于另一个预测变量的取值时，即存在交互作用。在R中用 * 拟合（例如 genotype * stratification），它会加入两个主效应以及它们的交互项。在对数几率尺度上，交互系数表示第二个预测变量改变对数几率的幅度在两组之间的差异。此处，冷层积对野生型种子萌发的提升远大于对休眠突变体的提升：冷处理的效应取决于基因型。

Gene × environment (G×E) 基因型 × 环境（G×E）

A specific kind of interaction in which the effect of an environmental treatment differs between genotypes (or vice versa). A genotype that responds strongly to a treatment while another barely responds is the signature of G×E. In this workshop, wild-type germination is highly sensitive to cold stratification whereas the dog1 dormancy mutant germinates readily regardless of temperature.

一种特定的交互作用，即环境处理的效应在不同基因型之间存在差异（反之亦然）。一个基因型对处理反应强烈而另一个几乎无反应，正是 G×E 的标志。在本次研讨课中，野生型的萌发对冷层积高度敏感，而 dog1 休眠突变体无论温度如何都能顺利萌发。

Learning Objectives

Learning objectives

By the end of this workshop you should be able to:

Explain why binary response data cannot be modelled with a Normal linear model and why a Bernoulli GLM is appropriate.
Describe how the logit link function connects probabilities to log-odds, and convert between the two using plogis().
Fit a Bernoulli GLM in R using glm() with family = binomial(link = "logit").
Diagnose a Bernoulli GLM using both the base residual plots and the binned residuals plot from check_model(), and explain why overdispersion is not a concern.
Interpret model coefficients on the log-odds scale, including an interaction between two categorical predictors alongside a separate continuous effect, and calculate predicted probabilities for specific combinations of predictor values.

Workshop structure

The Analytical Workflow

As always, we work through the same seven-step workflow:

Formulate a research question.
Perform exploratory data analysis (EDA).
Identify any hidden assumptions.
Fit an appropriate model [Focus of today].
Diagnose the model and check assumptions [Focus of today].
Summarise the results [Focus of today].
Interpret and draw conclusions.

In Workshop 9 we went from Normal linear models to Poisson GLMs for count data. Today we take one more step. The response variable is not a count but a binary outcome: each observation is either a 0 or a 1. You can’t model that with a Normal distribution without all kinds of ugly consequences, which we’ll get to shortly. The appropriate distribution is the Bernoulli, the appropriate link is the logit, and the whole thing is more commonly known as logistic regression.

1. Formulate a research question

The dataset for this workshop comes from a seed-germination experiment in the model plant Arabidopsis thaliana. Freshly harvested seeds are often dormant: they will not germinate even in good conditions until dormancy is broken, usually by a spell of cold (winter). The gene DOG1 (DELAY OF GERMINATION 1) is a key regulator of this dormancy. The experiment compares wild-type seeds (WT) with a dog1 loss-of-function mutant, each given either a warm or a cold pre-treatment, to ask how genotype and cold interact to control germination. Download the data below and save it to your BI3010/data/ folder.

Download germination.txt

Each row is one seed. The variables are:

germinated: whether the seed germinated within the test window (1 = yes, 0 = no).
seed_mass: the mass of the seed in micrograms (µg), a continuous measure of seed provisioning; heavier seeds tend to germinate more readily.
genotype: whether the seed was wild-type ("WT") or the dormancy mutant ("dog1").
stratification: whether the seed received a warm pre-treatment ("Warm") or a cold, dormancy-breaking pre-treatment ("Cold").

The research question is predictive: how do seed mass, genotype, and cold stratification combine to determine the probability that a seed germinates, and does the effect of cold differ between the two genotypes?

Load the data:

germination <- read.table("Data/germination.txt", header = TRUE)
germination$genotype <- factor(germination$genotype,
                               levels = c("WT", "dog1"))
germination$stratification <- factor(germination$stratification,
                                     levels = c("Warm", "Cold"))

The last two lines set the factor levels explicitly so that "WT" and "Warm" are the reference (baseline) categories. This means the intercept will represent a warm-treated wild-type seed, the genotype coefficient will represent the difference for the dog1 mutant, and the stratification coefficient will represent the difference for cold treatment, each relative to that baseline.

2. Perform exploratory data analysis

The response variable is binary, so some of the usual EDA tools behave a bit differently. geom_histogram() on a 0/1 variable produces a two-bar chart (not very illuminating, but fine for checking the split between germinated and not). table() gives counts; mean() gives the proportion that germinated.

str(germination)
summary(germination)
table(germination$germinated)
mean(germination$germinated)   # overall germination rate

table(germination$genotype, germination$germinated)
table(germination$stratification, germination$germinated)

A raw scatter plot of 0s and 1s is deeply unhelpful: every point sits on one of two horizontal lines and you get no sense of how the data are distributed. Adding a small amount of vertical jitter helps enormously:

ggplot(germination, aes(x = seed_mass, y = germinated)) +
  geom_jitter(height = 0.05, alpha = 0.4) +
  labs(x = "Seed mass (µg)", y = "Germinated (0/1)") +
  theme_minimal()

Because the research question is really about how cold interacts with genotype, the most useful EDA plot shows the germination rate for each genotype under each treatment. If the effect of cold is larger for one genotype than the other, the two lines won’t be parallel; that non-parallelism is the visual signature of an interaction:

ggplot(germination, aes(x = stratification, y = germinated,
                        colour = genotype, group = genotype)) +
  geom_point(stat = "summary", fun = "mean", size = 3) +
  geom_line(stat = "summary", fun = "mean") +
  labs(x = "Stratification", y = "Proportion germinated", colour = "Genotype") +
  theme_minimal()

3. Identify hidden assumptions

The usual hidden assumptions apply: validity (does germinated capture what we mean, a seed that produces a viable seedling rather than just a burst seed coat?), representativeness (do these lab-grown seeds generalise to field conditions?), and independence (each seed must be a separate, independent observation; seeds from the same maternal plant or the same Petri dish may be more alike than independent seeds, which would violate this).

Dispersion isn’t a concern for Bernoulli GLMs. The variance of a Bernoulli outcome is p(1 − p), which is completely determined by p, the probability of germination. There’s no free dispersion parameter to estimate or check.

4. Fit an appropriate model

Why not use a linear model?

The response variable germinated can only be 0 or 1. A Normal linear model fitted to binary data has two problems. First, it will happily predict values below 0 and above 1, and you can’t have a probability of −0.2 or 1.3. Second, the spread of 0/1 residuals is largest when p ≈ 0.5 and shrinks near 0 or 1, so the equal variance of errors assumption is violated by design. The linear model will run without complaint. It’s just wrong.

Question 1. If you fitted lm(germinated ~ seed_mass, data = germination), which of the following would be a problem?

The linear model cannot handle a continuous predictor when the response is binary. The linear model would refuse to run and return an error. The linear model can predict values below 0 and above 1, which are not valid probabilities, and the equal variance of errors assumption is violated because residuals from a binary response are not uniformly spread. The linear model assumes the response is a count, which conflicts with a binary response.

The linear model will run without error and produce coefficient estimates, but its predictions for seed mass at the extremes of the range will fall outside [0, 1]. For example, a predicted germination probability of 1.3 or −0.2 has no meaning. The model also assumes equal spread of errors at all fitted values, but for a binary response the spread depends on the probability itself: it’s largest when p ≈ 0.5 and smallest near p = 0 or p = 1. Both problems disappear when we use a GLM with the logit link.

The Bernoulli distribution and logit link

The Bernoulli distribution models a single binary trial: each observation is either 0 or 1, with probability p of being 1. The model for the distributional assumption is:

\[germinated_i \sim Bernoulli(p_i)\]

where $p_i$ is the probability that seed $i$ germinates.

The logit link connects the linear predictor to $p_i$ by transforming probabilities to log-odds:

\[logit(p_i) = \log\!\left(\frac{p_i}{1-p_i}\right) = \beta_0 + \beta_1 \times seedmass_i\]

Log-odds can range from $-\infty$ to $+\infty$, so they work as a linear predictor without risk of predictions outside [0, 1]. To get predicted probabilities from the log-odds, we apply plogis().

Question 2. A log-odds of 0 corresponds to which probability of success?

0 (impossible) 1 (certain) 0.5 (equally likely to succeed or fail) The probability cannot be determined from the log-odds alone.

Log-odds of 0 means the odds are exp(0) = 1, which means success and failure are equally likely: p = 0.5. You can verify this with plogis(0) in R. More generally: negative log-odds → probability below 0.5, positive log-odds → probability above 0.5.

Fitting the model

In R, a Bernoulli GLM uses family = binomial(link = "logit"). The family name is binomial rather than bernoulli because Bernoulli is actually a special case of the Binomial (one trial per observation), and R covers both with the same family. As always, fitting the model is the easiest part. Start with seed_mass as the only predictor:

germ_mod <- glm(germinated ~ seed_mass,
                family = binomial(link = "logit"),
                data = germination)

5. Diagnose the model and check assumptions

par(mfrow = c(2, 2))
plot(germ_mod)
par(mfrow = c(1, 1))

The residual plots for a Bernoulli GLM look different from those of a linear model or Poisson GLM. Because each observation is either 0 or 1, the raw residuals cluster into two distinct bands in the residuals vs fitted plot: one band for seeds that germinated (germinated = 1) and one for seeds that didn’t (germinated = 0). This two-banded appearance is expected and doesn’t indicate a problem.

Question 3. You notice the residuals vs fitted plot shows two distinct bands of points rather than a single cloud. What does this mean?

The model has failed to converge and the estimates are unreliable. There are two separate subgroups in the data that should be modelled separately. This is expected for a binary response. The two bands are the residuals from seeds that germinated (1) and those that did not (0). It does not indicate a violation of assumptions. The equal variance of errors assumption is violated and the model should be refitted with a transformation.

Think of it this way: every observation is exactly 0 or 1. There’s nowhere else for the points to go. The two bands are the residuals from the germinated seeds (1s) and the non-germinated seeds (0s), not a sign of trouble, just an unavoidable consequence of a binary response. What you’re looking for is a systematic trend within each band. A curve or funnel inside the bands is a concern. Two flat parallel bands with no internal pattern is fine.

This two-banded plot is honestly not much use for spotting real problems, which is exactly why we’ll turn to a better diagnostic (binned residuals) below.

Question 4. Should you use the normal Q-Q plot to assess this Bernoulli GLM?

Yes: the Q-Q plot shows whether residuals are normally distributed, which is an assumption of all GLMs. Yes, but only if the sample size is small. No: a Bernoulli GLM does not assume normally distributed residuals. The Q-Q plot is checking an assumption that was never made, so it tells us nothing useful about model fit. No: Q-Q plots can only be used for continuous response variables.

As with the Poisson GLM in Workshop 9, the Q-Q plot is largely uninformative here. Normality of residuals is an assumption of the Normal linear model, not of GLMs. A Bernoulli GLM assumes residuals follow the Bernoulli distribution, and the Q-Q plot doesn’t check that. Focus your diagnostic attention on the residuals vs fitted and scale-location plots, and on the residuals vs leverage plot for influential observations.

Question 5. Do you need to check for overdispersion in this Bernoulli GLM?

Yes: calculate residual deviance divided by residual degrees of freedom, just as in Workshop 9. Yes: binary data is even more prone to overdispersion than count data. No: the Bernoulli distribution has no free dispersion parameter. The variance is fixed as p(1 − p) once p is estimated, so there is nothing to check. No: dispersion only matters when the response variable has more than two possible values.

Overdispersion is a concern for Poisson GLMs but not for Bernoulli. Once you know p, the variance is exactly p(1 − p): there’s no free dispersion parameter to worry about. Honestly, after the doomsday dispersion values we saw in Workshop 9, it’s a relief. The residual deviance / degrees of freedom check doesn’t apply here.

Better diagnostics with `performance`

The two-banded base plot is honest but hard to read: with every point pinned to 0 or 1, a real problem is easy to miss. Back in Workshop 4 we met the performance package and its check_model() function for prettier, better-labelled diagnostics. It’s even more useful here, because for a Bernoulli GLM it automatically swaps in diagnostics designed for binary data:

install.packages(c("performance", "see"))  # Run once
library(performance)
check_model(germ_mod)

Two panels matter most:

Binned residuals. This is the fix for the two-banded plot. Instead of plotting all 200 individual 0/1 residuals, check_model() sorts the seeds by their fitted probability, groups them into bins, and plots the average residual of each bin against the average fitted value. Averaging cancels the 0/1 noise, so any systematic trend the raw plot hid becomes visible. The grey band is an approximate error bound: if the model fits, most points fall inside it, scattered around zero with no slope. A run of points drifting above or below zero, or a clear curve, would signal that the relationship with the predictor isn’t quite right (for example, non-linearity on the log-odds scale).
Posterior predictive check. check_model() uses the fitted model to simulate many fake datasets and overlays them on the observed data. If the model is reasonable, the simulated proportion of 1s and 0s should look like the real split.

As in Workshop 4, each panel prints guidance text describing what a passing plot looks like. Read it in context: not every assumption is equally important, and minor wobbles in real data are normal.

Question 6. On the binned residuals plot for germ_mod, almost all the points fall inside the grey error bounds and scatter around zero with no obvious slope. What should you conclude?

The two-banded pattern has been hidden, so this plot is less trustworthy than the base residuals vs fitted plot. The model fits acceptably: binning has removed the 0/1 noise, and with points inside the bounds and no systematic trend there is no sign of misfit with seed mass. The seeds should be re-binned until every point lies exactly on the zero line. Points inside the bounds indicate overdispersion, which must now be corrected.

Points scattered inside the error bounds with no slope is exactly what a well-fitting Bernoulli GLM looks like. Binning is what makes this readable: it averages away the 0/1 noise that produced the two useless bands, leaving a plot in which a genuine problem (a curve, or a run of points drifting off zero) would actually show up. Here nothing does, so seed mass enters the model on a sensible scale. And overdispersion still doesn’t apply; points inside the bounds are a good sign, not a dispersion diagnostic.

6. Summarise the results

summary(germ_mod)

The coefficient table gives estimates on the log-odds scale. From the output, the fitted equation is approximately:

\[logit(\hat{p}_i) = -1.63 + 0.089 \times seedmass_i\]

Predictions on the log-odds scale

Substitute a value of seed_mass directly into the equation to get the predicted log-odds.

Question 7. What is the predicted log-odds of germination for a seed of mass 20 µg? And for a seed of mass 30 µg?

-1.63 + 0.089 * 20   # =  0.15  (seed_mass = 20)
-1.63 + 0.089 * 30   # =  1.04  (seed_mass = 30)

A log-odds of +0.15 is just above zero, so the probability of germination is just above 0.5. A log-odds of +1.04 is well above zero, so the probability is clearly above 0.5. To get exact probabilities, back-transform with plogis().

Predictions on the probability scale

To convert log-odds to a probability, use plogis(). You can also work it out from the definition: p = exp(log-odds) / (1 + exp(log-odds)).

Question 8. Using plogis(), convert your log-odds predictions from Question 7 to probabilities. At what seed mass does the model predict a 50% chance of germination?

plogis(0.15)   # ≈ 0.54  (seed_mass = 20: 54% chance of germination)
plogis(1.04)   # ≈ 0.74  (seed_mass = 30: 74% chance of germination)

For a 50% chance of germination, we need the log-odds to equal 0:

-1.63 + 0.089 * seed_mass = 0
seed_mass = 1.63 / 0.089 ≈ 18.3

The model predicts that a seed mass of around 18.3 µg gives a 50% probability of germination. Below that the odds favour no germination; above it they favour germination.

Plotting the probability curve

Create a fake dataset, generate predictions on the log-odds scale, then back-transform with plogis() to get predicted probabilities and confidence interval limits:

fake <- data.frame(seed_mass = seq(from = min(germination$seed_mass),
                                   to   = max(germination$seed_mass),
                                   length.out = 50))

preds       <- predict(germ_mod, newdata = fake, se.fit = TRUE)
fake$fit    <- plogis(preds$fit)
fake$low    <- plogis(preds$fit - 1.96 * preds$se.fit)
fake$upp    <- plogis(preds$fit + 1.96 * preds$se.fit)

Question 9. Using the fake dataset, create a figure showing predicted probability of germination against seed mass, with a 95% confidence band. Add the raw data points using geom_jitter().

ggplot() +
  geom_ribbon(data = fake,
              aes(x = seed_mass, ymin = low, ymax = upp),
              fill = "grey80") +
  geom_line(data = fake,
            aes(x = seed_mass, y = fit)) +
  geom_jitter(data = germination,
              aes(x = seed_mass, y = germinated),
              height = 0.03, alpha = 0.35) +
  labs(x = "Seed mass (µg)",
       y = "Probability of germination") +
  theme_minimal()

The result is a gently S-shaped (logistic) curve: probability rises across the mass range while staying within [0, 1] at both ends. The S-shape is a direct consequence of the logit link. Because seed mass is only a modest predictor on its own, the rise here is gradual rather than steep; most of the action, as we’ll see next, is in genotype and stratification.

7. Interpret and draw conclusions

The model predicts that germination probability increases with seed mass. The positive slope means that each additional microgram of seed mass increases the log-odds of germination by approximately 0.089. On the probability scale the relationship isn’t linear, but seed mass alone is a fairly weak predictor: it nudges the probability rather than driving it. The real drivers are the ones we’ve ignored so far, genotype and cold stratification, so let’s add them.

Adding genotype and stratification

The research question isn’t really about seed mass; it’s about how genotype and cold stratification combine to control germination, and specifically whether cold matters more for one genotype than the other. That “whether the effect of one predictor depends on another” is exactly what an interaction captures. We fit it with *, which adds both main effects and their interaction, while keeping seed_mass as a separate additive term:

germ_mod2 <- glm(germinated ~ seed_mass + genotype * stratification,
                 family = binomial(link = "logit"),
                 data = germination)

Writing genotype * stratification is shorthand for genotype + stratification + genotype:stratification. The seed_mass term sits outside the interaction: it has one slope that applies to every seed, which is what we mean by a “separate continuous effect.”

Diagnose the new model

par(mfrow = c(2, 2))
plot(germ_mod2)
par(mfrow = c(1, 1))

library(performance)
check_model(germ_mod2)

The diagnostics should look similar to the first model: two-banded base residuals are still expected, the Q-Q plot is still uninformative, and the binned residuals should sit inside their bounds. One new panel appears in check_model() output: a Collinearity check. Because the model now contains an interaction, you may see high collinearity flagged for the terms involved in it. This is expected and not a problem: interaction terms are built from their main effects, so they’re correlated with them by construction. It’s only worth worrying about when two separate predictors are highly collinear.

Summarise and interpret

summary(germ_mod2)

From the output, the fitted equation on the log-odds scale is approximately:

\[logit(\hat{p}_i) = -4.43 + 0.11 \times seedmass_i + 2.71 \times dog1_i + 4.08 \times Cold_i - 2.46 \times (dog1_i \times Cold_i)\]

where dog1 and Cold are 1 when the seed is the mutant genotype and cold-treated respectively, and 0 otherwise.

Question 10. The stratificationCold coefficient is about +4.08 and the interaction genotypedog1:stratificationCold is about −2.46. What does the interaction term tell you?

Cold stratification reduces germination in dog1 seeds, because the interaction coefficient is negative. The boost from cold is much smaller for dog1 than for WT. Cold adds +4.08 to the log-odds of germination in WT but only 4.08 − 2.46 = +1.62 in dog1; in odds terms cold multiplies the odds by ~59 in WT but only ~5 in dog1. dog1 seeds are 2.46 times less likely to germinate than WT seeds overall. Seed mass has a weaker effect on germination in dog1 seeds than in WT seeds.

The interaction coefficient modifies the effect of cold depending on genotype. For WT (the reference), cold adds +4.08 to the log-odds. For dog1, cold adds 4.08 + (−2.46) = +1.62. Exponentiating: cold multiplies the odds of germination by exp(4.08) ≈ 59 in WT but only exp(1.62) ≈ 5 in dog1. Note cold is still positive for dog1; it just does far less. That makes biological sense: the dog1 mutant has lost dormancy, so its seeds germinate readily even when warm and cold has little left to add. This is a gene × environment interaction: the effect of the environment (cold) depends on the genotype.

Question 11. Using the fitted equation from germ_mod2, calculate the predicted probability of germination at seed mass = 25 µg for all four genotype × treatment combinations.

Start from the base term for a warm WT seed at 25 µg: −4.43 + 0.11 × 25 = −1.68. Then add the relevant coefficients:

plogis(-1.68)                  # WT,   Warm ≈ 0.16
plogis(-1.68 + 4.08)           # WT,   Cold ≈ 0.92
plogis(-1.68 + 2.71)           # dog1, Warm ≈ 0.74
plogis(-1.68 + 2.71 + 4.08 - 2.46)  # dog1, Cold ≈ 0.93

Cold transforms germination in WT (0.16 → 0.92) but barely moves dog1 (0.74 → 0.93). That gap between the two genotypes’ response to cold is the interaction, now expressed as probabilities. Notice these match the raw cell germination rates you saw in the EDA, a good sign the model is capturing the data.

Plot the model fit

Let’s use ggeffects and have it do the tedious coding for us. With an interaction, pass all three predictors to ggpredict(): the first goes on the x-axis, the second becomes coloured curves, and the third splits into panels:

library(ggeffects)

preds2 <- ggpredict(germ_mod2, terms = c("seed_mass", "genotype", "stratification"))
plot(preds2)

You’ll get two panels, one per stratification treatment. In the Warm panel the WT curve sits far below the dog1 curve (dormant wild-type vs non-dormant mutant). In the Cold panel both curves are high and close together. That change, a big genotype gap under warm conditions that nearly closes under cold, is the visual signature of the interaction.

A new challenger: amphibian chytrid infection

If time allows, work through the following dataset. It comes from a field study of chytridiomycosis, a disease caused by the fungus Batrachochytrium dendrobatidis (Bd) that has driven amphibian declines worldwide. Researchers swabbed frogs at a set of ponds and tested each for Bd infection, recording the water temperature at capture and whether the pond was permanent or temporary (temporary ponds dry out seasonally, which interrupts the fungus’s life cycle). Download the data below and save it to your BI3010/data/ folder.

Download chytrid.txt

The dataset contains 170 swabbed frogs. Each row is one frog, with the following variables:

infected: whether the frog tested positive for Bd (1 = yes, 0 = no).
water_temp: water temperature at the capture site, in °C (continuous, 12 to 24). Bd grows best in cool water and is inhibited when it warms up.
pond_type: whether the frog was caught at a "Permanent" pond or a "Temporary" one.

The research question is predictive: how are water temperature and pond type associated with the probability that a frog is infected with Bd?

chytrid <- read.table("Data/chytrid.txt", header = TRUE)
chytrid$pond_type <- factor(chytrid$pond_type, levels = c("Temporary", "Permanent"))

This dataset has the simpler additive shape (one continuous predictor plus one two-level factor, no interaction), so a model like infected ~ water_temp + pond_type is a sensible place to start. Go through the full seven-step workflow. A few things to keep in mind:

The response is binary, so use jitter plots rather than raw scatter plots to visualise it. mean(chytrid$infected) gives the overall infection rate; tapply(chytrid$infected, chytrid$pond_type, mean) gives it separately by pond type.
Diagnose the model with both plot() and check_model(). The binned residuals plot is again the one to focus on; there is no interaction here, so don’t expect a collinearity warning.
Think carefully about hidden assumptions. A single swab is an imperfect test: what does a negative result actually guarantee? And is water temperature at the moment of capture a good stand-in for the conditions the fungus experienced?
Use plogis() to back-transform any predictions from the log-odds scale to the probability scale.

Question 12. After completing your chytrid analysis, summarise your main findings. How do water temperature and pond type affect the probability of infection?

Warmer water is associated with a lower probability of infection: the water_temp coefficient is negative on the log-odds scale (about −0.23), consistent with Bd being inhibited at higher temperatures. Frogs at permanent ponds are more likely to be infected than those at temporary ponds: the pond_typePermanent coefficient is positive (about +1.57), which exponentiates to roughly a 4.8-fold increase in the odds of infection. In probability terms, the model predicts a 50% infection risk at around 14 °C at temporary ponds but not until about 21 °C at permanent ponds; permanent water shifts the whole curve towards higher infection because the fungus is never flushed out by the pond drying.

Key hidden assumptions worth flagging: a single swab can miss a genuine infection (a negative doesn’t guarantee the frog is Bd-free), and a one-off temperature reading may not represent the conditions the fungus actually experienced over the frog’s recent history. Both sit in the validity column of our assumption table.

What you can take away from this course

The skills you’ve learnt on this course aren’t actually “biology” skills. They’re data skills, and data skills go anywhere. Finance, healthcare, tech, government, sport, NGOs; any organisation sitting on data and trying to make sense of it needs people who can do what you can now do. That’s what “data analyst”, “quantitative researcher”, and “data scientist” actually mean in most job ads.

Here’s what to put on your CV:

Statistical methods

Linear regression
Interaction modelling
Generalised Linear Models: Poisson regression for count data, logistic regression for binary outcomes
Model diagnostics and assumption checking
Prediction and uncertainty quantification

Software

R and RStudio
ggplot2 for data visualisation

One-liner for applications:

Quantitative data analysis in R, including linear and generalised linear models (Poisson and logistic regression), model diagnostics, and data visualisation.

Logistic regression is worth mentioning specifically. It’s everywhere outside academia: customer churn, clinical outcomes, loan defaults, election forecasting. You can now run them.

You don’t need every detail memorised. You can always come back to this material for the details. You’ve fitted the models, read the diagnostics, and interpreted results in context.

Fin

That’s a wrap. This is the last workshop of the course, so you’re free to go.

Two final requests from me (Deon):

I know from past years that the majority of you will forget the details of the course and methods a few months after the course. That’s fine. What I do want you to try and hold onto is this: data and models (statistical, machine learning or AI) aren’t magic boxes that unravel the Truth or reveal the secrets of the universe. A model is just a bunch of code and equations. Don’t be taken in by people throwing data and results at you, whether that’s on the internet, in the news, or in peer-reviewed papers. Ask what assumptions were made. Ask whether the data can actually answer the question being claimed. Be sceptical.
If you found any typos, confusing sections, or parts of the workshops or lectures that just didn’t land: please tell me. I genuinely think this material matters, and if I’m not communicating it well, I want to know. Feel free to email or message however you prefer. I’d really appreciate any feedback (just please don’t be needlessly cruel about it).

Assumption table

Assumption	Description	Example of violation
1. Validity	The data can answer the research question.	Scoring a burst seed coat as “germinated” without checking that a viable seedling actually emerged.
2. Representativeness	The sample reflects the population of interest.	Testing only large, well-filled seeds because they were easiest to handle, then generalising to all seeds.
3. Additivity	Effects are additive on the log-odds scale unless an interaction is fitted. Here we deliberately relax this: `genotype * stratification` lets the effect of cold differ between genotypes.	Leaving out the interaction when cold genuinely helps WT far more than `dog1`, forcing both genotypes to share one cold effect.
4. Linearity	The relationship between each continuous predictor (seed mass) and the log-odds of germination is linear. Checked with the binned residuals plot.	Binned residuals curving away from zero, showing seed mass acts non-linearly on the log-odds scale.
5. Independence of errors	Each seed is an independent observation.	Seeds from the same maternal plant or the same Petri dish being more alike than independent seeds.
6. Equal variance of errors	Not assumed for Bernoulli GLMs. The variance is p(1 − p), fixed once p is known.	Not applicable: heteroscedasticity is built into the Bernoulli distribution and is not a violation.
7. Normality of errors	Not assumed for Bernoulli GLMs.	Not applicable: the Q-Q plot is uninformative and should not be used to assess model fit.
8. Dispersion	Not applicable. The Bernoulli distribution has no free dispersion parameter to check.	Unlike Poisson GLMs, overdispersion cannot occur with a true Bernoulli response.