Глава 9. Проверка гипотез

В датасетах, которые мы до сих пор исследовали в этой книге, мы наблюдали различия между группами людей (и пингвинов), корреляции между переменными и углы наклона линий регрессии. Подобные результаты называются наблюдаемыми эффектами, поскольку они проявляются в выборке — в отличие от реальных эффектов в генеральной совокупности, которые обычно не поддаются непосредственному наблюдению. Когда мы видим такой эффект, стоит задуматься, может ли он присутствовать в более многочисленной генеральной совокупности или его появление в выборке случайно.

Существуют разные способы ответить на этот вопрос, такие как проверка нулевой гипотезы Фишера, теория принятия решений Неймана — Пирсона или байесовская проверка гипотез. Здесь я представлю комбинацию этих подходов, часто используемую на практике.

Эксперимент с монетой

Начнем с простого примера. Когда в 2002 году были введены в обращение монеты евро, один пытливый нумизмат-любитель провел эксперимент с бельгийской монетой номиналом в один евро. Он раскрутил монету на ребре 250 раз, из которых она упала орлом вверх 140 раз, а решкой — 110 раз. Если монета идеально сбалансирована, то стоит ожидать только 125 орлов, так что эти данные говорят о несимметричности монеты. С другой стороны, не стоит ожидать, что в каждом таком эксперименте выпадет ровно 125 орлов, поэтому вполне возможно, что монета на самом деле справедливая, а видимое отклонение от ожидаемого значения обусловлено случайностью. Чтобы убедиться в правдоподобности такого объяснения, проведем проверку гипотезы (hypothesis test).

Для симметричной монеты будем использовать следующую функцию, вычисляющую абсолютную разность между наблюдаемым числом выпавших орлов (heads) и ожидаемым (expected):

n = 250

p = 0.5

 

def abs_deviation(heads):

    expected = n * p

    return np.abs(heads - expected)

В наблюдаемых данных это отклонение составляет 15:

heads = 140

tails = 110

 

observed_stat = abs_deviation(heads)

observed_stat

15.0

Если монета действительно справедливая, то эксперимент с ее вращением можно смоделировать путем генерации последовательности случайных строк — равновероятно либо H, либо T — и подсчитав количество выпадений H:

def simulate_flips():

    flips = np.random.choice(["H", "T"], size=n)

    heads = np.sum(flips == "H")

    return heads

Каждый раз при вызове этой функции мы получаем исход модельного эксперимента:

simulate_flips()

119

Цикл ниже повторяет этот эксперимент множество раз, вычисляет отклонение для каждого из них и с помощью спискового включения сохраняет результаты в виде списка:

simulated_stats = [abs_deviation(simulate_flips()) for i in range(10001)]

В предположении, что монета является справедливой, мы получаем выборку из распределения отклонений. Вот как выглядит это распределение:

from empiricaldist import Pmf

 

pmf_effects = Pmf.from_seq(simulated_stats)

pmf_effects.bar()

 

decorate(xlabel="Абсолютное отклонение", ylabel="Вероятность")

Значения вблизи 0 — самые распространенные; значения больше 10 встречаются реже. Вспоминая, что отклонение в наблюдаемых данных равно 15, можно прийти к выводу, что отклонения такого порядка редки, но не исключены. В этом примере результаты моделирования больше или равны 15 примерно в 7.1 % случаев:

(np.array(simulated_stats) >= 15).mean() * 100

7.079292070792921

Итак, если монета справедливая, то можно ожидать отклонения такой же величины, как наблюдаемое, примерно в 7.1 % случаев, просто в силу случайности.

Это приводит к выводу, что эффект такого масштаба встречается нечасто, но он, безусловно, возможен даже с симметричной монетой. На базе проведенного эксперимента мы не можем исключить возможность того, что монета справедливая.

На этом примере видна логика проверки статистических гипотез:

• Мы начали с наблюдения — 140 выпадений орла из 250 вращений монеты — и гипотезы, что монета несимметричная, другими словами, что вероятность выпадения орла отличается от 50 %.

