20  Comparative Exploration

20.1 Learning Objectives

By the end of this chapter, you will be able to:

  • Frame a comparison as a null hypothesis vs. an alternative hypothesis and explain what a p-value does (and does not) mean.
  • Run an independent and paired t-test using scipy.stats, and compute Cohen’s d as an effect-size measure.
  • Run a one-way ANOVA to compare means across three or more groups, interpret the F-statistic, and compute eta-squared as an effect size.
  • Run a chi-square test of independence on categorical variables, and compute Cramér’s V.
  • Compute Pearson and Spearman correlation coefficients and interpret the difference between them.
  • Report statistical results with four elements: what you tested, the result, the test statistic + p-value, and the effect size + interpretation.
  • Recognize when statistical significance does and does not imply practical significance.

20.2 Why Comparative Exploration

Descriptive and qualitative exploration tell you what your data looks like and what it means in context. Comparative exploration answers the next question: when you compare two groups, three groups, or two variables, is the difference you see real, or could it have arisen by chance? This chapter introduces a small set of statistical comparisons — t-tests for two groups, ANOVA for three or more groups, chi-square for categorical contingency, and correlation for pairs of numeric variables — that cover the vast majority of the comparisons an analyst makes with tabular data.

The framework is deliberately simple. Each test answers a version of the same question: is the pattern I see distinguishable from random variation? You’ll read each result the same way: name the test, report the statistic and p-value, compute an effect size, and interpret the magnitude in plain language. Once you can read one test’s output, you can read them all.

If you have never worked with p-values before, they can feel counterintuitive at first — the interpretation flips in a way most introductions gloss over. Re-read the interpretation paragraphs slowly, and come back to them after running a few tests yourself. The intuition settles with practice.

20.3 Setup

This chapter uses pandas, numpy, matplotlib, seaborn, and scipy.stats.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats

pd.set_option('display.max_columns', None)
pd.set_option('display.float_format', '{:.3f}'.format)
sns.set_style('whitegrid')

Install missing packages with conda install pandas numpy matplotlib seaborn scipy (or pip install pandas numpy matplotlib seaborn scipy) and restart the Jupyter kernel. scipy is the most likely miss if you are coming from a minimal pandas environment. See Chapter 1.

20.4 Datasets

