Po co sprawdzać założenia regresji, zanim przejdziesz do wniosków?
„Model działa” kontra „model daje wiarygodne wnioski”
Regresja liniowa w Pythonie potrafi dopasować się do danych prawie zawsze. Wywołujesz metodę .fit(), dostajesz ładną tabelkę ze współczynnikami, wartościami p, R² i teoretycznie wszystko wygląda poprawnie. Tylko czy faktycznie możesz na tej podstawie podejmować decyzje biznesowe lub naukowe?
Kluczowa różnica leży między stwierdzeniem „model dobrze przewiduje”, a „wnioski statystyczne z modelu są poprawne”. Możesz mieć model z wysokim R², który świetnie „uczy się” struktury danych, ale łamie podstawowe założenia regresji: ma heteroskedastyczne błędy, kilka obserwacji całkowicie dominuje wynik, a reszty pokazują silną nieliniowość. Wtedy:
- wartości p są przekłamane,
- przedziały ufności są zbyt wąskie lub zbyt szerokie,
- testy istotności mogą wskazywać „brak efektu” tam, gdzie on faktycznie istnieje (lub odwrotnie).
Jeśli korzystasz z regresji po to, aby odpowiedzieć na pytanie „czy budżet reklamowy wpływa na sprzedaż?” lub „o ile średnio rośnie sprzedaż po zwiększeniu ceny o 1 jednostkę?”, potrzebujesz poprawnej inferencji, a nie tylko dobrego dopasowania.
Skutki łamania założeń w praktyce
Co się dzieje, gdy ignorujesz diagnostykę regresji w Pythonie i opierasz się tylko na gotowym outputcie z statsmodels lub scikit-learn? Typowe konsekwencje:
- Heteroskedastyczność – estymatory współczynników nadal są (w wielu przypadkach) nieobciążone, ale wariancje są błędnie oszacowane. W rezultacie błędy standardowe, wartości p i przedziały ufności są niewiarygodne.
- Nieliniowość – model „na siłę” próbuje dopasować prostą linię tam, gdzie relacja jest krzywoliniowa. Interpretacja pojedynczego współczynnika przestaje mieć sens, bo efekt zmiennej zależy od jej poziomu.
- Obserwacje wpływowe – kilka punktów o dużym lewarze (leverage) może całkowicie przekręcić regresję, zmieniając znak i wielkość współczynników. Biznesowo: jeden klient „olbrzym” zaburza wnioski o całej populacji.
- Zależność obserwacji – w danych szeregów czasowych czy panelowych klasyczne OLS nie uwzględnia korelacji w resztach, co znów psuje inferencję.
Jeżeli Twoim celem jest tylko „wyplucie” prognozy i oceniasz model wyłącznie po RMSE/MAE na zbiorze walidacyjnym, wiele z tych problemów może być mniej dotkliwych. Jeżeli jednak opierasz na współczynnikach decyzje typu: „tnijmy budżet TV, bo nie jest istotny statystycznie”, diagnostyka przestaje być opcjonalna.
Predykcja czy interpretacja parametrów – jakie masz oczekiwanie?
Kiedy siadasz do regresji, zadaj sobie krótkie pytanie: jaki masz główny cel?
- Jeśli interesuje Cię przede wszystkim dokładna prognoza (np. przewidzenie sprzedaży w następnym miesiącu), kluczowe będą: jakość dopasowania, walidacja krzyżowa, stabilność prognoz. Założenia klasycznej regresji są nadal ważne, ale część naruszeń da się obejść, korzystając z innych metod (np. drzew, XGBoost) lub walidacji.
- Jeśli Twoim celem jest interpretacja parametrów („o ile rośnie sprzedaż przy zwiększeniu budżetu digital o 1 jednostkę?”), wtedy diagnostyka regresji w Pythonie staje się krytyczna: potrzebujesz poprawnych błędów standardowych, rozsądnego rozkładu reszt i braku dominacji pojedynczych punktów.
Warto tę decyzję nazwać wprost jeszcze przed dopasowaniem modelu. Dzięki temu łatwiej będzie świadomie dobrać narzędzia: testy heteroskedastyczności, robust standard errors, transformacje, czy może zupełnie inny typ modelu.
Gdzie najbardziej boli ignorowanie założeń?
Najczęściej problemy wychodzą w obszarach, gdzie na podstawie regresji podejmujesz decyzje finansowe:
- Budżety marketingowe – jeżeli model budżet → sprzedaż ma heteroskedastyczne reszty, a do tego kilka kampanii o gigantycznym zasięgu dominuje wyniki, to wnioski typu „kanał A jest nieistotny, przesuńmy budżet do kanału B” mogą być zupełnie błędne.
- Modele cenowe – jeśli przy zmianach cen reszty rosną nieliniowo, prosta regresja może sugerować stałą elastyczność cenową, która w praktyce zmienia się w zależności od poziomu ceny i segmentu klienta.
Zanim więc wyślesz „ostateczny” raport, zatrzymaj się na moment i zadaj sobie pytanie: czy sprawdziłeś reszty, heteroskedastyczność i wpływ kluczowych obserwacji? Jeśli nie – kolejne sekcje pokazują, jak to zrobić w Pythonie krok po kroku.
Środowisko pracy w Pythonie i przykładowy zestaw danych
Statsmodels czy scikit-learn – co wybrać do diagnostyki regresji?
Jeśli Twoim celem jest diagnostyka regresji liniowej, w praktyce lepiej sprawdza się statsmodels niż scikit-learn. Dlaczego?
- statsmodels oferuje:
- pełną tabelkę wyników (współczynniki, błędy standardowe, wartości t, p-value),
- łatwy dostęp do reszt, wartości dopasowanych, miar wpływu (Cook’s distance, lewar),
- wbudowane testy: Breuscha–Pagana, White’a, Durbin–Watsona, Jarque–Bera i wiele innych.
- scikit-learn skupia się na uczeniu maszynowym:
- prostsza, spójna API do wielu algorytmów,
- dużo narzędzi do walidacji i grid search,
- ale brak pełnej diagnostyki statystycznej out-of-the-box.
Rozsądne podejście w wielu projektach jest takie: modelowanie i diagnostyka w statsmodels, finalny pipeline predykcyjny w scikit-learn (jeśli potrzebujesz później klasycznego workflow ML). Reszty, heteroskedastyczność i wpływ obserwacji zdecydowanie wygodniej analizuje się w statsmodels.
Minimalny setup: biblioteki i konfiguracja
Do pracy z regresją i jej diagnostyką przydadzą się:
- pandas – wczytanie i przygotowanie danych,
- numpy – operacje numeryczne,
- statsmodels – estymacja i testy,
- matplotlib i seaborn – wizualizacje diagnostyczne.
Przykładowy setup:
import pandas as pd
import numpy as np
import statsmodels.api as sm
import statsmodels.stats.api as sms
from statsmodels.graphics.gofplots import qqplot
from statsmodels.graphics.regressionplots import plot_partregress_grid
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(style="whitegrid", context="notebook")
Jeżeli używasz Jupyter Notebooka, pamiętaj o:
%matplotlib inline
Przykładowy model: sprzedaż ~ budżet reklamowy + sezonowość
Załóżmy prostą sytuację: analizujesz sprzedaż miesięczną, a jako predyktory masz:
- budżet reklamowy (
budget), - sezonowość (np. miesiąc, zakodowany jako zmienne zero-jedynkowe albo indeks
month), - opcjonalnie: liczba aktywnych kampanii (
campaigns).
Dane mogą wyglądać tak:
df = pd.read_csv("sales_marketing.csv")
# zmienna zależna
y = df["sales"]
# zmienne niezależne
X = df[["budget", "campaigns"]]
# dodanie stałej do modelu (intercept)
X = sm.add_constant(X)
model = sm.OLS(y, X)
results = model.fit()
print(results.summary())
Na tym etapie masz już dopasowany model OLS i tabelę z parametrami. Pytanie brzmi: czy można jej ufać? Żeby to ocenić, trzeba przejść przez analizę reszt, heteroskedastyczności i wpływu obserwacji.
Jakie założenia regresji liniowej faktycznie są istotne w praktyce?
Kiedy linearyzacja relacji jest wystarczająca?
Klasyczna regresja liniowa zakłada, że relacja między zmienną zależną a predyktorami jest liniowa (w parametrach). To nie oznacza, że wykres musi być „sztywną” prostą – do modelu można dołożyć:
- składniki kwadratowe (
x**2), - logarytmy (
np.log(x)), - interakcje (
x1 * x2).
Pytanie brzmi: czy istnieje taka transformacja zmiennych, że model liniowy jest rozsądny? Jeśli na wykresie wartości dopasowanych vs reszty widzisz wyraźny kształt U lub S, oznacza to, że sama prosta jest niewystarczająca. W praktyce, zanim zaczniesz bawić się zaawansowanymi algorytmami, spróbuj:
- dodać składnik kwadratowy,
- przelogarytmować zmienną zależną lub predyktory,
- pomodelować interakcje między głównymi predyktorami.
Jaką masz tu motywację: zależy Ci na prostocie interpretacji, czy na możliwie najlepszym dopasowaniu? To wpływa na decyzję, jak bardzo komplikować model.
Niezależność obserwacji – gdzie się sypie najczęściej?
Założenie niezależności reszt bywa łamane przede wszystkim w:
- danych szeregów czasowych – sprzedaż w lutym jest silnie związana ze sprzedażą w styczniu, nawet po uwzględnieniu znanych predyktorów,
- danych panelowych – wiele obserwacji dla tych samych jednostek (klientów, sklepów, regionów) powoduje korelację wewnątrz grup.
Jeśli w takich sytuacjach zastosujesz prosty OLS, to:
- współczynniki mogą być nadal nieobciążone,
- ale błędy standardowe będą niedoszacowane,
- wartości p będą zbyt „optymistyczne” – więcej zmiennych wyjdzie istotnych, niż faktycznie jest.
Wyjścia są różne: modele z błędami skorelowanymi (np. Newey–West robust errors), modele mieszane, regresja panelowa, ARIMA/VAR dla czystych szeregów. Klucz jest jeden: jeśli dane są czasowe lub panelowe, klasyczna diagnostyka regresji w Pythonie to za mało – trzeba dodatkowo zbadać autokorelację reszt.
Homoskedastyczność i normalność reszt – co jest ważniejsze?
Założenia:
- Homoskedastyczność – wariancja błędów jest stała w całym zakresie predyktorów.
- Normalność reszt – reszty mają rozkład normalny (lub zbliżony).
Które z nich jest praktycznie ważniejsze? W zdecydowanej większości analiz kluczowa jest homoskedastyczność, ponieważ:
- wpływa bezpośrednio na poprawność błędów standardowych i testów istotności,
- jest testowalna i mamy dobre poprawki (robust standard errors, WLS).
Normalność reszt jest istotna głównie przy małych próbach i wtedy, gdy zależy Ci na dokładnym spełnieniu założeń testów opartych na rozkładzie t. Przy dużych próbach (kilkadziesiąt, kilkaset obserwacji i więcej) centralne twierdzenie graniczne „ratuje” wiele sytuacji: nawet jeśli rozkład reszt jest lekko skośny, inferencja często działa znośnie, o ile nie ma skrajnych outlierów.
Brak wielokolinearności i brak punktów dominujących
Wielokolinearność i obserwacje wpływowe często niszczą interpretowalność modelu:
- wielokolinearność – jeśli dwie zmienne są mocno skorelowane, trudno rozdzielić ich indywidualny wpływ; błędy standardowe rosną, a wartości p robią się dziwnie wysokie mimo dobrego dopasowania,
- obserwacje wpływowe – punkty o dużym lewarze mogą zmieniać znak i wielkość współczynników; to tak jakby jedna firma w zbiorze B2B decydowała za cały rynek.
Dobra diagnostyka regresji w Pythonie powinna więc uwzględniać:
- sprawdzenie współczynników VIF (Variance Inflation Factor) dla kluczowych predyktorów,
- analizę lewaru, Cook’s distance, DFFITS z
results.get_influence()w statsmodels.
Jeżeli budujesz model głównie po to, by uzasadnić decyzję typu „odcinamy jedną zmienną z modelu”, to ignorowanie tych kwestii może prowadzić do bardzo kosztownych pomyłek.

