Lab Overview

Time: ~30 minutes

Goal: Practice one-way ANOVA analysis from start to finish using real public health data

Learning Objectives:

Understand when and why to use ANOVA instead of multiple t-tests
Set up hypotheses for ANOVA
Conduct and interpret the F-test
Perform post-hoc tests when appropriate
Check ANOVA assumptions
Calculate and interpret effect size (η²)

Structure:

Part A: Guided Example (follow along)
Part B: Your Turn (independent practice)

Submission: Upload your completed .Rmd file and published to Brightspace by the end of class.

PART A: GUIDED EXAMPLE

Example: Blood Pressure and BMI Categories

Research Question: Is there a difference in mean systolic blood pressure (SBP) across three BMI categories (Normal weight, Overweight, Obese)?

Why ANOVA? We have one continuous outcome (SBP) and one categorical predictor with THREE groups (BMI category). Using multiple t-tests would inflate our Type I error rate.

Step 1: Setup and Data Preparation

# Load necessary libraries
library(tidyverse)   # For data manipulation and visualization
library(knitr)       # For nice tables
library(car)         # For Levene's test
library(NHANES)      # NHANES dataset

# Load the NHANES data
data(NHANES)

Create analysis dataset:

# Set seed for reproducibility
set.seed(553)

# Create BMI categories and prepare data
bp_bmi_data <- NHANES %>%
  filter(Age >= 18 & Age <= 65) %>%  # Adults 18-65
  filter(!is.na(BPSysAve) & !is.na(BMI)) %>%
  mutate(
    bmi_category = case_when(
      BMI < 25 ~ "Normal",
      BMI >= 25 & BMI < 30 ~ "Overweight",
      BMI >= 30 ~ "Obese",
      TRUE ~ NA_character_
    ),
    bmi_category = factor(bmi_category, 
                         levels = c("Normal", "Overweight", "Obese"))
  ) %>%
  filter(!is.na(bmi_category)) %>%
  select(ID, Age, Gender, BPSysAve, BMI, bmi_category)

# Display first few rows
head(bp_bmi_data) %>% 
  kable(caption = "Blood Pressure and BMI Dataset (first 6 rows)")

Blood Pressure and BMI Dataset (first 6 rows)
ID	Age	Gender	BPSysAve	BMI	bmi_category
51624	34	male	113	32.22	Obese
51624	34	male	113	32.22	Obese
51624	34	male	113	32.22	Obese
51630	49	female	112	30.57	Obese
51647	45	female	118	27.24	Overweight
51647	45	female	118	27.24	Overweight

# Check sample sizes
table(bp_bmi_data$bmi_category)

## 
##     Normal Overweight      Obese 
##       1939       1937       2150

Interpretation: We have 6026 adults with complete BP and BMI data across three BMI categories.

Step 2: Descriptive Statistics

# Calculate summary statistics by BMI category
summary_stats <- bp_bmi_data %>%
  group_by(bmi_category) %>%
  summarise(
    n = n(),
    Mean = mean(BPSysAve),
    SD = sd(BPSysAve),
    Median = median(BPSysAve),
    Min = min(BPSysAve),
    Max = max(BPSysAve)
  )

summary_stats %>% 
  kable(digits = 2, 
        caption = "Descriptive Statistics: Systolic BP by BMI Category")

Descriptive Statistics: Systolic BP by BMI Category
bmi_category	n	Mean	SD	Median	Min	Max
Normal	1939	114.23	15.01	113	78	221
Overweight	1937	118.74	13.86	117	83	186
Obese	2150	121.62	15.27	120	82	226

Observation: The mean SBP appears to increase from Normal (114.2) to Overweight (118.7) to Obese (121.6).

But is this difference statistically significant?

We can not answer if the difference of the groups means is statistically significant as it is merely descriptive statistics.

Step 3: Visualize the Data

