Чем больше я узнаю людей, тем больше мне нравится моя собака.
—Марк Твен

В предыдущих сериях постов для начинающих из ремикса книги Генри Гарнера «Clojure для исследования данных» (Clojure for Data Science) на языке Python мы рассмотрели методы описания выборок с точки зрения сводных статистик и методов статистического вывода из них параметров популяции. Такой анализ сообщает нам нечто о популяции в целом и о выборке в частности, но он не позволяет нам делать очень точные утверждения об их отдельных элементах. Это связано с тем, что в результате сведения данных всего к двум статистикам - среднему значению и стандартному отклонению - теряется огромный объем информации.

Нам часто требуется пойти дальше и установить связь между двумя или несколькими переменными либо предсказать одну переменную при наличии другой. И это подводит нас к теме данной серии из 5 постов - исследованию корреляции и регрессии. Корреляция имеет дело с силой и направленностью связи между двумя или более переменными. Регрессия определяет природу этой связи и позволяет делать предсказания на ее основе.

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

О данных

В этой серии постов используются данные, любезно предоставленные компанией Guardian News and Media Ltd., о спортсменах, принимавших участие в Олимпийских Играх 2012 г. в Лондоне. Эти данные изначально были взяты из блога газеты Гардиан.

Обследование данных

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

Файл all-london-2012-athletes.tsv достаточно небольшой. Мы можем обследовать данные при помощи pandas, как мы делали в первой серии постов «Python, исследование данных и выборы», воспользовавшись функцией read_csv:

def load_data():
    return pd.read_csv('data/ch03/all-london-2012-athletes-ru.tsv', '\t')
                                            
def ex_3_1():
    '''Загрузка данных об участниках 
       олимпийских игр в Лондоне 2012 г.'''
    return load_data()

Если выполнить этот пример в консоли интерпретатора Python либо в блокноте Jupyter, то вы должны увидеть следующий ниже результат:

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

ФИО атлета
страна, за которую он выступает
возраст, лет
рост, см.
вес, кг.
пол «М» или «Ж»
дата рождения в виде строки
место рождения в виде строки (со страной)
число выигранных золотых медалей
число выигранных серебряных медалей
число выигранных бронзовых медалей
всего выигранных золотых, серебряных и бронзовых медалей
вид спорта, в котором он соревновался
состязание в виде списка, разделенного запятыми

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

Визуализация данных

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

def ex_3_2():
    '''Визуализация разброса значений 
       роста спортсменов на гистограмме'''
    df = load_data()
    df['Рост, см'].hist(bins=20)
    plt.xlabel('Рост, см.')
    plt.ylabel('Частота')
    plt.show()

Этот пример сгенерирует следующую ниже гистограмму:

Как мы и ожидали, данные приближенно нормально распределены. Средний рост спортсменов составляет примерно 177 см. Теперь посмотрим на распределение веса олимпийских спортсменов:

def ex_3_3():
    '''Визуализация разброса значений веса спортсменов'''
    df = load_data()
    df['Вес'].hist(bins=20)
    plt.xlabel('Вес')
    plt.ylabel('Частота')
    plt.show()

Приведенный выше пример сгенерирует следующую ниже гистограмму:

Данные показывают четко выраженную асимметрию. Хвост с правой стороны намного длиннее, чем с левой, и поэтому мы говорим, что асимметрия - положительная. Мы можем оценить асимметрию данных количественно при помощи функции библиотеки pandas skew:

def ex_3_4():
    '''Вычисление асимметрии веса спортсменов'''
    df = load_data()
    swimmers = df[ df['Вид спорта'] == 'Swimming']
    return swimmers['Вес'].skew()

0.23441459903001483

К счастью, эта асимметрия может быть эффективным образом смягчена путем взятия логарифма веса при помощи функции библиотеки numpy np.log:

def ex_3_5():
    '''Визуализация разброса значений веса спортсменов на
       полулогарифмической гистограмме с целью удаления 
       асимметрии'''
    df = load_data()
    df['Вес'].apply(np.log).hist(bins=20)
    plt.xlabel('Логарифмический вес')
    plt.ylabel('Частота')
    plt.show()

Этот пример сгенерирует следующую ниже гистограмму:

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

Логнормальное распределение

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

Логарифм показывает степень, в которую должно быть возведено фиксированное число (основание) для получения данного числа. Изобразив логарифмы на графике в виде гистограммы, мы показали, что эти степени приближенно нормально распределены. Логарифмы обычно берутся по основанию 10 или основанию e, трансцендентному числу, приближенно равному 2.718. В функции библиотеки numpy np.log и ее инверсии np.exp используется основание e. Выражение log_e также называется натуральным логарифмом, или ln, из-за свойств, делающих его особенно удобным в исчислении.

