import numpy as np
import pandas as pd
from sklearn.metrics import r2_score, root_mean_squared_error
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNet
from sklearn.preprocessing import StandardScaler
from sklearn.impute import SimpleImputer

import matplotlib.pyplot as plt
import seaborn as sns

sns.set_theme("notebook", font_scale=1.4, palette="Set2")

import warnings

warnings.simplefilter(action="ignore")

random_state = 42

df = pd.read_csv("WDICSV.csv")
print(f"Размер таблицы: {df.shape}")
df.head(3)

Размер таблицы: (397936, 68)

countries = pd.Series(
    df["Country Name"].values, index=df["Country Code"], name="Country Name"
)
countries = countries[~countries.index.duplicated(keep="first")]
countries.head(3)

indicators = pd.Series(
    df["Indicator Name"].values, index=df["Indicator Code"], name="Indicator Name"
)
indicators = indicators[~indicators.index.duplicated(keep="first")]
indicators.sample(3, random_state=random_state)

df_groupped = df.groupby("Country Name").count()
df_groupped.sample(7, random_state=random_state)

plt.figure(figsize=(15, 6))

# Обычно смотрят не на количество записей, а на количество пропущенных ячеек
sns.heatmap(
    len(df["Indicator Name"].unique()) - df_groupped.iloc[:, 3:],
    cbar_kws={"label": "Количество пропусков"},
)
plt.title("Распределение пропусков в данных");

df22 = df.pivot_table(
    index="Country Code", columns="Indicator Code", values="2022", aggfunc="first"
)
df22.sample(5, random_state=random_state)

not_countries = ['AFE', 'ARB', 'CEB', 'EAP', 'EAR', 'EAS', 'ECA', 'ECS', 'EMU', 'EUU', 'FCS', 'HIC', 'HPC', 'IBD', 'IBT', 'IDA', 'IDB', 'IDX', 'LAC', 'LCN', 'LDC', 'LIC', 'LMC', 'LMY', 'LTE', 'MEA', 'MNA', 'NAC', 'OED', 'OSS', 'PRE', 'PST', 'SSA', 'SSF', 'SST', 'TEA', 'TEC', 'TLA', 'TMN', 'TSA', 'TSS', 'UMC', 'WLD']

print("Всего аггрегированных регионов:", len(not_countries))
countries[not_countries].sample(5, random_state=random_state)

Всего аггрегированных регионов: 43

target_like_features = df22.filter(like="CO2").columns.to_list() + [
    "EN.GHG.TOT.ZG.AR5",
    "EN.GHG.ALL.MT.CE.AR5",
    "EN.GHG.ALL.PC.CE.AR5",
    "EN.GHG.ALL.LU.MT.CE.AR5",
]
y = df22["EN.GHG.CO2.PC.CE.AR5"].dropna()
X = df22.drop(columns=target_like_features, errors="ignore").loc[y.index]
X.shape, y.shape

((251, 1047), (251,))

nans = X.isna().sum()

plt.figure(figsize=(8, 4))
nans.hist()
plt.xlabel("Количество пропусков")
plt.ylabel("Число признаков в бине")
plt.title("Распределение пропусков по индикаторам");

threshold = 50  # порог
X.drop(columns=X.columns[nans > threshold], inplace=True)
print(
    f"Осталось {len(X.columns)} признаков с менее чем {threshold} пропущенными значениями"
)

Осталось 301 признаков с менее чем 50 пропущенными значениями

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.4, random_state=random_state
)
X_train.shape, X_test.shape, y_train.shape, y_test.shape

((150, 301), (101, 301), (150,), (101,))

# Инициализируем класс
scaler = StandardScaler().set_output(transform="pandas")

# Преобразуем обучающую и валидационную выборки
X_train = scaler.fit_transform(X_train)
X_test = scaler.transform(X_test)

imputer = SimpleImputer(strategy="median").set_output(transform="pandas")

X_train = imputer.fit_transform(X_train)
X_test = imputer.transform(X_test)

np.any(X_train.isna()), np.any(X_test.isna()),

(False, False)

# Сюда будем записывать R2
r2 = pd.Series(index=X_train.columns, dtype=float)

for feature in X_train.columns:

    model_1d = LinearRegression()
    model_1d.fit(X_train[feature].values.reshape(-1, 1), y_train)
    r2[feature] = r2_score(
        y_test, model_1d.predict(X_test[feature].values.reshape(-1, 1))
    )

