{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9c36a211",
   "metadata": {},
   "source": [
    "# Analiza wariancji (ANOVA)\n",
    "\n",
    "Na tym wykładzie zajmiemy się **jednoczynnikową analizą wariancji** (*one-way ANOVA*). Zobaczymy:\n",
    "\n",
    "- skąd bierze się statystyka testowa $F$ i jak policzyć ją **„na piechotę”**,\n",
    "- jak czytać p-wartość w analizie wariancji (i czego ona *nie* mówi),\n",
    "- jak sprawdzić założenia w praktyce,\n",
    "- jak zrobić analizę *post-hoc* (Tukey HSD),\n",
    "- jak połączyć ANOVA z modelem liniowym (OLS): to ta sama rodzina metod."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f23089d1",
   "metadata": {},
   "source": [
    "## Cele\n",
    "\n",
    "Po zajęciach:\n",
    "\n",
    "- Umiesz sformułować hipotezy $H_0/H_A$ dla jednoczynnikowej ANOVA.\n",
    "- Umiesz policzyć składniki ANOVA „na piechotę”: $SSE$, $SSA$, $MSE$, $MSA$, a następnie $F$.\n",
    "- Umiesz policzyć p-wartość z rozkładu $F$ i wyznaczyć region krytyczny.\n",
    "- Umiesz wykonać ANOVA w bibliotece (weryfikacja wyniku), wykonać test Tukeya i zinterpretować wyniki.\n",
    "- Potrafisz raportować wynik: $F(df_1, df_2)$, p-wartość oraz miarę wielkości efektu (np. $\\eta^2$).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53ea2556",
   "metadata": {},
   "outputs": [],
   "source": [
    "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.formula.api as smf\n",
    "from statsmodels.stats.anova import anova_lm\n",
    "from statsmodels.stats.multicomp import pairwise_tukeyhsd\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": "117d2b7f",
   "metadata": {},
   "source": [
    "Na początek kilka ważnych faktów o analizie wariancji.\n",
    "\n",
    "- Jest powszechnie stosowaną metodą statystyczną pozwalającą na ocenę istotności różnic **wielu średnich**.\n",
    "- Opracował ją i upowszechnił Ronald Fisher w latach 20. XX wieku.\n",
    "- Jest to technika badania wyników (doświadczeń, obserwacji), które zależą od jednego (jednoczynnikowa analiza wariancji) lub kilku czynników działających równocześnie (dwuczynnikowa lub wieloczynnikowa analiza wariancji).\n",
    "  - przykłady takich czynników: metody leczenia, płeć, sposoby nawożenia.\n",
    "  - czynniki takie nazywamy *czynnikami grupującymi, klasyfikacyjnymi* lub *zmiennymi manipulacyjnymi*.\n",
    "  - każdy czynnik ma kilka poziomów (np. mało–średnio–dużo czegoś).\n",
    "- Założenia ANOVA są podobne jak w regresji liniowej: (1) homogeniczność wariancji w grupach, (2) normalność w grupach (częściej: normalność reszt w modelu) oraz (3) niezależność obserwacji.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8450dce",
   "metadata": {},
   "source": [
    "## Model teoretyczny\n",
    "\n",
    "Model strukturalny leżący u podstaw analizy wariancji ma postać:\n",
    "\n",
    "$$X_{ij} = \\mu + \\tau_j + \\epsilon_{ij}$$\n",
    "\n",
    "- $X_{ij}$ – indywidualna obserwacja numer $i$ należąca do grupy $j$,\n",
    "- $\\mu$ – ogólna średnia,\n",
    "- $\\tau_j$ – efekt przynależenia do grupy $j$,\n",
    "- $\\epsilon_{ij}$ – składnik losowy (błąd).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "249d411d",
   "metadata": {},
   "source": [
    "## Układ hipotez przy analizie wariancji\n",
    "\n",
    "W jednoczynnikowej ANOVA testujemy układ hipotez:\n",
    "\n",
    "$$H_0: \\mu_1 = \\mu_2 = \\dots = \\mu_k$$\n",
    "\n",
    "W tym zapisie $\\mu_1, \\dots, \\mu_k$ to średnie w populacji dla poszczególnych poziomów czynnika grupującego.\n",
    "\n",
    "$$H_A: \\neg H_0$$\n",
    "\n",
    "Hipoteza alternatywna mówi więc: *nie wszystkie średnie są równe* (co najmniej jedna różni się od pozostałych).\n",
    "\n",
    "Warto zauważyć, że $H_0$ jest sformułowana bardzo „ściśle” i może być fałszywa na wiele sposobów: wszystkie średnie mogą być różne, albo tylko niektóre.\n",
    "\n",
    "Sam test ANOVA jest **testem globalnym (omnibus)**: jeżeli odrzucimy $H_0$, to wiemy, że *gdzieś* są różnice, ale test nie mówi jeszcze *między którymi grupami*.\n",
    "\n",
    "Żeby dowiedzieć się, w jaki sposób $H_0$ jest fałszywa, wykonujemy analizę **post-hoc** (np. test Tukeya). Uwaga: sensownie robimy to dopiero wtedy, gdy wynik ANOVA jest istotny statystycznie.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0adff711",
   "metadata": {},
   "source": [
    "## Wariancja międzygrupowa i wewnątrzgrupowa\n",
    "\n",
    "Logika analizy wariancji polega na porównywaniu wariancji **międzygrupowej** i **wewnątrzgrupowej**.\n",
    "\n",
    "Korzystamy z faktu, że — przy założeniu hipotezy zerowej — na podstawie obu możemy skonstruować estymator wariancji populacji. Jeśli $H_0$ jest prawdziwa, oba estymatory „mierzą to samo” i ich iloraz powinien być bliski 1.\n",
    "\n",
    "Jeśli $H_0$ jest fałszywa, jeden z estymatorów przestaje być sensownym estymatorem tej samej wariancji i iloraz robi się większy niż 1.\n",
    "\n",
    "Mechanizm stanie się jaśniejszy na przykładzie.\n",
    "\n",
    "### Przykład: 5 grup, ale bez efektu (sytuacja z $H_0$)\n",
    "\n",
    "Stwórzmy dane, w których czynnik ma 5 poziomów (tu: „bardzo mało”, „mało”, „średnio”, „dużo”, „bardzo dużo”), ale **we wszystkich grupach losujemy z tej samej populacji** (ta sama średnia i to samo odchylenie standardowe). To jest sytuacja, w której $H_0$ jest prawdziwa.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c63979d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "levels = [\"bardzo mało\", \"mało\", \"średnio\", \"dużo\", \"bardzo dużo\"]\n",
    "\n",
    "n_per_group = 200\n",
    "k = len(levels)\n",
    "\n",
    "condition = np.repeat(levels, n_per_group)\n",
    "value = rng.normal(loc=10.0, scale=5.0, size=k * n_per_group)\n",
    "\n",
    "df = pd.DataFrame({\"condition\": condition, \"value\": value})\n",
    "\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "270cecee",
   "metadata": {},
   "source": [
    "### Wariancja wewnątrz grup\n",
    "\n",
    "Do obliczeń grupowych użyjemy `groupby`:\n",
    "\n",
    "- najpierw grupujemy dane po `condition`,\n",
    "- potem w każdej grupie liczymy wariancję,\n",
    "- na końcu (w tym przykładzie) bierzemy średnią z wariancji grupowych.\n",
    "\n",
    "To da nam pierwszy estymator wariancji populacji (intuicyjnie: „typowa” wariancja w grupie).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "613f0d9d",
   "metadata": {},
   "outputs": [],
   "source": [
    "group_variances = df.groupby(\"condition\")[\"value\"].var(ddof=1)\n",
    "group_variances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5643b88",
   "metadata": {},
   "outputs": [],
   "source": [
    "mse = float(group_variances.mean())\n",
    "mse"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ac5c8f5",
   "metadata": {},
   "source": [
    "Częściej spotkamy się ze sposobem obliczania, w którym najpierw liczymy **sumę kwadratów**, a następnie **średni kwadrat**.\n",
    "\n",
    "$$\n",
    "SSE  = \\sum_{i=1}^{k} \\sum_{j=1}^{n_i}(X_{ij} - \\bar{X_i})^2\n",
    "$$\n",
    "\n",
    "$$\n",
    "MSE = \\frac{SSE}{n-k}\n",
    "$$\n",
    "\n",
    "Ten sposób jest równoważny temu, co zrobiliśmy wcześniej (przy równych licznościach w grupach).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d621c49c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Liczebności i liczba grup\n",
    "n_total = len(df)\n",
    "k = df[\"condition\"].nunique()\n",
    "\n",
    "# SSE: suma kwadratów wewnątrz grup\n",
    "sse_terms = []\n",
    "for level_name, group_df in df.groupby(\"condition\"):\n",
    "    group_values = group_df[\"value\"].to_numpy()\n",
    "    group_mean = float(group_values.mean())\n",
    "    deviations = group_values - group_mean\n",
    "    sse_terms.append(float(np.sum(deviations**2)))\n",
    "\n",
    "sse = float(np.sum(sse_terms))\n",
    "\n",
    "sse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d806fca",
   "metadata": {},
   "outputs": [],
   "source": [
    "mse_from_sse = sse / (n_total - k)\n",
    "mse_from_sse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d0fc6a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# (Sanity check) Dwie wersje MSE powinny być prawie identyczne\n",
    "mse, mse_from_sse, abs(mse - mse_from_sse)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5959e7fb",
   "metadata": {},
   "source": [
    "### Wariancja między grupami\n",
    "\n",
    "Drugi estymator wariancji w populacji skonstruujemy korzystając z wiedzy o wariancji średnich z próby.\n",
    "\n",
    "Z Centralnego Twierdzenia Granicznego (dla średniej) wiemy, że:\n",
    "\n",
    "$$\n",
    "\\frac{\\sigma^2}{n} = s^2_{\\bar{X}}\n",
    "$$\n",
    "\n",
    "czyli równoważnie:\n",
    "\n",
    "$$\n",
    "\\sigma^2 = n s^2_{\\bar{X}}\n",
    "$$\n",
    "\n",
    "Widzimy więc, że wariancja w populacji ($\\sigma^2$) to wariancja średnich z próby ($s^2_{\\bar{X}}$) pomnożona przez liczebność próby.\n",
    "\n",
    "W takim razie, żeby wyestymować wariancję w populacji, możemy:\n",
    "\n",
    "1) policzyć średnią w każdej grupie, \n",
    "2) policzyć wariancję tych średnich, \n",
    "3) pomnożyć ją przez $n$ (liczebność w grupie).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89688ce9",
   "metadata": {},
   "outputs": [],
   "source": [
    "group_means = df.groupby(\"condition\")[\"value\"].mean()\n",
    "\n",
    "group_means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68c26716",
   "metadata": {},
   "outputs": [],
   "source": [
    "variance_of_means = float(group_means.var(ddof=1))\n",
    "variance_of_means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6250fa6",
   "metadata": {},
   "outputs": [],
   "source": [
    "mstreat = variance_of_means * n_per_group\n",
    "mstreat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2403dce",
   "metadata": {},
   "source": [
    "Częściej spotkamy się ze sposobem, w którym najpierw liczymy sumę kwadratów, a potem średni kwadrat:\n",
    "\n",
    "$$\n",
    "SSA = \\sum_{i=1}^{k} (\\bar{X_i} - \\bar{X})^2 n_i\n",
    "$$\n",
    "\n",
    "$$\n",
    "MSA = \\frac{SSA}{k-1}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2ec9884",
   "metadata": {},
   "outputs": [],
   "source": [
    "overall_mean = float(df[\"value\"].mean())\n",
    "\n",
    "ssa_terms = []\n",
    "for mean_value in group_means.to_list():\n",
    "    ssa_terms.append((float(mean_value) - overall_mean) ** 2 * n_per_group)\n",
    "\n",
    "ssa = float(np.sum(ssa_terms))\n",
    "msa = ssa / (k - 1)\n",
    "\n",
    "ssa, msa"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f718b15b",
   "metadata": {},
   "source": [
    "Zwróćmy uwagę: obliczona wartość będzie dobrym estymatorem wariancji w populacji **tylko wtedy**, gdy próby pochodzą z tej samej populacji.\n",
    "\n",
    "Jeżeli w populacji średnie rzeczywiście różnią się między poziomami czynnika, to „wariancja między grupami” będzie zawierać dodatkową składową związaną z efektem czynnika — i właśnie ten fakt wykorzystamy do skonstruowania testu statystycznego.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "febe2885",
   "metadata": {},
   "source": [
    "## Statystyka testowa\n",
    "\n",
    "Statystyka testowa $F$ będzie ilorazem dwóch obliczonych wyżej wartości:\n",
    "\n",
    "$$\n",
    "F = \\frac{MSA}{MSE}\n",
    "$$\n",
    "\n",
    "W przypadku, kiedy grupa nie ma znaczenia (czyli w uproszczeniu: $H_0$ jest prawdziwa), spodziewamy się, że statystyka ta będzie bliska 1 — bo licznik i mianownik są estymatorami tego samego parametru.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2810bee",
   "metadata": {},
   "outputs": [],
   "source": [
    "f_stat = msa / mse\n",
    "f_stat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2670cea1",
   "metadata": {},
   "source": [
    "### Widżet: skąd bierze się statystyka $F$?\n",
    "\n",
    "Poniższy widżet pokazuje ten sam mechanizm, który przed chwilą policzyliśmy ręcznie.\n",
    "\n",
    "Manipulujemy trzema rzeczami:\n",
    "\n",
    "- odstępem między średnimi grup,\n",
    "- odchyleniem standardowym wewnątrz grup,\n",
    "- liczebnością w każdej grupie.\n",
    "\n",
    "Oś pionowa jest stała. Dzięki temu widać, czy zmienia się faktyczny rozrzut danych, a nie tylko automatyczne skalowanie wykresu."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bab3a5cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def formatuj_p(p_value: float) -> str:\n",
    "    \"\"\"Zapisz p-wartość w krótkiej formie do opisu wykresu.\"\"\"\n",
    "\n",
    "    if p_value < 0.001:\n",
    "        return \"p < .001\"\n",
    "    return f\"p = {p_value:.3f}\"\n",
    "\n",
    "\n",
    "def anova_components_equal_n(values_by_group: np.ndarray) -> dict[str, float]:\n",
    "    \"\"\"Policz składniki jednoczynnikowej ANOVA dla równych liczebności grup.\"\"\"\n",
    "\n",
    "    n_groups, n_per_group_widget = values_by_group.shape\n",
    "    group_means_widget = values_by_group.mean(axis=1)\n",
    "    overall_mean_widget = float(values_by_group.mean())\n",
    "\n",
    "    sse_widget = float(np.sum((values_by_group - group_means_widget[:, None]) ** 2))\n",
    "    ssa_widget = float(\n",
    "        n_per_group_widget * np.sum((group_means_widget - overall_mean_widget) ** 2)\n",
    "    )\n",
    "\n",
    "    df_between_widget = n_groups - 1\n",
    "    df_within_widget = n_groups * (n_per_group_widget - 1)\n",
    "\n",
    "    msa_widget = ssa_widget / df_between_widget\n",
    "    mse_widget = sse_widget / df_within_widget\n",
    "    f_widget = msa_widget / mse_widget\n",
    "    p_widget = float(stats.f.sf(f_widget, dfn=df_between_widget, dfd=df_within_widget))\n",
    "    eta_squared_widget = ssa_widget / (ssa_widget + sse_widget)\n",
    "\n",
    "    return {\n",
    "        \"ssa\": ssa_widget,\n",
    "        \"sse\": sse_widget,\n",
    "        \"msa\": msa_widget,\n",
    "        \"mse\": mse_widget,\n",
    "        \"f\": f_widget,\n",
    "        \"p\": p_widget,\n",
    "        \"eta_squared\": eta_squared_widget,\n",
    "        \"df_between\": df_between_widget,\n",
    "        \"df_within\": df_within_widget,\n",
    "        \"overall_mean\": overall_mean_widget,\n",
    "    }\n",
    "\n",
    "\n",
    "anova_widget_rng = np.random.default_rng(SEED + 101)\n",
    "anova_widget_noise = anova_widget_rng.normal(size=(3, 60))\n",
    "anova_widget_jitter = anova_widget_rng.uniform(-0.08, 0.08, size=(3, 60))\n",
    "anova_widget_offsets = np.array([-1.0, 0.0, 1.0])\n",
    "anova_widget_labels = [\"grupa A\", \"grupa B\", \"grupa C\"]\n",
    "anova_widget_colors = sns.color_palette(\"colorblind\", n_colors=3)\n",
    "\n",
    "\n",
    "def rysuj_widget_f(odstep_srednich: float = 2.0, sigma: float = 5.0, n: int = 30) -> None:\n",
    "    group_values = 50 + odstep_srednich * anova_widget_offsets[:, None]\n",
    "    group_values = group_values + sigma * anova_widget_noise[:, :n]\n",
    "    components = anova_components_equal_n(group_values)\n",
    "\n",
    "    fig, (ax_data, ax_ratio) = plt.subplots(\n",
    "        ncols=2,\n",
    "        figsize=(11, 4.2),\n",
    "        gridspec_kw={\"width_ratios\": [1.35, 1.0]},\n",
    "    )\n",
    "\n",
    "    x_positions = np.arange(3)\n",
    "    for group_index, label in enumerate(anova_widget_labels):\n",
    "        x = x_positions[group_index] + anova_widget_jitter[group_index, :n]\n",
    "        y = group_values[group_index]\n",
    "        group_mean = float(y.mean())\n",
    "\n",
    "        ax_data.scatter(\n",
    "            x,\n",
    "            y,\n",
    "            s=34,\n",
    "            alpha=0.65,\n",
    "            color=anova_widget_colors[group_index],\n",
    "            edgecolor=\"white\",\n",
    "            linewidth=0.5,\n",
    "        )\n",
    "        ax_data.hlines(\n",
    "            group_mean,\n",
    "            group_index - 0.22,\n",
    "            group_index + 0.22,\n",
    "            color=\"black\",\n",
    "            lw=2.2,\n",
    "        )\n",
    "\n",
    "    ax_data.axhline(\n",
    "        components[\"overall_mean\"],\n",
    "        color=\"0.25\",\n",
    "        ls=\"--\",\n",
    "        lw=1.4,\n",
    "        label=\"średnia ogólna\",\n",
    "    )\n",
    "    ax_data.set_title(\"Dane grupowe: średnie vs rozrzut w grupach\")\n",
    "    ax_data.set_xlim(-0.5, 2.5)\n",
    "    ax_data.set_ylim(25, 75)\n",
    "    ax_data.set_xticks(x_positions)\n",
    "    ax_data.set_xticklabels(anova_widget_labels)\n",
    "    ax_data.set_ylabel(\"wynik\")\n",
    "    ax_data.legend(loc=\"upper left\")\n",
    "\n",
    "    bars = ax_ratio.barh(\n",
    "        [\"MSA\\n(między)\", \"MSE\\n(wewnątrz)\"],\n",
    "        [components[\"msa\"], components[\"mse\"]],\n",
    "        color=[\"#4477AA\", \"#CCBB44\"],\n",
    "        alpha=0.9,\n",
    "    )\n",
    "    ax_ratio.set_xlim(0, 1200)\n",
    "    ax_ratio.set_xlabel(\"średni kwadrat (stała skala)\")\n",
    "    ax_ratio.set_title(\"Licznik i mianownik statystyki F\")\n",
    "\n",
    "    for bar in bars:\n",
    "        width = bar.get_width()\n",
    "        ax_ratio.text(\n",
    "            min(width + 20, 1150),\n",
    "            bar.get_y() + bar.get_height() / 2,\n",
    "            f\"{width:.1f}\",\n",
    "            va=\"center\",\n",
    "            fontsize=10,\n",
    "        )\n",
    "\n",
    "    opis = (\n",
    "        f\"F({components['df_between']:.0f}, {components['df_within']:.0f}) = \"\n",
    "        f\"{components['f']:.2f}\\n\"\n",
    "        f\"{formatuj_p(components['p'])}\\n\"\n",
    "        f\"η² = {components['eta_squared']:.3f}\"\n",
    "    )\n",
    "    ax_ratio.text(\n",
    "        0.98,\n",
    "        0.05,\n",
    "        opis,\n",
    "        transform=ax_ratio.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_f = {\n",
    "        \"odstep_srednich\": widgets.FloatSlider(\n",
    "            value=2.0,\n",
    "            min=0.0,\n",
    "            max=4.0,\n",
    "            step=0.25,\n",
    "            description=\"odstęp średnich\",\n",
    "            continuous_update=False,\n",
    "            style=widget_style,\n",
    "        ),\n",
    "        \"sigma\": widgets.FloatSlider(\n",
    "            value=5.0,\n",
    "            min=2.0,\n",
    "            max=8.0,\n",
    "            step=0.25,\n",
    "            description=\"σ w grupach\",\n",
    "            continuous_update=False,\n",
    "            style=widget_style,\n",
    "        ),\n",
    "        \"n\": widgets.IntSlider(\n",
    "            value=30,\n",
    "            min=10,\n",
    "            max=60,\n",
    "            step=5,\n",
    "            description=\"n w grupie\",\n",
    "            continuous_update=False,\n",
    "            style=widget_style,\n",
    "        ),\n",
    "    }\n",
    "\n",
    "    output_f = widgets.interactive_output(rysuj_widget_f, controls_f)\n",
    "    display(widgets.VBox([widgets.HBox(list(controls_f.values())), output_f]))\n",
    "else:\n",
    "    rysuj_widget_f()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46ebc781",
   "metadata": {},
   "source": [
    "Czy ta wartość jest „duża” czy „mała”? Musimy odnieść ją do rozkładu statystyki pod $H_0$.\n",
    "\n",
    "W ANOVA (przy założeniach modelu) statystyka ma rozkład **$F$** z parametrami:\n",
    "\n",
    "- $df_1 = k - 1$ (stopnie swobody w liczniku),\n",
    "- $df_2 = n - k$ (stopnie swobody w mianowniku).\n",
    "\n",
    "Warto zobaczyć to również symulacyjnie.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c7a087f",
   "metadata": {},
   "outputs": [],
   "source": [
    "df1 = k - 1\n",
    "df2 = n_total - k\n",
    "\n",
    "df1, df2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4aad0f69",
   "metadata": {},
   "source": [
    "### Symulacja rozkładu $F$ pod $H_0$\n",
    "\n",
    "Powtórzymy nasz „eksperyment” 10 000 razy, a w każdym powtórzeniu:\n",
    "\n",
    "- losujemy dane w 5 grupach z tej samej populacji,\n",
    "- liczymy $F$,\n",
    "- a na końcu rysujemy histogram i porównujemy go z gęstością teoretycznego rozkładu $F(df_1, df_2)$.\n",
    "\n",
    "Żeby kod w pętli był krótki, zrobimy małą funkcję pomocniczą liczącą $F$ dla danych w postaci macierzy (wiersz = grupa, kolumny = obserwacje).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8bfad7a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_statistic_oneway(values_by_group: np.ndarray) -> float:\n",
    "    '''Policz statystykę F dla jednoczynnikowej ANOVA.\n",
    "\n",
    "    Parametry\n",
    "    ---------\n",
    "    values_by_group:\n",
    "        Tablica o kształcie (k, n_i): wiersze to grupy, kolumny to obserwacje.\n",
    "\n",
    "    Zwraca\n",
    "    ------\n",
    "    float\n",
    "        Wartość statystyki F = MSA / MSE.\n",
    "    '''\n",
    "\n",
    "    n_groups, n_per_group = values_by_group.shape\n",
    "\n",
    "    group_means = values_by_group.mean(axis=1)\n",
    "    overall_mean = values_by_group.mean()\n",
    "\n",
    "    # SSE i MSE (wewnątrz grup)\n",
    "    sse = np.sum((values_by_group - group_means[:, None]) ** 2)\n",
    "    mse = sse / (n_groups * n_per_group - n_groups)\n",
    "\n",
    "    # SSA i MSA (między grupami)\n",
    "    ssa = np.sum((group_means - overall_mean) ** 2) * n_per_group\n",
    "    msa = ssa / (n_groups - 1)\n",
    "\n",
    "    return float(msa / mse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ac3e33a",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_rep = 10_000\n",
    "\n",
    "f_values = np.empty(n_rep)\n",
    "\n",
    "for i in range(n_rep):\n",
    "    sample = rng.normal(loc=10.0, scale=5.0, size=k * n_per_group)\n",
    "    sample = sample.reshape(k, n_per_group)\n",
    "    f_values[i] = f_statistic_oneway(sample)\n",
    "\n",
    "f_values[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4726aea2",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "sns.histplot(f_values, stat=\"density\", bins=40, ax=ax)\n",
    "\n",
    "x = np.linspace(0, np.quantile(f_values, 0.995), 400)\n",
    "pdf = stats.f.pdf(x, dfn=df1, dfd=df2)\n",
    "ax.plot(x, pdf, color=\"black\", lw=2, label=f\"F(df1={df1}, df2={df2})\")\n",
    "\n",
    "ax.set_title(\"Rozkład statystyki F pod H0 (symulacja vs teoria)\")\n",
    "ax.set_xlabel(\"F\")\n",
    "ax.set_ylabel(\"Gęstość\")\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae713191",
   "metadata": {},
   "source": [
    "### Region krytyczny i p-wartość (pierwsza ANOVA krok-po-kroku)\n",
    "\n",
    "Teraz przeprowadzimy test w klasycznym schemacie:\n",
    "\n",
    "- **Historia/pytanie:** czy średnie 5 grup różnią się? (tu wiemy, że dane generowaliśmy tak, aby $H_0$ była prawdziwa).\n",
    "- **Hipotezy:**\n",
    "  - $H_0: \\mu_1 = \\dots = \\mu_5$,\n",
    "  - $H_A: \\neg H_0$.\n",
    "- **Poziom istotności:** przyjmijmy $\\alpha = 0.05$.\n",
    "- **Statystyka pod $H_0$:** $F \\sim F(df_1=k-1, df_2=n-k)$.\n",
    "- **Region krytyczny:** odrzucamy $H_0$, gdy $F \\ge F_{1-\\alpha}(df_1, df_2)$ (prawy ogon rozkładu).\n",
    "- **p-wartość „na piechotę”:** $p = P(F \\ge f_{obs} \\mid H_0)$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be05b6d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 0.05\n",
    "\n",
    "f_crit = stats.f.ppf(1 - alpha, dfn=df1, dfd=df2)\n",
    "p_value = stats.f.sf(f_stat, dfn=df1, dfd=df2)\n",
    "\n",
    "f_stat, f_crit, p_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6da9f19",
   "metadata": {},
   "outputs": [],
   "source": [
    "reject_h0 = f_stat >= f_crit\n",
    "reject_h0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78316705",
   "metadata": {},
   "source": [
    "Na podstawie p-wartości podejmujemy decyzję:\n",
    "\n",
    "- jeśli $p \\le \\alpha$, odrzucamy $H_0$,\n",
    "- jeśli $p > \\alpha$, nie odrzucamy $H_0$.\n",
    "\n",
    "W tym przykładzie spodziewamy się typowo $p > 0.05$, bo dane generowaliśmy bez efektu czynnika.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c4d0749",
   "metadata": {},
   "source": [
    "### Weryfikacja funkcją z biblioteki\n",
    "\n",
    "Zamiast ręcznie liczyć składniki ANOVA, możemy dopasować model liniowy z predyktorem kategorialnym i poprosić bibliotekę o tabelę ANOVA.\n",
    "\n",
    "W Pythonie najczęściej robi się to przez `statsmodels`:\n",
    "\n",
    "- dopasowujemy OLS: `value ~ C(condition)`,\n",
    "- a następnie liczymy tabelę ANOVA.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4719a211",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = smf.ols(\"value ~ C(condition)\", data=df).fit()\n",
    "anova_table = anova_lm(model, typ=1)\n",
    "\n",
    "anova_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb78106d",
   "metadata": {},
   "outputs": [],
   "source": [
    "f_lib = float(anova_table.loc[\"C(condition)\", \"F\"])\n",
    "p_lib = float(anova_table.loc[\"C(condition)\", \"PR(>F)\"])\n",
    "\n",
    "f_stat, f_lib, p_value, p_lib"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ddfa9ff",
   "metadata": {},
   "source": [
    "Widzimy, że wynik „na piechotę” i wynik z funkcji bibliotecznej są zgodne (różnice, jeśli występują, wynikają tylko z zaokrągleń).\n",
    "\n",
    "To była nasza pierwsza analiza wariancji: policzyliśmy statystykę $F$, p-wartość i porównaliśmy wynik z tym, co daje biblioteka.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a12cd624",
   "metadata": {},
   "source": [
    "## Analiza wariancji na przykładzie\n",
    "\n",
    "Prześledźmy analizę wariancji na przykładzie.\n",
    "\n",
    "Wykorzystamy niewielki zbiór danych `PlantGrowth`, dotyczący wpływu pewnych warunków na wagę roślin.\n",
    "\n",
    "- w jednej kolumnie mamy wagę roślin,\n",
    "- w drugiej – grupę (czynnik): dwie eksperymentalne (`trt1`, `trt2`) oraz kontrolna (`ctrl`).\n",
    "\n",
    "Za pomocą ANOVA badamy, czy warunki eksperymentalne miały wpływ na wagę roślin.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2adf86ab",
   "metadata": {},
   "source": [
    "### Historia → pytanie → hipotezy\n",
    "\n",
    "- **Historia:** mamy trzy grupy roślin (kontrola i dwa warianty „zabiegu”).\n",
    "- **Pytanie badawcze:** czy średnia waga roślin różni się między grupami?\n",
    "- **Hipotezy:**\n",
    "\n",
    "$$H_0: \\mu_{ctrl} = \\mu_{trt1} = \\mu_{trt2}$$\n",
    "\n",
    "$$H_A: \\neg H_0$$\n",
    "\n",
    "To jest hipoteza alternatywna „nieskierowana” (test globalny): nie wskazujemy, *która* grupa jest większa, tylko pytamy, czy **jakiekolwiek** średnie różnią się między sobą.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "804ca1f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "plantgrowth = pd.DataFrame(\n",
    "    {\n",
    "        \"weight\": [\n",
    "            4.17, 5.58, 5.18, 6.11, 4.50, 4.61, 5.17, 4.53, 5.33, 5.14,\n",
    "            4.81, 4.17, 4.41, 3.59, 5.87, 3.83, 6.03, 4.89, 4.32, 4.69,\n",
    "            6.31, 5.12, 5.54, 5.50, 5.37, 5.29, 4.92, 6.15, 5.80, 5.26,\n",
    "        ],\n",
    "        \"group\": [\n",
    "            \"ctrl\", \"ctrl\", \"ctrl\", \"ctrl\", \"ctrl\", \"ctrl\", \"ctrl\", \"ctrl\", \"ctrl\", \"ctrl\",\n",
    "            \"trt1\", \"trt1\", \"trt1\", \"trt1\", \"trt1\", \"trt1\", \"trt1\", \"trt1\", \"trt1\", \"trt1\",\n",
    "            \"trt2\", \"trt2\", \"trt2\", \"trt2\", \"trt2\", \"trt2\", \"trt2\", \"trt2\", \"trt2\", \"trt2\",\n",
    "        ],\n",
    "    }\n",
    ")\n",
    "\n",
    "plantgrowth"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9f1e304",
   "metadata": {},
   "source": [
    "### Szybka estymacja: średnie i przedziały ufności\n",
    "\n",
    "Zanim przejdziemy do testu, policzmy:\n",
    "\n",
    "- liczebności,\n",
    "- średnie,\n",
    "- odchylenia standardowe,\n",
    "- 95% przedziały ufności dla średniej w każdej grupie.\n",
    "\n",
    "To jest podejście „estimation-first”: najpierw pytamy, **jak duży** jest efekt i jak niepewny.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6967c11",
   "metadata": {},
   "outputs": [],
   "source": [
    "summary_rows = []\n",
    "\n",
    "for group_name, group_df in plantgrowth.groupby(\"group\"):\n",
    "    values = group_df[\"weight\"].to_numpy()\n",
    "    n = len(values)\n",
    "    mean = float(values.mean())\n",
    "    sd = float(values.std(ddof=1))\n",
    "    se = sd / np.sqrt(n)\n",
    "\n",
    "    t_crit = stats.t.ppf(0.975, df=n - 1)\n",
    "    ci_low = mean - t_crit * se\n",
    "    ci_high = mean + t_crit * se\n",
    "\n",
    "    summary_rows.append(\n",
    "        {\n",
    "            \"group\": group_name,\n",
    "            \"n\": n,\n",
    "            \"mean\": mean,\n",
    "            \"sd\": sd,\n",
    "            \"ci95_low\": ci_low,\n",
    "            \"ci95_high\": ci_high,\n",
    "        }\n",
    "    )\n",
    "\n",
    "summary_table = pd.DataFrame(summary_rows).sort_values(\"group\")\n",
    "summary_table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1cb84c0",
   "metadata": {},
   "source": [
    "### Wizualizacja danych i założeń\n",
    "\n",
    "Na początek warto stworzyć wizualizację danych — m.in. po to, aby wstępnie ocenić założenia ANOVA:\n",
    "\n",
    "- **homogeniczność wariancji** (czy rozrzut w grupach wygląda podobnie),\n",
    "- **normalność w grupach / normalność reszt** (czy rozkłady w grupach nie są wyraźnie „dziwne”),\n",
    "- **outliery**.\n",
    "\n",
    "Na wykresie pudełkowym możemy „na oko” ocenić, czy dane nie wyglądają na skrajnie problematyczne.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7206a0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "sns.boxplot(data=plantgrowth, x=\"group\", y=\"weight\", ax=ax)\n",
    "sns.stripplot(\n",
    "    data=plantgrowth,\n",
    "    x=\"group\",\n",
    "    y=\"weight\",\n",
    "    color=\"black\",\n",
    "    alpha=0.6,\n",
    "    jitter=True,\n",
    "    ax=ax,\n",
    ")\n",
    "\n",
    "ax.set_title(\"PlantGrowth: rozkład wagi roślin w grupach\")\n",
    "ax.set_xlabel(\"Grupa\")\n",
    "ax.set_ylabel(\"Waga\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65c0ecf5",
   "metadata": {},
   "source": [
    "Normalność w grupach można „na oko” ocenić za pomocą wykresów Q–Q dla każdej grupy z osobna.\n",
    "\n",
    "Przy małych próbach może to być bardziej pomocne niż formalny test normalności.\n",
    "Nie chcemy, żeby punkty w sposób wyraźny i systematyczny odbiegały od przekątnej.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "358d891c",
   "metadata": {},
   "outputs": [],
   "source": [
    "groups = [\"ctrl\", \"trt1\", \"trt2\"]\n",
    "\n",
    "fig, axes = plt.subplots(nrows=len(groups), ncols=1, figsize=(6, 10), dpi=120)\n",
    "\n",
    "for ax, g in zip(axes, groups, strict=True):\n",
    "    values = plantgrowth.loc[plantgrowth[\"group\"] == g, \"weight\"].to_numpy()\n",
    "    stats.probplot(values, dist=\"norm\", plot=ax)\n",
    "    ax.set_title(f\"Q–Q plot: grupa {g}\")\n",
    "    ax.set_xlabel(\"Kwantyle teoretyczne\")\n",
    "    ax.set_ylabel(\"Kwantyle empiryczne\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1261bcfc",
   "metadata": {},
   "source": [
    "Na tym etapie często robi się też szybki, formalny test równości wariancji (np. Levene’a) – ale pamiętajmy:\n",
    "\n",
    "- testy na założenia też mają p-wartości i też potrafią „wariować” przy małych próbach,\n",
    "- wykresy i sensowna wiedza domenowa zwykle są równie ważne.\n",
    "\n",
    "Poniżej tylko jako pomocnicza informacja:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efd4b3b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "levene_stat, levene_p = stats.levene(\n",
    "    plantgrowth.loc[plantgrowth[\"group\"] == \"ctrl\", \"weight\"],\n",
    "    plantgrowth.loc[plantgrowth[\"group\"] == \"trt1\", \"weight\"],\n",
    "    plantgrowth.loc[plantgrowth[\"group\"] == \"trt2\", \"weight\"],\n",
    ")\n",
    "\n",
    "float(levene_stat), float(levene_p)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "560c3597",
   "metadata": {},
   "source": [
    "### ANOVA „na piechotę” (krok po kroku)\n",
    "\n",
    "Teraz przeprowadzimy analizę wariancji ręcznie.\n",
    "\n",
    "1. Ustalamy $\\alpha$.\n",
    "2. Liczymy składniki: $SSE$, $MSE$, $SSA$, $MSA$.\n",
    "3. Liczymy statystykę $F = MSA/MSE$.\n",
    "4. Pod $H_0$: $F \\sim F(df_1, df_2)$.\n",
    "5. Wyznaczamy region krytyczny i p-wartość.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a13de7a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 0.05\n",
    "\n",
    "k = plantgrowth[\"group\"].nunique()\n",
    "n_total = len(plantgrowth)\n",
    "\n",
    "df1 = k - 1\n",
    "df2 = n_total - k\n",
    "\n",
    "df1, df2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20387cf7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# SSE: suma kwadratów odchyleń od średnich grupowych\n",
    "sse = 0.0\n",
    "\n",
    "for group_name, group_df in plantgrowth.groupby(\"group\"):\n",
    "    values = group_df[\"weight\"].to_numpy()\n",
    "    group_mean = float(values.mean())\n",
    "    deviations = values - group_mean\n",
    "    sse += float(np.sum(deviations**2))\n",
    "\n",
    "sse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1af6b983",
   "metadata": {},
   "outputs": [],
   "source": [
    "mse = sse / df2\n",
    "mse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "212d5678",
   "metadata": {},
   "outputs": [],
   "source": [
    "# SSA: międzygrupowa suma kwadratów\n",
    "overall_mean = float(plantgrowth[\"weight\"].mean())\n",
    "\n",
    "# W PlantGrowth w każdej grupie jest po 10 obserwacji.\n",
    "n_per_group = int(plantgrowth.groupby(\"group\").size().iloc[0])\n",
    "\n",
    "means = plantgrowth.groupby(\"group\")[\"weight\"].mean()\n",
    "ssa = float(np.sum(((means - overall_mean) ** 2) * n_per_group))\n",
    "\n",
    "ssa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa772b38",
   "metadata": {},
   "outputs": [],
   "source": [
    "msa = ssa / df1\n",
    "msa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54fe652e",
   "metadata": {},
   "outputs": [],
   "source": [
    "f_stat = msa / mse\n",
    "f_stat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0baeea1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Region krytyczny (prawy ogon rozkładu F)\n",
    "f_crit = stats.f.ppf(1 - alpha, dfn=df1, dfd=df2)\n",
    "\n",
    "f_crit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9222c9cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# p-wartość „na piechotę”: P(F >= f_obs | H0)\n",
    "p_value = stats.f.sf(f_stat, dfn=df1, dfd=df2)\n",
    "\n",
    "float(p_value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6526f10f",
   "metadata": {},
   "outputs": [],
   "source": [
    "reject_h0 = f_stat >= f_crit\n",
    "reject_h0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c709056",
   "metadata": {},
   "source": [
    "**Interpretacja:**\n",
    "\n",
    "- Otrzymana p-wartość mówi, jak bardzo „zaskakująca” byłaby statystyka $F$ równa co najmniej $f_{obs}$, **gdyby** wszystkie średnie w populacji były równe.\n",
    "- Ponieważ $p < 0.05$, odrzucamy $H_0$ na poziomie istotności $\\alpha = 0.05$.\n",
    "- To jeszcze nie mówi, *które* grupy różnią się między sobą — do tego wrócimy w analizie post-hoc.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a48c9576",
   "metadata": {},
   "source": [
    "### Weryfikacja funkcją z biblioteki\n",
    "\n",
    "Teraz wykonamy to samo przy pomocy biblioteki:\n",
    "\n",
    "- dopasujemy model `weight ~ C(group)`,\n",
    "- policzymy tabelę ANOVA.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1b914d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = smf.ols(\"weight ~ C(group)\", data=plantgrowth).fit()\n",
    "anova_table = anova_lm(model, typ=1)\n",
    "\n",
    "anova_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d44d3cd3",
   "metadata": {},
   "outputs": [],
   "source": [
    "f_lib = float(anova_table.loc[\"C(group)\", \"F\"])\n",
    "p_lib = float(anova_table.loc[\"C(group)\", \"PR(>F)\"])\n",
    "\n",
    "f_stat, f_lib, float(p_value), p_lib"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b13dddce",
   "metadata": {},
   "source": [
    "Widzimy zgodność: ręczne obliczenia i `statsmodels` dają to samo $F$ i tę samą p-wartość.\n",
    "\n",
    "Możemy więc stwierdzić, że dane dostarczają (przy założeniach ANOVA) podstaw do odrzucenia hipotezy, że wszystkie trzy grupy mają identyczną średnią w populacji.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4591b785",
   "metadata": {},
   "source": [
    "## Testy post-hoc\n",
    "\n",
    "Często interesuje nas nie tylko fakt, że odrzuciliśmy hipotezę „wszystkie średnie są równe”, ale również **między którymi grupami** występują statystycznie istotne różnice.\n",
    "\n",
    "W tym miejscu przydatne są testy *post-hoc*. Ich zaletą jest to, że uwzględniają wielokrotne porównania (w przeciwieństwie do „serii testów t” bez korekty).\n",
    "\n",
    "Jednym z najpopularniejszych testów *post-hoc* jest test Tukeya HSD (*honest significant difference*).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ef8ecc6",
   "metadata": {},
   "outputs": [],
   "source": [
    "tukey = pairwise_tukeyhsd(\n",
    "    endog=plantgrowth[\"weight\"],\n",
    "    groups=plantgrowth[\"group\"],\n",
    "    alpha=alpha,\n",
    ")\n",
    "\n",
    "tukey.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a084d95",
   "metadata": {},
   "source": [
    "Wydruk informuje nas o:\n",
    "\n",
    "- różnicy średnich (*meandiff*),\n",
    "- skorygowanej p-wartości (*p-adj*),\n",
    "- 95% przedziale ufności dla różnicy średnich (*lower*, *upper*),\n",
    "- decyzji (*reject*).\n",
    "\n",
    "Istotne statystycznie są te różnice, których przedziały ufności **nie obejmują zera**.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d832202",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = tukey.plot_simultaneous(\n",
    "    figsize=(8, 3),\n",
    "    xlabel=\"Średnia grupy (jednoczesny 95% CI)\",\n",
    ")\n",
    "plt.title(\"PlantGrowth: jednoczesne przedziały ufności dla średnich grup\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91ab737f",
   "metadata": {},
   "source": [
    "Powyższy wykres nie pokazuje bezpośrednio kontrastów. Pokazuje jednoczesne przedziały ufności dla **średnich grup**.\n",
    "\n",
    "Jeżeli chcemy wykres bezpośrednio pokazujący kontrasty Tukeya, bardziej naturalny jest wykres poniżej: każda linia to jedna **różnica średnich** wraz z 95% przedziałem ufności po korekcie Tukeya."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d375784",
   "metadata": {},
   "outputs": [],
   "source": [
    "pair_labels = []\n",
    "for first_index in range(len(tukey.groupsunique)):\n",
    "    for second_index in range(first_index + 1, len(tukey.groupsunique)):\n",
    "        first_group = tukey.groupsunique[first_index]\n",
    "        second_group = tukey.groupsunique[second_index]\n",
    "        pair_labels.append(f\"{second_group} - {first_group}\")\n",
    "\n",
    "y_positions = np.arange(len(pair_labels))\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 3.4))\n",
    "\n",
    "for comparison_index, y_position in enumerate(y_positions):\n",
    "    ci_low, ci_high = tukey.confint[comparison_index]\n",
    "    mean_diff = tukey.meandiffs[comparison_index]\n",
    "    is_rejected = bool(tukey.reject[comparison_index])\n",
    "    color = \"#228833\" if is_rejected else \"0.45\"\n",
    "\n",
    "    ax.plot([ci_low, ci_high], [y_position, y_position], color=color, lw=3)\n",
    "    ax.scatter(mean_diff, y_position, color=color, s=55, zorder=3)\n",
    "    ax.text(\n",
    "        1.65,\n",
    "        y_position,\n",
    "        f\"p = {tukey.pvalues[comparison_index]:.3f}\",\n",
    "        va=\"center\",\n",
    "        fontsize=9,\n",
    "    )\n",
    "\n",
    "ax.axvline(0, color=\"black\", lw=1.2, ls=\"--\")\n",
    "ax.set_xlim(-1.3, 2.1)\n",
    "ax.set_ylim(-0.6, len(pair_labels) - 0.4)\n",
    "ax.set_yticks(y_positions)\n",
    "ax.set_yticklabels(pair_labels)\n",
    "ax.set_xlabel(\"Różnica średnich (95% CI Tukeya)\")\n",
    "ax.set_title(\"PlantGrowth: kontrasty Tukeya\")\n",
    "\n",
    "ax.text(\n",
    "    0.03,\n",
    "    0.05,\n",
    "    \"zielone: CI nie obejmuje 0\\nszare: CI obejmuje 0\",\n",
    "    transform=ax.transAxes,\n",
    "    va=\"bottom\",\n",
    "    bbox={\"boxstyle\": \"round,pad=0.4\", \"fc\": \"white\", \"ec\": \"0.35\"},\n",
    ")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5809e3ad",
   "metadata": {},
   "source": [
    "### Widżet: ANOVA mówi „co najmniej jedna grupa”, Tukey pokazuje gdzie\n",
    "\n",
    "Globalny test ANOVA odpowiada na pytanie, czy wszystkie średnie można rozsądnie traktować jako równe. Nie mówi jednak, **które** grupy różnią się od których.\n",
    "\n",
    "W tym widżecie wybieramy prosty układ średnich i patrzymy jednocześnie na wynik ANOVA oraz porównania Tukeya. Skale osi są stałe, żeby było widać, kiedy różnice są małe, duże albo giną w szumie."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6452ee9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "posthoc_widget_rng = np.random.default_rng(SEED + 202)\n",
    "posthoc_widget_noise = posthoc_widget_rng.normal(size=(3, 45))\n",
    "posthoc_widget_jitter = posthoc_widget_rng.uniform(-0.08, 0.08, size=(3, 45))\n",
    "posthoc_widget_groups = np.array([\"A\", \"B\", \"C\"])\n",
    "posthoc_widget_labels = [\"A\", \"B\", \"C\"]\n",
    "posthoc_widget_colors = sns.color_palette(\"colorblind\", n_colors=3)\n",
    "\n",
    "posthoc_mean_patterns = {\n",
    "    \"wszystkie podobne\": np.array([50.0, 50.0, 50.0]),\n",
    "    \"jedna grupa odstaje\": np.array([50.0, 50.0, 56.0]),\n",
    "    \"gradient\": np.array([47.0, 50.0, 53.0]),\n",
    "    \"dwie podobne + jedna niższa\": np.array([46.0, 54.0, 54.0]),\n",
    "}\n",
    "\n",
    "\n",
    "def rysuj_widget_anova_tukey(\n",
    "    uklad_srednich: str = \"jedna grupa odstaje\",\n",
    "    sigma: float = 4.0,\n",
    "    n: int = 20,\n",
    ") -> None:\n",
    "    means = posthoc_mean_patterns[uklad_srednich]\n",
    "    values_by_group = means[:, None] + sigma * posthoc_widget_noise[:, :n]\n",
    "\n",
    "    values = values_by_group.ravel()\n",
    "    groups = np.repeat(posthoc_widget_groups, n)\n",
    "\n",
    "    f_stat_widget, p_value_widget = stats.f_oneway(\n",
    "        values_by_group[0],\n",
    "        values_by_group[1],\n",
    "        values_by_group[2],\n",
    "    )\n",
    "    tukey_widget = pairwise_tukeyhsd(endog=values, groups=groups, alpha=0.05)\n",
    "\n",
    "    fig, (ax_data, ax_tukey) = plt.subplots(\n",
    "        ncols=2,\n",
    "        figsize=(11, 4.2),\n",
    "        gridspec_kw={\"width_ratios\": [1.0, 1.15]},\n",
    "    )\n",
    "\n",
    "    x_positions = np.arange(3)\n",
    "    for group_index, label in enumerate(posthoc_widget_labels):\n",
    "        x = x_positions[group_index] + posthoc_widget_jitter[group_index, :n]\n",
    "        y = values_by_group[group_index]\n",
    "\n",
    "        ax_data.scatter(\n",
    "            x,\n",
    "            y,\n",
    "            s=34,\n",
    "            alpha=0.65,\n",
    "            color=posthoc_widget_colors[group_index],\n",
    "            edgecolor=\"white\",\n",
    "            linewidth=0.5,\n",
    "        )\n",
    "        ax_data.hlines(\n",
    "            float(y.mean()),\n",
    "            group_index - 0.22,\n",
    "            group_index + 0.22,\n",
    "            color=\"black\",\n",
    "            lw=2.2,\n",
    "        )\n",
    "\n",
    "    decyzja = \"odrzucamy H₀\" if p_value_widget <= 0.05 else \"nie odrzucamy H₀\"\n",
    "    opis_anova = f\"ANOVA: F = {f_stat_widget:.2f}\\n{formatuj_p(p_value_widget)}\\n{decyzja}\"\n",
    "    ax_data.text(\n",
    "        0.03,\n",
    "        0.97,\n",
    "        opis_anova,\n",
    "        transform=ax_data.transAxes,\n",
    "        va=\"top\",\n",
    "        bbox={\"boxstyle\": \"round,pad=0.45\", \"fc\": \"white\", \"ec\": \"0.35\"},\n",
    "    )\n",
    "    ax_data.set_title(\"Dane grupowe\")\n",
    "    ax_data.set_xlim(-0.5, 2.5)\n",
    "    ax_data.set_ylim(35, 65)\n",
    "    ax_data.set_xticks(x_positions)\n",
    "    ax_data.set_xticklabels(posthoc_widget_labels)\n",
    "    ax_data.set_ylabel(\"wynik\")\n",
    "\n",
    "    pair_labels = []\n",
    "    for first_index in range(len(tukey_widget.groupsunique)):\n",
    "        for second_index in range(first_index + 1, len(tukey_widget.groupsunique)):\n",
    "            first = tukey_widget.groupsunique[first_index]\n",
    "            second = tukey_widget.groupsunique[second_index]\n",
    "            pair_labels.append(f\"{second} - {first}\")\n",
    "\n",
    "    y_positions = np.arange(len(pair_labels))\n",
    "    for comparison_index, y_position in enumerate(y_positions):\n",
    "        ci_low, ci_high = tukey_widget.confint[comparison_index]\n",
    "        mean_diff = tukey_widget.meandiffs[comparison_index]\n",
    "        is_rejected = bool(tukey_widget.reject[comparison_index])\n",
    "        color = \"#228833\" if is_rejected else \"0.45\"\n",
    "\n",
    "        ax_tukey.plot([ci_low, ci_high], [y_position, y_position], color=color, lw=3)\n",
    "        ax_tukey.scatter(mean_diff, y_position, color=color, s=55, zorder=3)\n",
    "\n",
    "    ax_tukey.axvline(0, color=\"black\", lw=1.2, ls=\"--\")\n",
    "    ax_tukey.set_xlim(-15, 15)\n",
    "    ax_tukey.set_ylim(-0.6, len(pair_labels) - 0.4)\n",
    "    ax_tukey.set_yticks(y_positions)\n",
    "    ax_tukey.set_yticklabels(pair_labels)\n",
    "    ax_tukey.set_xlabel(\"różnica średnich (95% CI Tukeya)\")\n",
    "    ax_tukey.set_title(\"Porównania post-hoc\")\n",
    "\n",
    "    ax_tukey.text(\n",
    "        0.98,\n",
    "        0.05,\n",
    "        \"zielone: istotne po korekcie\\nszare: brak istotności\",\n",
    "        transform=ax_tukey.transAxes,\n",
    "        ha=\"right\",\n",
    "        va=\"bottom\",\n",
    "        bbox={\"boxstyle\": \"round,pad=0.4\", \"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_tukey = {\n",
    "        \"uklad_srednich\": widgets.Dropdown(\n",
    "            options=list(posthoc_mean_patterns.keys()),\n",
    "            value=\"jedna grupa odstaje\",\n",
    "            description=\"układ średnich\",\n",
    "            style=widget_style,\n",
    "        ),\n",
    "        \"sigma\": widgets.FloatSlider(\n",
    "            value=4.0,\n",
    "            min=2.0,\n",
    "            max=7.0,\n",
    "            step=0.25,\n",
    "            description=\"σ w grupach\",\n",
    "            continuous_update=False,\n",
    "            style=widget_style,\n",
    "        ),\n",
    "        \"n\": widgets.IntSlider(\n",
    "            value=20,\n",
    "            min=10,\n",
    "            max=45,\n",
    "            step=5,\n",
    "            description=\"n w grupie\",\n",
    "            continuous_update=False,\n",
    "            style=widget_style,\n",
    "        ),\n",
    "    }\n",
    "\n",
    "    output_tukey = widgets.interactive_output(rysuj_widget_anova_tukey, controls_tukey)\n",
    "    display(widgets.VBox([widgets.HBox(list(controls_tukey.values())), output_tukey]))\n",
    "else:\n",
    "    rysuj_widget_anova_tukey()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61831caf",
   "metadata": {},
   "source": [
    "## Wielkość efektu\n",
    "\n",
    "W przypadku analizy wariancji standardowo używa się miary $\\eta^2$.\n",
    "\n",
    "Definicja:\n",
    "\n",
    "$$\n",
    "\\eta^2 = \\frac{SSA}{SST}\n",
    "$$\n",
    "\n",
    "gdzie $SST$ to całkowita suma kwadratów.\n",
    "\n",
    "Interpretacja jest analogiczna jak w przypadku znanego już $R^2$: to część zmienności zmiennej zależnej wyjaśniana przez czynnik.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87b3f758",
   "metadata": {},
   "outputs": [],
   "source": [
    "# SST: całkowita suma kwadratów\n",
    "sst = float(np.sum((plantgrowth[\"weight\"] - overall_mean) ** 2))\n",
    "\n",
    "eta_squared = ssa / sst\n",
    "eta_squared"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e0855d8",
   "metadata": {},
   "source": [
    "### Widżet: liczebność, p-wartość, $\\eta^2$ i moc\n",
    "\n",
    "Ten widżet oddziela trzy rzeczy, które łatwo pomylić:\n",
    "\n",
    "- **wielkość efektu**: jak dużą część zmienności wyjaśnia czynnik,\n",
    "- **p-wartość**: jak nietypowy byłby taki wynik przy $H_0$,\n",
    "- **moc testu**: jak często test wykryłby taki efekt przy danej liczebności i zmienności.\n",
    "\n",
    "Najważniejsza obserwacja: zwiększanie `n` zwykle silnie zmienia p-wartość i moc, ale nie musi oznaczać większego efektu."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "239c8b28",
   "metadata": {},
   "outputs": [],
   "source": [
    "power_widget_rng = np.random.default_rng(SEED + 303)\n",
    "power_widget_noise = power_widget_rng.normal(size=(3, 80))\n",
    "power_widget_colors = sns.color_palette(\"colorblind\", n_colors=3)\n",
    "power_widget_x = np.linspace(30, 70, 500)\n",
    "\n",
    "\n",
    "def anova_power_balanced(means: np.ndarray, sigma: float, n: int, alpha: float = 0.05) -> dict[str, float]:\n",
    "    \"\"\"Policz teoretyczną moc jednoczynnikowej ANOVA dla zbalansowanego układu.\"\"\"\n",
    "\n",
    "    n_groups = len(means)\n",
    "    grand_mean = float(means.mean())\n",
    "    deviations = means - grand_mean\n",
    "\n",
    "    df_between_power = n_groups - 1\n",
    "    df_within_power = n_groups * (n - 1)\n",
    "    noncentrality = float(n * np.sum(deviations**2) / sigma**2)\n",
    "    f_critical = float(stats.f.ppf(1 - alpha, dfn=df_between_power, dfd=df_within_power))\n",
    "    power = float(stats.ncf.sf(f_critical, dfn=df_between_power, dfd=df_within_power, nc=noncentrality))\n",
    "\n",
    "    between_variance = float(np.mean(deviations**2))\n",
    "    eta_squared_population = between_variance / (between_variance + sigma**2)\n",
    "\n",
    "    return {\n",
    "        \"df_between\": df_between_power,\n",
    "        \"df_within\": df_within_power,\n",
    "        \"noncentrality\": noncentrality,\n",
    "        \"f_critical\": f_critical,\n",
    "        \"power\": power,\n",
    "        \"eta_squared_population\": eta_squared_population,\n",
    "    }\n",
    "\n",
    "\n",
    "def rysuj_widget_moc(odstep_srednich: float = 2.0, sigma: float = 5.0, n: int = 30) -> None:\n",
    "    means = np.array([50.0 - odstep_srednich, 50.0, 50.0 + odstep_srednich])\n",
    "    power_info = anova_power_balanced(means=means, sigma=sigma, n=n, alpha=0.05)\n",
    "\n",
    "    values_by_group = means[:, None] + sigma * power_widget_noise[:, :n]\n",
    "    sample_components = anova_components_equal_n(values_by_group)\n",
    "\n",
    "    fig, (ax_curves, ax_info) = plt.subplots(\n",
    "        ncols=2,\n",
    "        figsize=(11, 4.2),\n",
    "        gridspec_kw={\"width_ratios\": [1.25, 1.0]},\n",
    "    )\n",
    "\n",
    "    for group_index, mean_value in enumerate(means):\n",
    "        density = stats.norm.pdf(power_widget_x, loc=mean_value, scale=sigma)\n",
    "        ax_curves.plot(\n",
    "            power_widget_x,\n",
    "            density,\n",
    "            color=power_widget_colors[group_index],\n",
    "            lw=2.2,\n",
    "            label=f\"μ{group_index + 1} = {mean_value:.1f}\",\n",
    "        )\n",
    "        ax_curves.fill_between(\n",
    "            power_widget_x,\n",
    "            density,\n",
    "            color=power_widget_colors[group_index],\n",
    "            alpha=0.12,\n",
    "        )\n",
    "\n",
    "    ax_curves.set_title(\"Założone rozkłady w populacji\")\n",
    "    ax_curves.set_xlim(30, 70)\n",
    "    ax_curves.set_ylim(0, 0.15)\n",
    "    ax_curves.set_xlabel(\"wynik\")\n",
    "    ax_curves.set_ylabel(\"gęstość\")\n",
    "    ax_curves.legend(loc=\"upper left\")\n",
    "\n",
    "    opis = (\n",
    "        \"Parametry populacji\\n\"\n",
    "        f\"n w grupie = {n}\\n\"\n",
    "        f\"η²_pop = {power_info['eta_squared_population']:.3f}\\n\"\n",
    "        f\"moc przy α = .05: {power_info['power']:.3f}\\n\\n\"\n",
    "        \"Jedna przykładowa próba\\n\"\n",
    "        f\"F({sample_components['df_between']:.0f}, {sample_components['df_within']:.0f}) = \"\n",
    "        f\"{sample_components['f']:.2f}\\n\"\n",
    "        f\"{formatuj_p(sample_components['p'])}\\n\"\n",
    "        f\"η²_próby = {sample_components['eta_squared']:.3f}\"\n",
    "    )\n",
    "\n",
    "    ax_info.axis(\"off\")\n",
    "    ax_info.text(\n",
    "        0.05,\n",
    "        0.95,\n",
    "        opis,\n",
    "        transform=ax_info.transAxes,\n",
    "        va=\"top\",\n",
    "        fontsize=12,\n",
    "        bbox={\"boxstyle\": \"round,pad=0.55\", \"fc\": \"white\", \"ec\": \"0.35\"},\n",
    "    )\n",
    "    ax_info.text(\n",
    "        0.05,\n",
    "        0.08,\n",
    "        \"Stała skala osi: zmiana kształtu krzywych nie wynika z autoskalowania.\",\n",
    "        transform=ax_info.transAxes,\n",
    "        va=\"bottom\",\n",
    "        fontsize=10,\n",
    "        color=\"0.25\",\n",
    "    )\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "if WIDGETY_DOSTEPNE:\n",
    "    widget_style = {\"description_width\": \"initial\"}\n",
    "    controls_power = {\n",
    "        \"odstep_srednich\": widgets.FloatSlider(\n",
    "            value=2.0,\n",
    "            min=0.0,\n",
    "            max=5.0,\n",
    "            step=0.25,\n",
    "            description=\"odstęp średnich\",\n",
    "            continuous_update=False,\n",
    "            style=widget_style,\n",
    "        ),\n",
    "        \"sigma\": widgets.FloatSlider(\n",
    "            value=5.0,\n",
    "            min=3.0,\n",
    "            max=8.0,\n",
    "            step=0.25,\n",
    "            description=\"σ w populacji\",\n",
    "            continuous_update=False,\n",
    "            style=widget_style,\n",
    "        ),\n",
    "        \"n\": widgets.IntSlider(\n",
    "            value=30,\n",
    "            min=5,\n",
    "            max=80,\n",
    "            step=5,\n",
    "            description=\"n w grupie\",\n",
    "            continuous_update=False,\n",
    "            style=widget_style,\n",
    "        ),\n",
    "    }\n",
    "\n",
    "    output_power = widgets.interactive_output(rysuj_widget_moc, controls_power)\n",
    "    display(widgets.VBox([widgets.HBox(list(controls_power.values())), output_power]))\n",
    "else:\n",
    "    rysuj_widget_moc()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b202015",
   "metadata": {},
   "source": [
    "## Analiza wariancji a model liniowy\n",
    "\n",
    "Między analizą wariancji i modelem liniowym istnieją bardzo bliskie związki.\n",
    "Można nawet powiedzieć, że jednoczynnikowa ANOVA to szczególny przypadek modelu liniowego.\n",
    "\n",
    "Historycznie te tematy bywały omawiane oddzielnie, ale matematycznie są spójne.\n",
    "\n",
    "Zobaczmy podsumowanie modelu liniowego z jednym kategorialnym predyktorem:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ab32653",
   "metadata": {},
   "outputs": [],
   "source": [
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93518357",
   "metadata": {},
   "source": [
    "Zwróćmy uwagę na kilka rzeczy.\n",
    "\n",
    "Dla analizy wariancji hipoteza zerowa wygląda tak:\n",
    "\n",
    "$$H_0: \\mu_1 = \\mu_2 = \\mu_3$$\n",
    "\n",
    "Dla modelu liniowego (w ujęciu „testu globalnego” dla modelu) hipoteza zerowa jest równoważna:\n",
    "\n",
    "$$H_0: R^2 = 0$$\n",
    "\n",
    "Te hipotezy, chociaż brzmią różnie, są matematycznie równoważne.\n",
    "\n",
    "Jeżeli średnie w grupach są w populacji identyczne, to informacja o tym, do której grupy należy obserwacja, nie pozwala przewidywać jej wartości lepiej niż sama średnia ogólna.\n",
    "Wtedy model nie wyjaśnia wariancji ($R^2 = 0$).\n",
    "\n",
    "W jednoczynnikowej ANOVA $R^2$ z modelu OLS jest równy $\\eta^2$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "350391cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "r2 = float(model.rsquared)\n",
    "\n",
    "r2, eta_squared, abs(r2 - eta_squared)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "312e0be3",
   "metadata": {},
   "source": [
    "## Ćwiczenia\n",
    "\n",
    "1. Wróć do pierwszego przykładu z 5 grupami i:\n",
    "   - zwiększ odchylenie standardowe w jednej grupie (np. zamiast `scale=5` użyj `scale=10` tylko dla jednej grupy),\n",
    "   - zobacz, jak zmienia się tabela ANOVA i wnioski.\n",
    "\n",
    "2. W danych `PlantGrowth`:\n",
    "   - policz wszystkie pary różnic średnich (ctrl–trt1, ctrl–trt2, trt1–trt2),\n",
    "   - porównaj surowe przedziały ufności (bez korekty) z wynikami Tukeya.\n",
    "\n",
    "3. Zastanów się: co w praktyce oznacza $\\eta^2 \\approx 0.26$ w tym przykładzie?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "813fb168",
   "metadata": {},
   "source": [
    "## Podsumowanie\n",
    "\n",
    "- ANOVA testuje hipotezę: „wszystkie średnie są równe” vs „nie wszystkie są równe”.\n",
    "- Statystyka $F$ porównuje zmienność między grupami do zmienności w grupach.\n",
    "- p-wartość to informacja warunkowa na $H_0$, a nie miara „prawdziwości” hipotez.\n",
    "- Po istotnym wyniku ANOVA często robimy analizę *post-hoc* (np. Tukey HSD), żeby ustalić, które grupy się różnią.\n",
    "- W jednoczynnikowej ANOVA $\\eta^2$ jest równy $R^2$ z modelu liniowego z predyktorem kategorialnym.\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
}