Основы программирования в Python

Алла Тамбовцева

Библиотека scipy для статистики

Библиотека scipy (сокращение от scientific Python) включает в себя разные модули, позволяющие выполнять научные расчеты, решать задачи оптимизации, генерировать выборки из (псевдо)случайных величин с заданными параметрами, реализовывать статистические тесты и создавать статистические модели.

Импортируем библиотеку pandas и модуль stats из библиотеки scipy.

In [1]:
import scipy.stats as st # сократим название
import pandas as pd

Загрузим базу данных (датафрейм) из csv-файла. Описание данных см.здесь.

In [2]:
df = pd.read_csv('https://raw.githubusercontent.com/allatambov/Py-programming-3/master/add/swiss.csv')
In [3]:
df.head()
Out[3]:
Unnamed: 0 Fertility Agriculture Examination Education Catholic Infant.Mortality
0 Courtelary 80.2 17.0 15 12 9.96 22.2
1 Delemont 83.1 45.1 6 9 84.84 22.2
2 Franches-Mnt 92.5 39.7 5 5 93.40 20.2
3 Moutier 85.8 36.5 12 7 33.77 20.3
4 Neuveville 76.9 43.5 17 15 5.16 20.6

Далее мы попробуем реализовать различные статистические тесты (применить статистические критерии) для сравнения средних значений или распределений в двух выборках. Функции, встроенные в модуль stats, в отличие от аналогичных функций в R, требуют в качестве аргументов две выборки в явном виде, то есть два массива данных. В R обычно достаточно указать группирующую переменную через ~, и набор значений автоматически разделится на две ожидаемые группы.

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

In [4]:
df[df.Catholic > 50] # для иллюстрации
Out[4]:
Unnamed: 0 Fertility Agriculture Examination Education Catholic Infant.Mortality
1 Delemont 83.1 45.1 6 9 84.84 22.2
2 Franches-Mnt 92.5 39.7 5 5 93.40 20.2
5 Porrentruy 76.1 35.3 9 7 90.57 26.6
6 Broye 83.8 70.2 16 7 92.85 23.6
7 Glane 92.4 67.8 14 8 97.16 24.9
8 Gruyere 82.4 53.3 12 7 97.67 21.0
9 Sarine 82.9 45.2 16 13 91.38 24.4
10 Veveyse 87.1 64.5 14 6 98.61 24.5
30 Conthey 75.5 85.9 3 2 99.71 15.1
31 Entremont 69.3 84.9 7 6 99.68 19.8
32 Herens 77.3 89.7 5 2 100.00 18.3
33 Martigwy 70.5 78.2 12 6 98.96 19.4
34 Monthey 79.4 64.9 7 3 98.22 20.2
35 St Maurice 65.0 75.9 9 9 99.06 17.8
36 Sierre 92.2 84.6 3 3 99.46 16.3
37 Sion 79.3 63.1 13 13 96.83 18.1
45 Rive Droite 44.7 46.6 16 29 50.43 18.2
46 Rive Gauche 42.8 27.7 22 29 58.33 19.3
In [5]:
sample1 = df[df.Catholic > 50]["Infant.Mortality"] # выборка 1
sample1
Out[5]:
1     22.2
2     20.2
5     26.6
6     23.6
7     24.9
8     21.0
9     24.4
10    24.5
30    15.1
31    19.8
32    18.3
33    19.4
34    20.2
35    17.8
36    16.3
37    18.1
45    18.2
46    19.3
Name: Infant.Mortality, dtype: float64
In [6]:
sample2 = df[df.Catholic <= 50]["Infant.Mortality"] # выборка 2
sample2
Out[6]:
0     22.2
3     20.3
4     20.6
11    16.5
12    19.1
13    22.7
14    18.7
15    21.2
16    20.0
17    20.2
18    10.8
19    20.0
20    18.0
21    22.4
22    16.7
23    15.3
24    21.0
25    23.8
26    18.0
27    16.3
28    20.9
29    22.5
38    20.3
39    20.5
40    18.9
41    23.0
42    20.0
43    19.5
44    18.0
Name: Infant.Mortality, dtype: float64

