This lab introduces the workflow involved in ANOVA, starting with how ANOVA relates to our previous use of lm().

\(~\)

ANOVA

In our previous lab you worked with the professor salaries dataset:

fac_sal = read.csv("https://remiller1450.github.io/data/Salaries.csv")

Additionally, you’ve already seen how to fit the one-way ANOVA model using the lm() function:

rank_model = lm(salary ~ rank, data = fac_sal)
summary(rank_model)
## 
## Call:
## lm(formula = salary ~ rank, data = fac_sal)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -68972 -16376  -1580  11755 104773 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     93876       2954  31.777  < 2e-16 ***
## rankAsstProf   -13100       4131  -3.171  0.00164 ** 
## rankProf        32896       3290   9.997  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23630 on 394 degrees of freedom
## Multiple R-squared:  0.3943, Adjusted R-squared:  0.3912 
## F-statistic: 128.2 on 2 and 394 DF,  p-value: < 2.2e-16

At the bottom of the summary output, notice the line:

F-statistic: 128.2 on 2 and 394 DF, p-value: < 2.2e-16

This is the test statistic, degrees of freedom, and \(p\)-value of the ANOVA F-test.

While we can perform one-way ANOVA using lm(), it is often more useful to fit the one-way ANOVA model using the aov() function, which uses the exact same syntax of lm():

rank_anova = aov(salary ~ rank, data = fac_sal)
summary(rank_anova)
##              Df    Sum Sq   Mean Sq F value Pr(>F)    
## rank          2 1.432e+11 7.162e+10   128.2 <2e-16 ***
## Residuals   394 2.201e+11 5.586e+08                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Notice the same F-statistic (128.2), degrees of freedom (2 and 394), and \(p\)-value appear in the ANOVA table in this summary.

Question #1: In this question you’ll use the NFL player value data set, which was introduced in our lecture on correlation. The code below loads these data and converts the Rnd variable (the round in which a player was drafted) to a categorical data type.

## Load Data
draft_data = read.csv("https://remiller1450.github.io/data/draft_value.csv")

## Make the round the player was drafted categorical
draft_data$Rnd = factor(draft_data$Rnd)
  • Part A: Use the filter() function in the dplyr library to filter these data to only include the defensive positions, “DE”, “DT”, “LB”, and “CB”. Then use side-by-side boxplots to display the distribution of G (the number of games started by the player) for each of these defensive positions. Hint: Filtering was introduced in Lab 5.
  • Part B: Using the plot from Part A, do the data appear to exhibit more between-group or within-group variability? Briefly explain.
  • Part C: ANOVA assumes that errors are Normally distributed around each group mean, with a similar amount of variability across groups. Based upon your plot from Part A, do these assumptions seem reasonable for these data?
  • Part D: Use aov() to fit the ANOVA model using Pos (position) to predict G (number of games started), then use the summary() function to print the model’s ANOVA table. Provide a 1-sentence conclusion summarizing the results of the test.
  • Part E: Relying solely on side-by-side boxplots to check the assumptions of the ANOVA model can be problematic, particularly when some groups have relatively small sample sizes. Use the plot() function on your fitted aov model to evaluate the Normality assumption using a QQ-plot (the second plot), and use the command aggregate(G ~ Pos, data = your_filtered_data, FUN = sd) to find the standard deviations in each group. Briefly comment on whether the QQ-plot and group standard deviations support the assumptions of the ANOVA model.

\(~\)

Kruskal-Wallis

When the assumptions of the one-way ANOVA model are questionable, a common alternative is the Kruskal-Wallis test, which replaces the original values of the outcome variable with their ranks, then compares the distribution of ranks by group. This approach is conceptually similar to the Wilcoxon Rank-sum test introduced at the end of Lab 5.

The test is performed using the kruskal.test() function, which uses the same syntax as aov():

kruskal.test(salary ~ rank, data = fac_sal)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  salary by rank
## Kruskal-Wallis chi-squared = 194.01, df = 2, p-value < 2.2e-16

The hypotheses of the Kruskal-Wallis test are essentially the same as ANOVA: \[H_0: \text{all groups come from the same population distribution} \\ H_a: \text{at least one group differs}\]

