Поиск причинно-следственных связей с помощью Python

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

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

Исследования показывают, что качество питьевой воды напрямую связано с уровнем смертности в населении. По данным Роспотребнадзора, в 2023 году с ненормативным качеством водоснабжения было ассоциировано 7,18 дополнительных случаев смертности на 100 тыс. населения. Это подчеркивает необходимость применения аналитических методов для изучения взаимосвязей между качеством воды и здоровьем людей, что позволит не только выявить проблемы, но и разработать эффективные меры по их устранению.

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

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

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

Начнем анализ с импорта библиотек и загрузки датасета.

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import plotly.express as px
import plotly.graph_objects as go
df = pd.read_csv('water.csv')
df.head()

Датасет со статистикой по городам Англии

Рис. 1: Датасет со статистикой по городам Англии

Оценим размеры датасета:

print('В датасете строк/колонок:', df.shape)

В датасете строк/колонок: (61, 5)

Посмотрим какие типы данных в столбцах.

df.info()

Информация о типах столбцах в датасете

Рис. 2: Информация о типах столбцах в датасете

Видим что в таблице типы данных оформлены корректно. Числовые столбцы – int64, строковые – object. Все 61 строк содержат данные.

В нашем датасете содержатся данные о смертности людей и данные о жесткости воды по городам Англии. Давайте взглянем на базовые статистики этих показателей.

df['mortality'].describe()

count 61.000000
mean 1524.147541
std 187.668754
min 1096.000000
25% 1379.000000
50% 1555.000000
75% 1668.000000
max 1987.000000
Name: mortality, dtype: float64

Мы видим что средняя смертность в английских городах составляет 1524 на 100000 населения. Медианное значение – 1555. Что близко к среднему. Что может говорить о том, что разброс значений небольшой.

Это же подтверждают значения минимума (1096) и максимума (1987). Разброс между ними составляет:

df.mortality.max() - df.mortality.min()

891

Построим гистограмму распределения этих значений.

df.mortality.hist(bins=10)

Рис. 3: Гистограмма смертности на 100000 чел.

Рис. 3: Гистограмма смертности на 100000 чел.

Видим, что основная часть сконцентрирована вокруг 1600. Мы также видим резкий спад графика в районе 1800. Очевидно что это самые неблагополучные города Англии в плане смертности.

Давайте взглянем на статистики жесткости воды.

df['hardness'].describe()

count 61.000000
mean 47.180328
std 38.093966
min 5.000000
25% 14.000000
50% 39.000000
75% 75.000000
max 138.000000
Name: hardness, dtype: float64

Средняя жесткость воды (концентрация кальция в питьевой воде) составляет 47.18, медианное значение – 39.0. Здесь разброс данных уже выше, минимум – 5, максимум – 138.

df.hardness.hist(bins=10)

Гистограмма распределения жесткости воды по городам

Рис. 4: Гистограмма распределения жесткости воды по городам

Судя по графику распределения в большинстве городов Англии вода не содержит много кальция. Однако мы также можем заметить и 2-ую моду в диапазоне 58-83. Очевидно, что у воды в этих городах есть какая-то общая черта, возможно одна и та же река.

Давайте теперь приступим к поиску причинно-следственных связей и проанализируем есть ли взаимосвязь между жесткостью питьевой воды и смертностью.

figure = px.scatter(data_frame=df,
                    x='hardness',
                    y='mortality',
                    color='mortality', 
                    trendline='ols', 
                    title='Взаимосвязь между жесткостью питьевой воды и смертностью')
figure.update_traces(marker={'size': 12})
figure.show()

Взаимосвязь между жесткостью питьевой воды и смертностью

Рис. 5: Взаимосвязь между жесткостью питьевой воды и смертностью

df['hardness'].corr(df['mortality']) 

-0.6548486232042463

-0.65 – это умеренно сильная обратная корреляция. Давайте также рассчитаем коэффициент корреляции Спирмена.

df['hardness'].corr(df['mortality'], method='spearman')

-0.6316646189166502

Получаем почти идентичное значение. Вывод: жесткость воды имеет обратную взаимосвязь со смертностью людей.

Давайте теперь посмотрим какие города Англии имеют высокую смертность.

df.sort_values(by='mortality', ascending=False).head(5)

Города Англии с самой высокой смертностью

Рис. 6: Города Англии с самой высокой смертностью

Топ-5 анти-лидеров: Salford, Manchester, Bootle, Blackburn, Liverpool. Все находятся на севере страны. Жесткость воды колеблется от 8 до 15.

Теперь давайте взглянем на города с самой низкой смертностью.

df.sort_values(by='mortality', ascending=True).head(5)

Города Англии с самой низкой смертностью

Рис. 7: Города Англии с самой низкой смертностью

Топ-5 лидеров: Ipswich, Oxford, Reading, Bath, Croydon. Все находятся на юге Англии. У всех высокая жесткость питьевой воды (от 96 до 138).

Давайте также взглянем таким же образом на топы по жесткости.

df.sort_values(by='hardness', ascending=False).head(5)

Топ-5 городов с самой жесткой водой