Теперь приступим к формальной проверке гипотез.

T-test (две выборки)

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

Нулевая гипотеза: средние значения двух генеральных совокупностей (откуда взяты выборки) равны.

Альтернативная гипотеза: средние значения двух генеральных совокупностей не равны.

В stats в t-тесте в качестве альтернативной гипотезы используется двусторонняя альтернатива (средние не равны), и всегда выводится соответствующее p-value (two-tailed). То же будет характерно для всех последующих тестов. Так как наши выборки независимы (не связаны, пример связанных выборок: значения уровня смертности в одних и тех же кантонах до и после какой-нибудь реформы), нам нужна функция ttest_ind(), от independent.

In [12]:
st.ttest_ind(sample1, sample2)
Out[12]:
Ttest_indResult(statistic=1.1297965130690064, pvalue=0.26454837328688746)

Что возвращает эта функция? Наблюдаемое значение t-статистики и p-value. Результат понятный, только более лакончиный по сравнению с тем, что выводит R.

Выводы: так как p-value больше любого конвенционального уровня значимости (1%, 5%, 10%), на имеющихся данных на любом разумном уровне значимости нет оснований отвергнуть нулевую гипотезу. Средний уровень детской смертности в католических и протестантских районах можно считать одинаковым.

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

In [13]:
st.ttest_ind(sample1, sample2, equal_var = False)
Out[13]:
Ttest_indResult(statistic=1.0863471703503398, pvalue=0.28551301767919196)

Принципиальных отличий в результатах не наблюдается.

Wilcoxon test & Mann-Whitney test (две выборки)

Если мы не можем считать распределение генеральных совокупностей, откуда взяты выборки, нормальным, то следует использовать методы, основанные не на самих наблюдениях в выборках, а на их рангах. Для сравнения распределений (иногда речь идет о сравнении медиан) используются тесты Уилкоксона и Манна-Уитни. Начнем с теста Уилкоксона (не проверяем, является ли распределение данных нормальным, просто для примера используем те же выборки).

In [15]:
st.wilcoxon(sample1, sample2)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-15-6bfff17a9e4f> in <module>()
----> 1 st.wilcoxon(sample1, sample2)

/usr/local/lib/python3.5/dist-packages/scipy/stats/morestats.py in wilcoxon(x, y, zero_method, correction)
   2374         x, y = map(asarray, (x, y))
   2375         if len(x) != len(y):
-> 2376             raise ValueError('Unequal N in wilcoxon.  Aborting.')
   2377         d = x - y
   2378 

ValueError: Unequal N in wilcoxon.  Aborting.

Неудача! Проблема в том, что реализация критерия Уилкоксона в stats требует, чтобы выборки были одинакового размера. Но это ограничение можно обойти, просто выбрав другой критерий ‒ критерий Манна-Уитни.

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

Альтернативная гипотеза: выборки взяты из разных распределений.

In [16]:
st.mannwhitneyu(sample1, sample2)
Out[16]:
MannwhitneyuResult(statistic=235.5, pvalue=0.2920645577220585)

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

ANOVA

Если выборок больше, чем две, то использовать указанные выше критерии нельзя. В предположении, что все выборки взяты из нормального распределения, для сравнения средних значений более чем в двух группах используется однофакторный дисперсионный анализ (ANOVA, analysis of variance).

Нулевая гипотеза: средние значения по всем группам (во всех генеральных совокупностях) равны.

Альтернативная гипотеза: средние значения по всем группам (во всех генеральных совокупностях) не равны.

Чтобы не создавать искусственные группы на основе данных в swiss, загрузим таблицу с весами цыплят, которых кормили разным кормом :) Описание см. здесь.

In [7]:
dat = pd.read_csv('https://raw.githubusercontent.com/allatambov/Py-programming-3/master/add/chickwts.csv')
In [8]:
dat.head()
Out[8]:
Unnamed: 0 weight feed
0 1 179 horsebean
1 2 160 horsebean
2 3 136 horsebean
3 4 227 horsebean
4 5 217 horsebean

