In earlier chapters we learned many different statistical tools:
This final chapter pulls those ideas together. The aim is not to memorize a list of tests. The aim is to move carefully from a biological question to a planned model, and then from a fitted model to a checked, interpreted, plotted, and reported result.
Core idea
Do not start with “Which test do I use?” Start with the biological question, sketch the expected figure, identify the response variable and data structure, then choose the model.
The chapter has two parts:
Statistics does not begin after data collection. In a good study, the statistical plan is part of the study design.
Before collecting data, we should already be able to say:
Here is my biological question.
Here is the figure I expect to make.
Here is the response variable.
Here is the model I plan to fit.
Here is roughly how many observations I need, and how many residual degrees of freedom I expect.
This does not mean the plan can never change. Real data can surprise us. But without a plan, we often discover too late that the design does not answer the question.
Biological question
|
v
Sketch expected figure
|
v
Define response variable
|
+-- Continuous -> LM / Gaussian model
+-- Binary / proportion -> binomial GLM
+-- Count -> Poisson/count GLM
+-- Multivariate -> ordination / multivariate methods
|
v
Check data structure
|
+-- independent observations -> ordinary LM / GLM
+-- grouped / repeated / blocked / nested -> mixed model / random effects
|
v
Define explanatory variables
|
+-- categorical -> group differences
+-- continuous -> slopes
+-- categorical + continuous -> ANCOVA-style model
+-- multiple predictors -> multiple model
+-- interactions -> one effect depends on another
|
v
Planned statistical model
|
v
Plan sample size and expected residual dfTraditional names like ANOVA, regression, and ANCOVA are useful, but they are not the deepest decision. They are mostly different fixed-effect structures inside the linear model.
The deeper questions are:
A useful planning sentence is:
We will test whether X affects (or is related with) Y using model Z, because Y is [response type], X is [predictor type], and the observations are [independent/grouped]. We plan to collect n = … observations, giving n minus the number of estimated parameters residual degrees of freedom.
This sentence is not just a reporting sentence. It is a design sentence. If you cannot complete it, the study plan (including question) is probably not clear enough yet.
Residual degrees of freedom are not magic numbers that appear later in R. They are already implied by the design and the model.
Biological question:
Does fertilizer treatment affect plant biomass?
| Planning step | Decision |
|---|---|
| Expected figure | Biomass on the y-axis, fertilizer treatment on the x-axis; probably a boxplot, dot plot, or means with confidence intervals |
| Response variable | Plant biomass, continuous |
| Data structure | Each plant is measured once and plants are independent |
| Explanatory variables | Fertilizer treatment, categorical: control, low fertilizer, high fertilizer |
| Planned model | Gaussian linear model with biomass as the response and fertilizer treatment as the explanatory variable |
| Traditional name | One-way ANOVA, also understood as a linear model with one categorical predictor |
Example model:
m1 <- lm(biomass ~ fertilizer, data = plants)Suppose we plan 20 plants per treatment and 3 treatments. That gives \(n = 60\) observations.
The model estimates 3 treatment means, so the expected residual df is approximately:
60 - 3 = 57Planning sentence:
We will test whether fertilizer treatment affects plant biomass using a Gaussian linear model, because biomass is continuous, fertilizer treatment is categorical, and each plant is measured independently. We plan to collect 60 observations, 20 per treatment, giving an expected residual df of about 57.
Biological question:
Does temperature affect growth rate?
| Planning step | Decision |
|---|---|
| Expected figure | Growth rate on the y-axis, temperature on the x-axis; a scatterplot with a fitted line and uncertainty band |
| Response variable | Growth rate, continuous |
| Data structure | Each experimental unit is measured once and observations are independent |
| Explanatory variables | Temperature, continuous |
| Planned model | Gaussian linear model with a continuous predictor |
| Traditional name | Linear regression |
Example model:
m1 <- lm(growth_rate ~ temperature, data = growth)Suppose we plan \(n = 50\) observations across the temperature gradient.
The model estimates an intercept and a slope, so the expected residual df is:
50 - 2 = 48Planning sentence:
We will test whether temperature affects growth rate using a Gaussian linear model, because growth rate is continuous, temperature is continuous, and observations are independent. We plan to collect 50 observations, giving an expected residual df of about 48.
The expected figure helps clarify the model. If the figure is a scatterplot with one fitted line, the model is probably a regression-type model.
Biological question:
Does a drug treatment affect heart rate when each animal is measured repeatedly over time?
| Planning step | Decision |
|---|---|
| Expected figure | Heart rate on the y-axis, time or treatment on the x-axis, with repeated measurements connected or grouped by animal |
| Response variable | Heart rate, continuous |
| Data structure | Observations are not independent, because each animal contributes multiple measurements |
| Explanatory variables | Treatment is categorical; time may be continuous or categorical, depending on the design |
| Planned model | Mixed model, because repeated measurements are grouped within animals |
| Random effect | Animal identity, because measurements are repeated within animal |
Example model:
m1 <- lmer(heart_rate ~ treatment + time + (1 | animal_id), data = heart)If the biological question is whether the treatment changes the trajectory through time, then plan an interaction:
m2 <- lmer(heart_rate ~ treatment * time + (1 | animal_id), data = heart)Here we need to plan both the number of animals and the number of measurements per animal. For example, 20 animals measured 5 times each gives 100 observations, but those 100 observations are not equivalent to 100 independent animals.
Planning sentence:
We will test whether drug treatment affects heart rate over time using a mixed model, because heart rate is continuous, treatment and time are explanatory variables, and repeated measurements are grouped within animals. We plan to measure 20 animals 5 times each, giving 100 total observations, while accounting for animal identity as a random effect.
The random-effects question comes early. Before asking “ANOVA or regression?”, ask whether the observations are independent.
Biological question:
Does habitat type affect whether individuals are infected?
| Planning step | Decision |
|---|---|
| Expected figure | Infection probability on the y-axis and habitat type on the x-axis, perhaps shown as proportions with confidence intervals |
| Response variable | Infection status: infected or not infected, binary |
| Data structure | Each individual is sampled once and observations are independent |
| Explanatory variables | Habitat type, categorical |
| Planned model | Binomial GLM |
Example model:
m1 <- glm(infected ~ habitat, family = binomial, data = infection)Suppose we plan 50 individuals in each of two habitats, giving \(n = 100\) observations. With two habitat groups, the expected residual df is roughly:
100 - 2 = 98Planning sentence:
We will test whether infection probability differs between habitats using a binomial GLM, because infection status is binary, habitat is categorical, and individuals are independent. We plan to sample 100 individuals, 50 per habitat, giving an expected residual df of about 98.
The explanatory variable could be categorical or continuous, just as in a linear model. What changes the model class is the response variable.
Once we have the data, the key implementation question is:
Did the model behave well enough that we can interpret it?
The workflow is:
Fit planned model
|
v
Check model assumptions / diagnostics
|
+-- okay -> interpret model output
+-- not okay -> revise model / family / structure
|
v
Interpret output
|
v
Make explicit predictions
|
v
Plot data + model predictions
|
v
Report statistical and biological conclusionThe order matters. We do not fit a model, look for a small p-value, and only then check whether the model was reasonable. We check first, then interpret.
Good practice is to make predictions explicitly from the fitted model and then plot those predictions. The figure should show the model you fitted, not a model quietly fitted by the plotting function.
Biological question:
Does body size affect metabolic rate, and does this body-size effect differ among three species?
This example connects several ideas from the course:
| Planning step | Decision |
|---|---|
| Expected figure | Scatterplot of metabolic rate against body size, with separate fitted lines for each species |
| Response variable | metabolic_rate, continuous |
| Data structure | Independent observations, assuming each individual is measured once |
| Explanatory variables | body_size, continuous; species, categorical |
| Planned model | Gaussian linear model with a species-by-body-size interaction |
| Approximate residual df | Number of observations minus number of estimated parameters |
Planning sentence:
We will test whether body size affects metabolic rate differently among species using a Gaussian linear model with a species-by-body-size interaction, because metabolic rate is continuous, body size is continuous, species is categorical, and observations are independent.
Reading in the dataset:
metabolism <- read_csv("datasets/metabolism_species_body_size.csv")Rows: 135 Columns: 3
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): species
dbl (2): body_size, metabolic_rate
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Here are the first few rows of the dataset:
metabolism |>
slice_head(n = 6)# A tibble: 6 × 3
species body_size metabolic_rate
<chr> <dbl> <dbl>
1 species_A 9.92 3.71
2 species_A 27.2 8.35
3 species_A 14.2 5.58
4 species_A 13.7 8.06
5 species_A 18.4 6.82
6 species_A 18.4 5.48Planned figure:
ggplot(metabolism, aes(body_size, metabolic_rate, colour = species)) +
geom_point() +
labs(
x = "Body size",
y = "Metabolic rate",
colour = "Species"
) +
theme_classic()
The planned model has:
metabolic_ratebody_sizespeciesbody_size:speciesm_metabolism <- lm(metabolic_rate ~ species * body_size, data = metabolism)The model formula:
metabolic_rate ~ species * body_sizemeans:
metabolic_rate ~ species + body_size + species:body_sizeThe * means main effects plus the interaction. We use it when the biological question is whether the effect of one predictor depends on another predictor.
For a linear model, the common diagnostic questions are:
| Diagnostic plot | Question |
|---|---|
| Residuals vs fitted | Is there obvious nonlinearity or unequal variance? |
| QQ plot | Are residuals approximately normal? |
| Scale-location | Does residual spread change with fitted values? |
| Residuals vs leverage | Are a few observations dominating the model? |
autoplot(m_metabolism)Warning: `fortify(<lm>)` was deprecated in ggplot2 4.0.0.
ℹ Please use `broom::augment(<lm>)` instead.
ℹ The deprecated feature was likely used in the ggfortify package.
Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
ℹ Please use tidy evaluation idioms with `aes()`.
ℹ See also `vignette("ggplot2-in-packages")` for more information.
ℹ The deprecated feature was likely used in the ggfortify package.
Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
ℹ The deprecated feature was likely used in the ggfortify package.
Please report the issue at <https://github.com/sinhrks/ggfortify/issues>.
We are not looking for perfect plots. We are looking for obvious problems that would make the model inappropriate.
Possible actions if diagnostics look bad:
summary(m_metabolism)
Call:
lm(formula = metabolic_rate ~ species * body_size, data = metabolism)
Residuals:
Min 1Q Median 3Q Max
-4.1918 -0.8628 -0.0026 0.9103 3.5532
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.83601 0.56726 6.762 4.25e-10 ***
speciesspecies_B -0.42823 0.80217 -0.534 0.5944
speciesspecies_C -0.68723 0.79494 -0.865 0.3889
body_size 0.17335 0.03001 5.777 5.40e-08 ***
speciesspecies_B:body_size 0.23073 0.04327 5.333 4.20e-07 ***
speciesspecies_C:body_size -0.07204 0.04127 -1.746 0.0832 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.426 on 129 degrees of freedom
Multiple R-squared: 0.8028, Adjusted R-squared: 0.7952
F-statistic: 105 on 5 and 129 DF, p-value: < 2.2e-16anova(m_metabolism)Analysis of Variance Table
Response: metabolic_rate
Df Sum Sq Mean Sq F value Pr(>F)
species 2 632.80 316.40 155.492 < 2.2e-16 ***
body_size 1 324.39 324.39 159.417 < 2.2e-16 ***
species:body_size 2 111.55 55.78 27.411 1.201e-10 ***
Residuals 129 262.50 2.03
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Use this table to orient yourself:
| Question level | Model/output | Common test statistic | What it asks |
|---|---|---|---|
| Model / term level | ANOVA table for LM | F |
Does this explanatory variable or term explain variation in the response? |
| Model / term level | Deviance comparison for GLM | Deviance / chi-square-like test | Does adding/removing this term improve the model? |
| Coefficient level | Linear model coefficients | t |
Is this estimated coefficient different from zero? |
| Coefficient level | Binomial/Poisson GLM coefficients | z |
Is this estimated coefficient different from zero on the link scale? |
| Coefficient level | Confidence interval | CI | What range of effect sizes is compatible with the data/model? |
summary() often gives coefficient-level tests. anova() often gives term-level tests. These are related, but they are not always asking exactly the same question.
Residual degrees of freedom:
df.residual(m_metabolism)[1] 129Residual degrees of freedom are the amount of information left for estimating error after fitting the model. Here they are roughly the number of observations minus the number of estimated parameters.
nrow(metabolism)[1] 135length(coef(m_metabolism))[1] 6nrow(metabolism) - length(coef(m_metabolism))[1] 129Do not use geom_smooth() as the model.
Instead:
predict()new_metabolism <- expand_grid(
species = c("species_A", "species_B", "species_C"),
body_size = seq(5, 30, length.out = 100)
)
pred_metabolism <- predict(
m_metabolism,
newdata = new_metabolism,
interval = "confidence"
)
plot_metabolism <- bind_cols(
new_metabolism,
as_tibble(pred_metabolism)
)
plot_metabolism |>
slice_head(n = 6)# A tibble: 6 × 5
species body_size fit lwr upr
<chr> <dbl> <dbl> <dbl> <dbl>
1 species_A 5 4.70 3.85 5.56
2 species_A 5.25 4.75 3.91 5.59
3 species_A 5.51 4.79 3.96 5.62
4 species_A 5.76 4.83 4.02 5.65
5 species_A 6.01 4.88 4.08 5.68
6 species_A 6.26 4.92 4.13 5.71ggplot(metabolism, aes(body_size, metabolic_rate, colour = species)) +
geom_point(alpha = 0.65) +
geom_ribbon(
data = plot_metabolism,
aes(x = body_size, ymin = lwr, ymax = upr, fill = species),
alpha = 0.18,
colour = NA,
inherit.aes = FALSE
) +
geom_line(data = plot_metabolism, aes(y = fit), linewidth = 1) +
labs(
x = "Body size",
y = "Metabolic rate",
colour = "Species",
fill = "Species"
) +
theme_classic()
The model and the figure now match. The lines and ribbons come from m_metabolism, the model we fitted and checked.
Example structure:
We fitted a Gaussian linear model testing whether metabolic rate depended on body size, species, and their interaction. Diagnostic plots did not show major violations of model assumptions. Metabolic rate changed with body size, but the slope differed among species: the relationship was strongest in species B and weakest in species C. This indicates that species changed the relationship between body size and metabolism.
If giving formal results, add the relevant coefficient or term-level test from the output.
Biological question:
Does habitat type affect whether individuals are infected?
| Planning step | Decision |
|---|---|
| Expected figure | Infection probability or observed proportion infected by habitat |
| Response variable | infected, binary: 0 = not infected, 1 = infected |
| Data structure | Independent observations, assuming each individual is sampled once |
| Explanatory variables | habitat, categorical |
| Planned model | Binomial GLM |
| Approximate residual df | Number of observations minus number of estimated habitat probabilities |
Planning sentence:
We will test whether infection probability differs between habitats using a binomial GLM, because infection status is binary, habitat is categorical, and observations are independent.
Reading in the dataset:
infection <- read_csv("datasets/infection_habitat.csv")Rows: 160 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): habitat
dbl (1): infected
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Here are the first few rows of the dataset:
infection |>
slice_head(n = 6)# A tibble: 6 × 2
habitat infected
<chr> <dbl>
1 forest 0
2 forest 1
3 forest 0
4 forest 0
5 forest 0
6 forest 0Planned figure:
infection_summary <- infection |>
group_by(habitat) |>
summarise(
infected = mean(infected),
.groups = "drop"
)
ggplot(infection_summary, aes(habitat, infected)) +
geom_col(fill = "grey70", colour = "grey25") +
labs(
x = "Habitat",
y = "Observed proportion infected"
) +
theme_classic()
m_infection <- glm(
infected ~ habitat,
family = binomial,
data = infection
)For a binomial GLM, ask:
A simple overdispersion check:
deviance(m_infection) / df.residual(m_infection)[1] 1.238983A ratio much larger than 1 can suggest overdispersion. This is a warning sign that the model may be underestimating uncertainty.
summary(m_infection)
Call:
glm(formula = infected ~ habitat, family = binomial, data = infection)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.2368 0.2677 -4.619 3.85e-06 ***
habitatgrassland 1.3870 0.3492 3.972 7.14e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 212.70 on 159 degrees of freedom
Residual deviance: 195.76 on 158 degrees of freedom
AIC: 199.76
Number of Fisher Scoring iterations: 4In a binomial GLM, coefficients are on the log-odds scale by default. That scale is useful for the model, but not always intuitive biologically.
Predicted probabilities are often easier to interpret.
new_infection <- tibble(
habitat = c("forest", "grassland")
) |>
mutate(
predicted_probability = predict(
m_infection,
newdata = pick(everything()),
type = "response"
)
)
new_infection# A tibble: 2 × 2
habitat predicted_probability
<chr> <dbl>
1 forest 0.225
2 grassland 0.537ggplot(infection, aes(habitat, infected)) +
geom_jitter(width = 0.08, height = 0.03, alpha = 0.25) +
geom_point(
data = new_infection,
aes(y = predicted_probability),
size = 4,
colour = "firebrick"
) +
labs(
x = "Habitat",
y = "Infection status / predicted probability"
) +
theme_classic()
For GLMs, type = "response" gives predictions on the biological scale. For a binomial GLM, that means predicted probabilities.
Example structure:
We fitted a binomial GLM testing whether infection probability differed between habitats. The model used infection status as a binary response and habitat as a categorical explanatory variable. After checking model fit, we interpreted the model on the probability scale: predicted infection probability was higher in grassland than forest.
Biological question:
Does food availability affect the number of offspring produced?
| Planning step | Decision |
|---|---|
| Expected figure | Scatterplot of offspring number against food availability, with a fitted prediction line |
| Response variable | n_offspring, count |
| Data structure | Independent observations, assuming each individual or reproductive unit is measured once |
| Explanatory variables | food, continuous |
| Planned model | Poisson GLM, with overdispersion checked after fitting |
| Approximate residual df | Number of observations minus number of estimated parameters |
Planning sentence:
We will test whether food availability affects offspring number using a Poisson GLM, because offspring number is a count, food availability is continuous, and observations are independent.
Reading in the dataset:
offspring <- read_csv("datasets/offspring_food.csv")Rows: 90 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (2): food, n_offspring
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Here are the first few rows of the dataset:
offspring |>
slice_head(n = 6)# A tibble: 6 × 2
food n_offspring
<dbl> <dbl>
1 1.97 4
2 8.89 7
3 3.67 3
4 3.47 4
5 5.35 5
6 5.37 5Planned figure:
ggplot(offspring, aes(food, n_offspring)) +
geom_point(alpha = 0.7) +
labs(
x = "Food availability",
y = "Number of offspring"
) +
theme_classic()
m_offspring <- glm(
n_offspring ~ food,
family = poisson,
data = offspring
)deviance(m_offspring) / df.residual(m_offspring)[1] 1.338865For count models, overdispersion is a common issue. If the residual deviance is much larger than the residual df, the Poisson model may be too simple.
summary(m_offspring)
Call:
glm(formula = n_offspring ~ food, family = poisson, data = offspring)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.61383 0.12057 5.091 3.56e-07 ***
food 0.17327 0.01782 9.721 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 216.70 on 89 degrees of freedom
Residual deviance: 117.82 on 88 degrees of freedom
AIC: 401.27
Number of Fisher Scoring iterations: 5In a Poisson GLM, coefficients are on the log scale. A positive coefficient means the expected count increases as the predictor increases.
new_offspring <- tibble(
food = seq(0, 10, length.out = 100)
) |>
mutate(
predicted_count = predict(
m_offspring,
newdata = pick(everything()),
type = "response"
)
)ggplot(offspring, aes(food, n_offspring)) +
geom_point(alpha = 0.6) +
geom_line(
data = new_offspring,
aes(y = predicted_count),
colour = "firebrick",
linewidth = 1
) +
labs(
x = "Food availability",
y = "Number of offspring"
) +
theme_classic()
Example structure:
We fitted a Poisson GLM testing whether offspring number changed with food availability. Model checks did not suggest major overdispersion. Predicted offspring number increased with food availability, suggesting that food availability was positively associated with reproductive output.
Planning checklist:
Implementation checklist:
The goal is not to memorize every possible test. The goal is to move carefully from biological question to planned model, and then from fitted model to checked, interpreted, plotted, and reported result.