# Create boxplots with individual points
ggplot(bp_bmi_data, 
  aes(x = bmi_category, y = BPSysAve, fill = bmi_category)) +
  geom_boxplot(alpha = 0.7, outlier.shape = NA) +
  geom_jitter(width = 0.2, alpha = 0.1, size = 0.5) +
  scale_fill_brewer(palette = "Set2") +
  labs(
    title = "Systolic Blood Pressure by BMI Category",
    subtitle = "NHANES Data, Adults aged 18-65",
    x = "BMI Category",
    y = "Systolic Blood Pressure (mmHg)",
    fill = "BMI Category"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none")

What the plot tells us:

There appears to be a trend: higher BMI categories have higher median SBP
The boxes overlap, but the obese group appears shifted upward
Variability (box heights) looks similar across groups

Step 4: Set Up Hypotheses

Null Hypothesis (H₀): μ_Normal = μ_Overweight = μ_Obese
(All three population means are equal)

Alternative Hypothesis (H₁): At least one population mean differs from the others

Significance level: α = 0.05

Step 5: Fit the ANOVA Model

# Fit the one-way ANOVA model
anova_model <- aov(BPSysAve ~ bmi_category, data = bp_bmi_data)

# Display the ANOVA table
summary(anova_model)

##                Df  Sum Sq Mean Sq F value Pr(>F)    
## bmi_category    2   56212   28106   129.2 <2e-16 ***
## Residuals    6023 1309859     217                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretation:

F-statistic: 129.24
Degrees of freedom: df₁ = 2 (k-1 groups), df₂ = 6023 (n-k)
p-value: < 2e-16 (very small)
Decision: Since p < 0.05, we reject H₀
Conclusion: There is statistically significant evidence that mean systolic BP differs across at least two BMI categories.

Step 6: Post-Hoc Tests (Tukey HSD)

Why do we need this? The F-test tells us that groups differ, but not which groups differ. Tukey’s Honest Significant Difference controls the family-wise error rate for multiple pairwise comparisons.

# Conduct Tukey HSD test
tukey_results <- TukeyHSD(anova_model)
print(tukey_results)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = BPSysAve ~ bmi_category, data = bp_bmi_data)
## 
## $bmi_category
##                       diff      lwr      upr p adj
## Overweight-Normal 4.507724 3.397134 5.618314     0
## Obese-Normal      7.391744 6.309024 8.474464     0
## Obese-Overweight  2.884019 1.801006 3.967033     0

# Visualize the confidence intervals
plot(tukey_results, las = 0)

Interpretation:

Comparison	Mean Diff	95% CI	p-value	Significant?
Overweight - Normal	4.51	[3.4, 5.62]	1.98e-13	Yes
Obese - Normal	7.39	[6.31, 8.47]	< 0.001	Yes
Obese - Overweight	2.88	[1.8, 3.97]	1.38e-09	Yes

Conclusion: All three pairwise comparisons are statistically significant. Obese adults have higher SBP than overweight adults, who in turn have higher SBP than normal-weight adults.

Step 7: Calculate Effect Size

# Extract sum of squares from ANOVA table
anova_summary <- summary(anova_model)[[1]]

ss_treatment <- anova_summary$`Sum Sq`[1]
ss_total <- sum(anova_summary$`Sum Sq`)

# Calculate eta-squared
eta_squared <- ss_treatment / ss_total

cat("Eta-squared (η²):", round(eta_squared, 4), "\n")

## Eta-squared (η²): 0.0411

cat("Percentage of variance explained:", round(eta_squared * 100, 2), "%")

## Percentage of variance explained: 4.11 %

Interpretation: BMI category explains 4.11% of the variance in systolic BP.

Effect size guidelines: Small (0.01), Medium (0.06), Large (0.14)
Our effect: Small

While statistically significant, the practical effect is modest—BMI category alone doesn’t explain most of the variation in blood pressure.

Step 8: Check Assumptions