Задание: разбить датафрейм на группы с помощью groupby по переменной feed и сохранить значения weight в словарь.

Решение:

In [9]:
wgt = {}
for name, d in dat.groupby('feed'):
    wgt[name] = d.weight
In [10]:
wgt
Out[10]:
{'casein': 59    368
 60    390
 61    379
 62    260
 63    404
 64    318
 65    352
 66    359
 67    216
 68    222
 69    283
 70    332
 Name: weight, dtype: int64, 'horsebean': 0    179
 1    160
 2    136
 3    227
 4    217
 5    168
 6    108
 7    124
 8    143
 9    140
 Name: weight, dtype: int64, 'linseed': 10    309
 11    229
 12    181
 13    141
 14    260
 15    203
 16    148
 17    169
 18    213
 19    257
 20    244
 21    271
 Name: weight, dtype: int64, 'meatmeal': 48    325
 49    257
 50    303
 51    315
 52    380
 53    153
 54    263
 55    242
 56    206
 57    344
 58    258
 Name: weight, dtype: int64, 'soybean': 22    243
 23    230
 24    248
 25    327
 26    329
 27    250
 28    193
 29    271
 30    316
 31    267
 32    199
 33    171
 34    158
 35    248
 Name: weight, dtype: int64, 'sunflower': 36    423
 37    340
 38    392
 39    339
 40    341
 41    226
 42    320
 43    295
 44    334
 45    322
 46    297
 47    318
 Name: weight, dtype: int64}

Теперь ANOVA (f_oneway от One-Way ANOVA):

In [11]:
st.f_oneway(wgt['casein'], wgt['horsebean'], wgt['linseed'], wgt['meatmeal'], wgt['soybean'], wgt['sunflower'])
Out[11]:
F_onewayResult(statistic=15.364799774712534, pvalue=5.936419853471331e-10)

Функция возвращает наблюдаемое значение F-статистики и p-value. В данном случае p-value близко к 0, поэтому гипотезу о равенстве средних генеральных совокупностей по группам можно отвергнуть на 1% уровне значимости. Средний вес цыплят, которых кормили разным кормом, отличается (еще бы, horsebean или sunflower!).

Kruskal-Wallis test

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

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

Альтернативная гипотеза: выборки взяты из разных распределений.

In [12]:
st.kruskal(wgt['casein'], wgt['horsebean'], wgt['linseed'], wgt['meatmeal'], wgt['soybean'], wgt['sunflower'])
Out[12]:
KruskalResult(statistic=37.34271769425624, pvalue=5.112829511937094e-07)

Выводы аналогичны полученным выше.

Линейная регрессия (парная)

К сожалению, построение обычных регрессионных моделей, не связанных с машинным обучением, в Python, является нетривиальной задачей. Конечно, построить парную или множественную регрессию в Python не составит особенного труда, но вот получить готовую и красивую выдачу с результатами, как в R, сразу не получится (речь идет о библиотеках scipy и sklearn). Для того чтобы построить более сложную модель, потребуется больше сил, и сделать это без четкого понимания, как эта модель устроена, хорошо не получится, так как многие вещи придется дописывать самостоятельно.

Но мы начнем с самого простого, с парной линейной регрессии. В качестве зависимой переменной выберем показатель Agriculture (процент мужчин, занятых в сельском хозяйстве), а в качестве независимой ‒ Catholic (процент католического населения в кантоне). Построим модель с помощью функции linregress:

In [24]:
st.linregress(df.Agriculture, df.Catholic)
Out[24]:
LinregressResult(slope=0.736535100477301, intercept=3.8312750162456837, rvalue=0.4010950530487398, pvalue=0.005204433539191572, stderr=0.25075675209727843)

Эта функция вернула следующие значения: коэффициент при независимой переменной (slope), константу (intercept), коэффициент корреляции (r), pvalue, стандартную ошибку коэффициента (stderr). Для удобства можно сохранить их отдельно, используя списковое присваивание:

In [13]:
slope, intercept, rvalue, pvalue, stderr = st.linregress(df.Agriculture, df.Catholic)

И спокойно вызывать по отдельности:

In [14]:
slope
Out[14]:
0.736535100477301
In [15]:
pvalue
Out[15]:
0.005204433539191572

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

In [16]:
d = {'slope': slope, 'intercept': intercept, 'pvalue': pvalue, 'stderr': stderr}
In [17]:
d
Out[17]:
{'intercept': 3.8312750162456837,
 'pvalue': 0.005204433539191572,
 'slope': 0.736535100477301,
 'stderr': 0.25075675209727843}
In [18]:
res = pd.DataFrame.from_dict(d, orient='index')
res
Out[18]:
0
slope 0.736535
intercept 3.831275
stderr 0.250757
pvalue 0.005204

Транспонируем (поменяем местами строки и столбцы):

In [19]:
new = res.transpose()
new
Out[19]:
slope intercept stderr pvalue
0 0.736535 3.831275 0.250757 0.005204

Поменяем местами столбцы, чтобы они шли в таком порядке:

intercept, slope, stderr, pvalue
In [20]:
cols = ['intercept', 'slope', 'stderr', 'pvalue']
new = new[cols]
new
Out[20]:
intercept slope stderr pvalue
0 3.831275 0.736535 0.250757 0.005204

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

  • 5% уровень значимости: $*$
  • 1% уровень значимости: $**$
  • 0.1% уровень значимости: $***$

Задача: написать функцию star, которая принимает на вход датафрейм с результатами, и возвращает строку из звездочек, число которых зависит от p-value.

Решение:

In [21]:
def star(D):
    if D.pvalue[0] > 0.01 and D.pvalue[0] <= 0.05:
        return '*'
    elif D.pvalue[0] > 0.001 and D.pvalue[0] <= 0.01:
        return '**'
    elif D.pvalue[0] <= 0.001:
        return '***'
    else:
        return ''
In [22]:
star(new) # верно
Out[22]:
'**'

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

In [23]:
def add_star(x, D=new):
    s = star(D)
    r = str(x) + s
    return r
In [24]:
new['coef'] = new['slope'].apply(add_star)
new
/home/oem/.local/lib/python3.5/site-packages/ipykernel_launcher.py:1: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  """Entry point for launching an IPython kernel.
Out[24]:
intercept slope stderr pvalue coef
0 3.831275 0.736535 0.250757 0.005204 0.736535100477301**

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

In [25]:
def add_star(x, D=new):
    s = star(D)
    r = str('{:.2f}'.format(x)) + s
    return r
In [26]:
new['coef'] = new['slope'].apply(add_star)
new
/home/oem/.local/lib/python3.5/site-packages/ipykernel_launcher.py:1: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  """Entry point for launching an IPython kernel.
Out[26]:
intercept slope stderr pvalue coef
0 3.831275 0.736535 0.250757 0.005204 0.74**

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

In [27]:
# без coef - он у нас текстовый
new[['intercept', 'slope', 'stderr', 'pvalue']].applymap("{0:.2f}".format)
Out[27]:
intercept slope stderr pvalue
0 3.83 0.74 0.25 0.01
In [28]:
final = new[['intercept', 'slope', 'stderr', 'pvalue']].applymap("{0:.2f}".format)
final['coef'] = new['coef']
In [29]:
final
Out[29]:
intercept slope stderr pvalue coef
0 3.83 0.74 0.25 0.01 0.74**

Допишем вместо индекса 0 строку с переменными в модели:

In [30]:
final.index = ['Agriculture ~ Catholic']
In [31]:
final
Out[31]:
intercept slope stderr pvalue coef
Agriculture ~ Catholic 3.83 0.74 0.25 0.01 0.74**

Теперь можем выгрузить эту таблицу в LaTeX (для тех, кто им пользуется):

