{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Основы программирования в Python\n", "\n", "*Алла Тамбовцева*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Библиотека `scipy` для статистики" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Библиотека `scipy` (сокращение от *scientific Python*) включает в себя разные модули, позволяющие выполнять научные расчеты, решать задачи оптимизации, генерировать выборки из (псевдо)случайных величин с заданными параметрами, реализовывать статистические тесты и создавать статистические модели." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Импортируем библиотеку `pandas` и модуль `stats` из библиотеки `scipy`. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import scipy.stats as st # сократим название\n", "import pandas as pd" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Загрузим базу данных (датафрейм) из csv-файла. Описание данных см.[здесь](https://vincentarelbundock.github.io/Rdatasets/doc/datasets/swiss.html)." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "df = pd.read_csv('https://raw.githubusercontent.com/allatambov/Py-programming-3/master/add/swiss.csv')" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Unnamed: 0FertilityAgricultureExaminationEducationCatholicInfant.Mortality
0Courtelary80.217.015129.9622.2
1Delemont83.145.16984.8422.2
2Franches-Mnt92.539.75593.4020.2
3Moutier85.836.512733.7720.3
4Neuveville76.943.517155.1620.6
\n", "
" ], "text/plain": [ " Unnamed: 0 Fertility Agriculture Examination Education Catholic \\\n", "0 Courtelary 80.2 17.0 15 12 9.96 \n", "1 Delemont 83.1 45.1 6 9 84.84 \n", "2 Franches-Mnt 92.5 39.7 5 5 93.40 \n", "3 Moutier 85.8 36.5 12 7 33.77 \n", "4 Neuveville 76.9 43.5 17 15 5.16 \n", "\n", " Infant.Mortality \n", "0 22.2 \n", "1 22.2 \n", "2 20.2 \n", "3 20.3 \n", "4 20.6 " ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Далее мы попробуем реализовать различные статистические тесты (применить статистические критерии) для сравнения средних значений или распределений в двух выборках. Функции, встроенные в модуль `stats`, в отличие от аналогичных функций в R, требуют в качестве аргументов две выборки в явном виде, то есть два массива данных. В R обычно достаточно указать группирующую переменную через `~`, и набор значений автоматически разделится на две ожидаемые группы. \n", "\n", "Предположим, что нам необходимо сравнить средние значения уровня детской смертности в кантонах Швейцарии, где преобладает католическое население и где преобладает протестантское население. Сформируем две выборки на основе имеющихся данных: выберем соответствующие строки в таблице и возьмем столбец *Infant.Mortality*." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Unnamed: 0FertilityAgricultureExaminationEducationCatholicInfant.Mortality
1Delemont83.145.16984.8422.2
2Franches-Mnt92.539.75593.4020.2
5Porrentruy76.135.39790.5726.6
6Broye83.870.216792.8523.6
7Glane92.467.814897.1624.9
8Gruyere82.453.312797.6721.0
9Sarine82.945.2161391.3824.4
10Veveyse87.164.514698.6124.5
30Conthey75.585.93299.7115.1
31Entremont69.384.97699.6819.8
32Herens77.389.752100.0018.3
33Martigwy70.578.212698.9619.4
34Monthey79.464.97398.2220.2
35St Maurice65.075.99999.0617.8
36Sierre92.284.63399.4616.3
37Sion79.363.1131396.8318.1
45Rive Droite44.746.6162950.4318.2
46Rive Gauche42.827.7222958.3319.3
\n", "
" ], "text/plain": [ " Unnamed: 0 Fertility Agriculture Examination Education Catholic \\\n", "1 Delemont 83.1 45.1 6 9 84.84 \n", "2 Franches-Mnt 92.5 39.7 5 5 93.40 \n", "5 Porrentruy 76.1 35.3 9 7 90.57 \n", "6 Broye 83.8 70.2 16 7 92.85 \n", "7 Glane 92.4 67.8 14 8 97.16 \n", "8 Gruyere 82.4 53.3 12 7 97.67 \n", "9 Sarine 82.9 45.2 16 13 91.38 \n", "10 Veveyse 87.1 64.5 14 6 98.61 \n", "30 Conthey 75.5 85.9 3 2 99.71 \n", "31 Entremont 69.3 84.9 7 6 99.68 \n", "32 Herens 77.3 89.7 5 2 100.00 \n", "33 Martigwy 70.5 78.2 12 6 98.96 \n", "34 Monthey 79.4 64.9 7 3 98.22 \n", "35 St Maurice 65.0 75.9 9 9 99.06 \n", "36 Sierre 92.2 84.6 3 3 99.46 \n", "37 Sion 79.3 63.1 13 13 96.83 \n", "45 Rive Droite 44.7 46.6 16 29 50.43 \n", "46 Rive Gauche 42.8 27.7 22 29 58.33 \n", "\n", " Infant.Mortality \n", "1 22.2 \n", "2 20.2 \n", "5 26.6 \n", "6 23.6 \n", "7 24.9 \n", "8 21.0 \n", "9 24.4 \n", "10 24.5 \n", "30 15.1 \n", "31 19.8 \n", "32 18.3 \n", "33 19.4 \n", "34 20.2 \n", "35 17.8 \n", "36 16.3 \n", "37 18.1 \n", "45 18.2 \n", "46 19.3 " ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df[df.Catholic > 50] # для иллюстрации" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1 22.2\n", "2 20.2\n", "5 26.6\n", "6 23.6\n", "7 24.9\n", "8 21.0\n", "9 24.4\n", "10 24.5\n", "30 15.1\n", "31 19.8\n", "32 18.3\n", "33 19.4\n", "34 20.2\n", "35 17.8\n", "36 16.3\n", "37 18.1\n", "45 18.2\n", "46 19.3\n", "Name: Infant.Mortality, dtype: float64" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sample1 = df[df.Catholic > 50][\"Infant.Mortality\"] # выборка 1\n", "sample1" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0 22.2\n", "3 20.3\n", "4 20.6\n", "11 16.5\n", "12 19.1\n", "13 22.7\n", "14 18.7\n", "15 21.2\n", "16 20.0\n", "17 20.2\n", "18 10.8\n", "19 20.0\n", "20 18.0\n", "21 22.4\n", "22 16.7\n", "23 15.3\n", "24 21.0\n", "25 23.8\n", "26 18.0\n", "27 16.3\n", "28 20.9\n", "29 22.5\n", "38 20.3\n", "39 20.5\n", "40 18.9\n", "41 23.0\n", "42 20.0\n", "43 19.5\n", "44 18.0\n", "Name: Infant.Mortality, dtype: float64" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sample2 = df[df.Catholic <= 50][\"Infant.Mortality\"] # выборка 2\n", "sample2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Теперь приступим к формальной проверке гипотез." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### T-test (две выборки)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "T-test используется для сравнения средних значений двух генеральных совокупностей в предположении, что обе выборки взяты из нормального распределения. \n", "\n", "*Нулевая гипотеза:* средние значения двух генеральных совокупностей (откуда взяты выборки) равны. \n", "\n", "*Альтернативная гипотеза:* средние значения двух генеральных совокупностей не равны. \n", "\n", "В `stats` в t-тесте в качестве альтернативной гипотезы используется двусторонняя альтернатива (средние *не равны*), и всегда выводится соответствующее p-value (*two-tailed*). То же будет характерно для всех последующих тестов. Так как наши выборки независимы (не связаны, пример связанных выборок: значения уровня смертности в одних и тех же кантонах до и после какой-нибудь реформы), нам нужна функция `ttest_ind()`, от *independent*." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Ttest_indResult(statistic=1.1297965130690064, pvalue=0.26454837328688746)" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "st.ttest_ind(sample1, sample2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Что возвращает эта функция? Наблюдаемое значение t-статистики и p-value. Результат понятный, только более лакончиный по сравнению с тем, что выводит R.\n", "\n", "*Выводы:* так как p-value больше любого конвенционального уровня значимости (1%, 5%, 10%), на имеющихся данных на любом разумном уровне значимости нет оснований отвергнуть нулевую гипотезу. Средний уровень детской смертности в католических и протестантских районах можно считать одинаковым." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "По умолчанию считается, что дисперсии генеральных овокупностей равны. Часто это бывает не так, и такое предположение без формальной проверки и без содержательных соображений может казаться нереалистичным. Если мы предполагаем, что дисперсии генеральных совокупностей не равны, то это можно учесть, добавив аргумент `equal_var`:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Ttest_indResult(statistic=1.0863471703503398, pvalue=0.28551301767919196)" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "st.ttest_ind(sample1, sample2, equal_var = False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Принципиальных отличий в результатах не наблюдается." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Wilcoxon test & Mann-Whitney test (две выборки)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Если мы не можем считать распределение генеральных совокупностей, откуда взяты выборки, нормальным, то следует использовать методы, основанные не на самих наблюдениях в выборках, а на их рангах. Для сравнения распределений (иногда речь идет о сравнении медиан) используются тесты Уилкоксона и Манна-Уитни. Начнем с теста Уилкоксона (не проверяем, является ли распределение данных нормальным, просто для примера используем те же выборки)." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "ename": "ValueError", "evalue": "Unequal N in wilcoxon. Aborting.", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mst\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwilcoxon\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msample1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msample2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/scipy/stats/morestats.py\u001b[0m in \u001b[0;36mwilcoxon\u001b[0;34m(x, y, zero_method, correction)\u001b[0m\n\u001b[1;32m 2374\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0masarray\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2375\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2376\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Unequal N in wilcoxon. Aborting.'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2377\u001b[0m \u001b[0md\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mx\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2378\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mValueError\u001b[0m: Unequal N in wilcoxon. Aborting." ] } ], "source": [ "st.wilcoxon(sample1, sample2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Неудача! Проблема в том, что реализация критерия Уилкоксона в `stats` требует, чтобы выборки были одинакового размера. Но это ограничение можно обойти, просто выбрав другой критерий ‒ критерий Манна-Уитни. \n", "\n", "*Нулевая гипотеза:* выборки взяты из одного и того же распределения.\n", "\n", "*Альтернативная гипотеза:* выборки взяты из разных распределений." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "MannwhitneyuResult(statistic=235.5, pvalue=0.2920645577220585)" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "st.mannwhitneyu(sample1, sample2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Опять же, на имеющихся данных на любом уровне значимости нет оснований отвергнуть нулевую гипотезу. Выборки взяты из одного и того же распределения." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### ANOVA" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Если выборок больше, чем две, то использовать указанные выше критерии нельзя. В предположении, что все выборки взяты из нормального распределения, для сравнения средних значений более чем в двух группах используется однофакторный дисперсионный анализ (ANOVA, *analysis of variance*). \n", "\n", "*Нулевая гипотеза*: средние значения по всем группам (во всех генеральных совокупностях) равны.\n", "\n", "*Альтернативная гипотеза*: средние значения по всем группам (во всех генеральных совокупностях) не равны.\n", "\n", "Чтобы не создавать искусственные группы на основе данных в *swiss*, загрузим таблицу с весами цыплят, которых кормили разным кормом :) Описание см. [здесь](https://vincentarelbundock.github.io/Rdatasets/doc/datasets/chickwts.html)." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "dat = pd.read_csv('https://raw.githubusercontent.com/allatambov/Py-programming-3/master/add/chickwts.csv')" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Unnamed: 0weightfeed
01179horsebean
12160horsebean
23136horsebean
34227horsebean
45217horsebean
\n", "
" ], "text/plain": [ " Unnamed: 0 weight feed\n", "0 1 179 horsebean\n", "1 2 160 horsebean\n", "2 3 136 horsebean\n", "3 4 227 horsebean\n", "4 5 217 horsebean" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dat.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Задание:** разбить датафрейм на группы с помощью `groupby` по переменной `feed` и сохранить значения `weight` в словарь." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Решение:*" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "wgt = {}\n", "for name, d in dat.groupby('feed'):\n", " wgt[name] = d.weight" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'casein': 59 368\n", " 60 390\n", " 61 379\n", " 62 260\n", " 63 404\n", " 64 318\n", " 65 352\n", " 66 359\n", " 67 216\n", " 68 222\n", " 69 283\n", " 70 332\n", " Name: weight, dtype: int64, 'horsebean': 0 179\n", " 1 160\n", " 2 136\n", " 3 227\n", " 4 217\n", " 5 168\n", " 6 108\n", " 7 124\n", " 8 143\n", " 9 140\n", " Name: weight, dtype: int64, 'linseed': 10 309\n", " 11 229\n", " 12 181\n", " 13 141\n", " 14 260\n", " 15 203\n", " 16 148\n", " 17 169\n", " 18 213\n", " 19 257\n", " 20 244\n", " 21 271\n", " Name: weight, dtype: int64, 'meatmeal': 48 325\n", " 49 257\n", " 50 303\n", " 51 315\n", " 52 380\n", " 53 153\n", " 54 263\n", " 55 242\n", " 56 206\n", " 57 344\n", " 58 258\n", " Name: weight, dtype: int64, 'soybean': 22 243\n", " 23 230\n", " 24 248\n", " 25 327\n", " 26 329\n", " 27 250\n", " 28 193\n", " 29 271\n", " 30 316\n", " 31 267\n", " 32 199\n", " 33 171\n", " 34 158\n", " 35 248\n", " Name: weight, dtype: int64, 'sunflower': 36 423\n", " 37 340\n", " 38 392\n", " 39 339\n", " 40 341\n", " 41 226\n", " 42 320\n", " 43 295\n", " 44 334\n", " 45 322\n", " 46 297\n", " 47 318\n", " Name: weight, dtype: int64}" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "wgt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Теперь ANOVA (`f_oneway` от *One-Way ANOVA*):" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "F_onewayResult(statistic=15.364799774712534, pvalue=5.936419853471331e-10)" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "st.f_oneway(wgt['casein'], wgt['horsebean'], wgt['linseed'], wgt['meatmeal'], wgt['soybean'], wgt['sunflower'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Функция возвращает наблюдаемое значение F-статистики и p-value. В данном случае p-value близко к 0, поэтому гипотезу о равенстве средних генеральных совокупностей по группам можно отвергнуть на 1% уровне значимости. Средний вес цыплят, которых кормили разным кормом, отличается (еще бы, horsebean или sunflower!)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Kruskal-Wallis test" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Критерий Краскела-Уоллиса используется, когда нам необходимо сравнить распределения более, чем в двух группах в предположении, что выборки взяты не из нормального распределения (распределения неизвестны). \n", "\n", "*Нулевая гипотеза*: выборки взяты из одного и того же распределения.\n", "\n", "*Альтернативная гипотеза*: выборки взяты из разных распределений." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "KruskalResult(statistic=37.34271769425624, pvalue=5.112829511937094e-07)" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "st.kruskal(wgt['casein'], wgt['horsebean'], wgt['linseed'], wgt['meatmeal'], wgt['soybean'], wgt['sunflower'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Выводы аналогичны полученным выше. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Линейная регрессия (парная)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "К сожалению, построение обычных регрессионных моделей, не связанных с машинным обучением, в Python, является нетривиальной задачей. Конечно, построить парную или множественную регрессию в Python не составит особенного труда, но вот получить готовую и красивую выдачу с результатами, как в R, сразу не получится (речь идет о библиотеках `scipy` и `sklearn`). Для того чтобы построить более сложную модель, потребуется больше сил, и сделать это без четкого понимания, как эта модель устроена, хорошо не получится, так как многие вещи придется дописывать самостоятельно.\n", "\n", "Но мы начнем с самого простого, с парной линейной регрессии. В качестве зависимой переменной выберем показатель *Agriculture* (процент мужчин, занятых в сельском хозяйстве), а в качестве независимой ‒ *Catholic* (процент католического населения в кантоне). Построим модель с помощью функции `linregress`:" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "LinregressResult(slope=0.736535100477301, intercept=3.8312750162456837, rvalue=0.4010950530487398, pvalue=0.005204433539191572, stderr=0.25075675209727843)" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "st.linregress(df.Agriculture, df.Catholic)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Эта функция вернула следующие значения: коэффициент при независимой переменной (*slope*), константу (*intercept*), коэффициент корреляции (*r*), pvalue, стандартную ошибку коэффициента (*stderr*). Для удобства можно сохранить их отдельно, используя списковое присваивание:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "slope, intercept, rvalue, pvalue, stderr = st.linregress(df.Agriculture, df.Catholic)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "И спокойно вызывать по отдельности:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.736535100477301" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "slope" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.005204433539191572" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pvalue" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Давайте попробуем оформить результаты в выдачу, хотя бы отдаленно напоминающую аккуратную и подробную выдачу в R. Сначала сохраним все полученные значения в словарь, а затем из него сделаем маленький датафрейм. " ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "d = {'slope': slope, 'intercept': intercept, 'pvalue': pvalue, 'stderr': stderr}" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'intercept': 3.8312750162456837,\n", " 'pvalue': 0.005204433539191572,\n", " 'slope': 0.736535100477301,\n", " 'stderr': 0.25075675209727843}" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
0
slope0.736535
intercept3.831275
stderr0.250757
pvalue0.005204
\n", "
" ], "text/plain": [ " 0\n", "slope 0.736535\n", "intercept 3.831275\n", "stderr 0.250757\n", "pvalue 0.005204" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "res = pd.DataFrame.from_dict(d, orient='index')\n", "res" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Транспонируем (поменяем местами строки и столбцы):" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
slopeinterceptstderrpvalue
00.7365353.8312750.2507570.005204
\n", "
" ], "text/plain": [ " slope intercept stderr pvalue\n", "0 0.736535 3.831275 0.250757 0.005204" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "new = res.transpose()\n", "new" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Поменяем местами столбцы, чтобы они шли в таком порядке:\n", "\n", " intercept, slope, stderr, pvalue" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
interceptslopestderrpvalue
03.8312750.7365350.2507570.005204
\n", "
" ], "text/plain": [ " intercept slope stderr pvalue\n", "0 3.831275 0.736535 0.250757 0.005204" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cols = ['intercept', 'slope', 'stderr', 'pvalue']\n", "new = new[cols]\n", "new" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Теперь хорошо бы добавить звездочки, которые обычно показывают, на каком уровне значимости коэффициент является значимым. В R звездочки ставятся так: \n", "\n", "* 5% уровень значимости: $*$\n", "* 1% уровень значимости: $**$\n", "* 0.1% уровень значимости: $***$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Задача:** написать функцию `star`, которая принимает на вход датафрейм с результатами, и возвращает строку из звездочек, число которых зависит от p-value." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Решение:* " ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "def star(D):\n", " if D.pvalue[0] > 0.01 and D.pvalue[0] <= 0.05:\n", " return '*'\n", " elif D.pvalue[0] > 0.001 and D.pvalue[0] <= 0.01:\n", " return '**'\n", " elif D.pvalue[0] <= 0.001:\n", " return '***'\n", " else:\n", " return ''" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'**'" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "star(new) # верно" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Теперь напишем функцию, которая будет \"доклеивать\" звездочки к коэффициенту при независимой переменной и которую можно будет применить к столбцу `slope` в таблице `new`." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "def add_star(x, D=new):\n", " s = star(D)\n", " r = str(x) + s\n", " return r" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/oem/.local/lib/python3.5/site-packages/ipykernel_launcher.py:1: SettingWithCopyWarning: \n", "A value is trying to be set on a copy of a slice from a DataFrame.\n", "Try using .loc[row_indexer,col_indexer] = value instead\n", "\n", "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n", " \"\"\"Entry point for launching an IPython kernel.\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
interceptslopestderrpvaluecoef
03.8312750.7365350.2507570.0052040.736535100477301**
\n", "
" ], "text/plain": [ " intercept slope stderr pvalue coef\n", "0 3.831275 0.736535 0.250757 0.005204 0.736535100477301**" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "new['coef'] = new['slope'].apply(add_star)\n", "new" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Новый столбец со \"звездным\" коэффициентом добавился. Однако выглядит это не очень хорошо ‒ много знаков после запятой. Сократим до двух, поправив нашу функцию, используя форматирование строк." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "def add_star(x, D=new):\n", " s = star(D)\n", " r = str('{:.2f}'.format(x)) + s\n", " return r" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/oem/.local/lib/python3.5/site-packages/ipykernel_launcher.py:1: SettingWithCopyWarning: \n", "A value is trying to be set on a copy of a slice from a DataFrame.\n", "Try using .loc[row_indexer,col_indexer] = value instead\n", "\n", "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n", " \"\"\"Entry point for launching an IPython kernel.\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
interceptslopestderrpvaluecoef
03.8312750.7365350.2507570.0052040.74**
\n", "
" ], "text/plain": [ " intercept slope stderr pvalue coef\n", "0 3.831275 0.736535 0.250757 0.005204 0.74**" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "new['coef'] = new['slope'].apply(add_star)\n", "new" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "При желании можно убрать лишние знаки после запятой во всех столбцах. И это можно сделать, не затрагивая старую таблицу." ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
interceptslopestderrpvalue
03.830.740.250.01
\n", "
" ], "text/plain": [ " intercept slope stderr pvalue\n", "0 3.83 0.74 0.25 0.01" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# без coef - он у нас текстовый\n", "new[['intercept', 'slope', 'stderr', 'pvalue']].applymap(\"{0:.2f}\".format)" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "final = new[['intercept', 'slope', 'stderr', 'pvalue']].applymap(\"{0:.2f}\".format)\n", "final['coef'] = new['coef']" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
interceptslopestderrpvaluecoef
03.830.740.250.010.74**
\n", "
" ], "text/plain": [ " intercept slope stderr pvalue coef\n", "0 3.83 0.74 0.25 0.01 0.74**" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "final" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Допишем вместо индекса 0 строку с переменными в модели:" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [], "source": [ "final.index = ['Agriculture ~ Catholic']" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
interceptslopestderrpvaluecoef
Agriculture ~ Catholic3.830.740.250.010.74**
\n", "
" ], "text/plain": [ " intercept slope stderr pvalue coef\n", "Agriculture ~ Catholic 3.83 0.74 0.25 0.01 0.74**" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "final" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Теперь можем выгрузить эту таблицу в LaTeX (для тех, кто им пользуется):" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'\\\\begin{tabular}{llllll}\\n\\\\toprule\\n{} & intercept & slope & stderr & pvalue & coef \\\\\\\\\\n\\\\midrule\\nAgriculture \\\\textasciitilde Catholic & 3.83 & 0.74 & 0.25 & 0.01 & 0.74** \\\\\\\\\\n\\\\bottomrule\\n\\\\end{tabular}\\n'" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "final.to_latex()" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\\begin{tabular}{llllll}\n", "\\toprule\n", "{} & intercept & slope & stderr & pvalue & coef \\\\\n", "\\midrule\n", "Agriculture \\textasciitilde Catholic & 3.83 & 0.74 & 0.25 & 0.01 & 0.74** \\\\\n", "\\bottomrule\n", "\\end{tabular}\n", "\n" ] } ], "source": [ "print(final.to_latex()) # код можно просто скопировать" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Более подробно об опциях, которые можно учесть при выгрузке, можно узнать, вызвав `help(final.to_latex)`. И да, выгрузить в LaTeX можно любой датафрейм `pandas`!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Если нужно построить множественную линейную регрессию, то тут уже модулем `stats` не обойтись. Нужна другая, более мощная библиотека, которая хорошо знакома тем, кто занимается машинным обучением, а именно `sklearn`.\n", "\n", "Импортируем функцию для линейнойй модели:" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [], "source": [ "import sklearn.linear_model as lm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Создадим модель (точнее, объект `model`, который относится к классу `LinearRegression`, но про классы мы поговорим чуть позже)." ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [], "source": [ "model = lm.LinearRegression()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Теперь выберем из датафрейма набор независимых переменных (`X`):" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [], "source": [ "X = df[['Catholic', 'Infant.Mortality']]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Теперь осталось построить модель, в которой в качестве зависимой переменной будет выступать *Agriculture*, а в качестве независимых ‒ *Catholic* и *Infant.Mortality*. Обратите внимание: непривычно, но сначала должен быть указан набор **независимых** переменных, а только потом **зависимая**." ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.fit(X, df.Agriculture)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Модель построена. Можем получить коэффициенты при независимых переменных:" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 0.23136647, -1.05591198])" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.coef_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Значение константы:" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "62.197852894848076" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.intercept_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Предсказанные значения *Agriculture*:" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([41.06101685, 58.3857378 , 62.47805871, 48.57608516, 41.63991698,\n", " 55.06545491, 58.76070641, 58.3852103 , 62.62126393, 57.57586812,\n", " 59.14305646, 46.74654744, 42.55513587, 39.2536043 , 43.10475222,\n", " 45.4115873 , 41.84312255, 43.67027871, 51.29144137, 41.73669397,\n", " 44.4014838 , 39.59120088, 48.06701105, 47.01413869, 40.57898075,\n", " 38.27719429, 43.78373533, 46.77263667, 44.40031738, 39.85116869,\n", " 69.32313223, 64.35340491, 66.01131015, 64.60918585, 63.59324507,\n", " 66.32178167, 67.99819623, 65.48906085, 42.06311915, 43.74220078,\n", " 44.83704814, 41.82659786, 42.22950455, 43.60888913, 52.98749334,\n", " 54.64806565, 55.31435754])" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.predict(X)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Можем сравнить их с реальными значениями *Agriculture*:" ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0 24.061017\n", "1 13.285738\n", "2 22.778059\n", "3 12.076085\n", "4 -1.860083\n", "5 19.765455\n", "6 -11.439294\n", "7 -9.414790\n", "8 9.321264\n", "9 12.375868\n", "10 -5.356944\n", "11 -15.253453\n", "12 -24.944864\n", "13 -21.446396\n", "14 -26.195248\n", "15 -27.188413\n", "16 7.843123\n", "17 24.270279\n", "18 36.091441\n", "19 -31.263306\n", "20 -15.398516\n", "21 -15.508799\n", "22 -2.832989\n", "23 -7.085861\n", "24 -30.621019\n", "25 -19.822806\n", "26 -19.716265\n", "27 -14.027363\n", "28 17.600317\n", "29 -9.648831\n", "30 -16.576868\n", "31 -20.546595\n", "32 -23.688690\n", "33 -13.590814\n", "34 -1.306755\n", "35 -9.578218\n", "36 -16.601804\n", "37 2.389061\n", "38 3.663119\n", "39 36.042201\n", "40 28.137048\n", "41 24.226598\n", "42 4.629505\n", "43 24.908889\n", "44 51.787493\n", "45 8.048066\n", "46 27.614358\n", "Name: Agriculture, dtype: float64" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.predict(X)-df.Agriculture" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Посчитать сумму квадратов остатков:" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "19487.961646821" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sum((model.predict(X)-df.Agriculture)**2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "И убедиться, что это значение, посчитанное вручную, верное:" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/oem/.local/lib/python3.5/site-packages/sklearn/utils/deprecation.py:70: DeprecationWarning: Function residues_ is deprecated; ``residues_`` is deprecated and will be removed in 0.19\n", " warnings.warn(msg, category=DeprecationWarning)\n" ] }, { "data": { "text/plain": [ "19487.961646820993" ] }, "execution_count": 44, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.residues_ # остатки" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Единственное, получить статистическую значимость коэффициентов так же просто и быстро не получится. Придется, например, подключать для расчетов `scipy`. Мораль: если нужно построить обычные регрессионные модели и работать с ними в контексте статистики и \"классической\" эконометрики, без машинного обучения, лучше использовать тот же R, где с помощью двух строк кода можно получить модель и всю сводную информацию по ней." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.5.2" } }, "nbformat": 4, "nbformat_minor": 2 }