Алла Тамбовцева
scipy
для статистики¶Библиотека scipy
(сокращение от scientific Python) включает в себя разные модули, позволяющие выполнять научные расчеты, решать задачи оптимизации, генерировать выборки из (псевдо)случайных величин с заданными параметрами, реализовывать статистические тесты и создавать статистические модели.
Импортируем библиотеку pandas
и модуль stats
из библиотеки scipy
.
import scipy.stats as st
import pandas as pd
Загрузим базу данных (датафрейм) из csv-файла. Описание данных см.здесь.
df = pd.read_csv('https://raw.githubusercontent.com/allatambov/Py-programming-3/master/add/swiss.csv')
df.head()
Unnamed: 0 | Fertility | Agriculture | Examination | Education | Catholic | Infant.Mortality | |
---|---|---|---|---|---|---|---|
0 | Courtelary | 80.2 | 17.0 | 15 | 12 | 9.96 | 22.2 |
1 | Delemont | 83.1 | 45.1 | 6 | 9 | 84.84 | 22.2 |
2 | Franches-Mnt | 92.5 | 39.7 | 5 | 5 | 93.40 | 20.2 |
3 | Moutier | 85.8 | 36.5 | 12 | 7 | 33.77 | 20.3 |
4 | Neuveville | 76.9 | 43.5 | 17 | 15 | 5.16 | 20.6 |
Далее мы попробуем реализовать различные статистические тесты (применить статистические критерии) для сравнения средних значений или распределений в двух выборках.
Предположим, что нам необходимо сравнить средние значения уровня детской смертности в кантонах Швейцарии, где преобладает католическое население и где преобладает протестантское население. Сформируем две выборки на основе имеющихся данных: выберем соответствующие строки в таблице и возьмем столбец Infant.Mortality
.
df[df.Catholic > 50] # для иллюстрации
Unnamed: 0 | Fertility | Agriculture | Examination | Education | Catholic | Infant.Mortality | |
---|---|---|---|---|---|---|---|
1 | Delemont | 83.1 | 45.1 | 6 | 9 | 84.84 | 22.2 |
2 | Franches-Mnt | 92.5 | 39.7 | 5 | 5 | 93.40 | 20.2 |
5 | Porrentruy | 76.1 | 35.3 | 9 | 7 | 90.57 | 26.6 |
6 | Broye | 83.8 | 70.2 | 16 | 7 | 92.85 | 23.6 |
7 | Glane | 92.4 | 67.8 | 14 | 8 | 97.16 | 24.9 |
8 | Gruyere | 82.4 | 53.3 | 12 | 7 | 97.67 | 21.0 |
9 | Sarine | 82.9 | 45.2 | 16 | 13 | 91.38 | 24.4 |
10 | Veveyse | 87.1 | 64.5 | 14 | 6 | 98.61 | 24.5 |
30 | Conthey | 75.5 | 85.9 | 3 | 2 | 99.71 | 15.1 |
31 | Entremont | 69.3 | 84.9 | 7 | 6 | 99.68 | 19.8 |
32 | Herens | 77.3 | 89.7 | 5 | 2 | 100.00 | 18.3 |
33 | Martigwy | 70.5 | 78.2 | 12 | 6 | 98.96 | 19.4 |
34 | Monthey | 79.4 | 64.9 | 7 | 3 | 98.22 | 20.2 |
35 | St Maurice | 65.0 | 75.9 | 9 | 9 | 99.06 | 17.8 |
36 | Sierre | 92.2 | 84.6 | 3 | 3 | 99.46 | 16.3 |
37 | Sion | 79.3 | 63.1 | 13 | 13 | 96.83 | 18.1 |
45 | Rive Droite | 44.7 | 46.6 | 16 | 29 | 50.43 | 18.2 |
46 | Rive Gauche | 42.8 | 27.7 | 22 | 29 | 58.33 | 19.3 |
sample1 = df[df.Catholic > 50]["Infant.Mortality"] # выборка 1
sample1
1 22.2 2 20.2 5 26.6 6 23.6 7 24.9 8 21.0 9 24.4 10 24.5 30 15.1 31 19.8 32 18.3 33 19.4 34 20.2 35 17.8 36 16.3 37 18.1 45 18.2 46 19.3 Name: Infant.Mortality, dtype: float64
sample2 = df[df.Catholic <= 50]["Infant.Mortality"] # выборка 2
sample2
0 22.2 3 20.3 4 20.6 11 16.5 12 19.1 13 22.7 14 18.7 15 21.2 16 20.0 17 20.2 18 10.8 19 20.0 20 18.0 21 22.4 22 16.7 23 15.3 24 21.0 25 23.8 26 18.0 27 16.3 28 20.9 29 22.5 38 20.3 39 20.5 40 18.9 41 23.0 42 20.0 43 19.5 44 18.0 Name: Infant.Mortality, dtype: float64
Теперь приступим к формальной проверке гипотез.
T-test используется для сравнения средних значений двух генеральных совокупностей в предположении, что обе выборки взяты из нормального распределения. Существует две разновидности t-теста: t-тест для независимых выборок и t-тест для связных (парных) выборок. В связных выборках объекты связаны друг с другом. Пример связных выборок: значения уровня смертности в одних и тех же кантонах до и после какой-нибудь реформы)
Нулевая гипотеза: средние значения двух генеральных совокупностей, откуда взяты выборки, равны, то есть $H_0: a_1=a_2$.
Альтернативная гипотеза: средние значения двух генеральных совокупностей не равны, то есть: $H_1: a_1 \ne a_2$
В stats
в t-тесте в качестве альтернативной гипотезы используется двусторонняя альтернатива (средние не равны) и всегда выводится соответствующее p-value (two-tailed). То же будет характерно для всех последующих тестов. Так как наши выборки независимы, нам нужна функция ttest_ind()
, от independent.
st.ttest_ind(sample1, sample2)
Ttest_indResult(statistic=1.1297965130690064, pvalue=0.26454837328688746)
Что возвращает эта функция? Наблюдаемое значение t-статистики и p-value. Результат понятный, только более лакончиный по сравнению с тем, что выводит R и другие статистические пакеты.
Выводы: так как p-value больше любого конвенционального уровня значимости (1%, 5%, 10%), на имеющихся данных на любом разумном уровне значимости нет оснований отвергнуть нулевую гипотезу. Средний уровень детской смертности в католических и протестантских районах можно считать одинаковым.
По умолчанию считается, что дисперсии генеральных овокупностей равны. Часто это бывает не так, и такое предположение без формальной проверки и без содержательных соображений может казаться нереалистичным. Если мы предполагаем, что дисперсии генеральных совокупностей не равны, то это можно учесть, добавив аргумент equal_var
:
st.ttest_ind(sample1, sample2, equal_var = False)
Ttest_indResult(statistic=1.0863471703503398, pvalue=0.28551301767919196)
Принципиальных отличий в данном случае в результатах не наблюдается. А формально проверить гипотезу о равенстве дисперсий двух генеральных совокупностей (которые описываются двумя случайными величинами) можно с помощью F-критерия.
Нулевая гипотеза: дисперсии двух генеральных совокупностей равны, то есть $H_0: \sigma_1^2 = \sigma_2^2$
Альтернативная гипотеза: дисперсии двух генеральных совокупностей не равны, то есть $H_1: \sigma_1^2 \ne \sigma_2^2$.
Реализовать F-тест, который нам нужен именно для этого случая, в Python сразу не получится: встроенная функция f_oneway()
используется для однофакторного дисперсионного анализа (ANOVA), речь о котором пойдёт далее. Можно попробовать реализовать этот тест «вручную», рассчитав частное выборочных дисперсий и поработав с F-распределением, но давайте пойдём другим путём и воспользуемся другими критериями для сравнения дисперсий, которые явно встроены в Python.
Например, тест Левена:
st.levene(sample1, sample2)
LeveneResult(statistic=1.0380526498527436, pvalue=0.3137214012581385)
На любом разумном уровне значимости нет оснований отвергнуть нулевую гипотезу о равенстве дисперсий.
Если мы не можем считать распределение генеральных совокупностей, откуда взяты выборки, нормальным, то следует использовать методы, основанные не на самих наблюдениях в выборках, а на их рангах. Для сравнения распределений (иногда речь идет о сравнении медиан) используются тесты Уилкоксона и Манна-Уитни. Начнем с теста Уилкоксона (не проверяем, является ли распределение данных нормальным, просто для примера используем те же выборки).
st.wilcoxon(sample1, sample2)
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-15-6bfff17a9e4f> in <module>() ----> 1 st.wilcoxon(sample1, sample2) /usr/local/lib/python3.5/dist-packages/scipy/stats/morestats.py in wilcoxon(x, y, zero_method, correction) 2374 x, y = map(asarray, (x, y)) 2375 if len(x) != len(y): -> 2376 raise ValueError('Unequal N in wilcoxon. Aborting.') 2377 d = x - y 2378 ValueError: Unequal N in wilcoxon. Aborting.
Неудача! Проблема в том, что реализация критерия Уилкоксона в stats
требует, чтобы выборки были одинакового размера. Но это ограничение можно обойти, просто выбрав другой критерий – критерий Манна-Уитни, который используется для аналогичных задач.
Нулевая гипотеза: выборки взяты из одного и того же распределения, то есть $H_0: F(x) = G(x)$
Альтернативная гипотеза: выборки взяты из разных распределений, то есть $H_1: F(x) \ne G(x)$.
st.mannwhitneyu(sample1, sample2)
MannwhitneyuResult(statistic=235.5, pvalue=0.2920645577220585)
Опять же, на имеющихся данных на любом уровне значимости нет оснований отвергнуть нулевую гипотезу. Выборки взяты из одного и того же распределения.
Если выборок больше, чем две, то использовать указанные выше критерии нельзя. В предположении, что все выборки взяты из нормального распределения, для сравнения средних значений более чем в двух группах используется однофакторный дисперсионный анализ (ANOVA, analysis of variance).
Нулевая гипотеза: средние значения по всем $k$ группам (во всех генеральных совокупностях) равны, то есть $H_0: a_1 = a_2 = \dots = a_k$
Альтернативная гипотеза: средние значения по всем группам (во всех генеральных совокупностях) не равны.
Чтобы не создавать искусственные группы на основе данных в swiss, загрузим таблицу с весами цыплят, которых кормили разным кормом :) Описание см. здесь.
dat = pd.read_csv('https://raw.githubusercontent.com/allatambov/Py-programming-3/master/add/chickwts.csv')
dat.head()
Unnamed: 0 | weight | feed | |
---|---|---|---|
0 | 1 | 179 | horsebean |
1 | 2 | 160 | horsebean |
2 | 3 | 136 | horsebean |
3 | 4 | 227 | horsebean |
4 | 5 | 217 | horsebean |
Задание: разбить датафрейм на группы с помощью groupby
по переменной feed
и сохранить значения weight
в словарь.
Решение:
wgt = {}
for name, d in dat.groupby('feed'):
wgt[name] = d.weight
wgt
{'casein': 59 368 60 390 61 379 62 260 63 404 64 318 65 352 66 359 67 216 68 222 69 283 70 332 Name: weight, dtype: int64, 'horsebean': 0 179 1 160 2 136 3 227 4 217 5 168 6 108 7 124 8 143 9 140 Name: weight, dtype: int64, 'linseed': 10 309 11 229 12 181 13 141 14 260 15 203 16 148 17 169 18 213 19 257 20 244 21 271 Name: weight, dtype: int64, 'meatmeal': 48 325 49 257 50 303 51 315 52 380 53 153 54 263 55 242 56 206 57 344 58 258 Name: weight, dtype: int64, 'soybean': 22 243 23 230 24 248 25 327 26 329 27 250 28 193 29 271 30 316 31 267 32 199 33 171 34 158 35 248 Name: weight, dtype: int64, 'sunflower': 36 423 37 340 38 392 39 339 40 341 41 226 42 320 43 295 44 334 45 322 46 297 47 318 Name: weight, dtype: int64}
Теперь ANOVA (f_oneway
от One-Way ANOVA):
st.f_oneway(wgt['casein'], wgt['horsebean'], wgt['linseed'], wgt['meatmeal'], wgt['soybean'], wgt['sunflower'])
F_onewayResult(statistic=15.364799774712534, pvalue=5.936419853471331e-10)
Функция возвращает наблюдаемое значение F-статистики и p-value. В данном случае p-value близко к 0, поэтому гипотезу о равенстве средних генеральных совокупностей по группам можно отвергнуть на 1% уровне значимости. Средний вес цыплят, которых кормили разным кормом, отличается (ещё бы, horsebean или sunflower!).
Критерий Краскела-Уоллиса используется, когда нам необходимо сравнить распределения более, чем в двух группах в предположении, что выборки взяты не из нормального распределения (распределения неизвестны).
Нулевая гипотеза: выборки взяты из одного и того же распределения, то есть: $H_0: F(x) = G(x) = \dots = H(x)$
Альтернативная гипотеза: выборки взяты из разных распределений.
st.kruskal(wgt['casein'], wgt['horsebean'], wgt['linseed'], wgt['meatmeal'], wgt['soybean'], wgt['sunflower'])
KruskalResult(statistic=37.34271769425624, pvalue=5.112829511937094e-07)