Логнормальное распределение обычно имеет место в процессах роста, где темп роста не зависит от размера. Этот феномен известен как закон Джибрэта, который был cформулирован в 1931 г. Робертом Джибрэтом, заметившим, что он применим к росту фирм. Поскольку темп роста пропорционален размеру, более крупные фирмы демонстрируют тенденцию расти быстрее, чем фирмы меньшего размера.

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

С тех пор выяснилось, что закон Джибрэта применим к большому числу ситуаций, включая размеры городов и, согласно обширному математическому ресурсу Wolfram MathWorld, к количеству слов в предложениях шотландского писателя Джорджа Бернарда Шоу.

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

Визуализация корреляции

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

def swimmer_data():
    '''Загрузка данных роста и веса только олимпийских пловцов'''
    df = load_data()
    return df[df['Вид спорта'] == 'Swimming'].dropna()

def ex_3_6():
    '''Визуализация корреляции между ростом и весом'''
    df = swimmer_data()
    xs = df['Рост, см']
    ys = df['Вес'].apply( np.log )
    pd.DataFrame(np.array([xs,ys]).T).plot.scatter(0, 1, s=12, grid=True)
    plt.xlabel('Рост, см.')
    plt.ylabel('Логарифмический вес')
    plt.show()

Этот пример сгенерирует следующий ниже график:

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

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

Генерирование джиттера

Поскольку каждое значение округлено до ближайшего сантиметра или килограмма, то значение, записанное как 180 см, на самом деле может быть каким угодно между 179.5 и 180.5 см, тогда как значение 80 кг на самом деле может быть каким угодно между 79.5 и 80.5 кг. Для создания случайных искажений, мы можем добавить случайные помехи в каждую точку данных роста в диапазоне между -0.5 и 0.5 и в том же самом диапазоне проделать с точками данных веса (разумеется, это нужно cделать до того, как мы возьмем логарифм значений веса):

def jitter(limit):
    '''Генератор джиттера (произвольного сдвига точек данных)'''
    return lambda x: random.uniform(-limit, limit) + x

def ex_3_7():
    '''Визуализация корреляции между ростом и весом с джиттером'''
    df = swimmer_data()
    xs = df['Рост, см'].apply(jitter(0.5))
    ys = df['Вес'].apply(jitter(0.5)).apply(np.log)
    pd.DataFrame(np.array([xs,ys]).T).plot.scatter(0, 1, s=12, grid=True)
    plt.xlabel('Рост, см.')
    plt.ylabel('Логарифмический вес')
    plt.show()

График с джиттером выглядит следующим образом:

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

Ковариация

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

Если у нас имеется два ряда чисел, X и Y, то их отклонения от среднего значения составляют:

Здесь x_i — это значение X с индексом i, y_i — значение Y с индексом i, x̅ — среднее значение X, и y̅ — среднее значение Y. Если X и Y проявляют тенденцию изменяться вместе, то их отклонения от среднего будет иметь одинаковый знак: отрицательный, если они — меньше среднего, положительный, если они больше среднего. Если мы их перемножим, то произведение будет положительным, когда у них одинаковый знак, и отрицательным, когда у них разные знаки. Сложение произведений дает меру тенденции этих двух переменных отклоняться от среднего значения в одинаковом направлении для каждой заданной выборки.

Ковариация определяется как среднее этих произведений:

На чистом Python ковариация вычисляется следующим образом:

def covariance(xs, ys):
    '''Вычисление ковариации (несмещенная, т.е. n-1)'''
    dx = xs - xs.mean() 
    dy = ys - ys.mean()
    return (dx * dy).sum() / (dx.count() - 1)

В качестве альтернативы, мы можем воспользоваться функцией pandas cov:

df['Рост, см'].cov(df['Вес'])

1.3559273321696459

Ковариация роста и логарифма веса для наших олимпийских пловцов равна 1.356, однако это число сложно интерпретировать. Единицы измерения здесь представлены произведением единиц на входе.

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

^{Стандартная оценка, англ. standard score, также z-оценка — это относительное число стандартных отклонений, на которые значение переменной отстоит от среднего значения. Положительная оценка показывает, что переменная находится выше среднего, отрицательная — ниже среднего. Это безразмерная величина, получаемая при вычитании популяционного среднего из индивидуальных значений и деления разности на популяционное стандартное отклонение.}

Корреляция Пирсона

Корреляция Пирсона часто обозначается переменной r и вычисляется следующим образом, где отклонения от среднего dx_i и dy_i вычисляются как и прежде:

Поскольку для переменных X и Y стандартные отклонения являются константными, уравнение может быть упрощено до следующего, где σ_x и σ_y — это стандартные отклонения соответственно X и Y:

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

Ранее мы уже написали функции для вычисления стандартного отклонения. В сочетании с нашей функцией с вычислением ковариации получится следующая ниже имплементация корреляции Пирсона:

def variance(xs):
    '''Вычисление корреляции,
       несмещенная дисперсия при n <= 30'''
    x_hat = xs.mean()
    n = xs.count()
    n = n - 1 if n in range( 1, 30 ) else n  
    return sum((xs - x_hat) ** 2) / n

def standard_deviation(xs):
    '''Вычисление стандартного отклонения'''
    return np.sqrt(variance(xs))

def correlation(xs, ys): 
    '''Вычисление корреляции'''
    return covariance(xs, ys) / (standard_deviation(xs) * 
                                 standard_deviation(ys))

В качестве альтернативы мы можем воспользоваться функцией pandas corr:

df['Рост, см'].corr(df['Вес'])

Поскольку стандартные оценки безразмерны, то и коэффициент корреляции r тоже безразмерен. Если r равен -1.0 либо 1.0, то переменные идеально антикоррелируют либо идеально коррелируют.

Правда, если r = 0, то с необходимостью вовсе не следует, что переменные не коррелируют. Корреляция Пирсона измеряет лишь линейные связи. Как продемонстрировано на следующих графиках, между переменными может существовать еще некая нелинейная связь, которую r не объясняет:

Отметим, что корреляция центрального примера не определена, потому что стандартное отклонение y = 0. Поскольку наше уравнение для r содержало бы деление ковариации на 0, то результат получается бессмысленным. В этом случае между переменными не может быть никакой корреляции; y всегда будет иметь среднее значение. Простое обследование стандартных отклонений это подтвердит.

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

def ex_3_8():
    '''Вычисление корреляции средствами pandas
       на примере данных роста и веса'''
    df = swimmer_data()
    return df['Рост, см'].corr( df['Вес'].apply(np.log))

0.86748249283924894

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

Выборочный r и популяционный ρ

Аналогично среднему значению и стандартному отклонению, коэффициент корреляции является сводной статистикой. Он описывает выборку; в данном случае, выборку спаренных значений: роста и веса. Коэффициент корреляции известной выборки обозначается буквой r, тогда как коэффициент корреляции неизвестной популяции обозначается греческой буквой ρ (рхо).

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

Даже в допустимой популяции — такой как пловцы, выступавшие на недавних Олимпийских играх, — наша выборка коэффициента корреляции является всего лишь одной из многих потенциально возможных. То, насколько мы можем доверять нашему r, как оценке параметра ρ, зависит от двух факторов:

Размера выборки
Величины r

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

Проверка статистических гипотез

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

В первую очередь, мы должны сформулировать две гипотезы, нулевую гипотезу и альтернативную:

$H_1∶ρ\ne 0$

H₀ - это гипотеза, что корреляция в популяции нулевая. Другими словами, наше консервативное представление состоит в том, что измеренная корреляция целиком вызвана случайной ошибкой при отборе.

H₁ - это альтернативная возможность, что корреляция в популяции не нулевая. Отметим, что мы не определяем направление корреляции, а только что она существует. Это означает, что мы выполняем двустороннюю проверку.

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

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

Мы можем снова воспользоваться t-распределением и вычислить t-статистику:

В приведенной формуле df — это степень свободы наших данных. Для проверки корреляции степень свободы равна n - 2, где n — это размер выборки. Подставив это значение в формулу, получим:

В итоге получим t-значение 102.21. В целях его преобразования в p-значение мы должны обратиться к t-распределению. Библиотека scipy предоставляет интегральную функцию распределения (ИФР) для t-распределения в виде функции stats.t.cdf, и комплементарной ей (1-cdf) функции выживания stats.t.sf. Значение функции выживания соответствует p-значению для односторонней проверки. Мы умножаем его на 2, потому что выполняем двустороннюю проверку:

def t_statistic(xs, ys):
    '''Вычисление t-статистики'''
    r = xs.corr(ys)  # как вариант, correlation(xs, ys)
    df = xs.count() - 2
    return r * np.sqrt(df / 1 - r ** 2)

def ex_3_9():
    '''Выполнение двухстороннего t-теста'''
    df = swimmer_data()
    xs = df['Рост, см']
    ys = df['Вес'].apply(np.log)
    t_value = t_statistic(xs, ys)
    df = xs.count() - 2 
    p = 2 * stats.t.sf(t_value, df)  # функция выживания 
    return {'t-значение':t_value, 'p-значение':p}

