09 - Chi-Square Test: Independence and Goodness-of-Fit

An A/B test gave 120 conversions out of 1000 visits for version A and 150 out of 1000 for B: a relative increase of 25%, or just luck? The data are not measurements but counts in categories (converted yes or no, version A or B), and the t-test from the previous article does not apply. The right test is the chi-square (written χ²), and on your table it gives p = 0.0496 without correction and p = 0.0577 with Yates’ continuity correction: a verdict that depends on a choice that textbooks present as obvious.

This article, the ninth in the Probability and Statistics series, explains what the test computes, when the approximation it rests on holds, and what to do when the data are few. You need the analysis of variance and the p-value from the article on hypothesis tests.

What the χ² statistic computes

You have k categories. For each one you count how many observations there are (O_i, the observed frequencies) and you compute how many you would expect if the null hypothesis were true (E_i, the expected frequencies). Pearson’s statistic is χ² = Σ(O_i − E_i)²/E_i. Each term is the squared deviation, divided by what you expected: a deviation of 10 on a cell with 1000 expected weighs little, the same deviation on a cell with 20 expected weighs a lot. If the statistic is small, the data resemble the model; if it is large, they do not. The definition is in the NIST statistics handbook (section 1.3.5.15).

The reference distribution is the chi-square with d degrees of freedom: the law of the sum of the squares of d independent standard Normal variables, with mean d and variance 2d. The number of degrees of freedom counts how many free quantities there are in the comparison.

With k categories and the total fixed, only k − 1 counts are free: hence d = k − 1 in the goodness-of-fit test. The reason why the statistic has this law is an approximation, valid when the counts are large (the central limit theorem applied to the counts). This is why the test is “asymptotic”: the rest of the article measures how well it holds up with small numbers.

Goodness of fit: is a die fair?

You roll a die 60 times and get faces 1-6 respectively 8, 12, 9, 15, 6 and 10 times. If the die is fair you expect 10 per face, so χ² = ((−2)² + 2² + (−1)² + 5² + (−4)² + 0²)/10 = (4 + 4 + 1 + 25 + 16 + 0)/10 = 5.0, with 6 − 1 = 5 degrees of freedom. The p-value is the probability that a chi-square with 5 degrees of freedom exceeds 5.0.

The program does the calculation and compares it with scipy.stats.chisquare. Then it measures two things. The first: the p-value obtained by simulating one million fair dice rolled 60 times, which does not depend on the approximation. The second: how far off the approximation is with little data, on a case with four categories of probability 0.90, 0.05, 0.03 and 0.02. With 4 categories you can enumerate all the possible outcomes of the rolls and add up their exact probabilities, so the “true size” of the test at 5% is exact, not simulated.

# Goodness of fit: is a die fair? Worked exercise, then the error of the chi-square approximation
import numpy as np
from scipy import stats

osservate = np.array([8, 12, 9, 15, 6, 10])          # 60 rolls
n, k = osservate.sum(), len(osservate)
attese = np.full(k, n / k)                           # 10 per face if the die is fair
chi2 = ((osservate - attese) ** 2 / attese).sum()
p = stats.chi2.sf(chi2, df=k - 1)
r = stats.chisquare(osservate, attese)
print(f"n = {n}, expected = {attese[0]:.0f}, chi2 = {chi2:.3f}, degrees of freedom = {k - 1}, p = {p:.4f}")
print(f"scipy.stats.chisquare: chi2 = {r.statistic:.3f}, p = {r.pvalue:.4f}")

# "Almost exact" p-value: fraction of 1,000,000 fair dice, rolled 60 times, with chi2 at least this large
rng = np.random.default_rng(7)
sim = rng.multinomial(n, np.full(k, 1 / k), size=1_000_000)
chi2_sim = ((sim - attese) ** 2 / attese).sum(axis=1)
print(f"p-value by simulation = {(chi2_sim >= chi2 - 1e-9).mean():.4f} +/- 0.0005 (chi-square approximation: {p:.4f})")

