Глава 2. Распределение данных

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

Таблицы частот

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

Для представления распределений мы будем использовать библиотеку empiri­caldist. В данном контексте «эмпирический» означает, что распределения основаны на данных, а не на математических моделях. Библиотека empiricaldist предоставляет класс FreqTab, который можно использовать для вычисления и графического отображения таблиц частот. Импортируем его следующим образом:

from empiricaldist import FreqTab

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

t = [1.0, 2.0, 2.0, 3.0, 5.0]

Класс FreqTab предоставляет метод from_seq, который принимает на вход последовательность и создает объект класса FreqTab:

ftab = FreqTab.from_seq(t)

ftab

	частота
1.0	1
2.0	2
3.0	1
5.0	1

Объект FreqTab — это разновидность Series из Pandas, в котором содержатся значения и их частоты. В этом примере значение 1.0 соответствует частоте 1, значение 2.0 — частоте 2 и т.д.

FreqTab также предоставляет метод bar, отображающий таблицу частот в виде столбчатой диаграммы:

ftab.bar()

decorate(xlabel="Значение", ylabel="Частота")

Поскольку класс FreqTab является производным от Series из Pandas, можно воспользоваться оператором квадратных скобок для поиска значения и получения его частоты:

ftab[2.0]

2

Но, в отличие от Series, для поиска значений объект FreqTab также можно вызывать как функцию:

ftab(2.0)

2

Если попытаться найти значение, которого нет в объекте FreqTab, то вызов функции вернет 0:

ftab(4.0)

0

Объект FreqTab имеет атрибут qs, который содержит массив значений. qs здесь означает «величины» (quantities), хотя, строго говоря, не все значения являются количественными величинами:

ftab.qs

array([1., 2., 3., 5.])

Во FreqTab также есть атрибут fs, который содержит массив частот (frequencies):

ftab.fs

array([1, 2, 1, 1])

Класс FreqTab предоставляет метод items, который можно использовать для циклического прохождения пар «величина — частота»:

for x, freq in ftab.items():

    print(x, freq)

1.0 1

2.0 2

3.0 1

5.0 1

В дальнейшем вы познакомитесь и с другими методами класса FreqTab.

Распределения датасета NSFG

Когда вы приступаете к работе с новым датасетом, я советую исследовать переменные, которые вы планируете использовать по одной за раз; и хорошая практика — начать с просмотра таблиц частот.

В качестве примера рассмотрим данные из датасета NSFG. В предыдущей главе мы загрузили его, считали в объект Pandas DataFrame и очистили некоторые переменные. Код, который мы использовали для загрузки и очистки данных, находится в модуле под названием nsfg.py. Инструкции по его установке приведены в Jupyter-блокноте к этой главе.

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

from nsfg import read_fem_preg

 

preg = read_fem_preg()

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

Чтобы выбрать строки, в которых переменная outcome равна 1, можно использовать метод query:

live = preg.query("outcome == 1")

В строке, передаваемой в query, имена переменных, такие как outcome, ссылаются на имена столбцов в объекте DataFrame. Она также может содержать операторы, подобные ==, и операнды, такие как 1.

Теперь, чтобы подсчитать, сколько раз каждая величина появляется в пере­менной birthwgt_lb, являющейся частью веса при рождении в фунтах, можно воспользоваться методом FreqTab.from_seq. Аргумент name присваивает объекту FreqTab имя, которое используется в качестве метки при построении графиков:

ftab_lb = FreqTab.from_seq(live["birthwgt_lb"], name="birthwgt_lb")

Вот как выглядит это распределение:

ftab_lb.bar()

decorate(xlabel="Фунты", ylabel="Частота")

Глядя на него, первое, что бросается в глаза, — это форма, которая напоминает знаменитую колоколообразную кривую, формально называемую нормальным распределением или распределением Гаусса. Другой примечательной особенностью этого распределения является мода (mode) — значение, которое встречается чаще всего. Чтобы найти моду, можно использовать метод idxmax, который находит величину, соответствующую самой высокой частоте:

ftab_lb.idxmax()

7.0

Класс FreqTab имеет метод mode, который дает тот же результат:

ftab_lb.mode()

7.0

В данном распределении мода составляет 7 фунтов.

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

ftab_oz = FreqTab.from_seq(live["birthwgt_oz"], name="birthwgt_oz")

ftab_oz.bar()

decorate(xlabel="Унции", ylabel="Частота")

