{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "176532a1",
   "metadata": {},
   "source": [
    "# Interakcje w modelu liniowym: zmienne ciągłe i nominalne\n",
    "\n",
    "W tej części zamykamy temat regresji liniowej, pokazując, co się dzieje, gdy efekt jednego predyktora **zależy od wartości drugiego predyktora**.\n",
    "\n",
    "To jest właśnie interakcja. W modelu bez interakcji efekty predyktorów są addytywne: każdy predyktor „dokłada” swoją część niezależnie od pozostałych. W modelu z interakcją ten sam predyktor może mieć inny efekt przy niskiej, średniej i wysokiej wartości innego predyktora.\n",
    "\n",
    "Najwięcej czasu poświęcimy interakcjom między zmiennymi ciągłymi, bo tam interpretacja współczynników jest najmniej intuicyjna. Na końcu wrócimy do zmiennych nominalnych i zobaczymy, że idea jest ta sama: interakcja to brak addytywności, tylko kodowany w inny sposób."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d3b675e",
   "metadata": {},
   "source": [
    "## Cele\n",
    "\n",
    "Po zajęciach:\n",
    "\n",
    "- Rozumiesz interakcję jako sytuację, w której efekt jednego predyktora zależy od drugiego predyktora.\n",
    "- Umiesz dopasować model z interakcją składnią `x * z`.\n",
    "- Rozumiesz, dlaczego przy interakcjach ciągłych zwykle centrujemy predyktory.\n",
    "- Umiesz odczytać proste nachylenia (*simple slopes*) dla różnych wartości moderatora.\n",
    "- Rozumiesz interakcję 2×2 jako różnicę różnic.\n",
    "- Umiesz połączyć regresyjne ujęcie interakcji z tabelą ANOVA."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1ef89d9",
   "metadata": {},
   "source": [
    "## Wymagania wstępne\n",
    "\n",
    "- Regresja liniowa prosta i wielokrotna.\n",
    "- Interpretacja współczynników regresji.\n",
    "- Testy t i F w modelu liniowym.\n",
    "- Podstawy pracy z formułami w `statsmodels`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "334db941",
   "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.anova import anova_lm\n",
    "from statsmodels.stats.outliers_influence import variance_inflation_factor\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": "66a1bc67",
   "metadata": {},
   "source": [
    "## Dane: wynagrodzenie, doświadczenie i staż\n",
    "\n",
    "Użyjemy zbioru `Salaries`, który opisuje wynagrodzenia akademickie w amerykańskim college'u w roku akademickim 2008/2009. Każdy wiersz opisuje jedną osobę.\n",
    "\n",
    "Najważniejsze zmienne:\n",
    "\n",
    "- `salary` — dziewięciomiesięczne wynagrodzenie w dolarach,\n",
    "- `yrs_since_phd` — liczba lat od uzyskania doktoratu,\n",
    "- `yrs_service` — staż pracy w obecnej instytucji,\n",
    "- `rank` — stanowisko akademickie,\n",
    "- `sex` — płeć osoby.\n",
    "\n",
    "Ten przykład jest dobry dydaktycznie, bo dwie zmienne doświadczenia są ze sobą silnie powiązane, ale nie znaczą dokładnie tego samego. Możemy więc zapytać: czy związek lat od doktoratu z wynagrodzeniem jest taki sam dla osób z krótkim i długim stażem w obecnej instytucji?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1adb343f",
   "metadata": {},
   "outputs": [],
   "source": [
    "DATA_DIR = Path(\"data\")\n",
    "\n",
    "salaries_raw = pd.read_csv(DATA_DIR / \"Salaries.csv\")\n",
    "\n",
    "salaries = salaries_raw.rename(\n",
    "    columns={\n",
    "        \"yrs.since.phd\": \"yrs_since_phd\",\n",
    "        \"yrs.service\": \"yrs_service\",\n",
    "    }\n",
    ").drop(columns=\"rownames\")\n",
    "\n",
    "rank_order = [\"AsstProf\", \"AssocProf\", \"Prof\"]\n",
    "sex_order = [\"Female\", \"Male\"]\n",
    "\n",
    "salaries[\"rank\"] = pd.Categorical(salaries[\"rank\"], categories=rank_order, ordered=True)\n",
    "salaries[\"sex\"] = pd.Categorical(salaries[\"sex\"], categories=sex_order)\n",
    "\n",
    "salaries.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da7277af",
   "metadata": {},
   "outputs": [],
   "source": [
    "salaries[[\"salary\", \"yrs_since_phd\", \"yrs_service\"]].describe().T"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70f9dd5e",
   "metadata": {},
   "source": [
    "## Relacje dwuzmiennowe\n",
    "\n",
    "Zacznijmy tak, jak zwykle zaczynamy analizę regresji: od wykresów rozrzutu.\n",
    "\n",
    "Na tym etapie patrzymy osobno na relację wynagrodzenia z latami od doktoratu i osobno na relację wynagrodzenia ze stażem. To jeszcze nie mówi nam, jak te dwie zmienne działają razem, ale daje pierwszy obraz danych."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "176b7fc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(ncols=2, figsize=(12, 4), sharey=True)\n",
    "\n",
    "sns.regplot(\n",
    "    data=salaries,\n",
    "    x=\"yrs_since_phd\",\n",
    "    y=\"salary\",\n",
    "    ax=axes[0],\n",
    "    scatter_kws={\"alpha\": 0.55, \"s\": 28},\n",
    "    line_kws={\"color\": \"black\", \"lw\": 2},\n",
    ")\n",
    "axes[0].set_title(\"Wynagrodzenie a lata od doktoratu\")\n",
    "axes[0].set_xlabel(\"Lata od doktoratu\")\n",
    "axes[0].set_ylabel(\"Wynagrodzenie\")\n",
    "\n",
    "sns.regplot(\n",
    "    data=salaries,\n",
    "    x=\"yrs_service\",\n",
    "    y=\"salary\",\n",
    "    ax=axes[1],\n",
    "    scatter_kws={\"alpha\": 0.55, \"s\": 28},\n",
    "    line_kws={\"color\": \"black\", \"lw\": 2},\n",
    ")\n",
    "axes[1].set_title(\"Wynagrodzenie a staż w instytucji\")\n",
    "axes[1].set_xlabel(\"Staż w obecnej instytucji\")\n",
    "axes[1].set_ylabel(\"\")\n",
    "\n",
    "for ax in axes:\n",
    "    ax.set_xlim(0, 62)\n",
    "    ax.set_ylim(50_000, 240_000)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65af8a43",
   "metadata": {},
   "source": [
    "## Regresja bez interakcji\n",
    "\n",
    "Najpierw dopasujemy model addytywny:\n",
    "\n",
    "$$\n",
    "\\widehat{salary} = b_0 + b_1 \\cdot yrs\\_since\\_phd + b_2 \\cdot yrs\\_service\n",
    "$$\n",
    "\n",
    "W takim modelu efekt `yrs_since_phd` jest taki sam niezależnie od wartości `yrs_service`, a efekt `yrs_service` jest taki sam niezależnie od wartości `yrs_since_phd`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "665c65f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "fit_additive = smf.ols(\n",
    "    \"salary ~ yrs_since_phd + yrs_service\",\n",
    "    data=salaries,\n",
    ").fit()\n",
    "\n",
    "fit_additive.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb4bd6d6",
   "metadata": {},
   "source": [
    "Dwie zmienne doświadczenia są ze sobą powiązane, więc warto sprawdzić VIF. Nie traktujemy VIF jako magicznego testu, ale jako informację o tym, jak mocno predyktor jest przewidywalny przez pozostałe predyktory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e974813c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def vif_table(fit) -> pd.DataFrame:\n",
    "    exog = pd.DataFrame(fit.model.exog, columns=fit.model.exog_names)\n",
    "\n",
    "    rows = []\n",
    "    for idx, name in enumerate(exog.columns):\n",
    "        if name == \"Intercept\":\n",
    "            continue\n",
    "\n",
    "        vif = variance_inflation_factor(exog.to_numpy(), idx)\n",
    "        rows.append(\n",
    "            {\n",
    "                \"predyktor\": name,\n",
    "                \"VIF\": vif,\n",
    "                \"tolerancja\": 1 / vif,\n",
    "            }\n",
    "        )\n",
    "\n",
    "    return pd.DataFrame(rows)\n",
    "\n",
    "\n",
    "vif_table(fit_additive)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25473732",
   "metadata": {},
   "source": [
    "## Regresja z interakcją\n",
    "\n",
    "Teraz dopuszczamy możliwość, że efekt lat od doktoratu zależy od stażu w instytucji.\n",
    "\n",
    "W formule `statsmodels` zapis:\n",
    "\n",
    "```python\n",
    "salary ~ yrs_since_phd * yrs_service\n",
    "```\n",
    "\n",
    "oznacza:\n",
    "\n",
    "```python\n",
    "salary ~ yrs_since_phd + yrs_service + yrs_since_phd:yrs_service\n",
    "```\n",
    "\n",
    "Czyli model zawiera dwa efekty główne i jeden składnik interakcyjny."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "137ac02f",
   "metadata": {},
   "outputs": [],
   "source": [
    "fit_interaction = smf.ols(\n",
    "    \"salary ~ yrs_since_phd * yrs_service\",\n",
    "    data=salaries,\n",
    ").fit()\n",
    "\n",
    "fit_interaction.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75d06356",
   "metadata": {},
   "source": [
    "Interakcję między dwiema zmiennymi ciągłymi można rozumieć jako dodanie do modelu nowej zmiennej, która jest iloczynem tych zmiennych. Składnia `*` jest więc wygodnym skrótem, ale nie magią."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2cd7f2e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "salaries_product = salaries.assign(\n",
    "    interaction_product=salaries[\"yrs_since_phd\"] * salaries[\"yrs_service\"]\n",
    ")\n",
    "\n",
    "fit_product = smf.ols(\n",
    "    \"salary ~ yrs_since_phd + yrs_service + interaction_product\",\n",
    "    data=salaries_product,\n",
    ").fit()\n",
    "\n",
    "pd.DataFrame(\n",
    "    {\n",
    "        \"model z *\": fit_interaction.params,\n",
    "        \"model z ręcznym iloczynem\": fit_product.params.rename(\n",
    "            {\n",
    "                \"interaction_product\": \"yrs_since_phd:yrs_service\",\n",
    "            }\n",
    "        ),\n",
    "    }\n",
    ").round(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3cdfc4a3",
   "metadata": {},
   "source": [
    "## Co oznaczają współczynniki w modelu z interakcją?\n",
    "\n",
    "Dla dwóch predyktorów ciągłych model ma postać:\n",
    "\n",
    "$$\n",
    "\\widehat{Y} = b_0 + b_1X + b_2Z + b_3XZ\n",
    "$$\n",
    "\n",
    "Wtedy efekt $X$ nie jest już jedną stałą liczbą. Nachylenie dla $X$ przy konkretnej wartości $Z$ wynosi:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial \\widehat{Y}}{\\partial X} = b_1 + b_3Z\n",
    "$$\n",
    "\n",
    "To jest najważniejszy fakt w tej części wykładu: **współczynnik przy jednym predyktorze mówi o jego efekcie wtedy, gdy drugi predyktor ma wartość 0**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6cca5e43",
   "metadata": {},
   "source": [
    "## Centrowanie zmiennych\n",
    "\n",
    "Wartość 0 często nie jest sensownym punktem odniesienia. W naszym przykładzie `yrs_service = 0` oznacza osobę bez stażu w obecnej instytucji. To może być możliwe, ale niekoniecznie jest punktem, wokół którego chcemy interpretować całą relację.\n",
    "\n",
    "Dlatego przy interakcjach ciągłych zwykle centrujemy predyktory: od każdej wartości odejmujemy średnią. Wtedy 0 oznacza „wartość przeciętną w próbie”."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b71f2b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "salaries = salaries.assign(\n",
    "    yrs_since_phd_c=salaries[\"yrs_since_phd\"] - salaries[\"yrs_since_phd\"].mean(),\n",
    "    yrs_service_c=salaries[\"yrs_service\"] - salaries[\"yrs_service\"].mean(),\n",
    ")\n",
    "\n",
    "fit_centered = smf.ols(\n",
    "    \"salary ~ yrs_since_phd_c * yrs_service_c\",\n",
    "    data=salaries,\n",
    ").fit()\n",
    "\n",
    "fit_centered.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5eadcf4",
   "metadata": {},
   "source": [
    "Centrowanie nie zmienia dopasowanych wartości ani $R^2$. Zmienia natomiast interpretację wyrazu wolnego i efektów głównych, bo zmienia punkt odniesienia dla wartości 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4dcaee27",
   "metadata": {},
   "outputs": [],
   "source": [
    "comparison = pd.DataFrame(\n",
    "    {\n",
    "        \"bez centrowania\": fit_interaction.params,\n",
    "        \"po centrowaniu\": fit_centered.params.rename(\n",
    "            {\n",
    "                \"yrs_since_phd_c\": \"yrs_since_phd\",\n",
    "                \"yrs_service_c\": \"yrs_service\",\n",
    "                \"yrs_since_phd_c:yrs_service_c\": \"yrs_since_phd:yrs_service\",\n",
    "            }\n",
    "        ),\n",
    "    }\n",
    ")\n",
    "\n",
    "comparison.loc[\"R^2\", :] = [fit_interaction.rsquared, fit_centered.rsquared]\n",
    "comparison.round(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6987deb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "vif_uncentered = vif_table(fit_interaction).assign(model=\"bez centrowania\")\n",
    "vif_centered = vif_table(fit_centered).assign(model=\"po centrowaniu\")\n",
    "\n",
    "pd.concat([vif_uncentered, vif_centered], ignore_index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90fa3fa5",
   "metadata": {},
   "source": [
    "## Proste nachylenia (*simple slopes*)\n",
    "\n",
    "Skoro efekt `yrs_since_phd` zależy od `yrs_service`, możemy policzyć ten efekt dla kilku konkretnych wartości moderatora.\n",
    "\n",
    "Typowy wybór dydaktyczny to:\n",
    "\n",
    "- niski staż: $-1SD$,\n",
    "- przeciętny staż: średnia,\n",
    "- wysoki staż: $+1SD$.\n",
    "\n",
    "Dla wycentrowanego moderatora są to po prostu wartości `-SD`, `0`, `+SD`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a91bbf4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def simple_slope_table(\n",
    "    fit,\n",
    "    moderator_values: pd.DataFrame,\n",
    "    pred_name: str,\n",
    "    mod_name: str,\n",
    "    interaction_name: str,\n",
    ") -> pd.DataFrame:\n",
    "    params = fit.params\n",
    "    cov = fit.cov_params()\n",
    "    df_resid = fit.df_resid\n",
    "\n",
    "    rows = []\n",
    "    for _, row in moderator_values.iterrows():\n",
    "        z_value = row[mod_name]\n",
    "        slope = params[pred_name] + params[interaction_name] * z_value\n",
    "\n",
    "        var_slope = (\n",
    "            cov.loc[pred_name, pred_name]\n",
    "            + z_value**2 * cov.loc[interaction_name, interaction_name]\n",
    "            + 2 * z_value * cov.loc[pred_name, interaction_name]\n",
    "        )\n",
    "        se = np.sqrt(var_slope)\n",
    "        t_value = slope / se\n",
    "        p_value = 2 * stats.t.sf(abs(t_value), df=df_resid)\n",
    "        t_crit = stats.t.isf(0.05 / 2, df=df_resid)\n",
    "\n",
    "        rows.append(\n",
    "            {\n",
    "                \"poziom moderatora\": row[\"poziom moderatora\"],\n",
    "                \"yrs_service\": row[\"yrs_service\"],\n",
    "                \"nachylenie yrs_since_phd\": slope,\n",
    "                \"SE\": se,\n",
    "                \"t\": t_value,\n",
    "                \"p\": p_value,\n",
    "                \"CI 95% low\": slope - t_crit * se,\n",
    "                \"CI 95% high\": slope + t_crit * se,\n",
    "            }\n",
    "        )\n",
    "\n",
    "    return pd.DataFrame(rows)\n",
    "\n",
    "\n",
    "service_mean = salaries[\"yrs_service\"].mean()\n",
    "service_sd = salaries[\"yrs_service\"].std(ddof=1)\n",
    "\n",
    "moderator_values = pd.DataFrame(\n",
    "    {\n",
    "        \"poziom moderatora\": [\"niski staż (-1 SD)\", \"średni staż\", \"wysoki staż (+1 SD)\"],\n",
    "        \"yrs_service\": [service_mean - service_sd, service_mean, service_mean + service_sd],\n",
    "        \"yrs_service_c\": [-service_sd, 0, service_sd],\n",
    "    }\n",
    ")\n",
    "\n",
    "simple_slope_table(\n",
    "    fit=fit_centered,\n",
    "    moderator_values=moderator_values,\n",
    "    pred_name=\"yrs_since_phd_c\",\n",
    "    mod_name=\"yrs_service_c\",\n",
    "    interaction_name=\"yrs_since_phd_c:yrs_service_c\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58f1def9",
   "metadata": {},
   "source": [
    "## Wizualizacja interakcji ciągłej\n",
    "\n",
    "Tabela współczynników mówi nam, czy składnik interakcyjny jest różny od zera. Wykres mówi nam, **jak** wygląda ta interakcja.\n",
    "\n",
    "Narysujemy przewidywane wynagrodzenie jako funkcję lat od doktoratu dla trzech poziomów stażu: niskiego, średniego i wysokiego."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee0e97ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_raw = np.linspace(0, 60, 140)\n",
    "x_centered = x_raw - salaries[\"yrs_since_phd\"].mean()\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(9, 5))\n",
    "\n",
    "ax.scatter(\n",
    "    salaries[\"yrs_since_phd\"],\n",
    "    salaries[\"salary\"],\n",
    "    color=\"0.70\",\n",
    "    alpha=0.45,\n",
    "    s=24,\n",
    "    label=\"obserwacje\",\n",
    ")\n",
    "\n",
    "colors = [\"#4477AA\", \"#228833\", \"#CC6677\"]\n",
    "\n",
    "for color, (_, row) in zip(colors, moderator_values.iterrows(), strict=True):\n",
    "    new_data = pd.DataFrame(\n",
    "        {\n",
    "            \"yrs_since_phd_c\": x_centered,\n",
    "            \"yrs_service_c\": row[\"yrs_service_c\"],\n",
    "        }\n",
    "    )\n",
    "    pred = fit_centered.get_prediction(new_data).summary_frame()\n",
    "\n",
    "    ax.plot(\n",
    "        x_raw,\n",
    "        pred[\"mean\"],\n",
    "        color=color,\n",
    "        lw=2.5,\n",
    "        label=f\"{row['poziom moderatora']} ({row['yrs_service']:.1f} lat)\",\n",
    "    )\n",
    "    ax.fill_between(\n",
    "        x_raw,\n",
    "        pred[\"mean_ci_lower\"],\n",
    "        pred[\"mean_ci_upper\"],\n",
    "        color=color,\n",
    "        alpha=0.12,\n",
    "    )\n",
    "\n",
    "ax.set_xlim(0, 60)\n",
    "ax.set_ylim(50_000, 240_000)\n",
    "ax.set_title(\"Interakcja: efekt lat od doktoratu zależy od stażu\")\n",
    "ax.set_xlabel(\"Lata od doktoratu\")\n",
    "ax.set_ylabel(\"Przewidywane wynagrodzenie\")\n",
    "ax.legend(loc=\"upper left\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e2efd7a",
   "metadata": {},
   "source": [
    "### Widżet: interakcja między zmiennymi ciągłymi\n",
    "\n",
    "Ten widżet pokazuje samą logikę interakcji. Zmieniamy trzy składniki modelu:\n",
    "\n",
    "- efekt `X`,\n",
    "- efekt moderatora `Z`,\n",
    "- interakcję `X × Z`.\n",
    "\n",
    "Skala osi jest stała. Dzięki temu widać, jak manipulacja parametrami zmienia układ linii, zamiast tylko przeskalowywać wykres."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05becc79",
   "metadata": {},
   "outputs": [],
   "source": [
    "def rysuj_widget_interakcji_ciaglej(\n",
    "    efekt_x: float = 900,\n",
    "    efekt_z: float = 250,\n",
    "    interakcja_xz: float = -20,\n",
    ") -> None:\n",
    "    x_raw = np.linspace(0, 60, 140)\n",
    "    x_centered = x_raw - 30\n",
    "\n",
    "    moderator_levels = {\n",
    "        \"niski Z (-1 SD)\": -15,\n",
    "        \"średni Z\": 0,\n",
    "        \"wysoki Z (+1 SD)\": 15,\n",
    "    }\n",
    "    colors = {\n",
    "        \"niski Z (-1 SD)\": \"#4477AA\",\n",
    "        \"średni Z\": \"#228833\",\n",
    "        \"wysoki Z (+1 SD)\": \"#CC6677\",\n",
    "    }\n",
    "\n",
    "    fig, ax = plt.subplots(figsize=(8.5, 4.8))\n",
    "\n",
    "    opis = []\n",
    "    for label, z_value in moderator_levels.items():\n",
    "        y_hat = 110_000 + efekt_x * x_centered + efekt_z * z_value + interakcja_xz * x_centered * z_value\n",
    "        slope = efekt_x + interakcja_xz * z_value\n",
    "        opis.append(f\"{label}: nachylenie X = {slope:.0f}\")\n",
    "        ax.plot(x_raw, y_hat, color=colors[label], lw=2.6, label=label)\n",
    "\n",
    "    ax.set_xlim(0, 60)\n",
    "    ax.set_ylim(50_000, 240_000)\n",
    "    ax.set_title(\"Interakcja ciągła × ciągła: zmiana nachylenia\")\n",
    "    ax.set_xlabel(\"X\")\n",
    "    ax.set_ylabel(\"Przewidywane Y\")\n",
    "    ax.legend(loc=\"upper left\")\n",
    "\n",
    "    ax.text(\n",
    "        0.98,\n",
    "        0.05,\n",
    "        \"\\n\".join(opis),\n",
    "        transform=ax.transAxes,\n",
    "        ha=\"right\",\n",
    "        va=\"bottom\",\n",
    "        bbox={\"boxstyle\": \"round,pad=0.45\", \"fc\": \"white\", \"ec\": \"0.35\"},\n",
    "    )\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "if WIDGETY_DOSTEPNE:\n",
    "    widget_style = {\"description_width\": \"initial\"}\n",
    "    controls_ciagle = {\n",
    "        \"efekt_x\": widgets.FloatSlider(\n",
    "            value=900,\n",
    "            min=-2500,\n",
    "            max=3500,\n",
    "            step=100,\n",
    "            description=\"efekt X\",\n",
    "            continuous_update=False,\n",
    "            style=widget_style,\n",
    "        ),\n",
    "        \"efekt_z\": widgets.FloatSlider(\n",
    "            value=250,\n",
    "            min=-2500,\n",
    "            max=2500,\n",
    "            step=100,\n",
    "            description=\"efekt Z\",\n",
    "            continuous_update=False,\n",
    "            style=widget_style,\n",
    "        ),\n",
    "        \"interakcja_xz\": widgets.FloatSlider(\n",
    "            value=-20,\n",
    "            min=-150,\n",
    "            max=150,\n",
    "            step=5,\n",
    "            description=\"interakcja X×Z\",\n",
    "            continuous_update=False,\n",
    "            style=widget_style,\n",
    "        ),\n",
    "    }\n",
    "    widget_ciagle = widgets.interactive_output(\n",
    "        rysuj_widget_interakcji_ciaglej,\n",
    "        controls_ciagle,\n",
    "    )\n",
    "    display(widgets.HBox([widgets.VBox(list(controls_ciagle.values())), widget_ciagle]))\n",
    "else:\n",
    "    rysuj_widget_interakcji_ciaglej()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15630c52",
   "metadata": {},
   "source": [
    "## Interakcje nominalne: ta sama idea, inne kodowanie\n",
    "\n",
    "Przy zmiennych nominalnych nie mamy naturalnej osi liczbowej, więc poziomy czynnika muszą zostać zakodowane. Sama idea interakcji pozostaje jednak taka sama: efekt jednego czynnika zależy od poziomu drugiego czynnika.\n",
    "\n",
    "W tym samym zbiorze danych możemy zapytać, czy różnice wynagrodzeń między stanowiskami wyglądają podobnie dla kobiet i mężczyzn. W schemacie 2×2 najprostsza interpretacja brzmi: **interakcja to różnica różnic**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc9f615e",
   "metadata": {},
   "outputs": [],
   "source": [
    "salaries.groupby([\"rank\", \"sex\"], observed=True)[\"salary\"].agg([\"count\", \"mean\", \"std\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36d281a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(9, 4.5))\n",
    "\n",
    "sns.pointplot(\n",
    "    data=salaries,\n",
    "    x=\"rank\",\n",
    "    y=\"salary\",\n",
    "    hue=\"sex\",\n",
    "    order=rank_order,\n",
    "    hue_order=sex_order,\n",
    "    errorbar=(\"ci\", 95),\n",
    "    dodge=True,\n",
    "    markers=[\"o\", \"s\"],\n",
    ")\n",
    "\n",
    "sns.stripplot(\n",
    "    data=salaries,\n",
    "    x=\"rank\",\n",
    "    y=\"salary\",\n",
    "    hue=\"sex\",\n",
    "    order=rank_order,\n",
    "    hue_order=sex_order,\n",
    "    dodge=True,\n",
    "    alpha=0.28,\n",
    "    size=3,\n",
    "    color=\"black\",\n",
    ")\n",
    "\n",
    "handles, labels = plt.gca().get_legend_handles_labels()\n",
    "plt.legend(handles[:2], labels[:2], title=\"sex\", loc=\"upper left\")\n",
    "plt.ylim(50_000, 240_000)\n",
    "plt.title(\"Średnie wynagrodzenie w układzie rank × sex\")\n",
    "plt.xlabel(\"rank\")\n",
    "plt.ylabel(\"salary\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a34c7376",
   "metadata": {},
   "source": [
    "### Widżet: interakcja 2×2 jako różnica różnic\n",
    "\n",
    "W schemacie 2×2 interakcja jest równa temu, że różnica między poziomami jednego czynnika nie jest taka sama dla obu poziomów drugiego czynnika.\n",
    "\n",
    "Oś Y ma stałą skalę wynagrodzeń. Gdy interakcja jest równa 0, linie są równoległe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbff27f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def rysuj_widget_interakcji_2x2(\n",
    "    efekt_rank: float = 22_000,\n",
    "    efekt_sex: float = 4_000,\n",
    "    interakcja: float = 0.0,\n",
    ") -> None:\n",
    "    rank_levels = [\"AsstProf\", \"Prof\"]\n",
    "    sex_levels = [\"Female\", \"Male\"]\n",
    "    rank_code = {\"AsstProf\": -1, \"Prof\": 1}\n",
    "    sex_code = {\"Female\": -1, \"Male\": 1}\n",
    "\n",
    "    grand_mean = 110_000\n",
    "    rows = []\n",
    "    for rank in rank_levels:\n",
    "        for sex in sex_levels:\n",
    "            mean_value = (\n",
    "                grand_mean\n",
    "                + efekt_rank * rank_code[rank]\n",
    "                + efekt_sex * sex_code[sex]\n",
    "                + interakcja * rank_code[rank] * sex_code[sex]\n",
    "            )\n",
    "            rows.append({\"rank\": rank, \"sex\": sex, \"mean\": mean_value})\n",
    "\n",
    "    widget_df = pd.DataFrame(rows)\n",
    "\n",
    "    fig, ax = plt.subplots(figsize=(7.5, 4.5))\n",
    "    colors = {\"Female\": \"#4477AA\", \"Male\": \"#CC6677\"}\n",
    "    markers = {\"Female\": \"o\", \"Male\": \"s\"}\n",
    "\n",
    "    for sex in sex_levels:\n",
    "        subset = widget_df[widget_df[\"sex\"] == sex]\n",
    "        ax.plot(\n",
    "            rank_levels,\n",
    "            subset[\"mean\"],\n",
    "            marker=markers[sex],\n",
    "            color=colors[sex],\n",
    "            lw=2.5,\n",
    "            ms=8,\n",
    "            label=sex,\n",
    "        )\n",
    "\n",
    "    diff_female = float(\n",
    "        widget_df.query(\"sex == 'Female' and rank == 'Prof'\")[\"mean\"].iloc[0]\n",
    "        - widget_df.query(\"sex == 'Female' and rank == 'AsstProf'\")[\"mean\"].iloc[0]\n",
    "    )\n",
    "    diff_male = float(\n",
    "        widget_df.query(\"sex == 'Male' and rank == 'Prof'\")[\"mean\"].iloc[0]\n",
    "        - widget_df.query(\"sex == 'Male' and rank == 'AsstProf'\")[\"mean\"].iloc[0]\n",
    "    )\n",
    "    diff_of_diffs = diff_male - diff_female\n",
    "\n",
    "    ax.set_ylim(50_000, 240_000)\n",
    "    ax.set_xlabel(\"rank\")\n",
    "    ax.set_ylabel(\"średnie salary\")\n",
    "    ax.set_title(\"Interakcja nominalna jako różnica różnic\")\n",
    "    ax.legend(title=\"sex\", loc=\"upper left\")\n",
    "\n",
    "    opis = (\n",
    "        f\"różnica rank u kobiet: {diff_female:,.0f}\\n\"\n",
    "        f\"różnica rank u mężczyzn: {diff_male:,.0f}\\n\"\n",
    "        f\"różnica różnic: {diff_of_diffs:,.0f}\"\n",
    "    )\n",
    "    ax.text(\n",
    "        0.98,\n",
    "        0.05,\n",
    "        opis,\n",
    "        transform=ax.transAxes,\n",
    "        ha=\"right\",\n",
    "        va=\"bottom\",\n",
    "        bbox={\"boxstyle\": \"round,pad=0.45\", \"fc\": \"white\", \"ec\": \"0.35\"},\n",
    "    )\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "if WIDGETY_DOSTEPNE:\n",
    "    widget_style = {\"description_width\": \"initial\"}\n",
    "    controls_nominalny = {\n",
    "        \"efekt_rank\": widgets.FloatSlider(\n",
    "            value=22_000,\n",
    "            min=-40_000,\n",
    "            max=60_000,\n",
    "            step=2_000,\n",
    "            description=\"efekt rank\",\n",
    "            continuous_update=False,\n",
    "            style=widget_style,\n",
    "        ),\n",
    "        \"efekt_sex\": widgets.FloatSlider(\n",
    "            value=4_000,\n",
    "            min=-25_000,\n",
    "            max=25_000,\n",
    "            step=1_000,\n",
    "            description=\"efekt sex\",\n",
    "            continuous_update=False,\n",
    "            style=widget_style,\n",
    "        ),\n",
    "        \"interakcja\": widgets.FloatSlider(\n",
    "            value=0.0,\n",
    "            min=-25_000,\n",
    "            max=25_000,\n",
    "            step=1_000,\n",
    "            description=\"interakcja\",\n",
    "            continuous_update=False,\n",
    "            style=widget_style,\n",
    "        ),\n",
    "    }\n",
    "    widget_nominalny = widgets.interactive_output(\n",
    "        rysuj_widget_interakcji_2x2,\n",
    "        controls_nominalny,\n",
    "    )\n",
    "    display(widgets.HBox([widgets.VBox(list(controls_nominalny.values())), widget_nominalny]))\n",
    "else:\n",
    "    rysuj_widget_interakcji_2x2()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64435b13",
   "metadata": {},
   "source": [
    "## Model nominalny z interakcją\n",
    "\n",
    "Dopasujemy najpierw model dla uproszczonego schematu 2×2: `AsstProf` vs `Prof` oraz `Female` vs `Male`. Użyjemy kodowania `Sum`. Dzięki temu efekty główne są interpretowane jako odchylenia od średniej ze średnich komórkowych, a nie jako porównania do jednej kategorii referencyjnej.\n",
    "\n",
    "Formuła:\n",
    "\n",
    "```python\n",
    "salary ~ C(rank, Sum) * C(sex, Sum)\n",
    "```\n",
    "\n",
    "zawiera efekty główne `rank`, `sex` oraz interakcję `rank:sex`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d2ddc41",
   "metadata": {},
   "outputs": [],
   "source": [
    "salaries_2x2 = salaries[salaries[\"rank\"].isin([\"AsstProf\", \"Prof\"])].copy()\n",
    "salaries_2x2[\"rank\"] = salaries_2x2[\"rank\"].cat.remove_unused_categories()\n",
    "\n",
    "fit_nominal = smf.ols(\n",
    "    \"salary ~ C(rank, Sum) * C(sex, Sum)\",\n",
    "    data=salaries_2x2,\n",
    ").fit()\n",
    "\n",
    "fit_nominal.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c82b52a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "anova_nominal = anova_lm(fit_nominal, typ=2)\n",
    "anova_nominal"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bd24e9a",
   "metadata": {},
   "source": [
    "W schemacie 2×2 możemy też policzyć interakcję bardzo bezpośrednio jako różnicę różnic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "698bfb50",
   "metadata": {},
   "outputs": [],
   "source": [
    "means_2x2 = salaries_2x2.pivot_table(\n",
    "    index=\"rank\",\n",
    "    columns=\"sex\",\n",
    "    values=\"salary\",\n",
    "    aggfunc=\"mean\",\n",
    "    observed=True,\n",
    ")\n",
    "\n",
    "diff_female = means_2x2.loc[\"Prof\", \"Female\"] - means_2x2.loc[\"AsstProf\", \"Female\"]\n",
    "diff_male = means_2x2.loc[\"Prof\", \"Male\"] - means_2x2.loc[\"AsstProf\", \"Male\"]\n",
    "\n",
    "pd.DataFrame(\n",
    "    {\n",
    "        \"porównanie\": [\n",
    "            \"Prof - AsstProf u kobiet\",\n",
    "            \"Prof - AsstProf u mężczyzn\",\n",
    "            \"różnica różnic\",\n",
    "        ],\n",
    "        \"wartość\": [diff_female, diff_male, diff_male - diff_female],\n",
    "    }\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8419bb5",
   "metadata": {},
   "source": [
    "## Gdy czynnik ma więcej niż dwa poziomy\n",
    "\n",
    "Dla pełnego schematu `rank × sex` logika jest ta sama, ale tabela ANOVA staje się wygodniejsza niż ręczne liczenie jednej różnicy różnic.\n",
    "\n",
    "W niezrównoważonych danych typ sum kwadratów ma znaczenie, bo różne typy odpowiadają różnym pytaniom. W praktyce trzeba wiedzieć, czy testujemy efekt sekwencyjnie, po uwzględnieniu innych efektów głównych, czy w pełnym modelu. Tutaj pokazujemy typ II i typ III jako dwa często spotykane warianty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a7d5a04",
   "metadata": {},
   "outputs": [],
   "source": [
    "fit_nominal_2x3 = smf.ols(\n",
    "    \"salary ~ C(rank, Sum) * C(sex, Sum)\",\n",
    "    data=salaries,\n",
    ").fit()\n",
    "\n",
    "display(anova_lm(fit_nominal_2x3, typ=2))\n",
    "display(anova_lm(fit_nominal_2x3, typ=3))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b0ca752",
   "metadata": {},
   "source": [
    "## Podsumowanie\n",
    "\n",
    "- Interakcja oznacza, że efekt jednego predyktora zależy od drugiego predyktora.\n",
    "- W modelu `Y ~ X * Z` składnik `X:Z` mówi, czy nachylenie dla `X` zmienia się wraz z `Z`.\n",
    "- Przy interakcjach ciągłych efekty główne odnoszą się do sytuacji, w której drugi predyktor ma wartość 0.\n",
    "- Centrowanie sprawia, że 0 oznacza wartość przeciętną, więc efekty główne są zwykle łatwiejsze do zinterpretowania.\n",
    "- Proste nachylenia pokazują efekt jednego predyktora przy konkretnych wartościach moderatora.\n",
    "- Przy zmiennych nominalnych interakcję 2×2 można czytać jako różnicę różnic.\n",
    "- W danych niezrównoważonych trzeba świadomie wybierać typ sum kwadratów w ANOVA."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}