Введение в линейные модели со statsmodels

Важно: примеры используют pandas, NumPy и Seaborn для иллюстрации данных; код предназначен для выполнения в интерактивной среде (Jupyter, IPython) или в скриптах.

Что такое statsmodels?

statsmodels — это библиотека Python для запуска стандартных статистических тестов и моделирования. Она ориентирована на регрессионный анализ и инструменты, знакомые по эконометрике, но полезна всем, кому нужны надежные статистические оценки и диагностические тесты.

Пояснение в одну строку: statsmodels даёт детализированные статистические сводки, тесты и диагностические инструменты, которые часто отсутствуют в более «машинно-обученческих» библиотеках.

Почему стоит выбрать statsmodels:

• Подробные таблицы результатов (коэффициенты, стандартные ошибки, t- и p-значения, CI, AIC/BIC).
• Совместимость с формульным синтаксисом, похожим на R (patsy).
• Набор диагностик: тесты на гетероскедастичность, автокорреляцию, нормальность остатков.
• Лёгкая проверка согласованности с R, Stata и SAS в академических и прикладных работах.

Советы по использованию: начните с pandas DataFrame, проверяйте типы колонок (числовые vs категориальные) и чистоту данных (пропуски, выбросы). Для визуализации удобно использовать Seaborn или Matplotlib.

Установка и базовый импорт

Если у вас ещё не установлены пакеты:

pip install statsmodels seaborn pandas numpy

Импорты, которые мы будем использовать в примерах:

import numpy as np
import pandas as pd
import seaborn as sns
import statsmodels.api as sm
import statsmodels.formula.api as smf

sns.set_theme()

Пояснение: smf предоставляет формульный API (R-подобный), а sm — более низкоуровневый API, где вы сами собираете матрицу признаков (design matrix).

Простой пример: линейная регрессия на примере чаевых

Мы будем использовать датасет tips из Seaborn. Это реальный учебный пример: в нём есть счёт (total_bill), чаевые (tip), размер партии (size) и другие признаки.

tips = sns.load_dataset('tips')
tips.head()

Построим рассеяние и визуально оценим связь tip vs total_bill:

sns.relplot(x='total_bill', y='tip', data=tips)

Регрессия с визуализацией:

sns.regplot(x='total_bill', y='tip', data=tips)

Чтобы получить численные коэффициенты, применим statsmodels с формулой «R-стиля» и OLS (Ordinary Least Squares — метод наименьших квадратов):

results = smf.ols('tip ~ total_bill', data=tips).fit()
print(results.summary())

Или в интерактивной сессии просто:

results.summary()

Ключевые строки таблицы результатов:

• coef — оценённые коэффициенты модели (наклон и интерсепт).
• std err — стандартная ошибка оценок.
• t — t-статистика для проверки нулевой гипотезы coef = 0.
• P>|t| — p-значение; обычно порог 0.05.
• [0.025 0.975] — 95% доверительный интервал для коэффициента.

Пояснение: формула ‘tip ~ total_bill’ задаёт tip как зависимую переменную, а total_bill как регрессор. Тильда (~) — синтаксис указания зависимости.

Прямое использование API (матрицы проектирования)

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

x = np.linspace(-10, 10, 100)
y = 2 * x - 3
X = sm.add_constant(x)  # добавляет столбец единиц как интерсепт
model = sm.OLS(y, X)
results = model.fit()
print(results.summary())

Заметьте: в низкоуровневом API функция называется OLS (вверхними буквами). Вы контролируете порядок колонок в матрице X.

Множественная регрессия

Чтобы добавить в модель дополнительные регрессоры, просто перечислите их через плюс:

results = smf.ols('tip ~ total_bill + size', data=tips).fit()
results.summary()

Формула ‘y ~ x1 + x2’ добавляет оба регрессора. Результат — модель с несколькими коэффициентами. Для взаимодействий используйте ‘:’ или ‘*’ (см. секцию про взаимодействия дальше).

Нелинейные регрессоры

Чтобы встроить полиномиальные или трансформированные признаки, используйте I() вокруг выражения, например для квадратичной модели:

x = np.linspace(-10, 10, 100)
y = 3 * x**2 + 2 * x + 5
df = pd.DataFrame({'x': x, 'y': y})

results = smf.ols('y ~ x + I(x**2)', data=df).fit()
results.summary()

Пояснение: I(x2) сообщает, что это арифметическая операция над колонкой x, а не отдельная переменная ‘x2’.

Интерпретация результатов регрессии: на что смотреть

