0. Setup

Load libraries

library(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.


1. The data

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, …

2. Data preparation

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.

Making missingness visible

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.

Analytic datasets

# 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")


3. The logic of statistical inference

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.

The problem inference solves

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:

  1. Diet 3 really does produce heavier chicks, and our sample reflects that.
  2. The diets are equivalent, and we happened to draw a sample where the Diet 3 chicks came out ahead.

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 null hypothesis

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\).”

The p-value

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:

  • the probability that \(H_0\) is true
  • the probability that your finding is due to chance
  • the probability you are wrong
  • a measure of effect size

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.

The decision rule

We compare the p-value to a cutoff, \(\alpha\), chosen before looking at the data, and do exactly one of two things:

  • Reject \(H_0\) — equality is not a good description of the relationship.
  • Fail to reject \(H_0\) — we cannot rule out equality.

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.

Two ways to be right, two ways to be wrong

\(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\))
  • \(\alpha = P(\text{reject } H_0 \mid H_0 \text{ true})\) — a false positive
  • \(\beta = P(\text{fail to reject } H_0 \mid H_0 \text{ false})\) — a false negative
  • Power \(= 1 - \beta\) — the chance of detecting a real effect

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.

Seeing \(\alpha\) do its job

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.


4. Choosing a test

Which test to run is determined by the variable types, not by what you’re hoping to find. Two questions get you there:

  1. What kind of variable is the outcome?
  2. What kind of variable is the exposure, and if it’s categorical, how many levels?
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.


5. Correlation: two continuous variables

What it measures

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:

  1. r measures linear association only. A perfect U-shaped relationship can produce r ≈ 0. Always plot first.
  2. The p-value and the effect size are different questions. With a large enough sample, an r of 0.05 gets a tiny p-value. Statistically significant and practically meaningless are compatible.
  3. Correlation is symmetric and says nothing about direction of causation. cor(x, y) is identical to cor(y, x).

Example A: weight and day

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.

Example B: day-16 weight and day-21 weight

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.

Example C: day-0 weight and day-21 weight

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.


6. The paired t-test

When to use it

Use a paired t-test when each observation in one group has a natural partner in the other:

  • the same subject measured before and after an intervention
  • the same subject measured under two conditions
  • matched pairs (twins, case-control pairs matched on age and sex)

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\]

Why it’s a different test

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.

Look at the data

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.

Run it

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.

What ignoring the pairing costs

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.


7. The independent-samples t-test (equal variances)

Setup

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:

  1. The outcome is roughly normal within each group.
  2. The two groups have equal population variances.
  3. Observations are independent, within and between groups.

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.

Check the assumptions first

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 reports grouping factor must have exactly 2 levels and 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.

Run it

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.

A comparison that does not reach significance

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.


8. The independent-samples t-test (unequal variances)

Welch’s correction

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.

Compare the two

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.


9. One-way ANOVA

Why not just run all the t-tests

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.

How it works

Despite the name, ANOVA compares means. It does it by partitioning total variability into two parts:

  • Between-group variance — how far the group means sit from the grand mean
  • Within-group variance — how much chicks vary around their own group’s mean

\[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.

Look at the data

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")

Check the assumptions

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.)

Run it

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.

Post-hoc: which groups differ?

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()


10. Type II error and power

Returning to Diet 1 vs Diet 2

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.

How many chicks would we have needed?

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.

The takeaway

\(\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?


11. Practice

Try these before checking the answers below. Full solutions are in the chunk that follows each question.

Question 1

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.

Question 2

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.

Question 3

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.

Question 4

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.


12. Summary

Test selection

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)

What we found

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

Nine things to carry forward

  1. The p-value is P(data at least this extreme | \(H_0\) true). Never the reverse.
  2. Under a true null, p-values are uniform — 5% land below 0.05 by construction.
  3. “Fail to reject” is not “accept.” It often just means the sample was too small.
  4. A confidence interval carries the p-value’s information plus a statement about magnitude. Report it.
  5. Variable types choose the test. Study design chooses between paired and independent.
  6. Ignoring pairing inflates the standard error and violates the independence assumption.
  7. var.equal = FALSE is R’s default. Welch’s is a reasonable default in general.
  8. A significant ANOVA says “somewhere,” not “where.” Post-hoc tests must adjust.
  9. Statistical significance is not effect size, and neither is practical importance.

Session info

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