Поскольку природа не осведомлена о фунтах и унциях, то можно было ожидать, что все значения birthwgt_oz будут равновероятными — то есть полученное распределение должно быть равномерным (uniform). Но оказывается, что 0 встречается чаще других величин, а 1 и 15 — реже; это говорит о том, что респонденты регулярно округляют те значения веса при рождении, которые близки целому числу фунтов.

Еще один пример: рассмотрим таблицу частот переменной agepreg, содержащей возраст матери в конце беременности:

ftab_age = FreqTab.from_seq(live["agepreg"], name="agepreg")

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

ftab_age.bar(width=0.1)

decorate(xlabel="Возраст", ylabel="Частота")

Распределение только отдаленно напоминает колокол и при этом скошено вправо, то есть его хвост тянется вправо дальше, чем влево.

В заключение взглянем на таблицу частот переменной prglngth, представляющей собой продолжительность беременности в неделях. Аргумент xlim устанавливает пределы по оси x в диапазоне от 20 до 50 недель: значений за пределами этого диапазона не так много, и, вероятно, это ошибки.

ftab_length = FreqTab.from_seq(live["prglngth"], name="prglngth")

ftab_length.bar()

decorate(xlabel="Недели", ylabel="Частота", xlim=[20, 50])

Срок 39 недель встречается заметно чаще других. Также левый хвост распределения длиннее правого: ранние роды случаются регулярно, но беременность редко длится дольше 43 недель, и в таких случаях врачи часто принимают меры.

Выбросы

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

Для выявления выбросов функция ниже принимает объект FreqTab и целое число n, чтобы с использованием индекса среза выбрать n наименьших (smallest) величин и их частот:

def smallest(ftab, n=10):

    return ftab[:n]

Вот 10 наименьших значений в таблице частот переменной prglngth:

smallest(ftab_length)

prglngth

0     1

4     1

9     1

13    1

17    2

18    1

19    1

20    1

21    2

22    7

Name: prglngth, dtype: int64

Поскольку мы отобрали случаи рождения живых детей, то продолжительность беременности менее 10 недель, безусловно, является ошибочной. Скорее всего, причина в том, что результат был указан неверно. Продолжительность беременности более 30 недель, судя по всему, является корректной. О значениях между 10 и 30 неделями трудно сказать определенно: какие-то из них, скорее всего, являются ошибочными, но другие корректно отражают случаи преждевременных родов.

Следующая функция выделяет самые большие (largest) значения из объекта FreqTab:

def largest(ftab, n=10):

    return ftab[-n:]

Вот самые продолжительные сроки беременности в датасете:

largest(ftab_length)

prglngth

40    1116

41     587

42     328

43     148

44      46

45      10

46       1

47       1

48       7

50       2

Name: prglngth, dtype: int64

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

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

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

Первые дети

Сравним теперь распределение продолжительности беременности для женщин, рожающих впервые, и всех остальных. Для выбора строк, представляющих первых (firsts) и других (others) детей, можно использовать метод query:

firsts = live.query("birthord == 1")

others = live.query("birthord != 1")

Создадим также объекты FreqTab для продолжительности беременности женщин из каждой группы:

ftab_first = FreqTab.from_seq(firsts["prglngth"], name="первые")

ftab_other = FreqTab.from_seq(others["prglngth"], name="другие")

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

def two_bar_plots(ftab1, ftab2, width=0.45):

    ftab1.bar(align="edge", width=-width)

    ftab2.bar(align="edge", width=width, alpha=0.5)

Вот как они выглядят:

two_bar_plots(ftab_first, ftab_other)

decorate(xlabel="Недели", ylabel="Частота", xlim=[20, 50])

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

firsts["prglngth"].count(), others["prglngth"].count()

(4413, 4735)

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

first_mean = firsts["prglngth"].mean()

other_mean = others["prglngth"].mean()

first_mean, other_mean

(38.60095173351461, 38.52291446673706)

Но разница составляет всего 0.078 недели, что равняется приблизительно 13 часам:

diff = first_mean - other_mean

diff, diff * 7 * 24

(0.07803726677754952, 13.11026081862832)

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

• Возможно, существует реальная разница в средней продолжительности беременности при рождении первых и других детей.

• Явная разница, наблюдаемая в этом датасете, может быть результатом предвзятости в процессе составления выборки, то есть при отборе респондентов для опроса.

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

• Очевидная разница может быть результатом случайных вариаций в процессе составления выборки.

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

Размер эффекта

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

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

diff / live["prglngth"].mean() * 100

0.20237586646738304

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

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

diff / live["prglngth"].std()

0.028877623375210403

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