ANOVA Assumptions:

Independence: Observations are independent (assumed based on study design)
Normality: Residuals are approximately normally distributed
Homogeneity of variance: Equal variances across groups

par(mfrow = c(2, 2), mar = c(2, 2, 1, 1))
plot(anova_model)

Diagnostic Plot Interpretation:

Residuals vs Fitted: Points show random scatter around zero with no clear pattern → Good!
Q-Q Plot: Points follow the diagonal line reasonably well → Normality assumption is reasonable
Scale-Location: Red line is relatively flat → Equal variance assumption is reasonable
Residuals vs Leverage: No points beyond Cook’s distance lines → No highly influential outliers

# Levene's test for homogeneity of variance
levene_test <- leveneTest(BPSysAve ~ bmi_category, data = bp_bmi_data)
print(levene_test)

## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value  Pr(>F)  
## group    2  2.7615 0.06328 .
##       6023                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Levene’s Test Interpretation:

p-value: 0.0633
If p < 0.05, we would reject equal variances
Here: Equal variance assumption is met

Overall Assessment: With n > 2000, ANOVA is robust to minor violations. Our assumptions are reasonably satisfied.

Step 9: Report Results

Example Results Section:

We conducted a one-way ANOVA to examine whether mean systolic blood pressure (SBP) differs across BMI categories (Normal, Overweight, Obese) among 6,026 adults aged 18-65 from NHANES. Descriptive statistics showed mean SBP of 114.2 mmHg (SD = 15) for normal weight, 118.7 mmHg (SD = 13.9) for overweight, and 121.6 mmHg (SD = 15.3) for obese individuals.

The ANOVA revealed a statistically significant difference in mean SBP across BMI categories, F(2, 6023) = 129.24, p < 0.001. Tukey’s HSD post-hoc tests indicated that all pairwise comparisons were significant (p < 0.05): obese adults had on average 7.4 mmHg higher SBP than normal-weight adults, and 2.9 mmHg higher than overweight adults.

The effect size (η² = 0.041) indicates that BMI category explains 4.1% of the variance in systolic blood pressure, representing a small practical effect. These findings support the well-established relationship between higher BMI and elevated blood pressure, though other factors account for most of the variation in SBP.

PART B: YOUR TURN - INDEPENDENT PRACTICE

Practice Problem: Physical Activity and Depression

Research Question: Is there a difference in the number of days with poor mental health across three physical activity levels (None, Moderate, Vigorous)?

Your Task: Complete the same 9-step analysis workflow you just practiced, but now on a different outcome and predictor.

Step 1: Data Preparation

# Prepare the dataset
set.seed(553)

mental_health_data <- NHANES %>%
  filter(Age >= 18) %>%
  filter(!is.na(DaysMentHlthBad) & !is.na(PhysActive)) %>%
  mutate(
    activity_level = case_when(
      PhysActive == "No" ~ "None",
      PhysActive == "Yes" & !is.na(PhysActiveDays) & PhysActiveDays < 3 ~ "Moderate",
      PhysActive == "Yes" & !is.na(PhysActiveDays) & PhysActiveDays >= 3 ~ "Vigorous",
      TRUE ~ NA_character_
    ),
    activity_level = factor(activity_level, 
                           levels = c("None", "Moderate", "Vigorous"))
  ) %>%
  filter(!is.na(activity_level)) %>%
  select(ID, Age, Gender, DaysMentHlthBad, PhysActive, activity_level)

# YOUR TURN: Display the first 6 rows and check sample sizes
# Display first few rows
head(mental_health_data) %>% 
  kable(caption = "Activity level and Days with Bad Mental Helth Dataset (first 6 rows)")

Activity level and Days with Bad Mental Helth Dataset (first 6 rows)
ID	Age	Gender	DaysMentHlthBad	PhysActive	activity_level
51624	34	male	15	No	None
51624	34	male	15	No	None
51624	34	male	15	No	None
51630	49	female	10	No	None
51647	45	female	3	Yes	Vigorous
51647	45	female	3	Yes	Vigorous

