{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "daf9dd14",
   "metadata": {},
   "source": [
    "# Korelacja i regresja liniowa\n",
    "\n",
    "Na tym wykładzie zajmujemy się **korelacją Pearsona** oraz **prostą regresją liniową z jednym predyktorem**.\n",
    "\n",
    "W praktyce te dwa tematy bardzo często idą razem:\n",
    "\n",
    "- korelacja odpowiada na pytanie: *„czy jest liniowy związek i jak silny?”*,\n",
    "- regresja odpowiada na pytanie: *„o ile zmieni się przeciętne $Y$, gdy $X$ wzrośnie o 1 jednostkę?”* oraz pozwala budować **predykcje**.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c94b5b0",
   "metadata": {},
   "source": [
    "## Cele\n",
    "\n",
    "Po zajęciach:\n",
    "\n",
    "- Umiesz narysować i opisać wykres punktowy (*scatterplot*).\n",
    "- Rozumiesz, czym różni się **kowariancja** od **korelacji**.\n",
    "- Umiesz policzyć współczynnik korelacji $r$ oraz zbudować dla niego test istotności.\n",
    "- Umiesz dopasować prostą regresję liniową, odczytać i zinterpretować *intercept* i *slope*.\n",
    "- Potrafisz policzyć p-wartość „na piechotę” (z rozkładu pod $H_0$) i porównać ją z wynikiem funkcji bibliotecznej.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78b6aaad",
   "metadata": {},
   "source": [
    "## Ważne uwagi (interpretacja)\n",
    "\n",
    "- **Korelacja nie implikuje przyczynowości**: nawet duże $|r|$ nie oznacza, że $X$ „powoduje” $Y$.\n",
    "- W analizie relacji między zmiennymi **zawsze zaczynamy od wykresu** (łatwo przeoczyć nieliniowość, outliery, grupy).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d805708b",
   "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",
    "import statsmodels.formula.api as smf\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import clear_output, display\n",
    "from scipy import stats\n",
    "from statsmodels.stats.outliers_influence import OLSInfluence\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"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ed510c0",
   "metadata": {},
   "source": [
    "W tej części kursu będziemy zajmować się korelacją oraz regresją. Na początek kilka uwag wstępnych.\n",
    "\n",
    "W statystyce zasadniczo odróżnia się te dwa pojęcia na poziomie konceptualnym:\n",
    "\n",
    "**Regresja** – metoda pozwalająca na zbadanie związku pomiędzy zmiennymi i wykorzystanie tej wiedzy do przewidywania nieznanych wartości jednych wielkości na podstawie znajomości wartości innych. Poszukuje się związku między jedną (lub więcej) zmienną objaśniającą (*predyktorem*) $X$ a zmienną objaśnianą (*wynikiem*) $Y$.\n",
    "\n",
    "W **regresji** często myślimy o $X$ jako o zmiennej „ustawionej” (kontrolowanej) przez badacza.\n",
    "\n",
    "W **korelacji** obie zmienne traktujemy jako zmienne losowe.\n",
    "\n",
    "W praktyce różnica ta bywa zatarta, ale warto o niej pamiętać.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7c2a465",
   "metadata": {},
   "source": [
    "# Wykres punktowy (*scatterplot*)\n",
    "\n",
    "Na początek przyjrzyjmy się sposobowi wizualizacji relacji między dwiema zmiennymi.\n",
    "\n",
    "Wykres punktowy świetnie nadaje się do ilustrowania tego typu danych: każdy punkt to jedna obserwacja (np. osoba, kraj, miasto)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31203ba6",
   "metadata": {},
   "source": [
    "## Podstawy\n",
    "\n",
    "Poniżej wykorzystujemy trzy małe zbiory danych (każdy punkt to jeden kraj/miasto):\n",
    "\n",
    "1. liczba lekarzy a śmiertelność noworodków,\n",
    "2. wydatki na służbę zdrowia a oczekiwana długość życia,\n",
    "3. promieniowanie słoneczne a zachorowania na nowotwory.\n",
    "\n",
    "Najpierw wczytajmy dane.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cedc09c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "DATA_DIR = Path(\"data\")\n",
    "if not DATA_DIR.exists():\n",
    "    DATA_DIR = Path(\"09_Korelacja_i_regresja_liniowa\") / \"data\"\n",
    "\n",
    "d1 = pd.read_csv(DATA_DIR / \"Fig9-1a.dat\", sep=r\"\\s+\")\n",
    "d2 = pd.read_csv(DATA_DIR / \"Fig9-1b.dat\", sep=r\"\\s+\")\n",
    "d3 = pd.read_csv(DATA_DIR / \"Fig9-1c.dat\", sep=r\"\\s+\")\n",
    "\n",
    "print(\"d1 (pierwsze wiersze):\")\n",
    "print(d1.head())\n",
    "print(\"\\n---\")\n",
    "print(\"d2 (pierwsze wiersze):\")\n",
    "print(d2.head())\n",
    "print(\"\\n---\")\n",
    "print(\"d3 (pierwsze wiersze):\")\n",
    "print(d3.head())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89e0de07",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 3, figsize=(15, 4), constrained_layout=True)\n",
    "\n",
    "axes[0].scatter(d1[\"Physicians\"], d1[\"InfMort\"], alpha=0.9)\n",
    "axes[0].set_title(\"Lekarze a śmiertelność noworodków\")\n",
    "axes[0].set_xlabel(\"Liczba lekarzy\")\n",
    "axes[0].set_ylabel(\"Śmiertelność noworodków\")\n",
    "\n",
    "axes[1].scatter(d2[\"Expend\"], d2[\"LifeExp\"], alpha=0.9)\n",
    "axes[1].set_title(\"Wydatki a długość życia\")\n",
    "axes[1].set_xlabel(\"Wydatki na zdrowie\")\n",
    "axes[1].set_ylabel(\"Oczekiwana długość życia\")\n",
    "\n",
    "axes[2].scatter(d3[\"Radiation\"], d3[\"Cancer\"], alpha=0.9)\n",
    "axes[2].set_title(\"Promieniowanie a zachorowania\")\n",
    "axes[2].set_xlabel(\"Promieniowanie\")\n",
    "axes[2].set_ylabel(\"Zachorowania (nowotwory)\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47b059b9",
   "metadata": {},
   "source": [
    "## Dlaczego zawsze należy wizualizować dane? (przykład „datasaurus”)\n",
    "\n",
    "„Datasaurus” to przykład, który bardzo dobrze pokazuje, dlaczego **nie wystarczy patrzeć na same statystyki opisowe**. Dane mogą mieć prawie takie same średnie, odchylenia standardowe i korelacje, a jednocześnie tworzyć zupełnie inne wzory na wykresie.\n",
    "\n",
    "Krótka historia:\n",
    "\n",
    "- pierwotnego dinozaura narysował Alberto Cairo w 2016 roku jako ilustrację hasła: nie ufaj samym statystykom podsumowującym, zawsze wizualizuj dane;\n",
    "- Justin Matejka i George Fitzmaurice z Autodesk Research wykorzystali później tę ideę w pracy *Same Stats, Different Graphs* (CHI 2017);\n",
    "- za pomocą symulowanego wyżarzania stworzyli zestaw kształtów o prawie identycznych statystykach opisowych, znany jako **Datasaurus Dozen**;\n",
    "\n",
    "W tej sekcji pokazujemy trzy rzeczy: oryginalnego dinozaura, pełny zestaw Datasaurus Dozen, a potem lokalny wariant `data beaver`.\n",
    "\n",
    "Źródła: OpenIntro `datasaurus`, Autodesk Research: *Same Stats, Different Graphs*, pakiet R `datasauRus`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb3c06a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "datasaurus = pd.read_csv(DATA_DIR / \"datasaurus.csv\")\n",
    "dino = datasaurus[datasaurus[\"dataset\"] == \"dino\"]\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6, 5))\n",
    "ax.scatter(dino[\"x\"], dino[\"y\"], s=24, alpha=0.85, color=\"#2a9d8f\")\n",
    "ax.set_title(\"Datasaurus: najpierw zobaczmy wykres\")\n",
    "ax.set_xlabel(\"x\")\n",
    "ax.set_ylabel(\"y\")\n",
    "ax.set_xlim(0, 100)\n",
    "ax.set_ylim(0, 100)\n",
    "ax.set_aspect(\"equal\", adjustable=\"box\")\n",
    "plt.show()\n",
    "\n",
    "pd.DataFrame(datasaurus[\"dataset\"].value_counts().sort_index()).rename(columns={\"count\": \"n\"})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17904fd7",
   "metadata": {},
   "source": [
    "### Pełny zestaw: Datasaurus Dozen\n",
    "\n",
    "Na końcu zobaczmy pełny zestaw z OpenIntro: 13 kształtów, w tym oryginalny `dino`, które mają bardzo podobne statystyki opisowe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e13becc4",
   "metadata": {},
   "outputs": [],
   "source": [
    "summary = (\n",
    "    datasaurus.groupby(\"dataset\")\n",
    "    .agg(\n",
    "        mean_x=(\"x\", \"mean\"),\n",
    "        mean_y=(\"y\", \"mean\"),\n",
    "        sd_x=(\"x\", \"std\"),\n",
    "        sd_y=(\"y\", \"std\"),\n",
    "        cov_xy=(\"x\", lambda s: np.nan),\n",
    "    )\n",
    ")\n",
    "\n",
    "# Kowariancję i korelację policzymy jawnie, żeby było czytelnie, bez \"magii\".\n",
    "rows = []\n",
    "for dataset_name, group in datasaurus.groupby(\"dataset\"):\n",
    "    x = group[\"x\"]\n",
    "    y = group[\"y\"]\n",
    "    cov_xy = np.cov(x, y, ddof=1)[0, 1]\n",
    "    r_xy = np.corrcoef(x, y)[0, 1]\n",
    "    rows.append({\"dataset\": dataset_name, \"cov_xy\": cov_xy, \"r_xy\": r_xy})\n",
    "\n",
    "extra = pd.DataFrame(rows).set_index(\"dataset\")\n",
    "summary = summary.drop(columns=[\"cov_xy\"]).join(extra)\n",
    "\n",
    "pd.options.display.float_format = \"{:.3f}\".format\n",
    "summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e686677",
   "metadata": {},
   "outputs": [],
   "source": [
    "datasets = sorted(datasaurus[\"dataset\"].unique())\n",
    "\n",
    "n_cols = 4\n",
    "n_rows = int(np.ceil(len(datasets) / n_cols))\n",
    "\n",
    "fig, axes = plt.subplots(n_rows, n_cols, figsize=(14, 3.5 * n_rows), constrained_layout=True)\n",
    "axes = np.array(axes).reshape(n_rows, n_cols)\n",
    "\n",
    "for i, dataset_name in enumerate(datasets):\n",
    "    row = i // n_cols\n",
    "    col = i % n_cols\n",
    "    ax = axes[row, col]\n",
    "\n",
    "    group = datasaurus[datasaurus[\"dataset\"] == dataset_name]\n",
    "    ax.scatter(group[\"x\"], group[\"y\"], s=12, alpha=0.85)\n",
    "    ax.set_title(dataset_name)\n",
    "    ax.set_xlim(0, 100)\n",
    "    ax.set_ylim(0, 100)\n",
    "    ax.set_aspect(\"equal\", adjustable=\"box\")\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "\n",
    "# Usuń puste osie, jeśli liczba paneli nie wypełnia całej siatki.\n",
    "for j in range(len(datasets), n_rows * n_cols):\n",
    "    row = j // n_cols\n",
    "    col = j % n_cols\n",
    "    fig.delaxes(axes[row, col])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f94342ec",
   "metadata": {},
   "source": [
    "### `databeaver`\n",
    "\n",
    "Ten zestaw jest lokalną wariacją na ten sam temat: kilka bardzo różnych kształtów ma prawie te same statystyki opisowe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ef0da84",
   "metadata": {},
   "outputs": [],
   "source": [
    "beaver = pd.read_csv(DATA_DIR / \"beaver_full.csv\", index_col=0)\n",
    "pd.DataFrame(beaver[\"dataset\"].value_counts()).rename(columns={\"count\": \"n\"})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e94442ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "beaver_summary = (\n",
    "    beaver.groupby(\"dataset\")\n",
    "    .agg(\n",
    "        mean_x=(\"x\", \"mean\"),\n",
    "        mean_y=(\"y\", \"mean\"),\n",
    "        sd_x=(\"x\", \"std\"),\n",
    "        sd_y=(\"y\", \"std\"),\n",
    "        cov_xy=(\"x\", lambda s: np.nan),\n",
    "    )\n",
    ")\n",
    "\n",
    "# Kowariancję i korelację liczymy jawnie, żeby było widać dokładnie, co porównujemy.\n",
    "rows = []\n",
    "for dataset_name, group in beaver.groupby(\"dataset\"):\n",
    "    x = group[\"x\"]\n",
    "    y = group[\"y\"]\n",
    "    cov_xy = np.cov(x, y, ddof=1)[0, 1]\n",
    "    r_xy = np.corrcoef(x, y)[0, 1]\n",
    "    rows.append({\"dataset\": dataset_name, \"cov_xy\": cov_xy, \"r_xy\": r_xy})\n",
    "\n",
    "extra = pd.DataFrame(rows).set_index(\"dataset\")\n",
    "beaver_summary = beaver_summary.drop(columns=[\"cov_xy\"]).join(extra)\n",
    "\n",
    "pd.options.display.float_format = \"{:.3f}\".format\n",
    "beaver_summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d979ae7",
   "metadata": {},
   "outputs": [],
   "source": [
    "beaver_datasets = sorted(beaver[\"dataset\"].unique())\n",
    "\n",
    "n_cols = 2\n",
    "n_rows = int(np.ceil(len(beaver_datasets) / n_cols))\n",
    "\n",
    "fig, axes = plt.subplots(n_rows, n_cols, figsize=(12, 5 * n_rows), constrained_layout=True)\n",
    "axes = np.array(axes).reshape(n_rows, n_cols)\n",
    "\n",
    "for i, dataset_name in enumerate(beaver_datasets):\n",
    "    row = i // n_cols\n",
    "    col = i % n_cols\n",
    "    ax = axes[row, col]\n",
    "\n",
    "    group = beaver[beaver[\"dataset\"] == dataset_name]\n",
    "    ax.scatter(group[\"x\"], group[\"y\"], s=14, alpha=0.9)\n",
    "    ax.set_title(dataset_name)\n",
    "    ax.set_xlabel(\"x\")\n",
    "    ax.set_ylabel(\"y\")\n",
    "\n",
    "# Usuń puste osie, jeśli liczba paneli jest nieparzysta.\n",
    "for j in range(len(beaver_datasets), n_rows * n_cols):\n",
    "    row = j // n_cols\n",
    "    col = j % n_cols\n",
    "    fig.delaxes(axes[row, col])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba08c1e4",
   "metadata": {},
   "source": [
    "# Regresja liniowa\n",
    "\n",
    "Korelacja mówi, że dwie zmienne współzmieniają się liniowo.\n",
    "Regresja idzie krok dalej: zapisuje tę zależność jako równanie, którego możemy używać do opisu i przewidywania.\n",
    "\n",
    "Najprostszy model regresji liniowej ma postać:\n",
    "\n",
    "$$\n",
    "Y = a + bX + e\n",
    "$$\n",
    "\n",
    "gdzie:\n",
    "\n",
    "- $Y$ — zmienna wyjaśniana,\n",
    "- $X$ — predyktor,\n",
    "- $a$ — wyraz wolny,\n",
    "- $b$ — nachylenie prostej,\n",
    "- $e$ — reszta, czyli to, czego model nie przewidział.\n",
    "\n",
    "W Pythonie będziemy tu używać formułowego API `statsmodels`, np. `smf.ols(\"Y ~ X\", data=df).fit()`.\n",
    "Taka składnia jest bliska zapisowi statystycznemu: po lewej stronie mamy zmienną wyjaśnianą, po prawej predyktory.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4b82d68",
   "metadata": {},
   "source": [
    "## Dopasowanie modelu (prosty przykład)\n",
    "\n",
    "Wróćmy do danych z pierwszego wykresu: śmiertelność niemowląt i liczba lekarzy na 100 000 mieszkańców.\n",
    "\n",
    "Zanim poprosimy `statsmodels` o dopasowanie modelu, zatrzymajmy się na chwilę przy tym, co właściwie znaczy „dopasować linię”.\n",
    "Model liniowy przewiduje wartość $Y$ według równania:\n",
    "\n",
    "$$\\hat{Y} = a + bX$$\n",
    "\n",
    "Dla każdej obserwacji możemy policzyć **resztę**:\n",
    "\n",
    "$$e_i = Y_i - \\hat{Y}_i$$\n",
    "\n",
    "czyli pionową odległość między punktem a linią modelu.\n",
    "Metoda najmniejszych kwadratów wybiera taki wyraz wolny $a$ i takie nachylenie $b$, aby suma kwadratów reszt była jak najmniejsza:\n",
    "\n",
    "$$SSE = \\sum e_i^2$$\n",
    "\n",
    "Poniżej można ręcznie przesuwać linię dla tych samych danych i zobaczyć, jak zmieniają się reszty oraz suma kwadratów reszt. W widgecie $a$ oznacza poziom linii dla przeciętnej wartości $X$, dzięki czemu jego skala odpowiada osi Y na wykresie.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29614b90",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_demo = d1[\"Physicians\"]\n",
    "y_demo = d1[\"InfMort\"]\n",
    "x_center_demo = x_demo.mean()\n",
    "\n",
    "ols_model_demo = smf.ols(\"InfMort ~ Physicians\", data=d1).fit()\n",
    "ols_slope_demo = ols_model_demo.params[\"Physicians\"]\n",
    "ols_intercept_centered_demo = ols_model_demo.predict(\n",
    "    pd.DataFrame({\"Physicians\": [x_center_demo]})\n",
    ").iloc[0]\n",
    "\n",
    "\n",
    "def show_manual_fit(intercept, slope, show_ols):\n",
    "    y_hat = intercept + slope * (x_demo - x_center_demo)\n",
    "    residuals = y_demo - y_hat\n",
    "    sse = np.sum(residuals**2)\n",
    "    rmse = np.sqrt(np.mean(residuals**2))\n",
    "\n",
    "    x_line = np.linspace(10, 21, 200)\n",
    "    y_line = intercept + slope * (x_line - x_center_demo)\n",
    "\n",
    "    fig, ax = plt.subplots(figsize=(8, 5))\n",
    "    ax.scatter(x_demo, y_demo, s=70, alpha=0.9, label=\"obserwacje\")\n",
    "    ax.plot(x_line, y_line, color=\"black\", linewidth=2.5, label=\"linia z suwaków\")\n",
    "    ax.vlines(\n",
    "        x_demo,\n",
    "        ymin=np.minimum(y_demo, y_hat),\n",
    "        ymax=np.maximum(y_demo, y_hat),\n",
    "        color=\"tomato\",\n",
    "        alpha=0.65,\n",
    "        linewidth=2,\n",
    "        label=\"reszty\",\n",
    "    )\n",
    "\n",
    "    if show_ols:\n",
    "        y_ols = ols_intercept_centered_demo + ols_slope_demo * (x_line - x_center_demo)\n",
    "        ax.plot(\n",
    "            x_line,\n",
    "            y_ols,\n",
    "            color=\"steelblue\",\n",
    "            linestyle=\"--\",\n",
    "            linewidth=2.5,\n",
    "            label=\"minimum SSE (OLS)\",\n",
    "        )\n",
    "\n",
    "    ax.axvline(x_center_demo, color=\"gray\", linestyle=\":\", linewidth=1)\n",
    "    ax.set_xlim(10, 21)\n",
    "    ax.set_ylim(-8, 10)\n",
    "    ax.set_xlabel(\"Liczba lekarzy na 100 000 mieszkańców\")\n",
    "    ax.set_ylabel(\"Śmiertelność niemowląt\")\n",
    "    ax.set_title(f\"SSE = {sse:.2f}; RMSE = {rmse:.2f}\")\n",
    "    ax.legend(loc=\"upper left\")\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "intercept_slider = widgets.FloatSlider(\n",
    "    value=0.0,\n",
    "    min=-8.0,\n",
    "    max=10.0,\n",
    "    step=0.1,\n",
    "    description=\"a\",\n",
    "    readout_format=\".1f\",\n",
    "    continuous_update=False,\n",
    ")\n",
    "slope_slider = widgets.FloatSlider(\n",
    "    value=0.0,\n",
    "    min=-2.0,\n",
    "    max=2.0,\n",
    "    step=0.02,\n",
    "    description=\"b\",\n",
    "    readout_format=\".2f\",\n",
    "    continuous_update=False,\n",
    ")\n",
    "show_ols_checkbox = widgets.Checkbox(value=False, description=\"pokaż linię OLS\")\n",
    "\n",
    "manual_fit_output = widgets.interactive_output(\n",
    "    show_manual_fit,\n",
    "    {\n",
    "        \"intercept\": intercept_slider,\n",
    "        \"slope\": slope_slider,\n",
    "        \"show_ols\": show_ols_checkbox,\n",
    "    },\n",
    ")\n",
    "\n",
    "display(\n",
    "    widgets.VBox(\n",
    "        [\n",
    "            widgets.HBox([intercept_slider, slope_slider]),\n",
    "            show_ols_checkbox,\n",
    "            manual_fit_output,\n",
    "        ]\n",
    "    )\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cc1717c",
   "metadata": {},
   "source": [
    "To, co przed chwilą robiliśmy ręcznie suwakiem, `statsmodels` robi automatycznie: znajduje parametry minimalizujące sumę kwadratów reszt.\n",
    "\n",
    "W składni formułowej `Y ~ X` wyraz wolny jest dodawany automatycznie.\n",
    "Gdybyśmy pracowali bezpośrednio na macierzy projektującej, musielibyśmy dodać kolumnę stałej ręcznie.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6c6f3b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = d1[\"Physicians\"]\n",
    "y = d1[\"InfMort\"]\n",
    "\n",
    "model_d1 = smf.ols(\"InfMort ~ Physicians\", data=d1).fit()\n",
    "\n",
    "print(\"Parametry modelu:\")\n",
    "model_d1.params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2f94e46",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_d1.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae6591d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_grid = np.linspace(x.min(), x.max(), 200)\n",
    "grid = pd.DataFrame({\"Physicians\": x_grid})\n",
    "y_grid = model_d1.predict(grid)\n",
    "\n",
    "plt.scatter(x, y, alpha=0.8)\n",
    "plt.plot(x_grid, y_grid, color=\"black\", linewidth=2)\n",
    "plt.xlabel(\"Liczba lekarzy na 100 000 mieszkańców\")\n",
    "plt.ylabel(\"Śmiertelność niemowląt\")\n",
    "plt.title(\"Regresja liniowa: InfMort ~ Physicians\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5985448",
   "metadata": {},
   "source": [
    "## Linia regresji na trzech przykładach\n",
    "\n",
    "Na tym etapie możemy dorysować linie regresji do wszystkich trzech wykresów punktowych.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "314338a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 3, figsize=(15, 4), constrained_layout=True)\n",
    "\n",
    "# 1) InfMort ~ Physicians\n",
    "x1 = d1[\"Physicians\"]\n",
    "y1 = d1[\"InfMort\"]\n",
    "model1 = smf.ols(\"InfMort ~ Physicians\", data=d1).fit()\n",
    "x1_grid = np.linspace(x1.min(), x1.max(), 200)\n",
    "y1_grid = model1.predict(pd.DataFrame({\"Physicians\": x1_grid}))\n",
    "axes[0].scatter(x1, y1, alpha=0.8)\n",
    "axes[0].plot(x1_grid, y1_grid, color=\"black\")\n",
    "axes[0].set_title(\"InfMort ~ Physicians\")\n",
    "axes[0].set_xlabel(\"Physicians\")\n",
    "axes[0].set_ylabel(\"InfMort\")\n",
    "\n",
    "# 2) LifeExp ~ Expend\n",
    "x2 = d2[\"Expend\"]\n",
    "y2 = d2[\"LifeExp\"]\n",
    "model2 = smf.ols(\"LifeExp ~ Expend\", data=d2).fit()\n",
    "x2_grid = np.linspace(x2.min(), x2.max(), 200)\n",
    "y2_grid = model2.predict(pd.DataFrame({\"Expend\": x2_grid}))\n",
    "axes[1].scatter(x2, y2, alpha=0.8)\n",
    "axes[1].plot(x2_grid, y2_grid, color=\"black\")\n",
    "axes[1].set_title(\"LifeExp ~ Expend\")\n",
    "axes[1].set_xlabel(\"Expend\")\n",
    "axes[1].set_ylabel(\"LifeExp\")\n",
    "\n",
    "# 3) Cancer ~ Radiation\n",
    "x3 = d3[\"Radiation\"]\n",
    "y3 = d3[\"Cancer\"]\n",
    "model3 = smf.ols(\"Cancer ~ Radiation\", data=d3).fit()\n",
    "x3_grid = np.linspace(x3.min(), x3.max(), 200)\n",
    "y3_grid = model3.predict(pd.DataFrame({\"Radiation\": x3_grid}))\n",
    "axes[2].scatter(x3, y3, alpha=0.8)\n",
    "axes[2].plot(x3_grid, y3_grid, color=\"black\")\n",
    "axes[2].set_title(\"Cancer ~ Radiation\")\n",
    "axes[2].set_xlabel(\"Radiation\")\n",
    "axes[2].set_ylabel(\"Cancer\")\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25a742fc",
   "metadata": {},
   "source": [
    "# Kowariancja i współczynnik korelacji *r* Pearsona\n",
    "\n",
    "## Kowariancja\n",
    "\n",
    "Wzór na współczynnik kowariancji między dwiema zmiennymi wygląda następująco:\n",
    "\n",
    "$$\n",
    "Cov(X,Y) = E\\big[(X-E(X))\\,(Y-E(Y))\\big]\n",
    "$$\n",
    "\n",
    "Intuicyjnie:\n",
    "\n",
    "- jeśli dużym wartościom $X$ (względem średniej) towarzyszą duże wartości $Y$ → kowariancja dodatnia,\n",
    "- jeśli dużym $X$ towarzyszą małe $Y$ → kowariancja ujemna,\n",
    "- jeśli raz tak, raz tak → znoszenie się składników i kowariancja blisko 0.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a36f1413",
   "metadata": {},
   "source": [
    "## Korelacja\n",
    "\n",
    "Problem z kowariancją polega na tym, że jej wielkość zależy od jednostek i skali zmiennych.\n",
    "\n",
    "Korelacja jest „wystandaryzowaną” kowariancją. Dla populacji:\n",
    "\n",
    "$$\n",
    "\\rho = \\frac{Cov(X,Y)}{\\sqrt{Var(X)Var(Y)}}\n",
    "$$\n",
    "\n",
    "Dla próby (korelacja Pearsona):\n",
    "\n",
    "$$\n",
    "r = \\frac{cov_{XY}}{s_X s_Y}\n",
    "$$\n",
    "\n",
    "Wartość $r$ jest w przedziale $[-1, 1]$ i mierzy **siłę liniowej zależności**.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38626f54",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sprawdźmy \"na piechotę\" kowariancję i korelację na danych d1.\n",
    "\n",
    "x = d1[\"Physicians\"]\n",
    "y = d1[\"InfMort\"]\n",
    "\n",
    "n = len(x)\n",
    "x_mean = x.mean()\n",
    "y_mean = y.mean()\n",
    "\n",
    "x_dev = x - x_mean\n",
    "y_dev = y - y_mean\n",
    "\n",
    "cov_xy_hand = (x_dev * y_dev).sum() / (n - 1)\n",
    "\n",
    "print(\"Kowariancja (ręcznie):\", cov_xy_hand)\n",
    "print(\"Kowariancja (np.cov): \", np.cov(x, y, ddof=1)[0, 1])\n",
    "print(\"Kowariancja (pandas):\", x.cov(pd.Series(y)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "287d2eaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Korelacja \"ręcznie\": r = cov / (sd_x * sd_y)\n",
    "\n",
    "sd_x = x.std(ddof=1)\n",
    "sd_y = y.std(ddof=1)\n",
    "\n",
    "r_hand = cov_xy_hand / (sd_x * sd_y)\n",
    "\n",
    "r_scipy, p_scipy = stats.pearsonr(x, y)\n",
    "\n",
    "print(\"Korelacja (ręcznie):\", r_hand)\n",
    "print(\"Korelacja (np.corrcoef):\", np.corrcoef(x, y)[0, 1])\n",
    "print(\"Korelacja (scipy.stats.pearsonr):\", r_scipy)\n",
    "print(\"p-wartość (scipy.stats.pearsonr):\", p_scipy)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9522ac10",
   "metadata": {},
   "source": [
    "## Interpretacja wartości *r*\n",
    "\n",
    "Interpretując współczynnik korelacji Pearsona pamiętajmy, że „duże/małe” zależy od dziedziny.\n",
    "\n",
    "Przykładowe (bardzo orientacyjne) zasady kciuka w psychologii:\n",
    "\n",
    "- $r = 0$ – brak liniowej współzależności,\n",
    "- $0 < |r| < 0.3$ – słaby związek,\n",
    "- $0.3 \\le |r| < 0.5$ – średni związek,\n",
    "- $0.5 \\le |r| < 0.7$ – znaczny związek,\n",
    "- $0.7 \\le |r| < 0.9$ – wysoki związek,\n",
    "- $|r| \\ge 0.9$ – bardzo wysoki związek.\n",
    "\n",
    "W praktyce warto zawsze podawać też **niepewność** (np. przedział ufności) i **kontekst**.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ae3396d",
   "metadata": {},
   "source": [
    "## Macierze korelacyjne – wizualizacje\n",
    "\n",
    "Jeśli mamy więcej niż dwie zmienne, często chcemy policzyć korelacje dla wszystkich par zmiennych naraz.\n",
    "\n",
    "Zrobimy to na dwóch przykładach:\n",
    "\n",
    "1. dane sztuczne (losowe) – żeby zobaczyć „jak wygląda” macierz korelacji, gdy zmienne są niezależne,\n",
    "2. dane rzeczywiste z wieloma zmiennymi\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88382a73",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1) Dane losowe: 35 obserwacji, 10 zmiennych\n",
    "\n",
    "n_obs = 35\n",
    "n_vars = 10\n",
    "\n",
    "m = rng.normal(size=(n_obs, n_vars))\n",
    "columns = [f\"X{i+1}\" for i in range(n_vars)]\n",
    "\n",
    "m_df = pd.DataFrame(m, columns=columns)\n",
    "\n",
    "corr_m = m_df.corr()\n",
    "print(corr_m.round(2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac54ab34",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 6))\n",
    "ax = sns.heatmap(\n",
    "    corr_m,\n",
    "    vmin=-1,\n",
    "    vmax=1,\n",
    "    center=0,\n",
    "    cmap=\"vlag\",\n",
    "    annot=True,\n",
    "    fmt=\".2f\",\n",
    "    square=True,\n",
    ")\n",
    "ax.set_title(\"Macierz korelacji (dane losowe)\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a577cca",
   "metadata": {},
   "source": [
    "### Dane rzeczywiste: przykład z Afryki subsaharyjskiej\n",
    "\n",
    "Teraz wykorzystamy realny zbiór danych z przykładu Howella dotyczący krajów Afryki subsaharyjskiej.\n",
    "Interesuje nas nie pojedyncza korelacja, ale **wzorzec współzmienności** kilku wskaźników demograficznych i społecznych.\n",
    "\n",
    "W tabeli mamy m.in.:\n",
    "\n",
    "- `infant_mortality` - śmiertelność niemowląt,\n",
    "- `income` - dochód,\n",
    "- `young_mothers`, `older_mothers`, `too_soon` - wskaźniki opisujące strukturę urodzeń,\n",
    "- `contraception` - użycie antykoncepcji,\n",
    "- `unmet_need` - niezaspokojona potrzeba antykoncepcji.\n",
    "\n",
    "Dwie obserwacje oznaczone `in_analysis == 0` zostawiamy w pliku dla przejrzystości źródła, ale w analizie korzystamy z głównego zestawu krajów."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db394fdc",
   "metadata": {},
   "outputs": [],
   "source": [
    "africa_raw = pd.read_csv(DATA_DIR / \"africa_subsaharan.csv\")\n",
    "africa = africa_raw[africa_raw[\"in_analysis\"] == 1]\n",
    "\n",
    "africa.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f12c144",
   "metadata": {},
   "outputs": [],
   "source": [
    "africa_variables = [\n",
    "    \"infant_mortality\",\n",
    "    \"income\",\n",
    "    \"young_mothers\",\n",
    "    \"older_mothers\",\n",
    "    \"too_soon\",\n",
    "    \"contraception\",\n",
    "    \"unmet_need\",\n",
    "]\n",
    "\n",
    "print(\"Liczba krajów w analizie:\", len(africa))\n",
    "print(\"Zmienne w macierzy korelacji:\")\n",
    "print(africa_variables)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92913ca5",
   "metadata": {},
   "outputs": [],
   "source": [
    "corr_africa = africa[africa_variables].corr()\n",
    "corr_africa.round(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "622843ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(9, 7))\n",
    "ax = sns.heatmap(\n",
    "    corr_africa,\n",
    "    vmin=-1,\n",
    "    vmax=1,\n",
    "    center=0,\n",
    "    cmap=\"vlag\",\n",
    "    annot=True,\n",
    "    fmt=\".2f\",\n",
    "    square=True,\n",
    ")\n",
    "ax.set_title(\"Macierz korelacji: kraje Afryki subsaharyjskiej\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "657bf9ce",
   "metadata": {},
   "source": [
    "# Widget: zgadnij korelację\n",
    "\n",
    "Dobra intuicja dotycząca tego, „jaki rozrzut punktów odpowiada jakiej korelacji”, bardzo pomaga.\n",
    "Zrobimy krótką grę: notebook losuje chmurę punktów, a my próbujemy oszacować **korelację w tej konkretnej próbie**.\n",
    "\n",
    "Nie zgadujemy ukrytego parametru populacyjnego. Po kliknięciu przycisku porównujemy odpowiedź z korelacją policzoną dla punktów widocznych na wykresie.\n",
    "\n",
    "Ważne: skala osi jest stała, więc zmienia się wyłącznie układ punktów, a nie „zoom” wykresu.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15918555",
   "metadata": {},
   "outputs": [],
   "source": [
    "correlation_game_state = {\n",
    "    \"rng\": np.random.default_rng(SEED),\n",
    "    \"xy\": None,\n",
    "    \"sample_r\": None,\n",
    "    \"answer_visible\": False,\n",
    "    \"last_guess\": 0.0,\n",
    "}\n",
    "\n",
    "\n",
    "def generate_correlation_data(n):\n",
    "    rng_local = correlation_game_state[\"rng\"]\n",
    "    rho = rng_local.uniform(-1.0, 1.0)\n",
    "    covariance = np.array([[1.0, rho], [rho, 1.0]])\n",
    "    xy = rng_local.multivariate_normal(mean=[0.0, 0.0], cov=covariance, size=n)\n",
    "\n",
    "    correlation_game_state[\"xy\"] = xy\n",
    "    correlation_game_state[\"sample_r\"] = np.corrcoef(xy[:, 0], xy[:, 1])[0, 1]\n",
    "    correlation_game_state[\"answer_visible\"] = False\n",
    "    correlation_game_state[\"last_guess\"] = 0.0\n",
    "\n",
    "\n",
    "def draw_correlation_game():\n",
    "    xy = correlation_game_state[\"xy\"]\n",
    "    sample_r = correlation_game_state[\"sample_r\"]\n",
    "\n",
    "    if correlation_game_state[\"answer_visible\"]:\n",
    "        shown_guess = correlation_game_state[\"last_guess\"]\n",
    "    else:\n",
    "        shown_guess = guess_slider.value\n",
    "\n",
    "    x = xy[:, 0]\n",
    "    y = xy[:, 1]\n",
    "\n",
    "    fig, ax = plt.subplots(figsize=(6, 6))\n",
    "    ax.scatter(x, y, alpha=0.85)\n",
    "    ax.axhline(0, color=\"gray\", linewidth=0.8, alpha=0.5)\n",
    "    ax.axvline(0, color=\"gray\", linewidth=0.8, alpha=0.5)\n",
    "    ax.set_xlim(-3.5, 3.5)\n",
    "    ax.set_ylim(-3.5, 3.5)\n",
    "    ax.set_aspect(\"equal\", adjustable=\"box\")\n",
    "    ax.set_xlabel(\"X\")\n",
    "    ax.set_ylabel(\"Y\")\n",
    "\n",
    "    if correlation_game_state[\"answer_visible\"]:\n",
    "        delta = sample_r - shown_guess\n",
    "        result_text = \"\\n\".join(\n",
    "            [\n",
    "                f\"r w próbie = {sample_r:.2f}\",\n",
    "                f\"odpowiedź = {shown_guess:.2f}\",\n",
    "                f\"delta = {delta:.2f}\",\n",
    "                f\"|delta| = {abs(delta):.2f}\",\n",
    "            ]\n",
    "        )\n",
    "        ax.set_title(\"Wynik\")\n",
    "        ax.text(\n",
    "            0.04,\n",
    "            0.96,\n",
    "            result_text,\n",
    "            transform=ax.transAxes,\n",
    "            va=\"top\",\n",
    "            ha=\"left\",\n",
    "            fontsize=11,\n",
    "            bbox={\n",
    "                \"boxstyle\": \"round,pad=0.45\",\n",
    "                \"facecolor\": \"white\",\n",
    "                \"edgecolor\": \"black\",\n",
    "                \"alpha\": 0.88,\n",
    "            },\n",
    "        )\n",
    "    else:\n",
    "        ax.set_title(f\"Moja odpowiedź: r ≈ {shown_guess:.2f}\")\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def refresh_correlation_output():\n",
    "    with correlation_output:\n",
    "        clear_output(wait=True)\n",
    "        draw_correlation_game()\n",
    "\n",
    "\n",
    "def reset_guess_slider():\n",
    "    if guess_slider.value != 0.0:\n",
    "        guess_slider.value = 0.0\n",
    "    refresh_correlation_output()\n",
    "\n",
    "\n",
    "def handle_correlation_action(_):\n",
    "    if correlation_game_state[\"answer_visible\"]:\n",
    "        generate_correlation_data(n_slider.value)\n",
    "        action_button.description = \"Sprawdź\"\n",
    "        reset_guess_slider()\n",
    "    else:\n",
    "        correlation_game_state[\"last_guess\"] = guess_slider.value\n",
    "        correlation_game_state[\"answer_visible\"] = True\n",
    "        action_button.description = \"Nowy wykres\"\n",
    "        reset_guess_slider()\n",
    "\n",
    "\n",
    "def update_guess(_):\n",
    "    if not correlation_game_state[\"answer_visible\"]:\n",
    "        refresh_correlation_output()\n",
    "\n",
    "\n",
    "def update_sample_size(_):\n",
    "    generate_correlation_data(n_slider.value)\n",
    "    action_button.description = \"Sprawdź\"\n",
    "    reset_guess_slider()\n",
    "\n",
    "\n",
    "n_slider = widgets.IntSlider(\n",
    "    value=120,\n",
    "    min=30,\n",
    "    max=300,\n",
    "    step=10,\n",
    "    description=\"n\",\n",
    "    continuous_update=False,\n",
    ")\n",
    "guess_slider = widgets.FloatSlider(\n",
    "    value=0.0,\n",
    "    min=-0.95,\n",
    "    max=0.95,\n",
    "    step=0.05,\n",
    "    description=\"zgaduję r\",\n",
    "    readout_format=\".2f\",\n",
    "    continuous_update=False,\n",
    ")\n",
    "action_button = widgets.Button(description=\"Sprawdź\")\n",
    "correlation_output = widgets.Output()\n",
    "\n",
    "n_slider.observe(update_sample_size, names=\"value\")\n",
    "guess_slider.observe(update_guess, names=\"value\")\n",
    "action_button.on_click(handle_correlation_action)\n",
    "\n",
    "generate_correlation_data(n_slider.value)\n",
    "refresh_correlation_output()\n",
    "display(\n",
    "    widgets.VBox(\n",
    "        [\n",
    "            widgets.HBox([n_slider, guess_slider]),\n",
    "            action_button,\n",
    "            correlation_output,\n",
    "        ]\n",
    "    )\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45d1dd59",
   "metadata": {},
   "source": [
    "# Linia regresji: stres a symptomy psychologiczne\n",
    "\n",
    "Howell omawia ten przykład jako klasyczną ilustrację związku między stresem a zdrowiem psychicznym.\n",
    "Punktem wyjścia jest badanie Wagnera, Compasa i Howella (1988) dotyczące studentów pierwszego roku.\n",
    "\n",
    "Badacze mierzyli stres jako wynik oparty na niedawnych wydarzeniach życiowych: liczyła się nie tylko liczba zdarzeń, ale też ich częstość, ważność i subiektywna ocena wpływu.\n",
    "Symptomy psychologiczne mierzono kwestionariuszem Hopkins Symptom Checklist, obejmującym 57 objawów.\n",
    "Surowa liczba symptomów była wyraźnie skośna, dlatego Howell pracuje z logarytmem naturalnym tej zmiennej (`lnSymptoms`).\n",
    "\n",
    "Nasze pytanie brzmi:\n",
    "\n",
    "> Czy osoby z wyższym poziomem stresu raportują przeciętnie więcej symptomów psychologicznych?\n",
    "\n",
    "To nadal analiza obserwacyjna, więc nie będziemy udawać, że sama linia regresji dowodzi przyczynowości.\n",
    "Używamy jej do opisu zależności i do nauki interpretacji parametrów modelu.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be5e6cc7",
   "metadata": {},
   "outputs": [],
   "source": [
    "stress_df = pd.read_csv(DATA_DIR / \"Tab9-2.dat\", sep=r\"\\s+\")\n",
    "stress_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d547d06",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_stress = stress_df[\"Stress\"]\n",
    "Y_ln = stress_df[\"lnSymptoms\"]\n",
    "\n",
    "plt.scatter(X_stress, Y_ln, alpha=0.9)\n",
    "plt.title(\"Stres a logarytm liczby symptomów\")\n",
    "plt.xlabel(\"Stress\")\n",
    "plt.ylabel(\"lnSymptoms\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0dd16dd5",
   "metadata": {},
   "source": [
    "## Dopasowanie regresji i odczyt parametrów\n",
    "\n",
    "Dopasujemy model:\n",
    "\n",
    "$$\\text{lnSymptoms} = a + b\\,\\text{Stress} + \\varepsilon$$\n",
    "\n",
    "Jeżeli nachylenie $b$ jest dodatnie, to wyższy poziom stresu wiąże się z wyższą przewidywaną wartością `lnSymptoms`.\n",
    "Jeżeli $b$ jest bliskie zeru, to sama zmienność stresu niewiele pomagałaby w przewidywaniu symptomów.\n",
    "\n",
    "Po dopasowaniu modelu będziemy czytać wynik bardzo konkretnie:\n",
    "\n",
    "- *intercept* ($a$): przewidywana wartość `lnSymptoms`, gdy `Stress = 0`,\n",
    "- *slope* ($b$): przeciętna zmiana przewidywanego `lnSymptoms` przy wzroście `Stress` o 1 punkt,\n",
    "- przedział ufności i test dla nachylenia: jak niepewne jest nasze oszacowanie związku.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c48fa535",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_stress = smf.ols(\"lnSymptoms ~ Stress\", data=stress_df).fit()\n",
    "model_stress.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa9f17e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_grid = np.linspace(X_stress.min(), X_stress.max(), 200)\n",
    "stress_grid = pd.DataFrame({\"Stress\": x_grid})\n",
    "y_grid = model_stress.predict(stress_grid)\n",
    "\n",
    "plt.scatter(X_stress, Y_ln, alpha=0.9)\n",
    "plt.plot(x_grid, y_grid, color=\"black\", linewidth=2)\n",
    "plt.xlabel(\"Stress\")\n",
    "plt.ylabel(\"ln(Symptoms)\")\n",
    "plt.title(\"Model regresji dla danych o stresie\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39604a3a",
   "metadata": {},
   "source": [
    "## Wyraz wolny (*intercept*)\n",
    "\n",
    "- **Definicja:** wartość $\\hat{Y}$ kiedy $X=0$.\n",
    "- Jego interpretacja zależy od tego, czy $X=0$ ma sens w danym kontekście.\n",
    "\n",
    "Często (dla sensowniejszej interpretacji) **centruje się predyktor** wokół średniej.\n",
    "Wtedy *intercept* odpowiada $\\hat{Y}$ dla „przeciętnego” $X$.\n",
    "\n",
    "Wycentrowanie nie zmienia *slope*.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35b5e835",
   "metadata": {},
   "outputs": [],
   "source": [
    "intercept = model_stress.params[\"Intercept\"]\n",
    "slope = model_stress.params[\"Stress\"]\n",
    "print(\"Intercept (X=0):\", intercept)\n",
    "print(\"Slope:\", slope)\n",
    "\n",
    "stress_df = stress_df.assign(\n",
    "    Stress_centered=stress_df[\"Stress\"] - stress_df[\"Stress\"].mean()\n",
    ")\n",
    "model_centered = smf.ols(\"lnSymptoms ~ Stress_centered\", data=stress_df).fit()\n",
    "\n",
    "intercept_c = model_centered.params[\"Intercept\"]\n",
    "slope_c = model_centered.params[\"Stress_centered\"]\n",
    "print()\n",
    "print(\"Po centrowaniu X:\")\n",
    "print(\"Intercept (X=średnia):\", intercept_c)\n",
    "print(\"Slope:\", slope_c)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c39d17c1",
   "metadata": {},
   "source": [
    "## Współczynnik kierunkowy (*slope*)\n",
    "\n",
    "- **Interpretacja:** o ile (średnio) zmieni się $\\hat{Y}$, gdy $X$ wzrośnie o 1 jednostkę.\n",
    "\n",
    "W prostej regresji liniowej z jednym predyktorem da się go policzyć „od podstaw”:\n",
    "\n",
    "$$\n",
    "\\hat{b} = \\frac{cov_{XY}}{var_X}\n",
    "$$\n",
    "\n",
    "a potem:\n",
    "\n",
    "$$\n",
    "\\hat{a} = \\bar{Y} - \\hat{b}\\,\\bar{X}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d91205a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parametry \"na piechotę\"\n",
    "\n",
    "x = X_stress\n",
    "y = Y_ln\n",
    "\n",
    "cov_xy = np.cov(x, y, ddof=1)[0, 1]\n",
    "var_x = np.var(x, ddof=1)\n",
    "\n",
    "b_hand = cov_xy / var_x\n",
    "a_hand = y.mean() - b_hand * x.mean()\n",
    "\n",
    "print(\"a (ręcznie):\", a_hand)\n",
    "print(\"b (ręcznie):\", b_hand)\n",
    "print(\"a (statsmodels):\", intercept)\n",
    "print(\"b (statsmodels):\", slope)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae30f0ca",
   "metadata": {},
   "source": [
    "## Test istotności współczynnika kierunkowego (*slope*) – krok po kroku\n",
    "\n",
    "**Pytanie.** Czy w populacji istnieje liniowy związek między stresem a (log) liczbą symptomów?\n",
    "\n",
    "Testujemy hipotezę o współczynniku kierunkowym:\n",
    "\n",
    "- $H_0: \\beta = 0$ (brak liniowego związku; prosta jest „pozioma”),\n",
    "- $H_A: \\beta \\neq 0$ (związek istnieje).\n",
    "\n",
    "Ustalamy poziom istotności, np. $\\alpha = 0.05$.\n",
    "\n",
    "Statystyka testowa ma postać:\n",
    "\n",
    "$$t = \\frac{\\hat{\\beta} - 0}{SE(\\hat{\\beta})}$$\n",
    "\n",
    "i przy $H_0$ ma rozkład t z $n-2$ stopniami swobody.\n",
    "\n",
    "Poniżej policzymy:\n",
    "\n",
    "- $t$,\n",
    "- region krytyczny,\n",
    "- p-wartość „na piechotę”,\n",
    "- przedział ufności dla $\\beta$,\n",
    "\n",
    "…a następnie porównamy z wynikiem `statsmodels`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0259c193",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 0.05\n",
    "n = len(X_stress)\n",
    "x_mean = X_stress.mean()\n",
    "residuals = model_stress.resid\n",
    "s_yx = np.sqrt(np.sum(residuals**2) / (n - 2))\n",
    "se_b = s_yx / np.sqrt(np.sum((X_stress - x_mean) ** 2))\n",
    "t_value = slope / se_b\n",
    "df = n - 2\n",
    "p_value = 2 * stats.t.sf(abs(t_value), df=df)\n",
    "crit = stats.t.ppf(1 - alpha / 2, df=df)\n",
    "ci_low = slope - crit * se_b\n",
    "ci_high = slope + crit * se_b\n",
    "\n",
    "print(\"b_hat =\", slope)\n",
    "print(\"SE(b_hat) =\", se_b)\n",
    "print(\"t =\", t_value)\n",
    "print(\"df =\", df)\n",
    "print(\"p =\", p_value)\n",
    "print(\"95% CI =\", (ci_low, ci_high))\n",
    "print()\n",
    "print(\"b_hat (statsmodels) =\", model_stress.params[\"Stress\"])\n",
    "print(\"SE(b_hat) (statsmodels) =\", model_stress.bse[\"Stress\"])\n",
    "print(\"t (statsmodels) =\", model_stress.tvalues[\"Stress\"])\n",
    "print(\"p (statsmodels) =\", model_stress.pvalues[\"Stress\"])\n",
    "print(\"CI 95% (statsmodels) =\", model_stress.conf_int(alpha=alpha).loc[\"Stress\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a18a7fe",
   "metadata": {},
   "source": [
    "## Standaryzowany współczynnik regresji\n",
    "\n",
    "Jeśli wystandaryzujemy obie zmienne (z-score), to *slope* też staje się „w skali odchyleń standardowych”.\n",
    "\n",
    "W prostej regresji z jednym predyktorem zachodzi ważna zależność:\n",
    "\n",
    "- standaryzowany *slope* = korelacja $r$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "353a1d9d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Standaryzowany slope: b_std = b * (sd_x / sd_y)\n",
    "\n",
    "sd_x = x.std(ddof=1)\n",
    "sd_y = y.std(ddof=1)\n",
    "\n",
    "b_std = slope * (sd_x / sd_y)\n",
    "r_xy = np.corrcoef(x, y)[0, 1]\n",
    "\n",
    "print(\"Standaryzowany slope:\", b_std)\n",
    "print(\"Korelacja r:\", r_xy)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64b1e6e0",
   "metadata": {},
   "source": [
    "# Testowanie hipotez dla współczynnika korelacji $r$ Pearsona\n",
    "\n",
    "Najczęściej sprawdzamy, czy dwie zmienne są skorelowane w populacji.\n",
    "\n",
    "Hipotezy (test dwustronny):\n",
    "\n",
    "$$\n",
    "H_A: \\rho \\neq 0\n",
    "$$\n",
    "\n",
    "$$\n",
    "H_0: \\rho = 0\n",
    "$$\n",
    "\n",
    "Statystyka testowa:\n",
    "\n",
    "$$\n",
    "T = \\frac{r}{\\sqrt{1-r^2}}\\sqrt{n-2}\n",
    "$$\n",
    "\n",
    "Przy prawdziwości $H_0$ ma rozkład t z $n-2$ stopniami swobody.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c16d386c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Krok-po-kroku: test dla korelacji między Stress i lnSymptoms\n",
    "\n",
    "alpha = 0.05\n",
    "\n",
    "x = X_stress\n",
    "y = Y_ln\n",
    "n = len(x)\n",
    "\n",
    "r, _ = stats.pearsonr(x, y)\n",
    "df = n - 2\n",
    "\n",
    "t_stat = (r / np.sqrt(1 - r**2)) * np.sqrt(df)\n",
    "\n",
    "t_crit = stats.t.ppf(1 - alpha / 2, df=df)\n",
    "\n",
    "p_value_hand = 2 * (1 - stats.t.cdf(abs(t_stat), df=df))\n",
    "\n",
    "print(f\"n = {n}, df = {df}\")\n",
    "print(\"r =\", r)\n",
    "print(\"t =\", t_stat)\n",
    "print(f\"t krytyczne (|t| >= {t_crit:.3f})\")\n",
    "print(\"p-wartość (na piechotę) =\", p_value_hand)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9cf24cbd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Weryfikacja wyniku funkcją biblioteczną\n",
    "\n",
    "r_scipy, p_scipy = stats.pearsonr(X_stress, Y_ln)\n",
    "\n",
    "print(\"r (scipy):\", r_scipy)\n",
    "print(\"p (scipy):\", p_scipy)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2cf7ef9",
   "metadata": {},
   "source": [
    "### Przedział ufności dla korelacji (Fisher z)\n",
    "\n",
    "Oprócz p-wartości warto podać niepewność estymacji efektu.\n",
    "\n",
    "Dla korelacji Pearsona często używa się przybliżenia opartego o transformację Fishera:\n",
    "\n",
    "- $z = \\operatorname{arctanh}(r)$,\n",
    "- $SE(z) = 1/\\sqrt{n-3}$,\n",
    "- przedział dla $z$ liczymy z rozkładu normalnego, a potem wracamy transformacją $\\tanh$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72eefebe",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 0.05\n",
    "\n",
    "r, _ = stats.pearsonr(X_stress, Y_ln)\n",
    "n = len(X_stress)\n",
    "\n",
    "z = np.arctanh(r)\n",
    "se_z = 1 / np.sqrt(n - 3)\n",
    "\n",
    "z_crit = stats.norm.ppf(1 - 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\"r = {r:.3f}\")\n",
    "print(f\"CI 95% dla \\u03c1 (Fisher z): [{r_low:.3f}, {r_high:.3f}]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0af8f07f",
   "metadata": {},
   "source": [
    "## Rozkład statystyki pod $H_0$ (symulacja)\n",
    "\n",
    "Możemy sprawdzić symulacyjnie, że statystyka t zbudowana z korelacji ma pod $H_0$ rozkład t z $n-2$ stopniami swobody.\n",
    "\n",
    "Losujemy niezależne $X$ i $Y$ (czyli prawdziwe $\\rho=0$), liczymy $r$ i przekształcamy do t.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f86ddf01",
   "metadata": {},
   "outputs": [],
   "source": [
    "replications = 10_000\n",
    "n = 35\n",
    "\n",
    "t_stats = np.empty(replications)\n",
    "for i in range(replications):\n",
    "    x = rng.normal(size=n)\n",
    "    y = rng.normal(size=n)\n",
    "\n",
    "    r = np.corrcoef(x, y)[0, 1]\n",
    "    t_stats[i] = (r / np.sqrt(1 - r**2)) * np.sqrt(n - 2)\n",
    "\n",
    "# Histogram + gęstość t\n",
    "x_grid = np.linspace(-5, 5, 400)\n",
    "\n",
    "plt.hist(t_stats, bins=50, density=True, alpha=0.7, label=\"symulacja\")\n",
    "plt.plot(x_grid, stats.t.pdf(x_grid, df=n - 2), linewidth=2.5, label=\"t(df=n-2)\")\n",
    "\n",
    "plt.title(\"Rozkład statystyki t dla korelacji pod $H_0$\")\n",
    "plt.xlabel(\"t\")\n",
    "plt.ylabel(\"gęstość\")\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea2bfa5e",
   "metadata": {},
   "source": [
    "## Wiele korelacji naraz i „fałszywe alarmy”\n",
    "\n",
    "Jeśli policzymy bardzo dużo korelacji i dla każdej wykonamy test przy $\\alpha=0.05$, to nawet gdy wszystkie prawdziwe korelacje w populacji wynoszą 0, część testów wyjdzie „istotna” przez przypadek.\n",
    "\n",
    "Poniżej robimy eksperyment:\n",
    "\n",
    "- generujemy 10 niezależnych zmiennych (każda z rozkładu normalnego),\n",
    "- liczymy macierz korelacji i odpowiadające jej statystyki t,\n",
    "- sprawdzamy, ile testów odrzuca $H_0$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50f82e33",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Dane: 35 obserwacji, 10 niezależnych zmiennych\n",
    "n = 35\n",
    "p = 10\n",
    "\n",
    "m = rng.normal(size=(n, p))\n",
    "\n",
    "corr = np.corrcoef(m, rowvar=False)\n",
    "\n",
    "# Statystyka t ma sens tylko dla par różnych zmiennych, nie na przekątnej.\n",
    "df = n - 2\n",
    "off_diagonal = ~np.eye(p, dtype=bool)\n",
    "t_matrix = np.full_like(corr, np.nan)\n",
    "t_matrix[off_diagonal] = (\n",
    "    corr[off_diagonal] / np.sqrt(1 - corr[off_diagonal] ** 2)\n",
    ") * np.sqrt(df)\n",
    "\n",
    "alpha = 0.05\n",
    "critical = stats.t.ppf(1 - alpha / 2, df=df)\n",
    "reject = np.zeros_like(corr, dtype=bool)\n",
    "reject[off_diagonal] = np.abs(t_matrix[off_diagonal]) >= critical\n",
    "\n",
    "print(\"Liczba odrzuceń H0 (poza przekątną):\", reject.sum())\n",
    "print(\"Oczekiwana liczba ~ alpha * liczba_testów =\", alpha * (p * (p - 1)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e4fb17b",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(6, 5))\n",
    "ax = sns.heatmap(reject.astype(int), cmap=\"viridis\", cbar=False, square=True)\n",
    "ax.set_title(\"Odrzucenia H0 dla korelacji (1 = odrzucono)\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61afec45",
   "metadata": {},
   "source": [
    "# R², diagnostyka i predykcja w regresji liniowej\n",
    "\n",
    "Wróćmy teraz do regresji jako modelu, a nie tylko linii narysowanej na wykresie.\n",
    "Po współczynniku kierunkowym i korelacji interesują nas jeszcze trzy pytania:\n",
    "\n",
    "- jaka część zmienności $Y$ jest opisywana przez model,\n",
    "- czy założenia regresji wyglądają rozsądnie na wykresach diagnostycznych,\n",
    "- jak odróżnić przewidywanie średniej od przewidywania pojedynczej obserwacji."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e2240b4",
   "metadata": {},
   "source": [
    "## Przykład: wynik SAT jako funkcja wyniku w teście\n",
    "\n",
    "Howell omawia tu dane z badania Katz, Lautenschlager, Blackburn i Harris (1990).\n",
    "Badacze interesowali się tym, na ile wynik w specjalnym teście opartym na pytaniach wielokrotnego wyboru z SAT wiąże się z rzeczywistym wynikiem SAT Verbal.\n",
    "\n",
    "> Czy osoby lepiej radzące sobie z tym typem zadań mają również wyższe wyniki SAT Verbal?\n",
    "\n",
    "Zmienna `Score` opisuje wynik w teście, a `SATV` - wynik SAT Verbal.\n",
    "Jeśli zależność jest dodatnia, to wyższy `Score` będzie wiązał się z wyższą przewidywaną wartością `SATV`.\n",
    "Jednocześnie nie chcielibyśmy interpretować tego naiwnie: część związku może odzwierciedlać ogólne kompetencje językowe, część obycie z formatem testu, a część inne cechy studentów.\n",
    "\n",
    "Najpierw obejrzymy dane i dopasujemy prosty model liniowy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4de6c599",
   "metadata": {},
   "outputs": [],
   "source": [
    "sat_df = pd.read_csv(DATA_DIR / \"Tab9-6.dat\", sep=\"\\t\")\n",
    "sat_df.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6af6b40",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_score = sat_df[\"Score\"]\n",
    "Y_satv = sat_df[\"SATV\"]\n",
    "\n",
    "model_sat = smf.ols(\"SATV ~ Score\", data=sat_df).fit()\n",
    "model_sat.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a002f35",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_grid = np.linspace(X_score.min(), X_score.max(), 200)\n",
    "y_grid = model_sat.predict(pd.DataFrame({\"Score\": x_grid}))\n",
    "\n",
    "plt.scatter(X_score, Y_satv, alpha=0.9)\n",
    "plt.plot(x_grid, y_grid, color=\"black\", linewidth=2)\n",
    "plt.xlabel(\"Wynik w teście wiedzy o muzyce\")\n",
    "plt.ylabel(\"SAT Verbal\")\n",
    "plt.title(\"Regresja liniowa: SATV ~ Score\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "815eae23",
   "metadata": {},
   "source": [
    "## $R^2$ „od podstaw”\n",
    "\n",
    "W prostej regresji liniowej:\n",
    "\n",
    "- $SST$ – całkowita zmienność $Y$ (suma kwadratów względem średniej),\n",
    "- $SSR$ – zmienność wyjaśniona przez model (suma kwadratów względem średniej),\n",
    "- $SSE$ – zmienność resztowa (błędy predykcji).\n",
    "\n",
    "Zachodzi:\n",
    "\n",
    "$$SST = SSR + SSE$$\n",
    "\n",
    "a:\n",
    "\n",
    "$$R^2 = \\frac{SSR}{SST} = 1 - \\frac{SSE}{SST}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04a9fca9",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = X_score\n",
    "y = Y_satv\n",
    "y_hat = model_sat.predict(pd.DataFrame({\"Score\": x}))\n",
    "resid = y - y_hat\n",
    "\n",
    "sst = np.sum((y - y.mean()) ** 2)\n",
    "sse = np.sum(resid**2)\n",
    "ssr = np.sum((y_hat - y.mean()) ** 2)\n",
    "r_squared_manual = ssr / sst\n",
    "\n",
    "print(\"SST (całkowita zmienność):\", sst)\n",
    "print(\"SSR (wyjaśniona zmienność):\", ssr)\n",
    "print(\"SSE (niewyjaśniona zmienność):\", sse)\n",
    "print(\"R² ręcznie:\", r_squared_manual)\n",
    "print(\"R² statsmodels:\", model_sat.rsquared)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff4afc97",
   "metadata": {},
   "source": [
    "## Założenia regresji (w wersji „na skróty”)\n",
    "\n",
    "W prostej regresji liniowej zwykle zakładamy m.in.:\n",
    "\n",
    "- relacja średniej $Y$ względem $X$ jest liniowa,\n",
    "- reszty mają stałą wariancję (homoscedastyczność),\n",
    "- reszty są w przybliżeniu normalne (przydatne dla wnioskowania o współczynnikach),\n",
    "- obserwacje są niezależne.\n",
    "\n",
    "Żeby nie polegać wyłącznie na wynikach testów, warto obejrzeć wykresy diagnostyczne.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ce7cb4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Podstawowe wykresy diagnostyczne: reszty vs dopasowania oraz Q-Q plot\n",
    "\n",
    "fitted = model_sat.fittedvalues\n",
    "resid = model_sat.resid\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(12, 4), constrained_layout=True)\n",
    "\n",
    "axes[0].scatter(fitted, resid, alpha=0.9)\n",
    "axes[0].axhline(0, color=\"black\", linewidth=1)\n",
    "axes[0].set_title(\"Residuals vs Fitted\")\n",
    "axes[0].set_xlabel(\"Fitted values\")\n",
    "axes[0].set_ylabel(\"Residuals\")\n",
    "\n",
    "stats.probplot(resid, dist=\"norm\", plot=axes[1])\n",
    "axes[1].set_title(\"Normal Q-Q (residuals)\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff854a9d",
   "metadata": {},
   "source": [
    "## Obserwacje odstające i wpływowe (leverage, Cook's distance)\n",
    "\n",
    "Niektóre punkty mogą mieć nieproporcjonalnie duży wpływ na dopasowaną linię.\n",
    "\n",
    "W prostym przykładzie sprawdzimy, które obserwacje mają największy dystans Cooka.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21d84d6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "infl = OLSInfluence(model_sat)\n",
    "cooks_d = infl.cooks_distance[0]\n",
    "\n",
    "influential = (\n",
    "    pd.DataFrame(\n",
    "        {\n",
    "            \"ID\": sat_df[\"ID\"],\n",
    "            \"Score\": X_score,\n",
    "            \"SATV\": Y_satv,\n",
    "            \"cooks_d\": cooks_d,\n",
    "        }\n",
    "    )\n",
    "    .sort_values(\"cooks_d\", ascending=False)\n",
    "    .head(8)\n",
    ")\n",
    "\n",
    "print(\"Największe Cook's distance (top 8):\")\n",
    "influential"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "457aef9e",
   "metadata": {},
   "source": [
    "## Predykcja\n",
    "\n",
    "Regresji możemy używać do przewidywania wartości $Y$.\n",
    "\n",
    "W praktyce pojawiają się dwa różne rodzaje przedziałów:\n",
    "\n",
    "- **przedział ufności dla średniej** (węższy): niepewność co do $E[Y|X]$,\n",
    "- **przedział predykcji dla obserwacji** (szerszy): dodatkowo uwzględnia losowy rozrzut pojedynczych obserwacji wokół średniej.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7be19932",
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_all = model_sat.get_prediction(\n",
    "    pd.DataFrame({\"Score\": X_score})\n",
    ").summary_frame(alpha=0.05)\n",
    "\n",
    "sat_plot = sat_df.copy()\n",
    "sat_plot[\"fitted\"] = pred_all[\"mean\"]\n",
    "sat_plot[\"mean_ci_low\"] = pred_all[\"mean_ci_lower\"]\n",
    "sat_plot[\"mean_ci_high\"] = pred_all[\"mean_ci_upper\"]\n",
    "sat_plot[\"obs_ci_low\"] = pred_all[\"obs_ci_lower\"]\n",
    "sat_plot[\"obs_ci_high\"] = pred_all[\"obs_ci_upper\"]\n",
    "\n",
    "sat_plot.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f76a0f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "score_value = 60\n",
    "new_score = pd.DataFrame({\"Score\": [score_value]})\n",
    "pred_60 = model_sat.get_prediction(new_score).summary_frame(alpha=0.05)\n",
    "pred_60"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0fcbacb3",
   "metadata": {},
   "outputs": [],
   "source": [
    "points = np.linspace(30, 70, 400)\n",
    "pred = model_sat.get_prediction(\n",
    "    pd.DataFrame({\"Score\": points})\n",
    ").summary_frame(alpha=0.05)\n",
    "\n",
    "plt.scatter(X_score, Y_satv, alpha=0.8, label=\"Dane\")\n",
    "plt.plot(points, pred[\"mean\"], color=\"black\", label=\"Średnia predykcja\")\n",
    "plt.fill_between(\n",
    "    points,\n",
    "    pred[\"mean_ci_lower\"],\n",
    "    pred[\"mean_ci_upper\"],\n",
    "    color=\"steelblue\",\n",
    "    alpha=0.25,\n",
    "    label=\"95% CI dla średniej\",\n",
    ")\n",
    "plt.fill_between(\n",
    "    points,\n",
    "    pred[\"obs_ci_lower\"],\n",
    "    pred[\"obs_ci_upper\"],\n",
    "    color=\"orange\",\n",
    "    alpha=0.18,\n",
    "    label=\"95% PI dla osoby\",\n",
    ")\n",
    "plt.axvline(score_value, color=\"gray\", linestyle=\"--\", linewidth=1)\n",
    "plt.xlabel(\"Wynik w teście wiedzy o muzyce\")\n",
    "plt.ylabel(\"SAT Verbal\")\n",
    "plt.title(\"Przedział ufności vs przedział predykcji\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff128989",
   "metadata": {},
   "source": [
    "# Dlaczego diagnostyka jest ważna? Dane łamiące założenia\n",
    "\n",
    "Zobaczmy dane, o których wiemy, że **nie spełniają założeń** regresji liniowej:\n",
    "\n",
    "- relacja jest nieliniowa (kwadratowa),\n",
    "- wariancja błędu rośnie z $X$,\n",
    "- rozkład błędów nie jest normalny.\n",
    "\n",
    "Mimo to prosta regresja może wyjść „istotna statystycznie” – i nie powinniśmy wtedy bezrefleksyjnie mówić o relacji liniowej.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17051e97",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Dane symulacyjne łamiące założenia\n",
    "x = rng.uniform(0, 10, size=120)\n",
    "noise = rng.uniform(0, 5 * x)\n",
    "y = x**2 + noise\n",
    "\n",
    "bad_df = pd.DataFrame({\"x\": x, \"y\": y})\n",
    "model_bad = smf.ols(\"y ~ x\", data=bad_df).fit()\n",
    "model_bad.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c53773f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_grid = np.linspace(x.min(), x.max(), 200)\n",
    "y_grid = model_bad.predict(pd.DataFrame({\"x\": x_grid}))\n",
    "\n",
    "plt.scatter(x, y, alpha=0.7)\n",
    "plt.plot(x_grid, y_grid, color=\"black\", linewidth=2)\n",
    "plt.xlabel(\"X\")\n",
    "plt.ylabel(\"Y\")\n",
    "plt.title(\"Model liniowy dopasowany do nieliniowej zależności\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d12d6b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Diagnostyka: reszty vs dopasowania i Q-Q\n",
    "\n",
    "fitted = model_bad.fittedvalues\n",
    "resid = model_bad.resid\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(12, 4), constrained_layout=True)\n",
    "\n",
    "axes[0].scatter(fitted, resid, alpha=0.7)\n",
    "axes[0].axhline(0, color=\"black\", linewidth=1)\n",
    "axes[0].set_title(\"Residuals vs Fitted\")\n",
    "axes[0].set_xlabel(\"Fitted values\")\n",
    "axes[0].set_ylabel(\"Residuals\")\n",
    "\n",
    "stats.probplot(resid, dist=\"norm\", plot=axes[1])\n",
    "axes[1].set_title(\"Normal Q-Q (residuals)\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "003e00d2",
   "metadata": {},
   "source": [
    "# Ćwiczenia (samodzielnie)\n",
    "\n",
    "Poniższe ćwiczenia korzystają z danych `data/africa_subsaharan.csv`.\n",
    "\n",
    "1. Narysuj wykres punktowy pokazujący relację między `income` i `infant_mortality`.\n",
    "2. Dodaj do wykresu linię regresji.\n",
    "3. Sprawdź, jak zmienia się dopasowanie, gdy porównasz wszystkie kraje z podzbiorem `in_analysis == 1`.\n",
    "4. Oblicz korelacje między wszystkimi zmiennymi z listy `africa_variables`.\n",
    "5. Zwizualizuj macierz korelacji.\n",
    "6. Jak duże musi być $|r|$, aby było istotne statystycznie dla $\\alpha=0.05$ i $\\alpha=0.01$ przy tej liczebności?\n",
    "7. Które zmienne są najsilniej związane z `infant_mortality`? Oceniaj po $|r|$, ale pamiętaj o wykresach punktowych.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05dd52bb",
   "metadata": {},
   "source": [
    "# Podsumowanie\n",
    "\n",
    "- Zaczynamy od wykresu: bez wizualizacji łatwo przegapić nieliniowość i obserwacje odstające.\n",
    "- Korelacja $r$ to wystandaryzowana kowariancja i mierzy siłę **liniowej** zależności.\n",
    "- Test dla korelacji opiera się na statystyce t z $n-2$ stopniami swobody.\n",
    "- W prostej regresji liniowej *slope* można policzyć „od podstaw” i zinterpretować w jednostkach $Y$ na 1 jednostkę $X$.\n",
    "- $R^2$ mówi, jaka część zmienności $Y$ jest wyjaśniana przez model (w prostej regresji: $R^2 = r^2$).\n",
    "- Predykcja ma dwa „poziomy niepewności”: CI dla średniej i PI dla obserwacji.\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
}