Как найти связь между переменными

Анализ связи между двумя переменными

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

1) показать, что существует эмпирическая
взаимосвязь между переменными;

2) исключить возможность
обратного влияния Y
на X;

3) убедиться, что взаимосвязь между
переменными не может быть объяснена
зависимостью этих переменных от какой-то
дополнительной переменной (или
переменных).

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

Пусть мы располагаем
совокупностью данных о занятиях
физзарядкой и образовании для выборки
горожан. Для простоты мы предположим,
что обе переменные имеют лишь два уровня:
высокий и низкий. Так как данные об
образовании исходно разбиты на большее
количество категорий, нам придется их
перегруппировать, разбив весь диапазон
значений на два класса. Предположим, мы
выберем в качестве граничного значения
10 лет обучения, так что люди, получившие
неполное среднее и среднее образование,
попадут в «низкую» градацию, а остальные
– в «высокую». (Это, конечно, большое
огрубление, но мы используем его из
соображений простоты.) Для занятий
физическим^ упражнениями мы соответственно
воспользуемся двумя категориями –
«делают физзарядку» и «не делают
физзарядку». Таблица
8.3 показывает, как могло
бы выглядеть совместное распределение
этих двух переменных.

Таблица 8.3 Взаимосвязь
между уровнем образования и занятиями
физкультурой

	Уровень образования	Всего
Высокий	Низкий
Делают зарядку	50	200	250
Не делают зарядку	205	45	250
Всего	255	245	500

В таблице 8.3
два столбца (для образования) и две
строки (для занятий физкультурой),
следовательно, размерность этой таблицы
2х2. Кроме того, имеются дополнительные
крайний столбец и крайняя строка
(маргиналы таблицы), указывающие общее
количество наблюдений в данной строке
или в столбце. В правом нижнем углу
указана общая сумма, т. е. общее число
наблюдений в выборке. Не давшие ответа
уже исключены (для реальных данных их
число также стоит указать, но не в
таблице, а в подтабличной сноске). Заметим
здесь, что многие исследователи при
построении таких таблиц пользуются
неписаным правилом: для той переменной,
которую полагают независимой, отводится
верхняя строка (горизонталь), а зависимую
располагают «сбоку», по вертикали
(разумеется, соблюдение этого правила
не является обязательным и ничего с
точки зрения анализа не меняет).

Обычно характер взаимоотношений
между переменными в небольшой таблице
можно определить даже «на глазок»,
сравнивая числа в столбцах или строках.
Еще легче это сделать, если вместо
абсолютных значений стоят проценты.
Чтобы перевести абсолютные частоты,
указанные в клетках таблицы, в проценты,
нужно разделить их на маргинальные
частоты и умножить на 100. Если делить на
маргинал столбца, мы получим процент
по столбцу. Например,

50/255х100 = 19,6%, т. е. 19,6% имеющих
низкий уровень образования делают
зарядку (но не наоборот!). Если делить
на маргинал строки, то мы получим другую
величину –процент по
строке. В частности,
можно заметить, что 80% делающих зарядку,
составляют люди с высоким уровнем
образования(200/250х100). Деление на общую
численность выборки дает общий процент.
Так, всего в выборке 50% людей, делающих
зарядку.

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

Отметим специально, что
элементарная таблица сопряженности
размерности 2х2 – это минимально необходимое
условие для вывода о наличии взаимосвязи
двух переменных. Знания о распределении
зависимой переменной недостаточно.
Нельзя, например, утверждать, будто из
того, что 75% детей-первенцев имеют
интеллект выше среднего, а 25% – средний
и более низкий, следует зависимость
между порядком рождения и интеллектом.
Необходимо проанализировать и
распределение показателей интеллекта
для детей-непервенцев. Варьировать
должна не только зависимая, но и
независимая переменная.

Таблица 8.4 Oбщая
форма таблицы сопряженности размерности
2х2

	Переменная Х	Всего
0	1
Переменная У	0	а	б	А+б
1	с	д	С+д
Всего	А+с	Б+д	n

довольно ясно показывает, что всякое
приращение в размерах партийной кассы
(сдвиг вправо по оси Х)влечет за
собой увеличение парламентского
представи­тельства (сдвиг вверх по
оси ординат). Между переменнымиХиYсуществуетлинейное отношение:если одна переменная возрастает по
величине, то это же происходит и с другой.
Помимо указания на природу отношения
двух перемен­ных, диаграмма нарисунке
21позволяет также сделать некоторые
предположе­ния обинтенсивности,силе этого отношения. Очевидно, что чем
более ком­пактно, «скученно»
располагаются точки-наблюдения вокруг
пунктирной пря­мой линии (описывающейидеальное линейное отношение Х иY), тем сильнее зависимость. Нарисунке
22приведены еще три диаграммы
рассеивания.

Очевидно, что на рисунке 22акакая-либо
связь между x
и
y попросту
отсут­ствует. Нарисунке
22бвоображаемая прямая линия
(отмечена пунктиром) пе­ресекла бы
диаграмму сверху вниз, из левого верхнего
в правый нижний угол. Иными словами,
линейная связь в этом случае имеет
обратное направление:

чем больше X,тем меньше зависимая переменная У.
Заметим также, что «куч­ность»
расположения точек вдоль воображаемой
прямойна рисунке 226не очень велика, а значит и связь
(корреляция) между переменными не только
обратная, отрицательная, но еще и не
очень сильная, умеренная. Наконец,на
рисунке 22в зависимую и независимую
переменную связывает явно нелинейное
отноше­ние: воображаемый график
нисколько не похож на прямую линию и
напомина­ет скорее параболу¹².
Отметим, что методы анализа, о которых
сейчас пойдет речь, не годятся для этого
нелинейного случая, так как обычная
формула для подсчета коэффициента
корреляции даст нулевое значение, хотя
связь между переменными существует.

Существует обобщенный показатель,
позволяющий оценить, насколько связь
между переменными приближается к
линейному функциональному отношению,
которое на диаграмме рассеивания
выглядит как прямая линия. Это коэффици­ент
корреляции,измеряющий тесноту связи
между переменными, т. е. их тен-денцию
изменяться совместно. Как и в рассмотренных
выше мерах связи каче-ственных признаков,
коэффициент корреляции позволяет
оценить возможность предсказания
значений зависимой переменной по
значениям независимой. Об­щая формула
для вычислениякоэффициента корреляции
Пирсонавключает в себя величинуковариациизначенийXиY.
Эта величина (s_xy)характеризует со­вместное изменение
значений двух переменных. Она задается
как сумма произ-

¹²
Именно так обычно выглядит зависимость
между благожелательностью установ­ки
по отношению к некоторому объекту
(X)
и интенсивностью установки (Y):
люди, занимающие крайне благожелательную
или крайне неблагожелательную позицию
в каком-то вопросе, обычно оценивают
свои убеждения как более выраженные и
ин­тенсивные, чем те люди, чьи установки
лежат в области середины, «нейтральных»
значений шкалы.ведений отклонений
наблюдаемых значений Х
и У от средних X
и Y
соответ­ственно, т. е. Eⁱ⁼¹_n
(X_i
– Х )(Y_i
–
Y
),
деленная на количество наблюдений.
Чтобы понять «физический смысл»
ковариации, достаточно обратить внимание
на сле­дующее ее свойство: если для
какого-то объекта
i
в выборке оба значения
—X,
и
Y,—
окажутся высокими, то и произведение
(Х_i
– Х)
на (Y_i
– Y)
будет боль­шим и положительным. Если
оба значения (по Х
и по Y)
низки, то произведение двух отклонений,
т. е. двух отрицательных чисел, также
будет положительным. Таким образом,
если линейная связь Х
и Y
положительна и велика, сумма таких
произведений для всех наблюдений также
будет положительна. Если связь между Х
и У обратная,
томногимположительным
отклонениям по Х
будут
соответ­ствовать отрицательные
отклонения по Y,
т. е. сумма отрицательных произведе­ний
отклонений будет отрицательной.

Наконец, при
отсутствии систематической связи
произведения будут иногда положительными,
иногда отрицательными, а их сумма (и,
следовательно, ковариация Х
и Y)
будет, в пределе, равна нулю. Таким
образом, ко вариация показы­вает
величину и направление связи, совместного
изменения Х
и У. Если разде­лить ковариацию
s_xy
на стандартные отклонения
s_x
и s_y
(чтобы избавиться от влияния масштаба
шкал, в которых измеряются X
и Y),
то мы получим искомую формулу коэффициента
корреляции Пирсона
(r_xy):

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

Несмотря на несколько устрашающий вид,
расчетная формула очень проста. Для
«ручного» вычисления r_xy,вам понадобятся лишь пять величин: суммы
значений поXиY(Eⁱ⁼¹_nХ_iиEⁱ⁼¹_nY_i),
суммы квадратов значений поХ и
Y (Eⁱ⁼¹_nХ²иEⁱ⁼¹_nY²),
суммы произведенийХ иYпо всем
объектам-«случаям» (Eⁱ⁼¹_n
X_iY_i).

В таблице 8.11приведены данные о максимальных дневных
и ночных температурах, зарегистрированных
в 10городах¹³.

Просуммировав значения в столбцах, мы
получим E¹⁰ _i=1Х_i
=258иE¹⁰ _i=1 Y_i=156.

Возведя каждое из значений Х и
Yв квадрат и просуммировав, мы
найдем, чтоE¹⁰ _i=1Х_i²
= 7180иE¹⁰ _i=1
Y_i²
= 2962.Сумма попарных произведенийХ_iиY_i(E¹⁰
_i=1 X_iY_i)
составит 4359.Вы можете
самостоятельно убедиться в том, что
подстановка всех значений в расчетную
формулу даст (надеюсь) величинуr_xy=0,91.
Иными словами, корреляция между дневными
и ночными температурами воздуха очень
высока, но все же отлична от
1,0(коэффициент корреляции может
меняться в пределах от-1,0 до
+1,0).Это отличие, вероятно, объясняется
влиянием других факторов (продолжительность
дня и ночи, облачность, географическое
положение и т. п.). Судя по полученной
величине корреляции, знание дневных
температур позволяетпредсказыватьночные температуры с очень высокой
точностью, но не безошибочно.

¹³Погода (Гидрометцентр
РФ) //Сегодня.
1994. 23авг.

Величина, которая равна квадрату
коэффициента корреляции Пирсона, т. е.
г, имеет ряд интересных статистических
свойств. Отметим сейчас, что r²являет­ся ПУО-мерой связи, подобной
обсуждавшимся выше (см. с.
176—179).Мож­но показать, что она
характеризует ту долю дисперсии значенийY, которая объяс­няется наличием
корреляции междуХ и Y.(Естественно, величинаr²будет
все­гда положительной и не может
превзойти по абсолютной величине
коэффициент корреляции.)¹⁴Та
часть разброса в значенияхY, которая
не может быть пред­сказана по значениям
X,—этодисперсия
ошибки нашего прогноза,т. е.
1 –r².Необъясненный
разброс в значенияхУприсутствует
в том случае, когда при равных уровнях
Х (ср., например, дневные температуры
в Варшаве и Бонне изтаблицы
8.11) сохраняются различия в
значенияхY.

Коэффициент корреляции позволяет
оценить степень связи между переменны­ми.
Однако этого недостаточно для того,
чтобы непосредственно преобразовы­вать
информацию, относящуюся к одной
переменной, в оценки другой пере­менной.
Допустим, мы выяснили, что коэффициент
корреляции между пере­менными «величина
партийного бюджета» и «число мест в
парламенте» равен 0,8.Можем ли мы теперь предсказать, сколько
мест в парламенте полу­чит партия,
годовой бюджет которой равен
100млн рублей? Похоже, что знание
величины коэффициента корреляции нам
здесь не поможет. Однако мы можем
вспомнить, что коэффициент корреляции
—это еще и оценка соответствия
раз­броса наших наблюдений той
идеальной модели линейного функционального
отношения, которое на рассмотренных
выше диаграммах рассеивания (см.рис.
21—22)представлено пунктирными
прямыми. Эти линии называютлиниями
регрессии.

Если бы все наблюдения аккуратно
«укладывались» на линию регрессии, то
для предсказания значения зависимой
переменной достаточно было бы восстано
вить перпендикуляр к оси
Y из той точки прямой, которая
соответствует известному значению
X.

¹⁴Подробный анализ
можно найти в большинстве руководств
по прикладной статисти­ке. Здесь
мы ограничимся обсуждением общей логикиоценки объясненной дисперсии.