Reszty w praktyce: jak je obliczyć i zrozumieć w Pythonie?
Rodzaje reszt: zwykłe, standaryzowane, studentyzowane
Najpierw pytanie: jakie reszty chcesz oglądać? W statsmodels masz kilka opcji:
Jak dobrać typ reszt do konkretnego zadania?
Zanim zaczniesz rysować wykresy, zadaj sobie jedno pytanie: czy zależy Ci bardziej na ocenie ogólnego dopasowania, czy na wyłapaniu pojedynczych, podejrzanych obserwacji? Od odpowiedzi zalega wybór rodzaju reszt:
- reszty zwykłe – dobre do szybkiego podglądu kształtu (np. nieliniowość, heteroskedastyczność),
- reszty standaryzowane – lepsze, gdy chcesz mieć porównywalną skalę błędów dla wszystkich obserwacji,
- reszty studentyzowane – kluczowe przy szukaniu outlierów i punktów wpływowych.
Jaki masz cel na teraz: „czy model w ogóle ma sens?” czy raczej „które punkty mi go rozwalają?”.
Jak wyciągnąć reszty z obiektu results w statsmodels?
Statsmodels trzyma interesujące Cię wielkości bezpośrednio w obiekcie results. Najważniejsze atrybuty to:
results.resid– reszty zwykłe,results.fittedvalues– wartości dopasowane,results.get_influence()– pełny pakiet informacji o wpływie i resztach znormalizowanych.
# podstawowe reszty i wartości dopasowane
residuals = results.resid
fitted = results.fittedvalues
# wpływ obserwacji
influence = results.get_influence()
# reszty standaryzowane (standardized)
standardized_resid = influence.resid_studentized_internal
# reszty studentyzowane (externally studentized)
studentized_resid = influence.resid_studentized_external
Jeśli zależy Ci na identyfikacji outlierów, używaj przede wszystkim studentized_resid. Próg? Zazwyczaj:
|reszta studentyzowana| > 2– obserwacja potencjalnie problematyczna,|reszta studentyzowana| > 3– poważny kandydat na outlier.
Wykres reszt vs wartości dopasowane – pierwszy test „na oko”
Jeśli masz dopasowany model, kolejny ruch jest prosty: wykres reszt vs fitted. Po co? Żeby sprawdzić:
- czy widać systematyczny kształt (U, S, wachlarz),
- czy wariancja reszt rośnie lub maleje wraz z poziomem predykcji,
- czy kilka punktów odstaje wyraźnie od reszty chmury.
fig, ax = plt.subplots(figsize=(8, 5))
sns.scatterplot(x=fitted, y=standardized_resid, ax=ax, alpha=0.7)
ax.axhline(0, color="red", linestyle="--", linewidth=1)
ax.set_xlabel("Wartości dopasowane")
ax.set_ylabel("Reszty standaryzowane")
ax.set_title("Reszty vs wartości dopasowane")
plt.show()
Na co patrzeć?
- brak kształtu – chmura przypadkowa, bez trendu – to dobry znak,
- kształt U/S – sugeruje nieliniowość, może brak składników kwadratowych/logarytmów,
- wachlarz (rosnący rozrzut) – klasyczny sygnał heteroskedastyczności.
Jeśli na tym wykresie widzisz kilka punków „wiszących” z boku, zadaj sobie pytanie: czy te obserwacje są reprezentatywne dla Twojej populacji, czy raczej są to nietypowe zdarzenia?
QQ-plot reszt – czy normalność to faktyczny problem?
Gdy już sprawdzisz kształt reszt, możesz zajrzeć na ich rozkład. Statsmodels oferuje funkcję qqplot:
fig = qqplot(standardized_resid, line="45", fit=True)
plt.title("QQ-plot reszt standaryzowanych")
plt.show()
Interpretacja jest prosta:
- punkty blisko linii 45° – rozkład reszt jest zbliżony do normalnego,
- mocno odgięte ogony – możliwe problemy z outlierami lub skośnością,
- systematyczne wygięcie – np. reszty bardziej „szpiczaste” lub „płaskie” niż normalne.
Pytanie kontrolne: jak duża jest Twoja próba? Jeżeli masz setki obserwacji, lekkie odchylenia od normalności zwykle nie przekreślają modelu, pod warunkiem, że:
- nie widzisz ewidentnych outlierów w ogonach,
- homoskedastyczność jest rozsądnie spełniona lub skorygowana.
Testy formalne: Jarque–Bera a praktyka decyzji
Gdy chcesz potwierdzić wrażenie z QQ-plota, sięgasz po testy formalne. W statsmodels wynik testu Jarque–Bera znajdziesz już w results.summary(), ale możesz też wywołać go ręcznie:
from statsmodels.stats.stattools import jarque_bera
jb_stat, jb_pvalue, skew, kurtosis = jarque_bera(residuals)
print("JB stat:", jb_stat)
print("p-value:", jb_pvalue)
print("Skośność:", skew)
print("Kurtoza:", kurtosis)
Co dalej, gdy p-value jest „złe” (np. < 0.05)? Najpierw odpowiedz sobie: czy Twoje wnioski opierają się głównie na testach istotności? Jeśli tak, rozważ:
- robust standard errors (np.
HC3), - transformację zmiennej zależnej (np. log-sprzedaży),
- usunięcie lub osobne modelowanie skrajnych outlierów.
Jeżeli model służy stricte do prognoz krótkoterminowych i weryfikujesz go głównie na zbiorach walidacyjnych, brak idealnej normalności reszt rzadko jest kluczowym problemem.
Heteroskedastyczność: jak ją wykryć i co z nią zrobić?
Wizualna diagnoza: wachlarz, „lejek” i zmienna wariancja
Heteroskedastyczność często widać jeszcze zanim odpalisz jakikolwiek test statystyczny. Jakie wzorce są najbardziej typowe?
- rozrzut reszt rośnie wraz z poziomem predykcji (efekt „lejka”),
- reszty są bardzo małe dla większości obserwacji, ale ogromne dla kilku segmentów (np. dużych klientów),
- w rysunku reszt względem konkretnego predyktora wyróżniają się „pasma” o różnej szerokości.
Możesz to sprawdzić nie tylko względem wartości dopasowanych, ale też osobno dla kluczowych zmiennych:
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
sns.scatterplot(x=fitted, y=standardized_resid, ax=axes[0], alpha=0.7)
axes[0].axhline(0, color="red", linestyle="--", linewidth=1)
axes[0].set_xlabel("Fitted")
axes[0].set_ylabel("Std. residuals")
axes[0].set_title("Reszty vs fitted")
sns.scatterplot(x=df["budget"], y=standardized_resid, ax=axes[1], alpha=0.7)
axes[1].axhline(0, color="red", linestyle="--", linewidth=1)
axes[1].set_xlabel("Budget")
axes[1].set_ylabel("Std. residuals")
axes[1].set_title("Reszty vs budget")
plt.tight_layout()
plt.show()
Zadaj sobie wtedy pytanie: czy wariancja błędu rośnie razem ze skalą zjawiska? W case’ach biznesowych (sprzedaż, przychody) to wręcz norma.
Test Breuscha–Pagana i White’a w statsmodels
Po wstępnej diagnozie wizualnej przychodzi czas na test formalny. W statsmodels masz wygodną funkcję het_breuschpagan i het_white:
from statsmodels.stats.diagnostic import het_breuschpagan, het_white
# zmienna zależna: reszty, regresory: macierz projektująca (X)
bp_test = het_breuschpagan(residuals, X)
white_test = het_white(residuals, X)
bp_stat, bp_pvalue, _, _ = bp_test
white_stat, white_pvalue, _, _ = white_test
print("Breusch-Pagan p-value:", bp_pvalue)
print("White p-value:", white_pvalue)
Interpretacja:
- p-value < 0.05 – odrzucasz hipotezę homoskedastyczności,
- p-value ≥ 0.05 – brak podstaw do odrzucenia, heteroskedastyczność nie jest statystycznie wykazana.
Jeśli testy wychodzą „źle”, zadaj sobie od razu drugie pytanie: czy zależy Ci bardziej na poprawnej inferencji (p-value, przedziały), czy na samej jakości predykcji? Od tego zależy strategia naprawcza.
Robust standard errors: HC0, HC1, HC2, HC3
Najprostsza odpowiedź na heteroskedastyczność to zmiana estymatora wariancji, bez ruszania samych współczynników. W statsmodels możesz to zrobić wprost:
# ponowne dopasowanie modelu z robust standard errors
results_hc3 = model.fit(cov_type="HC3")
print(results_hc3.summary())
Dlaczego HC3? W praktyce:
HC0/HC1– bardziej „klasyczne”,HC2/HC3– lepiej radzą sobie przy mniejszych próbach i punktach wpływowych.
Jeżeli Twoim głównym celem jest odpowiedź na pytanie „czy budżet jest istotny po kontrolowaniu kampanii?”, zastosowanie cov_type="HC3" jest często najrozsądniejszym pierwszym krokiem, zanim zaczniesz przeprojektowywać cały model.
WLS i transformacje – kiedy zmieniać model zamiast samych błędów?
Czasem jednak nie chcesz tylko „naprawiać” błędów standardowych, ale też poprawić samo dopasowanie i strukturę błędu. Tu pojawiają się dwie główne ścieżki:
- Weighted Least Squares (WLS) – jeśli masz sensowną hipotezę, jak zmienia się wariancja błędu,
- transformacje – np. logarytm zmiennej zależnej, gdy błąd rośnie proporcjonalnie do poziomu zjawiska.
# przykład: wariancja błędu rośnie ~ kwadrat fitted
weights = 1 / (fitted**2 + 1e-6)
wls_model = sm.WLS(y, X, weights=weights)
wls_results = wls_model.fit()
print(wls_results.summary())
W praktyce biznesowej bardzo często lepiej zachowuje się model:
df["log_sales"] = np.log(df["sales"])
y_log = df["log_sales"]
X = sm.add_constant(df[["budget", "campaigns"]])
log_model = sm.OLS(y_log, X)
log_results = log_model.fit()
Jeśli po takiej transformacji wykres reszt staje się znacznie spokojniejszy, a testy heteroskedastyczności przestają „krzyczeć”, masz mocny argument, że to skala zmiennej była głównym źródłem problemu.
Wpływ obserwacji: lewar, Cook’s distance i DFFITS
Dlaczego kilka punktów może sterować całym modelem?
Zadaj sobie pytanie: czy w Twoich danych są obserwacje, które reprezentują zupełnie inny świat niż reszta? Duży klient, nietypowy miesiąc kryzysowy, jednostkowa awaria systemu – to wszystko przykłady obserwacji, które:
- mogą mieć duży lewar (są „daleko” w przestrzeni predyktorów),
- mogą mieć dużą resztę (model źle je opisuje),
- mogą łącznie mieć ogromny wpływ na współczynniki.
W statsmodels ten pakiet informacji dostarcza get_influence().
Lewar (hat values) – które obserwacje mają „głos decydujący”?
Lewar (hat value) opisuje, na ile dana obserwacja „ciągnie” dopasowanie w swoją stronę. Wyciągniesz go wprost:
influence = results.get_influence()
# lewar (hat values)
leverage = influence.hat_matrix_diag
Jak go czytać? Często używa się progu:
h_i > 2 * (k + 1) / n– gdziekto liczba predyktorów, anliczba obserwacji.
Możesz to policzyć i zidentyfikować podejrzane punkty:
n = len(df)
k = X.shape[1] - 1 # bez stałej
leverage_threshold = 2 * (k + 1) / n
high_leverage_points = np.where(leverage > leverage_threshold)[0]
print("Indeksy punktów o wysokim lewarze:", high_leverage_points)
Pytanie diagnostyczne: czy te punkty pochodzą z grupy, którą naprawdę chcesz modelować, czy są to np. jednorazowe duże zamówienia spoza „standardowego” rynku?
Cook’s distance – łączna miara wpływu obserwacji
Sama informacja o lewarze nie wystarcza. Potrzebujesz też reszty, żeby zrozumieć, jak bardzo punkt „psuje” model. Do tego służy Cook’s distance:
cooks_d, pvals = influence.cooks_distance
# sortowanie po największym wpływie
influential_points = np.argsort(cooks_d)[::-1]
print("Top 5 wpływowych punktów:")
for idx in influential_points[:5]:
print(idx, cooks_d[idx])
Często używa się orientacyjnego progu:
D_i > 4 / n– obserwacja potencjalnie wpływowa.
DFFITS, DFBETAS i reszty studentyzowane – jak mierzyć wpływ na prognozy i współczynniki?
Cook’s distance daje ogląd ogólny, ale czasem potrzebujesz precyzyjniejszej odpowiedzi: które obserwacje najmocniej zmieniają prognozy, a które same współczynniki? Tu pojawiają się trzy narzędzia:
- DFFITS – wpływ punktu na jego przewidywaną wartość,
- DFBETAS – wpływ punktu na każdy współczynnik z osobna,
- reszty studentyzowane – reszty skalowane własnym błędem standardowym.
W statsmodels większość tego masz „pod ręką” w obiekcie influence:
# DFFITS
dffits, dffits_threshold = influence.dffits
# DFBETAS
dfbetas = influence.dfbetas
# reszty studentyzowane (zewnętrzne)
studentized_resid = influence.resid_studentized_external
Jak to czytać? Zadaj sobie pytanie: czy patrzysz na wpływ globalny, czy na konkretny współczynnik? Jeśli np. kluczowy jest „budget”, liczysz głównie na DFBETAS dla tego predyktora:
# załóżmy, że kolumna 1 w X to "budget" (0 to stała)
budget_idx = 1
dfbetas_budget = dfbetas[:, budget_idx]
# prosty próg: |DFBETAS| > 2 / sqrt(n)
threshold_dfbetas = 2 / np.sqrt(n)
influential_on_budget = np.where(np.abs(dfbetas_budget) > threshold_dfbetas)[0]
print("Punkty mocno wpływające na współczynnik budget:", influential_on_budget)
Reszty studentyzowane przydają się przy wychwytywaniu dużych odchyleń względem modelu. Typowe reguły:
|t_i| > 2– obserwacja podejrzana (zwłaszcza przy większych próbach),|t_i| > 3– potencjalny outlier pod względem wartości y.
outliers = np.where(np.abs(studentized_resid) > 3)[0]
print("Potencjalne outliery:", outliers)
Przed wycięciem czegokolwiek zawsze zadaj sobie pytanie biznesowe: czy to błąd danych, czy rzadki, ale realny przypadek? Wycinanie realnych ekstremów bywa gorsze niż ich jawne modelowanie.
Wizualizacja wpływu: wykres lewaru i Cook’s distance w Pythonie
Same liczby bywają mało intuicyjne. Często szybciej zadziała prosty wykres: lewar vs reszty studentyzowane z wielkością punktu zależną od Cook’s distance.
fig, ax = plt.subplots(figsize=(8, 6))
scatter = ax.scatter(
leverage,
studentized_resid,
s=100 * np.sqrt(cooks_d + 1e-8), # wielkość ~ Cook's distance
alpha=0.7
)
ax.axhline(0, color="grey", linestyle="--", linewidth=1)
ax.set_xlabel("Leverage")
ax.set_ylabel("Studentized residuals")
ax.set_title("Leverage vs studentized residuals (rozmiar = Cook's D)")
# orientacyjny próg lewaru
ax.axvline(leverage_threshold, color="red", linestyle="--", linewidth=1)
plt.show()
Zadaj sobie wtedy dwa krótkie pytania:
- które punkty mają jednocześnie wysoki lewar i dużą resztę?
- czy te punkty to specyficzne grupy (np. pojedynczy kraj, grupa VIP-klientów)?
Jeśli tak, często bardziej sensowne jest wydzielenie ich do osobnego modelu niż „siłowe” czyszczenie.
Scenariusze postępowania z obserwacjami wpływowymi
Co zrobić, gdy kilka punktów przeciąga model, a Cook’s distance i DFBETAS migają ostrzegawczo? Zanim usuniesz cokolwiek z danych, odpowiedz sobie: jaki masz cel?
- Prognoza dla „typowych” przypadków – np. chcesz przewidywać miesięczną sprzedaż standardowych klientów.
- Pełny obraz rynku – interesują Cię też outliery, bo niosą kluczową wartość biznesową (np. duże kontrakty).
Od tego uzależniasz strategię:
- Czyszczenie błędów danych – gdy masz ewidentne literówki, błędne jednostki, duplikaty.
- Osobne modele dla segmentów – np. oddzielnie dla „typu klienta = Enterprise” i „SME”.
- Modele odporniejsze na outliery – np. RLM (Robust Linear Model) w statsmodels.
from statsmodels.robust.robust_linear_model import RLM
# robust linear model z funkcją Hampela / Huber
rlm_model = RLM(y, X) # można dodać argument M=sm.robust.norms.HuberT()
rlm_results = rlm_model.fit()
print(rlm_results.summary())
Gdy porównujesz OLS i RLM, zadaj sobie pytanie: czy wnioski o kluczowych współczynnikach są stabilne? Jeśli po przejściu na model odporny na outliery sygnał „budżet działa” nadal jest podobny (znak, rząd wielkości), masz dodatkowe potwierdzenie wiarygodności wniosków.
Porównanie modelu z i bez obserwacji wpływowych
Jedno z najbardziej praktycznych ćwiczeń: dopasuj model na pełnych danych, potem usuń wąską grupę najbardziej wpływowych obserwacji i sprawdź, co się zmieniło. Nie po to, by od razu usuwać, ale by ocenić wrażliwość.
# przyjmijmy próg dla Cook's distance
influential_mask = cooks_d > 4 / n
df_influential = df[influential_mask]
df_clean = df[~influential_mask]
print("Liczba obserwacji wpływowych:", df_influential.shape[0])
# model na pełnych danych
y_full = df["sales"]
X_full = sm.add_constant(df[["budget", "campaigns"]])
results_full = sm.OLS(y_full, X_full).fit()
# model na danych "oczyszczonych"
y_clean = df_clean["sales"]
X_clean = sm.add_constant(df_clean[["budget", "campaigns"]])
results_clean = sm.OLS(y_clean, X_clean).fit()
print("Współczynniki pełny:")
print(results_full.params)
print("Współczynniki bez punktów wpływowych:")
print(results_clean.params)
Kluczowe pytanie: czy wnioski jakościowe się zmieniają? Jeśli po usunięciu kilku punktów:
- znaki współczynników odwracają się,
- a p-value kluczowych zmiennych skaczą „przez granicę” 0.05,
to sygnał, że Twoje wnioski są mocno zależne od konkretnych obserwacji. Wtedy zamiast ukrywać ten fakt w kodzie, lepiej go jawnie komunikować i szukać modelu, który lepiej opisuje różne segmenty.
Łączenie diagnoz: reszty, lewar i heteroskedastyczność razem
Regresja liniowa rzadko „psuje się” tylko w jednym wymiarze. Często masz mieszankę:
- heteroskedastyczność – wariancja rośnie z poziomem zjawiska,
- nienormalne reszty – np. długi ogon po prawej,
- kilka obserwacji o ogromnym lewarze – np. nietypowe budżety marketingowe.
Jak to wszystko poukładać w jedną, praktyczną procedurę? Zadaj sobie sekwencję pytań:
- Czy model w ogóle „łapie” główny wzorzec? — sprawdź wykres y vs y_hat i proste metryki (R², RMSE).
- Czy reszty są losowe względem fitted i głównych predyktorów? — szukasz brakujących nieliniowości lub segmentów.
- Czy wariancja błędu zmienia się z poziomem zjawiska? — heteroskedastyczność i ewentualne transformacje.
- Czy kilka punktów nie „przestawia” całego modelu? — lewar, Cook’s distance, DFBETAS.
Na tej podstawie możesz zbudować realny workflow w Pythonie, zamiast skakać między pojedynczymi testami bez planu.
Propozycja workflow diagnostycznego w Pythonie
Jeśli chcesz poukładać to w logiczną sekwencję kroków, możesz wykorzystać prosty schemat:
- Dopasowanie modelu bazowego (OLS) oraz zapisanie reszt i fitted:
y = df["sales"] X = sm.add_constant(df[["budget", "campaigns"]]) model = sm.OLS(y, X) results = model.fit() residuals = results.resid fitted = results.fittedvalues - Diagnoza wzrokowa: reszty vs fitted, histogram, QQ-plot:
fig, axes = plt.subplots(1, 2, figsize=(12, 5)) sns.scatterplot(x=fitted, y=residuals, ax=axes[0], alpha=0.7) axes[0].axhline(0, color="red", linestyle="--", linewidth=1) axes[0].set_xlabel("Fitted") axes[0].set_ylabel("Residuals") axes[0].set_title("Reszty vs fitted") sm.qqplot(residuals, line="45", ax=axes[1]) axes[1].set_title("QQ-plot reszt") plt.tight_layout() plt.show() - Testy formalne: normalność, heteroskedastyczność, autokorelacja (jeśli dane czasowe):
from statsmodels.stats.diagnostic import het_breuschpagan, het_white from statsmodels.stats.stattools import jarque_bera, durbin_watson bp_stat, bp_pvalue, _, _ = het_breuschpagan(residuals, X) white_stat, white_pvalue, _, _ = het_white(residuals, X) jb_stat, jb_pvalue, skew, kurtosis = jarque_bera(residuals) dw_stat = durbin_watson(residuals) print("Breusch-Pagan p-value:", bp_pvalue) print("White p-value:", white_pvalue) print("Jarque-Bera p-value:", jb_pvalue) print("Durbin-Watson stat:", dw_stat) - Wpływ obserwacji: lewar, Cook’s distance, DFFITS, DFBETAS:
influence = results.get_influence() leverage = influence.hat_matrix_diag cooks_d, _ = influence.cooks_distance dffits, _ = influence.dffits dfbetas = influence.dfbetas studentized_resid = influence.resid_studentized_external - Decyzje korekcyjne:
- jeśli celem są wnioski z p-value – zastosuj
cov_type="HC3", - jeśli dominują problemy skali – rozważ transformacje (np. log) lub WLS,
- jeśli kilka punktów ma ogromny wpływ – sprawdź ich naturę (błąd vs inny segment), przetestuj RLM lub modele segmentowe.
- jeśli celem są wnioski z p-value – zastosuj
Przy każdym kroku wracaj do pytania: czemu ma służyć ten model? Jeśli ma wspierać decyzje o budżecie i prognozach przychodów, ważniejsze od „idealnych” testów jest to, by:
- współczynniki były stabilne na rozsądnym zakresie danych,
- reszty nie zdradzały systematycznych pominięć (np. silnej nieliniowości),
- kilka anomalii nie sterowało polityką dla całego rynku.
Z takim kompasem diagnostyka regresji w Pythonie przestaje być zbiorem oderwanych testów, a staje się narzędziem do świadomego projektowania modelu pod konkretny cel.
Najczęściej zadawane pytania (FAQ)
Po co w ogóle sprawdzać założenia regresji liniowej w Pythonie?
Jeśli używasz regresji tylko do prognozowania i oceniasz model po RMSE/MAE na zbiorze walidacyjnym, naruszenia części założeń mogą być mniej dotkliwe. Pytanie brzmi: czy chcesz tylko przewidywać, czy także interpretować współczynniki i wyciągać wnioski przyczynowo‑skutkowe?
Gdy opierasz decyzje biznesowe lub naukowe na wartościach p i przedziałach ufności, łamanie założeń oznacza fałszywie „istotne” lub „nieistotne” efekty. Heteroskedastyczność, nieliniowość czy kilka obserwacji o ogromnym wpływie mogą całkowicie wypaczyć wnioski typu „kanał marketingowy A nic nie daje”.
Jak sprawdzić reszty regresji liniowej w Pythonie (statsmodels)?
Najpierw zadaj sobie pytanie: co chcesz ocenić – liniowość, normalność, stałą wariancję, czy obecność punktów odstających? Od tego zależy, jakie wykresy i testy zbudujesz na podstawie reszt.
W statsmodels po dopasowaniu modelu masz wprost dostęp do results.resid i results.fittedvalues. Typowe kroki to:
- wykres reszt vs wartości dopasowane (szukasz losowej „chmury”, bez kształtu U/S),
- histogram lub KDE reszt,
- wykres QQ‑plot (funkcja
qqplot), który porównuje rozkład reszt z rozkładem normalnym.
Jeśli na wykresie reszt widzisz strukturę (np. rosnącą „tubę”), to sygnał, że model łamie założenia i trzeba coś zmienić w specyfikacji.
Jak wykryć heteroskedastyczność w regresji liniowej w Pythonie?
Najpierw spójrz na prosty wykres: reszty vs wartości dopasowane. Czy rozrzut reszt rośnie lub maleje wraz z poziomem przewidywanej zmiennej? Jeśli widzisz kształt „lejka”, masz wizualną wskazówkę heteroskedastyczności.
Potem możesz użyć testów formalnych ze statsmodels, np. Breuscha–Pagana lub White’a (moduł statsmodels.stats.api). Wynik w postaci niskiego p‑value sugeruje odrzucenie hipotezy o stałej wariancji błędów. W takiej sytuacji rozważ:
- robust standard errors (np.
HC3), - transformację zmiennej zależnej (np. log),
- modelowanie zmiennej pomocniczej (np. różne wariancje w różnych segmentach).
Kluczowe pytanie: czy zależy ci na poprawnej inferencji, czy tylko na jak najlepszej prognozie?
Kiedy lepiej użyć statsmodels zamiast scikit-learn do regresji?
Zastanów się, czego oczekujesz: pełnej diagnostyki statystycznej czy wygodnego pipeline’u ML? Jeśli chcesz testów istotności, błędów standardowych, miar wpływu i gotowych testów na heteroskedastyczność czy autokorelację, statsmodels będzie naturalnym wyborem.
Scikit‑learn świetnie nadaje się do uczenia wielu modeli, walidacji krzyżowej i integracji w produkcyjnym pipeline’ie, ale nie daje „z pudełka” tabeli współczynników z wartościami p ani miar wpływu. Częsta praktyka: dopasowanie i diagnostyka modelu w statsmodels, a po ustaleniu finalnej specyfikacji odtworzenie jej w scikit‑learn do celów stricte predykcyjnych.
Jak sprawdzić obserwacje wpływowe (Cook’s distance, leverage) w Pythonie?
Najpierw odpowiedz sobie: czy masz podejrzenie, że kilka obserwacji „ciągnie” model w jedną stronę? Na przykład pojedynczy klient z ekstremalnym budżetem reklamowym może mocno zmienić znak współczynnika.
W statsmodels możesz skorzystać z OLSInfluence (np. sm.OLS(y, X).fit().get_influence()). Z obiektu wpływu wyciągniesz:
- Cook’s distance – obserwacje o najwyższych wartościach są potencjalnie najbardziej wpływowe,
- leverage (hat values) – informują, jak „odległy” jest punkt w przestrzeni predyktorów,
- standardized/ studentized residuals – do wykrywania odstających reszt.
Obserwacje o dużym lewarze i wysokim Cook’s distance wymagają osobnego przeglądu: czy to błędne dane, inny segment, czy realny, ale rzadki przypadek?
Co zrobić, gdy zależność nie jest liniowa, a ja chcę zostać przy regresji liniowej?
Sprawdź najpierw wykres reszt vs wartości dopasowane oraz wykresy częściowe (np. plot_partregress_grid). Jeśli widzisz kształty U, S albo inne systematyczne wzory, to sygnał nieliniowości, a niekoniecznie powód, by od razu porzucać regresję liniową.
Masz kilka prostych opcji:
- dodać składniki kwadratowe lub wyższych rzędów (np.
budget,budget**2), - przetestować transformacje typu log, sqrt po stronie zmiennych lub odpowiedzi,
- uwzględnić interakcje (np.
budget * season), jeśli efekt zależy od poziomu innej zmiennej.
Zadaj sobie pytanie: czy po takiej transformacji interpretacja wciąż odpowiada na twoje biznesowe lub naukowe pytanie? Jeśli tak – zlinearyzowana wersja modelu bywa bardzo praktycznym kompromisem.
Kiedy naruszenia założeń regresji są naprawdę krytyczne, a kiedy można je „przełknąć”?
Najlepszy filtr to pytanie: czy na podstawie współczynników podejmujesz decyzje finansowe lub strategiczne? W modelach budżetów marketingowych, polityki cenowej czy oceny skuteczności kampanii błędne wartości p i przedziały ufności oznaczają ryzyko kosztownych pomyłek.
Jeżeli budujesz model wyłącznie po to, by wygenerować prognozę i porównujesz metryki na zbiorze walidacyjnym, drobne odchylenia od normalności reszt mogą być mniej istotne. Krytyczne są jednak sytuacje, w których:
- heteroskedastyczność mocno psuje estymację wariancji,
- kilka punktów o ogromnym lewarze zmienia znak i wielkość kluczowych współczynników,
- masz dane szeregów czasowych z silną autokorelacją reszt, a używasz „gołego” OLS.
W takich przypadkach albo poprawiasz model i błędy standardowe, albo uczciwie rezygnujesz z silnych wniosków przyczynowych.