# Check sample sizes
table(mental_health_data$activity_level)

## 
##     None Moderate Vigorous 
##     3139      768     1850

YOUR TURN - Answer these questions:

How many people are in each physical activity group?
- None: 3139
- Moderate: 768
- Vigorous: 1850

Step 2: Descriptive Statistics

# Calculate summary statistics by BMI category
summary_stats <- mental_health_data %>%
  group_by(activity_level) %>%
  summarise(
    n = n(),
    Mean = mean(DaysMentHlthBad),
    SD = sd(DaysMentHlthBad),
    Median = median(DaysMentHlthBad),
    Min = min(DaysMentHlthBad),
    Max = max(DaysMentHlthBad)
  )

summary_stats %>% 
  kable(digits = 2, 
        caption = "Descriptive Statistics: Days Mental Health is Bad by activity level categories")

Descriptive Statistics: Days Mental Health is Bad by activity level categories
activity_level	n	Mean	SD	Max
None	3139	5.08	9.01	30
Moderate	768	3.81	6.87	30
Vigorous	1850	3.54	7.17	30

YOUR TURN - Interpret:

Which group has the highest mean number of bad mental health days? Group with activity_level NONE has the highest mean number of bad mental health days
Which group has the lowest? Group with activity_level Vigorous has the Vigorous has the lowest mean number of bad mental health days

Step 3: Visualization

ggplot(mental_health_data, 
  aes(x = activity_level, y = DaysMentHlthBad, fill = activity_level)) +
  geom_boxplot(alpha = 0.7, outlier.shape = NA) +
  geom_jitter(width = 0.2, alpha = 0.1, size = 0.5) +
  scale_fill_brewer(palette = "Set2") +
  labs(
    title = "Poor mental health outcomes by activity level",
    subtitle = "NHANES Data, Adults aged 18-65",
    x = "Activity level",
    y = "Poor Meantal Health (days)",
    fill = "BMI Category"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none")

YOUR TURN - Describe what you see:

Do the groups appear to differ? First of all the data is not normally distributed (central tendency theorem), as the boxplots are showing that the data is heavily skewed to the right and has no values in Q1 to the median. So the group “Moderate” activity level group has the highest frequency of bad mental health days. IQR of “Vigorous” activity level group is smaller, which shows that the middle values cluster more tightly (from median to Q3);Also “vigorous” activity level group has the smaller extreme value except outliers in comparison with other two groups.
Are the variances similar across groups? Variances are not similar across the groups because Q3 of the vigorous group is “closer” to the median and upper whisker is shorter in comparison with other two groups. —

Step 4: Set Up Hypotheses

YOUR TURN - Write the hypotheses:

Null Hypothesis (H₀): μ_None = μ_Moderate = μ_Vigorous (All three population means are equal)

Alternative Hypothesis (H₁): At least one population mean differs from the others

Significance level: α = 0.05

Step 5: Fit the ANOVA Model

# Fit the one-way ANOVA model
anova_model <- aov(DaysMentHlthBad ~ activity_level, data = mental_health_data)

# Display the ANOVA table
summary(anova_model)

##                  Df Sum Sq Mean Sq F value   Pr(>F)    
## activity_level    2   3109  1554.6   23.17 9.52e-11 ***
## Residuals      5754 386089    67.1                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

YOUR TURN - Extract and interpret the results:

F-statistic: 23.17
Degrees of freedom: 2
p-value: 9.52e-11
Decision (reject or fail to reject H₀): reject H₀
Statistical conclusion in words: We reject H₀, at least one population mean differs from the others

Step 6: Post-Hoc Tests

# Conduct Tukey HSD test
tukey_results <- TukeyHSD(anova_model)
print(tukey_results)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = DaysMentHlthBad ~ activity_level, data = mental_health_data)
## 
## $activity_level
##                         diff       lwr        upr     p adj
## Moderate-None     -1.2725867 -2.045657 -0.4995169 0.0003386
## Vigorous-None     -1.5464873 -2.109345 -0.9836298 0.0000000
## Vigorous-Moderate -0.2739006 -1.098213  0.5504114 0.7159887