На рисунке 23показано, как можно было бы графически
определить ожидае­мые значенияYдля гипотетического примера с партийной
кассой и местами в парламенте. (Разумеется,
найти искомое значение У можно и без
линейки, с помощью вычислений, если
известен угол наклона регрессионной
прямой и точка пересечения с осью
ординат.)

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

Существуют специальные статистические
процедуры, которые позволяют най­ти
регрессионную прямую, максимально
соответствующую реальным данным.
Регрессионный анализ, таким образом,
дает возможность предсказывать
зна­ченияYпо значениям Х с
минимальным количеством ошибок.В
общем виде уравнение, описывающее прямую
линию регрессии Y
no X,выглядит
так:

Y=a_yx
+b_yx
X,

где Y—то предсказываемое значение по переменнойY(в только что рассмот­ренном
примере —количество
мест в парламенте),а —это точка, в которой прямая пересекает
осьY(т. е. значение Y для случая,
когдаX=0), и b—коэффи­циент
регрессии, т. е. наклон прямой. Часто
удобно измерять обе переменные не в
«сырых» шкалах, а в единицах отклонения
от среднего. Процедурастан­дартизации,т. е. перевода исходной шкалы в стандартные
Z-оценки,вам уже известна (см. с.
169).Преимущество использования
стандартизированных пе­ременных в
регрессионном анализе заключается в
том, что линия регрессии в этом случае
проходит через начало координат.
Стандартизованный коэффици­ент
регрессии (наклон прямой) обозначают
обычно греческой буквой b(либо лат. b*).Правда, в отличие от b-коэффициента,
bне позволяет
прямо заклю­чить, на какое количество
исходных единиц возрастет У при увеличенииХна одну единицу (например, насколько
увеличится число депутатских мандатов
при увеличении бюджета на 1млн рублей или насколько увеличивается
зара­ботная плата при увеличении
стажа работы на один год). С другой
стороны, bпозволяет сопоставить
влияние на независимую переменную
контрольных переменных, измеренных в
разных шкалах.

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

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

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

Далее, нужно помнить о том, что регрессионные
коэффициенты в общем слу­чае
асимметричны.Если мы решим, что это
У, а неХявляется независимой
переменной, то вполне можем рассчитать
другую по величине пару коэффициентов
—а_xуиb_xу.(Заметьте, что порядок букв в подстрочном
индексе значим: первой всегда идет
предсказываемая переменная, а второй
—предсказываю­щая.) Разумеется,
при выборе кандидатов в зависимые и
независимые перемен­ные также важны
не статистические, а содержательные
соображения.

Если вернуться к затронутой выше
взаимосвязи между линейной регрессией
и корреляцией,то здесь мы можем
сделать следующие дополнения. Пусть
все точ­ки-наблюдения аккуратно
размещены на регрессионной прямой.
Перед нами почти невероятный случай
абсолютной линейной зависимости. Зная,
например, что коэффициентb(нестандартизованный) равен
313,мы можем утверждать, что именно
такова величина воздействия переменнойXна зависимую перемен­нуюY.
Кроме того, мы можем точно сказать, что
единичная прибавка в величи­неXвызовет увеличениеYна ту же величину,
313(если, допустим,Х—стаж работы, а У —зарплата,
то с увеличением стажа на год зарплата
растет на313рублей). В
этом случае коэффициент корреляции
будет равен в точности 1,0, что
свидетельствует о сильном, «абсолютном»,
характере связи переменных. Различие
между предсказанными и наблюдаемыми
значениями в этом случае отсутствует.
Корреляция как мера точности прогноза
показывает, что ошибок предсказания
просто нет.

¹⁵Более детальные
сведения можно найти в статистической
литературе. Очень доступ­но проблема
излагается, в частности, в кн.:Гласе
Дж., Стенли Дж.Указ. соч. С.
123—141.Для тех же, кто захочет
осуществить «ручную» регрессию для
какого-либо из использованных примеров,
просто приведем формулы для вычисления
нестан­дартизированных коэффициентов
(обозначения те же, что и выше):

В действительности, однако, из-за влияния
других переменных и случайной выборочной
ошибки точки-наблюдения обычно лежат
выше или ниже прямой, которая, как
говорилось, является лишь наилучшим
приближением реальных данных. Коэффициент
корреляции Пирсона rи величинаr²по-прежнему слу­жат оценкой точности
прогноза, основанного на линии регрессии.
Вполне воз­можны ситуации, когда
коэффициент регрессии очень велик,
воздействиеXнаYпросто громадно, но корреляция низка
и, следовательно, точность прогноза
невелика. Нет ничего необычного и в
обратной ситуации: воздействиеХ на
Y относительно мало, а коэффициент
корреляции и объясненная дисперсия
очень велики. Посмотрев на приведенные
выше диаграммы рассеивания, можно лег­ко
уяснить себесмысл отношения между
корреляцией и регрессией:первая
имеет прямое отношение к «разбросанности»
точек наблюдения (чем выше «разбро­санность»,
тем нижеr²и ненадежнее
прогноз), тогда как коэффициент регрес­сии
описывает наклон, «крутизну» линии.
Однако существующее здесь разли­чие
не стоит и преувеличивать: регрессионный
коэффициент (наклон прямой) для
стандартизованных данных в точности
равен коэффициенту корреляции Пирсонаr¹⁶. ¹⁶
Легко понять, что при измерении в единицах
стандартного отклонения максимальная
связь (Р=
1,0)
соответствует ситуации, когда сдвигу
от начала координат в
1 ед.
стан­дартного отклонения по Х
соответствует увеличение У также на
1 ед.
стандартного от­клонения. Важно
заметить, что в случае стандартизированных
переменных (и только в этом случае)
коэффициенты регрессии
Y
no
X и
X noY
будут совпадать.

Предположим, что исследователь изучает
зависимость между образованием матери
(X)и образованием
детей (Y). Обе переменные измерены
как количе­ство лет, затраченных на
получение образования. Найдя достаточно
высокую корреляцию междуХи У—скажем, равную 0,71, —он
также находит коэффи­циенты регрессииаиbи устанавливает, чтоr²(называемый также коэффици­ентом
детерминации) в данном случае приближенно
равен 0,5.Это значит, что
доля вариации в значениях переменнойY(образование детей), объясненная
воздействием переменной-предиктораХ(материнское образование), составля­ет
около 50%общей дисперсии
предсказываемой переменной. Коэффициент
корреляции между переменными достаточно
велик и статистически значим даже для
не очень большой выборки. Следовательно,
обнаруженная взаимосвязь пе­ременных
не может быть объяснена случайными
погрешностями выборки. В пользу
предложенной исследователем причинной
гипотезы говорит и то об­стоятельство,
что альтернативная гипотеза
—образование детей влияет на
об­разовательный статус родителей
—крайне неправдоподобна и может
быть от­вергнута на основании
содержательных представлений о временной
упорядо­ченности событий. Однако все
еще не исключены те возможности, которые
мы обсуждали в параграфе, посвященном
методу уточнения. Иными словами, нам
следует считаться с вероятностью того,
что какая-то другая переменная (или
несколько переменных) определяют и
образование родителей и образование
детей (например, финансовые возможности
либо интеллект). Чтобы проверить такую
конкурирующую гипотезу, следует
рассчитать так называемуючастную
корреляцию.Логика расчета частной
корреляции совпадает с логикой построе­ния
частных таблиц сопряженности при
использовании метода уточнения.
По­строить частные таблицы сопряженности
для различных уровней контрольной
переменной в случае, когда переменные
измерены на интервальном уровне,—
это практически неразрешимая задача.

Чтобы убедиться в этом, достаточно
подсчитать, каким должно быть количество
таблиц уже при десяти-двенадцати
категориях каждой переменной. Расчет
коэффициента частной корреляции—
это простейшее средство уточнения
исходной причинной модели при введении
дополнительной переменной. Интерпретация
коэффициента частной корреля­ции не
отличается от интерпретации частных
таблиц сопряженности:частной
корреляцией называют корреляцию между
двумя переменными, когда статис­тически
контролируется, или «поддерживается
на постоянном уровне», тре­тья
переменная (набор переменных).

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

Мюллер и соавторы ¹⁷приводят
интересный пример использования
коэффици­ента частной корреляции. В
исследовании П. Риттербэнда и Р.
Силберстайна изучались студенческие
беспорядки 1968—1969гг. Одна
из гипотез заключа­лась в том, что
число нарушений дисциплины и демонстраций
протеста в стар­ших классах учебных
заведений связано с различиями показателей
академи­ческой успеваемости учащихся.
Корреляция между частотой «политических»
беспорядков и средней успеваемостью
оказалась отрицательной (хуже успева­емость
—больше беспорядков) и статистически
значимой(r = -0,36).Однако
еще более высокой была корреляция между
частотой беспорядков и долей чер­нокожих
учащихся (r=0,54). Исследователи решили
проверить, сохранится ли связь между
беспорядками и успеваемостью, если
статистически проконтроли­ровать
влияние расового состава учащихся.
Коэффициент частной корреляции частоты
беспорядков и успеваемостипри контролерасового состава учащихся оказался
равным нулю. Исходная корреляция между
беспорядками и успевае­мостью в данном
случае может быть описана причинной
моделью«ложной взаимосвязи» (см. рис.
19):наблюдаемые значения этих двух
переменных скор-релированы лишь потому,
что обе они зависят от третьей переменной
—доли чернокожих в общем количестве
учащихся. Чернокожие студенты, как
замети­ли исследователи, оказались
восприимчивее к предложенным самыми
актив­ными «политиканами» образцам
участия в политических беспорядках.
Кроме того, их успеваемость, помимо
всяких политических событий, была
устойчиво ниже, чем средняя успеваемость
белых.

Коэффициент частной корреляции между
переменными Х и. Y
при контроле дополнительной переменной
Z(т. е. при поддержании Z«на постоянном уров­не») обозначают
какr_xy.z.Для его
вычисления достаточно знать, величины
на­блюдаемых попарных корреляций
между переменнымиX, Yи Z
(т. е.простых корреляций —r_xy, r_yz, r_xz):

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

¹⁷ Mueller
J., Schuessler К., Costner
H. Statistical Reasoning in Sociology.
3rd ed. Boston:

Haighton Mifflin Co,
1977. P. 279—281.

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

Прежде всего вероятностная оценка
коэффициента корреляции подразумевает
оценку отношения к его случайной ошибке.
Удобная, хотя и не вполне надеж­ная
формула для вычисления ошибки коэффициента
корреляции (т_r),выглядит
так¹⁸:

Всегда полезно вычислить отношение
полученной величины rк его ошибке
(т. е.r/т_i).В использовавшемся
нами примере данных о погоде коэффициент
корреляции оказался равен
0,91,а его выборочная ошибка составляет:

Отношение rк m_rобозначаемое как t,составит (0,91/0,573)= 15,88.Разумеется, коэффициент, превосходящий
свою случайную ошибку почти в
16раз, может быть признан значимым
даже без построения доверительных
интервалов.

Когда значение г не столь близко к
единице и выборка невелика, нужно все
же проверить статистическую гипотезу
о равенстве rнулю в генеральной
совокуп­ности. Для этого нужно
определить tпо формуле:

где t —это величина так называемого
t-критерияСтьюдента (см. также гла­ву
4),г —выборочный
коэффициент корреляции,п
—объем выборки. Для ус­тановления
значимости вычисленной величины
t-критерияпользуются табли­цами
t-распределениядля (n – 2)степеней свободы (см.табл.
4.1).Во многих пособиях по статистике
можно найти и готовые таблицы критических
значе­ний коэффициента корреляцииrдля данного уровня значимостиа. В
этом слу­чае отпадает необходимость
в каких-либо вычислениях t:достаточно сравнить полученную величину
коэффициента корреляции с табличным
значением¹⁹. (Например, величина
коэффициента корреляцииr
= 0,55будет существенной на уровне
значимостир = 0,01даже
для выборки объемом 105,так как крити­ческое значение составляет
0,254.)

¹⁸См.:Дружинин Н.
К.Логика оценки статистических
гипотез. М.: Статистика, 1973.
С. 112—114.

“См., в частности:Ликеш
И., Ляга И.Основные таблицы математической
статистики. М.: Финансы и статистика,
1985.(Табл. 14.)

Множественная регрессия
и путевой анализ