Альтернативой является использование стандартного отклонения только одной из групп, но неясно, какой именно. Еще мы могли бы взять среднее значение двух стандартных отклонений, но если группы различаются по размеру, то это придало бы слишком большой вес одной группе и недостаточный — другой.

Распространенное решение — использование объединенного стандартного отклонения (pooled standard deviation), которое представляет собой квадратный корень из объединенной дисперсии. Последняя равна взвешенной сумме дисперсий в отдельных группах. Чтобы вычислить этот показатель, начнем с дисперсий:

group1, group2 = firsts["prglngth"], others["prglngth"]

 

v1, v2 = group1.var(), group2.var()

Вот взвешенная сумма с размерами групп в качестве весов:

n1, n2 = group1.count(), group2.count()

pooled_var = (n1 * v1 + n2 * v2) / (n1 + n2)

Наконец, вот объединенное стандартное отклонение:

np.sqrt(pooled_var)

2.7022108144953862

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

firsts["prglngth"].std(), others["prglngth"].std()

(2.7919014146687204, 2.6158523504392375)

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

def cohen_effect_size(group1, group2):

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

 

    v1, v2 = group1.var(), group2.var()

    n1, n2 = group1.count(), group2.count()

    pooled_var = (n1 * v1 + n2 * v2) / (n1 + n2)

 

    return diff / np.sqrt(pooled_var)

А вот размер эффекта для разницы в средней продолжительности беременности:

cohen_effect_size(firsts["prglngth"], others["prglngth"])

0.028879044654449834

В данном примере разница составляет 0.029 стандартного отклонения, это мало. Для сравнения: разница в росте между мужчинами и женщинами составляет около 1.7 стандартного отклонения.

Представление результатов

Мы рассмотрели несколько способов описания разницы в продолжительности беременности (если таковая имеется) при рождении первых и последующих детей. Как представить полученные результаты?

Ответ зависит от того, кто задает вопрос. Ученого может интересовать любой (реальный) эффект, каким бы незначительным он ни был. Врача могут интересовать только те эффекты, которые практически значимы, то есть различия, существенные для практики. Беременную женщину могут заинтересовать результаты, имеющие к ней отношение, например вероятность преждевременных или поздних родов.

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

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

Глоссарий

Распределение (distribution)

Набор значений и частота появления каждого значения в датасете.

Таблица частот (frequency table)

Сопоставление значений и частот.

Частота (frequency)

Количество появлений значения в выборке.

Скошенный (skewed)

Распределение называется скошенным, если оно асимметрично. При этом экстремальные величины распространяются дальше в одном направлении, чем в другом.

Мода (mode)

Наиболее часто встречающаяся в выборке величина или одна из наиболее часто встречающихся величин.

Равномерное распределение (uniform distribution)

Распределение, в котором все величины имеют одинаковую частоту.

Выброс (outlier)

Экстремальная величина в распределении.

Стандартизированный (standardized)

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

Объединенное стандартное отклонение (pooled standard deviation)

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

Размер эффекта Коэна (Cohen’s effect size)

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

Практически значимый (practically significant)

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

Упражнения

Для выполнения упражнений этой главы загрузите файл данных NSFG с ответами женщин-респондентов, по одной строке на каждую респондентку. Инструкции по загрузке данных и кодбука приведены в Jupyter-блокноте для этой главы.

Модуль nsfg.py содержит функцию, которая считывает файл ответов респонденток, очищает некоторые переменные и возвращает объект DataFrame:

from nsfg import read_fem_resp

 

resp = read_fem_resp()

resp.shape

(7643, 3092)

Данный DataFrame содержит 3092 столбца, но мы будем использовать лишь некоторые из них.

Упражнение 2.1

Начнем с переменной totincr, которая содержит общий доход семьи респондентки, закодированный значением от 1 до 14. Обратитесь к кодбуку для файла рес­пондентов, чтобы узнать, какой уровень дохода соответствует каждому значению.

Создайте объект FreqTab для распределения этой переменной и отобразите его при помощи столбчатой диаграммы.

Упражнение 2.2

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

Используйте функцию largest, чтобы найти наибольшие значения переменной parity. Существуют ли значения, которые вам кажутся ошибочными?

Упражнение 2.3

Выясним, рожают ли женщины с более высоким доходом больше детей. Используйте метод query, чтобы выбрать респонденток с самым высоким доходом (уровень 14). Постройте диаграмму частот переменной parity только для респонденток с высоким доходом.

Сравните среднее значение переменной parity для респонденток с высоким уровнем дохода и всех остальных респонденток.

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

Означают ли эти результаты, что у людей с более высоким доходом больше детей, или вы можете предложить другое объяснение имеющемуся различию?

---