Рис. 8: Топ-5 городов с самой жесткой водой

Видим снова в топе южные города. Смертность во всех ниже среднего (1524).

А теперь давайте посмотрим на города с самой мягкой питьевой водой.

df.sort_values(by='hardness', ascending=True).head(5)

Топ-5 городов с самой мягкой водой

Рис. 9: Топ-5 городов с самой мягкой водой

Здесь уже видим что в топ попали не только северные, но и южные города. Также у этих городов есть еще одна общая черта – смертность около (на юге) или выше среднего (на севере страны).

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

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

sns.pairplot(df[['hardness', 'mortality']],
             kind='reg')

Диаграмма парных сравнений жесткости воды и смертности

Рис. 10: Диаграмма парных сравнений жесткости воды и смертности

Давайте углубимся в наши исследования. Попробуем построить модель линейной регрессии и посмотреть на ее значения.

import statsmodels.api as sm
from statsmodels.sandbox.regression.predstd import wls_prediction_std
from sklearn.model_selection import train_test_split
X = df[['hardness']]
y = df['mortality']

Разделим данные на обучаемую и тестовую выборки.

X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                    test_size=0.2)
X_const = sm.add_constant(X_train)

И посмотрим какие получились размеры выборок.

print('X_train:', X_train.shape)
print('X_test:', X_test.shape)
print('y_train:', y_train.shape)
print('y_test:', y_test.shape)

X_train: (48, 1)
X_test: (13, 1)
y_train: (48,)
y_test: (13,)

Теперь укажем на тип модели и посмотрим на результаты регрессии.

model = sm.OLS(y_train, X_const)
results = model.fit()
print(results.summary())

Таблица статистик после обучения линейной регрессии

Рис. 11: Таблица статистик после обучения линейной регрессии

Итак, мы получили модель регрессии для жесткости воды и смертности. В таблице выше обратите внимание на показатели R-squared и Adj. R-squared. Это коэффициент детерминации и взвешенный КД. Оба имеют значения больше 0.45, что можно интерпретировать как зависимость данных на среднем уровне. Таким образом мы можем сделать вывод что эти 2 величины имеют некоторую взаимосвязь.

Давайте теперь построим график линейной регрессии и посмотрим какой получился доверительный интервал.

prstd, iv_l, iv_u = wls_prediction_std(results) 

fig, ax = plt.subplots(figsize=(8,6))
ax.plot(X_const.iloc[:, 1], y_train, 'o', label='data')
ax.plot(X_const.iloc[:, 1], results.fittedvalues, 'r--.', label='OLS')
ax.plot(X_const.iloc[:, 1], iv_u, 'g--')
ax.plot(X_const.iloc[:, 1], iv_l, 'g--')
ax.legend(loc='best')

График линейной регрессии с доверительными интервалами

Рис. 12: График линейной регрессии с доверительными интервалами

Проанализируем как распределились остатки регрессионной модели.

plt.scatter(x=X_const.iloc[:, 1], y=results.resid)
plt.xlabel('Прогноз')
plt.ylabel('Остатки')

График остатков модели линейной регрессии

Рис. 13: График остатков модели линейной регрессии

Видим что в расположении остатков какой-либо системности нет. Это значит что наша регрессионная модель построена достаточно хорошо.

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

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

north = df[df['location'] == 'North']
south = df[df['location'] == 'South']

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

north[['mortality', 'hardness']].describe()

Описательная статистика столбцов по городам на севере Англии

Рис. 14: Описательная статистика столбцов по городам на севере Англии

south[['mortality', 'hardness']].describe()

Описательная статистика столбцов по городам на юге Англии

Рис. 15: Описательная статистика столбцов по городам на юге Англии

Мы видим что в южных городах в среднем меньше смертность (1376 против 1633), при том что вода в кранах на юге в среднем жестче (69.76 против 30.40).

Давайте посмотрим на графики северных городов.

sns.pairplot(north[['hardness', 'mortality']],
             kind='reg')

Диаграмма парных сравнений показателей столбцов северных городов

Рис. 16: Диаграмма парных сравнений показателей столбцов северных городов

А теперь южных городов.

sns.pairplot(south[['hardness', 'mortality']],
             kind='reg')

Диаграмма парных сравнений показателей столбцов южных городов

Рис. 17: Диаграмма парных сравнений показателей столбцов южных городов

Мы видим что на севере Англии в кранах течет в основном мягкая вода, а на юге – жесткая.

Мы также видим что отличаются графики смертности. У юга Англии он заканчивается на отметке 1600-1650, в то время как на севере эти значения являются вершиной распределения с наиболее высокой частотой и идут дальше.

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

Для наглядности наложим графики друг на друга.

north['mortality'].plot(kind='hist', 
                        alpha=0.5,
                        bins=15,
                        label='Север Англии',
                        density=True)
south['mortality'].plot(kind='hist', 
                        alpha=0.5,
                        bins=15,
                        label='Юг Англии',
                        density=True)
plt.legend(loc='upper left')
plt.title('Сравнение смертности на 100000 южных и северных городов')

Графики сравнения смертности на 100000 жителей южных и северных городов