# Lahman batting — qualified player-seasons with salary
batting = pd.read_csv('Batting.csv')
batting['BA'] = batting['H'] / batting['AB'].replace(0, float('nan'))
batting['TB'] = batting['H'] + 2*batting['2B'] + 3*batting['3B'] + 4*batting['HR']
batting['SLG'] = batting['TB'] / batting['AB'].replace(0, float('nan'))
batting['HBP'] = batting['HBP'].fillna(0)
batting['SF'] = batting['SF'].fillna(0)
batting['OBP'] = (
    (batting['H'] + batting['BB'] + batting['HBP'])
    / (batting['AB'] + batting['BB'] + batting['HBP'] + batting['SF']).replace(0, float('nan'))
)
batting['OPS'] = batting['OBP'] + batting['SLG']
batting['decade'] = (batting['yearID'] // 10) * 10
qualified = batting[batting['AB'] >= 100].copy()

salaries = pd.read_csv('Salaries.csv')
bat_sal = pd.merge(
    qualified, salaries[['playerID', 'yearID', 'salary']],
    on=['playerID', 'yearID'], how='inner'
)
# Cannabis sales
cannabis = pd.read_csv(
    'colorado_cannabis_sales_by_county.csv',
    header=2, skipfooter=8, engine='python', na_values=['NR']
)
cannabis.columns = ['month', 'year', 'county', 'medical', 'retail']
cannabis['county'] = cannabis['county'].str.strip()
for col in ['medical', 'retail']:
    cannabis[col] = pd.to_numeric(
        cannabis[col].str.replace('$', '', regex=False).str.replace(',', '', regex=False),
        errors='coerce'
    )
cannabis = cannabis[cannabis['county'] != 'Total']

20.5 The Logic of Hypothesis Testing

Every hypothesis test follows the same structure:

  1. State the null hypothesis (H₀): There is no difference between the groups, or there is no relationship between the variables. The observed pattern is due to chance.
  2. State the alternative hypothesis (H₁): There is a real difference or relationship.
  3. Compute a test statistic that measures how far the observed data deviates from what the null hypothesis would predict.
  4. Compute a p-value — the probability of seeing a test statistic this extreme (or more extreme) if the null hypothesis were true.
  5. Interpret the result. A small p-value (conventionally below 0.05) means the observed pattern is unlikely to have arisen by chance alone. A large p-value means you do not have enough evidence to reject the null hypothesis.

A p-value does not tell you the probability that the null hypothesis is true. It tells you how surprising your data would be if the null hypothesis were true. This distinction matters: a p-value of 0.03 means “if there were truly no difference, we’d see data this extreme only 3% of the time” — not “there is a 97% chance the difference is real.”

WarningStatistical significance is not practical significance

A large dataset can produce a tiny p-value for a difference that is too small to matter. Always report and interpret the effect size — the magnitude of the difference — alongside the p-value. A statistically significant but tiny difference is not analytically useful.

20.6 Comparing Two Groups: The t-Test

The independent samples t-test asks: Is the mean of one group different from the mean of another group? Use it when you have one numeric outcome and two categorical groups.

20.6.1 Example: Batting average before and after the steroid era

The late 1990s through early 2000s are sometimes called baseball’s “steroid era.” Did batting averages differ between the pre-steroid period (1980–1993) and the steroid era (1994–2004)?

pre_steroid = qualified[
    (qualified['yearID'] >= 1980) & (qualified['yearID'] <= 1993)
]['BA'].dropna()

steroid_era = qualified[
    (qualified['yearID'] >= 1994) & (qualified['yearID'] <= 2004)
]['BA'].dropna()

print(f'Pre-steroid: n={len(pre_steroid)}, mean={pre_steroid.mean():.4f}, std={pre_steroid.std():.4f}')
print(f'Steroid era: n={len(steroid_era)}, mean={steroid_era.mean():.4f}, std={steroid_era.std():.4f}')

Visualize the distributions before testing:

fig, ax = plt.subplots(figsize=(8, 4))
sns.kdeplot(pre_steroid, label='Pre-steroid (1980–93)', fill=True, alpha=0.3, ax=ax)
sns.kdeplot(steroid_era, label='Steroid era (1994–2004)', fill=True, alpha=0.3, ax=ax)
ax.set_xlabel('Batting Average')
ax.set_title('Batting Average Distribution: Pre-Steroid vs. Steroid Era')
ax.legend()
plt.tight_layout()
plt.show()

Now run the t-test:

t_stat, p_value = stats.ttest_ind(pre_steroid, steroid_era)

print(f't-statistic: {t_stat:.3f}')
print(f'p-value: {p_value:.4f}')

20.6.2 Interpreting the result

The mean difference between the two eras can be computed directly:

diff = steroid_era.mean() - pre_steroid.mean()
print(f'Mean difference: {diff:.4f}')

Here is what you can conclude: if the p-value is below 0.05, the difference in mean batting average between the two eras is statistically significant — unlikely to have arisen by chance. The effect size (the mean difference) tells you whether the difference is large enough to be practically meaningful. A difference of 0.005 in batting average is detectable with thousands of player-seasons but may not be substantively interesting. A difference of 0.015 would be more notable.

What you cannot conclude: that steroids caused the difference. These are different groups of players measured in different years with different rules, ballparks, and training methods. The comparison is between-subjects and observational — all the caveats from Chapter 3 apply.

20.6.3 Effect size: Cohen’s d

Cohen’s d measures the difference in means in terms of standard deviations. It provides a standardized measure of effect size that is comparable across different variables and studies.

pooled_std = np.sqrt((pre_steroid.std()**2 + steroid_era.std()**2) / 2)
cohens_d = (steroid_era.mean() - pre_steroid.mean()) / pooled_std

print(f"Cohen's d: {cohens_d:.3f}")

Rules of thumb for interpreting Cohen’s d: 0.2 is a small effect, 0.5 is medium, and 0.8 is large. These are guidelines, not thresholds — context determines what counts as meaningful in your domain.

20.8 Comparing Three or More Groups: One-Way ANOVA

A t-test compares two groups. When you have three or more groups and you want to ask whether any of them differ from the others, use a one-way ANOVA (Analysis of Variance). The test produces an F-statistic — the ratio of between-group variance to within-group variance. A large F (and correspondingly small p-value) means the group means are farther apart than you would expect under the null hypothesis that every group was drawn from the same distribution.

20.8.1 Example: Batting average across decades

Did mean batting average differ across the 1980s, 1990s, 2000s, and 2010s?

groups = {}
for decade in [1980, 1990, 2000, 2010]:
    groups[decade] = qualified[qualified['decade'] == decade]['BA'].dropna()
    print(f"{decade}s: n={len(groups[decade])}, mean={groups[decade].mean():.4f}")

f_stat, p_value = stats.f_oneway(*groups.values())
print(f'\nF-statistic: {f_stat:.2f}')
print(f'p-value: {p_value:.6f}')

Here is what you can conclude: a small p-value means at least one decade’s mean batting average differs from the others. ANOVA does not tell you which. For that, you need a post-hoc test — conventionally Tukey’s HSD (Honestly Significant Difference), which compares every pair while controlling the family-wise error rate.

from statsmodels.stats.multicomp import pairwise_tukeyhsd

decade_ba = qualified[qualified['decade'].isin([1980, 1990, 2000, 2010])].dropna(subset=['BA'])
tukey = pairwise_tukeyhsd(
    endog=decade_ba['BA'],
    groups=decade_ba['decade'].astype(int),
    alpha=0.05,
)
print(tukey.summary())

The Tukey output lists every pair with an adjusted p-value and whether the difference is significant at the specified alpha. This is the table to report alongside the ANOVA.

20.8.2 Effect size: eta-squared

The F-statistic itself depends heavily on sample size, so for effect size you use eta-squared (η²) — the fraction of total variance explained by group membership. Interpretively, it’s the ANOVA analog of R² in regression.

# Manual eta-squared from the ANOVA decomposition
grand_mean = decade_ba['BA'].mean()
ss_total = ((decade_ba['BA'] - grand_mean) ** 2).sum()
ss_between = sum(
    len(g) * (g.mean() - grand_mean) ** 2 for g in groups.values()
)
eta_sq = ss_between / ss_total
print(f"η² = {eta_sq:.3f}")

Conventional interpretation: η² ≈ 0.01 is small, 0.06 is medium, 0.14 is large. In large samples, even tiny η² can produce a significant F — report η² explicitly and interpret its magnitude.

20.8.3 Assumptions and robustness

ANOVA assumes roughly normal distributions within each group and roughly equal variances across groups. Both are often violated in real data. If the group sizes are reasonably balanced, the test is fairly robust. If variances differ substantially, use Welch’s ANOVA (stats.alexandergovern or, more commonly, scipy.stats.f_oneway with data that has been checked via Levene’s test first). For small samples with non-normal data, use the non-parametric Kruskal–Wallis test (stats.kruskal) as a robust alternative.

20.9 Comparing Categorical Variables: The Chi-Square Test

The chi-square test of independence asks: Are two categorical variables associated, or are they independent? Use it when both your variables are categorical.

20.9.1 Example: Has the distribution of home run tiers changed across eras?

# Bin HR into categories
hr_data = qualified[qualified['decade'].isin([1980, 2000])].copy()
hr_data['hr_tier'] = pd.cut(
    hr_data['HR'], bins=[-1, 5, 15, 30, 100],
    labels=['0-5', '6-15', '16-30', '31+']
)
hr_data['era'] = hr_data['decade'].map({1980: '1980s', 2000: '2000s'})

# Build contingency table
contingency = pd.crosstab(hr_data['era'], hr_data['hr_tier'])
print('Contingency table:')
print(contingency)
chi2, p_value, dof, expected = stats.chi2_contingency(contingency)

print(f'\nChi-square statistic: {chi2:.2f}')
print(f'Degrees of freedom: {dof}')
print(f'p-value: {p_value:.6f}')

Here is what you can conclude: if the p-value is small, the distribution of home run tiers differs significantly between the two eras — the two variables (era and HR tier) are not independent. To understand how they differ, compare the observed counts to the expected counts (what you would see if era and HR tier were independent):

expected_df = pd.DataFrame(
    expected, index=contingency.index, columns=contingency.columns
).round(1)

print('Expected counts (if independent):')
print(expected_df)
print('\nObserved - Expected:')
print((contingency - expected_df).round(1))

Positive residuals indicate cells with more observations than expected; negative residuals indicate fewer. This tells you which specific categories drive the overall association.

20.9.2 Effect size: Cramér’s V

Cramér’s V is the chi-square analog of Cohen’s d — a measure of how strongly the two categorical variables are associated, scaled from 0 (independent) to 1 (perfectly associated):

n = contingency.sum().sum()
min_dim = min(contingency.shape) - 1
cramers_v = np.sqrt(chi2 / (n * min_dim))

print(f"Cramér's V: {cramers_v:.3f}")

20.10 Measuring Linear Relationships: Correlation

The Pearson correlation coefficient (r) measures the strength and direction of a linear relationship between two numeric variables. It ranges from -1 (perfect negative linear relationship) to +1 (perfect positive linear relationship), with 0 indicating no linear relationship.

20.10.1 Example: Salary and performance

# Filter to a single recent season with enough data
season_2015 = bat_sal[bat_sal['yearID'] == 2015].dropna(subset=['salary', 'OPS'])

r, p_value = stats.pearsonr(season_2015['salary'], season_2015['OPS'])

print(f'Pearson r: {r:.3f}')
print(f'p-value: {p_value:.4f}')
print(f'n: {len(season_2015)}')

Visualize alongside the statistic:

fig, ax = plt.subplots(figsize=(8, 5))
ax.scatter(season_2015['salary'] / 1e6, season_2015['OPS'], alpha=0.4, s=20)
ax.set_xlabel('Salary (millions $)')
ax.set_ylabel('OPS')
ax.set_title(f'Salary vs. OPS (2015) — r = {r:.3f}, p = {p_value:.4f}')
plt.tight_layout()
plt.show()

Here is what you can conclude: the correlation coefficient tells you the strength and direction of the linear relationship. A positive r means higher salary is associated with higher OPS. The p-value tells you whether this correlation is statistically significant. The magnitude of r tells you how tightly the points cluster around a line — an r of 0.3 indicates a weak-to-moderate positive relationship, meaning salary explains about 9% (r² = 0.09) of the variation in OPS.

What you cannot conclude: that paying a player more would improve their performance. As discussed in Chapter 3, the causal arrow is ambiguous, and the relationship is confounded by age, experience, and roster construction.

20.10.2 Spearman correlation for non-linear monotonic relationships

Pearson’s r measures linear relationships. If the relationship is monotonic but curved (both variables increase together, but not in a straight line), Spearman’s rank correlation is more appropriate:

rho, p_value = stats.spearmanr(season_2015['salary'], season_2015['OPS'])

print(f'Spearman rho: {rho:.3f}')
print(f'p-value: {p_value:.4f}')

Spearman’s ρ (rho) works by ranking both variables and computing the Pearson correlation on the ranks. It is robust to outliers and does not assume the relationship is linear.

20.11 Correlation Across Many Variables

When you have several numeric columns, compute the full pairwise correlation matrix and identify which pairs have the strongest relationships:

metrics = bat_sal[['AB', 'HR', 'BB', 'BA', 'SLG', 'OBP', 'OPS', 'salary']].dropna()
corr = metrics.corr()

fig, ax = plt.subplots(figsize=(8, 6))
sns.heatmap(corr, annot=True, fmt='.2f', cmap='coolwarm', center=0,
            square=True, linewidths=0.5, ax=ax)
ax.set_title('Correlation Matrix — Batting Metrics and Salary')
plt.tight_layout()
plt.show()

Look for the cells with the strongest colors. Correlations between derived metrics (like OBP and OPS) will be high by construction. The interesting correlations are between conceptually distinct variables — salary and performance metrics, counting stats and rate stats.

20.12 Summary: Choosing a Test

Question Test Input
Is the mean different between two independent groups? Independent t-test (stats.ttest_ind) One numeric + one binary categorical
Is there a systematic difference between two paired measurements? Paired t-test (stats.ttest_rel) Two numeric columns from the same entity
Is the mean different across three or more groups? One-way ANOVA (stats.f_oneway), followed by Tukey HSD for pairwise comparisons One numeric + one categorical with 3+ levels
Are two categorical variables associated? Chi-square test (stats.chi2_contingency) Two categorical variables (contingency table)
How strong is the linear relationship between two variables? Pearson correlation (stats.pearsonr) Two numeric variables
How strong is the monotonic (but possibly non-linear) relationship? Spearman correlation (stats.spearmanr) Two numeric variables

20.13 Reporting Results

When reporting the result of a statistical test, include four elements:

  1. What you tested: “We compared mean batting average between the pre-steroid (1980–93) and steroid (1994–2004) eras.”
  2. The result: “The steroid era had a slightly higher mean batting average (0.271 vs. 0.267).”
  3. The test statistic and p-value: “t(8452) = -3.21, p = 0.001”
  4. The effect size and interpretation: “Cohen’s d = 0.07, a small effect. While the difference is statistically significant due to the large sample size, the practical difference of 4 points of batting average is modest.”

This pattern — what, result, statistics, interpretation — keeps your reporting honest and complete. It prevents the common mistake of reporting a p-value without context about the magnitude of the effect.

20.14 Common Issues in This Chapter

WarningYour p-value is tiny but the effect is invisible

What it means: With a large sample, almost any non-zero difference is “statistically significant.” This is a feature of the test, not evidence that the difference matters.

How to fix it: Always report effect size (Cohen’s d, Cramér’s V, r) alongside the p-value. Interpret the magnitude in plain language — “a 4-point difference in batting average” — not the p-value in isolation.

Warningttest_ind complains about unequal variances or NaN values

What it means: By default ttest_ind assumes equal variances and breaks on NaN. If one group has much more spread than the other, or if either input has missing values, you will get misleading results or errors.

How to fix it: Drop NaNs with .dropna() before passing to the test. If variances differ substantially, use stats.ttest_ind(a, b, equal_var=False) to run Welch’s t-test.

WarningChi-square warns about low expected counts

What it means: The chi-square approximation breaks down when any expected cell count is below 5. The test still runs but the p-value may be unreliable.

How to fix it: Collapse rare categories before testing, or use stats.fisher_exact for 2×2 tables. The expected array returned by chi2_contingency tells you which cells are problematic.

WarningPearson r is near zero but the scatter plot shows a clear pattern

What it means: Pearson measures linear relationships. A U-shaped or curved pattern can have near-zero Pearson correlation but obvious structure.

How to fix it: Visualize first, always. If the pattern is monotonic but curved, switch to Spearman. If the pattern is not monotonic, correlation is not the right tool — consider binning the predictor and comparing means across bins.

20.15 Key Takeaways

  • Every test follows the same structure: state the hypotheses, compute the statistic, interpret the p-value and the effect size.
  • A p-value tells you how surprising your data would be under the null hypothesis — not the probability that the null is true.
  • Statistical significance ≠ practical significance. Always report effect sizes. Tiny differences can be “significant” with a large enough sample.
  • Pearson for linear, Spearman for monotonic-but-curved, t-test for two groups, chi-square for two categorical variables, paired t-test for repeated measurements on the same entity.
  • Visualize before you test. A scatter plot or histogram catches the problems a p-value cannot.

20.16 Further Reading

  • scipy.stats reference — docs.scipy.org/doc/scipy/reference/stats.html covers every test introduced in this chapter and many more.
  • scipy.stats tutorial — docs.scipy.org/doc/scipy/tutorial/stats.html walks through distributions, tests, and descriptive statistics.
  • “Statistical Thinking for the 21st Century” by Russell Poldrack — free online at statsthinking21.org. A modern introduction that emphasizes interpretation over mechanics.
  • ASA Statement on p-values — amstat.org is the American Statistical Association’s short, authoritative statement on what p-values mean and do not mean. Two pages; worth a careful read.
  • statsmodels — For regression, ANOVA, and generalized linear models beyond the tests covered here, see statsmodels. Its formula-based interface will feel familiar if you have used R or Stata.
  • scikit-learn — When your question shifts from inference (what explains my outcome?) to prediction (what outcome should I expect?), scikit-learn provides a consistent API for classical machine learning algorithms.

20.17 Exercises

Compare the mean number of home runs per qualified player-season between the 1980s and the 2010s. Run an independent t-test and compute Cohen’s d. Is the difference statistically significant? Is it practically meaningful?

hr_80s = qualified[qualified['decade'] == 1980]['HR'].dropna()
hr_10s = qualified[qualified['decade'] == 2010]['HR'].dropna()

print(f'1980s: n={len(hr_80s)}, mean={hr_80s.mean():.2f}')
print(f'2010s: n={len(hr_10s)}, mean={hr_10s.mean():.2f}')

t, p = stats.ttest_ind(hr_80s, hr_10s)
pooled = np.sqrt((hr_80s.std()**2 + hr_10s.std()**2) / 2)
d = (hr_10s.mean() - hr_80s.mean()) / pooled

print(f'\nt = {t:.3f}, p = {p:.4f}')
print(f"Cohen's d = {d:.3f}")

Create a contingency table of Colorado counties by year (pick two years at least 5 years apart) and whether they reported any retail sales (yes/no). Run a chi-square test. Has the proportion of counties with retail sales changed significantly?

for_chi = cannabis[cannabis['year'].isin([2016, 2023])].copy()
for_chi['has_retail'] = for_chi['retail'].notna() & (for_chi['retail'] > 0)
yearly = for_chi.groupby(['year', 'county'])['has_retail'].any().reset_index()

ct = pd.crosstab(yearly['year'], yearly['has_retail'])
ct.columns = ['No Retail', 'Has Retail']
print(ct)

chi2, p, dof, exp = stats.chi2_contingency(ct)
print(f'\nChi-square = {chi2:.2f}, p = {p:.4f}')

Compute the Pearson and Spearman correlations between at-bats (AB) and salary for the 2010s. Which is larger, and what does the difference tell you about the shape of the relationship? Produce a scatter plot to visually confirm your interpretation.

recent = bat_sal[bat_sal['decade'] == 2010].dropna(subset=['AB', 'salary'])

r_pearson, p_pearson = stats.pearsonr(recent['AB'], recent['salary'])
r_spearman, p_spearman = stats.spearmanr(recent['AB'], recent['salary'])

print(f'Pearson r = {r_pearson:.3f} (p = {p_pearson:.4f})')
print(f'Spearman ρ = {r_spearman:.3f} (p = {p_spearman:.4f})')

fig, ax = plt.subplots(figsize=(8, 5))
ax.scatter(recent['AB'], recent['salary']/1e6, alpha=0.3, s=15)
ax.set_xlabel('At Bats')
ax.set_ylabel('Salary (millions $)')
ax.set_title(f'AB vs. Salary — Pearson r={r_pearson:.3f}, Spearman ρ={r_spearman:.3f}')
plt.tight_layout()
plt.show()

If Spearman is larger than Pearson, the relationship is monotonic but nonlinear — likely because salary increases steeply at high AB counts but flattens at low counts.