plt.figure(figsize=(8, 3))
sns.histplot(r2, bins=10, binrange=(-1, 1))
plt.title("Распределение $R^2$ для моделей по отдельным фичам")
plt.xlabel("$R^2$")
plt.yscale("log")

# Порог установим с небольшим запасом
threshold = 0.1
features_to_exclude = X.columns[r2 < threshold]

X_train.drop(columns=features_to_exclude, inplace=True, errors="ignore")
X_test.drop(columns=features_to_exclude, inplace=True, errors="ignore")
X_train.shape, X_test.shape

((150, 68), (101, 68))

# Инициализация модели и обучение
linear_model = LinearRegression()
linear_model.fit(X_train, y_train)

# Предсказание на тестовой выборке
y_pred = linear_model.predict(X_test)

print(
    f"R2: {r2_score(y_test, y_pred):.2f} \nRMSE: {root_mean_squared_error(y_test, y_pred):.2f}"
)

R2: -3.70 
RMSE: 9.45

# Предсказание средним значением таргета на обучающей выборке
y_pred_mean = y_train.mean() * np.ones_like(y_test)

print(
    f"R2: {r2_score(y_test, y_pred_mean):.2f} \nRMSE: {root_mean_squared_error(y_test, y_pred_mean):.2f}"
)

R2: -0.07 
RMSE: 4.50

plt.figure(figsize=(8, 4), tight_layout=True)

err = y_test.values - y_pred
std = y_test.values - y_pred_mean

sns.kdeplot(err, fill=True, alpha=0.35, label="Ошибки модели")
sns.kdeplot(std, fill=True, alpha=0.35, label="Отклонение от\nсреднего на трейне")
plt.legend()
plt.xlabel("Величина ошибки");

np.linalg.cond(X_train)

2.0233202454319468e+16

fig, (ax_low_corr, ax_high_corr) = plt.subplots(
    1, 2, figsize=(15, 6), sharey=True, tight_layout=True
)

# Выберем некоторый индикатор для сравнения
compared = "NY.GDP.PCAP.PP.CD"

# Два других индикатора
high_corr = "NY.GNP.PCAP.PP.CD"  # Сильно коррелирует с первым
low_corr = "BX.TRF.PWKR.DT.GD.ZS"  # Слабо коррелирует с первым

# Объединим их в один массив и удалим выбросы
sub_X = X[[compared, high_corr, low_corr]].dropna()

# Отображаем зависимость между значениями двух слабо коррелирующих индикаторов
sns.regplot(
    x=sub_X[low_corr],
    y=sub_X[compared],
    scatter_kws={"s": 100},
    line_kws={"lw": 5},
    ax=ax_low_corr,
)

# Настраиваем подписи
ax_low_corr.set_ylabel(indicators[compared])
ax_low_corr.set_xlabel(indicators[low_corr])

# Вычисляем коэффициент корреляции для слабо коррелирующих индикаторов
weak_corr = np.corrcoef(sub_X[low_corr], sub_X[compared])[0, 1]
ax_low_corr.set_title(f"Слабая корреляция r = {weak_corr:.2f}")  # Заголовок графика

# Отображаем зависимость между значениями двух сильно коррелирующих индикаторов
sns.regplot(
    x=sub_X[high_corr],
    y=sub_X[compared],
    scatter_kws={"s": 100},
    line_kws={"lw": 5},
    ax=ax_high_corr,
)

# Настраиваем подписи
ax_high_corr.set_ylabel("")
ax_high_corr.set_xlabel(indicators[high_corr])

# Вычисляем коэффициент корреляции для сильно коррелирующих индикаторов
strong_corr = np.corrcoef(sub_X[high_corr], sub_X[compared])[0, 1]

plt.title(f"Сильная корреляция r = {strong_corr:.2f}");

def cut_str(s: str, max_length: int = 30) -> str:
    """Обрезает строку `s`, если число символов в ней превосходит `max_length`"""

    return s[:max_length] + "..." if len(s) > max_length else s


# Считаем корреляции между индикаторами и визуализируем что получилось
corr_matrix = (
    X_train.corr().abs()
)  # Знак коэффициентов корреляции нам будет, скорее, мешать

# Переименуем индексы массива в полные названия
corr_matrix.index = indicators[corr_matrix.index].apply(cut_str)

# Визуализация
plt.figure(figsize=(12, 7))
sns.heatmap(corr_matrix).set(title="Корреляции между признаками");