# Visualize the confidence intervals
plot(tukey_results, las = 0)

YOUR TURN - Complete the table:

Comparison	Mean Difference	95% CI Lower	95% CI Upper	p-value	Significant?
Moderate - None
Vigorous - None
Vigorous - Moderate

Interpretation:

Which specific groups differ significantly? The difference between the vigorous and moderate groups is not statistically significant because p=0.71 and the CI of the difference crosses the null value. —

Step 7: Calculate Effect Size

# Extract sum of squares from ANOVA table
anova_summary <- summary(anova_model)[[1]]

ss_treatment <- anova_summary$`Sum Sq`[1]
ss_total <- sum(anova_summary$`Sum Sq`)

# Calculate eta-squared
eta_squared <- ss_treatment / ss_total

cat("Eta-squared (η²):", round(eta_squared, 4), "\n")

## Eta-squared (η²): 0.008

cat("Percentage of variance explained:", round(eta_squared * 100, 2), "%")

## Percentage of variance explained: 0.8 %

YOUR TURN - Interpret:

η² = 0.008
Percentage of variance explained: 0.8% Effect size classification (small/medium/large): η²=0.01 indicates a small effect, so we have very smaller effect as it equals to 0.008
What does this mean practically? Even though the result is statistically significant (p-value) in practice it has small to insufficient effect (η²-value).

Step 8: Check Assumptions

par(mfrow = c(2, 2))
plot(anova_model)

YOUR TURN - Evaluate each plot:

Residuals vs Fitted: You see three distinct vertical “columns” of points. Each column represents one of our groups (Nono, moderate, and vigorous). In ANOVA, every person in a group has the same “fitted value” (the group mean), red line is declining → Linearity assumption is not satisfied
Q-Q Plot: Points do not follow the diagonal line (S-shape)→ Normality assumption is not satisfied
Scale-Location: Red line is increasing → Equal variance assumption is not satisfied
Residuals vs Leverage: There are no points beyond Cook’s distance lines → There are no influential outliers

# YOUR TURN: Conduct Levene's test
levene_test <- leveneTest(DaysMentHlthBad ~ activity_level, data = mental_health_data)
print(levene_test)

## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value    Pr(>F)    
## group    2  23.168 9.517e-11 ***
##       5754                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

YOUR TURN - Overall assessment:

Are assumptions reasonably met? The Anova assumptions are not fully met. Overall we do not have balanced design (similar group sizes) and distributions are not similar across groups (non-normality), as figures show that data is not normally distributed and Levene’s test being statistically significance means that there is no homogeneity of the variances. In other words, the spread of values around each group’s mean are not roughly the same across all groups.
Do any violations threaten your conclusions? As ANOVA is robust to moderate violations of these assumptions, particularly in our case with with larger sample sizes (typically n > 30 per group). So the validity of our p-values and inference should be interpreted with caution, because the Assumptions of Anova are not fully supported by the data. As η² = 0.008, statistically significant result is achieved by large sample size rather than strong clinical effect —

Step 9: Write Up Results

YOUR TURN - Write a complete 2-3 paragraph results section:

Include: 1. Sample description and descriptive statistics 2. F-test results 3. Post-hoc comparisons (if applicable) 4. Effect size interpretation 5. Public health significance

Your Results Section:

Results