1. R-squared и Adj. R-squared. R^2 показывает долю дисперсии зависимой переменной, объяснённую моделью. Adjusted R^2 корректирует на число регрессоров — полезно в множественной регрессии.
2. Коэффициенты (coef). Для линейной модели y = b0 + b1x1 + b2x2 + … коэффициенты показывают ожидаемое изменение y при изменении x на единицу (при прочих равных).
3. Стандартные ошибки и p-значения. Низкое p-значение (обычно <0.05) указывает на статистическую значимость коэффициента.
4. AIC/BIC и F-статистика. Помогают сравнивать модели по качеству и сложности.
5. Остатки (residuals). Проверьте их распределение (QQ-plot), автокорреляцию и гетероскедастичность.

Пример интерпретации: если coef(total_bill) = 0.11 и p < 0.001, это означает, что при увеличении счёта на 1 единицу чаевые в среднем увеличиваются на 0.11, и эффект статистически значим.

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

Хорошая модель — это не только значимые коэффициенты, но и корректные допущения:

• Нормальность остатков: визуально (QQ-plot) или тестом Jarque–Bera (sm.stats.jarque_bera(resid)).
• Гомоскедастичность (однородность дисперсий): тест Бреуша–Пагана (sm.stats.diagnostic.het_breuschpagan).
• Автокорреляция остатков: тест Дарбина–Уотсона (sm.stats.stattools.durbin_watson).
• Мультиколлинеарность: VIF (variance inflation factor) из statsmodels.stats.outliers_influence.

Примеры кода:

# Остатки
resid = results.resid

# Jarque-Bera
import statsmodels.stats.api as sms
jb = sms.jarque_bera(resid)

# Breusch-Pagan
from statsmodels.stats.diagnostic import het_breuschpagan
bp_test = het_breuschpagan(resid, results.model.exog)

# Durbin-Watson
from statsmodels.stats.stattools import durbin_watson
dw = durbin_watson(resid)

# VIF
from statsmodels.stats.outliers_influence import variance_inflation_factor
X = results.model.exog  # матрица признаков, включая константу
vif = pd.DataFrame({'feature': results.model.exog_names,
                    'VIF': [variance_inflation_factor(X, i) for i in range(X.shape[1])]})

Интерпретация: VIF > 10 часто указывает на сильную мультиколлинеарность; Durbin-Watson около 2 — отсутствие автокорреляции; значимый BP-тест указывает на гетероскедастичность.

Коррекция стандартных ошибок (робастные ошибки)

Если присутствует гетероскедастичность, можно использовать робастные стандартные ошибки (heteroscedasticity-consistent):

robust_res = results.get_robustcov_results(cov_type='HC3')
print(robust_res.summary())

Пояснение: HC0/HC1/HC2/HC3 — разные варианты корректировок; HC3 более консервативен при наличии влиятельных наблюдений.

Работа с категориальными переменными и взаимодействиями

В формулах statsmodels вы можете явно указать категорию через C():

penguins = sns.load_dataset('penguins')

penguin_lm = smf.ols('bill_length_mm ~ C(species)', data=penguins).fit()
print(penguin_lm.summary())

Если переменная уже имеет тип ‘category’ в pandas, statsmodels корректно создаст фиктивные переменные (dummy). Для взаимодействий используйте * или :

• x1 * x2 = x1 + x2 + x1:x2
• x1:x2 = только произведение взаимодействия

Пример двухфакторного (двухпутевого) ANOVA:

penguin_multi_lm = smf.ols('bill_length_mm ~ species * island', data=penguins).fit()
results_anova = sm.stats.anova_lm(penguin_multi_lm)
print(results_anova)

Однофакторный ANOVA (one-way)

ANOVA отвечает на вопрос: отличается ли среднее значение числовой переменной между группами категориальной переменной?

penguins = sns.load_dataset('penguins')
penguin_lm = smf.ols('bill_length_mm ~ species', data=penguins).fit()
results = sm.stats.anova_lm(penguin_lm)
print(results)

Интерпретация: смотрим на p-value у фактора. Низкое p-value указывает на статистически значимое различие средних между группами.

Когда statsmodels не подходит (ограничения и альтернативы)

Counterexamples / когда это не лучший выбор:

• Для скоростного обучения и большой инфраструктуры машинного обучения (pipelines, кросс-валидация, регуляризация) предпочтителен scikit-learn (sklearn.linear_model).
• Для очень больших наборов данных, где требуется распределённое вычисление, лучше рассмотреть Spark ML или специализированные решения.
• Для нейронных сетей и нелинейных моделей с большим числом признаков — PyTorch/TF.

Альтернативные подходы:

• scikit-learn для простых и быстрых пайплайнов и кросс-валидации.
• statsmodels — для интерпретации и статистической отчётности.
• R (lm, glm) — если нужно воспроизводимость с уже существующим R-кодом или пакетом.

Практический мини-метод (SOP) для регрессии и ANOVA