Выше описывалась модель линейной
регрессии для двух переменных. В
дей­ствительности социолог довольно
редко сталкивается со столь простыми
моде­лями данных. Влияние одного
фактора обычно может объяснить лишь
часть разброса наблюдаемых значений
независимой переменной. Метод частной
кор­реляции позволяет нам проконтролировать
эффекты воздействия любых дру­гих
контрольных переменных, которые мы в
состоянии измерить. (Стоит снова
подчеркнуть здесь, что статистические
методы изучения причинных взаимо­связей,
в отличие от экспериментальных,позволяют нам контролировать лишь те
источники вариации, которые мы способны
концептуализировать и измерить.) Однако
еще более интересной задачей является
контроль одновременного воз­действиянесколькихнезависимых наоднузависимую переменную, а также срав­нение
эффекта воздействия разных независимых
переменных и предсказание «отклика»
независимой переменной. Именно эти
задачи решают методы анали­за, о
которых пойдет речь в данном параграфе.
Наше изложение будет непол­ным, так
как более детальное обсуждение требует
дополнительной математи­ческой
подготовки. Мы будем ориентироваться
на сравнительно скромные цели понимания
общей логики и интерпретации результатов
соответствующих ста­тистических
процедур.

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

Множественная регрессия во многом
аналогична простой (бивариантной)
рег­рессии. Отличие состоит в том, что
регрессия осуществляется по двум и
более независимым переменным одновременно,причем каждая из них входит в рег­рессионное
уравнение с коэффициентом, позволяющим
предсказать значения зависимой переменной
с минимальным количеством ошибок
(критерием здесь снова является метод
наименьших квадратов). Частные коэффициенты
в урав­нении множественной регрессии
показывают, какой будет величина
воздействия соответствующей независимой
переменной на зависимуюпри контролевлия­ния других независимых переменных.
Если воспользоваться простейшей
сис­темой обозначений, то уравнение
множественной регрессии для трех
независи­мых переменных можно записать
как:

где Y—это предсказываемое
значение зависимой переменной,X₁
… Х₃,—неза­висимые переменные, аb,
…b₃,
—частные коэффициенты регрессии
для каж­дой из зависимых переменных.

Коэффициенты bмогут быть
интерпретированы как показатели влияния
каж­дой из независимых переменных на
зависимую при контроле всех других
неза­висимых переменных в уравнении.
В отличие от коэффициентов частной
кор­реляции коэффициенты регрессии
обладают размерностью. Они показывают,на
сколько единиц изменится зависимая
переменная при увеличении независи­мой
на одну единицу (при контроле всех
остальных переменных модели). Пусть,
например, мы построили уравнение
множественной регрессии, описываю­щее
зависимость дохода от интеллекта(Х₁)и стажа работы(Х₂).Если вели­чинаb₁оказалась
равной 100,это означает,
что каждый дополнительный балл по шкале
интеллекта увеличивает доход на
100рублей. Значениеb₂=
950 говорит нам, что год стажа прибавляет
950рублей. Однако «сырые» оценки
интеллекта и стажа измерены в разных
единицах. Для определениясравни­тельной
значимостинезависимых переменных,
входящих в уравнение мно­жественной
регрессии, мы должны подвергнуть все
переменныестандар­тизации(т. е.
перевести их в Z-оценки,см. выше). Стандартизованные ко­эффициенты
множественной регрессии, которые удобнее
всего обозначать какb*(либо греч.
«бета» — B),меняются в пределах от -1,0до +1,0.Они сохраняют свою
величину при изменении масштаба шкалы:
переход от измерения возраста в годах
к измерению в днях не изменит
соответству­ющийb*.

Стандартизованные коэффициенты позволяют
оценить «вклад» каждой из
переменных-предикторов в предсказание
значений независимой перемен­ной.
Если в примере с влиянием интеллекта и
стажа работы на доход ока­жется, что
b₁*= 0,25,
a b₂*=0,30,то можно заключить, что сравнительная
значимость «веса» интеллекта и стажа
в предсказании дохода различаются
незначительно. Если же для одной
переменнойb₁*
=0,80,тогда какb₂*
=0,40, мы можем сказать, что эффект
воздействия второй переменной в два
раза меньше эффекта первой.

Чтобы определить ожидаемые значения
зависимой переменной для отдельных
индивидов, достаточно подставить в
уравнение множественной регрессии
со­ответствующие значения
переменных-предикторов и вычисленных
коэффици­ентов Ь.Пусть, например,
мы хотим рассчитать прогнозное значение
величины дохода для человека, чей
коэффициент интеллекта равен
110,а стаж работы — 20годам. Еслиb₁,
как в вышеприведенном примере,
составляет 100,
b₂
= 950, а слагаемое а
= 50000,то мы получим:

ожидаемый доход = 50000 +100х 110 + 950х 20 =
80000руб.

Множественную регрессию можно использовать
и для предсказания средних групповых
значений,например среднего дохода
мужчин-врачей. Единственное различие
в данном случае заключается в использованиисредних значений неза­висимых
переменныхдля подстановки в уравнение
множественной регрессии. В качестве
независимой переменной множественной
регрессии могут исполь­зоваться и
дихотомические переменные, которым
приписывают значения 0и
1 (например, пол). Для того чтобы
включить в уравнение номинальную
перемен­ную с более чем двумя
категориями, нужно создать соответствующее
число новых,«фиктивных» переменных,каждая из которых будет кодироваться
как О или 1в зависимости
от наличия или отсутствия категории-признака.
Скажем, состоящую из трех категорий
переменную «цвет глаз» можно представить
с помощью трех переменных:Х₁
—«голубые глаза»,Х₂
—«карие глаза»,Х₃,
— «зеленые глаза». (Человек с
голубыми глазами получит 1поХ₁и 0по
двум другим переменным.)

Метод множественной регрессии очень
популярен среди социологов. Вот,
на­пример, как выглядели результаты
его применения в исследовании Л. Бэрона
и М. Строса, изучавших факторы, влияющие
на статистику изнасилований²⁰.
Использованная в планировании этого
исследования матрица данных включа­ла
в себя в качестве объектов («случаев»)
различные штаты США. Признаками, по
которым описывались штаты, служили
около десятка независимых и соб­ственно
контрольных переменных, предположительно
воздействующих на за­висимую
переменную,—количество
зарегистрированных полицией изнасило­ваний
на 100000населения в год
для данного штата (по данным ежегодных
статистических отчетов ФБР). Предполагалось,
что существующие различия между штатами
в уровне изнасилований можно будет
объяснить различиями в уровнях независимых
переменных. Нужно отметить, что разброс
«случаев» по зависимой переменной был
весьма велик —от
71,9на Аляске до 8,2в
Север­ной Дакоте (1979).Из десятка переменных, включенных в
уравнение множе­ственной регрессии,
девять оказались статистически значимы.
Основные ре­зультаты регрессионного
анализа для семи переменных представлены
втаб­лице 8.12.

Из таблицы видно, что индекс совокупного
тиража порнографических журна­лов
(интегральный показатель, учитывающий
уровни продаж восьми популяр­ных
изданий) имеет коэффициент регрессии
6,99.Это означает, что рост индек­са
на единицу в среднем увеличивает
количество изнасилований почти на7случаев (в расчете на 100000населения). Весьма значительно и влияние
чис­ла убийств, что особенно заметно
при сравнении стандартизованных
коэффи­циентов (b*), не зависящих
от шкалы измерения признака. Фактически
количе­ство убийств вносит самый
значительный «вклад» в предсказание
значений за­висимой переменной (b*
= 0,55).Интересно отметить, что одна
из независимых переменных в описываемом
исследовании —индекс
положения женщин, рас­считанный на
основании 22-х политических, экономических
и социальных ин­дикаторов,—при
анализе простых взаимосвязей
продемонстрировала практи­чески
нулевую корреляцию с количеством
изнасилований (г = 0,17),причем

²⁰ Baron
L., Strauss M. A. Sexual
Stratification, Pornography, and Rape in the United States //
Malamuth
N.. Donnerstein E. (eds.)
Pornography and Sexual Aggression. Orlando etal.:

Academic Press, 1984.
P. 185—209.

²¹Таблица приводится
в сокращении по источнику:
Baron L, Strauss V. A.
Sexual Stratification, Pornography, and Rape…

результаты анализа диаграммы рассеивания
также не дали никаких свидетельств в
пользу гипотезы о нелинейной связи.
Множественная регрессия позволила
уточнить первоначальные выводы: при
контроле прочих переменных модели, чем
вышестатус женщин, темвышеуровень изнасилований (результат,
которо­му довольно трудно найти
теоретическое объяснение). Использование
девяти независимых переменных позволило
объяснить 83%дисперсии в
показателях количества изнасилований
(квадрат коэффициента множественной
корреляцииr²составил
0,83).

При интерпретации результатов
множественной регрессии стандартизован­ные
коэффициенты, как уже говорилось,
используют в качестве показателей
значимости, «вклада» соответствующих
переменных. Эта трактовка верна лишь в
определенных пределах. При нарушении
некоторых условий сравне­ние абсолютных
величин стандартизованных коэффициентов
может вести к неверным выводам. Дело в
том, что коэффициенты регрессии подвержены
влиянию случайных ошибок измерения(см. с. 116).Использованиененадеж­ных индикаторов«сдвигает»
регрессионные коэффициенты к нулю²².
Ины­ми. словами,более надежныеиндикаторы даютболее высокиеоценки
коэф­фициентов. Пусть, например, для
предсказания риска сердечно-сосудистых
заболеваний использовались две
независимые переменные индивидуально­го
уровня—«ориентация на
достижения» и «склонность подавлять
агрес­сию»,—причем
шкала для измерения первой обладалаболее высокимкоэф­фициентом
надежности. Если стандартизованный
коэффициент регрессии для достиженческой
мотивации окажетсявыше,чем для
подавления агрес­сии, это может
рассматриваться как следствие таких
содержательных раз­личий между
переменными, которые важны с точки
зрения теории психосо­циальных
факторов заболеваемости. Но нельзя
исключить и альтернатив­ное объяснение,
связывающее более высокий регрессионный
коэффициент первой переменной с побочными
эффектами методов измерения: влияние
ориентации на достижения не превосходит
влияния, оказываемого на риск инфаркта
склонностью подавлять агрессию, а
наблюдаемые различия рег­рессионных
коэффициентов связаны лишь с ненадежностью
использован­ных индикаторов склонности
к подавлению.

Другая проблема, требующая некоторой
осторожности в интерпретации ко­эффициентов
регрессии, возникает вследствие того,
что модель множествен­ной регрессии
не обязывает нас ни к каким строгим
предположениям о при­чинных связях
между независимыми переменными.Регрессионное уравне­ние, образно
говоря, не делает никаких различий между
собственно независимыми, т. е. теоретически
специфицированными, переменными и
дополнительными—контрольными,
опосредующими и т. п.—факторами,
вводимыми в модель с целью уточнения.
В тех случаях, когда теоретическая
гипотеза, проверяемая в ходе исследования,
допускает: 1)существование
взаимосвязей между независимыми
переменными, 2)наличие
прямых и кос­венных (опосредованных)
влияний, а также 3)использование нескольких индикаторов
для каждого латентного фактора, могут
понадобиться более совершенные
статистические методы. Одна из возможностей
здесь—это использованиепутевого анализа.

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

Путевой анализ —
один из основных способов построения
и проверки причин­ных моделей в
социологии.Многие более продвинутые
статистические техники основаны на
сходной исследовательской методологии.

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

Путевой анализ включает в себя технику
представления прямых и косвен­ных
причинных влияний при помощи специальных
диаграмм (потоковых графов). Эти диаграммы
часто называют просто причинными
(структурны­ми) моделями.

Последовательно «считывая» такую
модель, можно легко определить все пути
влияния одной переменной на другую и
соответственно оценить величину чис­того
эффекта. Во многих разделах этой книги
причинные модели уже исполь­зовались
для представления сравнительно сложных
причинных гипотез (см., например, с.
122),поэтому общая логика их построения
не требует детального обсуждения.
Порядок представления переменных на
диаграмме отражает пред­полагаемое
направление причинной связи, а диапазон
включенных в диаграм­му переменных
и отношения между ними зависят от
принятых исследователем теоретических
гипотез. Так называемыепутевые
коэффициенты,описыва­ющие связи
между переменными (связям соответствуют
стрелочки на диаг­рамме),равны
стандартизованным коэффициентам
множественной рег­рессии (b*)²³.

Обычно путевую диаграмму рисуют слева
направо —от самых «ранних»
по порядку следования независимых
переменных до зависимой. Путевые
коэффи­циенты часто обозначают
латинскими«р»с подстрочными
индексами(р₂₁
— это путевой коэффициент для
связи между переменнымиX₁
—>Х₂).Нарисунке
24в качестве примера
изображена путевая диаграмма, отражающая
гипотети­ческие отношения между
интеллектом(Х₁),образованием(Х₂),социально-эко­номическим статусом(Х₃),доходом(Х₄)и размерами сбережений(Х₅).