Question #2: For this question you will continue using the NFL player value subset you created in Question 1, but with different variables to investigate whether starter seasons differ by draft round.

  • Part A: Create side-by-side boxplots showing the distribution of St (number of seasons as a primary starter) by Rnd. Briefly comment on whether the assumptions of ANOVA appear reasonable for these data.
  • Part B: Perform the Kruskal-Wallis test and provide a proper one-sentence conclusion that includes context and the strength of evidence.
  • Part C: Now perform one-way ANOVA on these data and report the \(p\)-value. How does it compare to the \(p\)-value of the Kruskal-Wallis test?

\(~\)

Post-hoc Testing

A key reason to use aov() is its compatibility with functions used for post-hoc testing, such as Tukey’s Honest Significant Difference (HSD) procedure, which helps identify which specific group means differ after the overall ANOVA finds statistical evidence of at least one difference. Tukey’s HSD procedure is demonstrated below:

TukeyHSD(rank_anova)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = salary ~ rank, data = fac_sal)
## 
## $rank
##                         diff       lwr       upr     p adj
## AsstProf-AssocProf -13100.45 -22818.71 -3382.195 0.0046514
## Prof-AssocProf      32895.67  25154.51 40636.836 0.0000000
## Prof-AsstProf       45996.12  38395.94 53596.307 0.0000000

This procedure is preferable to the \(T\)-tests we could obtain if we fit this model using lm() because it exhaustively covers every pairing of groups (rather than each group vs. the reference group) and it adjusts the \(p\)-values for multiple comparisons, ensuring a 5% family-wise error rate when these adjusted \(p\)-values are compared to \(\alpha = 0.05\).

In situations where we opt for the Kruskal-Wallis test after deciding the assumptions of ANOVA are questionable post-hoc testing can be done using the Dunn test, which is implemented in the dunn.test package via the dunn.test() function that is demonstrated below:

# install.packages("dunn.test")
library(dunn.test)

dunn.test(x = fac_sal$salary, g = fac_sal$rank, method="bonferroni")

Question #3: Answer these questions using the example provided in this section.

  • Part A: Based upon the Tukey pairwise test results, which faculty ranks have the largest difference means? Report the difference in means and the corresponding adjusted \(p\)-value.
  • Part B: Outside of this specific application, is it possible for a pairing of groups to have the largest observed differences in means but a non-significant adjusted \(p\)-value? Briefly explain your reasoning.
  • Part C: In your own words, what does it mean to control the family-wise error rate at 5% using adjusted \(p\)-values?

\(~\)

ANOVA Workflow

Summarizing the earlier sections of this lab and our most recent lecture, we can consider the complete ANOVA workflow:

  1. Visualize groups and preliminary check of assumptions.
  2. Perform a global test, either ANOVA or the Kruskal-Wallis, based upon the preliminary assumptions check.
  3. If one-way ANOVA was used, confirm the model assumptions using a QQ-plot and comparison of within group standard deviations.
  4. If the global test is significant, perform post-hoc pairwise testing to evaluate which groups are different.

When reporting a conclusion you should highlight both the results of the global test, as well as any notable pairwise differences, such as the “most significant” difference or comparisons involving the largest/highest interest groups.

Question #4: The code below loads data on the members of the 118th US Congress, which served until Jan 3, 2025, and creates a new variable Group that describes each member’s combination of chamber and political party.

  • Part A: Use side-by-side boxplots to compare the distribution of Age by Group. Do the assumptions of the one-way ANOVA model seem reasonable?
  • Part B: Based upon your assessment in Part A, perform the appropriate global hypothesis test.
  • Part C: If your global hypothesis test in Part B was statistically significant, perform post-hoc pairwise testing.
  • Part D: Summarize the findings of your global and pairwise test results.
  • Part E: In this question we used one-way ANOVA, which treated “House Democrat”, “House Republican”, “Senate Democrat”, and “Senate Republican” as four unrelated groups. A better alternative would have been to use a related approach known as two-way ANOVA to separately evaluate the effects of Chamber and Party. In your own words, why is it difficult to evaluate the effect of Party in our one-way ANOVA approach?
congress_data = read.csv("https://remiller1450.github.io/data/congress_2024.csv")
congress_data$Group = paste0(congress_data$Chamber, congress_data$Party)