Как обнаружить и устранить мультиколлинеарность с помощью Statsmodels в Питоне

Привет, Хабр!

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

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

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

В статье рассмотрим как обнаружить и устранить мультиколлинеарность с помощью Statsmodels в Питоне.

Коэффициент инфляции дисперсии

Коэффициент инфляции дисперсии или VIF — это мера, которая помогает обнаружить степень мультиколлинеарности в множественной регрессионной модели. Он показывает, насколько увеличивается дисперсия коэффициента оценки из-за линейной зависимости с другими предикторами.

Коэффициент VIF для переменной Xj рассчитывается как:

где Rj2​ — это коэффициент детерминации регрессии переменной Xj​ на все остальные предикторы в модели. VIF начинается от 1 и может увеличиваться до бесконечности. Значение VIF, равное 1, указывает на отсутствие корреляции между данной переменной и другими предикторами, а значительно большие значения свидетельствуют о серьезной мультиколлинеарности.

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

import pandas as pd
import numpy as np

data = pd.DataFrame({
    'advertising': np.random.rand(100) * 100,
    'price_level': np.random.rand(100) * 10,
    'store_design': np.random.rand(100) * 5,
    'sales': np.random.rand(100) * 500
})

Для расчета VIF каждого предиктора в модели юзаем функцию variance_inflation_factor из модуля statsmodels.stats.outliers_influence.

from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.tools.tools import add_constant

# добавление константы для перехвата
X = add_constant(data.drop('sales', axis=1))

# расчет VIF для каждого предиктора
VIFs = pd.DataFrame()
VIFs['Variable'] = X.columns
VIFs['VIF'] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
print(VIFs)

Результатом будет DataFrame с VIF для каждой переменной.

       Variable       VIF
0         const  7.894204
1   advertising  1.035176
2   price_level  1.024471
3  store_design  1.035239

Высокие значения VIF указывают на потенциальные проблемы с мультиколлинеарностью. Как правило, значение VIF выше 5 требует внимания, а выше 10 — серьезного рассмотрения изменений в модели.

Методы решения мультиколлинеарности

Один из самых простых способов борьбы с мультиколлинеарностью — удаление переменных, которые сильно коррелируют с другими предикторами. Для этого потребуется провести анализ корреляционной матрицы или коэффициентов VIF.

import pandas as pd
from statsmodels.stats.outliers_influence import variance_inflation_factor

def calculate_vif(data):
    vif_data = pd.DataFrame()
    vif_data["feature"] = data.columns
    vif_data["VIF"] = [variance_inflation_factor(data.values, i) for i in range(data.shape[1])]
    return vif_data

# предположим, data — это DataFrame с предикторами
vif_info = calculate_vif(data)
high_vif = vif_info[vif_info['VIF'] > 10]  # Установите порог, например 10
data = data.drop(columns=high_vif['feature'])

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

from sklearn.cross_decomposition import PLSRegression
from sklearn.datasets import make_regression

X, y = make_regression(n_samples=1000, n_features=10, noise=0.1)
pls = PLSRegression(n_components=3)
pls.fit(X, y)

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

from sklearn.linear_model import Lasso, Ridge

# LASSO регрессия
lasso = Lasso(alpha=0.1)
lasso.fit(X, y)
print("LASSO Coefficients:", lasso.coef_)

# Ridge регрессия
ridge = Ridge(alpha=1.0)
ridge.fit(X, y)
print("Ridge Coefficients:", ridge.coef_)

alpha контролирует силу регуляризации. Подбор оптимального значения alpha часто выполняется с помощью кросс-валидации.

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

import statsmodels.api as sm

def forward_selection(data, response):
    remaining = set(data.columns)
    remaining.remove(response)
    selected = []
    current_score, best_new_score = 0.0, 0.0
    while remaining and current_score == best_new_score:
        scores_with_candidates = []
        for candidate in remaining:
            formula = "{} ~ {}".format(response, ' + '.join(selected + [candidate]))
            score = sm.OLS.from_formula(formula, data).fit().aic
            scores_with_candidates.append((score, candidate))
        scores_with_candidates.sort()
        best_new_score, best_candidate = scores_with_candidates.pop(0)
        if current_score > best_new_score:
            remaining.remove(best_candidate)
            selected.append(best_candidate)
            current_score = best_new_score
    formula = "{} ~ {}".format(response, ' + '.join(selected))
    model = sm.OLS.from_formula(formula, data).fit()
    return model

model = forward_selection(data, 'dependent_variable')

---

Важно всегда учитывать мультиколлинеарность при интерпретации результатов регрессионных анализов.

А про регрессионный анализ и не только вы можете узнать в рамках практических онлайн-курсов по машинному обучению и Data Science. У коллег из OTUS есть большая линейка курсов по данным направлениям. Переходите в каталог и убедитесь сами.