Проектное задание¶

import numpy as np
import pandas as pd
import scipy

import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns

import scipy.stats as stats
from statsmodels.stats.multitest import multipletests
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn import tree

Правила игры:¶

Вам разрешается пользоваться любыми источниками информации (Python Docs, StackOverflow, etc.), кроме ваших коллег.

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

Задание не должно совпадать с анализом, уже выложенным в Интернет.

Дедлайн:¶

Датасет надо выбрать на странице до 23:59:59 10 мая 2021

Проектное домашнее задание нужно отправить на почту python.msu@mail.ru до 23:59:59 22 мая 2021.

Тема письма должна быть "Фамилия Имя - Проект".

К письму нужно прикрепить эту Jupyter-тетрадь с названием "Family_Name_project.ipynb" и ваш файл с датасетом с тем же названием, как вы используете в коде. Файлы должны лежать в zip-архиве в папке "Family_Name". Пример:

Ivanov_Ivan.zip
|-- Ivanov_Ivan
    |-- Ivanov_Ivan_project.ipynb
    |-- titanic_train.csv
    |-- titanic_test.csv        # не нужен, если не используется в коде

Проектное задание (EDA). (до 30 баллов, минимум 10)¶

Творческое задание на навыки работы с библиотекой Pandas.

Что необходимо сделать:

Зайдите на сайт https://www.kaggle.com/ или другой ресурс и скачайте там какой-либо понравившийся вам датасет. Вы также можете взять собственные данные, если дадите их подробное описание. Отметьте выбранный вами датасет в таблице. Не выбирайте датасет, который уже кем-то выбран! Помните, что часто данные делятся на train и test, вам обычно нужна только часть train. Данные должны удовлетворять следующим условиям:
- В датасете присутствует как минимум 1000 записей,
- Для записей присутствует как минимум один категориальный признак (желательно два или три),
- Как минимум три числовых признака,
- Данные имеют понятный смысл (это не должны быть данные с БАК со странными буквенными кодами, понятные только специалистам).
Выбрав датасет, загрузите его и проведите EDA (exploratory data analysis). За проведение красивого анализа данных будут начислены дополнительные баллы.
Если вы не знаете, что делать, вот определенный минимум информации, которую нужно будет извлечь:
1. Посмотрите на параметры распределений числовых признаков - их средние, максимумы, минимумы, медианы, дисперсии и т.д. Постройте графики (подумайте, какие) наиболее интересных распределений. [до 3 баллов]
2. Посмотрите на наиболее крупные/мелкие/значительные/редкие/и т.п. записи по какому-либо из критериев (не увлекайтесь). [до 2 баллов]
3. Посмотрите на параметры этих распределений в группировке по категориальным переменным. Постройте графики этих распределений (какие?). [до 2 баллов]
4. Составьте наиболее осмысленную сводную таблицу по нескольким признакам. [до 2 баллов]
5. Покажите совместное распределение нескольких признаков (таблицей/графиками/и тем, и другим). [до 2 баллов]
6. Примените какой-либо разумный фильтр по одному из признаков в виде функции и проведите анализ отфильтрованных данных. [до 2 баллов]
7. Сформулируйте несколько разумных статистических гипотез и приведите свидетельства в их пользу/покажите недостаточность свидетельств в каждом случае. -- [до 10 баллов] (зависит от умения применять мат. аппарат и количества грамотно использованных тестов)
8. ** Задание для интересующихся: Для многих датасетов часто присутствует часть test. В этой части не указан один или несколько признаков (обычно категориальных или числовых). Попробуйте предсказать значения этого пропущенного признака по другой части датасета (train) с помощью математической модели. Помните, что выполнение этого задания потребует больших временных затрат! Библиотека моделей для обучения называется scikit-learn. [много баллов :^)]

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

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

EDA без графиков не принимаются.

Если вы решили выполнять пункт с ML, можете обращаться к преподавателям за советом по поводу выбора модели для ваших данных.

Начинаем здесь¶

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

Basic analysis¶

df = pd.read_csv("./penguins_size.csv", encoding = 'utf-8', sep = ',')
df = df.rename(columns={"species": "Species",
                   "island": "Island",
                   "culmen_length_mm": "Culmen length (mm)",
                   "culmen_depth_mm": "Culmen depth (mm)",
                   "flipper_length_mm": "Flipper length (mm)",
                   "body_mass_g": "Body mass (g)",
                   "sex": "Sex"})
df.head()

df.shape

(344, 7)

Всего в датафрейме 343 записи и 7 колонок.

Посчитаем число строчек с NaN и значением NA в столбце "пол", и если их немного, удалим их.

