{ "cells": [ { "cell_type": "markdown", "id": "9d0adffe", "metadata": {}, "source": [ "# Testy chi-kwadrat\n", "\n", "W tym notatniku poznasz dwa podstawowe zastosowania testów opartych o statystykę $\\chi^2$:\n", "\n", "- **test zgodności** — czy rozkład zliczeń w kategoriach jest zgodny z rozkładem oczekiwanym?\n", "- **test niezależności** — czy dwie zmienne kategoryczne są ze sobą powiązane (tabela kontyngencji)?\n", "\n", "W kognitywistyce i psychologii testy $\\chi^2$ pojawiają się np. w analizie odpowiedzi kategorialnych (poprawnie/niepoprawnie, wybór etykiety, kategoria percepcyjna), a także w porównaniach częstości między warunkami eksperymentu.\n" ] }, { "cell_type": "markdown", "id": "2e237757", "metadata": {}, "source": [ "## Cele\n", "\n", "- Umiesz przeprowadzić test $\\chi^2$ **krok po kroku**: pytanie → $H_0/H_1$ → $\\alpha$ → rozkład pod $H_0$ → region krytyczny → p-wartość obliczona ręcznie → weryfikacja funkcją → interpretacja.\n", "- Potrafisz policzyć statystykę testową $\\chi^2$ ze wzoru oraz wyznaczyć stopnie swobody.\n", "- Umiesz policzyć wartości oczekiwane w tabeli kontyngencji: $E_{ij} = \\frac{R_i C_j}{N}$.\n", "- Rozumiesz, co znaczy p-wartość: $P(\\text{wynik co najmniej tak skrajny jak zaobserwowany} \\mid H_0)$.\n", "- Raportujesz też wielkość efektu (np. $w$ Cohena, $V$ Craméra) i odróżniasz ją od istotności statystycznej.\n" ] }, { "cell_type": "markdown", "id": "b1f23b60", "metadata": {}, "source": [ "## Wymagania wstępne\n", "\n", "- Podstawy testowania hipotez: $H_0/H_1$, poziom istotności $\\alpha$, p-wartość.\n", "- Podstawy rozkładów i CTG (wystarczy intuicja).\n", "- Podstawy pracy z `numpy` i `pandas`.\n", "- (Opcjonalnie) Uruchamianie widżetów w notatniku Jupyter (pakiet `ipywidgets`)." ] }, { "cell_type": "markdown", "id": "a883a32d", "metadata": {}, "source": [ "## Wprowadzenie\n", "\n", "Testy oparte o statystykę $\\chi^2$ pojawiają się wszędzie tam, gdzie pracujemy ze **zliczeniami** (liczebnościami) w kategoriach.\n", "\n", "Ważna obserwacja: statystyka $\\chi^2$ jest sumą kwadratów, więc jest **nieujemna**.\n", "\n", "Z tego powodu testy $\\chi^2$ są w praktyce **jednostronne (prawostronne)** — interesują nas „duże” wartości statystyki.\n" ] }, { "cell_type": "markdown", "id": "ebcf1dbc", "metadata": {}, "source": [ "### Gdzie w kognitywistyce spotkasz test chi-kwadrat?\n", "\n", "Testy $\\chi^2$ stosuje się wtedy, gdy wynik jest **kategoryczny**, a my zliczamy odpowiedzi w kategoriach. Przykłady (bardzo uproszczone):\n", "\n", "- zadanie **dwualternatywnego wyboru wymuszonego** (często zapisuje się skrótem *2AFC*): w każdej próbie wybór jednej z dwóch opcji, np. „lewo/prawo”,\n", "- wybór etykiety emocji (np. Szczęście/Złość/Strach/Smutek),\n", "- kategoryzacja percepcyjna (np. „słyszałem /ba/” vs „słyszałem /da/” w bodźcu audiowizualnym),\n", "- obszary zainteresowania w okulografii (np. pierwsza fiksacja: oczy/usta/tło),\n", "- typ odpowiedzi w zadaniu wykrywania (np. trafienie/fałszywy alarm/chybienie/poprawne odrzucenie).\n", "\n", "**Uwaga dydaktyczna.** W testach $\\chi^2$ zakładamy niezależność obserwacji. W realnych eksperymentach często mamy wiele prób na osobę; wtedy „surowe zliczenia prób” nie są niezależne. W tym notatniku dla prostoty myślimy o zliczeniach na poziomie osób albo o pojedynczej próbie na osobę.\n" ] }, { "cell_type": "markdown", "id": "c1d63c60", "metadata": {}, "source": [ "## Test $\\chi^2$\n", "\n", "Przypomnijmy sobie na początek podstawowy schemat statystycznego testowania hipotez.\n" ] }, { "cell_type": "markdown", "id": "97ee3155", "metadata": {}, "source": [ "### Testowanie hipotez (przypomnienie)\n", "\n", "0. Wymyślamy naszą hipotezę alternatywną (ta nas interesuje!).\n", "1. Przyjmujemy hipotezę zerową (i ew. dodatkowe założenia).\n", "2. Decydujemy się na określony poziom istotności $\\alpha$.\n", "3. Konstruujemy rozkład interesującej nas statystyki z próby (przy założeniu $H_0$ i ew. innych założeniach) (dziś będzie to statystyka testowa $\\chi^2$).\n", "4. Określamy region krytyczny (biorąc pod uwagę $\\alpha$ i ewentualnie kierunkowość hipotezy) (przy teście $\\chi^2$ kierunkowość nie ma znaczenia — test ten jest zawsze prawostronny).\n", "5. Losujemy próbę, zbieramy dane, liczymy statystykę.\n", "6. Odrzucamy bądź nie $H_0$.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "e91a695f", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "from scipy import stats\n", "\n", "from IPython.display import display\n", "\n", "SEED = 20250116\n", "\n", "sns.set_theme(style=\"whitegrid\")\n", "\n", "try:\n", " import ipywidgets as widgets\n", "\n", " WIDGETY_DOSTEPNE = True\n", "except ImportError:\n", " widgets = None\n", " WIDGETY_DOSTEPNE = False\n" ] }, { "cell_type": "markdown", "id": "ed37ec0d", "metadata": {}, "source": [ "### Statystyka testowa $\\chi^2$\n", "\n", "Poniżej podany jest wzór na statystykę testową $\\chi^2$:\n", "\n", "$$\n", "\\chi^2 = \\sum \\frac{(O-E)^2}{E}\n", "$$\n", "\n", "gdzie:\n", "\n", "- $O$ — zaobserwowana liczba wystąpień określonego zdarzenia,\n", "- $E$ — oczekiwana/teoretyczna liczba wystąpień określonego zdarzenia.\n", "\n", "Kształt rozkładu teoretycznego $\\chi^2$ zależy *wyłącznie od liczby stopni swobody*.\n" ] }, { "cell_type": "markdown", "id": "95fd95e4", "metadata": {}, "source": [ "Rozkład $\\chi^2$ będzie dla nas bardzo przydatny w testowaniu hipotez ze względu na to, że rozkład $\\chi^2$ to rozkład zmiennej losowej, która jest sumą $k$ kwadratów niezależnych zmiennych losowych o standardowym rozkładzie normalnym.\n", "\n", "W tej chwili może wydawać się to niezbyt przydatna informacja, jak jednak zobaczymy, jest ona kluczowa dla logiki działania testów $\\chi^2$.\n", "\n", "Wrócimy do tego jeszcze za chwilę.\n", "\n", "Teraz wyobraźmy sobie, że rzucamy 100 razy (uczciwą) monetą.\n", "\n", "Spodziewamy się, że orzeł wypadnie około 50 razy. Jest to teoretyczna/oczekiwana liczba orłów w naszym eksperymencie.\n", "\n", "Obliczmy różnicę pomiędzy faktyczną liczbą orłów a teoretyczną liczbą orłów i podzielmy ją przez pierwiastek z teoretycznej liczby orłów.\n", "\n", "Można dowieść z Centralnego Twierdzenia Granicznego (w bardziej „zaawansowanym” sformułowaniu, niż omawialiśmy), że taka zmienna wraz ze wzrostem liczebności próby będzie dążyć do rozkładu normalnego $N(0, 1-p)$.\n", "\n", "My zamiast dowodzić, pokażemy, że ta zależność działa, korzystając z prostej symulacji.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "1ff29656", "metadata": {}, "outputs": [], "source": [ "rng = np.random.default_rng(SEED)\n", "\n", "n = 100\n", "p = 0.5\n", "expected_heads = n * p\n", "\n", "repeats = 100_000\n", "heads = rng.binomial(n=n, p=p, size=repeats)\n", "\n", "standardized_diff = (heads - expected_heads) / np.sqrt(expected_heads)\n", "\n", "print(\"Średnia (symulacja):\", float(standardized_diff.mean()))\n", "print(\"Odchylenie standardowe (symulacja):\", float(standardized_diff.std(ddof=0)))\n", "\n", "plt.figure(figsize=(9, 4))\n", "plt.hist(\n", " standardized_diff,\n", " bins=60,\n", " density=True,\n", " color=\"lightgrey\",\n", " edgecolor=\"white\",\n", ")\n", "\n", "x = np.linspace(standardized_diff.min(), standardized_diff.max(), 400)\n", "normal_sd = np.sqrt(1 - p)\n", "plt.plot(\n", " x,\n", " stats.norm.pdf(x, loc=0.0, scale=normal_sd),\n", " color=\"red\",\n", " linewidth=2,\n", " label=r\"$N(0,\\sqrt{1-p})$\",\n", ")\n", "\n", "sns.kdeplot(\n", " standardized_diff,\n", " color=\"blue\",\n", " linewidth=2,\n", " label=\"gęstość empiryczna (jądrowa)\",\n", ")\n", "\n", "plt.title(\"Różnice między liczbą trafień a wartością oczekiwaną (n=100)\")\n", "plt.xlabel(r\"$(k-np)/\\sqrt{np}$\")\n", "plt.ylabel(\"gęstość\")\n", "plt.legend()\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "c0a096a5", "metadata": {}, "source": [ "To „pofalowanie” krzywej gęstości naszych symulowanych danych wynika oczywiście z faktu, że nasza zmienna nie jest tak naprawdę ciągła — może przyjmować pewne wartości (liczby naturalne od $-50$ do $50$ podzielone przez pierwiastek z $50$), a pewnych nie może.\n", "\n", "Można łatwo stwierdzić to, patrząc na wzór — jeśli $O$ może być tylko liczbą całkowitą, to $\\frac{(O-E)^2}{E}$ nie będzie mogło przyjąć wszystkich wartości ze zbioru dodatnich liczb rzeczywistych.\n", "\n", "Wraz ze wzrostem liczebności próby może przyjmować coraz więcej wartości, więc krzywa w sposób naturalny „wygładzi” się.\n", "\n", "Przy $N=10000$ już wygląda zupełnie jak rozkład normalny.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "cc15cc04", "metadata": {}, "outputs": [], "source": [ "rng = np.random.default_rng(SEED + 1)\n", "\n", "n = 10_000\n", "p = 0.5\n", "expected_heads = n * p\n", "\n", "repeats = 10_000\n", "heads = rng.binomial(n=n, p=p, size=repeats)\n", "standardized_diff = (heads - expected_heads) / np.sqrt(expected_heads)\n", "\n", "plt.figure(figsize=(9, 4))\n", "plt.hist(\n", " standardized_diff,\n", " bins=60,\n", " density=True,\n", " color=\"lightgrey\",\n", " edgecolor=\"white\",\n", ")\n", "\n", "x = np.linspace(standardized_diff.min(), standardized_diff.max(), 400)\n", "normal_sd = np.sqrt(1 - p)\n", "plt.plot(\n", " x,\n", " stats.norm.pdf(x, loc=0.0, scale=normal_sd),\n", " color=\"red\",\n", " linewidth=2,\n", " label=r\"$N(0,\\sqrt{1-p})$\",\n", ")\n", "\n", "sns.kdeplot(\n", " standardized_diff,\n", " color=\"blue\",\n", " linewidth=2,\n", " label=\"gęstość empiryczna (jądrowa)\",\n", ")\n", "\n", "plt.title(\"Różnice między liczbą trafień a wartością oczekiwaną (n=10 000)\")\n", "plt.xlabel(r\"$(k-np)/\\sqrt{np}$\")\n", "plt.ylabel(\"gęstość\")\n", "plt.legend()\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "4e7aaec2", "metadata": {}, "source": [ "W testach opartych na statystyce testowej $\\chi^2$ będzie nas interesował rozkład sumy kwadratów kilku takich zmiennych.\n", "\n", "Można dowieść, że dla $r$ takich zmiennych rozkład taki ma postać $\\chi^2(r-1)$ — tzn. jest to rozkład $\\chi^2$ o $r-1$ stopniach swobody.\n", "\n", "Nie będziemy tego dowodzić; dowód można znaleźć np. tu:\n", "\n", "https://ocw.mit.edu/courses/mathematics/18-443-statistics-for-applications-fall-2003/lecture-notes/lec23.pdf\n", "\n", "(prawdopodobnie nigdy nie będą Państwo tego dowodzić).\n", "\n", "Przetestujmy tę „prawdę objawioną” na przykładzie dwóch zmiennych.\n", "\n", "Jak widzimy, rozkład empiryczny $\\chi^2$, który otrzymamy za pomocą symulacji (histogram i niebieska krzywa), odpowiada rozkładowi teoretycznemu (czerwona krzywa).\n" ] }, { "cell_type": "code", "execution_count": null, "id": "ae344521", "metadata": {}, "outputs": [], "source": [ "rng = np.random.default_rng(SEED + 2)\n", "\n", "repeats = 100_000\n", "n = 1_000\n", "p = 0.5\n", "expected = n * p\n", "\n", "heads = rng.binomial(n=n, p=p, size=repeats)\n", "tails = n - heads\n", "\n", "chi2_stats = ((heads - expected) ** 2) / expected + ((tails - expected) ** 2) / expected\n", "\n", "plt.figure(figsize=(9, 4))\n", "plt.hist(chi2_stats, bins=80, density=True, color=\"lightgrey\", edgecolor=\"white\")\n", "\n", "x = np.linspace(0, np.quantile(chi2_stats, 0.995), 400)\n", "plt.plot(x, stats.chi2.pdf(x, df=1), color=\"red\", linewidth=2, label=r\"$\\chi^2(1)$\")\n", "\n", "sns.kdeplot(chi2_stats, color=\"blue\", linewidth=2, label=\"gęstość empiryczna (jądrowa)\")\n", "\n", "plt.title(\"Kwadraty różnic (symulacja) a rozkład $\\chi^2$\")\n", "plt.xlabel(r\"$\\chi^2$\")\n", "plt.ylabel(\"gęstość\")\n", "plt.legend()\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "93404c83", "metadata": {}, "source": [ "Na poniższej ilustracji przedstawione mamy wykresy rozkładu $\\chi^2$ dla różnej liczby stopni swobody.\n", "\n", "Na każdy z wykresów nałożona została również pionowa linia odcinająca górne 5% rozkładu.\n", "\n", "Tam będą wartości statystyki testowej uprawniające nas do odrzucenia hipotezy zerowej (oczywiście tylko wtedy, gdy przyjmiemy $\\alpha = 0.05$).\n" ] }, { "cell_type": "code", "execution_count": null, "id": "99cc1469", "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(2, 2, figsize=(10, 6), sharex=True, sharey=True)\n", "axes = axes.ravel()\n", "\n", "x = np.linspace(0, 20, 500)\n", "alpha = 0.05\n", "\n", "for df, ax in zip([1, 2, 3, 4], axes, strict=True):\n", " y = stats.chi2.pdf(x, df=df)\n", " ax.plot(x, y, color=\"black\")\n", "\n", " critical = stats.chi2.ppf(1 - alpha, df=df)\n", " ax.axvline(critical, color=\"grey\", linestyle=\":\", linewidth=2)\n", "\n", " ax.set_title(f\"Liczba stopni swobody: {df}\")\n", " ax.text(critical + 0.4, 0.10, f\"{critical:.3f}\")\n", "\n", "fig.suptitle(\"Rozkład $\\chi^2$ i wartość krytyczna (górne 5%)\")\n", "for ax in axes:\n", " ax.set_xlabel(\"$\\chi^2$\")\n", " ax.set_ylabel(\"gęstość\")\n", "\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "05acddbf", "metadata": {}, "source": [ "### Interaktywna ilustracja: rozkład $\\chi^2$ i ogon krytyczny\n", "\n", "Poniższy widżet pozwala zmieniać liczbę stopni swobody, poziom istotności $\\alpha$ oraz obserwowaną wartość statystyki $\\chi^2$.\n", "\n", "Zwróć uwagę, że p-wartość jest **polem pod prawym ogonem** rozkładu (dla testów $\\chi^2$ interesuje nas prawy ogon).\n" ] }, { "cell_type": "code", "execution_count": null, "id": "73c75179", "metadata": {}, "outputs": [], "source": [ "if not WIDGETY_DOSTEPNE:\n", " print(\n", " \"Widżety nie są dostępne. Jeśli pracujesz w Jupyterze, doinstaluj pakiet ipywidgets, \"\n", " \"aby korzystać z części interaktywnych.\"\n", " )\n", "else:\n", " df_slider = widgets.IntSlider(\n", " value=3,\n", " min=1,\n", " max=20,\n", " step=1,\n", " description=\"df\",\n", " continuous_update=False,\n", " )\n", " alpha_slider = widgets.FloatSlider(\n", " value=0.05,\n", " min=0.001,\n", " max=0.20,\n", " step=0.001,\n", " description=\"α\",\n", " readout_format=\".3f\",\n", " continuous_update=False,\n", " )\n", " chi2_slider = widgets.FloatSlider(\n", " value=5.0,\n", " min=0.0,\n", " max=float(stats.chi2.ppf(0.995, df=df_slider.value)),\n", " step=0.1,\n", " description=\"χ²_obs\",\n", " readout_format=\".1f\",\n", " continuous_update=False,\n", " )\n", "\n", " def aktualizuj_chi2_max(change=None):\n", " chi2_slider.max = float(stats.chi2.ppf(0.995, df=int(df_slider.value)))\n", " if chi2_slider.value > chi2_slider.max:\n", " chi2_slider.value = chi2_slider.max\n", "\n", " df_slider.observe(aktualizuj_chi2_max, names=\"value\")\n", " aktualizuj_chi2_max()\n", "\n", " def pokaz_rozklad_chi2(df, alpha, chi2_obs):\n", " df = int(df)\n", " alpha = float(alpha)\n", " chi2_obs = float(chi2_obs)\n", "\n", " critical = float(stats.chi2.ppf(1 - alpha, df=df))\n", " p_value = float(stats.chi2.sf(chi2_obs, df=df))\n", "\n", " print(f\"df = {df}\")\n", " print(f\"α = {alpha:.3f}\")\n", " print(f\"χ²_obs = {chi2_obs:.3f}\")\n", " print(f\"Wartość krytyczna χ²_(1-α, df) = {critical:.3f}\")\n", " print(f\"p-wartość = {p_value:.6f}\")\n", "\n", " x_max = float(stats.chi2.ppf(0.999, df=df))\n", " x = np.linspace(0, x_max, 600)\n", " y = stats.chi2.pdf(x, df=df)\n", "\n", " plt.figure(figsize=(9, 4))\n", " plt.plot(x, y, color=\"black\", linewidth=2)\n", "\n", " mask = x >= chi2_obs\n", " plt.fill_between(x[mask], y[mask], color=\"skyblue\", alpha=0.6, label=\"p-wartość\")\n", "\n", " plt.axvline(chi2_obs, color=\"navy\", linestyle=\"-\", linewidth=2)\n", " plt.axvline(critical, color=\"grey\", linestyle=\":\", linewidth=2, label=\"wartość krytyczna\")\n", "\n", " plt.title(\"Rozkład $\\\\chi^2$ i prawy ogon\")\n", " plt.xlabel(r\"$\\\\chi^2$\")\n", " plt.ylabel(\"gęstość\")\n", " plt.legend()\n", " plt.tight_layout()\n", " plt.show()\n", "\n", " out = widgets.interactive_output(\n", " pokaz_rozklad_chi2,\n", " {\n", " \"df\": df_slider,\n", " \"alpha\": alpha_slider,\n", " \"chi2_obs\": chi2_slider,\n", " },\n", " )\n", "\n", " display(widgets.VBox([df_slider, alpha_slider, chi2_slider]), out)\n" ] }, { "cell_type": "markdown", "id": "da1ad57f", "metadata": {}, "source": [ "## Testowanie hipotez za pomocą testu $\\chi^2$\n", "\n", "Na tych zajęciach będziemy zajmować się dwoma testami opartymi na statystyce testowej $\\chi^2$:\n", "\n", "- **test zgodności** — czy rozkład zliczeń jest zgodny z rozkładem oczekiwanym?\n", "- **test niezależności** — czy dwie zmienne kategoryczne są ze sobą powiązane (tabela kontyngencji)?\n", "\n", "Z obliczeniowego punktu widzenia testy te działają podobnie, ale stosuje się je w innych sytuacjach i inaczej formułuje hipotezę zerową." ] }, { "cell_type": "markdown", "id": "89cb9a67", "metadata": {}, "source": [ "### Przykład I (test zgodności): zadanie dwualternatywnego wyboru wymuszonego\n", "\n", "W zadaniu **dwualternatywnego wyboru wymuszonego** (często zapisuje się skrótem *2AFC*) osoba badana w każdej próbie musi wybrać **jedną z dwóch** opcji, np. „lewo/prawo”.\n", "\n", "Jeśli w bodźcu nie ma użytecznej informacji (albo jest on zbyt słaby), uczestnik zgaduje. Wtedy oczekujemy trafień z prawdopodobieństwem $p=0{,}5$.\n", "\n", "Załóżmy, że w prostym eksperymencie zebraliśmy 280 prób i zapisaliśmy, czy odpowiedź była poprawna.\n", "\n", "- poprawnie: 123,\n", "- niepoprawnie: 157.\n", "\n", "**Pytanie:** czy taki rozkład zliczeń jest zgodny z losowym zgadywaniem?\n", "\n", "**Krok po kroku (NHST):**\n", "\n", "1. Hipotezy:\n", " - $H_0$: rozkład zgodny z losowym zgadywaniem ($p=0{,}5$),\n", " - $H_1$: rozkład nie jest zgodny z losowym zgadywaniem.\n", "2. Poziom istotności: $\\alpha = 0{,}05$.\n", "3. Statystyka: $\\chi^2 = \\sum \\frac{(O-E)^2}{E}$.\n", "4. Rozkład pod $H_0$: $\\chi^2(\\text{df}=r-1)$, gdzie $r$ to liczba kategorii (tu $r=2$, więc df=1).\n", "5. Region krytyczny: $\\chi^2 \\ge \\chi^2_{1-\\alpha,\\,\\text{df}}$.\n", "6. p-wartość: $p = P(\\chi^2 \\ge \\chi^2_{\\text{obs}} \\mid H_0)$.\n", "\n", "Pamiętaj: p-wartość nie jest miarą wielkości efektu ani „prawdopodobieństwem, że $H_0$ jest prawdziwa”.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "25231c01", "metadata": {}, "outputs": [], "source": [ "observed = np.array([123, 157])\n", "expected = np.array([140, 140])\n", "\n", "labels = [\"Poprawnie\", \"Niepoprawnie\"]\n", "\n", "components = (observed - expected) ** 2 / expected\n", "chi2_stat = float(components.sum())\n", "\n", "tabela = pd.DataFrame(\n", " {\n", " \"O (zaobserwowane)\": observed,\n", " \"E (oczekiwane)\": expected,\n", " \"(O-E)^2/E\": components,\n", " },\n", " index=labels,\n", ")\n", "\n", "print(tabela)\n", "print(\"Chi-kwadrat:\", chi2_stat)\n" ] }, { "cell_type": "markdown", "id": "9a181b0e", "metadata": {}, "source": [ "Następnie, używając teoretycznego rozkładu $\\chi^2$, sprawdźmy jakie jest prawdopodobieństwo uzyskania takiego lub bardziej skrajnego wyniku przy założeniu prawdziwości hipotezy zerowej.\n", "\n", "Liczba stopni swobody jaka nas interesuje (tzn. 1) dana jest przez wzór:\n", "\n", "$$\n", "\\text{df} = r - 1,\n", "$$\n", "\n", "gdzie $r$ to liczba kategorii.\n", "\n", "Generalnie zasadą jest to, że liczba stopni swobody równa jest liczbie kategorii, którą zmniejszamy o liczbę stopni swobody, które „utraciliśmy” dopasowując rozkład do danych.\n", "\n", "Obrazowo: gdy mamy 3 kategorie i znamy tylko sumę zliczeń, to jeśli ustalimy wartości w dwóch komórkach, wartość trzeciej jest zdeterminowana.\n", "\n", "W innych kategoriach można myśleć o tym tak: wyobraźmy sobie że rzucamy monetą 100 razy. Czy liczba orłów i liczba reszek są zmiennymi niezależnymi? Oczywiście nie — jeśli znamy liczbę orłów, to wiemy już, jaka jest liczba reszek. Stąd nasza redukcja stopni swobody.\n", "\n", "Liczba stopni swobody przy teście zgodności jest równa liczbie zmiennych, które możemy swobodnie zmieniać zachowując taką samą sumę wszystkich zliczeń.\n", "\n", "(W przypadku tabeli 2x2 i testu niezależności sprawa wygląda troszkę inaczej. Miejsce naszej sumy wszystkich zliczeń zajmują sumy brzegowe. Tracimy stopnie swobody z „wiersza” i z „kolumny”.)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "588b05d0", "metadata": {}, "outputs": [], "source": [ "alpha = 0.05\n", "r = observed.size\n", "\n", "df = r - 1\n", "critical_value = float(stats.chi2.ppf(1 - alpha, df=df))\n", "p_value = float(stats.chi2.sf(chi2_stat, df=df))\n", "\n", "print(\"df =\", df)\n", "print(\"Wartość krytyczna (górne 5%):\", critical_value)\n", "print(\"p-wartość (górny ogon):\", p_value)\n", "\n", "if chi2_stat >= critical_value:\n", " print(\"Decyzja: odrzucamy H0 przy alfa=0.05\")\n", "else:\n", " print(\"Decyzja: nie odrzucamy H0 przy alfa=0.05\")\n" ] }, { "cell_type": "markdown", "id": "328c12eb", "metadata": {}, "source": [ "Dla poziomu istotności statystycznej $\\alpha = 0.05$ odrzucamy (albo nie) hipotezę zerową na podstawie p-wartości lub porównania z wartością krytyczną.\n", "\n", "Niezależnie od decyzji warto od razu zapytać: **jak duża** jest obserwowana różnica? Dlatego obliczymy także wielkość efektu (oraz, opcjonalnie, przedział ufności dla odsetka poprawnych odpowiedzi).\n" ] }, { "cell_type": "code", "execution_count": null, "id": "e2356351", "metadata": {}, "outputs": [], "source": [ "n_total = int(observed.sum())\n", "\n", "cohens_w = float(np.sqrt(chi2_stat / n_total))\n", "\n", "print(\"n =\", n_total)\n", "print(\"Wielkość efektu (w Cohena):\", cohens_w)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "f613b6f7", "metadata": {}, "outputs": [], "source": [ "test_result = stats.chisquare(f_obs=observed, f_exp=expected)\n", "print(test_result)\n" ] }, { "cell_type": "markdown", "id": "6ca5b565", "metadata": {}, "source": [ "Dla dwóch kategorii test zgodności $\\chi^2$ jest ściśle powiązany z testem dwumianowym (dla odsetka trafień).\n", "\n", "Poniżej policzymy test dwumianowy oraz 95% przedział ufności dla $\\hat p$ (odsetka poprawnych odpowiedzi).\n" ] }, { "cell_type": "code", "execution_count": null, "id": "756a7234", "metadata": {}, "outputs": [], "source": [ "n = int(observed.sum())\n", "k_obs = int(observed[0])\n", "p0 = 0.5\n", "\n", "p_hat = k_obs / n\n", "\n", "binom_res = stats.binomtest(k_obs, n=n, p=p0, alternative=\"two-sided\")\n", "ci = binom_res.proportion_ci(confidence_level=0.95, method=\"exact\")\n", "\n", "print(f\"Oszacowanie z próby: p̂ = {p_hat:.3f}\")\n", "print(f\"95% przedział ufności dla p: [{ci.low:.3f}, {ci.high:.3f}]\")\n", "print(f\"Test dwumianowy (dokładny): p = {float(binom_res.pvalue):.6f}\")\n", "\n", "z = (k_obs - n * p0) / np.sqrt(n * p0 * (1 - p0))\n", "print(f\"Aproksymacja normalna: z = {float(z):.3f}, z² = {float(z**2):.3f}\")\n", "print(f\"Statystyka χ² z testu zgodności: {chi2_stat:.3f}\")\n" ] }, { "cell_type": "markdown", "id": "9905128a", "metadata": {}, "source": [ "Rozkład teoretyczny $\\chi^2$ nie zawsze idealnie odpowiada rozkładowi statystyki $\\chi^2$ z próby (np. dla małych liczebności rozkład z próby bywa „poszarpany”, bo pracujemy na zliczeniach całkowitych).\n", "\n", "Czasem można uzyskać dokładniejszy wynik, szacując p-wartość metodą symulacyjną (Monte Carlo).\n", "\n", "Dla ilustracji zrobimy to „ręcznie”, symulując dane pod hipotezą zerową.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "ce1d8e20", "metadata": {}, "outputs": [], "source": [ "rng = np.random.default_rng(SEED + 10)\n", "\n", "repeats = 100_000\n", "n = int(observed.sum())\n", "expected_probs = np.array([0.5, 0.5])\n", "expected_counts = n * expected_probs\n", "\n", "sim_counts = rng.multinomial(n=n, pvals=expected_probs, size=repeats)\n", "\n", "sim_chi2 = ((sim_counts - expected_counts) ** 2 / expected_counts).sum(axis=1)\n", "\n", "p_mc = float((sim_chi2 >= chi2_stat).mean())\n", "print(\"Monte Carlo p ≈\", p_mc)\n" ] }, { "cell_type": "markdown", "id": "a4955ed2", "metadata": {}, "source": [ "### Przykład II (test zgodności): wybór etykiety emocji\n", "\n", "Wyobraźmy sobie bardzo proste zadanie: 80 osób widzi tę samą neutralną twarz i **musi** wybrać jedną z czterech etykiet: Szczęście, Złość, Strach, Smutek.\n", "\n", "To oczywiście uproszczenie, ale dobrze nadaje się do ćwiczenia testu zgodności.\n", "\n", "Hipoteza zerowa: przy braku preferencji każda etykieta powinna pojawiać się równie często (po 25%).\n", "\n", "**Krok po kroku (NHST):**\n", "\n", "1. Hipotezy: $H_0$ — wszystkie kategorie mają prawdopodobieństwo $1/4$; $H_1$ — przynajmniej jedna ma inne prawdopodobieństwo.\n", "2. Poziom istotności: $\\alpha = 0.05$.\n", "3. Statystyka: $\\chi^2 = \\sum \\frac{(O-E)^2}{E}$.\n", "4. Rozkład pod $H_0$: $\\chi^2(\\text{df}=r-1)=\\chi^2(3)$.\n", "5. Region krytyczny: $\\chi^2 \\ge \\chi^2_{1-\\alpha,3}$.\n", "6. p-wartość: $p = P(\\chi^2 \\ge \\chi^2_{\\text{obs}} \\mid H_0)$.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "428f7889", "metadata": {}, "outputs": [], "source": [ "labels = [\"Szczęście\", \"Złość\", \"Strach\", \"Smutek\"]\n", "\n", "observed = np.array([32, 18, 14, 16])\n", "expected_probs = np.array([0.25, 0.25, 0.25, 0.25])\n", "expected = observed.sum() * expected_probs\n", "\n", "components = (observed - expected) ** 2 / expected\n", "chi2_stat = float(components.sum())\n", "\n", "alpha = 0.05\n", "df = len(labels) - 1\n", "\n", "critical_value = float(stats.chi2.ppf(1 - alpha, df=df))\n", "p_value = float(stats.chi2.sf(chi2_stat, df=df))\n", "\n", "tabela = pd.DataFrame(\n", " {\n", " \"O (zaobserwowane)\": observed,\n", " \"E (oczekiwane)\": expected,\n", " \"(O-E)^2/E\": components,\n", " },\n", " index=labels,\n", ")\n", "\n", "print(tabela)\n", "print(\"Chi-kwadrat:\", chi2_stat)\n", "print(\"df =\", df)\n", "print(\"Wartość krytyczna:\", critical_value)\n", "print(\"p-wartość:\", p_value)\n", "\n", "n_total = int(observed.sum())\n", "print(\"w Cohena =\", float(np.sqrt(chi2_stat / n_total)))\n", "\n", "x = np.arange(len(labels))\n", "width = 0.42\n", "\n", "plt.figure(figsize=(8, 4))\n", "plt.bar(\n", " x - width / 2,\n", " observed,\n", " width=width,\n", " color=\"skyblue\",\n", " edgecolor=\"white\",\n", " label=\"O (zaobserwowane)\",\n", ")\n", "plt.bar(\n", " x + width / 2,\n", " expected,\n", " width=width,\n", " color=\"lightgrey\",\n", " edgecolor=\"white\",\n", " label=\"E (oczekiwane)\",\n", ")\n", "plt.xticks(x, labels)\n", "plt.title(\"Test zgodności: zliczenia w kategoriach\")\n", "plt.xlabel(\"kategoria\")\n", "plt.ylabel(\"liczebność\")\n", "plt.legend()\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "code", "execution_count": null, "id": "27f97248", "metadata": {}, "outputs": [], "source": [ "test_result = stats.chisquare(f_obs=observed, f_exp=expected)\n", "print(test_result)\n" ] }, { "cell_type": "markdown", "id": "22cf4b4f", "metadata": {}, "source": [ "Ponownie możemy obliczyć wartość *p* za pomocą metody Monte Carlo (czyli symulacji).\n", "\n", "Zasadniczo nie będziemy robić tego w praktyce, ale warto zrobić to kilka razy żeby uzmysłowić sobie, że rozkład teoretyczny $\\chi^2$ i rozkład statystyki testowej $\\chi^2$ nie zawsze są dokładnie tym samym rozkładem.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "71005f84", "metadata": {}, "outputs": [], "source": [ "rng = np.random.default_rng(SEED + 11)\n", "\n", "repeats = 100_000\n", "n = int(observed.sum())\n", "expected_probs = np.array([0.25, 0.25, 0.25, 0.25])\n", "expected_counts = n * expected_probs\n", "\n", "sim_counts = rng.multinomial(n=n, pvals=expected_probs, size=repeats)\n", "\n", "sim_chi2 = ((sim_counts - expected_counts) ** 2 / expected_counts).sum(axis=1)\n", "\n", "p_mc = float((sim_chi2 >= chi2_stat).mean())\n", "print(\"p (Monte Carlo) ≈\", p_mc)\n" ] }, { "cell_type": "markdown", "id": "6b1d2be5", "metadata": {}, "source": [ "Warto zwrócić uwagę na to, że nie zawsze teoretyczne prawdopodobieństwa będą równe i wobec tego nie zawsze będą równe częstości teoretyczne.\n", "\n", "W zadaniach spotyka się sytuacje, gdzie pod $H_0$ mamy np. $p = (0.2, 0.5, 0.3)$; wtedy $E_i = n p_i$.\n" ] }, { "cell_type": "markdown", "id": "285fbf66", "metadata": {}, "source": [ "### Interaktywna symulacja: kiedy statystyka $\\chi^2$ ma rozkład zbliżony do teoretycznego?\n", "\n", "W teście zgodności często korzystamy z faktu, że (przy spełnieniu założeń i dostatecznie dużych liczebnościach) rozkład statystyki z próby jest zbliżony do rozkładu teoretycznego $\\chi^2$.\n", "\n", "Poniżej możesz zobaczyć, jak zmienia się histogram statystyki $\\chi^2$ (dla czterech kategorii) wraz ze wzrostem liczebności $n$.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "ecd1532d", "metadata": {}, "outputs": [], "source": [ "if not WIDGETY_DOSTEPNE:\n", " print(\n", " \"Widżety nie są dostępne. Jeśli pracujesz w Jupyterze, doinstaluj pakiet ipywidgets, \"\n", " \"aby korzystać z części interaktywnych.\"\n", " )\n", "else:\n", " n_slider = widgets.IntSlider(\n", " value=40,\n", " min=10,\n", " max=500,\n", " step=10,\n", " description=\"n\",\n", " continuous_update=False,\n", " )\n", " repeats_slider = widgets.IntSlider(\n", " value=2000,\n", " min=200,\n", " max=20000,\n", " step=200,\n", " description=\"Powtórzenia\",\n", " continuous_update=False,\n", " )\n", " seed_input = widgets.IntText(value=SEED, description=\"Ziarno\")\n", "\n", " def symuluj_rozklad_chi2(n, repeats, seed):\n", " n = int(n)\n", " repeats = int(repeats)\n", " seed = int(seed)\n", "\n", " rng_local = np.random.default_rng(seed)\n", "\n", " probs = np.array([0.25, 0.25, 0.25, 0.25])\n", " expected = n * probs\n", "\n", " sim_counts = rng_local.multinomial(n=n, pvals=probs, size=repeats)\n", " sim_chi2 = ((sim_counts - expected) ** 2 / expected).sum(axis=1)\n", "\n", " df = len(probs) - 1\n", "\n", " plt.figure(figsize=(9, 4))\n", " plt.hist(\n", " sim_chi2,\n", " bins=40,\n", " density=True,\n", " color=\"lightgrey\",\n", " edgecolor=\"white\",\n", " label=\"symulacja\",\n", " )\n", "\n", " x = np.linspace(0, np.quantile(sim_chi2, 0.99), 400)\n", " plt.plot(\n", " x,\n", " stats.chi2.pdf(x, df=df),\n", " color=\"crimson\",\n", " linewidth=2,\n", " label=f\"gęstość $\\\\chi^2$(df={df})\",\n", " )\n", "\n", " plt.title(\"Rozkład statystyki $\\\\chi^2$ w teście zgodności (gdy H0 jest prawdziwa)\")\n", " plt.xlabel(r\"$\\\\chi^2$\")\n", " plt.ylabel(\"gęstość\")\n", " plt.legend()\n", " plt.tight_layout()\n", " plt.show()\n", "\n", " print(f\"Średnia symulowana: {float(sim_chi2.mean()):.3f} (teoria: {df})\")\n", " print(f\"Odchylenie standardowe (symulacja): {float(sim_chi2.std(ddof=0)):.3f}\")\n", "\n", " out = widgets.interactive_output(\n", " symuluj_rozklad_chi2,\n", " {\n", " \"n\": n_slider,\n", " \"repeats\": repeats_slider,\n", " \"seed\": seed_input,\n", " },\n", " )\n", "\n", " display(widgets.VBox([n_slider, repeats_slider, seed_input]), out)\n" ] }, { "cell_type": "markdown", "id": "fbebdaa2", "metadata": {}, "source": [ "### Przykład III (test niezależności): bodźce audiowizualne i kategoria odpowiedzi\n", "\n", "W teście niezależności $\\chi^2$ sprawdzamy, czy dwie zmienne kategoryczne są powiązane.\n", "\n", "Typowy scenariusz w kognitywistyce: **warunek bodźcowy** (np. bodziec zgodny vs niezgodny) oraz **kategoria odpowiedzi** (np. „słyszałem /ba/” vs „słyszałem /da/”).\n", "\n", "Efekt McGurka (w bardzo uproszczonej wersji) pokazuje, że gdy informacja słuchowa i wzrokowa jest niezgodna, rozkład odpowiedzi może się zmieniać.\n", "\n", "Załóżmy, że 100 osób zobaczyło po jednym bodźcu:\n", "\n", "- warunek 1: bodziec zgodny (n = 50),\n", "- warunek 2: bodziec niezgodny (n = 50),\n", "\n", "a odpowiedź zaklasyfikowaliśmy do jednej z dwóch kategorii.\n" ] }, { "cell_type": "markdown", "id": "62c3fe2c", "metadata": {}, "source": [ "Tabela kontyngencji (zliczenia):\n", "\n", "| Warunek | „/ba/” | „/da/” | Suma |\n", "|---|---:|---:|---:|\n", "| Zgodny | 45 | 5 | 50 |\n", "| Niezgodny | 20 | 30 | 50 |\n", "| Suma | 65 | 35 | 100 |\n", "\n", "Hipotezy:\n", "\n", "- $H_0$: warunek i odpowiedź są niezależne,\n", "- $H_1$: warunek i odpowiedź są zależne.\n", "\n", "**Krok po kroku (NHST):**\n", "\n", "1. Poziom istotności: $\\alpha = 0{,}05$.\n", "2. Rozkład pod $H_0$: $\\chi^2(\\text{df}=(r-1)(c-1))=\\chi^2(1)$.\n", "3. Region krytyczny: $\\chi^2 \\ge \\chi^2_{1-\\alpha,1}$.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "6c6c2ad5", "metadata": {}, "outputs": [], "source": [ "observed = np.array([[45, 5], [20, 30]])\n", "\n", "observed_df = pd.DataFrame(\n", " observed,\n", " index=[\"Bodziec zgodny\", \"Bodziec niezgodny\"],\n", " columns=[\"Odpowiedź: /ba/\", \"Odpowiedź: /da/\"],\n", ")\n", "\n", "print(\"Zaobserwowane liczebności (O):\")\n", "print(observed_df)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "0c621b1f", "metadata": {}, "outputs": [], "source": [ "row_totals = observed.sum(axis=1, keepdims=True)\n", "col_totals = observed.sum(axis=0, keepdims=True)\n", "\n", "total = int(observed.sum())\n", "expected = row_totals * col_totals / total\n", "\n", "expected_df = pd.DataFrame(\n", " expected,\n", " index=observed_df.index,\n", " columns=observed_df.columns,\n", ")\n", "\n", "print(\"Liczebności oczekiwane (E) = R_i C_j / N:\")\n", "print(expected_df)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "371adaa7", "metadata": {}, "outputs": [], "source": [ "chi2_stat = float((((observed - expected) ** 2) / expected).sum())\n", "\n", "alpha = 0.05\n", "df = (observed.shape[0] - 1) * (observed.shape[1] - 1)\n", "\n", "critical_value = float(stats.chi2.ppf(1 - alpha, df=df))\n", "p_value = float(stats.chi2.sf(chi2_stat, df=df))\n", "\n", "print(\"Wartość krytyczna:\", critical_value)\n", "print(\"Chi-kwadrat:\", chi2_stat)\n", "print(\"p-wartość:\", p_value)\n", "\n", "n_total = int(observed.sum())\n", "phi = float(np.sqrt(chi2_stat / n_total))\n", "print(\"V Craméra (dla 2×2 = φ):\", phi)\n", "\n", "n_zgodny = int(observed[0].sum())\n", "n_niezgodny = int(observed[1].sum())\n", "\n", "k_da_zgodny = int(observed[0, 1])\n", "k_da_niezgodny = int(observed[1, 1])\n", "\n", "p_hat_zgodny = k_da_zgodny / n_zgodny\n", "p_hat_niezgodny = k_da_niezgodny / n_niezgodny\n", "\n", "ci_zgodny = stats.binomtest(k_da_zgodny, n=n_zgodny).proportion_ci(\n", " confidence_level=0.95, method=\"exact\"\n", ")\n", "ci_niezgodny = stats.binomtest(k_da_niezgodny, n=n_niezgodny).proportion_ci(\n", " confidence_level=0.95, method=\"exact\"\n", ")\n", "\n", "print(\"\\nOdsetek odpowiedzi '/da/':\")\n", "print(\n", " f\"- bodziec zgodny: {p_hat_zgodny:.3f} (95% PU: [{ci_zgodny.low:.3f}, {ci_zgodny.high:.3f}])\"\n", ")\n", "print(\n", " f\"- bodziec niezgodny: {p_hat_niezgodny:.3f} (95% PU: [{ci_niezgodny.low:.3f}, {ci_niezgodny.high:.3f}])\"\n", ")\n", "print(\"Różnica odsetków (niezgodny − zgodny):\", float(p_hat_niezgodny - p_hat_zgodny))\n" ] }, { "cell_type": "markdown", "id": "f6af2f48", "metadata": {}, "source": [ "Do przeprowadzania testu $\\chi^2$ możemy wykorzystać funkcję `chi2_contingency`.\n", "\n", "W przypadku tabeli 2×2 często rozważa się tzw. poprawkę ciągłości (poprawkę Yatesa). Poniżej porównujemy wynik „bez poprawki” i „z poprawką”.\n", "\n", "(Uwaga: poprawka Yatesa bywa krytykowana; na tych zajęciach traktujemy ją jako element praktyki raportowania i okazję do dyskusji.)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "6eb231bf", "metadata": {}, "outputs": [], "source": [ "chi2_corr, p_corr, df_corr, expected_corr = stats.chi2_contingency(\n", " observed, correction=True\n", ")\n", "chi2_nocorr, p_nocorr, df_nocorr, expected_nocorr = stats.chi2_contingency(\n", " observed, correction=False\n", ")\n", "\n", "print(\"Z poprawką Yatesa (ciągłości):\")\n", "print(\"chi2 =\", float(chi2_corr), \"df =\", int(df_corr), \"p =\", float(p_corr))\n", "\n", "print(\"\\nBez poprawki:\")\n", "print(\"chi2 =\", float(chi2_nocorr), \"df =\", int(df_nocorr), \"p =\", float(p_nocorr))\n" ] }, { "cell_type": "markdown", "id": "3cabc581", "metadata": {}, "source": [ "Możemy też oszacować p-wartość metodą Monte Carlo, symulując tabele o tych samych sumach brzegowych.\n", "\n", "Dla tabeli 2×2 (i ustalonych sum brzegowych) rozkład wewnątrz jest powiązany z rozkładem hipergeometrycznym, co ułatwia symulację.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "dae17dcb", "metadata": {}, "outputs": [], "source": [ "rng = np.random.default_rng(SEED + 12)\n", "\n", "repeats = 200_000\n", "\n", "row_sums = observed.sum(axis=1)\n", "col_sums = observed.sum(axis=0)\n", "\n", "total = int(observed.sum())\n", "expected = (row_sums[:, None] * col_sums[None, :]) / total\n", "\n", "r1 = int(row_sums[0])\n", "c1 = int(col_sums[0])\n", "\n", "a_sim = rng.hypergeometric(ngood=c1, nbad=total - c1, nsample=r1, size=repeats)\n", "\n", "b_sim = r1 - a_sim\n", "c_sim = c1 - a_sim\n", "d_sim = total - a_sim - b_sim - c_sim\n", "\n", "sim_tables = np.stack([a_sim, b_sim, c_sim, d_sim], axis=1).reshape(repeats, 2, 2)\n", "\n", "sim_chi2 = ((sim_tables - expected) ** 2 / expected).sum(axis=(1, 2))\n", "\n", "p_mc = float((sim_chi2 >= chi2_stat).mean())\n", "print(\"Zasymulowana wartość p ≈\", p_mc)\n" ] }, { "cell_type": "markdown", "id": "8258cde7", "metadata": {}, "source": [ "### Interaktywna ilustracja: p-wartość a liczebność próby (tabela 2×2)\n", "\n", "Ten widżet symuluje dane w prostym schemacie 2×2: dwa warunki (np. bodziec zgodny/niezgodny) i odpowiedź „/da/” (tak/nie).\n", "\n", "Możesz zmieniać liczebność oraz prawdopodobieństwa odpowiedzi w każdym warunku i obserwować, jak zmieniają się:\n", "\n", "- p-wartość,\n", "- wielkość efektu ($\\varphi$ / $V$ Craméra),\n", "- przedziały ufności dla odsetków.\n", "\n", "To dobry moment, aby zauważyć, że przy bardzo dużych próbach nawet małe różnice mogą dawać małe p-wartości.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "f6849742", "metadata": {}, "outputs": [], "source": [ "if not WIDGETY_DOSTEPNE:\n", " print(\n", " \"Widżety nie są dostępne. Jeśli pracujesz w Jupyterze, doinstaluj pakiet ipywidgets, \"\n", " \"aby korzystać z części interaktywnych.\"\n", " )\n", "else:\n", " n_slider = widgets.IntSlider(\n", " value=50,\n", " min=10,\n", " max=300,\n", " step=5,\n", " description=\"n/warunek\",\n", " continuous_update=False,\n", " )\n", " p_zgodny_slider = widgets.FloatSlider(\n", " value=0.10,\n", " min=0.00,\n", " max=1.00,\n", " step=0.01,\n", " description=\"P('/da/' | zgodny)\",\n", " readout_format=\".2f\",\n", " continuous_update=False,\n", " )\n", " p_niezgodny_slider = widgets.FloatSlider(\n", " value=0.60,\n", " min=0.00,\n", " max=1.00,\n", " step=0.01,\n", " description=\"P('/da/' | niezgodny)\",\n", " readout_format=\".2f\",\n", " continuous_update=False,\n", " )\n", " alpha_slider = widgets.FloatSlider(\n", " value=0.05,\n", " min=0.001,\n", " max=0.20,\n", " step=0.001,\n", " description=\"α\",\n", " readout_format=\".3f\",\n", " continuous_update=False,\n", " )\n", " seed_input = widgets.IntText(value=SEED, description=\"Ziarno\")\n", " correction_checkbox = widgets.Checkbox(\n", " value=False,\n", " description=\"Poprawka Yatesa\",\n", " indent=False,\n", " )\n", "\n", " def symuluj_tabele_2x2(n_per_condition, p_da_zgodny, p_da_niezgodny, alpha, seed, correction):\n", " n_per_condition = int(n_per_condition)\n", " p_da_zgodny = float(p_da_zgodny)\n", " p_da_niezgodny = float(p_da_niezgodny)\n", " alpha = float(alpha)\n", " seed = int(seed)\n", " correction = bool(correction)\n", "\n", " rng_local = np.random.default_rng(seed)\n", "\n", " da_zgodny = int(rng_local.binomial(n=n_per_condition, p=p_da_zgodny))\n", " da_niezgodny = int(rng_local.binomial(n=n_per_condition, p=p_da_niezgodny))\n", "\n", " table = np.array(\n", " [\n", " [n_per_condition - da_zgodny, da_zgodny],\n", " [n_per_condition - da_niezgodny, da_niezgodny],\n", " ]\n", " )\n", "\n", " df_table = pd.DataFrame(\n", " table,\n", " index=[\"Bodziec zgodny\", \"Bodziec niezgodny\"],\n", " columns=[\"Odpowiedź: nie '/da/'\", \"Odpowiedź: '/da/'\"],\n", " )\n", "\n", " chi2, p_value, df, expected = stats.chi2_contingency(table, correction=correction)\n", "\n", " n_total = int(table.sum())\n", " phi = float(np.sqrt(float(chi2) / n_total))\n", "\n", " p_hat_zgodny = da_zgodny / n_per_condition\n", " p_hat_niezgodny = da_niezgodny / n_per_condition\n", "\n", " ci_zgodny = stats.binomtest(da_zgodny, n=n_per_condition).proportion_ci(\n", " confidence_level=0.95, method=\"exact\"\n", " )\n", " ci_niezgodny = stats.binomtest(da_niezgodny, n=n_per_condition).proportion_ci(\n", " confidence_level=0.95, method=\"exact\"\n", " )\n", "\n", " print(df_table)\n", " print(\"\\nchi2 =\", float(chi2), \"df =\", int(df), \"p =\", float(p_value))\n", " print(\"V Craméra / φ =\", phi)\n", "\n", " decyzja = \"Odrzucamy H0\" if float(p_value) < alpha else \"Nie odrzucamy H0\"\n", " print(f\"Decyzja przy α={alpha:.3f}: {decyzja}\")\n", "\n", " print(\"\\nOdsetek odpowiedzi '/da/':\")\n", " print(\n", " f\"- zgodny: {p_hat_zgodny:.3f} (95% PU: [{ci_zgodny.low:.3f}, {ci_zgodny.high:.3f}])\"\n", " )\n", " print(\n", " f\"- niezgodny: {p_hat_niezgodny:.3f} (95% PU: [{ci_niezgodny.low:.3f}, {ci_niezgodny.high:.3f}])\"\n", " )\n", " print(\"Różnica odsetków (niezgodny − zgodny):\", float(p_hat_niezgodny - p_hat_zgodny))\n", "\n", " plt.figure(figsize=(7, 3.8))\n", " props = [p_hat_zgodny, p_hat_niezgodny]\n", " ci_low = [p_hat_zgodny - ci_zgodny.low, p_hat_niezgodny - ci_niezgodny.low]\n", " ci_high = [ci_zgodny.high - p_hat_zgodny, ci_niezgodny.high - p_hat_niezgodny]\n", "\n", " plt.bar([0, 1], props, color=[\"skyblue\", \"steelblue\"], edgecolor=\"white\")\n", " plt.errorbar([0, 1], props, yerr=[ci_low, ci_high], fmt=\"none\", color=\"black\", capsize=4)\n", " plt.xticks([0, 1], [\"zgodny\", \"niezgodny\"])\n", " plt.ylim(0, 1)\n", " plt.title(\"Odsetek odpowiedzi '/da/' (z 95% PU)\")\n", " plt.ylabel(\"odsetek\")\n", " plt.tight_layout()\n", " plt.show()\n", "\n", " out = widgets.interactive_output(\n", " symuluj_tabele_2x2,\n", " {\n", " \"n_per_condition\": n_slider,\n", " \"p_da_zgodny\": p_zgodny_slider,\n", " \"p_da_niezgodny\": p_niezgodny_slider,\n", " \"alpha\": alpha_slider,\n", " \"seed\": seed_input,\n", " \"correction\": correction_checkbox,\n", " },\n", " )\n", "\n", " display(\n", " widgets.VBox(\n", " [\n", " n_slider,\n", " p_zgodny_slider,\n", " p_niezgodny_slider,\n", " alpha_slider,\n", " seed_input,\n", " correction_checkbox,\n", " ]\n", " ),\n", " out,\n", " )\n" ] }, { "cell_type": "markdown", "id": "1ebf9cff", "metadata": {}, "source": [ "### Przykład IV (test niezależności): dane w postaci ramki danych\n", "\n", "W praktyce dane często przechowujemy w ramce danych: jeden wiersz to jedna osoba (albo jedna obserwacja), a kolumny opisują warunek i wynik.\n", "\n", "Czasem jednak mamy zliczenia (tabelę) i chcemy je zamienić na format „długi” (wiersz = obserwacja), np. żeby użyć funkcji, która oczekuje danych w postaci przypadków.\n", "\n", "Załóżmy bardzo prosty eksperyment: dwie grupy wykonują to samo zadanie pamięciowe, ale w jednym warunku w tle jest cisza, a w drugim — hałas. Zapisujemy, czy wynik był poprawny.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "705b78ee", "metadata": {}, "outputs": [], "source": [ "cell_counts = pd.DataFrame(\n", " {\n", " \"Wynik\": [\"Poprawnie\", \"Poprawnie\", \"Niepoprawnie\", \"Niepoprawnie\"],\n", " \"Warunek\": [\"Cisza\", \"Hałas\", \"Cisza\", \"Hałas\"],\n", " \"Liczebność\": [34, 28, 16, 22],\n", " }\n", ")\n", "\n", "print(cell_counts)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "dc905dfb", "metadata": {}, "outputs": [], "source": [ "def counts_to_cases(df_counts: pd.DataFrame, count_col: str = \"Liczebność\") -> pd.DataFrame:\n", " # Zamienia tabelę zliczeń na format przypadków (wiersz = obserwacja).\n", " idx = df_counts.index.repeat(df_counts[count_col])\n", " cases = df_counts.drop(columns=[count_col]).loc[idx]\n", " return cases.reset_index(drop=True)\n", "\n", "\n", "df = counts_to_cases(cell_counts)\n", "print(df.head())\n" ] }, { "cell_type": "code", "execution_count": null, "id": "689a0594", "metadata": {}, "outputs": [], "source": [ "c_table = pd.crosstab(df[\"Wynik\"], df[\"Warunek\"])\n", "print(\"Tabela kontyngencji:\")\n", "print(c_table)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "1a431b2e", "metadata": {}, "outputs": [], "source": [ "observed = c_table.to_numpy()\n", "\n", "chi2_stat, p_value, df, expected = stats.chi2_contingency(observed, correction=False)\n", "\n", "print(\"chi2 =\", float(chi2_stat))\n", "print(\"df =\", int(df))\n", "print(\"p =\", float(p_value))\n", "\n", "n_total = int(observed.sum())\n", "phi = float(np.sqrt(float(chi2_stat) / n_total))\n", "print(\"V Craméra / φ =\", phi)\n", "\n", "expected_df = pd.DataFrame(expected, index=c_table.index, columns=c_table.columns)\n", "print(\"\\nWartości oczekiwane (E):\")\n", "print(expected_df)\n" ] }, { "cell_type": "markdown", "id": "425a0006", "metadata": {}, "source": [ "Jeśli chcesz wypisać więcej informacji, zwykle przydają się:\n", "\n", "- `pd.crosstab` (tabela kontyngencji),\n", "- `chi2_contingency` (test + wartości oczekiwane),\n", "- własne obliczenia (np. wkład komórek do statystyki $\\chi^2$, reszty Pearsona).\n" ] }, { "cell_type": "markdown", "id": "f6be99f0", "metadata": {}, "source": [ "### Przykład V (test niezależności): czas ekspozycji bodźca a poprawność\n", "\n", "W zadaniach percepcyjnych często manipulujemy czasem ekspozycji bodźca. Załóżmy (dla prostoty), że każda osoba wykonała tylko jedną próbę i została przydzielona do jednego z czterech czasów ekspozycji: 50 ms, 100 ms, 150 ms, 200 ms.\n", "\n", "Zmienna $X$: czas ekspozycji (4 kategorie).\n", "\n", "Zmienna $Y$: odpowiedź (poprawnie / niepoprawnie).\n", "\n", "Poniżej są zliczenia w tabeli kontyngencji.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "d02889a4", "metadata": {}, "outputs": [], "source": [ "observed = pd.DataFrame(\n", " {\n", " \"Poprawnie\": [1, 18, 28, 30],\n", " \"Niepoprawnie\": [5, 22, 12, 10],\n", " },\n", " index=[\"50 ms\", \"100 ms\", \"150 ms\", \"200 ms\"],\n", ")\n", "\n", "observed.index.name = \"Czas ekspozycji\"\n", "observed.columns.name = \"Wynik\"\n", "\n", "print(observed)\n" ] }, { "cell_type": "markdown", "id": "a0877aa1", "metadata": {}, "source": [ "Spróbujmy przeprowadzić test niezależności $\\chi^2$.\n", "\n", "Standardowo, najpierw zrobimy to ręcznie (obliczając wartości oczekiwane i statystykę testową ze wzoru).\n", "\n", "Potem sprawdzimy, czy nasz wynik zgadza się z wynikiem zwracanym przez funkcję `chi2_contingency`.\n", "\n", "**Krok po kroku (NHST):**\n", "\n", "1. Hipotezy: $H_0$ — niezależność; $H_1$ — zależność.\n", "2. Poziom istotności: $\\alpha = 0.05$.\n", "3. Statystyka: $\\chi^2 = \\sum \\frac{(O-E)^2}{E}$.\n", "4. Rozkład pod $H_0$: $\\chi^2(\\text{df}=(4-1)(2-1))=\\chi^2(3)$.\n", "5. Region krytyczny: $\\chi^2 \\ge \\chi^2_{1-\\alpha,3}$.\n", "6. p-wartość: $p = P(\\chi^2 \\ge \\chi^2_{\\text{obs}} \\mid H_0)$.\n", "\n", "Dla interpretacji przydaje się też wielkość efektu: $V$ Craméra.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "0ec26992", "metadata": {}, "outputs": [], "source": [ "obs = observed.to_numpy()\n", "\n", "row_totals = obs.sum(axis=1, keepdims=True)\n", "col_totals = obs.sum(axis=0, keepdims=True)\n", "\n", "n_total = int(obs.sum())\n", "expected = row_totals * col_totals / n_total\n", "\n", "components = (obs - expected) ** 2 / expected\n", "chi2_stat = float(components.sum())\n", "\n", "df = (obs.shape[0] - 1) * (obs.shape[1] - 1)\n", "alpha = 0.05\n", "\n", "critical_value = float(stats.chi2.ppf(1 - alpha, df=df))\n", "p_value = float(stats.chi2.sf(chi2_stat, df=df))\n", "\n", "print(\"Wartości oczekiwane (E):\")\n", "print(pd.DataFrame(expected, index=observed.index, columns=observed.columns))\n", "\n", "print(\"\\nChi-kwadrat:\", chi2_stat)\n", "print(\"df =\", int(df))\n", "print(\"Wartość krytyczna:\", critical_value)\n", "print(\"p-wartość:\", p_value)\n", "\n", "cramers_v = float(np.sqrt(chi2_stat / n_total))\n", "print(\"V Craméra:\", cramers_v)\n", "\n", "max_idx = np.unravel_index(np.argmax(components), components.shape)\n", "max_row = observed.index[int(max_idx[0])]\n", "max_col = observed.columns[int(max_idx[1])]\n", "print(f\"Największy wkład do χ² ma komórka: ({max_row}, {max_col})\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "79800353", "metadata": {}, "outputs": [], "source": [ "chi2_lib, p_lib, df_lib, expected_lib = stats.chi2_contingency(obs)\n", "\n", "print(\"chi2_contingency:\")\n", "print(\"chi2 =\", float(chi2_lib), \"df =\", int(df_lib), \"p =\", float(p_lib))\n", "\n", "min_expected = float(expected_lib.min())\n", "print(\"Najmniejsza oczekiwana liczebność:\", min_expected)\n", "\n", "n_total = int(obs.sum())\n", "cramers_v = float(np.sqrt(float(chi2_lib) / n_total))\n", "print(\"V Craméra:\", cramers_v)\n", "\n", "if min_expected < 5:\n", " print(\n", " \"Uwaga: w tabeli są małe wartości oczekiwane (< 5). \"\n", " \"Warto rozważyć symulację p-wartości (Monte Carlo).\"\n", " )\n" ] }, { "cell_type": "markdown", "id": "85c30f7f", "metadata": {}, "source": [ "Jeśli w tabeli pojawiają się małe wartości oczekiwane (np. < 5), aproksymacja rozkładem $\\chi^2$ może być słabsza.\n", "\n", "Jednym z prostych sposobów jest obliczenie p-wartości metodą Monte Carlo, symulując tabele o tych samych sumach brzegowych.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "2a777896", "metadata": {}, "outputs": [], "source": [ "rng = np.random.default_rng(SEED + 13)\n", "\n", "repeats = 10_000\n", "\n", "row_sums = obs.sum(axis=1)\n", "col_sums = obs.sum(axis=0)\n", "\n", "n_total = int(obs.sum())\n", "\n", "col0_total = int(col_sums[0])\n", "expected = (row_sums[:, None] * col_sums[None, :]) / n_total\n", "chi2_obs = float((((obs - expected) ** 2) / expected).sum())\n", "\n", "chi2_sim = np.empty(repeats, dtype=float)\n", "\n", "for i in range(repeats):\n", " remaining_good = col0_total\n", " remaining_total = n_total\n", "\n", " col0_counts = np.zeros_like(row_sums)\n", "\n", " for r in range(len(row_sums) - 1):\n", " draw = rng.hypergeometric(\n", " ngood=remaining_good,\n", " nbad=remaining_total - remaining_good,\n", " nsample=int(row_sums[r]),\n", " )\n", " col0_counts[r] = draw\n", "\n", " remaining_good -= draw\n", " remaining_total -= int(row_sums[r])\n", "\n", " col0_counts[-1] = remaining_good\n", "\n", " simulated = np.column_stack([col0_counts, row_sums - col0_counts])\n", " chi2_sim[i] = (((simulated - expected) ** 2) / expected).sum()\n", "\n", "p_mc = float((chi2_sim >= chi2_obs).mean())\n", "\n", "print(\"chi2_obs =\", chi2_obs)\n", "print(\"p (Monte Carlo) ≈\", p_mc)\n" ] }, { "cell_type": "markdown", "id": "bbc450d9", "metadata": {}, "source": [ "### Wizualizacje danych z tabeli kontyngencji\n", "\n", "Dla testu niezależności $\\chi^2$ przydają się tzw. **reszty Pearsona**:\n", "\n", "$$\n", "R_{ij} = \\frac{O_{ij}-E_{ij}}{\\sqrt{E_{ij}}}.\n", "$$\n", "\n", "Im większa (co do modułu) reszta, tym bardziej „zaskakujący” jest wynik w danej komórce (w porównaniu z tym, czego oczekiwalibyśmy przy niezależności).\n", "\n", "Poniżej pokazujemy dwa proste wykresy:\n", "\n", "1) tabela + mapa cieplna reszt,\n", "2) wykres mozaikowy (kolor = reszty).\n" ] }, { "cell_type": "code", "execution_count": null, "id": "a358d1f1", "metadata": {}, "outputs": [], "source": [ "standardized_residuals = (obs - expected) / np.sqrt(expected)\n", "\n", "residuals_df = pd.DataFrame(\n", " standardized_residuals,\n", " index=observed.index,\n", " columns=observed.columns,\n", ")\n", "\n", "print(residuals_df)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "e3e9542e", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(6, 3.5))\n", "sns.heatmap(\n", " residuals_df,\n", " center=0,\n", " cmap=\"RdBu_r\",\n", " annot=True,\n", " fmt=\".2f\",\n", " cbar_kws={\"label\": \"reszta (standaryzowana)\"},\n", ")\n", "plt.title(\"Reszty Pearsona: gdzie tabela odbiega od niezależności?\")\n", "plt.xlabel(\"Wynik\")\n", "plt.ylabel(\"Czas ekspozycji\")\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "code", "execution_count": null, "id": "2ca4b09a", "metadata": {}, "outputs": [], "source": [ "from statsmodels.graphics.mosaicplot import mosaic\n", "import matplotlib as mpl\n", "import matplotlib.colors as mcolors\n", "import matplotlib.cm as cm\n", "\n", "mosaic_data: dict[tuple[str, str], int] = {}\n", "for row_label in observed.index:\n", " for col_label in observed.columns:\n", " mosaic_data[(row_label, col_label)] = int(observed.loc[row_label, col_label])\n", "\n", "cmap = mpl.colormaps.get_cmap(\"RdBu_r\")\n", "res_min = float(residuals_df.to_numpy().min())\n", "res_max = float(residuals_df.to_numpy().max())\n", "\n", "norm = mcolors.TwoSlopeNorm(vmin=res_min, vcenter=0.0, vmax=res_max)\n", "\n", "\n", "def props(key: tuple[str, str]) -> dict[str, object]:\n", " row_label, col_label = key\n", " r = float(residuals_df.loc[row_label, col_label])\n", " return {\"facecolor\": cmap(norm(r)), \"edgecolor\": \"white\"}\n", "\n", "\n", "fig, ax = plt.subplots(figsize=(8, 4))\n", "mosaic(\n", " mosaic_data,\n", " properties=props,\n", " title=\"Wykres mozaikowy (kolor = reszty Pearsona)\",\n", " gap=0.02,\n", " label_rotation=0,\n", " ax=ax,\n", ")\n", "\n", "sm = cm.ScalarMappable(norm=norm, cmap=cmap)\n", "sm.set_array([])\n", "fig.colorbar(sm, ax=ax, shrink=0.8, label=\"reszta (standaryzowana)\")\n", "\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "b45b8a35", "metadata": {}, "source": [ "## Podsumowanie\n", "\n", "- W testach $\\chi^2$ pracujemy na **zliczeniach** w kategoriach.\n", "- Statystyka testowa to suma $\\sum \\frac{(O-E)^2}{E}$ — im większa, tym większe odejście od $H_0$.\n", "- Ponieważ $\\chi^2 \\ge 0$, testy $\\chi^2$ są praktycznie zawsze **prawostronne**.\n", "- Test zgodności: df zwykle $r-1$ (liczba kategorii minus 1).\n", "- Test niezależności: df = $(r-1)(c-1)$.\n", "- p-wartość mówi o zgodności danych z $H_0$, ale nie mówi „jak duży” jest efekt — dlatego warto raportować też wielkość efektu (np. $w$ Cohena, $V$ Craméra).\n", "- Gdy oczekiwane liczebności są małe, warto rozważyć symulację (Monte Carlo) albo test dokładny (np. Fisher dla tabeli 2×2).\n", "- Do interpretacji przydają się reszty Pearsona: które komórki najbardziej „ciągną” wynik.\n" ] }, { "cell_type": "markdown", "id": "4ada09df", "metadata": {}, "source": [ "## Ćwiczenia\n", "\n", "1) (Test zgodności) W przykładzie II zmień zliczenia na: `Szczęście=24, Złość=18, Strach=20, Smutek=18` (nadal $n=80$). Policz $\\chi^2$, p-wartość, wielkość efektu $w$ i podejmij decyzję przy $\\alpha=0.05$.\n", "\n", "2) (Test zgodności) W przykładzie I sprawdź, jak zmieni się p-wartość, gdy liczba prób wzrośnie dwukrotnie, ale proporcja trafień pozostanie taka sama. Czy p-wartość musi mierzyć „siłę” efektu?\n", "\n", "3) (Test niezależności) W przykładzie III policz wkład każdej komórki do statystyki $\\chi^2$, czyli wartości $(O-E)^2/E$. Która komórka wnosi najwięcej?\n", "\n", "4) (Test niezależności) W przykładzie V zwiększ liczbę osób w warunku `50 ms` o 10 (dodaj 10 obserwacji do tej samej proporcji poprawnie/niepoprawnie). Jak zmieni się $\\chi^2$ i p-wartość? Czy wnioski są wrażliwe na liczebność próby?\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "pygments_lexer": "ipython3" } }, "nbformat": 4, "nbformat_minor": 5 }