Специальные правила позволяют перевести
отношения, изображенные на ди­аграмме,
в совокупность структурных уравнений,описывающих механизмы прямого и
опосредованного воздействия одних
переменных на другие. Нари­сунке
24,в частности, видно, что не
существует пути дляпрямого воздействия
интеллекта на размеры сбережений,
однако общий эффект воздействия
интел­лекта будет включать в себясовокупность непрямых эффектов:
X₁воздейству­ет наX₅ и
через образование(Х₂),и через достигнутый статус(Х₃),и через доход(Х₄).Иными
словами, хотя и нельзя утверждать, что
склонность откладывать деньги «в
кубышку» зависит от умственных
способностей, последние влияют

²³В оценивании также
используется метод наименьших квадратов.

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

В общем случае, полный эффект влияния
переменной равен сумме ее непосред­ственного
эффекта и всех косвенных эффектов
влияния.Величины возмуще­ний (е₂
— e₄)
Haрисунке позволяют оценить, насколько
хорошо работает мо­дель, показывая,
какая часть дисперсии соответствующей
переменной осталась необъясненной. В
результате путевой анализ позволяет
пересматривать и уточ­нять исходную
теоретическую модель, сравнивать
«эффективность» несколь­ких
конкурирующих теорий для объяснения
существующей совокупности эм­пирических
наблюдений. Существуют даже компьютерные
программы, осу­ществляющие автоматический
поиск наилучшей структурной модели, т.
е. процедуру, сходную с отбором из
нескольких существующих теорий та­кой,
которая максимально соответствовала
бы полученным в исследовании дан­ным²⁴.
Важно, однако, осознавать, что сами по
себе результаты применения регрессионных
методов и причинных моделей (регрессионные
коэфициенты, линии регрессии, путевые
диаграммы) решают прежде всего задачуобобщен­ного описанияуже полученных
эмпирических данных. Они могут служить
на­дежной основой для интерполяции,
оценки положения гипотетических «точек»
в пределах ряда наблюдавшихся значений,
однако их использование в целяхэкстраполяции и прогнозаможет
вести к существенным ошибкам в тех
случа­ях, когда такой прогноз не
подкреплен более широкой теорией, не
сводимой к отдельной модели для конечной
совокупности данных. (Достаточно указать
в качестве примера на многочисленные
ошибочные прогнозы в экономике
—на уке, где количество эмпирических
данных и описывающих их структурных
моделей многократно превзошло количество
существующих теорий.)

²⁴Подробнее см.:И.
Ф. Девятко.Диагностическая процедура
в социологии: очерк ис­тории и теории.
М.: Наука, 1993.С.121—136.

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

Дополнительная литература

Вайнберг Дж., Шумекер Дж..Статистика. М.: Финансы и статис­тика,
1979.

Гласе Дж., Стэнли Дж.Статистические
методы в педагогике и пси­хологии.
М.: Прогресс, 1976.

Интерпретация и анализ данныхв
социологическом исследовании. М.: Наука,
1987.

Татарова Г. Г.Типологический анализ
в социологии. М.: Наука,

1993.

Типология и классификацияв
социологических исследованиях. М.:

Наука, 1982.

Толстова Ю. Н.Логика математического
анализа социологических

данных. М.: Наука, 1991.

Хейс Д.Причинный анализ в статистических
исследованиях. М.:

Финансы и статистика, 1981.

Флейс Дж.Статистические методы для
изучения таблиц долей и пропорций. М.:
Финансы и статистика, 1989.

Ядов В. А.Социологическое исследование:
методология, програм­ма, методы. 2-е
изд. М.: Наука, 1987.Гл.
5.

Исследуем отношение между переменными¶

В этой главе исследуются отношения между переменными.

• Мы будем визуализировать отношения с помощью диаграмм рассеяния (scatter plots), диаграмм размаха (box plots) и скрипичных диаграмм (violin plots),

• И мы будем количественно определять отношения, используя корреляцию (correlation) и простую регрессию (simple regression).

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

In [1]:

from os.path import basename, exists

def download(url):
    filename = basename(url)
    if not exists(filename):
        from urllib.request import urlretrieve
        local, _ = urlretrieve(url, filename)
        print('Downloaded ' + local)

download('https://github.com/AllenDowney/' +
         'ElementsOfDataScience/raw/master/brfss.hdf5')

Изучение отношений¶

В качестве первого примера мы рассмотрим взаимосвязь между ростом и весом.

Мы будем использовать данные из Системы наблюдения за поведенческими факторами риска (BRFSS), которая находится в ведении Центров по контролю за заболеваниями по адресу https://www.cdc.gov/brfss.

В опросе приняли участие более 400 000 респондентов, но, чтобы произвести анализ, я выбрал случайную подвыборку из 100 000 человек.

In [2]:

import pandas as pd

brfss = pd.read_hdf('brfss.hdf5', 'brfss')
brfss.shape

Out[3]:

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

Мы начнем с HTM4, который записывает рост каждого респондента в см, и WTKG3, который записывает вес в кг.

In [4]:

height = brfss['HTM4']
weight = brfss['WTKG3']

Чтобы визуализировать взаимосвязь между этими переменными, мы построим диаграмму рассеяния (scatter plot).

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

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

см. документацию по plot

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

Первый из них – перекрытие (overplotted), то есть точки данных накладываются друг на друга, поэтому вы не можете сказать, где много точек, а где только одна.

Когда это происходит, результаты могут вводить в заблуждение.

Один из способов улучшить график – использовать прозрачность (transparency), что мы можем сделать с помощью ключевого аргумента alpha. Чем ниже значение alpha, тем прозрачнее каждая точка данных.

Вот как это выглядит с alpha=0.02.

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

При markersize=1 и низком значении alpha диаграмма рассеяния будет менее насыщенной.

Вот как это выглядит.

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

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

Такое добавление случайного шума называется дрожанием (jittering).

Дрожание – это добавление случайного шума к данным для предотвращения перекрытия статистических графиков. Если непрерывное измерение округлено до некоторой удобной единицы, может произойти перекрытие. Это приводит к превращению непрерывной переменной в дискретную порядковую переменную. Например, возраст измеряется в годах, а масса тела – в фунтах или килограммах. Если вы построите диаграмму разброса веса в зависимости от возраста для достаточно большой выборки людей, там может быть много людей, записанных, скажем, с 29 годами и 70 кг, и, следовательно, в этой точке будет нанесено много маркеров (29, 70).

Чтобы уменьшить перекрытие, вы можете добавить к данным небольшой случайный шум. Размер шума часто выбирается равным ширине единицы измерения. Например, к значению 70 кг вы можете добавить количество u , где u – равномерная случайная величина в интервале [-0,5, 0,5]. Вы можете обосновать дрожание, предположив, что истинный вес человека весом 70 кг с равной вероятностью находится в любом месте интервала [69,5, 70,5].

Контекст данных важен при принятии решения о дрожании. Например, возраст обычно округляется в меньшую сторону: 29-летний человек может праздновать свой 29-й день рождения сегодня или, возможно, ему исполнится 30 завтра, но ей все равно 29 лет. Следовательно, вы можете изменить возраст, добавив величину v , где v – равномерная случайная величина в интервале [0,1]. (Мы игнорируем статистически значимый случай женщин, которым остается 29 лет в течение многих лет!)

Источник: Jittering to prevent overplotting in statistical graphics

Мы можем использовать NumPy для добавления шума из нормального распределения со средним 0 и стандартным отклонением 2.

см. документацию NumPy

In [8]:

import numpy as np

noise = np.random.normal(0, 2, size=len(brfss))
height_jitter = height + noise

Вот как выглядит график с дрожащими (jittered) высотами.

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

In [10]:

noise = np.random.normal(0, 2, size=len(brfss))
weight_jitter = weight + noise

Наконец, давайте увеличим масштаб области, где находится большинство точек данных.

Функции xlim и ylim устанавливают нижнюю и верхнюю границы для осей $x$ и $y$; в данном случае мы наносим рост от 140 до 200 сантиметров и вес до 160 килограмм.

Вот как это выглядит.

Теперь у нас есть достоверная картина взаимосвязи между ростом и весом.

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

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

Упражнение: Набирают ли люди вес с возрастом? Мы можем ответить на этот вопрос, визуализировав взаимосвязь между весом и возрастом.

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

Набор данных BRFSS включает столбец AGE, который представляет возраст каждого респондента в годах. Чтобы защитить конфиденциальность респондентов, возраст округляется до пятилетних интервалов. AGE содержит середину интервалов (bins).

• Извлеките переменную 'AGE' из фрейма данных brfss и присвойте ее age.

• Постройте функцию вероятности (Probability mass function, PMF) для age в виде гистограммы, используя Pmf из empiricaldist.

empiricaldist – библиотека Python, представляющая эмпирические функции распределения.

In [14]:

try:
    import empiricaldist
except ImportError:
    !pip install empiricaldist

In [15]:

from empiricaldist import Pmf

Упражнение: Теперь давайте посмотрим на распределение веса.

Столбец, содержащий вес в килограммах, – это WTKG3. Поскольку этот столбец содержит много уникальных значений, отображение его как функции вероятности (PMF) работает плохо.

Чтобы получить лучшее представление об этом распределении, попробуйте построить график функции распределения (Cumulative distribution function, CDF).

Вычислите функцию распределения (CDF) нормального распределения с тем же средним значением и стандартным отклонением и сравните его с распределением веса.

Подходит ли нормальное распределение для этих данных? А как насчет логарифмического преобразования весов?

Упражнение: Теперь давайте построим диаграмму разброса (scatter plot) для weight и age.

Отрегулируйте alpha и markersize, чтобы избежать наложения (overplotting). Используйте ylim, чтобы ограничить ось y от 0 до 200 килограммов.

Упражнение: В предыдущем упражнении возрасты указаны в столбцах, потому что они были округлены до 5-летних интервалов (bins). Если мы добавим дрожание (jitter), диаграмма рассеяния покажет взаимосвязь более четко.

• Добавьте случайный шум к age со средним значением 0 и стандартным отклонением 2.5.
• Создайте диаграмму рассеяния и снова отрегулируйте alpha и markersize.

Визуализация отношений¶

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

Я начну с диаграммы разброса веса в зависимости от возраста.

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

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

С этой точки зрения кажется, что вес увеличивается до 40-50 лет, а затем начинает уменьшаться.

Если мы пойдем дальше, то сможем использовать ядерную оценку плотности (Kernel Density Estimation, KDE) для оценки функции плотности в каждом столбце и построения графика. И для этого есть название – скрипичная диаграмма (violin plot).

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

Вот так:

In [24]:

data = brfss.dropna(subset=['AGE', 'WTKG3'])
data.shape

dropna() создает новый фрейм данных, который удаляет строки из brfss, где AGE или WTKG3 равны NaN.

Теперь мы можем вызвать функцию violinplot.

см. документацию по violinplot

Аргументы x и y означают, что нам нужно AGE на оси x и WTKG3 на оси y.

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

Аргумент inner=None немного упрощает график.

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

Другой, связанный с этим способ просмотра данных, называется диаграмма размаха (ящик с усами, box plot).

Код для создания диаграммы размаха очень похож.

см. документацию по boxplot

Я включил аргумент whis=10, чтобы отключить функцию, которая нам не нужна.

Каждый прямоугольник представляет распределение веса в возрастной группе. Высота каждого прямоугольника представляет собой диапазон от 25-го до 75-го процентиля. Линия в середине каждого прямоугольника – это медиана. Шипы, торчащие сверху и снизу, показывают минимальное и максимальное значения.

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

• Глядя на медианы, кажется, что люди в возрасте от 40 лет являются самыми тяжелыми; люди младшего и старшего возраста легче.

• Глядя на размеры ящиков, кажется, что люди в возрасте от 40 также имеют наибольший разброс в весе.

• Эти графики также показывают, насколько искажено распределение веса; то есть самые тяжелые люди намного дальше от медианы, чем самые легкие.

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

Мы можем сделать это с помощью Pyplot-функции yscale.

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

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

Упражнение: Ранее мы рассмотрели диаграмму рассеяния (scatter plot) по росту и весу и увидели, что более высокие люди, как правило, тяжелее. Теперь давайте более подробно рассмотрим диаграмму размаха (box plot).

Фрейм данных brfss содержит столбец с именем _HTMG10, который представляет высоту в сантиметрах, разбитую на группы по 10 см.

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

