from __future__ import division

import numpy as np
import pandas as pd

from scipy import stats
import statsmodels.formula.api as smf
import statsmodels.stats.api as sms

%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns

from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

"""
Explanation: Тест. Практика построения регрессии
End of explanation
"""


botswana = pd.read_csv('botswana.tsv', delimiter='\t')
botswana.info()

botswana.columns
botswana.head()

print(botswana.religion.value_counts())

"""
Explanation: Давайте проанализируем данные опроса 4361 женщин из Ботсваны:
botswana.tsv
О каждой из них мы знаем:
сколько детей она родила (признак ceb)
возраст (age)
длительность получения образования (educ)
религиозная принадлежность (religion)
идеальное, по её мнению, количество детей в семье (idlnchld)
была ли она когда-нибудь замужем (evermarr)
возраст первого замужества (agefm)
длительность получения образования мужем (heduc)
знает ли она о методах контрацепции (knowmeth)
использует ли она методы контрацепции (usemeth)
живёт ли она в городе (urban)
есть ли у неё электричество, радио, телевизор и велосипед (electric, radio, tv, bicycle)
Давайте научимся оценивать количество детей ceb по остальным признакам.
Загрузите данные и внимательно изучите их. Сколько разных значений принимает признак religion?
End of explanation
"""


bw_wo_nan = botswana.dropna()
bw_wo_nan.shape

"""
Explanation: Во многих признаках есть пропущенные значения. Сколько объектов из 4361 останется, если выбросить все, содержащие пропуски?
End of explanation
"""


botswana['nevermarr'] = [1 if botswana.loc[i, 'evermarr'] == 0 else 0 for i in range(botswana.shape[0])]
botswana.head()

np.unique(botswana.agefm[botswana.agefm.notnull()].values)

botswana.agefm[botswana.agefm.isnull()] = 0
botswana.head()

del botswana['evermarr']

botswana.heduc[botswana.heduc.isnull() & botswana.nevermarr.values == 1] = -1
botswana.head()

botswana.info()

botswana.heduc.isnull().value_counts()

"""
Explanation: В разных признаках пропуски возникают по разным причинам и должны обрабатываться по-разному.
Например, в признаке agefm пропуски стоят только там, где evermarr=0, то есть, они соответствуют женщинам, никогда не выходившим замуж. Таким образом, для этого признака NaN соответствует значению "не применимо".
В подобных случаях, когда признак x1 на части объектов в принципе не может принимать никакие значения, рекомендуется поступать так:
создать новый бинарный признак
x2={1,0,x1='не применимо',иначе;
заменить "не применимо" в x1 на произвольную константу c, которая среди других значений x1 не встречается.
Теперь, когда мы построим регрессию на оба признака и получим модель вида
y=β0+β1x1+β2x2,
на тех объектах, где x1 было измерено, регрессионное уравнение примет вид
y=β0+β1x,
а там, где x1 было "не применимо", получится
y=β0+β1c+β2.
Выбор c влияет только на значение и интерпретацию β2, но не β1.
Давайте используем этот метод для обработки пропусков в agefm и heduc.
Создайте признак nevermarr, равный единице там, где в agefm пропуски.
Удалите признак evermarr — в сумме с nevermarr он даёт константу, значит, в нашей матрице X будет мультиколлинеарность.
Замените NaN в признаке agefm на cagefm=0.
У объектов, где nevermarr = 1, замените NaN в признаке heduc на cheduc1=−1 (ноль использовать нельзя, так как он уже встречается у некоторых объектов выборки).
Сколько осталось пропущенных значений в признаке heduc?
End of explanation
"""


botswana['idlnchld_noans'] = 0
botswana.loc[botswana.idlnchld.isnull(), 'idlnchld_noans'] = 1

botswana['heduc_noans'] = 0
botswana.loc[botswana.heduc.isnull(), 'heduc_noans'] = 1

botswana['usemeth_noans'] = 0
botswana.loc[botswana.usemeth.isnull(), 'usemeth_noans'] = 1

botswana.head()

botswana.idlnchld[botswana.idlnchld.isnull()] = -1
botswana.heduc[botswana.heduc.isnull()] = -2
botswana.usemeth[botswana.usemeth.isnull()] = -1

botswana.info()

botswana = botswana.dropna()

botswana.info()

elem_num = botswana.shape[0] * botswana.shape[1]
print('Array size: %d' % elem_num)

