{ "cells": [ { "cell_type": "markdown", "id": "d78847d0", "metadata": {}, "source": [ "# Wielokrotna (wieloraka) regresja liniowa\n", "\n", "Na dzisiejszych zajęciach zajmiemy się **wielokrotną regresją liniową** (*multiple linear regression*).\n", "\n", "Będziemy pracować na klasycznym przykładzie dotyczącym edukacji w Stanach Zjednoczonych (Howell, rozdział o regresji). Kod w notatniku ma być bazą do ćwiczeń i ilustracją, jak różne rzeczy możemy policzyć w Pythonie.\n", "\n", "Najważniejsza idea, którą chcemy zrozumieć:\n", "\n", "- w regresji wielokrotnej współczynnik przy danym predyktorze opisuje efekt **„przy kontroli pozostałych predyktorów”**,\n", "- dlatego współczynniki (ich znak i wielkość) mogą się zmieniać, gdy do modelu dodamy kolejną zmienną.\n" ] }, { "cell_type": "markdown", "id": "9de50240", "metadata": {}, "source": [ "## Cele\n", "\n", "Po zajęciach:\n", "\n", "- Umiesz dopasować regresję liniową z kilkoma predyktorami (OLS).\n", "- Rozumiesz, jak interpretować współczynniki „przy kontroli”.\n", "- Potrafisz policzyć p-wartość „na piechotę” dla testów t i F w regresji i porównać z wynikiem funkcji.\n", "- Umiesz policzyć i zinterpretować $R^2$ oraz $R^2_{adj}$.\n", "- Umiesz sprawdzić współliniowość predyktorów (VIF i tolerancja).\n", "- Rozumiesz ideę korelacji cząstkowej / półcząstkowej i ich związek z regresją.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "4a691f00", "metadata": {}, "outputs": [], "source": [ "from pathlib import Path\n", "\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "\n", "from scipy import stats\n", "import statsmodels.api as sm\n", "import statsmodels.formula.api as smf\n", "from statsmodels.stats.outliers_influence import variance_inflation_factor\n", "from statsmodels.graphics.gofplots import ProbPlot\n", "from statsmodels.stats.power import FTestPowerF2\n", "\n", "from IPython.display import display\n", "\n", "SEED = 20260117\n", "rng = np.random.default_rng(SEED)\n", "\n", "sns.set_theme(style=\"whitegrid\")\n", "plt.rcParams[\"figure.figsize\"] = (9, 4)\n", "plt.rcParams[\"figure.dpi\"] = 120\n", "\n", "try:\n", " import ipywidgets as widgets\n", "\n", " WIDGETY_DOSTEPNE = True\n", "except ImportError:\n", " widgets = None\n", " WIDGETY_DOSTEPNE = False" ] }, { "cell_type": "markdown", "id": "28743ab4", "metadata": {}, "source": [ "## Dane\n", "\n", "Wczytajmy dane z pliku `Tab15-1.dat`. Plik jest rozdzielany tabulatorami i ma nagłówek.\n", "\n", "To klasyczny przykład omawiany przez Howella: dane opisują stany USA i wybrane wskaźniki edukacyjne. Nie traktujemy go jako realistycznego modelu przyczynowego polityki edukacyjnej. To raczej wygodny przykład pokazujący, jak bardzo interpretacja relacji między dwiema zmiennymi może się zmienić po uwzględnieniu trzeciej zmiennej.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "fceb2a24", "metadata": {}, "outputs": [], "source": [ "DATA_DIR = Path(\"data\")\n", "\n", "data = pd.read_csv(DATA_DIR / \"Tab15-1.dat\", sep=\"\t\")\n", "\n", "pd.DataFrame(\n", " {\n", " \"liczba_obserwacji\": [data.shape[0]],\n", " \"liczba_zmiennych\": [data.shape[1]],\n", " }\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "4247e283", "metadata": {}, "outputs": [], "source": [ "data.head()" ] }, { "cell_type": "markdown", "id": "696b07e1", "metadata": {}, "source": [ "W każdym wierszu mamy jeden stan USA.\n", "\n", "W tej lekcji będziemy najczęściej korzystać z:\n", "\n", "- `SATcombined` – średni wynik SAT,\n", "- `Expend` – wydatki na edukację,\n", "- `PTratio` – ilu uczniów przypada na jednego nauczyciela,\n", "- `PctSAT` i `LogPctSAT` – popularność testu SAT w danym stanie.\n", "\n", "Ważny kontekst: SAT nie był w każdym stanie równie popularny. Tam, gdzie do testu podchodzi niewielki odsetek uczniów, często są to uczniowie bardziej selektywni i mocniej nastawieni na rekrutację akademicką. To może tworzyć mylący obraz relacji między wydatkami na edukację a średnim wynikiem SAT.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "c874ad4d", "metadata": {}, "outputs": [], "source": [ "cols = [\"SATcombined\", \"Expend\", \"PTratio\", \"PctSAT\", \"LogPctSAT\"]\n", "data[cols].describe().T" ] }, { "cell_type": "markdown", "id": "1b29f2f2", "metadata": {}, "source": [ "## Rozkłady zmiennych i wykresy rozrzutu\n", "\n", "Zanim dopasujemy model, obejrzymy dane. Dla wybranych zmiennych narysujemy:\n", "\n", "- histogram,\n", "- wykres Q–Q,\n", "- wykres rozrzutu z `SATcombined` na osi Y.\n", "\n", "Nie chodzi tu o mechaniczne „zaliczenie diagnostyki”. Chcemy zobaczyć, czy relacje wyglądają mniej więcej liniowo, czy rozkłady nie mają skrajnych obserwacji oraz czy transformacja `PctSAT` do `LogPctSAT` rzeczywiście upraszcza relację z wynikiem SAT.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "6cd6edfb", "metadata": {}, "outputs": [], "source": [ "def qqplot(ax: plt.Axes, sample: np.ndarray, title: str) -> None:\n", " (osm, osr), (slope, intercept, _r) = stats.probplot(sample, dist=\"norm\")\n", " ax.scatter(osm, osr, s=20, alpha=0.8)\n", " x = np.array([osm.min(), osm.max()])\n", " ax.plot(x, intercept + slope * x, color=\"black\", lw=2)\n", " ax.set_title(title)\n", " ax.set_xlabel(\"Teoretyczne kwantyle (N(0,1))\")\n", " ax.set_ylabel(\"Kwantyle z próby\")" ] }, { "cell_type": "code", "execution_count": null, "id": "3de6ea01", "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(nrows=4, ncols=3, figsize=(12, 12))\n", "\n", "# 1) Expend\n", "sns.histplot(data=data, x=\"Expend\", bins=8, ax=axes[0, 0])\n", "axes[0, 0].set_title(\"Expend: histogram\")\n", "axes[0, 0].set_xlabel(\"Expend\")\n", "axes[0, 0].set_ylabel(\"Liczność\")\n", "\n", "qqplot(axes[0, 1], data[\"Expend\"].to_numpy(), \"Expend: Q–Q\")\n", "\n", "sns.regplot(\n", " data=data,\n", " x=\"Expend\",\n", " y=\"SATcombined\",\n", " ax=axes[0, 2],\n", " scatter_kws={\"alpha\": 0.8, \"s\": 30},\n", ")\n", "axes[0, 2].set_title(\"SATcombined vs Expend\")\n", "axes[0, 2].set_xlabel(\"Expend\")\n", "axes[0, 2].set_ylabel(\"SATcombined\")\n", "\n", "# 2) PTratio\n", "sns.histplot(data=data, x=\"PTratio\", bins=8, ax=axes[1, 0])\n", "axes[1, 0].set_title(\"PTratio: histogram\")\n", "axes[1, 0].set_xlabel(\"PTratio\")\n", "axes[1, 0].set_ylabel(\"Liczność\")\n", "\n", "qqplot(axes[1, 1], data[\"PTratio\"].to_numpy(), \"PTratio: Q–Q\")\n", "\n", "sns.regplot(\n", " data=data,\n", " x=\"PTratio\",\n", " y=\"SATcombined\",\n", " ax=axes[1, 2],\n", " scatter_kws={\"alpha\": 0.8, \"s\": 30},\n", ")\n", "axes[1, 2].set_title(\"SATcombined vs PTratio\")\n", "axes[1, 2].set_xlabel(\"PTratio\")\n", "axes[1, 2].set_ylabel(\"SATcombined\")\n", "\n", "# 3) PctSAT\n", "sns.histplot(data=data, x=\"PctSAT\", bins=8, ax=axes[2, 0])\n", "axes[2, 0].set_title(\"PctSAT: histogram\")\n", "axes[2, 0].set_xlabel(\"PctSAT\")\n", "axes[2, 0].set_ylabel(\"Liczność\")\n", "\n", "qqplot(axes[2, 1], data[\"PctSAT\"].to_numpy(), \"PctSAT: Q–Q\")\n", "\n", "sns.regplot(\n", " data=data,\n", " x=\"PctSAT\",\n", " y=\"SATcombined\",\n", " ax=axes[2, 2],\n", " scatter_kws={\"alpha\": 0.8, \"s\": 30},\n", ")\n", "axes[2, 2].set_title(\"SATcombined vs PctSAT\")\n", "axes[2, 2].set_xlabel(\"PctSAT\")\n", "axes[2, 2].set_ylabel(\"SATcombined\")\n", "\n", "# 4) SATcombined + scatter z LogPctSAT\n", "sns.histplot(data=data, x=\"SATcombined\", bins=8, ax=axes[3, 0])\n", "axes[3, 0].set_title(\"SATcombined: histogram\")\n", "axes[3, 0].set_xlabel(\"SATcombined\")\n", "axes[3, 0].set_ylabel(\"Liczność\")\n", "\n", "qqplot(axes[3, 1], data[\"SATcombined\"].to_numpy(), \"SATcombined: Q–Q\")\n", "\n", "sns.regplot(\n", " data=data,\n", " x=\"LogPctSAT\",\n", " y=\"SATcombined\",\n", " ax=axes[3, 2],\n", " scatter_kws={\"alpha\": 0.8, \"s\": 30},\n", ")\n", "axes[3, 2].set_title(\"SATcombined vs LogPctSAT\")\n", "axes[3, 2].set_xlabel(\"LogPctSAT\")\n", "axes[3, 2].set_ylabel(\"SATcombined\")\n", "\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "851baa79", "metadata": {}, "source": [ "Rzeczą, która powinna zwrócić naszą uwagę jest to, że relacja między `Expend` i `SATcombined` wygląda na ujemną. To jest dość nieintuicyjne: spodziewalibyśmy się raczej, że wyższe wydatki na edukację idą w parze z lepszymi wynikami.\n", "\n", "To nie znaczy jeszcze, że większe wydatki „pogarszają” wyniki. Możliwe, że patrzymy na relację zanieczyszczoną trzecią zmienną. W tym przykładzie naturalnym kandydatem jest popularność SAT: jeśli w danym stanie do SAT podchodzi mały, selektywny odsetek uczniów, średni wynik może być wysoki z powodów niezwiązanych bezpośrednio z wydatkami.\n", "\n", "Spróbujmy najpierw policzyć prostą korelację i od razu przeprowadzić test statystyczny dla hipotezy $H_0: \\rho = 0$.\n" ] }, { "cell_type": "markdown", "id": "03b9cc57", "metadata": {}, "source": [ "## Przykład 1: korelacja `SATcombined` i `Expend`\n", "\n", "**Historia.** Czy wydatki na edukację (`Expend`) są liniowo powiązane z wynikiem SAT (`SATcombined`)?\n", "\n", "**Pytanie.** Czy korelacja w populacji jest różna od zera?\n", "\n", "- $H_0: \\rho = 0$\n", "- $H_A: \\rho \\neq 0$ (dwustronnie)\n", "- $\\alpha = 0.05$\n", "\n", "**Statystyka testowa.**\n", "\n", "$$\n", " t = r\\sqrt{\\frac{n-2}{1-r^2}}\n", "$$\n", "\n", "Przy $H_0$: $t \\sim t_{n-2}$.\n", "\n", "**Region krytyczny.** Odrzucamy $H_0$ gdy $|t| > t_{1-\\alpha/2,\\,n-2}$.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "a53d7d96", "metadata": {}, "outputs": [], "source": [ "alpha = 0.05\n", "\n", "x = data[\"Expend\"].to_numpy()\n", "y = data[\"SATcombined\"].to_numpy()\n", "\n", "n = len(data)\n", "df = n - 2\n", "\n", "# r „na piechotę”\n", "r = np.corrcoef(x, y)[0, 1]\n", "\n", "# statystyka t „na piechotę”\n", "t_stat = r * np.sqrt(df / (1 - r**2))\n", "\n", "t_crit = stats.t.isf(alpha / 2, df=df)\n", "p_value = 2 * stats.t.sf(abs(t_stat), df=df)\n", "\n", "print(f\"n = {n}, df = {df}\")\n", "print(f\"r = {r:.4f}\")\n", "print(f\"t = {t_stat:.4f}\")\n", "print(f\"t_crit = ±{t_crit:.4f}\")\n", "print(f\"p (na piechotę) = {p_value:.6f}\")\n", "\n", "# Sprawdzenie funkcją\n", "r_lib, p_lib = stats.pearsonr(x, y)\n", "print()\n", "print(\"scipy.stats.pearsonr:\")\n", "print(f\"r = {r_lib:.4f}, p = {p_lib:.6f}\")" ] }, { "cell_type": "markdown", "id": "ff575000", "metadata": {}, "source": [ "Dla kompletności policzmy 95% CI dla $\\rho$ metodą Fishera.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "39e0296b", "metadata": {}, "outputs": [], "source": [ "z = np.arctanh(r)\n", "se_z = 1 / np.sqrt(n - 3)\n", "z_crit = stats.norm.isf(alpha / 2)\n", "\n", "z_low = z - z_crit * se_z\n", "z_high = z + z_crit * se_z\n", "\n", "r_low = np.tanh(z_low)\n", "r_high = np.tanh(z_high)\n", "\n", "print(f\"95% CI dla rho: [{r_low:.3f}, {r_high:.3f}]\")" ] }, { "cell_type": "markdown", "id": "b7dae1ac", "metadata": {}, "source": [ "Krótka uwaga: gdybyśmy chcieli policzyć korelacje dla wielu par zmiennych, ręczne wywoływanie testu po kolei byłoby uciążliwe.\n", "\n", "Poniżej prosta funkcja tworząca tabelkę wyników (r i p) dla wszystkich par.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "3bcd1e73", "metadata": {}, "outputs": [], "source": [ "def pairwise_pearson_table(df: pd.DataFrame, columns: list[str]) -> pd.DataFrame:\n", " rows: list[dict[str, float | str]] = []\n", "\n", " for i, col1 in enumerate(columns):\n", " for col2 in columns[i + 1 :]:\n", " r_val, p_val = stats.pearsonr(df[col1].to_numpy(), df[col2].to_numpy())\n", " rows.append({\"var1\": col1, \"var2\": col2, \"r\": r_val, \"p\": p_val})\n", "\n", " out = pd.DataFrame(rows)\n", " out = out.sort_values(\"p\", ascending=True).reset_index(drop=True)\n", " return out\n", "\n", "cols_many = [\"Expend\", \"PTratio\", \"Salary\", \"PctSAT\", \"SATcombined\", \"LogPctSAT\"]\n", "\n", "corr_table = pairwise_pearson_table(data, cols_many)\n", "corr_table.head(12)" ] }, { "cell_type": "markdown", "id": "e944adf7", "metadata": {}, "source": [ "## Wielokrotna regresja: kontrolowanie `LogPctSAT`\n", "\n", "Przechodzimy teraz do właściwej idei regresji wielokrotnej.\n", "\n", "Chcemy zbadać związek między wydatkami na edukację i wynikiem SAT **kontrolując** odsetek osób podchodzących do testu. W praktyce pytamy:\n", "\n", "> Jaki jest związek między `Expend` i `SATcombined`, jeśli porównujemy stany podobne pod względem popularności SAT?\n", "\n", "Zaczniemy od dwóch modeli:\n", "\n", "1. Model 1: `SATcombined ~ Expend`\n", "2. Model 2: `SATcombined ~ Expend + LogPctSAT`\n", "\n", "Jeżeli współczynnik przy `Expend` mocno zmieni się po dodaniu `LogPctSAT`, będzie to znak, że prosta korelacja mieszała ze sobą kilka relacji naraz.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "ec9317ff", "metadata": {}, "outputs": [], "source": [ "fit1 = smf.ols(\"SATcombined ~ Expend\", data=data).fit()\n", "fit1.summary()" ] }, { "cell_type": "code", "execution_count": null, "id": "cab7a836", "metadata": {}, "outputs": [], "source": [ "fit2 = smf.ols(\"SATcombined ~ Expend + LogPctSAT\", data=data).fit()\n", "fit2.summary()" ] }, { "cell_type": "markdown", "id": "a83f469a", "metadata": {}, "source": [ "Widzimy, że dodanie `LogPctSAT` radykalnie zmienia obraz relacji. W modelu prostym `Expend` wyglądał na związany ujemnie z `SATcombined`; po kontroli `LogPctSAT` współczynnik staje się dodatni.\n", "\n", "To dobry moment, żeby mocno podkreślić interpretację: współczynnik przy `Expend` w Modelu 2 nie opisuje już surowej relacji między wydatkami i wynikiem SAT. Opisuje relację **przy stałym poziomie `LogPctSAT`**, czyli po porównaniu stanów podobnych pod względem odsetka uczniów podchodzących do SAT.\n", "\n", "Teraz zróbmy test istotności dla współczynnika `Expend` w Modelu 2 krok po kroku.\n" ] }, { "cell_type": "markdown", "id": "0a8bbb08", "metadata": {}, "source": [ "## Przykład 2: test t dla współczynnika w regresji (Model 2)\n", "\n", "**Historia.** Czy `Expend` ma związek z `SATcombined`, gdy kontrolujemy `LogPctSAT`?\n", "\n", "**Pytanie.** Czy współczynnik w populacji jest różny od zera?\n", "\n", "- $H_0: \\beta_{Expend} = 0$\n", "- $H_A: \\beta_{Expend} \\neq 0$\n", "- $\\alpha = 0.05$\n", "\n", "**Statystyka testowa.**\n", "\n", "$$\n", " t = \\frac{b}{SE(b)}\n", "$$\n", "\n", "Przy $H_0$: $t \\sim t_{n-p-1}$ (tu: $p=2$).\n", "\n", "**Region krytyczny.** $|t| > t_{1-\\alpha/2,\\,n-p-1}$.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "38927936", "metadata": {}, "outputs": [], "source": [ "alpha = 0.05\n", "\n", "b = float(fit2.params[\"Expend\"])\n", "se_b = float(fit2.bse[\"Expend\"])\n", "\n", "df_resid = int(fit2.df_resid)\n", "\n", "t_stat = b / se_b\n", "p_value = 2 * stats.t.sf(abs(t_stat), df=df_resid)\n", "\n", "t_crit = stats.t.isf(alpha / 2, df=df_resid)\n", "ci_low = b - t_crit * se_b\n", "ci_high = b + t_crit * se_b\n", "\n", "print(f\"b = {b:.4f}\")\n", "print(f\"SE(b) = {se_b:.4f}\")\n", "print(f\"t = {t_stat:.4f}, df = {df_resid}\")\n", "print(f\"p (na piechotę) = {p_value:.6f}\")\n", "print(f\"95% CI: [{ci_low:.3f}, {ci_high:.3f}]\")\n", "\n", "print()\n", "print(\"Sprawdzenie (statsmodels):\")\n", "print(f\"t = {fit2.tvalues['Expend']:.4f}, p = {fit2.pvalues['Expend']:.6f}\")" ] }, { "cell_type": "markdown", "id": "caa68d1a", "metadata": {}, "source": [ "## Interpretacja wielokrotnej regresji (reszty)\n", "\n", "Najważniejsza intuicja: „kontrolować zmienną” znaczy usunąć z rozważanej relacji tę część zmienności, którą da się przypisać zmiennej kontrolowanej.\n", "\n", "Żeby zobaczyć to ręcznie:\n", "\n", "1. dopasujemy `SATcombined ~ LogPctSAT` i weźmiemy reszty,\n", "2. dopasujemy `Expend ~ LogPctSAT` i weźmiemy reszty,\n", "3. dopasujemy regresję reszt `resid_sat ~ resid_expend`.\n", "\n", "Reszty `resid_sat` to ta część wyniku SAT, której nie przewiduje `LogPctSAT`. Reszty `resid_expend` to ta część wydatków, której nie przewiduje `LogPctSAT`. Jeżeli skorelujemy lub zregresujemy te dwie „oczyszczone” części, dostaniemy tę samą logikę, którą model wielokrotny wykonuje automatycznie.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "9be41111", "metadata": {}, "outputs": [], "source": [ "sat_on_log = smf.ols(\"SATcombined ~ LogPctSAT\", data=data).fit()\n", "expend_on_log = smf.ols(\"Expend ~ LogPctSAT\", data=data).fit()\n", "\n", "resid_sat = sat_on_log.resid\n", "resid_expend = expend_on_log.resid\n", "\n", "df_resid = pd.DataFrame({\"resid_sat\": resid_sat, \"resid_expend\": resid_expend})\n", "fit_resid = smf.ols(\"resid_sat ~ resid_expend\", data=df_resid).fit()\n", "\n", "pd.DataFrame(\n", " {\n", " \"źródło\": [\"regresja na resztach\", \"Model 2\"],\n", " \"b dla Expend\": [\n", " fit_resid.params[\"resid_expend\"],\n", " fit2.params[\"Expend\"],\n", " ],\n", " }\n", ")" ] }, { "cell_type": "markdown", "id": "f6f15aab", "metadata": {}, "source": [ "### Widżet: co znaczy „kontrolować” zmienną?\n", "\n", "Poniższy widżet pokazuje ideę kontroli statystycznej na prostych danych syntetycznych.\n", "\n", "Lewy panel pokazuje surową relację między $X$ i $Y$. Prawy panel pokazuje relację między resztami: z $X$ usuwamy część przewidywaną przez $Z$, z $Y$ usuwamy część przewidywaną przez $Z$, a potem patrzymy, czy te dwie „oczyszczone” części nadal są powiązane.\n", "\n", "Osie są stałe, więc zmianę relacji widać bez autoskalowania wykresu." ] }, { "cell_type": "code", "execution_count": null, "id": "c90785f8", "metadata": {}, "outputs": [], "source": [ "control_widget_rng = np.random.default_rng(SEED + 120)\n", "control_z_base = control_widget_rng.normal(size=120)\n", "control_x_noise = control_widget_rng.normal(scale=0.8, size=120)\n", "control_y_noise = control_widget_rng.normal(scale=0.8, size=120)\n", "\n", "\n", "def standaryzuj_array(values: np.ndarray) -> np.ndarray:\n", " return (values - values.mean()) / values.std(ddof=1)\n", "\n", "\n", "def residualizuj_array(values: np.ndarray, control: np.ndarray) -> np.ndarray:\n", " x_design = sm.add_constant(control)\n", " fit = sm.OLS(values, x_design).fit()\n", " return np.asarray(fit.resid)\n", "\n", "\n", "def rysuj_widget_kontrolowania(\n", " efekt_x_na_y: float = 0.6,\n", " zwiazek_z_x: float = 0.8,\n", " zwiazek_z_y: float = -1.2,\n", ") -> None:\n", " z = standaryzuj_array(control_z_base)\n", " x = standaryzuj_array(zwiazek_z_x * z + control_x_noise)\n", " y = standaryzuj_array(efekt_x_na_y * x + zwiazek_z_y * z + control_y_noise)\n", "\n", " x_res = standaryzuj_array(residualizuj_array(x, z))\n", " y_res = standaryzuj_array(residualizuj_array(y, z))\n", "\n", " raw_slope, raw_intercept, raw_r, raw_p, _ = stats.linregress(x, y)\n", " res_slope, res_intercept, res_r, res_p, _ = stats.linregress(x_res, y_res)\n", "\n", " fig, axes = plt.subplots(ncols=2, figsize=(11, 4.4), sharex=True, sharey=True)\n", "\n", " axes[0].scatter(x, y, s=35, alpha=0.72, color=\"#4477AA\", edgecolor=\"white\", linewidth=0.4)\n", " x_line = np.linspace(-3.2, 3.2, 100)\n", " axes[0].plot(x_line, raw_intercept + raw_slope * x_line, color=\"black\", lw=2)\n", " axes[0].set_title(\"Surowa relacja X–Y\")\n", " axes[0].set_xlabel(\"X\")\n", " axes[0].set_ylabel(\"Y\")\n", " axes[0].text(\n", " 0.04,\n", " 0.96,\n", " f\"r = {raw_r:.2f}\\nb = {raw_slope:.2f}\\np = {raw_p:.3f}\",\n", " transform=axes[0].transAxes,\n", " va=\"top\",\n", " bbox={\"boxstyle\": \"round,pad=0.4\", \"fc\": \"white\", \"ec\": \"0.35\"},\n", " )\n", "\n", " axes[1].scatter(x_res, y_res, s=35, alpha=0.72, color=\"#228833\", edgecolor=\"white\", linewidth=0.4)\n", " axes[1].plot(x_line, res_intercept + res_slope * x_line, color=\"black\", lw=2)\n", " axes[1].set_title(\"Po kontroli Z: reszty X i Y\")\n", " axes[1].set_xlabel(\"X_res\")\n", " axes[1].text(\n", " 0.04,\n", " 0.96,\n", " f\"r = {res_r:.2f}\\nb = {res_slope:.2f}\\np = {res_p:.3f}\",\n", " transform=axes[1].transAxes,\n", " va=\"top\",\n", " bbox={\"boxstyle\": \"round,pad=0.4\", \"fc\": \"white\", \"ec\": \"0.35\"},\n", " )\n", "\n", " for ax in axes:\n", " ax.axhline(0, color=\"0.7\", lw=1)\n", " ax.axvline(0, color=\"0.7\", lw=1)\n", " ax.set_xlim(-3.2, 3.2)\n", " ax.set_ylim(-3.2, 3.2)\n", "\n", " fig.suptitle(\"Kontrolowanie zmiennej jako praca na resztach\", y=1.03)\n", " plt.tight_layout()\n", " plt.show()\n", "\n", "\n", "if WIDGETY_DOSTEPNE:\n", " widget_style = {\"description_width\": \"initial\"}\n", " controls_control = {\n", " \"efekt_x_na_y\": widgets.FloatSlider(\n", " value=0.6,\n", " min=-1.5,\n", " max=1.5,\n", " step=0.1,\n", " description=\"bezpośredni efekt X→Y\",\n", " continuous_update=False,\n", " style=widget_style,\n", " ),\n", " \"zwiazek_z_x\": widgets.FloatSlider(\n", " value=0.8,\n", " min=-1.5,\n", " max=1.5,\n", " step=0.1,\n", " description=\"związek Z→X\",\n", " continuous_update=False,\n", " style=widget_style,\n", " ),\n", " \"zwiazek_z_y\": widgets.FloatSlider(\n", " value=-1.2,\n", " min=-1.5,\n", " max=1.5,\n", " step=0.1,\n", " description=\"związek Z→Y\",\n", " continuous_update=False,\n", " style=widget_style,\n", " ),\n", " }\n", "\n", " output_control = widgets.interactive_output(rysuj_widget_kontrolowania, controls_control)\n", " display(widgets.VBox([widgets.HBox(list(controls_control.values())), output_control]))\n", "else:\n", " rysuj_widget_kontrolowania()\n" ] }, { "cell_type": "markdown", "id": "6bea45c2", "metadata": {}, "source": [ "Możemy też myśleć o regresji wielokrotnej jako o tworzeniu nowej zmiennej $\\hat{Y}$ będącej najlepszą liniową kombinacją predyktorów.\n", "\n", "Dla Modelu 2 ta kombinacja ma postać:\n", "\n", "$$\n", "\\widehat{SAT} = b_0 + b_1 \\cdot Expend + b_2 \\cdot LogPctSAT.\n", "$$\n", "\n", "„Najlepsza” znaczy tutaj: taka, która minimalizuje sumę kwadratów reszt. Wtedy $R = corr(Y, \\hat{Y})$ i $R^2 = corr(Y, \\hat{Y})^2$.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "42e146cc", "metadata": {}, "outputs": [], "source": [ "y = data[\"SATcombined\"]\n", "y_hat = fit2.fittedvalues\n", "\n", "r_y_yhat = np.corrcoef(y, y_hat)[0, 1]\n", "\n", "pd.DataFrame(\n", " {\n", " \"miara\": [\"corr(Y, Y_hat)\", \"corr(Y, Y_hat)^2\", \"R^2 z modelu\"],\n", " \"wartość\": [r_y_yhat, r_y_yhat**2, fit2.rsquared],\n", " }\n", ")" ] }, { "cell_type": "markdown", "id": "652b5460", "metadata": {}, "source": [ "## Tolerancja i VIF\n", "\n", "W regresji wielokrotnej predyktory mogą być ze sobą skorelowane. Nie jest to automatycznie błąd, ale silna współliniowość sprawia, że współczynniki stają się niestabilne: mała zmiana danych może mocno zmienić estymaty, a błędy standardowe rosną.\n", "\n", "VIF (*Variance Inflation Factor*) mówi, jak bardzo wariancja estymatora współczynnika jest „napompowana” przez współliniowość z innymi predyktorami. Tolerancja to po prostu $1/VIF$.\n", "\n", "Dopasujmy model z trzema predyktorami: `Expend`, `LogPctSAT`, `PTratio`. Następnie policzymy VIF i tolerancję.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "1fc3ecea", "metadata": {}, "outputs": [], "source": [ "fit3 = smf.ols(\"SATcombined ~ Expend + LogPctSAT + PTratio\", data=data).fit()\n", "\n", "predictors = [\"Expend\", \"LogPctSAT\", \"PTratio\"]\n", "X = sm.add_constant(data[predictors])\n", "\n", "rows = []\n", "for idx, name in enumerate(X.columns):\n", " if name == \"const\":\n", " continue\n", "\n", " vif = variance_inflation_factor(X.to_numpy(), idx)\n", " rows.append({\"predyktor\": name, \"VIF\": vif, \"tolerancja\": 1 / vif})\n", "\n", "pd.DataFrame(rows)" ] }, { "cell_type": "markdown", "id": "fe684988", "metadata": {}, "source": [ "## Standaryzowane współczynniki regresji\n", "\n", "Częste pytanie brzmi: który predyktor jest „ważniejszy”? Nie możemy odpowiedzieć na nie, porównując surowe współczynniki regresji, bo zależą od jednostek pomiaru. Jeden predyktor może być mierzony w dolarach, inny w procentach, a jeszcze inny w logarytmach.\n", "\n", "Jednym prostym rozwiązaniem jest standaryzacja zmiennych (z-score) i dopasowanie modelu na standaryzowanych danych. Wtedy współczynniki mówią o zmianie w $Y$ (w SD), gdy $X$ wzrośnie o 1 SD.\n", "\n", "To nadal nie jest pełna teoria „ważności predyktorów”, ale daje uczciwszy punkt startu niż porównywanie surowych jednostek.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "2803f44c", "metadata": {}, "outputs": [], "source": [ "def zscore(x: pd.Series) -> pd.Series:\n", " return (x - x.mean()) / x.std(ddof=1)\n", "\n", "std = pd.DataFrame(\n", " {\n", " \"SATcombined_z\": zscore(data[\"SATcombined\"]),\n", " \"Expend_z\": zscore(data[\"Expend\"]),\n", " \"LogPctSAT_z\": zscore(data[\"LogPctSAT\"]),\n", " \"PTratio_z\": zscore(data[\"PTratio\"]),\n", " }\n", ")\n", "\n", "fit3_std = smf.ols(\n", " \"SATcombined_z ~ Expend_z + LogPctSAT_z + PTratio_z\",\n", " data=std,\n", ").fit()\n", "\n", "fit3_std.params.rename(\"beta standaryzowane\").to_frame()\n" ] }, { "cell_type": "markdown", "id": "bf475da8", "metadata": {}, "source": [ "## Skorygowane $R^2$ (adjusted $R^2$)\n", "\n", "Zwykłe $R^2$ prawie zawsze rośnie, gdy dodajemy predyktory (nawet słabe).\n", "\n", "Dlatego często pokazuje się też wersję skorygowaną:\n", "\n", "$$\n", "R^2_{adj} = 1 - \\frac{(1-R^2)(N-1)}{N - p - 1},\n", "$$\n", "\n", "gdzie $N$ to liczba obserwacji, a $p$ to liczba predyktorów.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "214cf58f", "metadata": {}, "outputs": [], "source": [ "n = int(fit3.nobs)\n", "p = 3\n", "\n", "r2 = fit3.rsquared\n", "r2_adj_manual = 1 - ((1 - r2) * (n - 1)) / (n - p - 1)\n", "\n", "pd.DataFrame(\n", " {\n", " \"miara\": [\"R^2\", \"R^2_adj (manual)\", \"R^2_adj (statsmodels)\"],\n", " \"wartość\": [r2, r2_adj_manual, fit3.rsquared_adj],\n", " }\n", ")" ] }, { "cell_type": "markdown", "id": "0d08c8ca", "metadata": {}, "source": [ "## Przykład 3: test F dla $R^2$ (Model 3)\n", "\n", "Testujemy:\n", "\n", "- $H_0: R^2 = 0$ (równoważnie: wszystkie nachylenia w populacji są 0)\n", "- $H_A: R^2 > 0$\n", "\n", "Statystyka:\n", "\n", "$$\n", "F = \\frac{(N-p-1)R^2}{p(1-R^2)}\n", "$$\n", "\n", "Pod $H_0$: $F \\sim F(p, N-p-1)$.\n", "\n", "Region krytyczny: $F > F_{1-\\alpha,\\,p,\\,N-p-1}$.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "c522ff1e", "metadata": {}, "outputs": [], "source": [ "alpha = 0.05\n", "\n", "n = int(fit3.nobs)\n", "p = 3\n", "\n", "r2 = fit3.rsquared\n", "\n", "f_stat_manual = ((n - p - 1) * r2) / (p * (1 - r2))\n", "f_crit = stats.f.isf(alpha, dfn=p, dfd=n - p - 1)\n", "p_manual = stats.f.sf(f_stat_manual, dfn=p, dfd=n - p - 1)\n", "\n", "pd.DataFrame(\n", " {\n", " \"miara\": [\"R^2\", \"F (manual)\", \"F_crit\", \"p (manual)\", \"F (statsmodels)\", \"p (statsmodels)\"],\n", " \"wartość\": [r2, f_stat_manual, f_crit, p_manual, fit3.fvalue, fit3.f_pvalue],\n", " }\n", ")" ] }, { "cell_type": "markdown", "id": "ca42feec", "metadata": {}, "source": [ "## Założenia regresji i wykresy diagnostyczne\n", "\n", "Diagnostyka regresji nie polega na jednym teście. Chodzi raczej o sprawdzenie, czy model nie łamie założeń w sposób, który zmienia interpretację wyników.\n", "\n", "Użyjemy układu zgodnego z oficjalnym przykładem `statsmodels` dla diagnostyki regresji liniowej:\n", "\n", "- **Residuals vs Fitted**: czy w resztach została nieliniowa struktura?\n", "- **Normal Q-Q**: czy rozkład reszt jest z grubsza zgodny z normalnością?\n", "- **Scale-Location**: czy wariancja reszt jest podobna w całym zakresie dopasowanych wartości?\n", "- **Residuals vs Leverage**: czy pojedyncze obserwacje mają duży wpływ na dopasowanie modelu?\n", "\n", "To są wykresy diagnostyczne. Nie „unieważniają” automatycznie modelu, ale pokazują, gdzie interpretacja powinna być ostrożniejsza." ] }, { "cell_type": "code", "execution_count": null, "id": "33e06845", "metadata": {}, "outputs": [], "source": [ "def regression_diagnostic_plots(fit) -> None:\n", " fitted = fit.fittedvalues\n", " residuals = fit.resid\n", "\n", " influence = fit.get_influence()\n", " standardized_residuals = influence.resid_studentized_internal\n", " leverage = influence.hat_matrix_diag\n", " cooks_distance = influence.cooks_distance[0]\n", "\n", " fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(11, 9))\n", "\n", " sns.residplot(\n", " x=fitted,\n", " y=residuals,\n", " lowess=True,\n", " scatter_kws={\"alpha\": 0.65, \"s\": 35},\n", " line_kws={\"color\": \"#CC6677\", \"lw\": 2},\n", " ax=axes[0, 0],\n", " )\n", " axes[0, 0].axhline(0, color=\"black\", lw=1, alpha=0.6)\n", " axes[0, 0].set_title(\"Residuals vs Fitted\")\n", " axes[0, 0].set_xlabel(\"Wartości dopasowane\")\n", " axes[0, 0].set_ylabel(\"Reszty\")\n", "\n", " qq = ProbPlot(standardized_residuals)\n", " qq.qqplot(line=\"45\", alpha=0.65, lw=1, ax=axes[0, 1])\n", " axes[0, 1].set_title(\"Normal Q-Q\")\n", " axes[0, 1].set_xlabel(\"Teoretyczne kwantyle\")\n", " axes[0, 1].set_ylabel(\"Standaryzowane reszty\")\n", "\n", " sqrt_abs_resid = np.sqrt(np.abs(standardized_residuals))\n", " axes[1, 0].scatter(fitted, sqrt_abs_resid, alpha=0.65, s=35, color=\"#4477AA\")\n", " sns.regplot(\n", " x=fitted,\n", " y=sqrt_abs_resid,\n", " scatter=False,\n", " lowess=True,\n", " ci=False,\n", " line_kws={\"color\": \"#CC6677\", \"lw\": 2},\n", " ax=axes[1, 0],\n", " )\n", " axes[1, 0].set_title(\"Scale-Location\")\n", " axes[1, 0].set_xlabel(\"Wartości dopasowane\")\n", " axes[1, 0].set_ylabel(r\"$\\sqrt{|\\mathrm{standaryzowane\\ reszty}|}$\")\n", "\n", " axes[1, 1].scatter(leverage, standardized_residuals, alpha=0.65, s=35, color=\"#4477AA\")\n", " sns.regplot(\n", " x=leverage,\n", " y=standardized_residuals,\n", " scatter=False,\n", " lowess=True,\n", " ci=False,\n", " line_kws={\"color\": \"#CC6677\", \"lw\": 2},\n", " ax=axes[1, 1],\n", " )\n", " axes[1, 1].axhline(0, color=\"black\", lw=1, alpha=0.6)\n", "\n", " n = int(fit.nobs)\n", " p = len(fit.params)\n", " leverage_threshold = 2 * p / n\n", " axes[1, 1].axvline(\n", " leverage_threshold,\n", " color=\"0.35\",\n", " ls=\"--\",\n", " lw=1,\n", " label=\"2p/n\",\n", " )\n", "\n", " x_max = max(float(leverage.max()) + 0.01, leverage_threshold * 1.15)\n", " x_values = np.linspace(0.001, x_max, 200)\n", " for cook_value in [0.5, 1.0]:\n", " cook_line = np.sqrt((cook_value * p * (1 - x_values)) / x_values)\n", " axes[1, 1].plot(\n", " x_values,\n", " cook_line,\n", " color=\"0.45\",\n", " ls=\"--\",\n", " lw=1,\n", " label=f\"Cook D={cook_value:g}\",\n", " )\n", " axes[1, 1].plot(x_values, -cook_line, color=\"0.45\", ls=\"--\", lw=1)\n", "\n", " axes[1, 1].set_xlim(0, x_max)\n", " axes[1, 1].set_title(\"Residuals vs Leverage\")\n", " axes[1, 1].set_xlabel(\"Dźwignia\")\n", " axes[1, 1].set_ylabel(\"Standaryzowane reszty\")\n", " axes[1, 1].legend(loc=\"best\")\n", "\n", " influential_idx = np.argsort(cooks_distance)[-3:]\n", " for idx in influential_idx:\n", " axes[1, 1].annotate(\n", " str(idx),\n", " xy=(leverage[idx], standardized_residuals[idx]),\n", " xytext=(5, 5),\n", " textcoords=\"offset points\",\n", " color=\"#882255\",\n", " )\n", "\n", " fig.suptitle(\"Diagnostyka modelu regresji\", y=1.01)\n", " fig.tight_layout()\n", " plt.show()\n", "\n", "\n", "regression_diagnostic_plots(fit3)" ] }, { "cell_type": "markdown", "id": "1bfcc580", "metadata": {}, "source": [ "Dodatkowo warto obejrzeć wykresy częściowej regresji. One nie zastępują klasycznej diagnostyki reszt, ale pomagają zobaczyć, jak każdy predyktor wygląda po uwzględnieniu pozostałych predyktorów w modelu." ] }, { "cell_type": "code", "execution_count": null, "id": "6f6d4d34", "metadata": {}, "outputs": [], "source": [ "fig = sm.graphics.plot_partregress_grid(fit3)\n", "fig.suptitle(\"Wykresy częściowej regresji\", y=1.02)\n", "fig.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "930f1ab6", "metadata": {}, "source": [ "## Wielkość próbki i moc testu dla $R^2$\n", "\n", "Dla testu globalnego w regresji wielokrotnej używamy miary Cohena:\n", "\n", "$$\n", "f^2 = \\frac{R^2}{1-R^2}\n", "$$\n", "\n", "Nie musimy implementować mocy testu F ręcznie. `statsmodels.stats.power.FTestPowerF2` liczy moc i rozwiązuje równanie mocy dla testu F opisanego przez $f^2$, licznikowe stopnie swobody i mianownikowe stopnie swobody." ] }, { "cell_type": "code", "execution_count": null, "id": "998f91a4", "metadata": {}, "outputs": [], "source": [ "def f2_from_r2(r2: float) -> float:\n", " return r2 / (1 - r2)\n", "\n", "\n", "def required_sample_size_for_r2(\n", " num_predictors: int,\n", " r2: float,\n", " alpha: float,\n", " target_power: float,\n", ") -> int:\n", " f2 = f2_from_r2(r2)\n", " df_denom = FTestPowerF2().solve_power(\n", " effect_size=f2,\n", " df_num=num_predictors,\n", " alpha=alpha,\n", " power=target_power,\n", " )\n", " return int(np.ceil(df_denom + num_predictors + 1))\n", "\n", "\n", "alpha = 0.05\n", "target_power = 0.80\n", "power_calculator = FTestPowerF2()\n", "\n", "scenarios = pd.DataFrame(\n", " {\n", " \"liczba_predyktorów\": [3, 1, 5],\n", " \"R^2\": [0.20, 0.09, 0.09],\n", " }\n", ")\n", "\n", "rows = []\n", "for _, scenario in scenarios.iterrows():\n", " num_predictors = int(scenario[\"liczba_predyktorów\"])\n", " r2_value = scenario[\"R^2\"]\n", " f2_value = f2_from_r2(r2_value)\n", " n_required = required_sample_size_for_r2(\n", " num_predictors=num_predictors,\n", " r2=r2_value,\n", " alpha=alpha,\n", " target_power=target_power,\n", " )\n", " df_denom = n_required - num_predictors - 1\n", " achieved_power = power_calculator.power(\n", " effect_size=f2_value,\n", " df_num=num_predictors,\n", " df_denom=df_denom,\n", " alpha=alpha,\n", " )\n", " rows.append(\n", " {\n", " \"liczba_predyktorów\": num_predictors,\n", " \"R^2\": r2_value,\n", " \"f^2\": f2_value,\n", " \"N dla mocy 0.80\": n_required,\n", " \"uzyskana moc\": achieved_power,\n", " }\n", " )\n", "\n", "pd.DataFrame(rows)" ] }, { "cell_type": "markdown", "id": "e071bd68", "metadata": {}, "source": [ "## Cząstkowa i półcząstkowa korelacja\n", "\n", "Korelacja cząstkowa $r_{XY\\cdot Z}$ to korelacja między resztami:\n", "\n", "- $X_{res}$ z modelu $X \\sim Z$,\n", "- $Y_{res}$ z modelu $Y \\sim Z$.\n", "\n", "Inaczej mówiąc: pytamy, czy ta część $X$, której nie przewiduje $Z$, jest związana z tą częścią $Y$, której nie przewiduje $Z$.\n", "\n", "Test istotności (dla jednej zmiennej kontrolowanej, $k=1$):\n", "\n", "- $t = r\\sqrt{\\frac{n-k-2}{1-r^2}}$, df = $n-k-2$.\n", "\n", "Korelacja półcząstkowa często oznacza $corr(Y, X_{res})$: usuwamy wpływ $Z$ tylko z predyktora, ale nie ze zmiennej zależnej. To rozróżnienie jest ważne, bo kwadraty tych korelacji odpowiadają nieco innym pytaniom o wyjaśnianą zmienność." ] }, { "cell_type": "code", "execution_count": null, "id": "b2c4c6c6", "metadata": {}, "outputs": [], "source": [ "x = data[\"Expend\"].to_numpy()\n", "y = data[\"SATcombined\"].to_numpy()\n", "\n", "x_res = smf.ols(\"Expend ~ LogPctSAT\", data=data).fit().resid.to_numpy()\n", "y_res = smf.ols(\"SATcombined ~ LogPctSAT\", data=data).fit().resid.to_numpy()\n", "\n", "r_partial = np.corrcoef(x_res, y_res)[0, 1]\n", "\n", "k = 1\n", "n = len(data)\n", "df = n - k - 2\n", "\n", "t_stat = r_partial * np.sqrt(df / (1 - r_partial**2))\n", "p_value = 2 * stats.t.sf(abs(t_stat), df=df)\n", "\n", "r_semipartial = np.corrcoef(y, x_res)[0, 1]\n", "\n", "print(f\"r_partial = {r_partial:.4f}\")\n", "print(f\"t = {t_stat:.4f}, df = {df}\")\n", "print(f\"p = {p_value:.6f}\")\n", "\n", "print()\n", "print(\"r_semipartial =\", round(float(r_semipartial), 4))\n", "print(\"r_partial^2 =\", round(float(r_partial**2), 4))\n" ] }, { "cell_type": "markdown", "id": "040d97bf", "metadata": {}, "source": [ "## Porównywanie modeli\n", "\n", "W praktyce często mamy więcej niż jeden sensowny model. Jedna teoria może sugerować zestaw predyktorów A, druga zestaw predyktorów B. Samo porównanie $R^2$ bywa mylące, bo dodawanie predyktorów prawie zawsze poprawia dopasowanie w próbie.\n", "\n", "### Modele zagnieżdżone\n", "\n", "Porównamy:\n", "\n", "- model pełny: `SATcombined ~ Salary + PTratio + LogPctSAT + Expend`\n", "- model zredukowany: `SATcombined ~ LogPctSAT + Expend`\n", "\n", "Modele są zagnieżdżone, bo model zredukowany jest szczególnym przypadkiem modelu pełnego. Testujemy, czy dodanie `Salary` i `PTratio` istotnie poprawia dopasowanie.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "5382bca2", "metadata": {}, "outputs": [], "source": [ "fit_full = smf.ols(\"SATcombined ~ Salary + PTratio + LogPctSAT + Expend\", data=data).fit()\n", "fit_reduced = smf.ols(\"SATcombined ~ LogPctSAT + Expend\", data=data).fit()\n", "\n", "f_stat, p_val, df_diff = fit_full.compare_f_test(fit_reduced)\n", "\n", "pd.DataFrame(\n", " {\n", " \"porównanie\": [\"model pełny vs model zredukowany\"],\n", " \"F\": [f_stat],\n", " \"df_diff\": [df_diff],\n", " \"p\": [p_val],\n", " }\n", ")" ] }, { "cell_type": "markdown", "id": "9295dd63", "metadata": {}, "source": [ "### AIC dla modeli niezagnieżdżonych\n", "\n", "Jeżeli modele nie są zagnieżdżone, klasyczny test porównujący modele nie jest już tak prosty. Możemy wtedy użyć kryteriów informacyjnych, np. AIC.\n", "\n", "Porównamy dwa modele różniące się tym, czy używamy `LogPctSAT` czy `PctSAT`.\n", "\n", "Mniejsze AIC sugeruje lepszy kompromis między dopasowaniem i złożonością. Nie jest to „test istotności” w takim sensie jak p-wartość; to narzędzie do porównywania modeli.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "17d67f31", "metadata": {}, "outputs": [], "source": [ "fit_log = smf.ols(\"SATcombined ~ LogPctSAT + Expend\", data=data).fit()\n", "fit_pct = smf.ols(\"SATcombined ~ PctSAT + Expend\", data=data).fit()\n", "\n", "compare = pd.DataFrame(\n", " {\n", " \"model\": [\"LogPctSAT + Expend\", \"PctSAT + Expend\"],\n", " \"AIC\": [fit_log.aic, fit_pct.aic],\n", " \"BIC\": [fit_log.bic, fit_pct.bic],\n", " \"R^2\": [fit_log.rsquared, fit_pct.rsquared],\n", " }\n", ")\n", "\n", "compare\n" ] }, { "cell_type": "markdown", "id": "3e7dc7df", "metadata": {}, "source": [ "## Podsumowanie\n", "\n", "- W regresji wielokrotnej współczynniki są **warunkowe** („przy kontroli”).\n", "- Kontrolowanie można rozumieć przez reszty: usuwamy z dwóch zmiennych część przewidywaną przez trzecią zmienną.\n", "- Testy t i F możemy policzyć ręcznie z rozkładów pod $H_0$.\n", "- VIF/tolerancja diagnozują współliniowość.\n", "- Korelacje cząstkowe są ściśle powiązane z interpretacją współczynników regresji.\n", "- Przy planowaniu badań można liczyć moc testu dla $R^2$.\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.14.2" } }, "nbformat": 4, "nbformat_minor": 5 }