• Постройте ось Y в логарифмическом масштабе.

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

plt.xticks(rotation='45')

Упражнение: В качестве второго примера давайте посмотрим на взаимосвязь между доходом (income) и ростом.

В BRFSS доход представлен как категориальная переменная; то есть респондентов относят к одной из 8 категорий доходов. Имя столбца – INCOME2.

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

• Извлеките INCOME2 из brfss и присвойте его income.

• Постройте функцию вероятности (PMF) для income в виде гистограммы (bar chart).

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

Упражнение: Создайте скрипичную диаграмму (violin plot), которая показывает распределение роста в каждой группе дохода.

Вы видите взаимосвязь между этими переменными?

Корреляция¶

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

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

Для демонстрации я выберу три столбца из набора данных BRFSS:

In [31]:

columns = ['HTM4', 'WTKG3', 'AGE']
subset = brfss[columns]

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

С этим подмножеством данных мы можем использовать метод corr, например:

Out[32]:

Результатом является корреляционная матрица. В первой строке корреляция HTM4 с самим собой равна 1. Это ожидаемо; корреляция чего-либо с самим собой равна 1.

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

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

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

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

Люди за сорок – самые тяжелые; люди младшего и старшего возраста легче. Итак, эта связь нелинейна.

Но корреляция измеряет только линейные отношения. Если связь нелинейная, корреляция обычно недооценивает ее силу.

Чтобы продемонстрировать, я сгенерирую несколько поддельных данных: xs содержит точки с равным интервалом между -1 и 1.

ys – это квадрат xs плюс некоторый случайный шум.

In [34]:

xs = np.linspace(-1, 1)
ys = xs**2 + np.random.normal(0, 0.05, len(xs))

Вот диаграмма рассеяния для xs и ys.

Понятно, что это сильная связь; если вам дано x, вы можете гораздо лучше догадаться о y.

Но вот корреляционная матрица:

Out[36]:

array([[1.        , 0.01135475],
       [0.01135475, 1.        ]])

Несмотря на то, что существует сильная нелинейная зависимость, вычисленная корреляция близка к 0.

В общем, если корреляция высока, то есть близка к 1 или -1, вы можете сделать вывод, что существует сильная линейная зависимость.
Но если корреляция близка к 0, это не означает, что связи нет; может быть связь нелинейная.

Это одна из причин, по которой я считаю, что корреляция не является хорошей статистикой.

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

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

Я создам два поддельных набора данных, чтобы продемонстрировать суть дела. В каждом наборе данных xs представляет возраст, а ys – вес.

Я использую np.random.seed для инициализации генератора случайных чисел, поэтому мы получаем одни и те же результаты при каждом запуске.

А вот и второй набор данных:

Я построил эти примеры так, чтобы они выглядели одинаково, но имели существенно разные корреляции:

In [39]:

rho1 = np.corrcoef(xs1, ys1)[0][1]
rho1

In [40]:

rho2 = np.corrcoef(xs2, ys2)[0][1]
rho2

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

В первом примере средняя прибавка в весе за 30 лет составляет менее 1 килограмма; во втором больше 5 килограммов!

Если нас беспокоит влияние увеличения веса на здоровье, второе соотношение, вероятно, более важно, даже если корреляция ниже.

Статистика, которая нас действительно волнует, – это наклон линии, а не коэффициент корреляции.

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

Упражнения: Цель BRFSS – изучить факторы риска для здоровья, поэтому в него включены вопросы о диете.

Столбец _VEGESU1 представляет количество порций овощей, которые респонденты ели в день.

Посмотрим, как эта переменная связана с возрастом и доходом.

• Во фрейме данных brfss выберите столбцы 'AGE', INCOME2 и _VEGESU1.
• Вычислите корреляционную матрицу для этих переменных.

Упражнение: В предыдущем упражнении корреляция между доходом и потреблением овощей составляет около 0.12. Корреляция между возрастом и потреблением овощей составляет примерно -0.01.

Что из следующего является правильной интерпретацией этих результатов?

• A: люди в этом наборе данных с более высоким доходом едят больше овощей.
• B: Связь между доходом и потреблением овощей линейна.
• C: Пожилые люди едят больше овощей.
• D: Между возрастом и потреблением овощей может быть сильная нелинейная зависимость.

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

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

Простая регрессия¶

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

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

Тот, что слева, имеет более высокую корреляцию, около 0,75 по сравнению с 0,5.

Но в этом контексте статистика, которая нас, вероятно, волнует, – это наклон линии, а не коэффициент корреляции.

Чтобы оценить наклон, мы можем использовать linregress из SciPy-библиотеки stats.

см. документацию по scipy.stats.linregress

Результатом является объект LinregressResult, содержащий пять значений: slope – наклон линии, наиболее подходящей для данных; intercept – это пересечение линии регрессии.

Для фальшивого набора данных 1 расчетный наклон составляет около 0,019 кг в год или около 0,56 кг за 30-летний период.

Вот результаты для фальшивого набора данных 2.

Расчетный наклон почти в 10 раз выше: около 0,18 килограмма в год или около 5,3 килограмма за 30 лет:

То, что здесь называется rvalue, – это корреляция, которая подтверждает то, что мы видели раньше; первый пример имеет более высокую корреляцию, около 0,75 по сравнению с 0,5.

Но сила эффекта, измеренная по наклону линии, во втором примере примерно в 10 раз выше.

Мы можем использовать результаты linregress для вычисления линии тренда: сначала мы получаем минимум и максимум наблюдаемых xs; затем мы умножаем на наклон и добавляем точку пересечения.

Вот как это выглядит для первого примера.

То же самое и со вторым примером.

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

Рост и вес¶

Теперь рассмотрим пример с реальными данными.
Вот еще раз диаграмма рассеяния для роста и веса.

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

In [52]:

subset = brfss.dropna(subset=['WTKG3', 'HTM4'])
height_clean = subset['HTM4']
weight_clean = subset['WTKG3']

Теперь мы можем вычислить линейную регрессию.

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

Как и раньше, мы можем вычислить линию тренда:

In [54]:

fx = np.array([height_clean.min(), height_clean.max()])
fy = res_hw.intercept + res_hw.slope * fx

Наклон этой линии соответствует диаграмме рассеяния.

Линейная регрессия имеет ту же проблему, что и корреляция; она только измеряет силу линейной связи.

Вот диаграмма рассеяния веса по сравнению с возрастом, которую мы видели ранее.

Люди в возрасте от 40 – самые тяжелые; люди младшего и старшего возраста легче. Так что отношения нелинейные.

Если мы не посмотрим на диаграмму рассеяния и вслепую вычислим линию регрессии, мы получим вот что.

Расчетный уклон составляет всего 0,02 килограмма в год или 0,6 килограмма за 30 лет.

А вот как выглядит линия тренда.

Прямая линия плохо отражает взаимосвязь между этими переменными.

Давайте попрактикуемся в простой регрессии.

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

Как мы видели ранее, столбец INCOME2 представляет уровень дохода, а _VEGESU1 представляет количество порций овощей, которые респонденты ели в день.

Постройте диаграмму рассеяния порций овощей в зависимости от дохода, то есть с порциями овощей по оси y и группой доходов по оси x.

Вы можете использовать ylim для увеличения нижней половины оси y.

Упражнение: Теперь давайте оценим наклон зависимости между потреблением овощей и доходом.

• Используйте dropna для выбора строк, в которых INCOME2 и _VEGESU1 не равны NaN.

• Извлеките INCOME2 и _VEGESU1 и вычислите простую линейную регрессию этих переменных.

Каков наклон линии регрессии? Что означает этот наклон в контексте изучаемого нами вопроса?

Упражнение: Наконец, постройте линию регрессии поверх диаграммы рассеяния.

Постановка задачи

При статистическом анализе зависимостей между количественными переменными возникают ситуации, когда представляет интерес расчет и анализ такого показателя как корреляционное отношение (η).

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

В данном разборе рассмотрим способы расчета и анализа η средствами Python.

Не будем углубляться в теорию (про η достаточно подробно написано, например, в [1, с.73], [2, с.412], [3, с.609]), но вспомним основные свойства η:

1. η характеризует степень тесноты любой корреляционной связи (как линейной, так и нелинейной), в отличие от коэффициента корреляции Пирсона r, который характеризует тесноту только линейной связи. Условие r=0 означает отсутствие линейной корреляционной связи между величинами, но при этом между ними может существовать нелинейная корреляционная связь (η>0).

2. η принимает значения от 0 до 1; при η=0 корреляционная связь отсутствует, при η=1 связь считается функциональной; степень тесноты связи можно оценивать по различным общепринятым шкалам, например, по шкале Чеддока и др.

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

4. η не может быть меньше абсолютной величины r: η ≥ |r|.

5. η несимметрично по отношению к исследуемым переменным, то есть η_XY ≠ η_YX.

6. Для расчета η необходимо иметь эмпирические данные эксперимента с повторностями; если же мы имеем просто два набора значений переменных X и Y, то данные нужно группировать. Этот вывод, в общем-то, очевиден – если предпринять попытку рассчитать η по негруппированным данным, получим результат η=1.

Группировка данных для расчета η заключается в разбиении области значений переменных X и Y на интервалы, подсчет частот попадания данных в интервалы и формирование корреляционной таблицы.

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

• когда нарушается условие η ≥ |r|;

• когда r окажется значим, а η нет;

• когда нижняя граница доверительного интервала для η окажется меньше 0 или верхняя граница – больше 1.

Это нужно учитывать при выполнении анализа.

Итак, перейдем к расчетам.

Формирование исходных данных

В качестве исходных данных рассмотрим зависимость расхода среднемесячного расхода топлива автомобиля (л/100 км) (FuelFlow) от среднемесячного пробега (км) (Mileage).

Загрузим исходные данные из csv-файла (исходные данные доступны в моем репозитории на GitHub):

fuel_df = pd.read_csv(
    filepath_or_buffer='data/fuel_df.csv',
    sep=';',
    index_col='Number')
dataset_df = fuel_df.copy()    # создаем копию исходной таблицы для работы
display(dataset_df.head())

Загруженный DataFrame содержит следующие столбцы:

• Month — месяц (в формате Excel)

• Mileage – месячный пробег (км)

• Temperature – среднемесячная температура (°C)

• FuelFlow – среднемесячный расход топлива (л/100 км)

Сохраним нужные нам переменные Mileage и FuelFlow в виде numpy.ndarray.

X = np.array(dataset_df['Mileage'])
Y = np.array(dataset_df['FuelFlow'])

Для удобства дальнейшей работы сформируем сформируем отдельный DataFrame из двух переменных – X и Y:

data_XY_df = pd.DataFrame({
    'X': X,
    'Y': Y})

Настройка заголовков отчета (для дальнейшего формирования графиков):

# Общий заголовок проекта
Task_Project = "Расчет и анализ корреляционного отношения средствами Python"
# Заголовок, фиксирующий момент времени
AsOfTheDate = ""
# Заголовок раздела проекта
Task_Theme = "Анализ расхода топлива автомобиля"
# Общий заголовок проекта для графиков
Title_String = f"{Task_Project}n{AsOfTheDate}"
# Наименования переменных
Variable_Name_X = "Среднемесячный пробег (км)"
Variable_Name_Y = "Среднемесячный расход топлива автомобиля (л/100 км)"

Визуализация и первичная обработка данных

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

mask1 = data_XY_df['X'] > 200
mask2 = data_XY_df['X'] < 2000
data_XY_df = data_XY_df[mask1 & mask2]

X = np.array(data_XY_df['X'])
Y = np.array(data_XY_df['Y'])

Описательная статистика исходных данных:

data_XY_df.describe()

Выполним визуализацию исходных данных:

fig, axes = plt.subplots(figsize=(297/INCH, 210/INCH))
fig.suptitle(Task_Theme)
axes.set_title('Зависимость расхода топлива от пробега')
data_df = data_XY_df
sns.scatterplot(
    data=data_df,
    x='X', y='Y',
    label='эмпирические данные',
    s=50,
    ax=axes)
axes.set_xlabel(Variable_Name_X)
axes.set_ylabel(Variable_Name_Y)
#axes.tick_params(axis="x", labelsize=f_size+4)
#axes.tick_params(axis="y", labelsize=f_size+4)
#axes.legend(prop={'size': f_size+6})
plt.show()
fig.savefig('graph/scatterplot_XY_sns.png', orientation = "portrait", dpi = 300)

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