{'p-значение': 1.8980236317815443e-106, 't-значение': 25.384018200627057}

P-значение настолько мало, что в сущности равно 0, означая, что шанс, что нулевая гипотеза является истинной, фактически не существует. Мы вынуждены принять альтернативную гипотезу о существовании корреляции.

Интервалы уверенности

Установив, что в более широкой популяции, безусловно, существует корреляция, мы, возможно, захотим количественно выразить диапазон значений, внутри которого, как мы ожидаем, будет лежать параметр ρ, вычислив для этого интервал уверенности. Как и в случае со средним значением в предыдущей серии постов, интервал уверенности для r выражает вероятность (выраженную в %), что параметр ρ популяции находится между двумя конкретными значениями.

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

Приведенный выше график показывает отрицательно скошенное распределение r-выборок для параметра ρ, равного 0.6.

К счастью, трансформация под названием z-преобразование Фишера стабилизирует дисперсию r по своему диапазону. Она аналогична тому, как наши данные о весе спортсменов стали нормально распределенными, когда мы взяли их логарифм.

Уравнение для z-преобразования следующее:

Стандартная ошибка z равна:

Таким образом, процедура вычисления интервалов уверенности состоит в преобразовании r в z с использованием z-преобразования, вычислении интервала уверенности в терминах стандартной ошибки SE_z и затем преобразовании интервала уверенности в r.

В целях вычисления интервала уверенности в терминах SE_z, мы можем взять число стандартных отклонений от среднего, которое дает нам требуемый уровень доверия. Обычно используют число 1.96, так как оно является числом стандартных отклонений от среднего, которое содержит 95% площади под кривой. Другими словами, 1.96 стандартных ошибок от среднего значения выборочного r содержит истинную популяционную корреляцию ρ с 95%-ой определенностью.

Мы можем убедиться в этом, воспользовавшись функцией scipy stats.norm.ppf. Она вернет стандартную оценку, связанную с заданной интегральной вероятностью в условиях односторонней проверки.

Однако, как показано на приведенном выше графике, мы хотели бы вычесть ту же самую величину, т.е. 2.5%, из каждого хвоста с тем, чтобы 95%-й интервал уверенности был центрирован на нуле. Для этого при выполнении двусторонней проверки нужно просто уменьшить разность наполовину и вычесть результат из 100%. Так что, требуемый уровень доверия в 95% означает, что мы обращаемся к критическому значению 97.5%:

def critical_value(confidence, ntails): # ДИ и число хвостов
    '''Расчет критического значения путем
       вычисления квантиля и получения 
       для него нормального значения'''
    lookup = 1 - ((1 - confidence) / ntails) 
    return stats.norm.ppf(lookup, 0, 1)  # mu=0, sigma=1

critical_value(0.95, 2)

1.959963984540054

Поэтому наш 95%-й интервал уверенности в z-пространстве для ρ задается следующей формулой:

Подставив в нашу формулу z_r и SE_z, получим:

Для r=0.867 и n=859 она даст нижнюю и верхнюю границу соответственно 1.137 и 1.722. В целях их преобразования из z-оценок в r-значения, мы используем следующее обратное уравнение z-преобразования:

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

def z_to_r(z):
    '''Преобразование z-оценки обратно в r-значение'''
    return (np.exp(z*2) - 1) / (np.exp(z*2) + 1)

def r_confidence_interval(crit, xs, ys): 
    '''Расчет интервала уверенности
       для критического значения и данных'''
    r   = xs.corr(ys)
    n   = xs.count()
    zr  = 0.5 * np.log((1 + r) / (1 - r)) 
    sez = 1 / np.sqrt(n - 3)
    return (z_to_r(zr - (crit * sez))), (z_to_r(zr + (crit * sez)))

def ex_3_10():
    '''Расчет интервала уверенности
       на примере данных роста и веса'''
    df = swimmer_data()
    X = df['Рост, см']
    y = df['Вес'].apply(np.log)
    interval = r_confidence_interval(1.96, X, y) 
    print('Интервал уверенности (95%):', interval)

Интервал уверенности (95%): (0.8499088588880347, 0.8831284878884087)

В результате получаем 95%-й интервал уверенности для ρ, расположенный между 0.850 и 0.883. Мы можем быть абсолютно уверены в том, что в более широкой популяции олимпийских пловцов существует сильная положительная корреляция между ростом и весом.

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

В следующем посте, посте №2, будет рассмотрена сама тема серии - регрессия и приемы оценивания ее качества.

Python, корреляция и регрессия: часть 1