print(df.isna().sum())
print('\n')
print(df.loc[df['Sex'] == '.'])
#Всего 10 строчек - можно удалить
#Всего одна строчка с точкой в столбце "пол"

Species                 0
Island                  0
Culmen length (mm)      2
Culmen depth (mm)       2
Flipper length (mm)     2
Body mass (g)           2
Sex                    10
dtype: int64


    Species  Island  Culmen length (mm)  Culmen depth (mm)  \
336  Gentoo  Biscoe                44.5               15.7   

     Flipper length (mm)  Body mass (g) Sex  
336                217.0         4875.0   .

df = df[~pd.isna(df).any(axis=1)]
df = df.loc[~df['Sex'].isin(['.'])]
df

print(df.isna().sum())
print('\n')
print(df.loc[df['Sex'] == '.'])
#Убедились, что ничего не осталось, исчезло всего 11 строчек

Species                0
Island                 0
Culmen length (mm)     0
Culmen depth (mm)      0
Flipper length (mm)    0
Body mass (g)          0
Sex                    0
dtype: int64


Empty DataFrame
Columns: [Species, Island, Culmen length (mm), Culmen depth (mm), Flipper length (mm), Body mass (g), Sex]
Index: []

df.shape

(333, 7)

df.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 333 entries, 0 to 343
Data columns (total 7 columns):
 #   Column               Non-Null Count  Dtype  
---  ------               --------------  -----  
 0   Species              333 non-null    object 
 1   Island               333 non-null    object 
 2   Culmen length (mm)   333 non-null    float64
 3   Culmen depth (mm)    333 non-null    float64
 4   Flipper length (mm)  333 non-null    float64
 5   Body mass (g)        333 non-null    float64
 6   Sex                  333 non-null    object 
dtypes: float64(4), object(3)
memory usage: 20.8+ KB

Почищенный датасет состоит из 333 строчек и 7 столбцов Столбцы:

Species - вид
- Adelie - Пингвин Адели́ (Pygoscelis adeliae)
- Chinstrap - Антарктический пингвин (Pygoscelis antarctica)
- Gentoo - Субантарктический пингвин, или папуанский пингвин (Pygoscelis papua)
Island - остров Антарктики, на котором проживает пингвин
- Torgersen
- Biscoe
- Dream
Culmen length - длина клюва (мм)
Culmen depth - высота клюва (мм)
Flipper length - длина крыла (мм)
Body mass - масса тела (г)
Sex - пол пингвина

df.describe()

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

mn = df['Body mass (g)'].min()
mx = df['Body mass (g)'].max()
mean = df['Body mass (g)'].mean()
med = df['Body mass (g)'].median()
print(f'Минимальная масса тела - {round(mn)/1000} кг., максимальная - {round(mx)/1000} кг.;')
print(f'Средняя масса - {round(mean)/1000} кг., а медиана - {round(med)/1000}.кг.')

Минимальная масса тела - 2.7 кг., максимальная - 6.3 кг.;
Средняя масса - 4.207 кг., а медиана - 4.05.кг.

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

df.groupby(['Species','Sex']).agg({'Body mass (g)':'mean'})

Представим сводную таблицу средней массы тела в зависимости от пола и острова и визуализируем её в виде теплокарты. Видно, что более крупные особи живут на острове Biscoe.

plt.figure(figsize=(5,5))
sns.heatmap(df.pivot_table('Body mass (g)', index='Sex', columns='Island', aggfunc = 'mean'),
            cmap='coolwarm_r', annot=True)

<AxesSubplot:xlabel='Island', ylabel='Sex'>

Информация о цветовой палитре (представлена в виде HEX и RGB-percent кода), которую будем использовать далее:

Красный - f72585, (0.97,0.15,0.52)
Фиолетовый - 560bad, (0.34,0.04,0.68)
Синий - 4895ef, (0.28,0.58,0.94)

pal = ((0.97,0.15,0.52), (0.34,0.04,0.68), (0.28,0.58,0.94))

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

plt.rcParams['figure.figsize']=5,5
sns.displot(x='Culmen length (mm)', y='Culmen depth (mm)', hue='Species', kind="kde", data=df, palette=(pal))

<seaborn.axisgrid.FacetGrid at 0x11d5b1610>

Видим, что кластеризуются, хотя есть и пересечения.

Интересно, а есть ли корреляция между признаками? Построим матрицу корреляции и представим её в виде теплокарты.

plt.figure(figsize=(7,7))
correlation_matrix = df.corr()
sns.heatmap(correlation_matrix, annot=True, cmap="coolwarm_r", center=0, linewidths=0)

<AxesSubplot:>

Видим хорошую положительную корреляцию между длиной крыла и массой тела.