-We conducted a one-way ANOVA to examine whether the mean number of days with poor mental health differs across physical activity levels (None, Moderate, Vigorous) among 5,757 adults aged 18–65 from NHANES. -Descriptive statistics indicated that adults with no physical activity reported the highest mean number of days with poor mental health (5.1 days, SD = 9), followed by those with moderate physical activity (3.8 days, SD = 6.9). Adults engaging in vigorous physical activity reported the lowest mean number of poor mental health days (3.5 days, SD = 7.2). -The one-way ANOVA revealed a statistically significant difference in mean days of poor mental health across activity level groups, F (2, 5754) = 23.17, p < 0.001. Post-hoc comparisons using Tukey’s HSD test showed that adults with no physical activity reported significantly more days of poor mental health compared to those with moderate and vigorous activity levels (p < 0.05). In contrast, the difference between the moderate and vigorous activity groups was not statistically significant. -The effect size was small (η² = 0.008), indicating that physical activity level explains approximately 0.8% of the variance in days of poor mental health. This suggests that although the association is statistically significant, the practical significance is limited, and poor mental health outcomes are likely influenced by multiple complex factors beyond physical activity alone. —

Reflection Questions

1. How does the effect size help you understand the practical vs. statistical significance?

Along with statistical significance, eta-squared values should be taken into consideration for interpreting clinical/practical effect of the independent variable on outcomes.

2. Why is it important to check ANOVA assumptions? What might happen if they’re violated?

Before making conclusions about validity of p-values and inferences ANOVA requires meeting three key assumptions. For my understanding it shows if we even should use ANOVA test analyzing data if the variables are not normally distributed (normality) if the variances of the groups are not roughly similar, if there are outliers that heavilty infulence our results then we should consider transformations or use Welch’s ANOVA, especially if the sample size is quite small.

3. In public health practice, when might you choose to use ANOVA?

We should choose AnOVA for comparison of more than 2 groups, to see if there are any differences between the means of these groups. And we should use one way ANOVA only when the data is normally distributed and variances are roughly equal, and Welch’s ANOVA when there is statistically significant deifference in groups variances.

4. What was the most challenging part of this lab activity?

The interpretation of plots for testing ANOVA’s assumptions and theoretical background for interpretation of the inferences.

Submission Checklist

Before submitting, verify you have:

Completed all code chunks in Part B
Filled in all interpretation questions
Completed the Tukey HSD comparison table
Written a complete results section
Answered all reflection questions
Rendered your .Rmd file to HTML without errors
Checked that all outputs display correctly

To submit: Upload both your .Rmd file and the HTML output to Brightspace.

Lab completed on: February 05, 2026

In-Class Lab Activity: Analysis of Variance (ANOVA)

EPI 553: Principles of Statistical Inference II

Aigerim Azamatova

2026-02-05

Lab Overview

PART A: GUIDED EXAMPLE

Example: Blood Pressure and BMI Categories

Research Question: Is there a difference in mean systolic blood pressure (SBP) across three BMI categories (Normal weight, Overweight, Obese)?

Step 1: Setup and Data Preparation

Step 2: Descriptive Statistics

But is this difference statistically significant?

Step 3: Visualize the Data

Step 4: Set Up Hypotheses

Step 5: Fit the ANOVA Model

Step 6: Post-Hoc Tests (Tukey HSD)

Step 7: Calculate Effect Size

Step 8: Check Assumptions

Step 9: Report Results

PART B: YOUR TURN - INDEPENDENT PRACTICE

Practice Problem: Physical Activity and Depression

Research Question: Is there a difference in the number of days with poor mental health across three physical activity levels (None, Moderate, Vigorous)?

Step 1: Data Preparation

Step 2: Descriptive Statistics

Step 3: Visualization

Step 4: Set Up Hypotheses

Step 5: Fit the ANOVA Model

Step 6: Post-Hoc Tests

Step 7: Calculate Effect Size

Step 8: Check Assumptions

Step 9: Write Up Results

Your Results Section:

Results

Reflection Questions

Submission Checklist