In [32]:
final.to_latex()
Out[32]:
'\\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'
In [33]:
print(final.to_latex()) # код можно просто скопировать
\begin{tabular}{llllll}
\toprule
{} & intercept & slope & stderr & pvalue &    coef \\
\midrule
Agriculture \textasciitilde Catholic &      3.83 &  0.74 &   0.25 &   0.01 &  0.74** \\
\bottomrule
\end{tabular}

Более подробно об опциях, которые можно учесть при выгрузке, можно узнать, вызвав help(final.to_latex). И да, выгрузить в LaTeX можно любой датафрейм pandas!

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

Импортируем функцию для линейнойй модели:

In [34]:
import sklearn.linear_model as lm

Создадим модель (точнее, объект model, который относится к классу LinearRegression, но про классы мы поговорим чуть позже).

In [36]:
model = lm.LinearRegression()

Теперь выберем из датафрейма набор независимых переменных (X):

In [37]:
X = df[['Catholic', 'Infant.Mortality']]

Теперь осталось построить модель, в которой в качестве зависимой переменной будет выступать Agriculture, а в качестве независимых ‒ Catholic и Infant.Mortality. Обратите внимание: непривычно, но сначала должен быть указан набор независимых переменных, а только потом зависимая.

In [38]:
model.fit(X, df.Agriculture)
Out[38]:
LinearRegression(copy_X=True, fit_intercept=True, n_jobs=1, normalize=False)

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

In [39]:
model.coef_
Out[39]:
array([ 0.23136647, -1.05591198])

Значение константы:

In [40]:
model.intercept_
Out[40]:
62.197852894848076

Предсказанные значения Agriculture:

In [41]:
model.predict(X)
Out[41]:
array([41.06101685, 58.3857378 , 62.47805871, 48.57608516, 41.63991698,
       55.06545491, 58.76070641, 58.3852103 , 62.62126393, 57.57586812,
       59.14305646, 46.74654744, 42.55513587, 39.2536043 , 43.10475222,
       45.4115873 , 41.84312255, 43.67027871, 51.29144137, 41.73669397,
       44.4014838 , 39.59120088, 48.06701105, 47.01413869, 40.57898075,
       38.27719429, 43.78373533, 46.77263667, 44.40031738, 39.85116869,
       69.32313223, 64.35340491, 66.01131015, 64.60918585, 63.59324507,
       66.32178167, 67.99819623, 65.48906085, 42.06311915, 43.74220078,
       44.83704814, 41.82659786, 42.22950455, 43.60888913, 52.98749334,
       54.64806565, 55.31435754])

Можем сравнить их с реальными значениями Agriculture:

In [42]:
model.predict(X)-df.Agriculture
Out[42]:
0     24.061017
1     13.285738
2     22.778059
3     12.076085
4     -1.860083
5     19.765455
6    -11.439294
7     -9.414790
8      9.321264
9     12.375868
10    -5.356944
11   -15.253453
12   -24.944864
13   -21.446396
14   -26.195248
15   -27.188413
16     7.843123
17    24.270279
18    36.091441
19   -31.263306
20   -15.398516
21   -15.508799
22    -2.832989
23    -7.085861
24   -30.621019
25   -19.822806
26   -19.716265
27   -14.027363
28    17.600317
29    -9.648831
30   -16.576868
31   -20.546595
32   -23.688690
33   -13.590814
34    -1.306755
35    -9.578218
36   -16.601804
37     2.389061
38     3.663119
39    36.042201
40    28.137048
41    24.226598
42     4.629505
43    24.908889
44    51.787493
45     8.048066
46    27.614358
Name: Agriculture, dtype: float64

Посчитать сумму квадратов остатков:

In [43]:
sum((model.predict(X)-df.Agriculture)**2)
Out[43]:
19487.961646821

И убедиться, что это значение, посчитанное вручную, верное:

In [44]:
model.residues_ # остатки
/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
  warnings.warn(msg, category=DeprecationWarning)
Out[44]:
19487.961646820993

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