Совместное распределение массы тела по полу и виду¶

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

plt.figure(figsize=(7,7))
sns.boxplot(x='Species', y='Body mass (g)', hue='Sex', data=df, palette=pal)

<AxesSubplot:xlabel='Species', ylabel='Body mass (g)'>

Статистика¶

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

Создадим разнообразные переменные-фильтры.

adelie_mass = df[df['Species'] == 'Adelie']['Body mass (g)']
chinstrap_mass = df[df['Species'] == 'Chinstrap']['Body mass (g)']
gentoo_mass = df[df['Species'] == 'Gentoo']['Body mass (g)']
#
male_mass = df[df['Sex'] == 'MALE']['Body mass (g)']
female_mass = df[df['Sex'] == 'FEMALE']['Body mass (g)']

fvalue, pvalue = stats.f_oneway(adelie_mass, chinstrap_mass, gentoo_mass)
print(f'f-value is {round(fvalue)}, p-value is {pvalue}')

f-value is 342, p-value is 3.74450512630046e-81

Затем проведём поправку Холма на множественное тестирование.

p_val_1 = stats.f_oneway(adelie_mass, chinstrap_mass, gentoo_mass)[1]
a = multipletests(p_val_1, alpha=0.05, method='holm')
print(f'reject hypothesis - {bool(a[0])}, adjusted p-value is {float(a[1])}')

reject hypothesis - True, adjusted p-value is 3.74450512630046e-81

На примере видов разберём post hoc анализ сделаем ANOVA попарно для каждого пингвина и уже на трёх p-value сделаем поправку Холма.

fvalue1, pvalue1 = stats.f_oneway(adelie_mass, chinstrap_mass)
fvalu2, pvalue2 = stats.f_oneway(chinstrap_mass, gentoo_mass)
fvalue3, pvalue3 = stats.f_oneway(adelie_mass, gentoo_mass)
print(f'adj.p-value for Adelie-Chinstrap is {pvalue1},\nadj.p-value for Chinstrap-Gentoo is {pvalue2},\nadj.p-value for Adelie-Gentoo is {pvalue3}')

adj.p-value for Adelie-Chinstrap is 0.6748289682757607,
adj.p-value for Chinstrap-Gentoo is 2.7817455797857023e-46,
adj.p-value for Adelie-Gentoo is 1.8806652580951526e-66

a = multipletests((pvalue1, pvalue2, pvalue3), alpha=0.05, method='holm')
print(f'reject hypothesis - {list(a[0])},\nadjusted p-value is {list(a[1])}')

reject hypothesis - [False, True, True],
adjusted p-value is [0.6748289682757607, 5.563491159571405e-46, 5.641995774285457e-66]

Видим, что гипотеза H₀ о том, что между тремя видами есть хотя бы одна пара, у которой вес не различается статистически значимо - отвергается. То есть в хотя бы одной паре видов различие есть. При попарном сравнении и множественной поправке, видим, что две пары видов различаются значимо (Chinstrap-Gentoo и Adelie-Gentoo), а вот пара Adelie-Chinstrap не различается значимо. Это соотносится с графиком.

a = sns.displot(x='Body mass (g)', hue='Species', kde=True, data=df, palette=(pal), alpha=0.3,
            bins=40, element='step').add_legend()
a.fig.set_size_inches(10,5)

Проверим, есть ли статистически значимое различие между полами непараметрическим тестом Манна-Уитни и ANOVA.

fvalue, pvalue = stats.mannwhitneyu(male_mass, female_mass)
print(f'f-value is {round(fvalue)}, p-value is {pvalue}')

f-value is 6874, p-value is 9.066671516230528e-16

fvalue, pvalue = stats.f_oneway(male_mass, female_mass)
print(f'f-value is {round(fvalue)}, p-value is {pvalue}')

f-value is 73, p-value is 4.897246751596325e-16

По результатам двух тестов отвергаем гипотезу о неразличимости полов и тут.

Машинное обучение¶

Попробуем разобраться с основами машинного обучения на примере двух задач: предсказание пола с помощью логистической регрессии и предсказание вида с помощью дерева решений. Логистическая регрессия подходит для определения бинарной категориальной переменной по численным значениям. А дерево решений подойдёт, когда признаков больше двух (три в нашем случае). Для обучения будем давать массу пингвина, длину и высоту клюва, а также длину крыльев.