1. Эксплоративный анализ (EDA): визуализация, проверка описательной статистики, поиск пропусков и выбросов.
2. Подготовка данных: типы колонок, кодирование категорий, масштабирование при необходимости.
3. Простая модель: начните с 1–2 регрессоров, проверьте качество и остатки.
4. Диагностика: проверка нормальности остатков, гетероскедастичности, автокорреляции, VIF.
5. Улучшение модели: добавьте трансформации (лог, квадратичные), взаимодействия, или используйте робастные ошибки.
6. Отчётность: представьте коэффициенты с доверительными интервалами, p-значениями и диагностикой.

Критерии приёмки модели:

• Остатки не демонстрируют сильных отклонений от нормальности без объективной причины.
• Гетероскедастичность отсутствует или скорректирована робастными ошибками.
• VIF в допустимых пределах или коллинеарные переменные устранены/объединены.
• Модель интерпретируема и соответствует задаче исследования.

Чеклисты по ролям

Data Scientist:

• Проверить мультиколлинеарность (VIF).
• Выполнить кросс-валидацию при сравнении моделей (если применимо).
• Использовать регуляризации (Ridge/Lasso) при большом числе признаков.

Аналитик / Исследователь:

• Подготовить таблицу коэффициентов и доверительных интервалов.
• Включить визуализацию остатков и QQ-plot.
• Проверить чувствительность результатов к включению/исключению переменных.

Преподаватель / репортёр:

• Пояснить значение p-value, R^2, и ограничений модели.
• Указать, какие допущения проверялись и какие коррекции применялись.

Советы по миграции и совместимости (R → statsmodels)

Ментальная модель: формула в statsmodels очень похожа на R, но есть различия в кодировке категорий и контрастных кодах. Если у вас есть R-скрипт:

• Проверьте, как в R кодируются факторы; в statsmodels можно явно указать тип pd.Categorical или использовать C() в формуле.
• Переобучение моделей может дать слегка отличающиеся числовые результаты из-за разных настроек по умолчанию — сравнивайте таблицы и тесты.

Рекомендации по хорошей практике

• Всегда визуализируйте данные перед моделью.
• Начинайте с простой модели и постепенно усложняйте.
• Документируйте предположения и тесты диагностики.
• Используйте робастные стандартные ошибки при наличии гетероскедастичности.

Примеры типичных ошибок и как их исправлять

1. Пропуски в данных приводят к уменьшению объёма выборки — используйте pd.dropna() аккуратно или методы имфутации.
2. Неправильная кодировка категорий — проверьте dtype и используйте C(var) для явного указания категориальной переменной.
3. Автокорреляция в временных рядах — вместо OLS используйте модели для временных рядов (ARIMA, GLSAR) или включайте лагающие переменные.

Факто-бокс: что помнить

• OLS = Ordinary Least Squares — базовый метод для линейной регрессии.
• p-value < 0.05 обычно используется как порог значимости.
• VIF > 10 — признак проблем мультиколлинеарности.
• Durbin–Watson ≈ 2 — отсутствие автокорреляции.

Визуальная помошь: решение «Какой метод использовать?»

flowchart TD
  A[Начать с данных] --> B{Цель исследования}
  B -->|Интерпретация коэффициентов| C[statsmodels 'OLS/GLM']
  B -->|Скорость и ML пайплайны| D[scikit-learn]
  B -->|Большие данные/распределение| E[Spark ML]
  C --> F{Остатки нормальны?}
  F -->|Да| G[Документировать модель]
  F -->|Нет| H[Диагностика и коррекции]
  H --> I[Трансформации, робастные ошибки, GLS]

Краткая галерея краевых случаев

• Очень маленькая выборка (<30): будьте осторожны с интерпретацией p-значений.
• Сильная мультиколлинеарность: задумайтесь об объединении переменных или регуляризации.
• Выбросы с большой взвешенностью: проверьте influential points (sm.graphics.influence_plot).

Глоссарий в одну строку

• OLS: метод наименьших квадратов;
• R^2: доля объяснённой дисперсии;
• VIF: фактор инфляции дисперсии (коллинеарность);
• ANOVA: анализ дисперсии для сравнения групп.

Итог и рекомендации

statsmodels — отличный инструмент, когда вам нужны интерпретируемые модели и детальная статистическая отчётность. Для продакшена и быстрых ML-пайплайнов чаще подходят scikit-learn, но при исследовательской работе и академической отчётности statsmodels остаётся предпочтительным выбором. Используйте формулы, тщательно диагностируйте модели и документируйте все шаги.

Краткие шаги для практики:

1. Подготовьте данные и проведите EDA.
2. Сформулируйте простую модель.
3. Проведите диагностику остатков и проверку допущений.
4. При необходимости примените робастные ошибки или трансформации.
5. Задокументируйте выводы и представьте графики и таблицы результатов.

Статья дала практическое введение и рабочий набор приёмов для начала работы с линейными моделями в statsmodels. Начните с простого, затем проверяйте допущения и расширяйте модель только при необходимости.