ChickWeightlibrary(dplyr) # wrangling verbs
library(tidyr) # pivot_wider, drop_na
library(ggplot2) # plots
library(forcats) # factor handling
library(labelled) # set_value_labels, to_factor, is.labelled
library(broom) # tidy(), glance() -- turns test output into data frames
broom is worth calling out. Every test we run today
returns a messy list object that prints nicely but is painful to work
with. tidy() converts any of them into a one-row data frame
with consistent column names (estimate,
statistic, p.value, conf.low,
conf.high). That makes results easy to store, compare side
by side, or drop into a table.
ChickWeight records the body weight of 50 newly hatched
chicks, each weighed repeatedly from day 0 through day 21, under one of
four experimental diets. It is a small, clean, genuinely experimental
dataset, and it happens to contain every structure we need this
week:
| What we need | Where it lives in ChickWeight |
|---|---|
| A continuous outcome | weight (grams) |
| Two continuous variables to correlate | weight at two different days |
| Repeated measures on the same unit | each chick weighed 12 times |
| A grouping variable with 2 levels | any two diets |
| A grouping variable with 4 levels | Diet |
| Real, non-random missingness | 5 chicks stop appearing before day 21 |
One structural note before we start. Each chick contributes 12 rows, not one. The 578 rows in this file are not 578 independent observations — they are 50 chicks measured repeatedly. Almost every mistake you can make with this dataset comes from forgetting that.
glimpse(ChickWeight)
## Rows: 578
## Columns: 4
## $ weight <dbl> 42, 51, 59, 64, 76, 93, 106, 125, 149, 171, 199, 205, 40, 49, 5…
## $ Time <dbl> 0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 21, 0, 2, 4, 6, 8, 10, 1…
## $ Chick <ord> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
## $ Diet <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
We will follow the same wrangling spine used in class: rename into readable names, attach value labels to numeric codes, convert labelled variables to factors, deal with missingness explicitly, then summarise.
chicks <- ChickWeight %>%
as_tibble() %>%
select(chick_id = Chick, day = Time, weight_g = weight, diet = Diet) %>%
mutate(
chick_id = as.character(chick_id),
day = as.numeric(day),
diet = as.numeric(as.character(diet)) # back to numeric codes 1-4
) %>%
set_value_labels(
diet = c("Diet 1" = 1, "Diet 2" = 2, "Diet 3" = 3, "Diet 4" = 4)
) %>%
mutate_if(is.labelled, to_factor)
glimpse(chicks)
## Rows: 578
## Columns: 4
## $ chick_id <chr> "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "…
## $ day <dbl> 0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 21, 0, 2, 4, 6, 8, 10,…
## $ weight_g <dbl> 42, 51, 59, 64, 76, 93, 106, 125, 149, 171, 199, 205, 40, 49,…
## $ diet <fct> Diet 1, Diet 1, Diet 1, Diet 1, Diet 1, Diet 1, Diet 1, Diet …
levels(chicks$diet)
## [1] "Diet 1" "Diet 2" "Diet 3" "Diet 4"
Two things to notice.
First, diet arrives as a factor whose levels are the
bare characters "1", "2", "3",
"4". That is exactly the situation you hit with survey data
coded 1 = Male, 2 = Female: the numbers are labels, not
quantities. We push it back to numeric, attach real labels with
set_value_labels(), and convert with
to_factor(). It is a round trip on this dataset, but it is
the same round trip you will do on every coded variable you meet, and
doing it here means the printed output says Diet 3 instead
of 3.
Second, as.numeric(as.character(diet)) is deliberate.
as.numeric() applied directly to a factor returns the
internal level codes, not the labels. On this dataset
the codes and labels happen to coincide, so nothing breaks — but on a
factor with levels c("10", "20", "30"),
as.numeric() returns 1, 2, 3. Always go through
as.character() first.
There are no NA values in ChickWeight.
There is still missing data.
Five chicks stop being weighed before day 21. They do not appear as
NA rows — their rows simply do not exist. This is
implicit missingness, and it is invisible until you
reshape the data so that every chick is guaranteed a slot for every
day.
any(is.na(x)) outputs false.
chicks_wide <- chicks %>%
filter(day %in% c(0, 16, 21)) %>%
pivot_wider(
id_cols = c(chick_id, diet),
names_from = day,
values_from = weight_g,
names_prefix = "day_"
)
nrow(chicks_wide) # 50 chicks, one row each
## [1] 50
chicks_wide %>% summarise(across(starts_with("day_"), ~ sum(is.na(.))))
## # A tibble: 1 × 3
## day_0 day_16 day_21
## <int> <int> <int>
## 1 0 3 5
Now the gaps are NA and countable: every chick has a day
0 weight, 3 are missing at day 16, and 5 are missing at day 21.
chicks_wide %>%
filter(is.na(day_21)) %>%
count(diet, name = "n_dropped")
## # A tibble: 2 × 2
## diet n_dropped
## <fct> <int>
## 1 Diet 1 4
## 2 Diet 4 1
Four of the five dropouts are on Diet 1 and one is on Diet 4. That is not random. Whatever caused those chicks to leave the study is related to their diet, and plausibly to their weight as well. When we compare day-21 weights across diets, we are comparing chicks that survived to day 21 — which is not quite the same population we started with. This does not invalidate the analysis, but it is the first thing to say out loud in the limitations section, and it is a much more common real-world problem than anything involving p-values.
# 45 chicks with both a day-16 and a day-21 weight (for the paired analysis)
chicks_paired <- chicks_wide %>%
drop_na(day_16, day_21) %>%
mutate(gain = day_21 - day_16)
# 45 chicks measured at day 21 (for the group comparisons)
chicks_d21 <- chicks %>%
filter(day == 21) %>%
mutate(diet = fct_relevel(diet, "Diet 1")) # Diet 1 = reference group
nrow(chicks_paired)
## [1] 45
nrow(chicks_d21)
## [1] 45
chicks_d21 %>%
group_by(diet) %>%
summarise(
n = n(),
mean_wt = round(mean(weight_g), 1),
sd_wt = round(sd(weight_g), 1),
se = round(sd(weight_g) / sqrt(n()), 1),
.groups = "drop"
)
## # A tibble: 4 × 5
## diet n mean_wt sd_wt se
## <fct> <int> <dbl> <dbl> <dbl>
## 1 Diet 1 16 178. 58.7 14.7
## 2 Diet 2 10 215. 78.1 24.7
## 3 Diet 3 10 270. 71.6 22.6
## 4 Diet 4 9 239. 43.3 14.4
Write these numbers down. Diet 1 averages 178 g, Diet 3 averages 270 g, and the standard deviations run between 43 and 78 g. Every test in this document is a formal way of asking whether differences of that size, in samples this small, are bigger than what noise alone would produce.
Look at the data before testing it.
ggplot(chicks, aes(x = day, y = weight_g, group = chick_id, colour = diet)) +
geom_line(alpha = 0.4) +
facet_wrap(~ diet) +
labs(title = "Every chick's growth trajectory, by diet",
subtitle = "One line per chick; lines ending early are dropouts",
x = "Day", y = "Weight (g)") +
theme_minimal() +
theme(legend.position = "none")
This is the section to slow down on. The tests are five lines of code each; the reasoning underneath them is the actual content of the week.
We have 45 chicks. We do not care about these 45 chicks. We care about chicks, in general, raised on these diets — a population we will never observe. The 45 are a sample, and any statement we make about the population has to be filtered through the fact that a different 45 chicks would have given us different numbers.
So when we observe that Diet 3 chicks average 92.6 g more than Diet 1 chicks, there are two possibilities on the table:
Inference is a procedure for deciding between these two, with a stated and controlled rate of being wrong. It does not tell you which one is true.
The move that makes this tractable is counterintuitive: we start by assuming the thing we probably don’t believe.
\[H_0: \mu_{\text{Diet 1}} = \mu_{\text{Diet 3}}\]
The null hypothesis states no difference in the population. We assume it, and then ask how strange our sample would look if it were true. Why work backwards like this? Because “there is no difference” is a single, precise claim — it pins the population down to one exact value, and from a pinned-down population you can calculate probabilities. “There is a difference” is infinitely many claims (a difference of 1 g? 50 g? 200 g?), and you cannot compute a probability from it.
The alternative hypothesis, \(H_A\), is simply “not \(H_0\).”
Assume the null is true. Assume our sample was drawn at random from that population. Now:
p-value = P(observing a result at least as extreme as ours | \(H_0\) is true in the population)
Read that conditional bar carefully, because nearly every misinterpretation of a p-value comes from reading it backwards.
The p-value is not:
It is one thing only: how surprising your data would be if the null were true. Small p-value: this sample would be strange in a world where the null holds. Large p-value: this sample is perfectly ordinary in a world where the null holds — which is not the same as evidence that the null holds.
We compare the p-value to a cutoff, \(\alpha\), chosen before looking at the data, and do exactly one of two things:
Note what is missing: “accept \(H_0\).” We never accept the null. Absence of evidence for a difference is not evidence of absence, and later we show exactly how badly that can go with samples this size.
| \(H_0\) true in population | \(H_0\) false in population | |
|---|---|---|
| Reject \(H_0\) | Type I error (prob = \(\alpha\)) | Correct |
| Fail to reject \(H_0\) | Correct | Type II error (prob = \(\beta\)) |
Both are conditional probabilities, conditioned on a state of the world you can never observe. That is why you cannot know whether this particular test made an error. You only control the long-run rate.
Setting \(\alpha = 0.05\) is a statement about how often you are willing to be fooled: if the null is true, I accept a 5% chance of rejecting it anyway.
That last claim is checkable. Let’s build a world where the null is exactly true — two groups drawn from the same distribution, guaranteed no real difference — and run 5,000 t-tests on it.
set.seed(2026)
null_p <- replicate(5000, {
group_a <- rnorm(20, mean = 200, sd = 60)
group_b <- rnorm(20, mean = 200, sd = 60) # identical population
t.test(group_a, group_b)$p.value
})
mean(null_p < 0.05) # our realised Type I error rate at alpha = 0.05
## [1] 0.05
mean(null_p < 0.01) # ... and at alpha = 0.01
## [1] 0.0104
Exactly 5.0% of tests rejected a null that was true by construction, and 1.04% at the stricter cutoff. Nothing was wrong with those 250 tests — no bad data, no coding error. That rate is the price of admission, and \(\alpha\) is where you set it.
ggplot(data.frame(p = null_p), aes(x = p)) +
geom_histogram(binwidth = 0.05, boundary = 0, fill = "grey40", colour = "white") +
geom_vline(xintercept = 0.05, linetype = "dashed", colour = "red", linewidth = 0.8) +
labs(title = "5,000 p-values from a world where the null is TRUE",
subtitle = "Flat, by construction. The bar left of the red line is your Type I error rate.",
x = "p-value", y = "Count") +
theme_minimal()
The histogram is flat. Under a true null, every p-value is equally likely — 0.03 is exactly as probable as 0.83. This is the single most useful picture in introductory inference, because it makes two things obvious at once: a small p-value is not rare under the null in any absolute sense, and if you run twenty tests looking for something interesting, you should expect one “significant” result for free.
Which test to run is determined by the variable types, not by what you’re hoping to find. Two questions get you there:
| Exposure (independent) | Outcome (dependent) | Test | In R |
|---|---|---|---|
| Continuous | Continuous | Pearson correlation | cor.test() |
| Categorical, 2 levels, independent groups | Continuous | Two-sample t-test | t.test(y ~ g) |
| Categorical, 2 levels, same units measured twice | Continuous | Paired t-test | t.test(y1, y2, paired = TRUE) |
| Categorical, 3+ levels | Continuous | One-way ANOVA | aov(y ~ g) |
The row that trips people up is the third. “Two groups” is not enough information — you have to know whether the two sets of numbers come from different units or the same units twice. That is a fact about the study design, and no amount of staring at the data will tell you. Section 6 shows what it costs to get it wrong.
All of these are parametric tests, which means they assume the outcome is roughly normally distributed within groups. We check that as we go.
Pearson’s r summarises the linear association between two continuous variables in one number between -1 and +1. The sign is the direction; the absolute value is the strength; 0 means no linear relationship.
The associated test has
\[H_0: \rho = 0\]
where \(\rho\) is the correlation in the population. So the p-value answers “how surprising would a correlation this far from zero be, if the true correlation were exactly zero?”
Three warnings, in order of how often they cause trouble:
r measures linear association only. A
perfect U-shaped relationship can produce r ≈ 0. Always plot first.cor(x, y) is identical to
cor(y, x).ggplot(chicks, aes(x = day, y = weight_g)) +
geom_point(alpha = 0.35) +
geom_smooth(method = "lm", se = TRUE, formula = y ~ x) +
labs(title = "Weight against day, all chicks pooled",
x = "Day", y = "Weight (g)") +
theme_minimal()
cor.test(~ weight_g + day, data = chicks) %>% tidy()
## # A tibble: 1 × 8
## estimate statistic p.value parameter conf.low conf.high method alternative
## <dbl> <dbl> <dbl> <int> <dbl> <dbl> <chr> <chr>
## 1 0.837 36.7 5.02e-153 576 0.811 0.860 Pearson… two.sided
r = 0.837, p < 0.001. Chicks get heavier over time — a strong positive linear association, and completely unsurprising.
Note the cloud fanning out to the right: variability in weight grows with age. That widening is worth remembering, because “equal variance” is an assumption we lean on later.
Now something more useful. Is a chick that was heavy on day 16 also heavy on day 21?
ggplot(chicks_paired, aes(x = day_16, y = day_21)) +
geom_point(alpha = 0.7, size = 2) +
geom_smooth(method = "lm", se = TRUE, formula = y ~ x) +
labs(title = "Day-16 weight predicts day-21 weight",
subtitle = "One point per chick (n = 45)",
x = "Weight at day 16 (g)", y = "Weight at day 21 (g)") +
theme_minimal()
cor.test(~ day_16 + day_21, data = chicks_paired) %>% tidy()
## # A tibble: 1 × 8
## estimate statistic p.value parameter conf.low conf.high method alternative
## <dbl> <dbl> <dbl> <int> <dbl> <dbl> <chr> <chr>
## 1 0.919 15.2 6.07e-19 43 0.856 0.955 Pearson'… two.sided
r = 0.919, p < 0.001. Very strong. Heavy chicks stay heavy.
Hold onto this result — it is the entire justification for the next section. The two measurements are not independent pieces of information about a chick; they are two looks at the same chick, and they move together.
cor.test(~ day_0 + day_21, data = chicks_paired) %>% tidy()
## # A tibble: 1 × 8
## estimate statistic p.value parameter conf.low conf.high method alternative
## <dbl> <dbl> <dbl> <int> <dbl> <dbl> <chr> <chr>
## 1 -0.302 -2.08 0.0437 43 -0.547 -0.00942 Pearson's… two.sided
r = -0.302, p = 0.044. Two things here are worth more discussion than the result itself.
The p-value is 0.044. At \(\alpha = 0.05\) we reject. At \(\alpha = 0.01\) we do not. The data did not change; the conclusion did. The 0.05 cutoff is a convention, not a law of nature, and results that land near it deserve to be reported as “weak evidence,” not waved through as “significant.”
The correlation is negative, and weak. Hatch weight tells you almost nothing about day-21 weight, and if anything, the slightly heavier hatchlings ended up marginally lighter. With r = -0.30, hatch weight explains about 9% of the variance in final weight (\(r^2 = 0.091\)). Reporting this as “a significant relationship between hatch weight and final weight” would be technically true and substantively misleading.
Use a paired t-test when each observation in one group has a natural partner in the other:
The test does not actually compare two groups. It computes the difference within each pair, then runs a one-sample test asking whether the mean of those differences is zero:
\[H_0: \mu_{\text{difference}} = 0\]
Section 5 showed that day-16 and day-21 weights correlate at r = 0.919. That correlation is exactly what the paired test exploits. A chick that is 60 g above average on day 16 is probably 60-ish g above average on day 21 too — so when we subtract, that chick-specific “bigness” cancels out, and what remains is just growth.
You can see it in the standard deviations:
chicks_paired %>%
summarise(
sd_day16 = round(sd(day_16), 1),
sd_day21 = round(sd(day_21), 1),
sd_gain = round(sd(gain), 1)
)
## # A tibble: 1 × 3
## sd_day16 sd_day21 sd_gain
## <dbl> <dbl> <dbl>
## 1 47.4 71.5 33.7
The individual weights vary by 47 and 72 g. The differences vary by only 34 g. Subtracting removed most of the between-chick variation and left the within-chick change, which is what we actually asked about. That is the payoff: a smaller denominator, and a more sensitive test.
chicks_paired %>%
select(chick_id, diet, day_16, day_21) %>%
pivot_longer(c(day_16, day_21), names_to = "timepoint", values_to = "weight_g") %>%
ggplot(aes(x = timepoint, y = weight_g, group = chick_id)) +
geom_line(alpha = 0.4, colour = "grey40") +
geom_point(aes(colour = timepoint), size = 2, alpha = 0.8) +
labs(title = "Each chick's day-16 and day-21 weight",
subtitle = "Lines connect the same chick. Nearly all slope upward.",
x = NULL, y = "Weight (g)") +
theme_minimal() +
theme(legend.position = "none")
The connecting lines are the point of this plot. If you cannot draw them, you cannot run a paired test.
ggplot(chicks_paired, aes(x = gain)) +
geom_histogram(bins = 12, fill = "grey40", colour = "white") +
geom_vline(xintercept = 0, linetype = "dashed", colour = "red") +
labs(title = "Distribution of within-chick weight gain (day 21 − day 16)",
subtitle = "This is the variable the paired t-test actually analyses",
x = "Gain (g)", y = "Count") +
theme_minimal()
Roughly symmetric, no severe skew, and sitting well to the right of zero. The normality assumption applies to this distribution — the differences — not to the original weights.
t.test(chicks_paired$day_21, chicks_paired$day_16, paired = TRUE) %>% tidy()
## # A tibble: 1 × 8
## estimate statistic p.value parameter conf.low conf.high method alternative
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 49.2 9.79 1.27e-12 44 39.0 59.3 Paired t… two.sided
Interpretation. Chicks gained an average of 49.2 g between day 16 and day 21 (95% CI: 39.0 to 59.3), t(44) = 9.79, p < 0.001. We reject the null of no change. The confidence interval is the more informative half of that sentence: it excludes zero, which is why we rejected, but it also says the gain is somewhere between about 39 and 59 g — a statement about magnitude that the p-value alone does not give you.
Suppose we forget the design and treat these as two independent groups of chicks.
t.test(chicks_paired$day_21, chicks_paired$day_16, var.equal = TRUE) %>% tidy()
## # A tibble: 1 × 10
## estimate estimate1 estimate2 statistic p.value parameter conf.low conf.high
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 49.2 219. 170. 3.84 0.000228 88 23.7 74.6
## # ℹ 2 more variables: method <chr>, alternative <chr>
Compare:
| Estimate | t | df | p-value | 95% CI | |
|---|---|---|---|---|---|
| Paired (correct) | 49.2 | 9.79 | 44 | 1.3e-12 | 39.0 to 59.3 |
| Unpaired (wrong) | 49.2 | 3.84 | 88 | 2.3e-04 | 23.7 to 74.6 |
The point estimate is identical — both compute the same 49.2 g difference in means. Everything about the precision changes. The unpaired test throws away the information that these are the same 45 chicks, so all the between-chick variability lands in the standard error. The t-statistic falls from 9.79 to 3.84 and the confidence interval more than doubles in width.
Here both tests reject, so the conclusion survives. That is luck, and it is the dangerous case: the error is invisible when the effect is large. With a subtler effect the unpaired test would have missed it entirely.
The unpaired test is also wrong, not merely inefficient — it assumes independent observations, and these observations are not independent. Note the df of 88, as though we had 90 chicks. We have 45.
Now a genuine two-group comparison: different chicks, different diets, no pairing available.
\[H_0: \mu_{\text{Diet 1}} = \mu_{\text{Diet 3}}\]
The classic (“pooled” or “Student’s”) t-test assumes:
Assumption 2 is the one that earns this test its name. Under it, both groups estimate the same underlying variance, so we can pool them into a single estimate — which buys degrees of freedom and therefore power.
d13 <- chicks_d21 %>%
filter(diet %in% c("Diet 1", "Diet 3")) %>%
droplevels()
ggplot(d13, aes(x = diet, y = weight_g, fill = diet)) +
geom_boxplot(alpha = 0.5, outlier.shape = NA) +
geom_jitter(width = 0.15, alpha = 0.8, size = 2) +
labs(title = "Day-21 weight: Diet 1 vs Diet 3",
x = NULL, y = "Weight (g)") +
theme_minimal() +
theme(legend.position = "none")
Watch out:
droplevels()is not optional.filter()removes the rows for Diets 2 and 4 but leaves their empty levels attached to the factor.t.test()then reportsgrouping factor must have exactly 2 levelsand refuses to run. This is one of the most common errors in the whole workflow.
d13 %>%
group_by(diet) %>%
summarise(n = n(), mean = round(mean(weight_g), 1), sd = round(sd(weight_g), 1),
var = round(var(weight_g), 0), .groups = "drop")
## # A tibble: 2 × 5
## diet n mean sd var
## <fct> <int> <dbl> <dbl> <dbl>
## 1 Diet 1 16 178. 58.7 3446
## 2 Diet 3 10 270. 71.6 5130
var.test(weight_g ~ diet, data = d13) %>% tidy()
## # A tibble: 1 × 9
## estimate num.df den.df statistic p.value conf.low conf.high method alternative
## <dbl> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 0.672 15 9 0.672 0.476 0.178 2.10 F tes… two.sided
An F-test with \(H_0\): the two population variances are equal. F = 0.672, p = 0.476 — we do not reject, so the equal-variance assumption is defensible.
Read this test with some scepticism, though. With n = 16 and n = 10, it has very little power to detect unequal variances, so “we failed to reject” mostly reflects small samples rather than positive evidence of equality. It is also sensitive to non-normality. Treat it as one input alongside the boxplot and what you know about the measurement.
ggplot(d13, aes(sample = weight_g)) +
stat_qq(alpha = 0.8) + stat_qq_line(colour = "red") +
facet_wrap(~ diet) +
labs(title = "Q-Q plots: is weight roughly normal within each diet?",
x = "Theoretical quantiles", y = "Sample quantiles") +
theme_minimal()
Points track the line reasonably well in both groups. Good enough for a t-test, which is fairly robust to mild departures from normality.
t.test(weight_g ~ diet, data = d13, var.equal = TRUE) %>% tidy()
## # A tibble: 1 × 10
## estimate estimate1 estimate2 statistic p.value parameter conf.low conf.high
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 -92.6 178. 270. -3.60 0.00145 24 -146. -39.4
## # ℹ 2 more variables: method <chr>, alternative <chr>
Interpretation. Chicks on Diet 3 weighed on average 92.6 g more at day 21 than chicks on Diet 1 (270.3 g vs 177.8 g), t(24) = -3.60, p = 0.001. We reject the null of equal means. The 95% CI for the difference runs from -145.7 to -39.4 g and excludes zero.
The sign is negative because R sorted the factor levels and computed
Diet 1 − Diet 3. The direction of the comparison is
determined by your factor ordering, not by anything meaningful, so
always check estimate1 and estimate2 before
writing up which group came out ahead.
d12 <- chicks_d21 %>%
filter(diet %in% c("Diet 1", "Diet 2")) %>%
droplevels()
t.test(weight_g ~ diet, data = d12, var.equal = TRUE) %>% tidy()
## # A tibble: 1 × 10
## estimate estimate1 estimate2 statistic p.value parameter conf.low conf.high
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 -36.9 178. 215. -1.38 0.182 24 -92.4 18.5
## # ℹ 2 more variables: method <chr>, alternative <chr>
Diet 2 chicks averaged 36.9 g more than Diet 1 chicks — a real difference in our sample, about 21% of Diet 1’s mean — and p = 0.182. We fail to reject.
Be careful with the sentence you write next. This is not evidence that the diets are equivalent. The 95% CI runs from -92.4 to +18.5 g, meaning the data are consistent with Diet 2 being anywhere from substantially better to slightly worse. That interval is compatible with a large effect; we simply do not have the sample size to distinguish it from zero. Section 10 quantifies exactly how badly.
Drop the equal-variance assumption and you get Welch’s t-test. It estimates each group’s variance separately instead of pooling, and pays for that with a fractional, downward-adjusted degrees of freedom.
This matters most when unequal variances line up with unequal sample sizes. If the smaller group also has the larger variance, the pooled test understates the standard error and the Type I error rate drifts above the nominal 5%.
In R this is one argument — and note the default:
t.test(weight_g ~ diet, data = d13, var.equal = FALSE) %>% tidy()
## # A tibble: 1 × 10
## estimate estimate1 estimate2 statistic p.value parameter conf.low conf.high
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 -92.6 178. 270. -3.43 0.00334 16.4 -150. -35.5
## # ℹ 2 more variables: method <chr>, alternative <chr>
var.equal = FALSE is R’s default. If
you type t.test(y ~ g) you have run Welch’s test, not
Student’s. That surprises people arriving from SPSS or Stata, and it is
worth stating explicitly in your methods section which one you used.
| t | df | p-value | 95% CI | |
|---|---|---|---|---|
Pooled (var.equal = TRUE) |
-3.60 | 24.0 | 0.0015 | -145.7 to -39.4 |
Welch (var.equal = FALSE) |
-3.43 | 16.4 | 0.0033 | -149.6 to -35.5 |
Same estimate, same conclusion, slightly wider interval and a p-value about twice as large. The degrees of freedom dropped from 24 to 16.4 — that fractional df is Welch’s signature, and it is the cost of not assuming equal variances.
The conclusion is unchanged here, which is the usual outcome. Because of that, and because the equal-variance assumption is rarely well supported by a low-powered F-test, a defensible default is to just use Welch’s. You give up a little power when variances really are equal, and you avoid an inflated error rate when they are not.
Now the second comparison:
t.test(weight_g ~ diet, data = d12, var.equal = FALSE) %>% tidy()
## # A tibble: 1 × 10
## estimate estimate1 estimate2 statistic p.value parameter conf.low conf.high
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 -36.9 178. 215. -1.29 0.218 15.3 -98.1 24.2
## # ℹ 2 more variables: method <chr>, alternative <chr>
p = 0.218 versus 0.182 pooled. Same conclusion, again.
With four diets there are six possible pairwise comparisons. Run all six at \(\alpha = 0.05\) and, if the null is true everywhere, the probability of at least one false positive is roughly \(1 - 0.95^6 = 0.26\). Not 5%. A one-in-four chance of “finding” something that isn’t there.
ANOVA asks a single omnibus question instead:
\[H_0: \mu_1 = \mu_2 = \mu_3 = \mu_4\]
\(H_A\) is not “all four differ.” It is “at least one differs from at least one other” — which means a significant ANOVA tells you that something is going on, never what.
Despite the name, ANOVA compares means. It does it by partitioning total variability into two parts:
\[F = \frac{\text{between-group variance}}{\text{within-group variance}}\]
If the diets do nothing, both pieces estimate the same noise and F sits near 1. If the diets matter, the numerator inflates and F grows. The p-value asks how often an F this large would arise by chance under the null.
Assumptions are the t-test’s, extended: normality within groups, equal variances across all groups (homoscedasticity), independent observations.
ggplot(chicks_d21, aes(x = diet, y = weight_g, fill = diet)) +
geom_boxplot(alpha = 0.5, outlier.shape = NA) +
geom_jitter(width = 0.15, alpha = 0.8, size = 2) +
stat_summary(fun = mean, geom = "point", shape = 18, size = 4, colour = "black") +
labs(title = "Day-21 weight across all four diets",
subtitle = "Black diamonds mark group means",
x = NULL, y = "Weight (g)") +
theme_minimal() +
theme(legend.position = "none")
fit <- aov(weight_g ~ diet, data = chicks_d21)
bartlett.test(weight_g ~ diet, data = chicks_d21) %>% tidy() # equal variances?
## # A tibble: 1 × 4
## statistic p.value parameter method
## <dbl> <dbl> <dbl> <chr>
## 1 3.05 0.384 3 Bartlett test of homogeneity of variances
shapiro.test(residuals(fit)) %>% tidy() # normal residuals?
## # A tibble: 1 × 3
## statistic p.value method
## <dbl> <dbl> <chr>
## 1 0.988 0.929 Shapiro-Wilk normality test
Bartlett’s test: p = 0.384, no evidence against equal variances. Shapiro-Wilk on the residuals: p = 0.929, no evidence against normality. Both assumptions look fine, so we can read the standard ANOVA table with a clear conscience. (Note we test normality on the residuals of the fitted model, not on the raw outcome — the raw weights are a mixture across four groups and needn’t be normal.)
tidy(fit)
## # A tibble: 2 × 6
## term df sumsq meansq statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 diet 3 57164. 19055. 4.65 0.00686
## 2 Residuals 41 167839. 4094. NA NA
glance(fit)
## # A tibble: 1 × 6
## logLik AIC BIC deviance nobs r.squared
## <dbl> <dbl> <dbl> <dbl> <int> <dbl>
## 1 -249. 508. 517. 167839. 45 0.254
Interpretation. Mean day-21 weight differed significantly across the four diets, F(3, 41) = 4.66, p = 0.007. We reject the null that all four diets produce the same mean weight. Diet explains about 25% of the variance in day-21 weight (\(R^2\) = 0.254).
And now we are stuck. We know at least one difference exists. We do not know which.
If you are worried about the equal-variance assumption,
oneway.test() is the ANOVA analogue of Welch’s t-test:
oneway.test(weight_g ~ diet, data = chicks_d21, var.equal = FALSE) %>% tidy()
## # A tibble: 1 × 5
## num.df den.df statistic p.value method
## <dbl> <dbl> <dbl> <dbl> <chr>
## 1 3 20.6 4.66 0.0122 One-way analysis of means (not assuming equal…
F = 4.66, p = 0.012, fractional denominator df = 20.6. Same conclusion.
Tukey’s Honest Significant Difference runs all six pairwise comparisons while holding the family-wise error rate at 5% — the probability of making at least one Type I error across the whole set, not per comparison.
TukeyHSD(fit) %>% tidy() %>% select(contrast, estimate, conf.low, conf.high, adj.p.value)
## # A tibble: 6 × 5
## contrast estimate conf.low conf.high adj.p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Diet 2-Diet 1 37.0 -32.1 106. 0.487
## 2 Diet 3-Diet 1 92.6 23.5 162. 0.00470
## 3 Diet 4-Diet 1 60.8 -10.6 132. 0.119
## 4 Diet 3-Diet 2 55.6 -21.0 132. 0.226
## 5 Diet 4-Diet 2 23.9 -54.9 103. 0.849
## 6 Diet 4-Diet 3 -31.7 -110. 47.0 0.704
Only one comparison survives: Diet 3 vs Diet 1, a difference of 92.6 g (95% CI: 23.5 to 161.6), adjusted p = 0.005. Every other pair, including Diet 4 vs Diet 1 at 60.8 g, has an adjusted p above 0.05.
Two things to take from this table.
The adjustment has teeth. In Section 7 the unadjusted Diet 1 vs Diet 2 test gave p = 0.182; Tukey reports 0.487. Same data, same comparison — the difference is entirely the correction for having looked at six comparisons instead of one. If you go hunting through pairwise tests after a significant ANOVA without adjusting, you are back to the 26% error rate we started this section with.
A significant omnibus test does not guarantee a significant pairwise one. The overall ANOVA detected structure that Tukey can only localise to a single pair. That is normal, and it is honest: the omnibus test pools evidence across all groups, so it can pick up a pattern too diffuse for any individual comparison to nail down.
chicks_d21 %>%
group_by(diet) %>%
summarise(mean_wt = mean(weight_g), se = sd(weight_g) / sqrt(n()), .groups = "drop") %>%
ggplot(aes(x = fct_reorder(diet, mean_wt), y = mean_wt)) +
geom_col(fill = "steelblue", alpha = 0.85, width = 0.6) +
geom_errorbar(aes(ymin = mean_wt - 1.96 * se, ymax = mean_wt + 1.96 * se), width = 0.15) +
coord_flip() +
labs(title = "Mean day-21 weight by diet, with 95% confidence intervals",
subtitle = "Only Diet 3 vs Diet 1 survives Tukey adjustment",
x = NULL, y = "Mean weight (g)") +
theme_minimal()
Section 7 gave a 36.9 g difference with p = 0.182, and we said “fail to reject” rather than “no difference.” Here is why that distinction is not pedantry.
Power is \(1 - \beta\): the probability of rejecting the null given that it is false. It depends on three things you can estimate before collecting data — effect size, variability, and sample size — plus your chosen \(\alpha\).
# Pooled SD across the two groups
s1 <- chicks_d21 %>% filter(diet == "Diet 1") %>% pull(weight_g)
s2 <- chicks_d21 %>% filter(diet == "Diet 2") %>% pull(weight_g)
sp <- sqrt(((length(s1) - 1) * var(s1) + (length(s2) - 1) * var(s2)) /
(length(s1) + length(s2) - 2))
round(sp, 1)
## [1] 66.7
# Power to detect a 36.9 g difference at our actual sample size
power.t.test(n = 13, delta = 36.95, sd = sp, sig.level = 0.05)
##
## Two-sample t test power calculation
##
## n = 13
## delta = 36.95
## sd = 66.65805
## sig.level = 0.05
## power = 0.2732614
## alternative = two.sided
##
## NOTE: n is number in *each* group
Power = 0.27. If Diet 2 genuinely produces chicks 36.9 g heavier, this study had a 27% chance of detecting it. Which means a 73% chance of a Type II error — of concluding “no significant difference” when a real difference exists.
That reframes the finding entirely. The correct reading of p = 0.182 is not “the diets are equivalent.” It is “this study was never capable of answering the question.” Failing to detect something with a 27%-powered test is roughly as informative as failing to find your keys with the lights off.
We can confirm it by simulation — build a world where the difference is real and see how often we catch it:
set.seed(2026)
power_sim <- replicate(5000, {
a <- rnorm(16, mean = 178, sd = 66) # Diet 1, n = 16
b <- rnorm(10, mean = 215, sd = 66) # Diet 2, n = 10 -- 37 g heavier, truly
t.test(a, b)$p.value
})
mean(power_sim < 0.05) # simulated power
## [1] 0.2672
27%, matching the analytic calculation. In 73% of these simulated studies, a real 37 g effect went undetected.
power.t.test(delta = 36.95, sd = sp, sig.level = 0.05, power = 0.80)
##
## Two-sample t test power calculation
##
## n = 52.06544
## delta = 36.95
## sd = 66.65805
## sig.level = 0.05
## power = 0.8
## alternative = two.sided
##
## NOTE: n is number in *each* group
About 52 chicks per group for the conventional 80% power — against the 16 and 10 we actually had.
This is why power calculations belong at the design stage. Run one after a null result and all you learn is that you already wasted the study.
\(\alpha\) and \(\beta\) trade off against each other. Lower \(\alpha\) to 0.01 and you make fewer false positives and more false negatives. The only way to reduce both at once is a larger sample or a more precise measurement.
So when you read “no significant difference,” always ask the second question: was this study capable of finding one?
Try these before checking the answers below. Full solutions are in the chunk that follows each question.
Do chicks on Diet 2 and Diet 4 differ in day-21 weight? State the null, choose a test and justify the choice, check the assumptions, and interpret.
d24 <- chicks_d21 %>%
filter(diet %in% c("Diet 2", "Diet 4")) %>%
droplevels()
d24 %>% group_by(diet) %>%
summarise(n = n(), mean = round(mean(weight_g), 1), sd = round(sd(weight_g), 1), .groups = "drop")
## # A tibble: 2 × 4
## diet n mean sd
## <fct> <int> <dbl> <dbl>
## 1 Diet 2 10 215. 78.1
## 2 Diet 4 9 239. 43.3
var.test(weight_g ~ diet, data = d24) %>% tidy() %>% select(statistic, p.value)
## # A tibble: 1 × 2
## statistic p.value
## <dbl> <dbl>
## 1 3.25 0.112
t.test(weight_g ~ diet, data = d24, var.equal = TRUE) %>% tidy()
## # A tibble: 1 × 10
## estimate estimate1 estimate2 statistic p.value parameter conf.low conf.high
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 -23.9 215. 239. -0.809 0.430 17 -86.1 38.3
## # ℹ 2 more variables: method <chr>, alternative <chr>
Answer. \(H_0: \mu_{\text{Diet 2}} = \mu_{\text{Diet 4}}\). Two independent groups of different chicks with a continuous outcome, so an independent-samples t-test. The F-test gives F = 3.25, p = 0.112 — the SDs do differ noticeably (78.1 vs 43.3), but not enough to reject equal variances at this sample size, and Welch gives the same conclusion (p = 0.418). Diet 4 chicks averaged 23.9 g more (238.6 vs 214.7), t(17) = -0.81, p = 0.430, 95% CI -86.0 to 38.3. We fail to reject. Given that a 37 g difference had only 27% power in Section 10, a 24 g difference at n = 10 and 9 has even less — this is uninformative rather than evidence of equivalence.
Does weight gain between day 16 and day 21 differ across the four diets? How is this different from the ANOVA in Section 9?
gain_fit <- aov(gain ~ diet, data = chicks_paired)
tidy(gain_fit)
## # A tibble: 2 × 6
## term df sumsq meansq statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 diet 3 10486. 3495. 3.64 0.0204
## 2 Residuals 41 39402. 961. NA NA
chicks_paired %>%
group_by(diet) %>%
summarise(n = n(), mean_gain = round(mean(gain), 1), sd = round(sd(gain), 1), .groups = "drop")
## # A tibble: 4 × 4
## diet n mean_gain sd
## <fct> <int> <dbl> <dbl>
## 1 Diet 1 16 31.9 23.4
## 2 Diet 2 10 50 35.5
## 3 Diet 3 10 72.9 38.8
## 4 Diet 4 9 52.4 28.1
Answer. F(3, 41) = 3.64, p = 0.020 — reject the null that mean gain is equal across diets. Gains ran from 31.9 g (Diet 1) to 72.9 g (Diet 3). This differs from Section 9 in a way that matters: the outcome is now a within-chick change, so each chick serves as its own control and differences in starting weight are removed. It combines the logic of Section 6 (pairing) with the logic of Section 9 (comparing more than two groups). When you have repeated measures on units in different groups, analysing the change rather than the endpoint is usually the more powerful and more interpretable choice.
A colleague runs all six pairwise t-tests, finds Diet 3 vs Diet 1 at p = 0.0015, and reports it as significant at \(\alpha = 0.05\) with no adjustment. What is wrong, and what is the fix?
Answer. Each test carries a 5% Type I error rate, so across six tests the family-wise error rate is about \(1 - 0.95^6 = 26\%\). Reporting the smallest of six p-values against an unadjusted 0.05 cutoff overstates the evidence. The fix is an omnibus ANOVA followed by an adjusted post-hoc procedure (Tukey, Bonferroni, Holm). Here the honest result is Tukey’s adjusted p = 0.005 — still significant, so the substantive conclusion survives, but that was not guaranteed in advance. Note that Diet 1 vs Diet 2 went from an unadjusted 0.182 to an adjusted 0.487.
In Section 2 we found that 5 chicks are missing day-21 weights, 4 of them on Diet 1. Why does that complicate the Section 9 conclusion?
Answer. The missingness is not random with respect to the exposure. Our day-21 analysis covers only chicks that survived to day 21, so the estimand quietly shifted from “the effect of diet on weight” to “the effect of diet on weight among chicks that made it.” If chicks dropped out because they were failing to thrive, Diet 1’s surviving average is biased upward — its weakest members were removed — which would make Diet 1 look better than it is and shrink the estimated gap to Diet 3. The direction of bias here works against our finding, which is reassuring, but the general lesson holds: report the dropout, state which population your estimate actually describes, and never let
drop_na()be a silent step.
| Exposure | Outcome | Test | R |
|---|---|---|---|
| Continuous | Continuous | Pearson correlation | cor.test(~ x + y, data) |
| 2 groups, independent | Continuous | Two-sample t-test | t.test(y ~ g, data, var.equal =) |
| 2 measures, same units | Continuous | Paired t-test | t.test(y1, y2, paired = TRUE) |
| 3+ groups, independent | Continuous | One-way ANOVA | aov(y ~ g, data) |
| Follow-up to ANOVA | Continuous | Tukey HSD | TukeyHSD(fit) |
| Question | Test | Result | Decision |
|---|---|---|---|
| Weight vs day | Pearson r | r = 0.84, p < 0.001 | Reject \(H_0\) |
| Day-16 vs day-21 weight | Pearson r | r = 0.92, p < 0.001 | Reject \(H_0\) |
| Day-0 vs day-21 weight | Pearson r | r = -0.30, p = 0.044 | Reject at 0.05, not at 0.01 |
| Weight change, day 16 to 21 | Paired t | +49.2 g, p < 0.001 | Reject \(H_0\) |
| Diet 1 vs Diet 3 | Pooled t | -92.6 g, p = 0.0015 | Reject \(H_0\) |
| Diet 1 vs Diet 3 | Welch t | -92.6 g, p = 0.0033 | Reject \(H_0\) |
| Diet 1 vs Diet 2 | Pooled t | -36.9 g, p = 0.182 | Fail to reject (power = 0.27) |
| All four diets | ANOVA | F(3,41) = 4.66, p = 0.007 | Reject \(H_0\) |
| Which pair? | Tukey | Only Diet 3 vs Diet 1, p = 0.005 | — |
var.equal = FALSE is R’s default. Welch’s is a
reasonable default in general.sessionInfo()
## R version 4.5.2 (2025-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/Denver
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] broom_1.0.12 labelled_2.16.0 forcats_1.0.1 ggplot2_4.0.2
## [5] tidyr_1.3.2 dplyr_1.2.0
##
## loaded via a namespace (and not attached):
## [1] Matrix_1.7-4 gtable_0.3.6 jsonlite_2.0.0 compiler_4.5.2
## [5] tidyselect_1.2.1 stringr_1.6.0 jquerylib_0.1.4 splines_4.5.2
## [9] scales_1.4.0 yaml_2.3.12 fastmap_1.2.0 lattice_0.22-9
## [13] R6_2.6.1 labeling_0.4.3 generics_0.1.4 knitr_1.51
## [17] backports_1.5.0 tibble_3.3.1 bslib_0.10.0 pillar_1.11.1
## [21] RColorBrewer_1.1-3 rlang_1.1.7 utf8_1.2.6 stringi_1.8.7
## [25] cachem_1.1.0 xfun_0.56 sass_0.4.10 S7_0.2.1
## [29] otel_0.2.0 cli_3.6.1 mgcv_1.9-4 withr_3.0.2
## [33] magrittr_2.0.4 digest_0.6.39 grid_4.5.2 rstudioapi_0.18.0
## [37] haven_2.5.5 hms_1.1.4 nlme_3.1-168 lifecycle_1.0.5
## [41] vctrs_0.7.1 evaluate_1.0.5 glue_1.8.0 farver_2.1.2
## [45] rmarkdown_2.30 purrr_1.2.1 tools_4.5.2 pkgconfig_2.0.3
## [49] htmltools_0.5.9