• Мы выбрали статистику критерия, или тестовую статистику (test statistic), которая количественно измеряет величину наблюдаемого эффекта. В данном примере статистикой критерия является абсолютное отклонение от ожидаемого исхода.

• Далее мы определили нулевую гипотезу, которая представляет собой модель, построенную на предположении, что наблюдаемый эффект обусловлен случайностью. В данном случае нулевая гипотеза заключается в том, что монета симметричная.

• Затем мы вычислили p-значение, которое представляет собой вероятность обнаружения наблюдаемого эффекта в условиях верности нулевой гипо­тезы. В этом примере p-значение — это вероятность отклонения, большего или равного 15.

На последнем шаге интерпретируется результат. Если p-значение мало, то делается вывод, что эффект вряд ли может быть обусловлен случайностью. Если оно велико, мы приходим к заключению, что эффект может быть правдоподобно объяснен случайностью. А если он находится где-то посередине, как в этом примере, то можно сделать вывод, что эффект вряд ли появился случайно, но нельзя исключать такую возможность.

Любая проверка гипотезы строится из этих элементов — статистики критерия, нулевой гипотезы и p-значения.

Проверка разницы средних

На материале датасета NSFG мы видели, что средняя продолжительность беременности первым ребенком немного дольше, чем последующими детьми. Теперь посмотрим, может ли эта разница быть случайной.

Функция get_nsfg_groups считывает данные, выбирает случаи рождения живых младенцев (live) и разделяет их на группы первенцев (firsts) и остальных новорожденных (others):

from nsfg import get_nsfg_groups

 

live, firsts, others = get_nsfg_groups()

Теперь можно выбрать продолжительность беременности в неделях для обеих групп:

data = firsts["prglngth"].values, others["prglngth"].values

Следующая функция принимает на вход данные в виде кортежа из двух последовательностей и вычисляет абсолютную разность средних значений:

def abs_diff_means(data):

    group1, group2 = data

    diff = np.mean(group1) - np.mean(group2)

    return np.abs(diff)

Наблюдаемая разница в продолжительности беременности первыми и остальными детьми составляет 0.078 недели:

observed_diff = abs_diff_means(data)

observed_diff

0.07803726677754952

Итак, гипотеза, которую мы проверим, заключается в том, что продолжительность беременности первым ребенком обычно дольше. Нулевая гипотеза заключается в том, что продолжительность беременности в обеих группах на самом деле одинакова, а кажущаяся разница обусловлена случайностью. Если продолжительность беременности в обеих группах одинакова, то эти группы можно объединить. Чтобы смоделировать такой эксперимент, можно при помощи функции библиотеки NumPy shuffle расположить объединенные значения в произвольном порядке, а затем разделить их на две группы первоначального размера с использованием индексов срезов:

def simulate_groups(data):

    group1, group2 = data

    n, m = len(group1), len(group2)

 

    pool = np.hstack(data)

    np.random.shuffle(pool)

    return pool[:n], pool[-m:]

Эта функция возвращает кортеж из двух последовательностей, который можно передать в функцию abs_diff_means:

abs_diff_means(simulate_groups(data))

0.031193045602279312

В следующем цикле многократно повторяется этот эксперимент и для каждого моделируемого датасета вычисляется абсолютная разность средних зна­чений:

simulated_diffs = [abs_diff_means(simulate_groups(data)) for i in range(1001)]

Чтобы отобразить результаты, воспользуемся функцией ниже, которая принимает выборку результатов моделирования и создает объект Pmf, аппроксимирующий ее распределение:

from scipy.stats import gaussian_kde

from empiricaldist import Pmf

 

def make_pmf(sample, low, high):

    kde = gaussian_kde(sample)

    qs = np.linspace(low, high, 201)

    ps = kde(qs)

    return Pmf(ps, qs)

Вот как выглядит распределение результатов моделирования. Заштрихованная область показывает случаи, когда при справедливости нулевой гипотезы разница в средних значениях превышает наблюдаемую разницу. Площадь этой области дает p-значение:

from thinkstats import fill_tail

 

pmf = make_pmf(simulated_diffs, 0, 0.2)

pmf.plot()

fill_tail(pmf, observed_diff, "right")

decorate(xlabel="Абсолютная разность средних (недели)", ylabel="Плотность

вероятности")

Следующая функция вычисляет p-значение, представляющее собой долю смоделированных значений, которые больше наблюдаемого или равны ему:

def compute_p_value(simulated, observed):

    """Доля смоделированных значений, которые больше наблюдаемого или равны ему."""

    return (np.asarray(simulated) >= observed).mean()

В нашем примере p-значение примерно равно 18 %, из чего следует, что разница в 0.078 недели может быть правдоподобно объяснена случайностью:

compute_p_value(simulated_diffs, observed_diff)

0.1838161838161838

Получив такой результат, мы не можем быть уверены, что продолжительность беременности первым ребенком обычно дольше. Возможно, что наблюдаемая разница в этом датасете обусловлена случайностью.

Обратите внимание, что в обоих примерах проверки гипотез присутствовали одни и те же элементы: статистика критерия, нулевая гипотеза и модель нулевой гипотезы. В последнем примере статистика критерия представляет собой абсолютную разность средних значений. Нулевая гипотеза заключается в том, что распределение продолжительности беременности в обеих группах на самом деле одинаковое. И мы смоделировали нулевую гипотезу, объединив данные из обеих групп, перетасовав и разделив их на две группы с теми же размерами, что и у исходных данных. Такой процесс перетасовки по-другому называется перестановкой (permutation).

Такой вычислительный подход к проверке гипотез позволяет легко комбинировать перечисленные выше элементы для проверки различных статистических гипотез.

Прочие статистики критерия

Возникает вопрос: может ли быть срок беременности первым ребенком не просто более продолжительным, но и более изменчивым? Чтобы проверить эту гипотезу, в качестве статистики критерия можно использовать абсолютную разность стандартных отклонений двух групп. Следующая функция вычисляет эту статистику:

def abs_diff_stds(data):

    group1, group2 = data

    diff = np.std(group1) - np.std(group2)

    return np.abs(diff)

В данных NSFG разница в стандартных отклонениях составляет около 0.18:

observed_diff = abs_diff_stds(data)

observed_diff

0.17600895913991677

Чтобы проверить, может ли это различие объясняться случайностью, снова используем перестановку. Следующий цикл многократно моделирует нулевую гипотезу и вычисляет абсолютную разность стандартных отклонений для каждого смоделированного датасета:

simulated_diffs = [abs_diff_stds(simulate_groups(data)) for i in range(1001)]

Вот как выглядит распределение результатов. Штриховка по-прежнему обозначает область, где статистика критерия, при условии истинности нулевой гипотезы, превышает наблюдаемую разницу:

pmf = make_pmf(simulated_diffs, 0, 0.5)

pmf.plot()

fill_tail(pmf, observed_diff, "right")

decorate(xlabel="Абсолютная разность стандартных отклонений (недели)",

ylabel="Плотность вероятности")

Площадь этой области можно оценить, вычислив долю результатов, которые равны наблюдаемой разнице или превышают ее:

compute_p_value(simulated_diffs, observed_diff)

0.17082917082917082

В данном случае p-значение составляет около 0.17, так что даже в случае, если эти две группы одинаковы, появление разницы такой величины вполне правдоподобно. Подводя итог, мы не можем быть уверены в том, что продолжительность беременности первым ребенком в целом более изменчива: разница, присутствующая в этом датасете, может быть случайной.

Проверка значимости корреляции

Такую же схему можно использовать для проверки значимости корреляции. Например, в датасете NSFG есть корреляция между весом новорожденного и возрастом матери: у матерей более старшего возраста в среднем рождаются более крупные дети. Но может ли это явление объясняться случайностью?

Давайте выяснять. Начнем с подготовки данных. Из всех случаев рождения живых детей выберем только те, когда известны возраст матери и вес новорожденного:

valid = live.dropna(subset=["agepreg", "totalwgt_lb"]) valid.shape

(9038, 244)

Теперь выберем нужные столбцы:

ages = valid["agepreg"]

birthweights = valid["totalwgt_lb"]

Следующая функция принимает на вход кортеж из двух последовательностей, xs и ys, и вычисляет модуль корреляции, которая может быть положительной или отрицательной:

def abs_correlation(data):

    xs, ys = data

    corr = np.corrcoef(xs, ys)[0, 1]

    return np.abs(corr)

В датасете NSFG корреляция составляет около 0.07:

data = ages, birthweights

observed_corr = abs_correlation(data)

observed_corr

0.0688339703541091

В качестве нулевой гипотезы примем отсутствие корреляции между возрастом матери и весом новорожденного. Перетасовав наблюдаемые значения, смоделируем ситуацию, в которой распределения возраста матери и веса новорожденного остаются теми же, но переменные перестают быть связаны.

Следующая функция принимает кортеж из последовательностей xs и ys, перетасовывает xs и возвращает кортеж, содержащий перетасованную xs и исходную ys. Мы бы добились того же результата, если бы перетасовали только ys или обе последовательности одновременно:

def permute(data):

    xs, ys = data

    new_xs = xs.values.copy()

    np.random.shuffle(new_xs)

    return new_xs, ys

Корреляция перетасованных значений обычно близка к 0:

abs_correlation(permute(data))

0.0019269515502894237

В цикле ниже генерируется множество перетасованных датасетов и в каждом из них вычисляется корреляция:

simulated_corrs = [abs_correlation(permute(data)) for i in range(1001)]

Вот как выглядит распределение полученных результатов. Вертикальная пунк­тирная линия обозначает наблюдаемую корреляцию:

pmf = make_pmf(simulated_corrs, 0, 0.07)

pmf.plot()

plt.axvline(observed_corr, color="gray", ls=":")

decorate(xlabel="Абсолютное значение корреляции", ylabel="Плотность

вероятности")

Обращает на себя внимание то, что наблюдаемая корреляция находится в хвосте распределения, где видимая область под кривой отсутствует. Если мы попытаемся вычислить p-значение, результат будет равен 0, это указывает на то, что ни в одной из моделей корреляция в перетасованных данных не превышает наблюдаемого значения:

compute_p_value(simulated_corrs, observed_corr)

0.0

Эти цифры говорят, что значение p, вероятно, меньше 1 на 1000, хотя оно на самом деле и не равно нулю. Маловероятно, что корреляция перетасованных данных превысит наблюдаемое значение, но такую возможность нельзя исключать.

Когда p-значение мало, по сложившейся традиции меньше 0.05, можно сказать, что результат статистически значим. Но такая интерпретация p-значения всегда имела свои недостатки, и постепенно от нее отказываются.

Один из таких недостатков — произвольность традиционного значения порога, которое к тому же подходит не для всех применений. Другая сложность заключается в том, что использование термина «значимый» может приводить к неверным выводам, поскольку оно намекает, что эффект важен на практике. Корреляция между возрастом матери и весом новорожденного — хороший пример. Она статистически значима, но настолько мала, что практически не важна.

С учетом вышеказанного можно прибегнуть к качественной интерпретации p-значения:

• если p-значение велико, то вполне вероятно, что наблюдаемый эффект мог возникнуть случайно;

• если p-значение мало, то, как правило, вероятность того, что эффект вызван случайностью, можно исключить. Но нужно помнить, что он все еще может объясняться нерепрезентативностью выборки или ошибками измерений.

Проверка пропорций

В качестве завершающего примера данной главы возьмем случай, когда выбор статистики критерия не так очевиден. Представьте, что вы управляете казино и подозреваете, что какой-то клиент жульничает и пользуется игральной костью, модифицированной таким образом, чтобы одна из граней выпадала с большей вероятностью, чем другие. Вы задерживаете предполагаемого мошенника и конфискуете игральную кость, но теперь вам нужно доказать, что она поддельная. Вы бросаете кубик 60 раз и записываете число повторений каждого исхода от 1 до 6. Вот результаты в виде объекта класса Hist:

from empiricaldist import Hist

 

qs = np.arange(1, 7)

freqs = [8, 9, 19, 5, 8, 11]

observed = Hist(freqs, qs)

observed.index.name = "исход"

observed

исход	частота
1	8
2	9
3	19
4	5
5	8
6	11

Вы ожидаете, что в среднем каждое значение должно выпадать 10 раз. В этих данных значение 3 появляется чаще, чем ожидалось, а значение 4 — реже. Но могут ли такие расхождения возникать по стечению обстоятельств?

Чтобы проверить эту гипотезу, начнем с вычисления ожидаемой частоты (expected) всех исходов (outcomes):

num_rolls = observed.sum()

outcomes = observed.qs

expected = Hist(num_rolls / 6, outcomes)

Следующая функция принимает на вход наблюдаемые частоты и вычисляет сумму их абсолютных разностей с ожидаемыми частотами:

def total_abs_deviation(observed):

    return np.sum(np.abs(observed - expected))

У наблюдаемых данных эта статистика критерия равна 20:

observed_dev = total_abs_deviation(observed)

observed_dev

20.0

Функция ниже принимает на вход наблюдаемые данные, моделирует бросание кости то же самое количество раз и возвращает объект класса Hist, содержащий смоделированные частоты:

def simulate_dice(observed):

    num_rolls = np.sum(observed)

    rolls = np.random.choice(observed.qs, num_rolls, replace=True)

    hist = Hist.from_seq(rolls)

    return hist

В следующем цикле многократно повторяется данный эксперимент и вычисляется общее абсолютное отклонение для каждого из них:

simulated_devs = [total_abs_deviation(simulate_dice(observed))

                  for i in range(1001)]

Вот как выглядит распределение статистики критерия в условиях нулевой гипотезы. Обратите внимание, что общая сумма всегда четная, потому что каждый раз, когда один из исходов выпадает чаще, чем ожидалось, какой-то другой исход автоматически выпадает реже:

pmf_devs = Pmf.from_seq(simulated_devs)

pmf_devs.bar()

 

decorate(xlabel="Общее абсолютное отклонение", ylabel="Вероятность")

Теперь становится понятно, что общее отклонение в 20 очков не является необычным. Его p-значение составляет около 13 %, поэтому нельзя с уверенностью заявить, что игральная кость поддельная:

compute_p_value(simulated_devs, observed_dev)

0.13086913086913088

Но выбранная нами статистика критерия — не единственно возможный вариант. Для решения подобной задачи чаще принято использовать статистику хи-квадрат, которую можно вычислить так:

def chi_squared_stat(observed):

    diffs = (observed - expected) ** 2

    return np.sum(diffs / expected)

Возведение отклонений в квадрат (вместо использования абсолютных значений) придает больший вес значительным отклонениям. Деление на значения из expected стандартизирует отклонения, хотя в данном случае это не влияет на результаты, поскольку все ожидаемые частоты равны:

observed_chi2 = chi_squared_stat(observed)

observed_chi2

11.6

Статистика хи-квадрат наблюдаемых данных равна 11.6. Эта цифра сама по себе не имеет большого значения, но ее можно сравнить с результатами моделирования. В следующем цикле генерируется множество модельных датасетов и вычисляется статистика хи-квадрат для каждого из них:

simulated_chi2 = [chi_squared_stat(simulate_dice(observed)) for i

in range(1001)]

Вот как выглядит распределение этой статистики критерия в условиях нулевой гипотезы.

Штриховкой обозначена область, где результаты превышают наблюдаемое значение:

pmf = make_pmf(simulated_chi2, 0, 20)

pmf.plot()

fill_tail(pmf, observed_chi2, "right")

decorate(xlabel="Статистика хи-квадрат", ylabel="Плотность вероятности")

Площадь заштрихованной области опять же дает p-значение:

compute_p_value(simulated_chi2, observed_chi2)

0.04495504495504495

Полученное с помощью статистики хи-квадрат p-значение составляет около 0.04, что значительно меньше, чем получилось по формуле общего отклонения — 0.13. Если серьезно относиться к пороговому значению в 5 %, то можно считать этот эффект статистически значимым. Но, принимая во внимание результаты обеих проверок, я бы сказал, что результаты неубедительны. Я бы не исключал возможности того, что игральная кость поддельная, но и не стал бы осуждать подозреваемого в мошенничестве.

Этот пример демонстрирует один важный момент: p-значение зависит от выбора как статистики критерия, так и модели нулевой гипотезы. И принятые решения иногда определяют, является эффект статистически значимым или нет.

Глоссарий

Проверка гипотез (hypothesis testing)

Набор методов для проверки правдоподобности того, что наблюдаемый эффект может быть следствием случайной выборки.

Статистика критерия (test statistic)

Статистический показатель, используемый при проверке гипотезы для количественной оценки величины наблюдаемого эффекта.

Нулевая гипотеза (null hypothesis)

Статистическая модель, основанная на предположении, что эффект, наблюдаемый в выборке, отсутствует в генеральной совокупности.

Перестановка (permutation)

Способ моделирования нулевой гипотезы путем случайного перемешивания набора данных.

p-значение (p-value)

Вероятность появления эффекта величиной с тот, который наблюдается в условиях нулевой гипотезы.

Статистически значимый (statistically significant)

Эффект статистически значим, если p-значение меньше выбранного порога, часто равного 5 %. При большом объеме данных наблюдаемый эффект может быть статистически значимым, даже если он слишком мал, чтобы иметь практическое значение.

Упражнения

Упражнение 9.1

Испытаем методику проверки гипотез на данных о пингвинах из раздела «Выборочное распределение» главы 8 на с. 151. Инструкции по загрузке данных приведены в Jupyter-блокноте этой главы.

Для начала считаем данные и выберем пингвинов вида антарктический пингвин:

penguins = pd.read_csv("penguins_raw.csv").dropna(subset=["Body Mass (g)"])

chinstrap = penguins.query('Species.str.startswith("Chinstrap")')

chinstrap.shape

(68, 17)

Теперь извлечем вес самцов (male) и самок (female) пингвинов в килограммах:

male = chinstrap.query("Sex == 'MALE'")

weights_male = male["Body Mass (g)"] / 1000

weights_male.mean()

3.9389705882352937

female = chinstrap.query("Sex == 'FEMALE'")

weights_female = female["Body Mass (g)"] / 1000

weights_female.mean()

3.5272058823529413

Используйте функции abs_diff_means и simulate_groups для моделирования большого количества датасетов в условиях нулевой гипотезы равенства распределений весов двух групп и вычислите разницу в средних значениях для каждого из датасетов. Сравните результаты моделирования с наблюдаемой разницей и вычислите p-значение. Насколько правдоподобно то, что видимая разница между группами обусловлена случайностью?

Упражнение 9.2

Извлечем из данных о пингвинах из предыдущего упражнения глубину (depth) и длину (length) верхнего края клюва (culmen) самок пингвинов:

data = female["Culmen Depth (mm)"], female["Culmen Length (mm)"]

Корреляция между этими переменными примерно равна 0.26:

observed_corr = abs_correlation(data)

observed_corr

0.2563170802728449

Проверим, может ли эта корреляция возникнуть случайно, когда ее в действительности нет. Используйте функцию permute, чтобы создать множество перестановок этих данных, и функцию abs_correlation, чтобы вычислить корреляцию для каждой такой перестановки. Постройте график распределения корреляций в условиях нулевой гипотезы и вычислите p-значение для наблюдаемой корреляции. Как можно интерпретировать этот результат?

---