# True size of the 5% test with 4 categories of probability (0.90, 0.05, 0.03, 0.02):
# I enumerate ALL the multinomial outcomes, so the result is exact
prob = np.array([0.90, 0.05, 0.03, 0.02])
crit = stats.chi2.ppf(0.95, df=3)
print("   n  minimum expected  exact size of the 5% test")
for n in (10, 20, 50, 100, 200):
    a, b, c = np.meshgrid(np.arange(n + 1), np.arange(n + 1), np.arange(n + 1), indexing="ij")
    ok = a + b + c <= n
    x = np.stack([a[ok], b[ok], c[ok], n - a[ok] - b[ok] - c[ok]], axis=1)
    pmf = stats.multinomial.pmf(x, n, prob)
    att = n * prob
    stat = ((x - att) ** 2 / att).sum(axis=1)
    print(f"{n:4d}  {att.min():16.2f}  {pmf[stat >= crit].sum():.4f}   (sum of probabilities = {pmf.sum():.6f})")

# Output:
# n = 60, expected = 10, chi2 = 5.000, degrees of freedom = 5, p = 0.4159
# scipy.stats.chisquare: chi2 = 5.000, p = 0.4159
# p-value by simulation = 0.4306 +/- 0.0005 (chi-square approximation: 0.4159)
#    n  minimum expected  exact size of the 5% test
#   10              0.20  0.0721   (sum of probabilities = 1.000000)
#   20              0.40  0.0750   (sum of probabilities = 1.000000)
#   50              1.00  0.0556   (sum of probabilities = 1.000000)
#  100              2.00  0.0515   (sum of probabilities = 1.000000)
#  200              4.00  0.0491   (sum of probabilities = 1.000000)

The calculation by hand coincides with SciPy: χ² = 5.000, p = 0.4159, and so the die gives no evidence of being loaded. The simulation gives 0.4306 (standard error 0.0005): more than the approximation, because the statistic takes few distinct values and “at least 5.0” also counts the rolls that give exactly 5.0.

In the table of the true size, with minimum expected 0.20 the “5%” test actually rejects 7.21% of the time, with 0.40 7.50%; with 1.00 5.56%, with 2.00 5.15% and with 4.00 4.91%. The practical threshold you find in textbooks (expected frequencies of at least 5; the NIST handbook and the SciPy documentation report it) is a prudent rule, not a sharp boundary: the damage grows gradually as the expected values fall below 2.

Parameters estimated from the data: the degrees of freedom drop

Often the null hypothesis is not “the die is fair” but “the data are Poisson”, without saying with which parameter. You have to estimate the parameter from the data, and this has a cost: the model has already adapted to the counts, so the deviations are smaller than those of a fixed model. The correct theorem, reported in the goodness-of-fit testing lecture of the MIT 18.650 course, says that the statistic tends to a chi-square with K − d − 1 degrees of freedom, where d is the number of estimated parameters. Each estimated parameter removes one degree of freedom.

The program verifies it: it generates 20,000 samples of 200 Poisson observations with parameter 3, groups them into 8 categories (0, 1, …, 6 and “7 or more”), estimates the parameter by maximum likelihood on the grouped counts (which is the estimate for which the theorem holds), computes χ² and counts how many times it exceeds the 5% threshold with K − 1 = 7 degrees of freedom and with K − 2 = 6.

# If you estimate a parameter from the data, the degrees of freedom drop: k - 1 - (estimated parameters)
import numpy as np
from scipy import stats, optimize

lam_vero, n, R = 3.0, 200, 20_000
k = 8                                          # categories: 0, 1, ..., 6 and "7 or more"
rng = np.random.default_rng(12)

def probabilita(lam):
    p = stats.poisson.pmf(np.arange(k - 1), lam)
    return np.append(p, 1 - p.sum())           # the last category collects the tail

def stima_ml_raggruppata(conteggi):
    """Maximum likelihood of lambda on the grouped counts (multinomial likelihood)."""
    f = lambda lam: -(conteggi * np.log(probabilita(lam))).sum()
    return optimize.minimize_scalar(f, bounds=(0.1, 10), method="bounded", options={"xatol": 1e-9}).x

chi2 = np.empty(R)
for r in range(R):
    x = np.minimum(rng.poisson(lam_vero, size=n), k - 1)       # data already grouped
    o = np.bincount(x, minlength=k)
    lam = stima_ml_raggruppata(o)
    e = n * probabilita(lam)
    chi2[r] = ((o - e) ** 2 / e).sum()

print(f"minimum expected with true lambda = {n * probabilita(lam_vero).min():.2f}; repetitions = {R}")
print("degrees of freedom  5% threshold  rejections with that threshold (H0 true)")
for df in (k - 1, k - 2):
    soglia = stats.chi2.ppf(0.95, df)
    print(f"{df:18d}  {soglia:13.3f}  {(chi2 > soglia).mean():.4f}")