fig = plt.figure(figsize=(420/INCH, 297/INCH))
ax1 = plt.subplot(2,2,1)
ax2 = plt.subplot(2,2,2)
ax3 = plt.subplot(2,2,3)
ax4 = plt.subplot(2,2,4)
fig.suptitle(Task_Theme)
ax1.set_title('X')
ax2.set_title('Y')

# инициализация данных
data_df = data_XY_df
X_mean = data_df['X'].mean()
X_std = data_df['X'].std(ddof = 1)
Y_mean = data_df['Y'].mean()
Y_std = data_df['Y'].std(ddof = 1)
bins_hist = 'sturges'    # выбор числа интервалов ('auto', 'fd', 'doane', 'scott', 'stone', 'rice', 'sturges', 'sqrt')

# данные для графика плотности распределения X
xmin = np.amin(data_df['X'])
xmax = np.amax(data_df['X'])
nx = 100
hx = (xmax - xmin)/(nx - 1)
x1 = np.linspace(xmin, xmax, nx)
xnorm1 = sps.norm.pdf(x1, X_mean, X_std)
kx = len(np.histogram(X, bins=bins_hist, density=False)[0])
xnorm2 = xnorm1*len(X)*(xmax-xmin)/kx

# данные для графика плотности распределения Y
ymin = np.amin(Y)
ymax = np.amax(Y)
ny = 100
hy = (ymax - ymin)/(ny - 1)
y1 = np.linspace(ymin, ymax, ny)
ynorm1 = sps.norm.pdf(y1, Y_mean, Y_std)
ky = len(np.histogram(Y, bins=bins_hist, density=False)[0])
ynorm2 = ynorm1*len(Y)*(ymax-ymin)/ky

# гистограмма распределения X
ax1.hist(
    data_df['X'],
    bins=bins_hist,
    density=False,
    histtype='bar',    # 'bar', 'barstacked', 'step', 'stepfilled'
    orientation='vertical',   # 'vertical', 'horizontal'
    color = "#1f77b4",
    label='эмпирическая частота')
ax1.plot(
    x1, xnorm2,
    linestyle = "-",
    color = "r",
    linewidth = 2,
    label = 'теоретическая нормальная кривая')
ax1.axvline(X_mean, color='magenta', label = 'среднее значение')
ax1.axvline(np.median(data_df['X']), color='orange', label = 'медиана')
ax1.legend(fontsize = f_size+4)

# гистограмма распределения Y
ax2.hist(
    data_df['Y'],
    bins=bins_hist,
    density=False,
    histtype='bar',    # 'bar', 'barstacked', 'step', 'stepfilled'
    orientation='vertical',   # 'vertical', 'horizontal'
    color = "#1f77b4",
    label='эмпирическая частота')    
ax2.plot(
    y1, ynorm2,
    linestyle = "-",
    color = "r",
    linewidth = 2,
    label = 'теоретическая нормальная кривая')
ax2.axvline(Y_mean, color='magenta', label = 'среднее значение')
ax2.axvline(np.median(data_df['Y']), color='orange', label = 'медиана')
ax2.legend(fontsize = f_size+4)

# коробчатая диаграмма X
sns.boxplot(
    #data=corn_yield_df,
    x=data_df['X'],
    orient='h',
    width=0.3,
    ax=ax3)

# коробчатая диаграмма Y
sns.boxplot(
    #data=corn_yield_df, 
    x=data_df['Y'], 
    orient='h', 
    width=0.3, 
    ax=ax4)

# подписи осей 
ax3.set_xlabel(Variable_Name_X)
ax4.set_xlabel(Variable_Name_Y)

plt.show()
fig.savefig('graph/scatterplot_boxplot_X_Y_sns.png', orientation = "portrait", dpi = 300)

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

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

 Для проверки нормальности распределения воспользуемся критерием Шапиро-Уилка:

# функция для обработки реализации теста Шапиро-Уилка
def Shapiro_Wilk_test(data):
    data = np.array(data)
    result = sci.stats.shapiro(data)
    s_calc = result.statistic    # расчетное значение статистики критерия
    a_calc = result.pvalue    # расчетный уровень значимости
    
    print(f"Расчетный уровень значимости: a_calc = {round(a_calc, DecPlace)}")
    print(f"Заданный уровень значимости: a_level = {round(a_level, DecPlace)}")
    
    if a_calc >= a_level:
        conclusion_ShW_test = f"Так как a_calc = {round(a_calc, DecPlace)} >= a_level = {round(a_level, DecPlace)}" + 
            ", то гипотеза о нормальности распределения по критерию Шапиро-Уилка ПРИНИМАЕТСЯ"
    else:
        conclusion_ShW_test = f"Так как a_calc = {round(a_calc, DecPlace)} < a_level = {round(a_level, DecPlace)}" + 
            ", то гипотеза о нормальности распределения по критерию Шапиро-Уилка ОТВЕРГАЕТСЯ"
    print(conclusion_ShW_test)

# проверка нормальности распределения переменной X
Shapiro_Wilk_test(X)

# проверка нормальности распределения переменной Y
Shapiro_Wilk_test(Y)

Итак, гипотеза о нормальном распределении исходных данных принимается, что позволяет нам в дальнейшем пользоваться статистическим инструментарием для интервального оценивания величины η, проверки гипотез и т.д.

Переходим собственно к расчету корреляционного отношения.

Переходим собственно к расчету корреляционного отношения.

Расчёт и анализ корреляционного отношения

1. Выполним группировку исходных данным по обоим признакам X и Y:

Создадим новую переменную matrix_XY_df для работы с группированными данными:

matrix_XY_df = data_XY_df.copy()

Определим число интервалов группировки (воспользуемся формулой Стерджесса); при этом минимальное число интервалов должно быть не менее 2:

# объем выборки для переменных X и Y
n_X = len(X)
n_Y = len(Y)

# число интервалов группировки
group_int_number = lambda n: round (3.31*log(n_X, 10)+1) if round (3.31*log(n_X, 10)+1) >=2 else 2

K_X = group_int_number(n_X)
K_Y = group_int_number(n_Y)
print(f"Число интервалов группировки для переменной X: {K_X}")
print(f"Число интервалов группировки для переменной Y: {K_Y}")

Выполним группировку данных средствами библиотеки pandas, для этого воспользуемся функцией pandas.cut. В результате получим новые признаки cut_X и cut_X, которые показывают, в какой из интервалов попадает конкретное значение X и Y. Полученные новые признаки добавим в DataFrame matrix_XY_df:

cut_X = pd.cut(X, bins=K_X)
cut_Y = pd.cut(Y, bins=K_Y)

matrix_XY_df['cut_X'] = cut_X
matrix_XY_df['cut_Y'] = cut_Y

display(matrix_XY_df.head())

Теперь мы можем получить корреляционную таблицу с помощью функции pandas.crosstab:

CorrTable_df = pd.crosstab(
    index=matrix_XY_df['cut_X'],
    columns=matrix_XY_df['cut_Y'],
    rownames=['cut_X'],
    colnames=['cut_Y'])

display(CorrTable_df)

# проверка правильности подсчета частот по интервалам
print(f"sum = {np.sum(np.array(CorrTable_df))}")

Функция pandas.crosstab также позволяет формировать более ровные и удобные для восприятия границы интервалов группировки путем задания их вручную. Для расчета η это принципиального значения не имеет, но в отдельных случаях может быть полезно.

Например, зададим вручную границы интервалов группировки для X и Y:

bins_X = pd.IntervalIndex.from_tuples([(200, 400), (400, 600), (600, 800), (800, 1000), (1000, 1200), (1200, 1400), (1400, 1600)])
cut_X = pd.cut(X, bins=bins_X)

bins_Y = pd.IntervalIndex.from_tuples([(6.0, 7.0), (7.0, 8.0), (8.0, 9.0), (9.0, 10.0), (10.0, 11.0), (11.0, 12.0), (12.0, 13.0)])
cut_Y = pd.cut(X, bins=bins_Y)

CorrTable_df2 = pd.crosstab(
    index=pd.cut(X, bins=bins_X),
    columns=pd.cut(Y, bins=bins_Y),
    rownames=['cut_X'],
    colnames=['cut_Y'])
display(CorrTable_df2)

# проверка правильности подсчета частот по интервалам
print(f"sum = {np.sum(np.array(CorrTable_df2))}")

Есть и другой способ получения корреляционной таблицы – с помощью pandas.pivot_table:

matrix_XY_df.pivot_table(
    values=['Y'],
    index='cut_X',
    columns='cut_Y',
    aggfunc=len,
    fill_value=0)

2. Выполним расчет корреляционного отношения:

Для дальнейших расчетов приведем корреляционную таблицу к типу numpy.ndarray:

CorrTable_np = np.array(CorrTable_df)
print(CorrTable_np, type(CorrTable_np))

Итоги корреляционной таблицы по строкам и столбцам:

# итоги по строкам
n_group_X = [np.sum(CorrTable_np[i]) for i in range(K_X)]
print(f"n_group_X = {n_group_X}")

# итоги по столбцам
n_group_Y = [np.sum(CorrTable_np[:,j]) for j in range(K_Y)]
print(f"n_group_Y = {n_group_Y}")

Также нам необходимо получить среднегрупповые значения X и Y для каждой группы (интервала). При этом нужно помнить, что функция pandas.crosstab при группировании расширяет крайние диапазоны на 0.1% с каждой стороны, чтобы включить минимальное и максимальное значения.

Для доступа к данным – границам интервалов, полученным с помощью pandas.cut – существуют методы right и left:

# Среднегрупповые значения переменной X
Xboun_mean = [(CorrTable_df.index[i].left + CorrTable_df.index[i].right)/2 for i in range(K_X)]
Xboun_mean[0] = (np.min(X) + CorrTable_df.index[0].right)/2    # исправляем значения в крайних интервалах
Xboun_mean[K_X-1] = (CorrTable_df.index[K_X-1].left + np.max(X))/2
print(f"Xboun_mean = {Xboun_mean}")

# Среднегрупповые значения переменной Y
Yboun_mean = [(CorrTable_df.columns[j].left + CorrTable_df.columns[j].right)/2 for j in range(K_Y)]
Yboun_mean[0] = (np.min(Y) + CorrTable_df.columns[0].right)/2    # исправляем значения в крайних интервалах
Yboun_mean[K_Y-1] = (CorrTable_df.columns[K_Y-1].left + np.max(Y))/2
print(f"Yboun_mean = {Yboun_mean}", 'n')

Находим средневзевешенные значения X и Y для каждой группы:

Xmean_group = [np.sum(CorrTable_np[:,j] * Xboun_mean) / n_group_Y[j] for j in range(K_Y)]
print(f"Xmean_group = {Xmean_group}")

Ymean_group = [np.sum(CorrTable_np[i] * Yboun_mean) / n_group_X[i] for i in range(K_X)]
print(f"Ymean_group = {Ymean_group}")

Общая дисперсия X и Y:

Sum2_total_X = np.sum(n_group_X * (Xboun_mean - np.mean(X))**2)
print(f"Sum2_total_X = {Sum2_total_X}")

Sum2_total_Y = np.sum(n_group_Y * (Yboun_mean - np.mean(Y))**2)
print(f"Sum2_total_Y = {Sum2_total_Y}")

Межгрупповая дисперсия X и Y (дисперсия групповых средних):

Sum2_between_group_X = np.sum(n_group_Y * (Xmean_group - np.mean(X))**2)
print(f"Sum2_between_group_X = {Sum2_between_group_X}")

Sum2_between_group_Y = np.sum(n_group_X * (Ymean_group - np.mean(Y))**2)
print(f"Sum2_between_group_Y = {Sum2_between_group_Y}")

Внутригрупповая дисперсия X и Y (возникает за счет других факторов – не связанных с другой переменной):

print(f"Sum2_within_group_X = {Sum2_total_X - Sum2_between_group_X}")
print(f"Sum2_within_group_Y = {Sum2_total_Y - Sum2_between_group_Y}")

Эмпирическое корреляционное отношение:

corr_ratio_XY = sqrt(Sum2_between_group_Y / Sum2_total_Y)
print(f"corr_ratio_XY = {corr_ratio_XY}")

corr_ratio_YX = sqrt(Sum2_between_group_X / Sum2_total_X)
print(f"corr_ratio_YX = {corr_ratio_YX}")

Итак, мы получили результат – значение корреляционного отношения.

Оценим тесноту корреляционной связи по шкале Чеддока, для удобства создадим пользовательскую функцию:

def Cheddock_scale_check(r, name='r'):
    # задаем шкалу Чеддока
    Cheddock_scale = {
        f'no correlation (|{name}| <= 0.1)':    0.1,
        f'very weak (0.1 < |{name}| <= 0.2)':   0.2,
        f'weak (0.2 < |{name}| <= 0.3)':        0.3,
        f'moderate (0.3 < |{name}| <= 0.5)':    0.5,
        f'perceptible (0.5 < |{name}| <= 0.7)': 0.7,
        f'high (0.7 < |{name}| <= 0.9)':        0.9,
        f'very high (0.9 < |{name}| <= 0.99)':  0.99,
        f'functional (|{name}| > 0.99)':        1.0}
    
    r_scale = list(Cheddock_scale.values())
    for i, elem in enumerate(r_scale):
        if abs(r) <= elem:
            conclusion_Cheddock_scale = list(Cheddock_scale.keys())[i]
            break
    return conclusion_Cheddock_scale

Шкала Чеддока изначально предназначалась для оценки тесноты линейно корреляционной связи (на основе коэффициента корреляции Пирсона r), но мы ее применим и для корреляционного отношения η (не забывая про свойство η ≥ r!). В выводе функции Cheddock_scale_check можно указать символ, обозначающий величину – аргумент name=chr(951) выводит η вместо r.

В современных исследованиях шкала Чеддока теряет популярность, в последние годы все чаще применяется шкала Эванса (в психосоциальных, медико-биологических и др.исследованиях) ( более подробно про шкалы Чеддока, Эванса и др. – см.[4]). Оценим тесноту корреляционной связи по шкале Эванса, для удобства также создадим пользовательскую функцию:

def Evans_scale_check(r, name='r'):
    # задаем шкалу Эванса
    Evans_scale = {
        f'very weak (|{name}| < 0.19)':         0.2,
        f'weak (0.2 < |{name}| <= 0.39)':       0.4,  
        f'moderate (0.4 < |{name}| <= 0.59)':   0.6,
        f'strong (0.6 < |{name}| <= 0.79)':     0.8,
        f'very strong (0.8 < |{name}| <= 1.0)': 1.0}
    
    r_scale = list(Evans_scale.values())
    for i, elem in enumerate(r_scale):
        if abs(r) <= elem:
            conclusion_Evans_scale = list(Evans_scale.keys())[i]
            break
    return conclusion_Evans_scale    

print(f"Оценка тесноты корреляции по шкале Эванса: {Evans_scale_check(corr_ratio_XY, name=chr(951))}")

Итак, степень тесноты корреляционной связи может быть оценена как высокая (по шкале Чеддока), сильная (по шкале Эванса).

3. Проверка значимости корреляционного отношения:

Рассмотрим нулевую гипотезу:

H0: ηXY = 0
H1: ηXY ≠ 0

Для проверки нулевой гипотезы воспользуемся критерием Фишера:

# расчетное значение статистики критерия Фишера
F_corr_ratio_calc = (n_X - K_X)/(K_X - 1) * corr_ratio_XY**2 / (1 - corr_ratio_XY**2)
print(f"Расчетное значение статистики критерия Фишера: F_calc = {round(F_corr_ratio_calc, DecPlace)}")
# табличное значение статистики критерия Фишера
dfn = K_X - 1
dfd = n_X - K_X
F_corr_ratio_table = sci.stats.f.ppf(p_level, dfn, dfd, loc=0, scale=1)
print(f"Табличное значение статистики критерия Фишера: F_table = {round(F_corr_ratio_table, DecPlace)}")
# расчетный уровень значимости
a_corr_ratio_calc = 1 - sci.stats.f.cdf(F_corr_ratio_calc, dfn, dfd, loc=0, scale=1)
print(f"Расчетный уровень значимости: a_calc = {round(a_corr_ratio_calc, DecPlace)}")
print(f"Заданный уровень значимости: a_level = {round(a_level, DecPlace)}")
# вывод
if F_corr_ratio_calc < F_corr_ratio_table:
    conclusion_corr_ratio_sign = f"Так как F_calc = {round(F_corr_ratio_calc, DecPlace)} < F_table = {round(F_corr_ratio_table, DecPlace)}" + 
        ", то гипотеза о равенстве нулю корреляционного отношения ПРИНИМАЕТСЯ, т.е. корреляционная связь НЕЗНАЧИМА"
else:
    conclusion_corr_ratio_sign = f"Так как F_calc = {round(F_corr_ratio_calc, DecPlace)} >= F_table = {round(F_corr_ratio_table, DecPlace)}" + 
        ", то гипотеза о равенстве нулю корреляционного отношения ОТВЕРГАЕТСЯ, т.е. корреляционная связь ЗНАЧИМА"
print(conclusion_corr_ratio_sign)

4. Доверительный интервал для корреляционного отношения:

# число степеней свободы
f1 = round ((K_X - 1 + n_X * corr_ratio_XY**2)**2 / (K_X - 1 + 2 * n_X * corr_ratio_XY**2))
f2 = n_X - K_X
# вспомогательные величины
z1 = (n_X - K_X) / n_X * corr_ratio_XY**2 / (1 - corr_ratio_XY**2) * 1/sci.stats.f.ppf(p_level, f1, f2, loc=0, scale=1) - (K_X - 1)/n_X
z2 = (n_X - K_X) / n_X * corr_ratio_XY**2 / (1 - corr_ratio_XY**2) * 1/sci.stats.f.ppf(1 - p_level, f1, f2, loc=0, scale=1) - (K_X - 1)/n_X
# доверительный интервал
corr_ratio_XY_low = sqrt(z1) if sqrt(z1) >= 0 else 0
corr_ratio_XY_high = sqrt(z2) if sqrt(z2) <= 1 else 1
print(f"{p_level*100}%-ный доверительный интервал для корреляционного отношения: {[round(corr_ratio_XY_low, DecPlace), round(corr_ratio_XY_high, DecPlace)]}")

Важное замечание: при значениях η близких к 0 или 1 левая или правая граница доверительного интервала может выходить за пределы отрезка [0; 1], теряя содержательный смысл (см. [1, с.80]). Причина этого – в аппроксимационном подходе к определению границ доверительного интервала. Подобные нежелательные явления возможны, и их нужно учитывать при выполнении анализа.

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

Оценим величину коэффициента линейной корреляции:

corr_coef = sci.stats.pearsonr(X, Y)[0]
print(f"Коэффициент линейной корреляции: r = {round(corr_coef, DecPlace)}")

print(f"Оценка тесноты линейной корреляции по шкале Чеддока: {Cheddock_scale_check(corr_coef)}")
print(f"Оценка тесноты линейной корреляции по шкале Эванса: {Evans_scale_check(corr_coef)}")

Проверим значимость коэффициента линейной корреляции:

# расчетный уровень значимости
a_corr_coef_calc = sci.stats.pearsonr(X, Y)[1]
print(f"Расчетный уровень значимости коэффициента линейной корреляции: a_calc = {a_corr_coef_calc}")
print(f"Заданный уровень значимости: a_level = {round(a_level, DecPlace)}")
# вывод
if a_corr_coef_calc >= a_level:
    conclusion_corr_coef_sign = f"Так как a_calc = {a_corr_coef_calc} >= a_level = {round(a_level, DecPlace)}" + 
        ", то гипотеза о равенстве нулю коэффициента линейной корреляции ПРИНИМАЕТСЯ, т.е. линейная корреляционная связь НЕЗНАЧИМА"
else:
    conclusion_corr_coef_sign = f"Так как a_calc = {a_corr_coef_calc} < a_level = {round(a_level, DecPlace)}" + 
        ", то гипотеза о равенстве нулю коэффициента линейной корреляции ОТВЕРГАЕТСЯ, т.е. линейная корреляционная связь ЗНАЧИМА"
print(conclusion_corr_coef_sign)

Теперь проверим значимость отличия линейной корреляционной связи от нелинейной. Для этого рассмотрим нулевую гипотезу:

H0: η² - r² = 0
H1: η² - r² ≠ 0

Для проверки нулевой гипотезы воспользуемся критерием Фишера:

print(f"Корреляционное отношение: {chr(951)} = {round(corr_ratio_XY, DecPlace)}")
print(f"Коэффициент линейной корреляции: r = {round(corr_coef, DecPlace)}")
# расчетное значение статистики критерия Фишера
F_line_corr_sign_calc = (n_X - K_X)/(K_X - 2) * (corr_ratio_XY**2 - corr_coef**2) / (1 - corr_ratio_XY**2)
print(f"Расчетное значение статистики критерия Фишера: F_calc = {round(F_line_corr_sign_calc, DecPlace)}")
# табличное значение статистики критерия Фишера
dfn = K_X - 2
dfd = n_X - K_X
F_line_corr_sign_table = sci.stats.f.ppf(p_level, dfn, dfd, loc=0, scale=1)
print(f"Табличное значение статистики критерия Фишера: F_table = {round(F_line_corr_sign_table, DecPlace)}")
# расчетный уровень значимости
a_line_corr_sign_calc = 1 - sci.stats.f.cdf(F_line_corr_sign_calc, dfn, dfd, loc=0, scale=1)
print(f"Расчетный уровень значимости: a_calc = {round(a_line_corr_sign_calc, DecPlace)}")
print(f"Заданный уровень значимости: a_level = {round(a_level, DecPlace)}")
# вывод
if F_line_corr_sign_calc < F_line_corr_sign_table:
    conclusion_line_corr_sign = f"Так как F_calc = {round(F_line_corr_sign_calc, DecPlace)} < F_table = {round(F_line_corr_sign_table, DecPlace)}" + 
        f", то гипотеза о равенстве {chr(951)} и r ПРИНИМАЕТСЯ, т.е. корреляционная связь ЛИНЕЙНАЯ"
else:
    conclusion_line_corr_sign = f"Так как F_calc = {round(F_line_corr_sign_calc, DecPlace)} >= F_table = {round(F_line_corr_sign_table, DecPlace)}" + 
        f", то гипотеза о равенстве {chr(951)} и r ОТВЕРГАЕТСЯ, т.е. корреляционная связь НЕЛИНЕЙНАЯ"
print(conclusion_line_corr_sign)

Создание пользовательской функции для корреляционного анализа

Для практической работы целесообразно все вышеприведенные расчеты реализовать в виде пользовательских функций:

• функция corr_coef_check – для расчета и анализа коэффициента линейной корреляции Пирсона

• функция corr_ratio_check – для расчета и анализа корреляционного отношения

• функция line_corr_sign_check – для проверка значимости линейной корреляционной связи

Данные функции выводят результаты анализа в виде DataFrame, что удобно для визуального восприятия и дальнейшего использования результатов анализа (впрочем, способ вывода – на усмотрение каждого исследователя).

# Функция для расчета и анализа коэффициента линейной корреляции Пирсона
def corr_coef_check(X, Y, p_level=0.95, scale='Cheddok'):
    a_level = 1 - p_level
    X = np.array(X)
    Y = np.array(Y)
    n_X = len(X)
    n_Y = len(Y)
    # оценка коэффициента линейной корреляции средствами scipy
    corr_coef, a_corr_coef_calc = sci.stats.pearsonr(X, Y)
    # несмещенная оценка коэффициента линейной корреляции (при n < 15) (см.Кобзарь, с.607)
    if n_X < 15:
        corr_coef = corr_coef * (1 + (1 - corr_coef**2) / (2*(n_X-3)))
    # проверка гипотезы о значимости коэффициента корреляции
    t_corr_coef_calc = abs(corr_coef) * sqrt(n_X-2) / sqrt(1 - corr_coef**2)
    t_corr_coef_table = sci.stats.t.ppf((1 + p_level)/2 , n_X - 2)
    conclusion_corr_coef_sign = 'significance' if t_corr_coef_calc >= t_corr_coef_table else 'not significance'
    # доверительный интервал коэффициента корреляции
    if t_corr_coef_calc >= t_corr_coef_table:
        z1 = np.arctanh(corr_coef) - sci.stats.norm.ppf((1 + p_level)/2, 0, 1) / sqrt(n_X-3) - corr_coef / (2*(n_X-1))
        z2 = np.arctanh(corr_coef) + sci.stats.norm.ppf((1 + p_level)/2, 0, 1) / sqrt(n_X-3) - corr_coef / (2*(n_X-1))
        corr_coef_conf_int_low = tanh(z1)
        corr_coef_conf_int_high = tanh(z2)
    else:
        corr_coef_conf_int_low = corr_coef_conf_int_high = '-'    
    # оценка тесноты связи
    if scale=='Cheddok':
        conclusion_corr_coef_scale = scale + ': ' + Cheddock_scale_check(corr_coef)
    elif scale=='Evans':
        conclusion_corr_coef_scale = scale + ': ' + Evans_scale_check(corr_coef)
    # формируем результат            
    result = pd.DataFrame({
        'notation': ('r'),
        'coef_value': (corr_coef),
        'coef_value_squared': (corr_coef**2),
        'p_level': (p_level),
        'a_level': (a_level),
        't_calc': (t_corr_coef_calc),
        't_table': (t_corr_coef_table),
        't_calc >= t_table': (t_corr_coef_calc >= t_corr_coef_table),
        'a_calc': (a_corr_coef_calc),
        'a_calc <= a_level': (a_corr_coef_calc <= a_level),
        'significance_check': (conclusion_corr_coef_sign),
        'conf_int_low': (corr_coef_conf_int_low),
        'conf_int_high': (corr_coef_conf_int_high),
        'scale': (conclusion_corr_coef_scale)
        },
        index=['Correlation coef.'])
        
    return result        