Рис. 18: Графики сравнения смертности на 100000 жителей южных и северных городов

Таким же образом оценим распределение жесткости воды.

north['hardness'].plot(kind='hist', 
                       alpha=0.5,
                       bins=12,
                       label='Север Англии',
                       density=True)
south['hardness'].plot(kind='hist', 
                       alpha=0.5,
                       bins=12,
                       label='Юг Англии',
                       density=True)
plt.legend(loc='upper left')
plt.title('Сравнение жесткости воды южных и северных городов')

Графики сравнения жесткости воды в южных и северных городах Англии

Рис. 19: Графики сравнения жесткости воды в южных и северных городах Англии

Давайте теперь сравним коэффициенты корреляций севера и юга Англии.

print('Север, Пирсон: ', north['hardness'].corr(north['mortality']))
print('Юг, Пирсон:    ', south['hardness'].corr(south['mortality']))
print('Север, Спирмен:', north['hardness'].corr(north['mortality'], method='spearman'))
print('Юг, Спирмен:   ', south['hardness'].corr(south['mortality'], method='spearman'))

Север, Пирсон: -0.3685978383288718
Юг, Пирсон: -0.6021532715484156
Север, Спирмен: -0.4042078956511175
Юг, Спирмен: -0.5957229185013566

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

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

Давайте также посмотрим как изменится модель линейной регрессии если рассматривать юг и север Англии по отдельности.

Начнем традиционно с севера.

X = north[['hardness']]
y = north['mortality']
X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                    test_size=0.2)
X_const = sm.add_constant(X_train)

Теперь укажем на тип модели и посмотрим на результаты регрессии.

model = sm.OLS(y_train, X_const)
results = model.fit()
print(results.summary())

Таблица результатов построения регрессии по северным городам

Рис. 20: Таблица результатов построения регрессии по северным городам

Видим что регрессия получилась намного хуже чем в общем случае, коэффициент детерминации упал с 0.45 до 0.23.

print('Parameters: ', results.params)
print('R2: ', results.rsquared)

Parameters: const 1689.945775
hardness -2.676408
dtype: float64
R2: 0.2366580947844319

Давайте теперь построим график линейной регрессии и посмотрим какой получился доверительный интервал.

prstd, iv_l, iv_u = wls_prediction_std(results) 

fig, ax = plt.subplots(figsize=(8,6))
ax.plot(X_const.iloc[:, 1], y_train, 'o', label='data')
ax.plot(X_const.iloc[:, 1], results.fittedvalues, 'r--.', label='OLS')
ax.plot(X_const.iloc[:, 1], iv_u, 'g--')
ax.plot(X_const.iloc[:, 1], iv_l, 'g--')
ax.legend(loc='best')

График линейной регрессии с доверительными интервалами для северных городов

Рис. 20: График линейной регрессии с доверительными интервалами для северных городов

Проанализируем как распределились остатки регрессионной модели.

plt.scatter(x=X_const.iloc[:, 1], y=results.resid)
plt.xlabel('Прогноз')
plt.ylabel('Остатки')

График остатков регрессии для северных городов

Рис. 21: График остатков регрессии для северных городов

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

Давайте проделаем аналогичный анализ с южными городами.

X = south[['hardness']]
y = south['mortality']
X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                    test_size=0.2)
X_const = sm.add_constant(X_train)

Теперь укажем на тип модели и посмотрим на результаты регрессии.

model = sm.OLS(y_train, X_const)
results = model.fit()
print(results.summary())

Таблица результатов построения регрессии по южным городам

Рис. 22: Таблица результатов построения регрессии по южным городам

Видим что регрессия получилась намного хуже чем в общем случае, коэффициент детерминации упал с 0.45 до 0.27.

print('Parameters: ', results.params)
print('R2: ', results.rsquared)

Parameters: const 1523.292086
hardness -1.842919
dtype: float64
R2: 0.2735659054615863

Давайте теперь построим график линейной регрессии и посмотрим какой получился доверительный интервал.

prstd, iv_l, iv_u = wls_prediction_std(results) 

fig, ax = plt.subplots(figsize=(8,6))
ax.plot(X_const.iloc[:, 1], y_train, 'o', label='data')
ax.plot(X_const.iloc[:, 1], results.fittedvalues, 'r--.', label='OLS')
ax.plot(X_const.iloc[:, 1], iv_u, 'g--')
ax.plot(X_const.iloc[:, 1], iv_l, 'g--')
ax.legend(loc='best')

График линейной регрессии с доверительными интервалами для южных городов

Рис. 23: График линейной регрессии с доверительными интервалами для южных городов

Проанализируем как распределились остатки регрессионной модели.

plt.scatter(x=X_const.iloc[:, 1], y=results.resid)
plt.xlabel('Прогноз')
plt.ylabel('Остатки')

График остатков регрессии для южных городов

Рис. 24: График остатков регрессии для южных городов

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

plt.hist(results.resid, bins=9)

Диаграмма распределения остатков

Рис. 25: Диаграмма распределения остатков

Однако как показывает диаграмма выше о нормальности распределения остатков говорить не приходится.

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

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