def calculate_single_vif(X_train: pd.DataFrame, feature: str):
    """Вычисляет коэффициент инфляции дисперсии (VIF) для отдельного признака.

    Args:
        X_train (pd.DataFrame): обучающая выборка
        feature (str): название признака

    Returns:
        float: значение VIF для указанного объекта.
    """
    # Один из признаков перемещаем в отдельную переменную
    X_VIF, y_VIF = X_train.drop(columns=feature), X_train[feature]

    # Делим выборки на трейн и тест
    X_train_VIF, X_test_VIF, y_train_VIF, y_test_VIF = train_test_split(
        X_VIF, y_VIF, test_size=0.5, random_state=random_state
    )

    # Создание модели, обучение и подсчёт метрики
    model_VIF = LinearRegression()
    model_VIF.fit(X_train_VIF, y_train_VIF)
    r2 = r2_score(y_test_VIF, model_VIF.predict(X_test_VIF))

    return np.inf if r2 == 1 else 1 / (1 - r2)

VIFs = []

for feature in X_train.columns:
    VIFs.append(calculate_single_vif(X_train, feature))

plt.figure(figsize=(10, 4))

sns.histplot(VIFs, bins=10, binrange=(0, 800))
plt.title("Гистограмма распределения VIF");

# Создаём сетку параметров регуляризации.
# Здесь удобнее работать с логарифмической шкалой
alphas = {"lasso": np.logspace(-5, 2, 500), "ridge": np.logspace(-8, 7, 500)}

# Создаём списки, куда будем записывать метрику качества
r2_train = {"lasso": [], "ridge": []}
r2_test = {"lasso": [], "ridge": []}

for model_name, model in {"lasso": Lasso(), "ridge": Ridge()}.items():

    for alpha in alphas[model_name]:

        # Устанавливаем очередное альфа и обучаем
        model.set_params(alpha=alpha)
        model.fit(X_train, y_train)

        # Считаем метрику и добавляем в списки
        r2_test[model_name].append(r2_score(y_test, model.predict(X_test)))
        r2_train[model_name].append(r2_score(y_train, model.predict(X_train)))

best_alpha_lasso = alphas["lasso"][np.argmax(r2_test["lasso"])]
best_alpha_ridge = alphas["ridge"][np.argmax(r2_test["ridge"])]
print(
    f'Для Lasso: R2_best = {np.max(r2_test["lasso"]):.2f} (test) при alpha = {best_alpha_lasso:.2f}\nДля Ridge: R2_best = {np.max(r2_test["ridge"]):.2f} (test) при alpha = {best_alpha_ridge:.2f}'
)

Для Lasso: R2_best = 0.54 (test) при alpha = 0.27
Для Ridge: R2_best = 0.59 (test) при alpha = 125.95

fig, axes = plt.subplots(1, 2, figsize=(15, 6), tight_layout=True)

for ax, model_name in zip(axes, alphas.keys()):

    ax.plot(alphas[model_name], r2_test[model_name], label="test", lw=3, ls="-")
    ax.plot(alphas[model_name], r2_train[model_name], label="train", lw=3, ls="--")
    ax.set_xscale("log")  # Логарифмический масштаб
    ax.set_xlabel("Значение $\\alpha$")
    ax.set_title(model_name)
    ax.legend()

axes[0].set_ylabel("$R^2$")
fig.suptitle("Зависимость качества модели от параметра регуляризации");

# Задаём сетки гиперпараметров
alphas = np.logspace(-2, 2, 50)
l1_ratios = np.linspace(0, 1.0, 50)

# Метрика качества - функция двух переменных alpha и l1_ratio,
# поэтому создаём пустой двумерный массив
r2_elastic = np.empty((len(l1_ratios), len(alphas)))

# Итерируемся по всем комбинациям alpha и l1_ratio
for alpha_idx, alpha in enumerate(alphas):
    for l1r_idx, l1_ratio in enumerate(l1_ratios):

        # Создаём модель ElasticNet
        model = ElasticNet(alpha=alpha, l1_ratio=l1_ratio, random_state=random_state)

        # Обучаем
        model.fit(X_train, y_train)

        # Делаем предсказание
        y_pred = model.predict(X_test)

        # Считаем метрику и сохраненяем в массив
        r2_elastic[l1r_idx, alpha_idx] = r2_score(y_test, y_pred)