# Функция для расчета и анализа корреляционного отношения
def corr_ratio_check(X, Y, p_level=0.95, orientation='XY', scale='Cheddok'):
    a_level = 1 - p_level
    X = np.array(X)
    Y = np.array(Y)
    n_X = len(X)
    n_Y = len(Y)    
    # запишем данные в DataFrame
    matrix_XY_df = pd.DataFrame({
        'X': X,
        'Y': Y})
    # число интервалов группировки
    group_int_number = lambda n: round (3.31*log(n_X, 10)+1) if round (3.31*log(n_X, 10)+1) >=2 else 2
    K_X = group_int_number(n_X)
    K_Y = group_int_number(n_Y)
    # группировка данных и формирование корреляционной таблицы
    cut_X = pd.cut(X, bins=K_X)
    cut_Y = pd.cut(Y, bins=K_Y)
    matrix_XY_df['cut_X'] = cut_X
    matrix_XY_df['cut_Y'] = cut_Y
    CorrTable_df = pd.crosstab(
        index=matrix_XY_df['cut_X'],
        columns=matrix_XY_df['cut_Y'],
        rownames=['cut_X'],
        colnames=['cut_Y'])
    CorrTable_np = np.array(CorrTable_df)
    # итоги корреляционной таблицы по строкам и столбцам
    n_group_X = [np.sum(CorrTable_np[i]) for i in range(K_X)]
    n_group_Y = [np.sum(CorrTable_np[:,j]) for j in range(K_Y)]
    # среднегрупповые значения переменной X
    Xboun_mean = [(CorrTable_df.index[i].left + CorrTable_df.index[i].right)/2 for i in range(K_X)]
    Xboun_mean[0] = (np.min(X) + CorrTable_df.index[0].right)/2    # исправляем значения в крайних интервалах
    Xboun_mean[K_X-1] = (CorrTable_df.index[K_X-1].left + np.max(X))/2
    # среднегрупповые значения переменной Y
    Yboun_mean = [(CorrTable_df.columns[j].left + CorrTable_df.columns[j].right)/2 for j in range(K_Y)]
    Yboun_mean[0] = (np.min(Y) + CorrTable_df.columns[0].right)/2    # исправляем значения в крайних интервалах
    Yboun_mean[K_Y-1] = (CorrTable_df.columns[K_Y-1].left + np.max(Y))/2
    # средневзевешенные значения X и Y для каждой группы
    Xmean_group = [np.sum(CorrTable_np[:,j] * Xboun_mean) / n_group_Y[j] for j in range(K_Y)]
    Ymean_group = [np.sum(CorrTable_np[i] * Yboun_mean) / n_group_X[i] for i in range(K_X)]
    # общая дисперсия X и Y
    Sum2_total_X = np.sum(n_group_X * (Xboun_mean - np.mean(X))**2)
    Sum2_total_Y = np.sum(n_group_Y * (Yboun_mean - np.mean(Y))**2)
    # межгрупповая дисперсия X и Y (дисперсия групповых средних)
    Sum2_between_group_X = np.sum(n_group_Y * (Xmean_group - np.mean(X))**2)
    Sum2_between_group_Y = np.sum(n_group_X * (Ymean_group - np.mean(Y))**2)
    # эмпирическое корреляционное отношение
    corr_ratio_XY = sqrt(Sum2_between_group_Y / Sum2_total_Y)
    corr_ratio_YX = sqrt(Sum2_between_group_X / Sum2_total_X)
    try:
        if orientation!='XY' and orientation!='YX':
            raise ValueError("Error! Incorrect orientation!")
        if orientation=='XY':
            corr_ratio = corr_ratio_XY
        elif orientation=='YX':
            corr_ratio = corr_ratio_YX
    except ValueError as err:
        print(err)
    # проверка гипотезы о значимости корреляционного отношения
    F_corr_ratio_calc = (n_X - K_X)/(K_X - 1) * corr_ratio**2 / (1 - corr_ratio**2)
    dfn = K_X - 1
    dfd = n_X - K_X
    F_corr_ratio_table = sci.stats.f.ppf(p_level, dfn, dfd, loc=0, scale=1)
    a_corr_ratio_calc = 1 - sci.stats.f.cdf(F_corr_ratio_calc, dfn, dfd, loc=0, scale=1)
    conclusion_corr_ratio_sign = 'significance' if F_corr_ratio_calc >= F_corr_ratio_table else 'not significance'
    # доверительный интервал корреляционного отношения
    if F_corr_ratio_calc >= F_corr_ratio_table:
        f1 = round ((K_X - 1 + n_X * corr_ratio**2)**2 / (K_X - 1 + 2 * n_X * corr_ratio**2))
        f2 = n_X - K_X
        z1 = (n_X - K_X) / n_X * corr_ratio**2 / (1 - corr_ratio**2) * 1/sci.stats.f.ppf(p_level, f1, f2, loc=0, scale=1) - (K_X - 1)/n_X
        z2 = (n_X - K_X) / n_X * corr_ratio**2 / (1 - corr_ratio**2) * 1/sci.stats.f.ppf(1 - p_level, f1, f2, loc=0, scale=1) - (K_X - 1)/n_X
        corr_ratio_conf_int_low = sqrt(z1) if sqrt(z1) >= 0 else 0
        corr_ratio_conf_int_high = sqrt(z2) if sqrt(z2) <= 1 else 1
    else:
        corr_ratio_conf_int_low = corr_ratio_conf_int_high = '-'    
    # оценка тесноты связи
    if scale=='Cheddok':
        conclusion_corr_ratio_scale = scale + ': ' + Cheddock_scale_check(corr_ratio, name=chr(951))
    elif scale=='Evans':
        conclusion_corr_ratio_scale = scale + ': ' + Evans_scale_check(corr_ratio, name=chr(951))
    # формируем результат            
    result = pd.DataFrame({
        'notation': (chr(951)),
        'coef_value': (corr_ratio),
        'coef_value_squared': (corr_ratio**2),
        'p_level': (p_level),
        'a_level': (a_level),
        'F_calc': (F_corr_ratio_calc),
        'F_table': (F_corr_ratio_table),
        'F_calc >= F_table': (F_corr_ratio_calc >= F_corr_ratio_table),
        'a_calc': (a_corr_ratio_calc),
        'a_calc <= a_level': (a_corr_ratio_calc <= a_level),
        'significance_check': (conclusion_corr_ratio_sign),
        'conf_int_low': (corr_ratio_conf_int_low),
        'conf_int_high': (corr_ratio_conf_int_high),
        'scale': (conclusion_corr_ratio_scale)
        },
        index=['Correlation ratio'])
    
    return result        

# Функция для проверка значимости линейной корреляционной связи
def line_corr_sign_check(X, Y, p_level=0.95, orientation='XY'):
    a_level = 1 - p_level
    X = np.array(X)
    Y = np.array(Y)
    n_X = len(X)
    n_Y = len(Y)    
    # коэффициент корреляции
    corr_coef = sci.stats.pearsonr(X, Y)[0]
    # корреляционное отношение
    try:
        if orientation!='XY' and orientation!='YX':
            raise ValueError("Error! Incorrect orientation!")
        if orientation=='XY':
            corr_ratio = corr_ratio_check(X, Y, orientation='XY', scale='Evans')['coef_value'].values[0]
        elif orientation=='YX':
            corr_ratio = corr_ratio_check(X, Y, orientation='YX', scale='Evans')['coef_value'].values[0]
    except ValueError as err:
        print(err)
    # число интервалов группировки
    group_int_number = lambda n: round (3.31*log(n_X, 10)+1) if round (3.31*log(n_X, 10)+1) >=2 else 2
    K_X = group_int_number(n_X)
    # проверка гипотезы о значимости линейной корреляционной связи
    if corr_ratio >= abs(corr_coef):
        F_line_corr_sign_calc = (n_X - K_X)/(K_X - 2) * (corr_ratio**2 - corr_coef**2) / (1 - corr_ratio**2)
        dfn = K_X - 2
        dfd = n_X - K_X
        F_line_corr_sign_table = sci.stats.f.ppf(p_level, dfn, dfd, loc=0, scale=1)
        comparison_F_calc_table = F_line_corr_sign_calc >= F_line_corr_sign_table
        a_line_corr_sign_calc = 1 - sci.stats.f.cdf(F_line_corr_sign_calc, dfn, dfd, loc=0, scale=1)
        comparison_a_calc_a_level = a_line_corr_sign_calc <= a_level
        conclusion_null_hypothesis_check = 'accepted' if F_line_corr_sign_calc < F_line_corr_sign_table else 'unaccepted'
        conclusion_line_corr_sign = 'linear' if conclusion_null_hypothesis_check == 'accepted' else 'non linear'
    else:
        F_line_corr_sign_calc = ''
        F_line_corr_sign_table = ''
        comparison_F_calc_table = ''
        a_line_corr_sign_calc = ''
        comparison_a_calc_a_level = ''
        conclusion_null_hypothesis_check = 'Attention! The correlation ratio is less than the correlation coefficient'
        conclusion_line_corr_sign = '-'
    
    # формируем результат            
    result = pd.DataFrame({
        'corr.coef.': (corr_coef),
        'corr.ratio.': (corr_ratio),
        'null hypothesis': ('ru00b2 = ' + chr(951) + 'u00b2'),
        'p_level': (p_level),
        'a_level': (a_level),
        'F_calc': (F_line_corr_sign_calc),
        'F_table': (F_line_corr_sign_table),
        'F_calc >= F_table': (comparison_F_calc_table),
        'a_calc': (a_line_corr_sign_calc),
        'a_calc <= a_level': (comparison_a_calc_a_level),
        'null_hypothesis_check': (conclusion_null_hypothesis_check),
        'significance_line_corr_check': (conclusion_line_corr_sign),
        },
        index=['Significance of linear correlation'])
    
    return result

display(corr_coef_check(X, Y, scale='Evans'))
display(corr_ratio_check(X, Y, orientation='XY', scale='Evans'))
display(line_corr_sign_check(X, Y, orientation='XY'))

Сделаем выводы по результатам расчетов:

1. Между величинами существует значимая (a_calc<0.05) корреляционная связь, корреляционное отношение η = 0.7936 (т.е. связь сильная по Эвансу).

2. Линейная корреляционная связь между величинами также значимая (a_calc<0.05), отрицательная, коэффициент корреляции r = -0.7189 (связь сильная по Эвансу); линейная корреляция между переменными объясняет 51.68% вариации.

3. Гипотеза о равенстве корреляционного отношения и коэффициента корреляции отвергается (a_calc<0.05), то есть отличие линейной формы связи от нелинейной значимо.

ИТОГИ

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

Исходный код находится в моем репозитории на GitHub (https://github.com/AANazarov/Statistical-methods).

ЛИТЕРАТУРА

1. Айвазян С.А. и др. Прикладная статистика: исследование зависимостей. – М.: Финансы и статистика, 1985. – 487 с.

2. Айвазян С.А., Мхитарян В.С. Прикладная статистика. Основы эконометрики: В 2 т. – Т.1: Теория вероятностей и прикладная статистика. – М.: ЮНИТИ-ДАНА, 2001. – 656 с.

3. Кобзарь А.И. Прикладная математическая статистика. Для инженеров и научных работников. – М.: ФИЗМАТЛИТ, 2006. – 816 с.

4. Котеров А.Н. и др. Сила связи. Сообщение 2. Градации величины корреляции. – Медицинская радиология и радиационная безопасность. 2019. Том 64. № 6. с.12–24 (https://medradiol.fmbafmbc.ru/journal_medradiol/abstracts/2019/6/12-24_Koterov_et_al.pdf).