#!/usr/bin/env python
# coding: utf-8

# # Основы программирования в Python
# 
# *Алла Тамбовцева*

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

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

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

# In[1]:


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


# Загрузим базу данных (датафрейм) из csv-файла. Описание данных см.[здесь](https://vincentarelbundock.github.io/Rdatasets/doc/datasets/swiss.html).

# In[2]:


df = pd.read_csv('https://raw.githubusercontent.com/allatambov/Py-programming-3/master/add/swiss.csv')


# In[3]:


df.head()


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

# In[4]:


df[df.Catholic > 50] # для иллюстрации


# In[5]:


sample1 = df[df.Catholic > 50]["Infant.Mortality"] # выборка 1
sample1


# In[6]:


sample2 = df[df.Catholic <= 50]["Infant.Mortality"] # выборка 2
sample2


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

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

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

# In[12]:


st.ttest_ind(sample1, sample2)


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

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

# In[13]:


st.ttest_ind(sample1, sample2, equal_var = False)


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

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

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

# In[15]:


st.wilcoxon(sample1, sample2)


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

# In[16]:


st.mannwhitneyu(sample1, sample2)


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

# ### ANOVA

# Если выборок больше, чем две, то использовать указанные выше критерии нельзя. В предположении, что все выборки взяты из нормального распределения, для сравнения средних значений более чем в двух группах используется однофакторный дисперсионный анализ (ANOVA, *analysis of variance*). 
# 
# *Нулевая гипотеза*: средние значения по всем группам (во всех генеральных совокупностях) равны.
# 
# *Альтернативная гипотеза*: средние значения по всем группам (во всех генеральных совокупностях) не равны.
# 
# Чтобы не создавать искусственные группы на основе данных в *swiss*, загрузим таблицу с весами цыплят, которых кормили разным кормом :) Описание см. [здесь](https://vincentarelbundock.github.io/Rdatasets/doc/datasets/chickwts.html).

# In[7]:


dat = pd.read_csv('https://raw.githubusercontent.com/allatambov/Py-programming-3/master/add/chickwts.csv')


# In[8]:


dat.head()


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

# *Решение:*

# In[9]:


wgt = {}
for name, d in dat.groupby('feed'):
    wgt[name] = d.weight


# In[10]:


wgt


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

# In[11]:


st.f_oneway(wgt['casein'], wgt['horsebean'], wgt['linseed'], wgt['meatmeal'], wgt['soybean'], wgt['sunflower'])


# Функция возвращает наблюдаемое значение 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'])


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

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

# К сожалению, построение обычных регрессионных моделей, не связанных с машинным обучением, в Python, является нетривиальной задачей. Конечно, построить парную или множественную регрессию в Python не составит особенного труда, но вот получить готовую и красивую выдачу с результатами, как в R, сразу не получится (речь идет о библиотеках `scipy` и `sklearn`). Для того чтобы построить более сложную модель, потребуется больше сил, и сделать это без четкого понимания, как эта модель устроена, хорошо не получится, так как многие вещи придется дописывать самостоятельно.
# 
# Но мы начнем с самого простого, с парной линейной регрессии. В качестве зависимой переменной выберем показатель *Agriculture* (процент мужчин, занятых в сельском хозяйстве), а в качестве независимой ‒ *Catholic* (процент католического населения в кантоне). Построим модель с помощью функции `linregress`:

# In[24]:


st.linregress(df.Agriculture, df.Catholic)


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

# In[13]:


slope, intercept, rvalue, pvalue, stderr = st.linregress(df.Agriculture, df.Catholic)


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

# In[14]:


slope


# In[15]:


pvalue


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

# In[16]:


d = {'slope': slope, 'intercept': intercept, 'pvalue': pvalue, 'stderr': stderr}


# In[17]:


d


# In[18]:


res = pd.DataFrame.from_dict(d, orient='index')
res


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

# In[19]:


new = res.transpose()
new


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

# In[20]:


cols = ['intercept', 'slope', 'stderr', 'pvalue']
new = new[cols]
new


# Теперь хорошо бы добавить звездочки, которые обычно показывают, на каком уровне значимости коэффициент является значимым. В 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) # верно


# Теперь напишем функцию, которая будет "доклеивать" звездочки к коэффициенту при независимой переменной и которую можно будет применить к столбцу `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


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

# 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


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

# In[27]:


# без coef - он у нас текстовый
new[['intercept', 'slope', 'stderr', 'pvalue']].applymap("{0:.2f}".format)


# In[28]:


final = new[['intercept', 'slope', 'stderr', 'pvalue']].applymap("{0:.2f}".format)
final['coef'] = new['coef']


# In[29]:


final


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

# In[30]:


final.index = ['Agriculture ~ Catholic']


# In[31]:


final


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

# In[32]:


final.to_latex()


# In[33]:


print(final.to_latex()) # код можно просто скопировать


# Более подробно об опциях, которые можно учесть при выгрузке, можно узнать, вызвав `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)


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

# In[39]:


model.coef_


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

# In[40]:


model.intercept_


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

# In[41]:


model.predict(X)


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

# In[42]:


model.predict(X)-df.Agriculture


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

# In[43]:


sum((model.predict(X)-df.Agriculture)**2)


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

# In[44]:


model.residues_ # остатки


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