Swing & A Deep Drive: Exploring the Statistician’s Heaven of the MLB
Aymane BELHADDAOUI
2026-08-10
Introduction: A Statistician’s Heaven
It is a well-worn cliché to say that Major League Baseball is a data analyst’s heaven—we have known that for decades. Instead, the motivation for this project stems from a specific, fascinating paradox in the sport: the reality that no one bats a thousand.
Think about the math for a second. In almost any other profession, succeeding 35% of the time is a failure. Yet, in baseball, when a player like Luis Arraez hits .354, it is considered an absolute masterclass in hitting. On the flip side, succeeding just 15% less — hitting .200 — drops a player right on the infamous “Mendoza Line,” fighting just to keep a job in the major leagues.
How can the margin between a perennial All-Star and a minor leaguer be so razor-thin? The answer lies in the extreme variance, luck, and probability inherent in the sport. When the gap between greatness and failure is that incredibly small, you cannot just rely on surface-level averages to tell the whole story. You need inferential statistics to separate the true signal from the noise.
To put my R skills to the test, I decided to dive into the famous Lahman Database — a comprehensive historical record of baseball statistics. In this notebook, I am walking through four fundamental statistical tests (T-Test, F-Test, ANOVA, and ANCOVA) to see what the data actually tells us about this game of inches. Finally, a special shoutout to the Toronto Blue Jays. Let’s go Blue Jays 🍁🫶!
1. The T-Test: Comparing Two Means
A two-sample t-test compares the averages of two distinct groups to see if they are significantly different from one another, or if the difference is just due to random chance.
The Question: Do Left-Handed Batters hit more Home Runs than Right-Handed Batters?
Note: We are looking at modern-era players (Year > 2000) with at least 300 At-Bats in a season.
batting_modern <- Batting %>%
filter(yearID > 2000, AB > 300) %>%
left_join(People %>% select(playerID, bats), by = "playerID") %>%
filter(bats %in% c("L", "R")) # Exclude switch hitters ("S") for a clean two-group test
# Run the T-Test
# Formula: Home Runs (HR) explained by Handedness (bats)
hr_ttest <- t.test(HR ~ bats, data = batting_modern)
hr_ttest##
## Welch Two Sample t-test
##
## data: HR by bats
## t = -0.87082, df = 3993.1, p-value = 0.3839
## alternative hypothesis: true difference in means between group L and group R is not equal to 0
## 95 percent confidence interval:
## -0.8387168 0.3228064
## sample estimates:
## mean in group L mean in group R
## 16.84227 17.10022The Verdict: Based on these results, there is no real difference in home run power between left-handed and right-handed batters. Both groups hit an average of about 17 home runs, and the high p-value (0.38) confirms the slight difference is just statistical noise. Essentially, it’s a tie.
2. The F-Test: Comparing Variances
While a t-test compares averages, an F-test compares variances (how spread out the data is). Averages lie. Two players can average 15 home runs a year, but one might hit exactly 15 every year (low variance), while the other hits 30 one year and 0 the next (high variance). The F-test helps analysts understand consistency, risk, and volatility.
The Question: Is the variance in stolen bases (SB) different between the American League (AL) and National League (NL)?
# 1. Prepare the data (Using modern era, players with at least 50 games)
speed_data <- Batting %>%
filter(yearID > 2010, G > 50, lgID %in% c("AL", "NL"))
# 2. Run the F-Test for variance
# Formula: Stolen Bases (SB) explained by League (lgID)
var_test_result <- var.test(SB ~ lgID, data = speed_data)
var_test_result##
## F test to compare two variances
##
## data: SB by lgID
## F = 0.93642, num df = 3601, denom df = 3784, p-value = 0.04613
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.8779158 0.9988728
## sample estimates:
## ratio of variances
## 0.9364189The Verdict: There’s a slight but real difference in the volatility of stolen bases between the two leagues. With a p-value just under 0.05, the test shows that the variation in stolen base numbers isn’t identical—specifically, the National League has slightly more spread or unpredictability in stolen bases than the American League. It’s a small difference, but statistically significant.
3. ANOVA: Analysis of Variance
ANOVA is essentially a t-test on steroids. Instead of comparing just two groups, it compares the means of three or more groups simultaneously. To compare 6 different divisions, running 15 separate t-tests is exhausting and mathematically dangerous. ANOVA gives you one clean “Yes” or “No” as to whether the group category matters at all.
The Question: Does run production (R) differ significantly across the baseball decades?
decades_data <- Teams %>%
filter(yearID >= 1980) %>%
mutate(Decade = case_when(
yearID < 1990 ~ "1980s",
yearID < 2000 ~ "1990s",
yearID < 2010 ~ "2000s",
yearID < 2020 ~ "2010s",
TRUE ~ "2020s"
))
# Run the ANOVA
# We wrap Decade in as.factor() so the subsequent Tukey test recognizes it as categorical
runs_anova <- aov(R ~ as.factor(Decade), data = decades_data)
summary(runs_anova)## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(Decade) 4 2244220 561055 47.88 <2e-16 ***
## Residuals 1313 15384224 11717
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The Verdict: Run scoring has definitely changed over the decades. With a p-value practically at zero, the results show a massive, undeniable difference in the average number of runs scored from the 1980s through the 2020s. Essentially, baseball has gone through distinct scoring eras—run production has been far from consistent over time.
3.5 The Follow-Up: Tukey HSD Test
Since our ANOVA told us that the Decade heavily influenced Run production, we need to know which decades were different. Was the Steroid Era (1990s/2000s) significantly different from the 1980s? We will run a Tukey Honest Significant Difference test to compare every decade head-to-head.
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = R ~ as.factor(Decade), data = decades_data)
##
## $`as.factor(Decade)`
## diff lwr upr p adj
## 1990s-1980s 54.12369 28.61430 79.633069 0.0000001
## 2000s-1980s 98.18359 73.13031 123.236870 0.0000000
## 2010s-1980s 38.30692 13.25364 63.360203 0.0003043
## 2020s-1980s -24.86752 -53.53710 3.802057 0.1244218
## 2000s-1990s 44.05990 19.44496 68.674844 0.0000112
## 2010s-1990s -15.81676 -40.43170 8.798177 0.4004921
## 2020s-1990s -78.99121 -107.27854 -50.703876 0.0000000
## 2010s-2000s -59.87667 -84.01861 -35.734722 0.0000000
## 2020s-2000s -123.05111 -150.92783 -95.174394 0.0000000
## 2020s-2010s -63.17444 -91.05116 -35.297727 0.0000000# The 'las = 1' argument flips the text labels so they are readable horizontally
par(mar = c(5, 6, 4, 2))
plot(tukey_results, las = 1, col = "#134A8E") 
The Verdict: This test pinpoints exactly when the
scoring shifts happened. The 2000s were the undisputed peak of baseball
offense, significantly outscoring every other decade. Since then,
scoring has plummeted—in fact, run production in the 2020s has dropped
so much that it is now statistically indistinguishable from the 1980s.
Let’s not forget however how the 2020 MLB season was shortened to 60
games due to Covid-19. —
4. ANCOVA: Analysis of Covariance
ANCOVA compares the means of different groups while controlling for a continuous third variable (a covariate) that might be unfairly influencing the results. Basically, It levels the playing field. It prevents you from making false conclusions based on unfair advantages.
The Question: Do American League (AL) teams score more runs than National League (NL) teams, assuming we control for the number of hits they get?
dh_era <- Teams %>%
filter(yearID >= 2010, yearID < 2020, lgID %in% c("AL", "NL"))
# Run the ANCOVA using a linear model (lm)
# We predict Runs (R) based on Hits (H) AND League (lgID)
ancova_model <- lm(R ~ H + lgID, data = dh_era)
anova(ancova_model)## Analysis of Variance Table
##
## Response: R
## Df Sum Sq Mean Sq F value Pr(>F)
## H 1 747143 747143 216.0134 < 2e-16 ***
## lgID 1 19584 19584 5.6622 0.01797 *
## Residuals 297 1027258 3459
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The Verdict:Unsurprisingly, getting hits is the biggest driver behind scoring runs. But the real takeaway here is the league effect: even if an American League and a National League team get the exact same number of hits, there is still a statistically significant difference in the number of runs they actually score. Because of the Designated Hitter rule during the 2010s which the NL only adopted in 2022, AL teams were simply more efficient at turning their hits into runs on the board.
5. Visuals
Raw console output is great for the analyst, but terrible for stakeholders. A massive part of data analysis is communication. Let’s filter down to the Toronto Blue Jays and show how to turn raw data into a beautiful, presentation-ready chart and table.
A. Using knitr::kable
Let’s look at the top 5 Blue Jays seasons by total Home Runs.
jays_top_hr <- Teams %>%
filter(teamID == "TOR") %>%
arrange(desc(HR)) %>%
head(5) %>%
select(Year = yearID, Wins = W, Losses = L, Runs = R, Hits = H, HomeRuns = HR)
kable(jays_top_hr, caption = "Top 5 Toronto Blue Jays Seasons by Home Runs")| Year | Wins | Losses | Runs | Hits | HomeRuns |
|---|---|---|---|---|---|
| 2021 | 91 | 71 | 846 | 1455 | 262 |
| 2010 | 85 | 77 | 755 | 1364 | 257 |
| 2019 | 67 | 95 | 726 | 1299 | 247 |
| 2000 | 83 | 79 | 861 | 1562 | 244 |
| 2015 | 93 | 69 | 891 | 1480 | 232 |
B. Using ggplot2
Let’s visualize the Blue Jays’ win trajectory over the last 20 years. We will use proper titles, minimal themes, and the official Blue Jays blue!
jays_modern <- Teams %>%
filter(teamID == "TOR", yearID >= 2000)
ggplot(jays_modern, aes(x = yearID, y = W)) +
geom_line(color = "#134A8E", linewidth = 1.2) +
geom_point(color = "#1D2D5C", size = 3) +
geom_hline(yintercept = 81, linetype = "dashed", color = "red", alpha = 0.6) +
labs(
title = "Toronto Blue Jays: Wins per Season (2000 - Present)",
subtitle = "Red dashed line indicates a .500 winning season",
x = "Season",
y = "Total Wins"
) +
theme_minimal() +
theme(
plot.title = element_text(face = "bold", size = 14),
axis.text = element_text(size = 10)
)
Conclusion
The MLB is proof that math doesn’t have to be dry. By applying
T-tests to bust myths, F-tests to find volatility, and ANOVA/ANCOVA to
dissect complex relationships, we can read the story behind the stats.
And with ggplot2 and knitr, we can tell that
story to anyone.
OK Blue Jays, Let’s Play Ball.