#Перепишем категориальные признаки в числа
df_learn = df
df_learn['Species'] = np.where(df_learn['Species'] == 'Adelie', '1', df_learn['Species'])
df_learn['Species'] = np.where(df_learn['Species'] == 'Chinstrap', '2', df_learn['Species'])
df_learn['Species'] = np.where(df_learn['Species'] == 'Gentoo', '3', df_learn['Species'])
#
df_learn['Island'] = np.where(df_learn['Island'] == 'Torgersen', '1', df_learn['Island'])
df_learn['Island'] = np.where(df_learn['Island'] == 'Biscoe', '2', df_learn['Island'])
df_learn['Island'] = np.where(df_learn['Island'] == 'Dream', '3', df_learn['Island'])
#
df_learn['Sex'] = np.where(df_learn['Sex'] == 'MALE', '1', df_learn['Sex'])
df_learn['Sex'] = np.where(df_learn['Sex'] == 'FEMALE', '2', df_learn['Sex'])
#
df_learn.head(2)

Логистическая регрессия для предсказания пола пингвина¶

Логистическая регрессия - аналог одного нейрона в нейросети, который настраивает веса для определения бинарного признака.

np.random.seed(999)
tr, tst  = train_test_split(df_learn, shuffle=True)
x_tr = tr.drop(['Species', 'Island', 'Sex'], axis=1)
y_tr = tr['Sex']
x_tst = tst.drop(['Species', 'Island', 'Sex'], axis=1)
y_tst = tst['Sex']
x_tr.head(1)

model = LogisticRegression(solver='liblinear', random_state=0).fit(x_tr, y_tr)
print(f'Исследуются бинарные классы: {model.classes_}: 1 - самцы и 2 - самки')
print(f'Коэффициенты регрессии для 4 параметров: {model.coef_}')
print(f'Пересечение с осью OX: {model.intercept_}')

Исследуются бинарные классы: ['1' '2']: 1 - самцы и 2 - самки
Коэффициенты регрессии для 4 параметров: [[-0.13028095 -0.77438928  0.20553451 -0.00543067]]
Пересечение с осью OX: [0.3855365]

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

model.score(x_tst, y_tst)

0.8571428571428571

pred = np.array(model.predict(x_tst), dtype=int)
test = np.array(y_tst, dtype=int)
dif = pred - test
print(f'Предсказанный массив отличается от истинного на {np.count_nonzero(dif)} элементов из {len(dif)}.')

Предсказанный массив отличается от истинного на 12 элементов из 84.

Мы получили модель, которая предсказывает пол пингвина в 85% случаев.

Дерево решений для предсказания вида пингвина¶

np.random.seed(999)
tr, tst = train_test_split(df_learn, shuffle=True)
x_tr = tr.drop(['Species', 'Island', 'Sex'], axis=1)
y_tr = tr['Species']
x_tst = tst.drop(['Species', 'Island', 'Sex'], axis=1)
y_tst = tst['Species']

model = tree.DecisionTreeClassifier().fit(x_tr, y_tr)
model.score(x_tst, y_tst)

0.9761904761904762

pred = np.array(model.predict(x_tst), dtype=int)
test = np.array(y_tst, dtype=int)
dif = pred - test
print(f'Предсказанный массив отличается от истинного на {np.count_nonzero(dif)} элемента из {len(dif)}.')

Предсказанный массив отличается от истинного на 2 элемента из 84.

А вот дерево решений предсказывает вид гораздо лучше, чем логистическая регрессия пол. Оно работает в 97% случаев.

Послесловие¶

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

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

	Culmen length (mm)	Culmen depth (mm)	Flipper length (mm)	Body mass (g)
count	333.000000	333.000000	333.000000	333.000000
mean	43.992793	17.164865	200.966967	4207.057057
std	5.468668	1.969235	14.015765	805.215802
min	32.100000	13.100000	172.000000	2700.000000
25%	39.500000	15.600000	190.000000	3550.000000
50%	44.500000	17.300000	197.000000	4050.000000
75%	48.600000	18.700000	213.000000	4775.000000
max	59.600000	21.500000	231.000000	6300.000000

	Species	Island	Culmen length (mm)	Culmen depth (mm)	Flipper length (mm)	Body mass (g)	Sex
0	1	1	39.1	18.7	181.0	3750.0	1
1	1	1	39.5	17.4	186.0	3800.0	2

	Species	Island	Culmen length (mm)	Culmen depth (mm)	Flipper length (mm)	Body mass (g)	Sex
0	Adelie	Torgersen	39.1	18.7	181.0	3750.0	MALE
1	Adelie	Torgersen	39.5	17.4	186.0	3800.0	FEMALE
2	Adelie	Torgersen	40.3	18.0	195.0	3250.0	FEMALE
3	Adelie	Torgersen	NaN	NaN	NaN	NaN	NaN
4	Adelie	Torgersen	36.7	19.3	193.0	3450.0	FEMALE