arg = np.argmax(r2_elastic)

best_elastic = {
    "alpha": alphas[arg % alphas.size],
    "l1_ratio": l1_ratios[arg // l1_ratios.size],
}

print(f"R2_best = {np.max(r2_elastic):.2f} при {best_elastic}")

R2_best = 0.58 при {'alpha': 0.7543120063354615, 'l1_ratio': 0.0}

models = {
    "Linear": LinearRegression().fit(X_train, y_train),
    "Ridge": Ridge(alpha=best_alpha_ridge).fit(X_train, y_train),
    "Lasso": Lasso(alpha=best_alpha_lasso).fit(X_train, y_train),
}

# Коэффициенты моделей сохраним в виде DataFrame
coef_df = pd.DataFrame(
    {name: model.coef_ for name, model in models.items()},
    index=indicators[X_train.columns],
)

(coef_df["Lasso"] != 0).sum()

17

# Выбираем признаки с ненулевыми коэффициентами у Lasso
top_features = coef_df["Lasso"][(coef_df["Lasso"] != 0)].index

# Вытаскиваем аналогичные коэффициенты у МНК и Ridge
top_coef_df = coef_df.loc[top_features].reset_index(names="Индикатор")

# Сортируем по алфавиту
top_coef_df.sort_values(by="Индикатор", inplace=True)

# Преобразуем таблицу в "длинный формат" для удобной визуализации
top_coef_df = top_coef_df.melt(
    id_vars="Индикатор", var_name="Модель", value_name="Значение коэффициента"
)

# Визуализация
plt.figure(figsize=(15, 15))
sns.barplot(
    data=top_coef_df, y="Индикатор", x="Значение коэффициента", hue="Модель", orient="h"
)
plt.grid("both")
plt.xscale("symlog")  # Логарифмический масштаб с учётом знака
plt.xlim((-10, 10))
plt.vlines(0, *plt.gca().get_ylim(), color="black")
plt.title("Визуализация некоторых коэффициентов линейных моделей");

	Country Name	Country Code	Indicator Name	Indicator Code	1960	1961	1962	1963	1964	1965	...	2014	2015	2016	2017	2018	2019	2020	2021	2022	2023
0	Africa Eastern and Southern	AFE	Access to clean fuels and technologies for coo...	EG.CFT.ACCS.ZS	NaN	NaN	NaN	NaN	NaN	NaN	...	17.401410	17.911234	18.463874	18.924037	19.437054	20.026254	20.647969	21.165877	21.863139	NaN
1	Africa Eastern and Southern	AFE	Access to clean fuels and technologies for coo...	EG.CFT.ACCS.RU.ZS	NaN	NaN	NaN	NaN	NaN	NaN	...	6.728819	7.005877	7.308571	7.547226	7.875917	8.243018	8.545483	8.906711	9.261320	NaN
2	Africa Eastern and Southern	AFE	Access to clean fuels and technologies for coo...	EG.CFT.ACCS.UR.ZS	NaN	NaN	NaN	NaN	NaN	NaN	...	38.080931	38.422282	38.722108	38.993157	39.337872	39.695279	40.137847	40.522209	41.011132	NaN

	Indicator Name
Indicator Code
DT.NFL.UNTA.CD	Net official flows from UN agencies, UNTA (cur...
GC.TAX.TOTL.GD.ZS	Tax revenue (% of GDP)
SP.POP.65UP.TO.ZS	Population ages 65 and above (% of total popul...

	Country Code	Indicator Name	Indicator Code	1960	1961	1962	1963	1964	1965	1966	...	2014	2015	2016	2017	2018	2019	2020	2021	2022	2023
Country Name
Northern Mariana Islands	1496	1496	1496	83	86	86	86	86	89	86	...	191	201	189	189	190	190	180	130	104	75
Jamaica	1496	1496	1496	181	203	217	221	212	220	232	...	1052	1030	1076	1057	1076	991	863	806	680	355
Liberia	1496	1496	1496	125	151	156	160	155	161	161	...	923	780	903	886	753	749	763	673	573	197
Small states	1496	1496	1496	117	127	129	128	129	133	130	...	766	771	753	746	756	725	654	605	524	303
Chile	1496	1496	1496	261	285	307	314	309	318	317	...	1072	1150	1074	1133	996	987	970	862	766	423
Uzbekistan	1496	1496	1496	86	89	89	89	89	91	89	...	1027	1030	1008	1024	1034	1028	986	930	759	416
Luxembourg	1496	1496	1496	135	133	136	135	135	139	144	...	1034	1034	1009	1004	971	975	954	857	721	369

Indicator Code	AG.LND.CREL.HA	AG.PRD.CREL.MT	AG.PRD.CROP.XD	AG.PRD.FOOD.XD	AG.PRD.LVSK.XD	AG.YLD.CREL.KG	BG.GSR.NFSV.GD.ZS	BM.GSR.CMCP.ZS	BM.GSR.FCTY.CD	BM.GSR.GNFS.CD	...	VA.EST	VA.NO.SRC	VA.PER.RNK	VA.PER.RNK.LOWER	VA.PER.RNK.UPPER	VA.STD.ERR	VC.BTL.DETH	VC.IDP.NWCV	VC.IDP.NWDS	VC.IDP.TOCV
Country Code
NZL	113518.0	967398.93	108.20	100.57	99.00	8522.0	11.645655	45.100731	1.409778e+10	7.134960e+10	...	1.639671	10.0	99.516907	94.685989	100.000000	0.137781	NaN	NaN	2800.0	NaN
ITA	3010620.0	14300570.00	93.95	96.97	102.34	4750.0	12.904907	44.904146	8.588524e+10	7.855490e+11	...	1.070306	10.0	82.608696	75.362320	89.371979	0.108861	NaN	NaN	4100.0	NaN
HKG	0.0	0.06	101.49	141.66	149.97	2079.6	40.787834	45.344914	2.071293e+11	6.828809e+11	...	-0.406678	9.0	35.265701	29.951691	40.579712	0.138847	NaN	NaN	330.0	NaN
BIH	207036.0	1431169.00	128.06	113.39	83.81	6912.7	15.812834	19.064109	1.058158e+09	1.516245e+10	...	-0.325713	9.0	37.681160	32.850243	42.028984	0.127861	NaN	NaN	78.0	91000.0
ARG	18063791.0	91584182.57	107.33	109.39	113.91	5070.0	5.648610	39.767467	1.540821e+10	9.739913e+10	...	0.543527	12.0	62.801933	58.454105	70.531403	0.127701	NaN	NaN	730.0	NaN

Phystech@DataScience ¶

Семинар. Регуляризация на реальных данных¶

Введение¶

Смотрим на данные¶

Постановка задачи¶

1. Подготовка данных¶

1.1 Выделение таргета и признаков¶

1.2 Выбор признаков¶

1.3 Стандартизация¶

1.4 Заполняем пропуски¶

1.5 Отбираем признаки¶

2. Что не так с методом наименьших квадратов?¶

3. Как распознать мультиколлинеарность?¶

3.1 Число обусловленности¶

3.1 Матрица корреляций¶

3.3 Коэффициент инфляции дисперсии (VIF)¶

4. Как бороться с мультиколлинеарностью?¶

4.1 Регуляризация¶

Краткое напоминание¶

Применяем методы Ridge и Lasso¶

Применяем ElasticNet¶

Визуализируем коэффициенты¶

5. Интерпретация результатов¶

Заключение¶

	Country Name
Country Code
AFE	Africa Eastern and Southern
AFW	Africa Western and Central
ARB	Arab World

	Country Name
Country Code
TLA	Latin America & the Caribbean (IDA & IBRD coun...
LTE	Late-demographic dividend
MEA	Middle East & North Africa
TEC	Europe & Central Asia (IDA & IBRD countries)
SST	Small states

Phystech@DataScience¶

Семинар. Регуляризация на реальных данных¶

Введение¶

Смотрим на данные¶

Постановка задачи¶

1. Подготовка данных¶

1.1 Выделение таргета и признаков¶

1.2 Выбор признаков¶

1.3 Стандартизация¶

1.4 Заполняем пропуски¶

1.5 Отбираем признаки¶

2. Что не так с методом наименьших квадратов?¶

3. Как распознать мультиколлинеарность?¶

3.1 Число обусловленности¶

3.1 Матрица корреляций¶

3.3 Коэффициент инфляции дисперсии (VIF)¶

4. Как бороться с мультиколлинеарностью?¶

4.1 Регуляризация¶

Краткое напоминание¶

Применяем методы Ridge и Lasso¶

Применяем ElasticNet¶

Визуализируем коэффициенты¶

5. Интерпретация результатов¶

Заключение¶

Phystech@DataScience ¶