print(f"simulated mean of chi2 = {chi2.mean():.3f} (k-1 = {k - 1}, k-2 = {k - 2}); standard error of the level: {np.sqrt(0.05 * 0.95 / R):.4f}")

# Output:
# minimum expected with true lambda = 6.70; repetitions = 20000
# degrees of freedom  5% threshold  rejections with that threshold (H0 true)
#                  7         14.067  0.0303
#                  6         12.592  0.0478
# simulated mean of chi2 = 5.987 (k-1 = 7, k-2 = 6); standard error of the level: 0.0015

With 7 degrees of freedom the test rejects 3.03% of the time, that is, it is too prudent: it loses power without telling you. With 6 it rejects 4.78%, close to the declared 5% (standard error 0.0015). The simulated mean of the statistic, 5.987, is practically 6, the expected value of a chi-square with 6 degrees of freedom, and not 7. The most common mistake is to compute the parameter from the data and then use K − 1 anyway: the test is too prudent and accepts wrong models more easily. The estimate must be made on the grouped counts; with the estimate on the raw data the limiting law is no longer exactly that.

Independence in a contingency table

A contingency table with r rows and c columns contains the counts for every combination of two categorical variables. The null hypothesis is that the variables are independent: in your A/B test, that the probability of converting does not depend on the version. Under independence, the expected frequency of cell (i, j) is E_ij = (row total i) × (column total j)/n, because the row and column probabilities multiply.

The degrees of freedom are (r − 1)(c − 1): the rc cells minus the total, minus the r − 1 row probabilities and the c − 1 column ones that you estimated from the margins. In your case, (2 − 1)(2 − 1) = 1. With two samples already separated (the two versions) one speaks of a test of homogeneity, but the statistic is the same.

The program applies the formula to the A/B test table, compares it with chi2_contingency with and without Yates’ correction (which the function applies only with one degree of freedom), shows that for a 2 × 2 the test coincides with the comparison of two proportions (z² = χ²), and measures the effect size with three indices: Cramér’s V (√(χ²/(n(min(r, c) − 1))), between 0 and 1), the odds ratio with the Woolf interval and the difference between the two proportions with its interval. Finally it verifies the degrees of freedom on a 3 × 4 table with true independence, simulating 20,000 tables.

# Independence test on a 2 x 2 table (outcome of an A/B test), by hand and with scipy
import numpy as np
from scipy import stats

tabella = np.array([[120, 880],      # version A: 120 conversions out of 1000
                    [150, 850]])     # version B: 150 conversions out of 1000
n = tabella.sum()
righe, colonne = tabella.sum(axis=1), tabella.sum(axis=0)
attese = np.outer(righe, colonne) / n                 # E_ij = (row total) * (column total) / n
chi2 = ((tabella - attese) ** 2 / attese).sum()
gdl = (tabella.shape[0] - 1) * (tabella.shape[1] - 1)
print(f"expected = {attese.tolist()}")
print(f"chi2 = {chi2:.4f}, degrees of freedom = {gdl}, p = {stats.chi2.sf(chi2, gdl):.4f}")
c, p, d, _ = stats.chi2_contingency(tabella, correction=False)
print(f"scipy (without Yates): chi2 = {c:.4f}, p = {p:.4f}")
cy, py_, _, _ = stats.chi2_contingency(tabella, correction=True)
print(f"scipy (with Yates):    chi2 = {cy:.4f}, p = {py_:.4f}")

# The same test as a comparison of two proportions: z^2 = chi2
pA, pB, pc = 120 / 1000, 150 / 1000, 270 / 2000
z = (pB - pA) / np.sqrt(pc * (1 - pc) * (1 / 1000 + 1 / 1000))
print(f"z = {z:.4f}, z^2 = {z ** 2:.4f}")

