{ "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** (*goodness-of-fit test*) — czy rozkład zliczeń w kategoriach jest zgodny z rozkładem oczekiwanym?\n", "- **test niezależności** (*test of independence*) — 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": "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 **kategorialny**, a my zliczamy odpowiedzi lub wystąpienia 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ę." ] }, { "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", "1. Wymyślamy naszą hipotezę alternatywną (ta nas interesuje!).\n", "2. Przyjmujemy hipotezę zerową (i ew. dodatkowe założenia).\n", "3. Decydujemy się na określony poziom istotności $\\alpha$.\n", "4. 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", "5. 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", "6. Losujemy próbę, zbieramy dane, liczymy statystykę.\n", "7. Odrzucamy bądź nie $H_0$." ] }, { "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 matplotlib.colors as mcolors\n", "import matplotlib.cm as cm\n", "import seaborn as sns\n", "from scipy import stats\n", "import cmcrameri.cm as cmc\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" ] }, { "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. 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$. Wrócimy do tego jeszcze za chwilę.\n", "\n", "Teraz wyobraźmy sobie, że rzucamy 100 razy (uczciwą) monetą.\n", "\n", "1. Spodziewamy się, że orzeł wypadnie około 50 razy. Jest to teoretyczna/oczekiwana liczba orłów w naszym eksperymencie.\n", "2. 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", "3. 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", "4. My zamiast dowodzić, pokażemy, że ta zależność działa, korzystając z prostej symulacji." ] }, { "cell_type": "code", "execution_count": null, "id": "1ff29656", "metadata": {}, "outputs": [], "source": [ "rng = np.random.default_rng(SEED + 2)\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):\", standardized_diff.mean())\n", "print(\"Odchylenie standardowe (symulacja):\", standardized_diff.std(ddof=0))\n", "\n", "plt.figure(figsize=(9, 4))\n", "plt.hist(\n", " standardized_diff,\n", " bins=20,\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()" ] }, { "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=20,\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()" ] }, { "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: https://ocw.mit.edu/courses/mathematics/18-443-statistics-for-applications-fall-2003/lecture-notes/lec23.pdf (prawdopodobnie nigdy nie będą Państwo tego dowodzić).\n", "\n", "Przetestujmy tę „prawdę objawioną” na przykładzie dwóch zmiennych.\n", "\n", "(Warto jednak pamiętać, że zwykłe KDE z gaussowskimi jądrami nie odtworzy idealnie kształtu $\\chi^2(1)$ tuż przy zerze: ten rozkład ma nośnik tylko dla wartości nieujemnych i bardzo wysoki *peak* w pobliżu $0$. Dlatego ograniczymy KDE do $x \\ge 0$ i lekko zmniejszymy szerokość pasma, żeby wizualnie było bliższe rozkładowi teoretycznemu.)\n", "\n", "Jak widzimy, rozkład empiryczny $\\chi^2$, który otrzymamy za pomocą symulacji (histogram i niebieska krzywa), odpowiada (w miarę...) rozkładowi teoretycznemu (czerwona krzywa)." ] }, { "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=30, 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(\n", " chi2_stats,\n", " color=\"blue\",\n", " linewidth=2,\n", " bw_adjust=0.7,\n", " cut=0,\n", " clip=(0, np.inf),\n", " label=\"gęstość empiryczna (jądrowa)\",\n", ")\n", "\n", "plt.title(r\"Kwadraty różnic (symulacja) a rozkład $\\chi^2$\")\n", "plt.xlabel(r\"$\\chi^2$\")\n", "plt.ylabel(\"gęstość\")\n", "plt.xlim(0, 5.5)\n", "plt.legend()\n", "plt.tight_layout()\n", "plt.show()" ] }, { "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=False)\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()" ] }, { "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, 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=\"$\\\\chi^2_{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 $\\\\chi^2_(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(\"$\\\\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)" ] }, { "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** (*goodness-of-fit test*) — czy rozkład zliczeń jest zgodny z rozkładem oczekiwanym?\n", "- **test niezależności** (*test of independence*) — 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ą.\n" ] }, { "cell_type": "markdown", "id": "89cb9a67", "metadata": {}, "source": [ "### Przykład I (test zgodności, *goodness-of-fit test*): zadanie dwualternatywnego wyboru wymuszonego\n", "\n", "Wyobraźmy sobie prosty eksperyment z percepcji wzrokowej: na ekranie przez bardzo krótki czas pojawia się słaby bodziec, na przykład strzałka skierowana w lewo albo w prawo, a osoba badana musi wskazać kierunek. Taki paradygmat służy do badania, czy w bodźcu jest w ogóle dość informacji, by odpowiedź była lepsza niż zgadywanie. Z punktu widzenia kognitywistyki pytamy więc, czy układ percepcyjny wydobywa sygnał z szumu, czy uczestnik odpowiada praktycznie losowo.\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 niezależnych obserwacji 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 = 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)" ] }, { "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 = stats.chi2.ppf(1 - alpha, df=df)\n", "p_value = 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\")" ] }, { "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 = observed.sum()\n", "\n", "cohens_w = np.sqrt(chi2_stat / n_total)\n", "\n", "print(\"n =\", n_total)\n", "print(\"Wielkość efektu (w Cohena):\", cohens_w)" ] }, { "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)" ] }, { "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 = observed.sum()\n", "k_obs = 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 = {binom_res.pvalue:.6f}\")\n", "\n", "z = (k_obs - n * p0) / np.sqrt(n * p0 * (1 - p0))\n", "print(f\"Aproksymacja normalna: z = {z:.3f}, z² = {z**2:.3f}\")\n", "print(f\"Statystyka χ² z testu zgodności: {chi2_stat:.3f}\")" ] }, { "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ą." ] }, { "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 = 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 = (sim_chi2 >= chi2_stat).mean()\n", "print(\"Monte Carlo p ≈\", p_mc)" ] }, { "cell_type": "markdown", "id": "a4955ed2", "metadata": {}, "source": [ "### Przykład II (test zgodności, *goodness-of-fit test*): wybór etykiety emocji\n", "\n", "Wyobraźmy sobie ćwiczenie z percepcji społecznej: uczestnicy oglądają tę samą niejednoznaczną twarz, na przykład fotografię balansującą między kilkoma emocjami, i muszą przypisać jej jedną etykietę. Taki eksperyment może badać, czy percepcja emocji jest równomierna, czy też badani mają systematyczny bias ku pewnym kategoriom. Z punktu widzenia kognitywistyki interesuje nas tu, jak niejednoznaczność bodźca siatka naszych ukrytych przekonań wpływają na decyzję percepcyjną.\n", "\n", "Wyobraźmy sobie bardzo proste zadanie: 80 osób widzi tę samą niejednoznaczną 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$ — nieprawda, że wszystkie kategorie mają prawdopodobieństwo $1/4$.\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 = components.sum()\n", "\n", "alpha = 0.05\n", "df = len(labels) - 1\n", "\n", "critical_value = stats.chi2.ppf(1 - alpha, df=df)\n", "p_value = 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 = observed.sum()\n", "print(\"w Cohena =\", 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()" ] }, { "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)" ] }, { "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." ] }, { "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 = 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 = (sim_chi2 >= chi2_stat).mean()\n", "print(\"p (Monte Carlo) ≈\", p_mc)" ] }, { "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(\"$\\\\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)" ] }, { "cell_type": "markdown", "id": "fbebdaa2", "metadata": {}, "source": [ "### Przykład III (test niezależności, *test of independence*): bodźce audiowizualne i kategoria odpowiedzi\n", "\n", "To przykład inspirowany klasycznymi badaniami nad integracją wielozmysłową. Uczestnik jednocześnie słyszy sylabę i widzi ruch ust mówiącej osoby; w warunku zgodnym oba kanały przekazują tę samą informację, a w niezgodnym nie. Kognitywistycznie pytamy, czy to, co widzimy, może systematycznie zmieniać to, co wydaje nam się, że słyszymy. Test niezależności pozwala więc sprawdzić, czy rodzaj bodźca audiowizualnego i kategoria odpowiedzi są ze sobą powiązane.\n", "\n", "W teście niezależności $\\chi^2$ sprawdzamy, czy dwie zmienne kategorialne 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)" ] }, { "cell_type": "markdown", "id": "3b7e38af", "metadata": {}, "source": [ "Zanim policzymy wartości oczekiwane ze wzoru, zatrzymajmy się na intuicji. Hipoteza zerowa mówi, że **warunek bodźcowy** i **kategoria odpowiedzi** są niezależne. To znaczy, że prawdopodobieństwo trafienia do danej komórki powinno być iloczynem prawdopodobieństwa wiersza i prawdopodobieństwa kolumny:\n", "\n", "$$\n", "P(i,j)=P(i)P(j).\n", "$$\n", "\n", "Policzmy jedną komórkę „na piechotę”, np. **(bodziec zgodny, odpowiedź /da/)**. Z danych mamy: $P(\\text{zgodny})=50/100=0.5$ oraz $P(\\text{/da/})=35/100=0.35$. Gdyby zmienne były niezależne, to\n", "\n", "$$\n", "P(\\text{zgodny i /da/}) = 0.5 \\cdot 0.35 = 0.175.\n", "$$\n", "\n", "A ponieważ wszystkich obserwacji jest $N=100$, oczekiwana liczebność w tej komórce wynosi\n", "\n", "$$\n", "E = 100 \\cdot 0.175 = 17.5.\n", "$$\n", "\n", "To jest dokładnie ten sam wynik co ze wzoru\n", "\n", "$$\n", "E_{ij} = N \\cdot P(i)P(j) = \\frac{R_i C_j}{N}.\n", "$$\n", "\n", "Czyli: bierzemy rozkłady brzegowe z danych, zakładamy niezależność, a potem przeliczamy oczekiwane prawdopodobieństwo z powrotem na liczebność. W tej tabeli oba wiersze mają po 50 osób, więc przy niezależności oczekujemy w każdym wierszu tego samego podziału: **32.5** odpowiedzi `/ba/` i **17.5** odpowiedzi `/da/`." ] }, { "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 = 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)" ] }, { "cell_type": "code", "execution_count": null, "id": "371adaa7", "metadata": {}, "outputs": [], "source": [ "chi2_stat = (((observed - expected) ** 2) / expected).sum()\n", "\n", "alpha = 0.05\n", "df = (observed.shape[0] - 1) * (observed.shape[1] - 1)\n", "\n", "critical_value = stats.chi2.ppf(1 - alpha, df=df)\n", "p_value = 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 = observed.sum()\n", "phi = np.sqrt(chi2_stat / n_total)\n", "print(\"V Craméra (dla 2×2 = φ):\", phi)\n", "\n", "n_zgodny = observed[0].sum()\n", "n_niezgodny = observed[1].sum()\n", "\n", "k_da_zgodny = observed[0, 1]\n", "k_da_niezgodny = 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):\", p_hat_niezgodny - p_hat_zgodny)" ] }, { "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 =\", chi2_corr, \"df =\", df_corr, \"p =\", p_corr)\n", "\n", "print(\"\\nBez poprawki:\")\n", "print(\"chi2 =\", chi2_nocorr, \"df =\", df_nocorr, \"p =\", p_nocorr)" ] }, { "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 = observed.sum()\n", "expected = (row_sums[:, None] * col_sums[None, :]) / total\n", "\n", "r1 = row_sums[0]\n", "c1 = 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 = (sim_chi2 >= chi2_stat).mean()\n", "print(\"Zasymulowana wartość p ≈\", p_mc)" ] }, { "cell_type": "markdown", "id": "8258cde7", "metadata": {}, "source": [ "### Interaktywna ilustracja: wiele możliwych badań, liczebność próby i wynik testu (tabela 2×2)\n", "\n", "Ten widżet nie pokazuje jednej tabeli, tylko **wiele możliwych badań** przy tych samych „prawdziwych” prawdopodobieństwach odpowiedzi w obu warunkach.\n", "\n", "Dla wybranych parametrów symulujemy wiele eksperymentów 2×2 i patrzymy:\n", "\n", "- jak zmienia się pojedyncza p-wartość od próby do próby,\n", "- jak często otrzymujemy wynik istotny ($p < \\alpha$),\n", "- jak ten odsetek rośnie wraz z liczebnością próby.\n", "\n", "Jeśli ustawisz takie same prawdopodobieństwa w obu warunkach, odsetek wyników istotnych powinien być bliski $\\alpha$. Jeśli ustawisz różne prawdopodobieństwa, to przy większym `n` będzie zwykle rósł. To jest intuicja stojąca za **mocą testu**, chociaż pełniej wrócimy do niej dopiero w osobnym notatniku.\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=10,\n", " description=\"n/warunek\",\n", " continuous_update=False,\n", " )\n", " p_zgodny_slider = widgets.FloatSlider(\n", " value=0.20,\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.35,\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", " repeats_slider = widgets.IntSlider(\n", " value=1000,\n", " min=200,\n", " max=3000,\n", " step=200,\n", " description=\"Symulacje\",\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 statystyki_chi2_2x2(a, b, c, d, correction=False):\n", " a = np.asarray(a, dtype=float)\n", " b = np.asarray(b, dtype=float)\n", " c = np.asarray(c, dtype=float)\n", " d = np.asarray(d, dtype=float)\n", "\n", " n_total = a + b + c + d\n", " row1 = a + b\n", " row2 = c + d\n", " col1 = a + c\n", " col2 = b + d\n", "\n", " denom = row1 * row2 * col1 * col2\n", " delta = a * d - b * c\n", "\n", " chi2 = np.zeros_like(n_total, dtype=float)\n", " valid = denom > 0\n", "\n", " if correction:\n", " chi2[valid] = n_total[valid] * (np.abs(delta[valid]) - n_total[valid] / 2) ** 2 / denom[valid]\n", " else:\n", " chi2[valid] = n_total[valid] * delta[valid] ** 2 / denom[valid]\n", "\n", " p_values = np.ones_like(chi2)\n", " p_values[valid] = stats.chi2.sf(chi2[valid], df=1)\n", "\n", " phi = np.zeros_like(chi2)\n", " phi[valid] = np.sqrt(chi2[valid] / n_total[valid])\n", "\n", " return chi2, p_values, phi\n", "\n", " def symuluj_wiele_badan_2x2(\n", " n_per_condition,\n", " p_da_zgodny,\n", " p_da_niezgodny,\n", " alpha,\n", " repeats,\n", " seed,\n", " correction,\n", " ):\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", " repeats = int(repeats)\n", " seed = int(seed)\n", " correction = bool(correction)\n", "\n", " rng_local = np.random.default_rng(seed)\n", "\n", " da_zgodny = rng_local.binomial(n=n_per_condition, p=p_da_zgodny, size=repeats)\n", " da_niezgodny = rng_local.binomial(n=n_per_condition, p=p_da_niezgodny, size=repeats)\n", "\n", " a = n_per_condition - da_zgodny\n", " b = da_zgodny\n", " c = n_per_condition - da_niezgodny\n", " d = da_niezgodny\n", "\n", " chi2_values, p_values, phi_values = statystyki_chi2_2x2(a, b, c, d, correction=correction)\n", " diff_props = d / n_per_condition - b / n_per_condition\n", "\n", " true_diff = p_da_niezgodny - p_da_zgodny\n", " significant_rate = (p_values < alpha).mean()\n", "\n", " example_table = np.array([[a[0], b[0]], [c[0], d[0]]])\n", " example_df = pd.DataFrame(\n", " example_table,\n", " index=[\"Bodziec zgodny\", \"Bodziec niezgodny\"],\n", " columns=[\"Odpowiedź: nie '/da/'\", \"Odpowiedź: '/da/'\"],\n", " )\n", "\n", " n_grid = np.arange(10, n_per_condition + 1, 10)\n", " curve_rng = np.random.default_rng(seed + 1)\n", " significant_curve = []\n", "\n", " for n_current in n_grid:\n", " b_curve = curve_rng.binomial(n=n_current, p=p_da_zgodny, size=repeats)\n", " d_curve = curve_rng.binomial(n=n_current, p=p_da_niezgodny, size=repeats)\n", " a_curve = n_current - b_curve\n", " c_curve = n_current - d_curve\n", "\n", " _, p_curve, _ = statystyki_chi2_2x2(a_curve, b_curve, c_curve, d_curve, correction=correction)\n", " significant_curve.append((p_curve < alpha).mean())\n", "\n", " significant_curve = np.array(significant_curve)\n", "\n", " print(\"Ustalony model generujący dane:\")\n", " print(f\"- P('/da/' | zgodny) = {p_da_zgodny:.2f}\")\n", " print(f\"- P('/da/' | niezgodny) = {p_da_niezgodny:.2f}\")\n", " print(f\"- prawdziwa różnica odsetków = {true_diff:.2f}\")\n", "\n", " print(\"\\nJedna przykładowa wylosowana tabela:\")\n", " print(example_df)\n", " print(\n", " f\"\\nDla tej jednej próby: chi2 = {chi2_values[0]:.3f}, p = {p_values[0]:.4f}, φ = {phi_values[0]:.3f}\"\n", " )\n", "\n", " print(\"\\nPodsumowanie po wielu symulowanych badaniach:\")\n", " print(f\"- odsetek prób z p < α = {significant_rate:.3f}\")\n", " print(f\"- mediana p-wartości = {np.median(p_values):.3f}\")\n", " print(f\"- średnia obserwowana różnica odsetków = {diff_props.mean():.3f}\")\n", " if np.isclose(true_diff, 0.0):\n", " print(\"- ponieważ prawdziwa różnica wynosi 0, odsetek p < α powinien być bliski α\")\n", " else:\n", " print(\"- ten odsetek można czytać jako przybliżenie mocy testu dla tych parametrów\")\n", "\n", " fig, axes = plt.subplots(1, 3, figsize=(15, 4))\n", "\n", " axes[0].hist(p_values, bins=np.linspace(0, 1, 21), color=\"lightgrey\", edgecolor=\"white\")\n", " axes[0].axvline(alpha, color=\"crimson\", linestyle=\"--\", linewidth=2, label=\"α\")\n", " axes[0].set_title(\"P-wartości z wielu prób\")\n", " axes[0].set_xlabel(\"p\")\n", " axes[0].set_ylabel(\"liczba symulacji\")\n", " axes[0].legend()\n", "\n", " axes[1].hist(diff_props, bins=20, color=\"skyblue\", edgecolor=\"white\")\n", " axes[1].axvline(true_diff, color=\"navy\", linestyle=\"--\", linewidth=2, label=\"prawdziwa różnica\")\n", " axes[1].set_title(\"Różnica odsetków\\n(niezgodny − zgodny)\")\n", " axes[1].set_xlabel(\"różnica odsetków\")\n", " axes[1].set_ylabel(\"liczba symulacji\")\n", " axes[1].legend()\n", "\n", " axes[2].plot(n_grid, significant_curve, marker=\"o\", color=\"black\")\n", " axes[2].scatter([n_per_condition], [significant_rate], color=\"crimson\", zorder=3, label=\"wybrane n\")\n", " axes[2].axhline(alpha, color=\"grey\", linestyle=\":\", linewidth=2, label=\"poziom α\")\n", " axes[2].set_ylim(0, 1)\n", " axes[2].set_title(\"Jak często p < α,\\ngdy rośnie liczebność\")\n", " axes[2].set_xlabel(\"n na warunek\")\n", " axes[2].set_ylabel(\"częstość p < α\")\n", " axes[2].legend()\n", "\n", " plt.tight_layout()\n", " plt.show()\n", "\n", " out = widgets.interactive_output(\n", " symuluj_wiele_badan_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", " \"repeats\": repeats_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", " repeats_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, *test of independence*): dane w postaci ramki danych\n", "\n", "To bardziej „warsztatowa” wersja pytania kognitywistycznego. Wyobraźmy sobie prosty eksperyment nad pamięcią roboczą: uczestnicy zapamiętują krótki ciąg cyfr, ale część wykonuje zadanie w ciszy, a część przy rozpraszającym hałasie w tle. Interesuje nas, czy warunek środowiskowy wiąże się z rozkładem wyników poprawnych i niepoprawnych. Przy okazji pokazujemy też ważną rzecz praktyczną: te same dane można zapisać jako tabelę zliczeń albo jako „długą” ramkę danych, gdzie każdy wiersz oznacza jedną obserwację.\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." ] }, { "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)" ] }, { "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)" ] }, { "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 =\", chi2_stat)\n", "print(\"df =\", df)\n", "print(\"p =\", p_value)\n", "\n", "n_total = observed.sum()\n", "phi = np.sqrt(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)" ] }, { "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, *test of independence*): czas ekspozycji bodźca a poprawność\n", "\n", "Wyobraźmy sobie mini-eksperyment nad progami percepcyjnymi. Badani widzą bodziec bardzo krótko, na przykład literę, cyfrę albo prosty obiekt, po czym muszą od razu podać odpowiedź. Z punktu widzenia kognitywistyki pytamy, czy wydłużenie czasu ekspozycji wiąże się ze zmianą jakości percepcji, czyli czy dodatkowe dziesiątki milisekund pozwalają zbudować stabilniejszą reprezentację bodźca. Test niezależności sprawdza tu, czy czas ekspozycji i poprawność odpowiedzi są ze sobą powiązane.\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)" ] }, { "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 = obs.sum()\n", "expected = row_totals * col_totals / n_total\n", "\n", "components = (obs - expected) ** 2 / expected\n", "chi2_stat = components.sum()\n", "\n", "df = (obs.shape[0] - 1) * (obs.shape[1] - 1)\n", "alpha = 0.05\n", "\n", "critical_value = stats.chi2.ppf(1 - alpha, df=df)\n", "p_value = 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 =\", df)\n", "print(\"Wartość krytyczna:\", critical_value)\n", "print(\"p-wartość:\", p_value)\n", "\n", "cramers_v = 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[max_idx[0]]\n", "max_col = observed.columns[max_idx[1]]\n", "print(f\"Największy wkład do χ² ma komórka: ({max_row}, {max_col})\")\n" ] }, { "cell_type": "markdown", "id": "f9e54766", "metadata": {}, "source": [ "Warto zatrzymać się na chwilę przy $V$ Craméra. Sama statystyka $\\chi^2$ rośnie nie tylko wtedy, gdy związek jest silniejszy, ale także wtedy, gdy rośnie liczebność próby. Dlatego wygodnie mieć dodatkową miarę, która bardziej opisuje **siłę związku** niż sam rozmiar badania.\n", "\n", "Dla tabeli kontyngencji używamy wzoru\n", "\n", "$$\n", "V = \\sqrt{\\frac{\\chi^2}{N\\,\\min(r-1,c-1)}}.\n", "$$\n", "\n", "W naszym przykładzie tabela ma rozmiar $4\\times 2$, więc $\\min(4-1,2-1)=1$. Dlatego tutaj wzór upraszcza się do\n", "\n", "$$\n", "V = \\sqrt{\\frac{\\chi^2}{N}}.\n", "$$\n", "\n", "To znaczy, że $V$ można czytać jako miarę tego, **jak daleko tabela odchodzi od niezależności po uwzględnieniu liczebności próby**. Gdy $V=0$, nie ma związku; im większe $V$, tym silniejszy związek między zmiennymi." ] }, { "cell_type": "code", "execution_count": null, "id": "0973ee7b", "metadata": {}, "outputs": [], "source": [ "obs_x2 = 2 * obs\n", "\n", "chi2_orig, p_orig, _, _ = stats.chi2_contingency(obs, correction=False)\n", "chi2_x2, p_x2, _, _ = stats.chi2_contingency(obs_x2, correction=False)\n", "\n", "v_orig = np.sqrt(chi2_orig / obs.sum())\n", "v_x2 = np.sqrt(chi2_x2 / obs_x2.sum())\n", "\n", "comparison = pd.DataFrame(\n", " {\n", " \"Chi-kwadrat\": [chi2_orig, chi2_x2],\n", " \"p-wartość\": [p_orig, p_x2],\n", " \"V Craméra\": [v_orig, v_x2],\n", " },\n", " index=[\"oryginalna tabela\", \"ta sama tabela × 2\"],\n", ")\n", "\n", "print(comparison)\n", "print(\"\\nProporcje w tabeli się nie zmieniły, więc V Craméra zostaje takie samo, \"\n", " \"ale chi-kwadrat rośnie, a p-wartość maleje.\")" ] }, { "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 =\", chi2_lib, \"df =\", df_lib, \"p =\", p_lib)\n", "\n", "min_expected = expected_lib.min()\n", "print(\"Najmniejsza oczekiwana liczebność:\", min_expected)\n", "\n", "n_total = obs.sum()\n", "cramers_v = np.sqrt(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 = obs.sum()\n", "\n", "col0_total = col_sums[0]\n", "expected = (row_sums[:, None] * col_sums[None, :]) / n_total\n", "chi2_obs = (((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=row_sums[r],\n", " )\n", " col0_counts[r] = draw\n", "\n", " remaining_good -= draw\n", " remaining_total -= 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 = (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", "abs_residual_max = np.abs(residuals_df.to_numpy()).max()\n", "residual_norm = mcolors.TwoSlopeNorm(vmin=-abs_residual_max, vcenter=0.0, vmax=abs_residual_max)\n", "\n", "print(residuals_df)" ] }, { "cell_type": "code", "execution_count": null, "id": "e3e9542e", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(6.2, 3.8))\n", "sns.heatmap(\n", " residuals_df,\n", " center=0,\n", " cmap=cmc.vik,\n", " norm=residual_norm,\n", " annot=True,\n", " fmt=\".2f\",\n", " linewidths=0.6,\n", " linecolor=\"white\",\n", " cbar_kws={\"label\": \"reszta (standaryzowana)\", \"shrink\": 0.9},\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()" ] }, { "cell_type": "code", "execution_count": null, "id": "2ca4b09a", "metadata": {}, "outputs": [], "source": [ "from statsmodels.graphics.mosaicplot import mosaic\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)] = observed.loc[row_label, col_label]\n", "\n", "cmap = cmc.vik\n", "\n", "\n", "def props(key: tuple[str, str]) -> dict[str, object]:\n", " row_label, col_label = key\n", " r = residuals_df.loc[row_label, col_label]\n", " return {\n", " \"facecolor\": cmap(residual_norm(r)),\n", " \"edgecolor\": \"white\",\n", " \"linewidth\": 1.2,\n", " }\n", "\n", "\n", "fig, ax = plt.subplots(figsize=(8.4, 4.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=residual_norm, cmap=cmap)\n", "sm.set_array([])\n", "fig.colorbar(sm, ax=ax, shrink=0.82, label=\"reszta (standaryzowana)\")\n", "\n", "plt.show()" ] }, { "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 (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 }