"""
Explanation: Избавимся от оставшихся пропусков.
Для признаков idlnchld, heduc и usemeth проведите операцию, аналогичную предыдущей: создайте индикаторы пропусков по этим признакам (idlnchld_noans, heduc_noans, usemeth_noans), замените пропуски на нехарактерные значения (cidlnchld=−1, cheduc2=−2 (значение -1 мы уже использовали), cusemeth=−1).
Остались только пропуски в признаках knowmeth, electric, radio, tv и bicycle. Их очень мало, так что удалите объекты, на которых их значения пропущены.
Какого размера теперь наша матрица данных? Умножьте количество строк на количество всех столбцов (включая отклик ceb).
End of explanation
"""


botswana.columns

formula = 'ceb ~ ' + ' + '.join(botswana.columns[1:])
formula

reg_m = smf.ols(formula, data=botswana)
fitted_m = reg_m.fit()
print(fitted_m.summary())

print(botswana.religion.value_counts())

"""
Explanation: Постройте регрессию количества детей ceb на все имеющиеся признаки методом smf.ols, как в разобранном до этого примере. Какой получился коэффициент детерминации R2? Округлите до трёх знаков после десятичной точки.
Если код из примера у вас не воспроизводится:
убедитесь, что вы сделали так:
import statsmodels.formula.api as smf
возможно, вам нужно обновить библиотеку patsy; выполните в командной строке
pip install -U patsy
End of explanation
"""


print('Breusch-Pagan test: p=%f' % sms.het_breuschpagan(fitted_m.resid, fitted_m.model.exog)[1])

reg_m2 = smf.ols(formula, data=botswana)
fitted_m2 = reg_m2.fit(cov_type='HC1')
print(fitted_m2.summary())

"""
Explanation: Проверьте критерием Бройша-Пагана гомоскедастичность ошибки в построенной модели. Выполняется ли она?
Если ошибка гетероскедастична, перенастройте модель, сделав поправку Уайта типа HC1.
End of explanation
"""


formula2 = 'ceb ~ age + educ + idlnchld + knowmeth + usemeth + agefm + heduc + urban + electric + bicycle \
+ nevermarr + idlnchld_noans + heduc_noans + usemeth_noans'
formula2

reg_m3 = smf.ols(formula2, data=botswana)
fitted_m3 = reg_m3.fit()
print(fitted_m3.summary())

print('Breusch-Pagan test: p=%f' % sms.het_breuschpagan(fitted_m3.resid, fitted_m3.model.exog)[1])

reg_m4 = smf.ols(formula2, data=botswana)
fitted_m4 = reg_m4.fit(cov_type='HC1')
print(fitted_m4.summary())

print('F=%f, p=%f, k1=%f' % reg_m2.fit().compare_f_test(reg_m4.fit()))

"""
Explanation: Удалите из модели незначимые признаки religion, radio и tv. Проверьте гомоскедастичность ошибки, при необходимости сделайте поправку Уайта.
Не произошло ли значимого ухудшения модели после удаления этой группы признаков? Проверьте с помощью критерия Фишера. Чему равен его достигаемый уровень значимости? Округлите до четырёх цифр после десятичной точки.
Если достигаемый уровень значимости получился маленький, верните все удалённые признаки; если он достаточно велик, оставьте модель без религии, тв и радио.
End of explanation
"""


formula3 = 'ceb ~ age + educ + idlnchld + knowmeth + agefm + heduc + urban + electric + bicycle \
+ nevermarr + idlnchld_noans + heduc_noans'
formula3

reg_m5 = smf.ols(formula3, data=botswana)
fitted_m5 = reg_m5.fit()
print(fitted_m5.summary())

print('F=%f, p=%.40f, k1=%f' % reg_m4.fit().compare_f_test(reg_m5.fit()))

"""
Explanation: Признак usemeth_noans значим по критерию Стьюдента, то есть, при его удалении модель значимо ухудшится. Но вообще-то отдельно его удалять нельзя: из-за того, что мы перекодировали пропуски в usemeth произвольно выбранным значением cusemeth=−1, удалять usemeth_noans и usemeth можно только вместе.
Удалите из текущей модели usemeth_noans и usemeth. Проверьте критерием Фишера гипотезу о том, что качество модели не ухудшилось. Введите номер первой значащей цифры в достигаемом уровне значимости (например, если вы получили 5.5×10−8, нужно ввести 8).
Если достигаемый уровень значимости получился маленький, верните удалённые признаки; если он достаточно велик, оставьте модель без usemeth и usemeth_noans.
End of explanation
"""


print(fitted_m4.summary())

"""
Explanation: Посмотрите на доверительные интервалы для коэффициентов итоговой модели (не забудьте использовать поправку Уайта, если есть гетероскедастичность ошибки) и выберите правильные выводы.
End of explanation
"""