# How large the effect is: Cramer's V, odds ratio with the Woolf interval, difference between proportions
V = np.sqrt(chi2 / (n * (min(tabella.shape) - 1)))
OR = tabella[1, 0] * tabella[0, 1] / (tabella[0, 0] * tabella[1, 1])
se_log = np.sqrt((1 / tabella).sum())
lo, hi = np.exp(np.log(OR) + np.array([-1.96, 1.96]) * se_log)
se_d = np.sqrt(pA * (1 - pA) / 1000 + pB * (1 - pB) / 1000)
print(f"Cramer's V = {V:.4f}")
print(f"odds ratio B versus A = {OR:.4f}, 95% Woolf CI = [{lo:.4f}, {hi:.4f}]")
print(f"difference between proportions = {pB - pA:.3f}, 95% CI = [{pB - pA - 1.96 * se_d:.4f}, {pB - pA + 1.96 * se_d:.4f}]")

# Check of the degrees of freedom with a 3 x 4 table and true independence: 20,000 simulated tables
rng = np.random.default_rng(3)
pr, pc_ = np.array([0.5, 0.3, 0.2]), np.array([0.4, 0.3, 0.2, 0.1])
celle = np.outer(pr, pc_).ravel()
R, m = 20_000, 300
stat = np.empty(R)
for r in range(R):
    t = rng.multinomial(m, celle).reshape(3, 4)
    stat[r] = stats.chi2_contingency(t, correction=False)[0]
soglia = stats.chi2.ppf(0.95, 6)
print(f"3 x 4 table, n = {m}: degrees of freedom 6, mean of chi2 = {stat.mean():.3f}, rejections at 5% = {(stat > soglia).mean():.4f}")

# Output:
# expected = [[135.0, 865.0], [135.0, 865.0]]
# chi2 = 3.8536, degrees of freedom = 1, p = 0.0496
# scipy (without Yates): chi2 = 3.8536, p = 0.0496
# scipy (with Yates):    chi2 = 3.6009, p = 0.0577
# z = 1.9630, z^2 = 3.8536
# Cramer's V = 0.0439
# odds ratio B versus A = 1.2941, 95% Woolf CI = [0.9999, 1.6748]
# difference between proportions = 0.030, 95% CI = [0.0001, 0.0599]
# 3 x 4 table, n = 300: degrees of freedom 6, mean of chi2 = 5.958, rejections at 5% = 0.0488

The expected frequencies are 135 and 865 in both rows, and χ² = 3.8536, p = 0.0496: with Yates’ correction, 3.6009 and p = 0.0577. The z statistic of the comparison between proportions is 1.9630 and its square is identical to χ². The verdict “significant” or “not significant” depends on a correction, and this is a signal: the result is borderline. The effect indices confirm it.

Cramér’s V is 0.0439, that is, almost nothing; the odds ratio is 1.2941 with interval [0.9999; 1.6748], which grazes 1; the difference between proportions is 3 percentage points with interval [0.0001; 0.0599]. The probable effect is small and uncertain: the data are not enough to decide, and a p-value of 0.0496 without these numbers is a poor summary. On the 3 × 4 table, the mean of χ² is 5.958, close to the 6 degrees of freedom (3 − 1)(4 − 1), and it rejects 4.88% of the time.

Little data: Yates and Fisher

When the expected frequencies are small, there are two alternatives. Yates’ correction reduces each deviation |O − E| by 0.5 before squaring it, to compensate for the fact that the counts are integers.

Fisher’s exact test does not use the approximation: it conditions on the row and column totals and computes the exact probability of every table with those totals, which is a hypergeometric distribution; scipy.stats.fisher_exact implements it for 2 × 2 tables (from version 1.15 it also accepts larger tables, with the method parameter). Neither is “the best”: they are two different answers to the same difficulty, and it is worth seeing what they do.

The program takes the table [[8, 2], [1, 5]], with minimum expected frequency 2.62, and computes the three p-values. Then, for two groups of 10 observations with success probability p₁ and p₂, it enumerates all the 121 possible tables, adds up the probabilities of those in which each test rejects at 5%, and obtains the exact probabilities of rejecting: the true level when p₁ = p₂ and the power when they differ. Tables with an empty column (no successes or no failures in total) give no information and do not reject.

# 2 x 2 tables with little data: true (exact) size of chi-square, chi-square with Yates and Fisher's test
import numpy as np
from scipy import stats

print("example table [[8, 2], [1, 5]]:")
t = np.array([[8, 2], [1, 5]])
print(f"  minimum expected = {stats.contingency.expected_freq(t).min():.2f}")
print(f"  chi2 without correction: p = {stats.chi2_contingency(t, correction=False)[1]:.4f}")
print(f"  chi2 with Yates:         p = {stats.chi2_contingency(t, correction=True)[1]:.4f}")
print(f"  Fisher:                  p = {stats.fisher_exact(t)[1]:.4f}")

