Введение в медицинскую статистику
Автор: old.medach.pro
Публикация: 18.02.2018
Может быть, вы видели, как в споре на околомедицинские темы собеседники обмениваются ссылками на Пабмед, не всегда ознакомившись хотя бы с абстрактами исследований. Забавно, что обычно такой подход себя оправдывает, ведь мало кто станет утруждать себя поиском и чтением полного текста публикаций; одно только наличие ссылки на некую работу может поколебать уверенность оппонента в своей правоте. Но, как вы понимаете, сама по себе ссылка на исследование не является полноценным аргументом; всегда стоит обратить внимание на то, как были собраны и обработаны данные. Для этого не нужно профессионально разбираться в статистике, достаточно знать основы. Кроме того, эти основы пригодятся, если вы соберетесь публиковать собственное исследование.
Сразу скажу, что весьма полезно (хотя и не обязательно) уметь использовать специальные инструменты для работы с данными. Это могут быть как языки программирования (R, Python), так и специальные программы вроде Excel или SPSS. В этой статье расчеты и построение графиков выполнены на Питоне. Если нажмете на эту кнопку, то увидите версию статьи с кодом:
#Импортируем библиотеки.
#Для работы с датасетами:
import pandas as pd
#Математика:
import numpy as np
#Чтобы строить графики:
import matplotlib.pyplot as plt
import matplotlib.patches as patches
get_ipython().magic('matplotlib inline')
#Чтобы построенные графики стали красивыми:
import seaborn as sns
#Собственно статистика:
from scipy import stats
#Отключим предупреждения:
import warnings
warnings.filterwarnings('ignore')
Выборка и генеральная совокупность
Должен ли повар, который варит кашу, съесть всю кастрюлю, чтобы оценить плод своих кулинарных усилий? Вероятно, нет: он может попробовать чуть-чуть, исходя из того простого соображения, что остальная каша в кастрюле похожа на ту, что оказалась в его ложке. В повседневной жизни мы постоянно экстраполируем представление о части на целое. Возможно, вы любите мандарины, но чтобы это понять, вам не пришлось съесть все мандарины на свете.
Множество объектов, о котором исследователь собирается делать выводы, называется генеральной совокупностью. Например, это могут быть все люди, страдающие диабетом. Или все женщины старше 65, перенесшие инсульт. Или все мандарины на свете. Собрать информацию обо всех объектах обычно невозможно, так что используются выборки. Выборка — это та часть генеральной совокупности, с которой работает исследователь. Если выборка отражает изучаемые свойства генеральной совокупности, она называется репрезентативной, а если не отражает, то нерепрезентативной (вряд ли пакет гнилых мандаринов убедит вас, что все мандарины на свете плохи на вкус; впрочем, не всегда нерепрезентативность выборки так очевидна).
Эталонный пример работы с нерепрезентативной выборкой — опрос, проведенный в 1936 году журналом "The Literary Digest": были разосланы бюллетени десяти миллионам человек, среди которых были подписчики журнала и люди, случайным образом выбранные из телефонных справочников и реестров владельцев автомобилей. Два с половиной миллиона человек вернули их заполненными, и большинство из них — 57% — в этом опросе поддерживали республиканца Альфа Лэндона. Демократ Франклин Рузвельт набрал 40%, но когда дело дошло до выборов, — победил. Позже стало понятно, почему так произошло: во-первых, большая часть подписчиков журнала были республиканцами, а во-вторых, телефонами и автомобилями владели главным образом состоятельные люди, которые тоже в основном были республиканцами.
Описательная статистика
Гистограмма
Первым делом нужно составить представление о выборке и описать ее. Наиболее удобное графическое представление выборки — гистограмма. Чтобы понять, как устроена гистограмма, возьмем для примера таблицу с данными из знаменитого исследования 1885 года английского ученого Френсиса Гальтона, в котором он изучал взаимосвязь между ростом взрослых людей и ростом их родителей. В таблицу вошли данные о росте детей (то есть уже взрослых людей) и их родителей из 197 семей. Давайте разберемся с ростом детей.
#Загружаем данные:
galton = pd.read_csv('Galton.csv')
#Выводим первые пять строк:
galton.head()
| Family | Father | Mother | Gender | Height | Kids | |
| 0 | 1 | 78.5 | 67.0 | M | 73.2 | 4 |
| 1 | 1 | 78.5 | 67.0 | F | 69.2 | 4 |
| 2 | 1 | 78.5 | 67.0 | F | 69.0 | 4 |
| 3 | 1 | 78.5 | 67.0 | F | 69.0 | 4 |
| 4 | 2 | 75.5 | 66.5 | M | 73.5 | 4 |
Как формально описать это распределение?
Для начала можно указать на область числовой оси, где единицы наблюдения имеют тенденцию к концентрации. А потом показать, насколько широко и симметрично они рассеиваются в стороны. Для этого используются специальные показатели, которые называются мерами центральной тенденции и мерами рассеяния.
Меры центральной тенденции
Это значения, которые показывают ту область числовой оси, где данные имеют тенденцию к концентрации.
Вообще говоря, мер центральной тенденции довольно много, но наиболее часто используются среднее арифметическое, мода и медиана, о них и поговорим.
Как найти среднее арифметическое для нашей выборки, все знают: \frac{1}{n} \sum_{i=1}^{n} x_{i} — сложить рост всех людей и поделить на количество людей. Это прекрасная характеристика выборки, но у нее есть один досадный недостаток: она неустойчива к выбросам. Если вы забыли, что такое выброс, представьте себе, что когда вы собирали выборку и заполняли таблицу, к одному значению приписали лишний нолик. Или, допустим, вы решили оценить средний месячный доход жителя Сомали, но в число респондентов затесался Билл Гейтс, приехавший туда на выходные. Если выборка небольшая, присутствие необычно выдающихся объектов сильно сместит среднее арифметическое, так что в некоторых случаях от выбросов стоит избавляться (то есть просто исключить из выборки).
Медиану найти тоже очень просто. Это значение, которое больше одной половины выборки и меньше другой половины выборки. Соответственно, чтобы ее найти, мы выстраиваем все элементы выборки по неубыванию (то есть по возрастанию, а одинаковые значения пусть просто стоят рядом) — получится вариационный ряд. То, что будет ровно в его центре — медиана. Например, вот вариационный ряд:
1, 2, 2, 4, 7, 8, 11
Ровно посередине стоит число 4, оно и будет медианой. Если в вариационном ряду четное число элементов, то медианой будет среднее арифметическое двух чисел, которые стоят посередине. Скажем, в вариационном ряду:
3, 3, 3, 4, 4, 4 — медиана равна 3,5. Ничего сложного. Однако если в выборке много единиц наблюдения, быстро найти ее середину бывает проблематично. Порядковый номер медианы в вариационном ряду можно быстро определить вот так: i_{m} = \frac{n+1}{2} — всего-навсего прибавить к числу элементов выборки единицу и поделить на два; для ряда с четным количеством значений результат получится дробным, например 2,5 — тогда медианой будет среднее арифметическое этих, в данном случае 2-го и 3-го, значений. Медиану стоит использовать и в тех случаях, когда среднее арифметическое смысла не имеет (например, нелепо считать среднее арифметическое для стадии гипертонической болезни).
Осталась мода. А мода - это самое часто встречающееся в выборке значение. Она интересует нас меньше, поскольку обычно является наименее информативным показателем из трех (при работе с непрерывными данными, особенно если значения в выборке довольно точные, — скажем, пара знаков после запятой, — часто оказывается, что моды нет вообще или она совершенно не показательна). Но это очень полезный показатель, когда распределение имеет больше одной области, где данные концентрируются — скоро вы увидите пример бимодального распределения. В таких случаях можно указать несколько мод.
Если среднее арифметическое и медиана примерно равны, а гистограмма симметрична и напоминает колокол, то такое распределение называют нормальным. Если один бок колокола становится более пологим, чем другой, распределение называют скошенным. На этом рисунке первое распределение — нормальное, а остальные — все более скошенные.
print('fig = plt.figure(figsize=(24, 12))
#Лучше сделаю через цикл, чтобы сто раз не писать:
for i in range(0, 4):
fig.add_subplot(2, 2, i+1)
#Просто вытаскиываю выборки объемом сто тысяч элементов из распределений разной степени скошенности
sample = stats.skewnorm.rvs(i*3, size=100000)
x = np.linspace(stats.skewnorm.ppf(0.01, i), stats.skewnorm.ppf(0.99, i), 100)
#И отрисовываю гистограммы
plt.plot(x, stats.skewnorm.pdf(x, i*3), 'r-', lw=4, alpha=0.6)
plt.hist(sample, 16, normed=True, alpha=0.8)
#Ну и да, линии для медианы и среднего арифметического:
line1 = plt.axvline(np.mean(sample), label=('Ср. арифм. = '+str(round(np.mean(sample), 2))), color='c', linestyle='dashed', linewidth=1.5)
line2 = plt.axvline(np.median(sample), label = ('Медиана = '+str(round(np.median(sample), 2))), color='r', linestyle='dashed', linewidth=1.5)
#Подписи добавить:
plt.legend(handles=)
plt.show()')
Как видите, медиана лучше подходит для описания центральной тенденции в скошенных распределениях, поскольку не так сильно съезжает к тяжелому хвосту, как среднее арифметическое.
Меры рассеяния
С основными мерами центральной тенденции — медианой, модой и средним арифметическим — мы разобрались. Теперь посмотрим на меры рассеяния. Нужны они, повторюсь, затем, чтобы описать, насколько сильно данные рассеяны вдоль оси x. Сравните две гистограммы (обратите внимание, что ось x на них имеет одинаковый масштаб):
#Есть два распределения:
norm_rv1 = stats.norm(loc=0, scale = 1)
norm_rv2 = stats.norm(loc=0, scale = 5)
#Вытаскиываю из них две выборки по десять тысяч единиц наблюдения:
sample_norm1 = norm_rv1.rvs(size=10000)
sample_norm2 = norm_rv2.rvs(size=10000)
#Ну тут все ясно, две гистограммы:
fig = plt.figure(figsize=(16, 7))
ax1 = plt.subplot(211)
plt.hist(sample_norm1, 12, alpha=0.8)
ax2 = plt.subplot(212, sharex=ax1)
plt.hist(sample_norm2, 20, alpha=0.8)
plt.show()
Как видите, в обоих случаях данные распределены вдоль оси x вокруг нуля, но вариабельность значений второй выборки выше. Там, очевидно, есть и объекты со значениями меньше -10 и больше 10, а в первой выборке таких нет вовсе, посмотрите на гистограмму.
Для описания рассеивания значений в более-менее симметричных колоколообразных распределениях используют среднеквадратическое отклонение, обозначается оно буквой σ (сигма). Для первой выборки σ=1, а для второй выборки σ=5, поэтому вторая гистограмма шире. Чем больше среднеквадратическое отклонение, тем сильнее данные распростираются в стороны от медианы и среднего арифметического. Но не стоит судить об этом показателе по внешнему виду гистограммы. На следующем рисунке две гистограммы одной и той же выборки, просто сделанные с разным масштабом. Обратите внимание на значения на оси x.
fig = plt.figure(figsize=(12, 4))
#Это вот чтобы масштаб был разным:
from matplotlib import gridspec
gs = gridspec.GridSpec(1, 2, width_ratios=)
#Две гистограммы для одной выборки:
sample_norm = np.random.randn(1000)
ax1 = fig.add_subplot(gs)
plt.hist(sample_norm, alpha=0.8)
plt.xticks(range(-3, 4))
ax2 = fig.add_subplot(gs)
plt.hist(sample_norm, bins=16, alpha=0.8)
plt.show()
Считается среднеквадратическое отклонение по вот такой формуле: \sigma = \sqrt{\frac{\sum_{i=1}^n (x_i - \overline{x})^2}{n}} — на первый взгляд довольно пугающей, но на самом деле очень простой. Смотрите на рисунок:
#Тут, конечно, ад. Надо было в пейнте нарисовать.
#Ну да ладно, берем четыре значения, считаем среднее.
a, b, c, d = 0, 1, 4, 7
mean = (a+b+c+d)/4
#И пусть будет прозрачненько.
alpha=0.3
#Рисунок будет крупным:
fig = plt.figure(figsize=(15, 8))
#Все сделано по одной схеме. Обозначаем место, где на картинке будет график:
ax1 = fig.add_subplot(222)
#Рисуем ось х и штриховую линию для среднего:
plt.axhline(0., color='k')
plt.axvline(mean, color='r', linestyle='dashed')
#И добавляем сколько надо квадратов. Рисуются они как-то вот так:
ax1.add_patch(
patches.Rectangle(
(a, 0.), # (x,y)
mean-a, # ширина
mean-a, # высота
alpha=alpha)
)
ax1.add_patch(patches.Rectangle((b, 0.), mean-b, mean-b, alpha=alpha, color='r'))
ax1.add_patch(patches.Rectangle((mean, 0.), c-mean, c-mean, alpha=alpha, color='c'))
ax1.add_patch(patches.Rectangle((mean, 0.), d-mean, d-mean, alpha=alpha, color='g'))
ax1.plot(, , ls='', marker='o', ms=8, mfc='r')
ax2 = fig.add_subplot(221, sharey=ax1)
plt.axhline(0., color='k')
plt.axvline(mean, color='r', linestyle='dashed')
plt.plot(, , ls='', marker='o', ms=8, mfc='r')
ax3 = fig.add_subplot(223)
sigma = np.std(np.array())
plt.axhline(0., color='k')
plt.axvline(mean, color='r', linestyle='dashed')
plt.plot(, , ls='', marker='o', ms=8, mfc='r')
ax3.add_patch(patches.Rectangle((a, 0.), mean-a, mean-a, alpha=alpha))
ax3.add_patch(patches.Rectangle((b, 0.), mean-b, mean-b, alpha=alpha, color='r'))
ax3.add_patch(patches.Rectangle((mean, 0.), c-mean, c-mean, alpha=alpha, color='c'))
ax3.add_patch(patches.Rectangle((mean, 0.), d-mean, d-mean, alpha=alpha, color='g'))
ax3.add_patch(patches.Rectangle((mean, 0.), sigma, sigma, alpha=alpha, color='m'))
ax3.plot(, , ls='', marker='o', ms=8, mfc='r')
ax4 = fig.add_subplot(224, sharey=ax3)
plt.axhline(0., color='k')
plt.axvline(mean, color='r', linestyle='dashed')
plt.plot(, , ls='', marker='o', ms=8, mfc='r')
ax4.add_patch(patches.Rectangle((mean, 0.), sigma, sigma, alpha=alpha, color='m'))
#А вот этот нелепый кусок кода для добавления подписи с сигмой к квадрату.
ax4.annotate('σ', xy=((mean+(mean+sigma))/2, sigma+0.2), \
xytext=((mean+(mean+sigma))/2, sigma+0.6), xycoords='data',
fontsize=18, ha='center', va='bottom',
bbox=dict(boxstyle='square', fc='white'),
arrowprops=dict(arrowstyle='-
Давайте сделаем все буквально по формуле. Возьмем для наглядной иллюстрации выборку из четырех элементов: 0,1,4,7. Среднее арифметическое этой выборки равно 3. Покажем это на первом графике: на оси x красными точками отмечены элементы нашей выборки, а вертикальная штриховая линия — среднее арифметическое.
Каждый элемент выборки отклоняется от среднего на некоторое расстояние. Например, 0 отклоняется от среднего на 3, а 1 отклоняется от среднего на 2. На втором рисунке мы достроили эти отклонения до квадратов. Вот смотрите, у зеленого квадрата сторона равна 4, потому что 7 отклоняется от среднего на 4.
На третьем рисунке появился сиреневый квадрат. Получился он так: мы просто нашли среднюю площадь первых четырех квадратов и построили пятый с такой площадью.
На четвертом рисунке видно, что сторона сиреневого квадрата равна среднеквадратическому отклонению.
Итак, среднеквадратическое отклонение — это среднее арифметическое квадратов абсолютных отклонений от выборочного среднего.
При этом, кстати говоря, площадь сиреневого квадрата является еще одной мерой рассеяния, дисперсией (среднее арифметическое квадратов отклонений от выборочного среднего). Поэтому вряд ли вас удивит, что обозначают дисперсию обычно вот так: \sigma^2, а высчитывается она следующим образом: \sigma^2=\frac{\sum_{i=1}^n (x_i - \overline{x})^2}{n} — ровно так же, как среднеквадратическое отклонение, только без корня.
Да, трудно понять, что такое дисперсия. И что такое среднеквадратическое отклонение. Их нельзя пощупать, даже представить довольно трудно. Просто думайте так: чем они больше, тем выше вариабельность значений в выборке, тем шире гистограмма.
Еще одна мера рассеяния — коэффициент вариации. Собственно, это то же самое среднеквадратическое отклонение, только выраженное в процентах от выборочного среднего. То есть чтобы его найти, нужно разделить σ на среднее и умножить на 100%. Коэффициент вариации удобно использовать для сравнения рассеяний разных выборок, ведь он, повторюсь, выражается в процентах.
Любопытная деталь, связанная с дисперсией и среднеквадратическим отклонением: посчитанные по выборке, они обычно немного меньше, чем в генеральной совокупности. И чем выборка меньше, тем хуже она отражает истинный разброс значений. Чтобы получить представление об этих показателях для генеральной совокупности, применяют так называемую поправку Бесселя. Она очень проста и представляет собой единицу с минусом в знаменателях обеих формул. То есть для оценки дисперсии и среднеквадратического отклонения в генеральной совокупности формулы станут вот такими:
s = \sqrt{\frac{\sum_{i=1}^n (x_i - \overline{x})^2}{n-1}}, s^2=\frac{\sum_{i=1}^n (x_i - \overline{x})^2}{n-1}
Такое скорректированное среднеквадратическое отклонение называют стандартным отклонением. Эту формулу следует использовать, когда вы не можете включить в свою выборку всю генеральную совокупность, то есть почти всегда. Но, очевидно, чем больше размер вашей выборки, тем меньший вклад будет вносить поправка.
Следующий занимательный и очень важный момент, связанный со стандартным отклонением, заключается в том, что для нормального распределения выполняется следующее правило: 68,27% единиц наблюдения попадают в интервал от μ−1σ до μ+1σ; 95,45% единиц наблюдения попадают в интервал от μ−2σ до μ+2σ; 99,7% единиц наблюдения попадают в интервал от μ−3σ до μ+3σ, где μ(мю) — генеральное среднее арифметическое. Правило называется, соответственно, правилом 68-95-99,7. Эти значения можно запомнить. Но можно и не запоминать: скоро я объясню, откуда взялись они взялись, как их получить и как оценить долю единиц наблюдения между двумя любыми точками, а не только между целыми значениями σ. Это необходимо для проверки гипотез и построения доверительных интервалов.
norm_rv = stats.norm(loc=0, scale = 1) #Стандартное нормальное распределение
fig_rule = plt.figure(figsize=(14, 5))
ax = fig_rule.add_subplot(111)
x = np.linspace(-4,5,100)
pdf = norm_rv.pdf(x) #График функции плотности вероятности
xticks =
xticks_s =
plt.xticks(xticks, xticks_s)
plt.ylim(0, 0.8)
plt.plot(x, pdf)
for i in : plt.axvline(i, ls='dotted', c='#2ecc71', ymax=0.9)
for i in : plt.axvline(i, ls='dotted', c="#34495e", ymax=0.78)
for i in : plt.axvline(i, ls='dotted', c="#e74c3c", ymax=0.66)
#Снова проклятые неудобные подписи:
ax.annotate('68%', xy=(-1, 0.5), \
xytext=(1, 0.5), xycoords='data',
fontsize=15, ha='left', va='center',
bbox=dict(boxstyle='square', fc='white'),
arrowprops=dict(arrowstyle='<->, head_length=0.4,head_width=0.2', lw=1.0, color="#e74c3c"))
ax.annotate('95%', xy=(-2, 0.6), \
xytext=(2, 0.6), xycoords='data',
fontsize=15, ha='left', va='center',
bbox=dict(boxstyle='square', fc='white'),
arrowprops=dict(arrowstyle='<->, head_length=0.4,head_width=0.2', lw=1.0, color="#34495e"))
ax.annotate('99,7%', xy=(-3, 0.7), \
xytext=(3, 0.7), xycoords='data',
fontsize=15, ha='left', va='center',
bbox=dict(boxstyle='square', fc='white'),
arrowprops=dict(arrowstyle='<->, head_length=0.4,head_width=0.2', lw=1.0, color='#2ecc71'))
plt.show()
Итак, для описания нормального распределения удобно указать среднее арифметическое и стандартное отклонение. Но для описания распределений, которые не похожи на нормальное, применять стандартное отклонение нельзя, ведь оно не отражает асимметрию. Для этого можно использовать квартили.
Как вы помните, если взять выборку и выстроить ее по возрастанию (получится вариационный ряд), а потом посмотреть, что находится посередине, то вы найдете медиану.
А если взять куски вариационного ряда слева и справа от медианы и для каждого снова найти медиану, это будут первый и третий квартили. То есть 1/4 часть выборки меньше первого квартиля, половина меньше второго (медианы), и 3/4 меньше третьего. Иными словами, квартили — это значения, которые делят вариационный ряд на четыре части. Расстояние между первым и третьим квартилями называется интерквартильным размахом — это еще одна мера рассеяния.
Помимо квартилей, есть еще и процентили, только они дробят вариационный ряд на сто кусков. Первый квартиль - это то же самое, что 25-й процентиль, а третий квартиль - это то же самое, что 75-й процентиль. Квартили и процентили это разновидности квантилей — значений, которые дробят вариационный ряд на некоторое количество частей. Но об остальных знать не обязательно.
#Тут ничего нового. Обычные прямоугольники.
fig = plt.figure(figsize=(15, 8))
ax = fig.add_subplot(111)
for p in :
ax.add_patch(p)
for x in range(0, 8, 2):
ax.annotate('25% данных', xy=(0.2+x/10, 0.25), \
xytext=(0.2+x/10, 0.25), xycoords='data',
fontsize=18, ha='center', va='center',
bbox=dict(boxstyle='round', fc='white'))
ax.annotate('Медиана\n(2-й квартиль)', xy=(0.5, 0.4), \
xytext=(0.5, 0.6), xycoords='data',
fontsize=18, ha='center', va='center',
bbox=dict(boxstyle='square', pad=0.3, fc='white'),
arrowprops=dict(arrowstyle=']-, widthA=4.2, lengthA=0.1', lw=1.5))
ax.annotate('1-й квартиль', xy=(0.3, 0.4), \
xytext=(0.3, 0.6), xycoords='data',
fontsize=18, ha='center', va='top',
bbox=dict(boxstyle='square', pad=0.3, fc='white'),
arrowprops=dict(arrowstyle=']-, widthA=3.8, lengthA=0.1', lw=1.5))
ax.annotate('3-й квартиль', xy=(0.7, 0.4), \
xytext=(0.7, 0.6), xycoords='data',
fontsize=18, ha='center', va='top',
bbox=dict(boxstyle='square', pad=0.3, fc='white'),
arrowprops=dict(arrowstyle=']-, widthA=3.8, lengthA=0.1', lw=1.5))
plt.ylim(0, 0.8)
plt.axis('off')
plt.show()
Давайте посмотрим, как квартили выглядят на скошенной гистограмме. Мы ожидаем, что первый и третий квартили будут на разном расстоянии от медианы, ведь именно эта фишка и заставляет нас их использовать.
plt.figure(figsize=(12, 6))
plt.subplot(211)
i = 6
sample = stats.skewnorm.rvs(i, size=10000)
x = np.linspace(stats.skewnorm.ppf(0.01, i), stats.skewnorm.ppf(0.99, i), 100)
plt.plot(x, stats.skewnorm.pdf(x, i), 'r-', lw=4, alpha=0.6)
plt.hist(sample, 16, normed=True, alpha=0.7)
l1 = plt.axvline(np.percentile(sample, 25), c="#e74c3c", alpha=0.8, linestyle='dashed', label='1-й квартиль')
l2 = plt.axvline(np.percentile(sample, 50), c="#34495e", alpha=0.8, linestyle='dashed', label='Медиана')
l3 = plt.axvline(np.percentile(sample, 75), c="#2ecc71", alpha=0.8, linestyle='dashed', label='3-й квартиль')
plt.legend(handles=)
plt.subplot(212)
i = -6
sample = stats.skewnorm.rvs(i, size=10000)
x = np.linspace(stats.skewnorm.ppf(0.01, i), stats.skewnorm.ppf(0.99, i), 100)
plt.plot(x, stats.skewnorm.pdf(x, i), 'r-', lw=4, alpha=0.6)
plt.hist(sample, 16, normed=True, alpha=0.7)
l1 = plt.axvline(np.percentile(sample, 25), c="#e74c3c", alpha=0.8, linestyle='dashed', label='1-й квартиль')
l2 = plt.axvline(np.percentile(sample, 50), c="#34495e", alpha=0.8, linestyle='dashed', label='Медиана')
l3 = plt.axvline(np.percentile(sample, 75), c="#2ecc71", alpha=0.8, linestyle='dashed', label='3-й квартиль')
plt.legend(handles=)
plt.show()
Смотрите, в первом случае правый склон гистограммы более пологий, и третий квартиль дальше от медианы, чем первый. Во втором случае первый квартиль дальше от медианы, чем третий, и у гистограммы более пологий левый склон.
Не всегда легко определить, что лучше использовать для описания рассеяния — стандартное отклонение или квартили. Я писал, что стандартное отклонение подходит для более-менее симметричных колоколообразных распределениях. Но что такое более-менее? И когда распределение становится слишком скошенным? На этот случай нет четкого правила, так что решайте сами. В статистике так бывает.
Ящик с усами
Есть один наигоднейший график, очень удобный для того, чтобы быстро оценить параметры выборки и сравнить их для нескольких выборок. Называется он ящиком с усами, или диаграммой размаха, или ящичковой диаграммой Тьюки.
Выглядят ящики с усами вот так:
plt.figure(figsize=(9, 6))
#Бокс-плоты из документации Seaborn -- каноничный датасет ирисов.
iris = sns.load_dataset("iris")
ax = sns.boxplot(data=iris, orient="v", palette="Set2")
ax.get_xaxis().set_visible(False)
ax.get_yaxis().set_visible(False)
Ориентированы ящики могут быть как горизонтально, так и вертикально, на усмотрение исследователя. Суть их при этом не меняется.
Устроены эти диаграммы довольно незатейливо. Границы самого ящика показывают первый и третий квартили, линия в центре — медиана. Максимальная длина каждого уса — полтора интерквартильных размаха. Если данные отстоят в меньшую сторону от первого квартиля или в большую от третьего на расстояние более полутора интерквартильных размахов (расстояние между 1-м и 3-м квартилями), то эти данные считаются выбросами и изображаются в виде точек. А если, например, единица наблюдения с наименьшим значением отстоит от первого квартиля только на половину интерквартильного размаха, то ус закончится на этой единице наблюдения (сравните оранжевый ящик на рисунке сверху и синий рядом с ним; длина усов оранжевого ящика максимальна — полтора интерквартильных размаха, а у синего ящика меньше, потому что данные так далеко не распростираются).
plt.figure(figsize=(15, 7))
plt.subplot(221)
i = 6
sample = stats.skewnorm.rvs(i, size=1000)
x = np.linspace(stats.skewnorm.ppf(0.01, i), stats.skewnorm.ppf(0.99, i), 100)
plt.plot(x, stats.skewnorm.pdf(x, i), 'r-', lw=4, alpha=0.6)
plt.hist(sample, 16, normed=True, alpha=0.7)
l1 = plt.axvline(np.percentile(sample, 25), c="#e74c3c", alpha=0.8, linestyle='dashed', label='1-й квартиль')
l2 = plt.axvline(np.percentile(sample, 50), c="#34495e", alpha=0.8, linestyle='dashed', label='Медиана')
l3 = plt.axvline(np.percentile(sample, 75), c="#2ecc71", alpha=0.8, linestyle='dashed', label='3-й квартиль')
plt.legend(handles=)
plt.subplot(223)
sns.boxplot(sample, width=0.5)
plt.subplot(222)
i = -6
sample = stats.skewnorm.rvs(i, size=1000)
x = np.linspace(stats.skewnorm.ppf(0.01, i), stats.skewnorm.ppf(0.99, i), 100)
plt.plot(x, stats.skewnorm.pdf(x, i), 'r-', lw=4, alpha=0.6)
plt.hist(sample, 16, normed=True, alpha=0.7)
l1 = plt.axvline(np.percentile(sample, 25), c="#e74c3c", alpha=0.8, linestyle='dashed', label='1-й квартиль')
l2 = plt.axvline(np.percentile(sample, 50), c="#34495e", alpha=0.8, linestyle='dashed', label='Медиана')
l3 = plt.axvline(np.percentile(sample, 75), c="#2ecc71", alpha=0.8, linestyle='dashed', label='3-й квартиль')
plt.legend(handles=)
plt.subplot(224)
sns.boxplot(sample, width=0.5)
plt.show()
Важно помнить, что не каждую выборку можно описать с помощью перечисленных показателей или отразить ящичковой диаграммой. Хороший тому пример — бимодальные распределения. Смотрите на картинку.
#Это функция для получения бимодально распределенной величины.
def bimodal( low1, high1, mode1, low2, high2, mode2 ):
toss = np.random.choice((1, 2))
if toss == 1:
return np.random.triangular( low1, high1, mode1 )
else:
return np.random.triangular( low2, high2, mode2 )
a =
for x in range(1000):
a.append(bimodal(1, 4, 4, 3, 6, 8))
fig = plt.figure(figsize=(12, 5))
ax1 = plt.subplot(211)
sns.distplot(a)
ax2 = plt.subplot(212, sharex=ax1)
sns.boxplot(a)
plt.show()
В подобных случаях разумным выходом будет просто показать гистограмму.
Пример описания выборки
Давайте для примера опишем мужчин и женщин из таблицы Гальтона, о которой мы говорили при обсуждении гистограмм.
fig = plt.figure(figsize=(12, 4))
fig.text(0.5, 0.04, "Рост", ha='center')
fig.text(0.08, 0.5, "Число людей", va='center', rotation='vertical')
ax1 = plt.subplot(121)
galton.Height.hist(bins=14, alpha=0.8)
plt.title('Рост женщин')
ax2 = plt.subplot(122)
galton.Height.hist(bins=12, alpha=0.8)
plt.title('Рост мужчин')
plt.show()
Гистограмма роста женщин достаточно похожа на нормальное распределение, так что для иллюстрации опишем первую группу с помощью среднего арифметического и стандартного отклонения. Вторую группу опишем с помощью медианы и квартилей.
df = pd.DataFrame(columns=, index=galton.describe().index)
df = galton.describe().Height
df = galton.describe().Height
df.ix
Рост испытуемых женского пола варьирует в пределах от 56 до 70 дюймов, среднее арифметическое составляет 64.1 дюйма (стандартное отклонение 2.4 дюйма).
Рост испытуемых мужского пола варьирует в пределах от 60 до 79 дюймов с медианой, равной 69.2 дюйма (интерквартильная широта от 67.5 до 71 дюйма).
sns.boxplot(x=galton.Gender, y=galton.Height, data=galton, width=0.4)
plt.ylabel('Рост')
plt.xlabel('Пол')
plt.show()
Многие исследователи, чтобы обозначить среднее арифметическое и стандартное отклонение, пишут так: 64.1±2.4, но делать этого не стоит, потому что при таком оформлении возникает путаница между стандартным отклонением и доверительными интервалами (о доверительных интервалах скоро тоже поговорим). Лучше вообще не используйте символ "±".
Кажется, вьетнамец победил, он находится ближе к правому краю колокола, чем серб. Как в этом убедиться?
Можно посчитать, на сколько стандартных отклонений рост серба и рост вьетнамца отличаются от средних значений. Это самое количество стандартных отклонений называется z-оценкой и вычисляется по такой формуле:
z=\frac{x-\mu}{\sigma}
То есть все, что надо сделать в случае с вьетнамцем — это вычесть 162 из 174 и поделить результат на 5,39.
Итак, рост вьетнамца отклоняется на 2,23 стандартного отклонения от среднего, а рост серба — на 1,78 стандартного отклонения от среднего, то есть вьетнамец у себя на родине воспринимается более высоким, чем серб у себя.
Таким образом можно стандартизировать любую нормально распределенную величину. По сути, мы наносим ее на график нормального распределения со средним, равным нулю и стандартным отклонением, равным единице.
x = np.linspace(-6,5,100)
pdf = norm_rv.pdf(x)
xticks =
plt.xticks(xticks)
plt.plot(x, pdf)
l1=plt.axvline(viet, color='g', ls='dashed', label='z-оценка\nроста вьетнамца')
l2=plt.axvline(serb, color='b', ls='dashed', label='z-оценка\nроста серба')
plt.legend(handles=)
plt.show()
Такое нормальное распределение со средним, равным 0, и среднеквадратическим отклонением, равным 1, называется стандартным нормальным распределением, или z-распределением.
У него есть одна любопытная и важная особенность. Зная z-оценку единицы наблюдения, мы можем также узнать, какой процент наблюдений превосходит значение этой единицы наблюдения. То есть мы можем прямо сейчас оценить, какая доля вьетнамцев выше или ниже ростом, чем наш вьетнамец. Для этого существует специальная таблица z-значений.
Z-оценка роста нашего вьетнамского приятеля составила 2,23 (с округлением до второго знака после запятой). В вертикальном столбце под буквой z находим 2,2, а в горизонтальном 0,03 (потому что 2,2 + 0,03 = 2,23), на пересечении стоит 0,9871. Это означает, что он выше 98,7% единиц наблюдения.
То есть если согнать всех-всех вьетнамцев на поляну и сделать из них гистограмму таким образом:
#International Statistical Review
#1975 vol: 43 (3) pp: 339
— только 2,2% людей справа будут выше него.
Для вычисления z-значения не обязательно использовать таблицу, можно получить более точную оценку с использованием специальных инструментов для работы с данными (это не обязательно должны быть R, SPSS или что-то подобное, можно поставить превосходное приложение "Probability Distributions" на телефон).
Зная про z-оценку, нетрудно догадаться, откуда взялось правило 68-95-99,7. Давайте посмотрим по таблице, какая доля единиц наблюдения лежит в интервале от -1 до 1.
Левее z=1 лежит 84,13% единиц наблюдения, левее z=-1 лежит 15,87% единиц наблюдения. Между ними, соответственно, 84,13-15,87 = 68,26%.
Можете проверить для 95 и 99,7, ну или попробовать для любых других точек.
#Этот график я уже строил выше, так что не буду копировать тонну кода и просто перетащу сам график.
fig_rule
Подведу итог. Z-оценка — это количество стандартных отклонений, на которое единица наблюдения отклоняется от среднего арифметического. С помощью таблицы z-оценок можно определить, какая доля единиц наблюдения лежит выше или ниже произвольного значения. Применима z-оценка только для нормально распределенных величин.
Что будет, если я извлеку много-много выборок из генеральной совокупности и посчитаю для каждой из этих выборок среднее арифметическое? Видимо, получится выборка из средних арифметических. Но как будет выглядеть гистограмма такой выборки?
Давайте проделаем этот трюк для каждой из наших генеральных совокупностей. Вытащим из каждой по десять тысяч выборок по 3 единицы наблюдения, посчитаем средние арифметические для всех них и получим две гистограммы, показывающие распределения средних арифметических.
#Вот так вытаскиваю выборки:
mean_list_norm = np.array()
mean_list_skewed = np.array()
fig = plt.figure(figsize=(12, 5))
fig.add_subplot(121)
plt.axvline(mean_list_norm.mean(), color='b', linestyle='dashed', linewidth=1)
plt.legend()
plt.hist(mean_list_norm, bins = 20)
fig.add_subplot(122)
plt.hist(mean_list_skewed, bins = 20)
plt.axvline(mean_list_skewed.mean(), color='b', linestyle='dashed', linewidth=1)
plt.legend()
plt.show()
Как видите, распределения выборочных средних похожи на распределения величин в генеральных совокупностях (хотя правое скошено несколько меньше). Средние арифметические примерно совпадают с истинными средними арифметическими, что неудивительно. Бросается в глаза, что эти гистограммы получились уже, чем гистограммы генеральных совокупностей (обратите внимание на ось x. Это тоже неудивительно, ведь, понятное дело, выборочные средние не так сильно варьируют, как единицы наблюдения во всей генеральной совокупности.
Давайте провернем этот трюк еще раз, но теперь будет извлекать выборки объемом 15.
mean_list_norm = np.array()
mean_list_skewed = np.array()
fig = plt.figure(figsize=(12, 5))
fig.add_subplot(121)
plt.axvline(mean_list_norm.mean(), color='b', linestyle='dashed', linewidth=1)
plt.legend()
plt.hist(mean_list_norm, bins = 20)
fig.add_subplot(122)
plt.hist(mean_list_skewed, bins = 20)
plt.axvline(mean_list_skewed.mean(), color='b', linestyle='dashed', linewidth=1)
plt.legend()
plt.show()
Гистограммы стали еще уже, а правая еще менее скошенной. Давайте еще раз. Извлечем из генеральных совокупностей по 10000 выборок, но теперь объемом 100.
mean_list_norm = np.array()
mean_list_skewed = np.array()
fig = plt.figure(figsize=(12, 5))
fig.add_subplot(121)
plt.axvline(mean_list_norm.mean(), color='b', linestyle='dashed', linewidth=1)
plt.legend()
plt.hist(mean_list_norm, bins = 20)
fig.add_subplot(122)
plt.hist(mean_list_skewed, bins = 20)
plt.axvline(mean_list_skewed.mean(), color='b', linestyle='dashed', linewidth=1)
plt.legend()
plt.show()
На вид две гистограммы нормально распределенной величины. Кстати, смотрите, они стали еще уже. Это абсолютно логично: чем больше размер выборки, тем более точную оценку среднего она может нам предоставить. Маленькие выборки хуже отражают параметры генеральной совокупности, чем большие.
И, как видите, чем больше выборка, тем сильнее распределение выборочных средних похоже на нормальное, и тем меньше будет рассеяние значений вокруг среднего (выше точность оценки среднего). Стало быть, чем больше объем выборки, тем меньше стандартное отклонение в выборке выборочных средних :)
Вот это стандартное отклонение в выборке, составленной из выборочных средних, называется стандартной ошибкой среднего. Приятная новость состоит в том, что для его определения не надо получать кучу выборок, есть простая маленькая формула:
SD_\overline{x} = \frac{\sigma}{\sqrt{n}} — истинная стандартная ошибка среднего.
Но, как вы понимаете, для оценки истинной стандартной ошибки среднего надо знать среднеквадратическое отклонение в генеральной совокупности, а это редкая роскошь, которой мы обычно не обладаем. Так что стандартную ошибку среднего мы оцениваем по выборочному стандартному отклонению:
\Huge SE_\overline{x} = \frac{s}{\sqrt{n}}
Некоторые исследователи используют стандартную ошибку среднего для описания выборки вместо стандартного отклонения (первая всегда меньше второй, поэтому иногда возникает такой соблазн). Этого делать не стоит, так как этот показатель нужен для оценки среднего в генеральной совокупности, он не является описательной статистикой.
Одной только стандартной ошибки для оценки среднего недостаточно. Но она позволит нам вычислять доверительные интервалы и проверять гипотезы, об этом мы сейчас поговорим.
Когда мы хотим найти среднее генеральной совокупности, нас почти всегда ждет горькое разочарование, ведь генеральные совокупности обычно огромны. Так что мы считаем выборочное среднее и надеемся, что оно не сильно отличается от генерального. Меру точности нашей оценки можно отразить с помощью доверительного интервала.
Мы знаем, что если извлечь из генеральной совокупности множество выборок и для каждой посчитать среднее арифметическое, полученные значения распределятся по нормальному закону со стандартным отклонением, равным \frac{\sigma}{\sqrt{n}}. Еще мы знаем, что, будь нам известно истинное среднее, то если бы мы отложили от него в обе стороны по три стандартных ошибки среднего, охватили бы примерно 99,7% всех возможных стандартных ошибок среднего (помните правило 68-95-99,7?).
Ну так вот, это работает и в другую сторону. Мы, конечно, не знаем истинного среднего, но запросто можем получить выборочное среднее. И можем от него в обе стороны отложить три стандартных ошибки среднего. И вот что интересно: если извлекать из генеральной совокупности выборки одну за другой много-много раз, а потом для каждой считать выборочное среднее и для него строить такой интервал: \overline{x} \pm 3SE_\overline{x}, то примерно 99,7% интервалов будут включать в себя истинное среднее.
Итак, когда вы находите выборочное среднее, а потом считаете, каким значениям соответствуют \overline{x} - 3SE_\overline{x} и \overline{x} + 3SE_\overline{x}, это и есть построение доверительного интервала. Совсем не обязательно откладывать именно по три стандартных ошибки среднего влево и вправо. Мы можем, например, захотеть построить 95%-й доверительный интервал. И мы знаем, что для этого надо отступать на две стандартных ошибки среднего, но это не вполне точное значение. Чтобы узнать более точное, надо заглянуть в таблицу z-оценок. Количество стандартных отклонений, которое нам требуется, добываем так: 1-(100-ДИ)/200, и полученное значение ищем в таблице, но не по индексам, а в куче z-значений.
Давайте со мной для 95%-го доверительного интервала. 1-(100-95)/200 = 1-0,025 = 0,975. Ищем 0,975 в таблице, и обнаруживаем это число на пересечении строки 1,9 и столбца 0,06. Чудно, это значит, что для построения нужного нам интервала надо отложить от выборочного среднего по 1,96 стандартного отклонения в обе стороны. Соответственно, формула 95%-го ДИ:
\overline{x} \pm 1,96SE_\overline{x}
Итак, 95% построенных таким образом интервалов будут содержать истинное среднее. В качестве иллюстрации я сгенерировал генеральную совокупность, извлек из нее 100 выборок и для каждой построил 95%-й доверительный интервал для среднего. Те интервалы, которые не включили в себя истинное среднее, отмечены жирным.
#Тут нужна специальная библиотека для статистики:
from statsmodels.stats.weightstats import zconfint
population = np.random.randn(100000)
plt.figure(figsize=(10, 10))
plt.xlim(xmin=-1)
for i in range(100):
sample = np.random.choice(population, size=60)
interval = zconfint(sample)
#Проверяю, включает ли интервал среднее, и если нет, то рисую его жирным.
if interval < 0 and interval > 0:
plt.plot((interval, interval), (i, i), color='#9A2EFE', alpha=0.4)
plt.scatter(x=, y=i, color='#9A2EFE', alpha=0.5, s=10)
else:
plt.plot((interval, interval), (i, i), color='#9A2EFE', alpha=1)
plt.scatter(x=, y=i, color='#9A2EFE', alpha=1, s=15)
plt.axvline(0, c='#04B45F', alpha=0.3, ls='dashed')
plt.show()
Ориентированы ящики могут быть как горизонтально, так и вертикально, на усмотрение исследователя. Суть их при этом не меняется.
Устроены эти диаграммы довольно незатейливо. Границы самого ящика показывают первый и третий квартили, линия в центре — медиана. Максимальная длина каждого уса — полтора интерквартильных размаха. Если данные отстоят в меньшую сторону от первого квартиля или в большую от третьего на расстояние более полутора интерквартильных размахов (расстояние между 1-м и 3-м квартилями), то эти данные считаются выбросами и изображаются в виде точек. А если, например, единица наблюдения с наименьшим значением отстоит от первого квартиля только на половину интерквартильного размаха, то ус закончится на этой единице наблюдения (сравните оранжевый ящик на рисунке сверху и синий рядом с ним; длина усов оранжевого ящика максимальна — полтора интерквартильных размаха, а у синего ящика меньше, потому что данные так далеко не распростираются).
plt.figure(figsize=(15, 7))
plt.subplot(221)
i = 6
sample = stats.skewnorm.rvs(i, size=1000)
x = np.linspace(stats.skewnorm.ppf(0.01, i), stats.skewnorm.ppf(0.99, i), 100)
plt.plot(x, stats.skewnorm.pdf(x, i), 'r-', lw=4, alpha=0.6)
plt.hist(sample, 16, normed=True, alpha=0.7)
l1 = plt.axvline(np.percentile(sample, 25), c="#e74c3c", alpha=0.8, linestyle='dashed', label='1-й квартиль')
l2 = plt.axvline(np.percentile(sample, 50), c="#34495e", alpha=0.8, linestyle='dashed', label='Медиана')
l3 = plt.axvline(np.percentile(sample, 75), c="#2ecc71", alpha=0.8, linestyle='dashed', label='3-й квартиль')
plt.legend(handles=)
plt.subplot(223)
sns.boxplot(sample, width=0.5)
plt.subplot(222)
i = -6
sample = stats.skewnorm.rvs(i, size=1000)
x = np.linspace(stats.skewnorm.ppf(0.01, i), stats.skewnorm.ppf(0.99, i), 100)
plt.plot(x, stats.skewnorm.pdf(x, i), 'r-', lw=4, alpha=0.6)
plt.hist(sample, 16, normed=True, alpha=0.7)
l1 = plt.axvline(np.percentile(sample, 25), c="#e74c3c", alpha=0.8, linestyle='dashed', label='1-й квартиль')
l2 = plt.axvline(np.percentile(sample, 50), c="#34495e", alpha=0.8, linestyle='dashed', label='Медиана')
l3 = plt.axvline(np.percentile(sample, 75), c="#2ecc71", alpha=0.8, linestyle='dashed', label='3-й квартиль')
plt.legend(handles=)
plt.subplot(224)
sns.boxplot(sample, width=0.5)
plt.show()
Важно помнить, что не каждую выборку можно описать с помощью перечисленных показателей или отразить ящичковой диаграммой. Хороший тому пример — бимодальные распределения. Смотрите на картинку.
#Это функция для получения бимодально распределенной величины.
def bimodal( low1, high1, mode1, low2, high2, mode2 ):
toss = np.random.choice((1, 2))
if toss == 1:
return np.random.triangular( low1, high1, mode1 )
else:
return np.random.triangular( low2, high2, mode2 )
a =
for x in range(1000):
a.append(bimodal(1, 4, 4, 3, 6, 8))
fig = plt.figure(figsize=(12, 5))
ax1 = plt.subplot(211)
sns.distplot(a)
ax2 = plt.subplot(212, sharex=ax1)
sns.boxplot(a)
plt.show()
В подобных случаях разумным выходом будет просто показать гистограмму.
Пример описания выборки
Давайте для примера опишем мужчин и женщин из таблицы Гальтона, о которой мы говорили при обсуждении гистограмм.
fig = plt.figure(figsize=(12, 4))
fig.text(0.5, 0.04, "Рост", ha='center')
fig.text(0.08, 0.5, "Число людей", va='center', rotation='vertical')
ax1 = plt.subplot(121)
galton.Height.hist(bins=14, alpha=0.8)
plt.title('Рост женщин')
ax2 = plt.subplot(122)
galton.Height.hist(bins=12, alpha=0.8)
plt.title('Рост мужчин')
plt.show()
Гистограмма роста женщин достаточно похожа на нормальное распределение, так что для иллюстрации опишем первую группу с помощью среднего арифметического и стандартного отклонения. Вторую группу опишем с помощью медианы и квартилей.
df = pd.DataFrame(columns=, index=galton.describe().index)
df = galton.describe().Height
df = galton.describe().Height
df.ix
| F | M | |
| mean | 64.110162 | 69.228817 |
| std | 2.370320 | 2.631594 |
| min | 56.000000 | 60.000000 |
| 25% | 62.500000 | 67.500000 |
| 50% | 64.000000 | 69.200000 |
| 75% | 65.500000 | 71.000000 |
| max | 70.500000 | 79.000000 |
Многие исследователи, чтобы обозначить среднее арифметическое и стандартное отклонение, пишут так: 64.1±2.4, но делать этого не стоит, потому что при таком оформлении возникает путаница между стандартным отклонением и доверительными интервалами (о доверительных интервалах скоро тоже поговорим). Лучше вообще не используйте символ "±".
Z-оценка
Встречаются как-то раз серб и вьетнамец. Рост серба 194 см, и у себя на родине он считается весьма высоким. Рост вьетнамца 174 см, и у себя на родине он тоже считается весьма высоким. Кто выше: серб среди сербов или вьетнамец среди вьетнамцев?
В Сербии средний рост мужчины 182 см. Так что рост нашего сербского господина на 12 см выше среднего.
Во Вьетнаме средний рост мужчины 162 см. Так что рост нашего вьетнамского господина тоже на 12 см выше среднего.
Стандартное отклонение для вьетнамцев — 5,39 см, а для сербов — 6,74 см (посчитаны по огромным выборкам, так что хорошо аппроксимируют среднеквадратические отклонения). Давайте посмотрим на график.
plt.figure(figsize=(10, 4))
norm_rv1 = stats.norm(loc=182, scale = 6.74)
norm_rv2 = stats.norm(loc=162, scale = 5.39)
x = np.linspace(140,210,100)
pdf1 = norm_rv1.pdf(x)
y = np.linspace(140,210,100)
pdf2 = norm_rv2.pdf(x)
plt.plot(x,pdf1)
plt.fill_between(y1=pdf1, x=x, alpha=0.4)
plt.plot(y,pdf2)
plt.fill_between(y1=pdf2, x=x, alpha=0.4)
l1=plt.axvline(174, color='g', ls='dashed', label='Наш вьетнамец')
l2=plt.axvline(194, color='b', ls='dashed', label='Наш серб')
plt.legend(handles=)
plt.show()
Кажется, вьетнамец победил, он находится ближе к правому краю колокола, чем серб. Как в этом убедиться?
Можно посчитать, на сколько стандартных отклонений рост серба и рост вьетнамца отличаются от средних значений. Это самое количество стандартных отклонений называется z-оценкой и вычисляется по такой формуле:
z=\frac{x-\mu}{\sigma}
То есть все, что надо сделать в случае с вьетнамцем — это вычесть 162 из 174 и поделить результат на 5,39.
Итак, рост вьетнамца отклоняется на 2,23 стандартного отклонения от среднего, а рост серба — на 1,78 стандартного отклонения от среднего, то есть вьетнамец у себя на родине воспринимается более высоким, чем серб у себя.
Таким образом можно стандартизировать любую нормально распределенную величину. По сути, мы наносим ее на график нормального распределения со средним, равным нулю и стандартным отклонением, равным единице.
x = np.linspace(-6,5,100)
pdf = norm_rv.pdf(x)
xticks =
plt.xticks(xticks)
plt.plot(x, pdf)
l1=plt.axvline(viet, color='g', ls='dashed', label='z-оценка\nроста вьетнамца')
l2=plt.axvline(serb, color='b', ls='dashed', label='z-оценка\nроста серба')
plt.legend(handles=)
plt.show()
Такое нормальное распределение со средним, равным 0, и среднеквадратическим отклонением, равным 1, называется стандартным нормальным распределением, или z-распределением.
У него есть одна любопытная и важная особенность. Зная z-оценку единицы наблюдения, мы можем также узнать, какой процент наблюдений превосходит значение этой единицы наблюдения. То есть мы можем прямо сейчас оценить, какая доля вьетнамцев выше или ниже ростом, чем наш вьетнамец. Для этого существует специальная таблица z-значений.
Z-оценка роста нашего вьетнамского приятеля составила 2,23 (с округлением до второго знака после запятой). В вертикальном столбце под буквой z находим 2,2, а в горизонтальном 0,03 (потому что 2,2 + 0,03 = 2,23), на пересечении стоит 0,9871. Это означает, что он выше 98,7% единиц наблюдения.
То есть если согнать всех-всех вьетнамцев на поляну и сделать из них гистограмму таким образом:
#International Statistical Review
#1975 vol: 43 (3) pp: 339
— только 2,2% людей справа будут выше него.
Для вычисления z-значения не обязательно использовать таблицу, можно получить более точную оценку с использованием специальных инструментов для работы с данными (это не обязательно должны быть R, SPSS или что-то подобное, можно поставить превосходное приложение "Probability Distributions" на телефон).
Зная про z-оценку, нетрудно догадаться, откуда взялось правило 68-95-99,7. Давайте посмотрим по таблице, какая доля единиц наблюдения лежит в интервале от -1 до 1.
Левее z=1 лежит 84,13% единиц наблюдения, левее z=-1 лежит 15,87% единиц наблюдения. Между ними, соответственно, 84,13-15,87 = 68,26%.
Можете проверить для 95 и 99,7, ну или попробовать для любых других точек.
#Этот график я уже строил выше, так что не буду копировать тонну кода и просто перетащу сам график.
fig_rule
Подведу итог. Z-оценка — это количество стандартных отклонений, на которое единица наблюдения отклоняется от среднего арифметического. С помощью таблицы z-оценок можно определить, какая доля единиц наблюдения лежит выше или ниже произвольного значения. Применима z-оценка только для нормально распределенных величин.
Центральная предельная теорема
Представьте: у меня есть две генеральные совокупности по миллиону объектов в каждой. В одной величина распределена нормально, в другой распределение скошенное.
#Делаю две совокупности по миллиону единиц наблюдения.
norm = norm_rv.rvs(size=1000000)
skewed = stats.skewnorm.rvs(5, size=1000000)
fig = plt.figure(figsize=(12, 5))
#Рисую их гистограммы вместе с графиками функций плотности вероятности.
ax = fig.add_subplot(121)
plt.hist(norm, normed=True, bins = 20)
x = np.linspace(-4,5,100)
pdf = norm_rv.pdf(x)
plt.axvline(norm.mean(), color='b', linestyle='dashed', linewidth=1)
plt.legend()
plt.plot(x, pdf)
ax = fig.add_subplot(122)
plt.hist(skewed, normed=True, bins = 20)
y = np.linspace(-1,4,100)
pdf = stats.skewnorm(5).pdf(y)
plt.axvline(skewed.mean(), color='b', linestyle='dashed', linewidth=1)
plt.legend()
plt.plot(y, pdf)
plt.show()
Что будет, если я извлеку много-много выборок из генеральной совокупности и посчитаю для каждой из этих выборок среднее арифметическое? Видимо, получится выборка из средних арифметических. Но как будет выглядеть гистограмма такой выборки?
Давайте проделаем этот трюк для каждой из наших генеральных совокупностей. Вытащим из каждой по десять тысяч выборок по 3 единицы наблюдения, посчитаем средние арифметические для всех них и получим две гистограммы, показывающие распределения средних арифметических.
#Вот так вытаскиваю выборки:
mean_list_norm = np.array()
mean_list_skewed = np.array()
fig = plt.figure(figsize=(12, 5))
fig.add_subplot(121)
plt.axvline(mean_list_norm.mean(), color='b', linestyle='dashed', linewidth=1)
plt.legend()
plt.hist(mean_list_norm, bins = 20)
fig.add_subplot(122)
plt.hist(mean_list_skewed, bins = 20)
plt.axvline(mean_list_skewed.mean(), color='b', linestyle='dashed', linewidth=1)
plt.legend()
plt.show()
Как видите, распределения выборочных средних похожи на распределения величин в генеральных совокупностях (хотя правое скошено несколько меньше). Средние арифметические примерно совпадают с истинными средними арифметическими, что неудивительно. Бросается в глаза, что эти гистограммы получились уже, чем гистограммы генеральных совокупностей (обратите внимание на ось x. Это тоже неудивительно, ведь, понятное дело, выборочные средние не так сильно варьируют, как единицы наблюдения во всей генеральной совокупности.
Давайте провернем этот трюк еще раз, но теперь будет извлекать выборки объемом 15.
mean_list_norm = np.array()
mean_list_skewed = np.array()
fig = plt.figure(figsize=(12, 5))
fig.add_subplot(121)
plt.axvline(mean_list_norm.mean(), color='b', linestyle='dashed', linewidth=1)
plt.legend()
plt.hist(mean_list_norm, bins = 20)
fig.add_subplot(122)
plt.hist(mean_list_skewed, bins = 20)
plt.axvline(mean_list_skewed.mean(), color='b', linestyle='dashed', linewidth=1)
plt.legend()
plt.show()
Гистограммы стали еще уже, а правая еще менее скошенной. Давайте еще раз. Извлечем из генеральных совокупностей по 10000 выборок, но теперь объемом 100.
mean_list_norm = np.array()
mean_list_skewed = np.array()
fig = plt.figure(figsize=(12, 5))
fig.add_subplot(121)
plt.axvline(mean_list_norm.mean(), color='b', linestyle='dashed', linewidth=1)
plt.legend()
plt.hist(mean_list_norm, bins = 20)
fig.add_subplot(122)
plt.hist(mean_list_skewed, bins = 20)
plt.axvline(mean_list_skewed.mean(), color='b', linestyle='dashed', linewidth=1)
plt.legend()
plt.show()
На вид две гистограммы нормально распределенной величины. Кстати, смотрите, они стали еще уже. Это абсолютно логично: чем больше размер выборки, тем более точную оценку среднего она может нам предоставить. Маленькие выборки хуже отражают параметры генеральной совокупности, чем большие.
И, как видите, чем больше выборка, тем сильнее распределение выборочных средних похоже на нормальное, и тем меньше будет рассеяние значений вокруг среднего (выше точность оценки среднего). Стало быть, чем больше объем выборки, тем меньше стандартное отклонение в выборке выборочных средних :)
Вот это стандартное отклонение в выборке, составленной из выборочных средних, называется стандартной ошибкой среднего. Приятная новость состоит в том, что для его определения не надо получать кучу выборок, есть простая маленькая формула:
SD_\overline{x} = \frac{\sigma}{\sqrt{n}} — истинная стандартная ошибка среднего.
Но, как вы понимаете, для оценки истинной стандартной ошибки среднего надо знать среднеквадратическое отклонение в генеральной совокупности, а это редкая роскошь, которой мы обычно не обладаем. Так что стандартную ошибку среднего мы оцениваем по выборочному стандартному отклонению:
\Huge SE_\overline{x} = \frac{s}{\sqrt{n}}
Некоторые исследователи используют стандартную ошибку среднего для описания выборки вместо стандартного отклонения (первая всегда меньше второй, поэтому иногда возникает такой соблазн). Этого делать не стоит, так как этот показатель нужен для оценки среднего в генеральной совокупности, он не является описательной статистикой.
Одной только стандартной ошибки для оценки среднего недостаточно. Но она позволит нам вычислять доверительные интервалы и проверять гипотезы, об этом мы сейчас поговорим.
Доверительные интервалы для среднего
Если выборка, для среднего арифметического которой мы строим ДИ, небольшая, — интервал будет более широким, а с увеличением ее размера будет все сужаться и сужаться.
Также надо помнить, что, например, 95%-й интервал для выборочного среднего будет уже, чем 99%-й, так как ему позволено чаще ошибаться и не включать истинное среднее. 99%-й более осторожен и дает более размытую оценку местонахождения истинного среднего.
В качестве примера построим доверительный интервал для среднего роста в датасете Гальтона.
Выборочное среднее составляет 66,76 дюймов, стандартная ошибка среднего — 0,12. Нас интересует 95%-й доверительный интервал, так что будем откладывать 1,96 стандартной ошибки среднего.
66,76 ± 1,96 * 0,12 = 66,76 ± 0,23
Таким образом, 95%-й ДИ от 66.52 до 67.0 дюймов.
Важный момент: поскольку доверительный интервал мы построили на основании z-распределения, распределение нашей выборки должно быть более-менее нормальным. И она не должна быть маленькой. Если распределение не особенно похоже на нормальное, а выборка маленькая, интервал следует строить на основе t-распределения, о котором мы поговорим немного позже. Знать, как устроены доверительные интервалы, конечно, полезно. Но вручную лучше их не считать. Во-первых это сопряжено с округлением, а во-вторых порождает больше простора для ошибок.
Проверка гипотез
Концепция
В разделе про z-оценку мы говорили о вьетнамце и сербе, оба они были весьма высокими. На самом деле там был был и третий человек, малайзиец. "Вообще-то я тоже серб," — сказал он. Однако это было довольно подозрительно, ведь, помимо всего прочего, рост малайзийца был 167 см.
Впрочем, хорошо, давайте на секунду предположим, что перед нами и правда серб. Какова вероятность встретить столь низкого (или еще более низкого) серба? На этот вопрос ответить нетрудно, нужно лишь узнать, насколько характерен такой рост для сербских мужчин. Если вероятность встретить такого или более низкого человека менее 5% (то есть если имеющееся значение отличается от среднего очень сильно), будем считать, что подозреваемый, скорее всего, пытается нас обмануть.
Иными словами, надо посмотреть, какая доля генеральной совокупности попадает в этот интервал:
plt.figure(figsize=(10, 4))
norm_rv = stats.norm(loc=182, scale = 6.74)
x = np.linspace(140,210,71)
y = norm_rv.pdf(x)
lol = plt.plot(x,y)
l1=plt.axvline(167, color='g', ls='dashed', label='Псевдосерб')
plt.fill_between(x], \
y], color='r', alpha=0.5, label='1.3% сербов')
plt.legend(handles=)
plt.show()
Сначала надо посчитать z-оценку (количество стандартных отклонений, на которое интересующее нас значение отклоняется от популяционного среднего). Как мы помним, средний рост сербского мужчины равен 182 см, стандартное отклонение 6,74 см.
Считаем z-оценку:
z=\frac{x-\mu}{\sigma} = \frac{167-182}{6,74} = -2.23
Итак, рост странного незнакомца отстоит от среднего роста серба на -2.23 стандартных отклонения. Посмотрим по таблице z-значений, какая доля сербов имеет такой же или еще меньший рост.
Итак, 1,3% сербов имеют такой или еще более низкий рост. Всего лишь 1,3% — не такое уж частое явление. Иными словами, у нас есть серьезный повод усомниться, что перед нами серб.
plt.figure(figsize=(10, 4))
x = np.linspace(-4,4,801, endpoint=True)
norm_rv = stats.norm(loc=0, scale = 1)
y1 = norm_rv.pdf(x)
plt.plot(x,y1)
l1=plt.axvline(-2.23, color='g', ls='dashed', label='Псевдосерб')
fill=plt.fill_between(x], \
y1], color='r', alpha=0.5, label='1.3% сербов')
plt.legend(handles=)
plt.show()
То, что мы сейчас сделали, называется проверкой гипотезы. Это невероятно важная, если не самая важная, концепция в доказательной медицине. Когда мы сравниваем уровень гликозилированного гемоглобина у диабетиков и здоровых людей и считаем уровень статистической значимости, мы проверяем гипотезу. Когда мы сравниваем доли умерших за 10 лет в группе больных гипертонической болезнью I стадии и в контрольной группе, мы тоже проверяем гипотезу. Мы все время этим занимаемся.
В ходе проверки гипотезы мы проходим несколько этапов:
- Формулируем нулевую гипотезу (H_0). Это обычно путь наименьшего сопротивления, самое унылое и серое утверждение: об отсутствии различия, отсутствии взаимосвязи, отсутствии изменений. В нашем случае мы предполагаем, что рост незнакомца относится к генеральной совокупности ростов сербов.
- Формулируем альтернативную гипотезу (H_1). Это уже более интересная гипотеза: есть различие, есть взаимосвязь, есть изменения, есть что-то любопытное. Ну или как в нашем случае: наш субъект слишком низкий, чтобы быть сербом.
- Решаем, что для нас такое "слишком низкий", "слишком большое различие" и т.п. То есть определяемся, в каком случае отвергнем нулевую гипотезу в пользу альтернативной. В нашем примере мы отвергаем H_0, если такие или еще более низкие сербы встречаются реже, чем в 5% случаев (то есть его слова неправдоподобны, и такой рост является редкостью). Выбранная нами вероятность называется уровнем α. В медицине, как правило, устанавливают пороговый уровень статистической значимости α = 0.05 (так мы ограничиваем вероятность ошибочно отвергнуть нулевую гипотезу; по исходному замыслу эта ошибка должна встречаться не чаще, чем в каждом 20 исследовании, — на практике, конечно, это не так, скоро это обсудим). Иными словами, величину считаем статистически значимой, если мала вероятность получить такое или еще более выдающееся ее значение (а что мы понимаем под "малой вероятностью" — как раз и есть α).
- Применяем статистический тест: рассчитываем некий показатель, который подчиняется своему закону распределения, в нашем случае это z-оценка, подчиняющаяся нормальному закону.
- Получаем уровень статистической значимости — p-значение. Если он меньше, чем выбранный нами уровень α (то есть вероятность получить такое значение меньше выбранного нами уровня), то мы отвергаем нулевую гипотезу. Если больше или равен, то мы не можем отвергнуть нулевую гипотезу. Потому что вы видите, что полученные вами данные плохо с ней уживаются.
- Низкое p-значение: ваши данные плохо уживаются с нулевой гипотезой.
- Высокое p-значение: ваши данные хорошо уживаются с нулевой гипотезой.
- Оно не показывает вероятность того, что имеющиеся результаты получены случайно.
- Оно не показывает вероятность того, что нулевая гипотеза верна. Или что альтернативная неверна.
- Если p-значение слишком велико, чтобы отвергнуть нулевую гипотезу, это не значит, что нулевая гипотеза верна.
- Маленькое p-значение не говорит о высокой эффективности лечения.
- Маленькое p-значение не означает, что наши результаты воспроизводимы.
- Маленькое p-значение не означает, что исследование автоматически имеет научную ценность.
Одновыборочный z-тест
А центральная предельная теорема говорит, что для генеральной совокупности распределение ее выборочных средних будет иметь такое же среднее, а стандартное отклонение будет меньше в \sqrt{n} раз (и будет называться стандартной ошибкой среднего, \frac{\sigma}{\sqrt{n}}), ну и самое главное — с увеличением размера выборок распределение выборочных средних будет все ближе к нормальному. Вот картинка для простоты:
scale = 1
norm_rv1 = stats.norm(loc=0, scale = 1)
norm_rv2 = stats.norm(loc=0, scale = scale/(30**0.5))
sample_norm1 = norm_rv1.rvs(size=10000)
sample_norm2 = norm_rv2.rvs(size=5000)
fig = plt.figure(figsize=(16, 7))
p = plt.hist(sample_norm1, 15, alpha=0.8, color='#A9E2F3')
sd = plt.hist(sample_norm2, 15, alpha=0.8, color='#A9A9F5')
#Красивые прямоугольнички для легенды:
import matplotlib.patches as mpatches
p = mpatches.Patch(color='#A9E2F3', label='Генеральная совокупность')
sd = mpatches.Patch(color='#A9A9F5', label='Выборочные средние, n=30')
plt.legend(handles=(p, sd))
plt.show()
Теперь давайте вообразим, что мы проводим опыт. Предположим, для эксперимента мы используем крыс из определенной линии, для которой известна средняя масса тела — 400 грамм. Стандартное отклонение (посчитано по огромной выборке и хорошо аппроксимирует генеральное среднеквадратическое отклонение) составляет 18.
Мы взяли 20 крысят и с раннего возраста накачивали их соматотропным гормоном. В итоге получили 20 взрослых крыс со средней массой тела 410.
Много это или мало? Можно ли сказать, что инъекция гормона статистически значимо повысила массу тела крыс?
Давайте проверять гипотезу.
H_0: \mu=\mu_0 — мы предполагаем, что средняя масса тела в генеральной совокупности крыс, которые получили инъекцию гормона, не отличается от средней массы тела крыс, которые не участвовали ни в каких экспериментах.
H_1: \mu>\mu_0 — альтернативная гипотеза: средняя масса тела крыс, которым делают инъекцию соматотропина, больше средней массы тела обычных крыс.
Различия будем считать достоверными при уровне значимости меньше 0,05.
Как вы помните, z-оценку для одного объекта можно получить по такой формуле: z=\frac{x-\mu}{\sigma}. Мы рассматриваем этот объект как часть генеральной совокупности и оцениваем, где он находится в ее распределении. Но сейчас мы хотим оценить выборочное среднее. Понять, насколько далеко оно отстоит от среднего в распределении выборочных средних (:D). А значит, нам нужны сведения о распределении выборочных средних (для выборок такого объема, как наша, то есть 20 единиц наблюдения). Его параметры мы отлично себе представляем: среднее как в генеральной совокупности, 400 грамм, а стандартное отклонение равно \frac{18}{\sqrt{20}}(оно же стандартная ошибка среднего). Вот теперь и посмотрим, как далеко отклоняется среднее по нашей выборке от среднего в распределении выборочных средних (да, абзац вышел довольно жестким).
Формула, очевидно, будет почти такая же, только в знаменателе вместо σ появится стандартная ошибка среднего, а вместо x будет \overline{x}.
z=\frac{\overline{x}-\mu_0}{\sigma/\sqrt{n}} = \frac{410-400}{18/\sqrt{20}} = 2,48
Смотрим в таблицу z-оценок: наше значение больше, чем 99,34% популяции.
plt.figure(figsize=(10, 4))
norm_rv = stats.norm(loc=400, scale = 18/(20**0.5))
x = np.linspace(380,420,401)
y = norm_rv.pdf(x, )
plt.plot(x,y)
l1=plt.axvline(410, color='g', ls='dashed', label='Наше выборочное среднее')
fill=plt.fill_between(x], \
y], color='r', alpha=0.5, label='99,34% выборочных средних')
plt.legend(handles=)
plt.show()
Но тут нам нужно посмотреть не это, нас интересует правый хвост распределения — какая доля выборочных средних отклоняются от генерального среднего так же или еще более значительно. Это сделать легко: если слева 99,34%, то справа 100-99,34%, или 0.66%, то есть p=0.0066, что явно меньше установленного нами уровня α.
plt.figure(figsize=(10, 4))
norm_rv = stats.norm(loc=400, scale = 18/(20**0.5))
x = np.linspace(380,420,401)
y = norm_rv.pdf(x, )
plt.plot(x,y)
l1=plt.axvline(410, color='g', ls='dashed', label='Наше выборочное среднее')
fill=plt.fill_between(x:], \
y:], color='r', alpha=0.5, label='0,066% выборочных средних')
plt.legend(handles=)
plt.show()
Итак, вероятность получить такое же или еще более экстремальное значение при условии справедливости нулевой гипотезы — 0,0066, это очень мало. То есть наблюдаемые нами данные сильно противоречат нулевой гипотезе, и мы ее отвергаем.
Тут надо заметить, что альтернативная гипотеза может быть не только в стиле "больше" или "меньше" — это, скорее, редкость. Обычно мы не знаем заранее, как повлияет вносимое нами различие, так что нужно использовать не одностороннюю, а двустороннюю альтернативу:
H_1: \mu\neq\mu_0 — то есть смотреть в таком случае следует на оба хвоста. В нашем случае z-оценка составила 2,48, но при такой альтернативной гипотезе нужно учесть и хвост левее -2.48. Ведь там единицы наблюдения будут отклоняться столь же сильно или сильнее, просто в другую сторону.
plt.figure(figsize=(10, 4))
norm_rv = stats.norm(loc=400, scale = 18/(20**0.5))
x = np.linspace(380,420,401)
y = norm_rv.pdf(x, )
plt.plot(x,y)
fill=plt.fill_between(x:], \
y:], color='r', alpha=0.5, label='(0,066 * 2)% выборочных средних')
plt.fill_between(x], \
y], color='r', alpha=0.5)
plt.legend(handles=)
plt.show()
Хорошая новость состоит в том, что, раз уж распределение симметрично, лазить по таблице не нужно: меньше z=-2.48 столько же единиц наблюдения, сколько и больше z=2.48, то есть всего 0.0066 * 2 = 0,0132.
Стало быть, двусторонний вариант альтернативной гипотезы мы также можем принять; средняя масса тела крыс, получавших инъекции соматотропина, статистически значимо отличается от популяционного среднего, p<0,05.
Вкратце
Одновыборочный z-критерий используется, чтобы оценить значимость различий между предполагаемым популяционным средним \mu_0 и выборочным средним \overline{x}.
Ограничения к применению критерия:
- Мы должны знать генеральную дисперсию (или среднеквадратическое отклонение) или располагать ее оценкой с высокой точностью. Если этого нет, то лучше использовать t-статистику, о ней скоро поговорим. Это нужно потому, что в формуле фигурирует стандартная ошибка среднего \sigma/\sqrt{n}
- Тест подходит для работы с нормально распределенными величинами. И в вашей работе необходимо указать, что это было проверено (как проверить, поговорим в разделе про двухвыборочный t-тест). Если распределение немного отличается от нормального, но выборка очень большая, центральная предельная теорема все равно дает нам возможность использовать тест (ведь нас интересует распределение выборочных средних). Но будьте осторожны, стройте вначале гистограмму, бывают очень сильно скошенные распределения.
Как вы понимаете, в таких случаях z-статистика перестает работать. Надо заменить распределение, на которое мы опирались до этого, на другое, рассчитанное для размера той выборки, которая у нас есть.
В игру вступает t-статистика, основанная, как следует из названия, не на z-распределении, а на t-распределении Стьюдента. "Стьюдент" — это псевдоним английского ученого Уильяма Сили Госсета. Он работал на всем известную пивоваренную компанию Гинесс. Из страха утечек конфиденциальной информации компания запретила своим сотрудникам публикацию каких бы то ни было исследовательских материалов. Однако Госсет опубликовал все равно, но из осторожности выступил под псевдонимом.
Итак, t-распределение — это на самом деле целое семейство распределений, напоминающих нормальное, но с более тяжелыми хвостами. У выборки любого размера есть свое распределение Стьюдента. Соответственно, если нормальное распределение описывалось двумя параметрами — средним значением и стандартным отклонением, — у t-распределений параметров три: среднее, стандартное отклонение и число степеней свободы.
Число степеней свободы - это довольно сложная концепция. Если не влезать в дебри, суть ее в том, что чем меньше объем выборки, тем меньше степеней свободы. А чем меньше степеней свободы, тем более жирные хвосты у распределения. Чаще всего число степеней свободы приравнивают к n - 1, то есть на один меньше объема выборки. О том, когда стоит так делать, и стоит ли вообще, поговорим немного позже.
from matplotlib import cm
cm = plt.get_cmap('flag')
plt.figure(figsize=(12, 14))
ax1 = plt.subplot(211)
ax1.set_color_cycle()
x = np.linspace(-5,5,1000)
y = np.linspace(-5,5,1000)
labels =
for i in range(2, 30, 3):
rv = stats.t(loc=0, scale = 1, df=i)
pdf = rv.pdf(x)
labels.append(plt.plot(x,pdf,linewidth=1, label=('Распределение Стьюдента, df= ' + str(i))))
rv_n = stats.norm(loc=0, scale = 1)
pdf_n = rv.pdf(x)
labels.append(plt.plot(x,pdf_n, label=('Нормальное распределение')))
plt.legend(handles= for x in labels])
plt.title('t-распределения с разным числом степеней свободы')
ax2 = plt.subplot(212)
ax2.set_color_cycle()
x = np.linspace(-5,5,1000)
y = np.linspace(-5,5,1000)
labels =
for i in range(2, 30, 3):
rv = stats.t(loc=0, scale = 1, df=i)
pdf = rv.pdf(x)
labels.append(plt.plot(x,pdf, label=('Распределение Стьюдента, df= ' + str(i))))
rv_n = stats.norm(loc=0, scale = 1)
pdf_n = rv.pdf(x)
labels.append(plt.plot(x,pdf_n, label=('Нормальное распределение')))
plt.legend(handles= for x in labels])
plt.title('Посмотрим поближе на верхушки распределений')
plt.ylim(0.33, 0.4)
plt.xlim(-0.25, 0.8)
plt.show()
Как видите, с ростом числа степеней свободы t-распределение все ближе к нормальному. Поэтому на огромных выборках нет принципиальной разницы, z- или t-статистику использовать. "Что такое огромная выборка? Как понять, когда еще нельзя использовать z, а когда уже можно?" — зададите вы резонный вопрос.
На самом деле проще всегда использовать t-статистику вместо z-статистики.
Таким образом, из-за того, что t-распределение для маленьких выборок обладает жирными хвостами, нам сложнее получить волшебное p < 0,05. Мы получаем большее p-значение, чем если бы использовали z-критерии. Это очень хорошо и разумно. За небольшой объем выборки надо расплачиваться.
plt.figure(figsize=(12, 6))
norm_rv = stats.norm(loc=0, scale = 1)
t_rv = stats.t(loc=0, scale = 1, df=3)
x = np.linspace(-5,5,1001)
pdf1 = norm_rv.pdf(x)
pdf2 = t_rv.pdf(x)
p1, = plt.plot(x,pdf1, c='b', alpha=0.4, label='Нормальное распределение')
p2, = plt.plot(x,pdf2, color='r', label='Распределение Стьдента (df=3)\nс более тяжелыми хвостами')
fill1=plt.fill_between(x:], \
pdf2:], alpha=0.5, color='r')
fill2=plt.fill_between(x:], \
pdf1:], alpha=0.5, color='b')
fill3=plt.fill_between(x], \
pdf2], alpha=0.5, color='r')
fill4=plt.fill_between(x], \
pdf1], alpha=0.5, color='b')
plt.legend(handles=)
plt.show()
На графике хорошо видно, что если величина выборки всего 4 элемента, вам будет намного сложнее попасть в закрашенную зону, где нулевая гипотеза отвергается.
Одновыборочный t-тест
Одновыборочный t-критерий используется, чтобы оценить значимость различий между предполагаемым популяционным средним \mu_0 и выборочным средним \overline{x}. Такая задача перед исследователем стоит редко, но этот критерий удобно рассматривать в качестве перехода от одновыборочного z-теста к двухвыборочному t-тесту.
Ограничения к применению критерия:
- В отличие от одновыборочного z-теста, тут генеральную дисперсию знать не требуется. Если так вышло, что вы ее откуда-то знаете или располагаете ее оценкой с высокой точностью, это тот редкий случай, когда z-статистика подходит лучше.
- Тест годится для работы с нормально распределенными величинами. Если распределение немного отличается от нормального, но выборка очень большая, центральная предельная теорема все равно дает нам возможность использовать тест (ведь нас интересует распределение выборочных средних). Но будьте осторожны, стройте вначале гистограмму, бывают и сильно скошенные распределения.
Нас традиционно будет интересовать уровень значимости 0,05. Смотрим на пересечение столбца 0,05 и ряда 24, там стоит 2,064.
plt.figure(figsize=(12, 6))
t_rv = stats.t(loc=0, scale = 1, df=24)
x = np.linspace(-5,5,1001)
pdf = t_rv.pdf(x)
p, = plt.plot(x,pdf, color='r', label='Распределение Стьдента (df=3)\nс более тяжелыми хвостами')
fill1=plt.fill_between(x:], \
pdf:], alpha=0.5, color='r')
fill3=plt.fill_between(x], \
pdf], alpha=0.5, color='r')
plt.legend(handles=)
plt.xticks()
plt.show()
Иными словами, t-значения, которые будут отстоять от среднего еще дальше, чем 2,064, будут указывать на статистически значимое отличие на уровне α=0,05.
Наше t-значение составляет 2.5, то есть различие статистически значимо.
Двухвыборочный t-тест
Наконец мы подобрались к тому, ради чего многие, возможно, и читают эту статью. Двухвыборочный тест Стьюдента применяется в медицинских исследованиях очень часто, но в большинстве случаев неправильно.
Здесь мы рассмотрим модификацию Уэлча для двухвыборочного t-теста, а потом я объясню, почему всегда стоит использовать именно ее.
Итак, тест нужен для того, чтобы оценить различие между средними значениями двух выборок.
У метода, как водится, есть ограничения.
- Сравниваемые выборки должны быть извлечены из совокупностей с нормальным распределением интересующей нас величины. Всегда нужно проверять выборку на нормальность и сообщать об этом в своей работе (как это проверить, скоро обсудим; просто посмотреть на гистограммы недостаточно). Однако если распределение немного скошено, а выборки крупные (например, по 100 единиц наблюдения), центральная предельная теорема все равно позволяет использовать этот критерий. Но будьте осторожны. Если сомневаетесь, воспользуйтесь критериями, которые не требуют нормальности (непараметрические критерии, например U-тест Манна—Уитни).
- Стандартный вариант двухвыборочного t-теста требует равенства генеральных дисперсий для сравниваемых групп. Если дисперсии различаются значимо, критерий использовать нельзя. Поэтому следует непременно проверить это условие и упомянуть об этом в работе. Для проверки можно, например, использовать тест Левена. Нужно все это из-за одной досадной статистической проблемы (проблемы Беренса—Фишера): невозможно оценить разницу между средними арифметическими совокупностей, если у них разная дисперсия.
- Хорошая новость состоит в том, что у проблемы Беренса—Фишера есть приближенные решения, одним из которых является t-тест Уэлча, который мы и собираемся рассмотреть. Поэтому нам можно не переживать о неравенстве генеральных дисперсий. Не забудьте, в ваших работах обязательно нужно написать о проверке на равенство генеральных дисперсий, или же сообщить, что вы используете модифицированный t-критерий.
Достаточно ли нормальны распределения? По гистограмме сказать невозможно. К счастью, есть специальные методы проверки нормальности. Их довольно много, наиболее часто применяются критерии Шапиро—Уилка и Колмогорова—Смирнова. Да, это отличные критерии, и вы можете их применять. Однако у них есть недостаток. На маленьких выборках они, как правило, не отвергают нормальность. А при проверке больших выборок вы, скорее всего, отклоните нормальность, даже если распределение отличается от нормального совсем чуть-чуть. Это особенно неуместно с учетом того, что с ростом выборки у t-критерия снижается требовательность к нормальности распределения (помните центральную предельную теорему?).
Но существует и более подходящий нашей задаче способ проверки нормальности. Это Q-Q-график. Выглядит он как квадрат с проведенной по диагонали линией, вдоль которой выстраиваются единицы наблюдения интересующей выборки. Не будем вдаваться в подробности его построения, а сразу посмотрим на пример:
plt.figure(figsize=(12,10))
sample1 = stats.skewnorm.rvs(0, size=100)
sample2 = stats.skewnorm.rvs(10, size=100)
#Ура, есть готовая реализация функции построения Q-Q-плота.
plt.subplot(2,2,1)
plt.hist(sample1, alpha=0.5)
plt.subplot(2,2,2)
plt.hist(sample2, alpha=0.5)
plt.subplot(2,2,3)
stats.probplot(sample1, dist="norm", plot=plt)
plt.subplot(2,2,4)
stats.probplot(sample2, dist="norm", plot=plt)
plt.show()
Если распределение нормальное, единицы наблюдения выстраиваются более-менее по диагональной линии, а если скошенное, то отклоняются. «Опять "более-менее",» — наверное, возмутитесь вы. Ну да. Это значит, что четкой границы между "еще можно" и "уже нельзя" нет. Решать вам. Это еще один пример того, что статистика может быть весьма субъективной.
Посмотрим на Q-Q-графики для мужчин и женщин из выборки Гальтона.
plt.figure(figsize=(12,5))
plt.subplot(1,2,1)
stats.probplot(galton.Height, dist="norm", plot=plt)
plt.title('Мужчины')
plt.subplot(1,2,2)
stats.probplot(galton.Height, dist="norm", plot=plt)
plt.title('Женщины')
plt.show()
Достаточно ли нормальны эти распределения для использования критерия Стьюдента? Да, безусловно. Редко такие увидишь. Да и в любом случае объем первой выборки 465, а второй 433. Они довольно крупные, и даже будь распределения немного скошены, мы бы все равно могли безопасно использовать t-статистику.
При этом критерий Шапиро—Уилка отвергает гипотезу о нормальности обеих выборок.
#sha = stats.shapiro(galton.Height)
#print('Критерий Шапиро - Уилка отвергает гипотезу о нормальности распределения на уровне значимости p =', round(sha, 2))
Итак.
H_0: \mu_1 = \mu_2
H_1: \mu_1 \neq \mu_2
Средний рост мужчины по нашей выборке составляет 69,23 дюйма (стандартное отклонение 2,63), женщины — 64,11 дюйма (стандартное отклонение 2,37). Мужчин 465 человек, женщин 433. Подставляем значения в формулу:
t = \frac{\overline{X}_1 - \overline{X}_2}{\sqrt{\frac{s^2_1}{n_1} + \frac{s^2_2}{n_2}}} = \frac{69,23 - 64,11}{\sqrt{\frac{2,63^2}{465} + \frac{2,37^2}{433}}} = 30,68
Теперь из семейства t-распределений надо выбрать то, которое нам подходит. Считаем число степеней свободы:
\nu \approx \frac{\left(\frac{s^2_1}{n_1} + \frac{s^2_2}{n_2}\right)^2}{\frac{s^4_1}{n^2_1(n_1 - 1)} + \frac{s^4_2}{n^2_2(n_2-1)}} =
\frac{\left(\frac{2,63^2}{465} + \frac{2,37^2}{433}\right)^2}{\frac{2,63^4}{465^2(465 - 1)} + \frac{2,37^4}{433^2(433-1)}} = 895
Что ж, славно, нам нужно распределение Стьюдента с 895 степенями свободы. В упрощенных таблицах, конечно, такого значения нет. Но есть удобные калькуляторы (например, приложение Probability Distributions, которое можно поставить из Эпп-Стора или Плэй-Маркета), которые говорят, что для такого распределения Стьюдента t-значение, равное 30,68, позволяет отвергнуть нулевую гипотезу с p < 0,001.
И еще раз агитирую вас использовать для вычислений готовые программные решения. В Экселе можно посчитать все это за две минуты, просто указав выборки, нужную функцию и тип альтернативы. Кроме того, вы избежите округлений в расчетах и получите более точный ответ (например, мы бы получили t = 30,66, а не t = 30,68; впрочем, для нас это не критично, потому что p-значение крайне мало).
#print(len(galton.Height))
#print(len(galton.Height))
#print(round(galton.Height.mean(), 2))
#print(round(galton.Height.mean(), 2))
#print(round(galton.Height.std(), 2))
#print(round(galton.Height.std(), 2))
# Проверим результат:
#stats.ttest_ind(galton.Height,\
# galton.Height,\
# equal_var = False)
Проверка гипотез и доверительные интервалы
Гипотезы часто проверяют с помощью доверительных интервалов. Например, вместо одновыборочного z-теста можно просто построить доверительный интервал для выборочного среднего и посмотреть, включает ли он генеральное среднее. Допустим, генеральное среднее равно нулю. А выборочное среднее равно 0,6 (95%-й ДИ 0,2-1,0). Как видите, доверительный интервал не включает 0, поэтому нулевую гипотезу можно отвергнуть.
Или, к примеру, взять относительный риск. Этот показатель используется для оценки риска развития заболевания у группы, подвергаемой некоторому воздействию, по сравнению с группой без этого воздействия (например, группа курильщиков и группа некурящих). Если относительный риск равен единице, это означает, что разницы между группами нет. Если относительный риск меньше единицы — в экспериментальной группе заболевание развивается реже, чем в контрольной. А если относительный риск больше единицы — наоборот, в контрольной группе заболевание развивается реже, чем в экспериментальной.
И если вы видите в научной работе, что относительный риск составил 0,9 (95%-й ДИ 0,8 - 1,0), вы понимаете, что на уровне значимости 0,05 нельзя отвергнуть нулевую гипотезу об отсутствии различий между группами. Потому что в интервал попала единица — значение относительного риска, при котором различий нет.
Также с помощью доверительных интервалов иногда сравнивают средние значения двух выборок. Сделать это можно двумя способами. Например, можно построить доверительные интервалы для средних каждой выборки и посмотреть, пересекаются они или нет. Проблема тут лишь в том, что этот способ не слишком чувствителен к различию средних. Например, вы убедились, что два 95%-х доверительных интервала не пересекаются. Но на самом деле этот факт обеспечивает статистическую значимость различий не 0,05, а в районе 0,005.
То есть если интервалы для средних не пересекаются, вы можете уверенно отвергать нулевую гипотезу с α=0,05. А если пересекаются, то это не значит, что различий нет, потому что на самом деле ошибка первого рода ограничена куда более жестким порогом.
Разумнее использовать доверительный интервал для разницы средних. Его формула сложнее, но это значительно более мощный метод.
Не забудьте, что доверительные интервалы наследуют ограничения, свойственные соответствующим критериям. Поэтому если вы, к примеру, хотите использовать t-интервал для разности средних, вам нужно сначала убедиться в равенстве генеральных дисперсий и нормальности распределений.
Проблема множественной проверки гипотез
Есть две игры. Правила их таковы:
- Игра №1. Я подкидываю монетку. Если выпадает орел, вы получаете тысячу рублей.
- Игра №2. Я подкидываю монетку 10 раз. Если хоть раз из этих десяти выпадет орел, вы получите тысячу рублей.
Наиболее часто используемая поправка на множественные сравнения — поправка Бонферрони.
Устроена она предельно просто. Все полученные p-значения домножаются на количество проверяемых гипотез, и потом уже сравниваются с α.
Иными словами, чтобы отвергнуть нулевую гипотезу на уровне значимости 0,05 при десяти проверяемых гипотезах, p должно быть меньше 0.005.
Поправка Бонферрони — чрезвычайно жесткий метод. При увеличении числа гипотез резко падает шансы отклонить хоть какие-то из них резко падают.
Во всех отношениях лучше метод Холма. Его стоит использовать вместо поправки Бонферрони абсолютно всегда. Он дает не слишком большой прирост мощности, но все же действует уже не столь радикально.
Работает он так.
Вначале p-значения выстраиваются по возрастанию. Давайте для примера возьмем 4 p-значения:
0,01; 0,02; 0,03; 0,04
Затем каждому значению мы присваиваем коэффициент, равный m + 1 - k, где m — длина нашего вариационного ряда, а k — номер конкретного p-значения в ряду. Выглядит это понятнее, чем звучит:
4; 3; 2; 1
Умножаем p-значения на коэффициенты:
0,04; 0,06; 0,06; 0,04
А теперь идем по списку и смотрим, чтобы каждое следующее значение было не меньше предыдущего. А если вдруг оно меньше, то заменяем его на предыдущее.
0,04; 0,06; 0,06; 0,06
Как видите, на уровне значимости 0,05 мы можем отвергнуть только первую нулевую гипотезу.
#Проверю на всякий.
#from statsmodels.sandbox.stats.multicomp import multipletests
#pvals =
#reject, p_corrected, a1, a2 = multipletests(pvals,
# alpha = 0.05,
# method = 'holm')
#p_corrected
Такое ограничение вероятности совершения по крайней мере одной ошибки первого рода называется ограничением FWER (family-wise error rate). И поправка Бонферрони, и метод Холма ограничивают FWER.
Но если проверяется большое количество гипотез, как, например, в биоинформатике, когда сравнивают геномы людей, или же в нейробиологии при сравнении данных фМРТ, — часто используется ограничение FDR (false discovery rate). При этом ограничивается не вероятность по крайней мере одного ложноположительного срабатывания, а средняя доля ложноположительных срабатываний. Например, мы можем жить с мыслью, что в наших выводах 5% ложноположительных заключений. Что ж, хорошо, тогда применим метод, ограничивающий FDR, например метод Бенджамини—Хохберга.
Необходимо всегда явно сообщать в своей работе об использовании коррекции на множественное сравнение и указывать метод, который вы применили.
p-хакинг
Отчего-то исследования, в которых не удалось совершить открытие, многие авторы считают неудавшимися и не публикуют. Для науки отрицательный результат не менее ценен, но он не так здорово смотрится в публикации. Поэтому часто исследователи делают ужасную вещь.
Они берут свои данные и начинают искать что-нибудь интересное. Какие-нибудь шокирующие различия или закономерности. А когда найдут, проверяют их значимость прямо на этих же данных. Если p меньше 0,05 — прекрасно, научная работа готова. Если не меньше, — что ж, ищем дальше! Главное, потом не указывать в работе, сколько гипотез было проверено.
Собственно, это и называется p-хакингом. Манера мучить данные, пока они не разродятся чем-то подходящим.
На самом деле проверять все найденные интересные закономерности нужно на новых данных. Проблема, конечно, в том, что их еще нужно собирать.
Надо понимать, что к моменту, когда вы начинаете анализ собранных данных, вы должны уже определиться с гипотезами, которые намереваетесь проверить. Иначе ошибка множественных сравнений становится частью вашего исследования.
Столь же безответственно перебирать критерии один за другим, пока наконец какой-то не выдаст приемлемое p. Или выбирать тип альтернативы уже после того, как вы увидели данные.
Существует огромный простор для манипуляции результатами, и исследователи часто пользуются этим безо всякого стеснения. Одна моя умная знакомая как-то раз пошутила, что статистика — это набор грязных хаков. Возможно, в области проверки гипотез сегодня так и есть.
Что еще почитать и посмотреть
В статье я затронул только несколько базовых моментов. Но они должны помочь вам двигаться дальше. Чтобы понять некоторые вещи (p-значение, как мне кажется, одна из них), их нужно прочитать или услышать много раз. Поэтому посоветую вам несколько источников:
Основы статистики — прекрасный курс Анатолия Карпова на Stepic. Ссылка ведет к первой его части, всего их три.
"Медико-биологическая статистика" Стентона Гланца — всем известная книга, горячо рекомендую.
"Statistics for the Behavioral Sciences", Frederick J Gravetter, Larry B. Wallnau — еще одна хорошая книга. На русский, к сожалению, не переведена.
"100 Statistical Tests", Gopal K. Kanji — небольшой справочник, в котором собраны 100 основных статистических критериев.
Statistical Inference — красивые интерактивные графики нескольких основных концепций статистики.
Misunderstanding of p-values — статья в Википедии про заблуждения относительно p-значений.
Голова профессора Бамблдорфа — вводящая в ступор заметка о неверной интерпретации доверительных интервалов.
Автор: Михаил Берегов
Редакция: Елизавета Зайцева
Обложка: Никита Родионов
Нашли опечатку? Выделите фрагмент и нажмите Ctrl+Enter.