# Two groups of 10, success probabilities p1 and p2: I enumerate all the 121 possible tables,
# so the probability of rejecting at 5% is EXACT (level if p1 = p2, power if p1 != p2)
n1 = n2 = 10

def prob_rifiuto(p1, p2):
    rif = np.zeros(3)
    for x1 in range(n1 + 1):
        for x2 in range(n2 + 1):
            if x1 + x2 in (0, n1 + n2):                    # an empty column: no evidence, no rejection
                continue
            tab = np.array([[x1, n1 - x1], [x2, n2 - x2]])
            w = stats.binom.pmf(x1, n1, p1) * stats.binom.pmf(x2, n2, p2)
            rif += w * np.array([stats.chi2_contingency(tab, correction=False)[1] < 0.05,
                                 stats.chi2_contingency(tab, correction=True)[1] < 0.05,
                                 stats.fisher_exact(tab)[1] < 0.05])
    return rif

print("situation            chi2    Yates   Fisher   (exact probability of rejecting at 5%)")
for etichetta, p1, p2 in (("H0, p = 0.1", 0.1, 0.1), ("H0, p = 0.3", 0.3, 0.3), ("H0, p = 0.5", 0.5, 0.5),
                          ("H1, 0.2 versus 0.8", 0.2, 0.8), ("H1, 0.1 versus 0.6", 0.1, 0.6)):
    r = prob_rifiuto(p1, p2)
    print(f"{etichetta:19s} {r[0]:7.4f} {r[1]:7.4f} {r[2]:7.4f}")

# Output:
# example table [[8, 2], [1, 5]]:
#   minimum expected = 2.62
#   chi2 without correction: p = 0.0134
#   chi2 with Yates:         p = 0.0510
#   Fisher:                  p = 0.0350
# situation            chi2    Yates   Fisher   (exact probability of rejecting at 5%)
# H0, p = 0.1          0.0090  0.0011  0.0011
# H0, p = 0.3          0.0371  0.0119  0.0119
# H0, p = 0.5          0.0422  0.0128  0.0128
# H1, 0.2 versus 0.8   0.8054  0.6353  0.6353
# H1, 0.1 versus 0.6   0.6590  0.4740  0.4740

On the example table the three tests do not agree at 5%: χ² without correction gives p = 0.0134, Yates 0.0510, Fisher 0.0350. With expected values of 2.62 none of the three is “the true one” by definition, but Fisher is exact conditionally on the totals. The exact comparison refutes two common ideas. First, at this size the χ² does not inflate false positives: it rejects 0.90%, 3.71% and 4.22% with p equal to 0.1, 0.3 and 0.5, so it is prudent.

Second, Fisher and Yates are even more prudent, 0.11%-1.28%, and the prudence is paid for in power: against an alternative 0.2 versus 0.8, the χ² rejects in 80.5% of cases, Yates and Fisher in 63.5%; with 0.1 versus 0.6, 65.9% against 47.4%. The choice is therefore between a test that errs little but makes little use of the data, and a more powerful one that declares a level that is not guaranteed.

Four mistakes that are really made

• Using K − 1 degrees of freedom after estimating a parameter from the data: the true level is 3.0% instead of 5%.
• Declaring a result “significant” with p = 0.0496 without looking at the effect: in the A/B test V = 0.0439 and the interval of the difference starts at 0.0001.
• Applying the χ² with expected values below 2 on rare categories: the true level rises to 7.5%.
• Choosing Yates or Fisher “to be safe” without noticing that power is lost: from 80.5% to 63.5% in one of the exact cases computed above.

The mistakes have in common trust in an unverified rule: the threshold of 5, the automatic correction, the degrees of freedom from memory. Each rule has a reason and a cost, and the cost is measured with a program that enumerates the tables or simulates the counts; if your data sit near a threshold, redo it with your own numbers.

You can now compute by hand and with SciPy the χ² statistic for a goodness of fit and for a contingency table, choose the degrees of freedom (even with estimated parameters), recognise when the approximation is weak, compare χ², Yates and Fisher on cases with little data, and report the effect together with the p-value. In the next article you move on to continuous variables linked to each other, with linear regression.

The same series also includes Nonparametric Statistics: Wilcoxon and Mann-Whitney and Exploratory Data Analysis (EDA) with Python.