Analiza tablic wielodzielczych i wnioskowanie symulacyjne

Author

Natalia Bąchor

Wstęp

Projekt dotyczy analizy zmiennych jakościowych przy użyciu tablic kontyngencji. Skupia się na metodach obliczeniowych w środowisku R, które zastępują zawodne testy asymptotyczne (jak $\chi^2$) w przypadku małych prób. Takie podejście pozwala na automatyzację i elastyczne badanie nietypowych zależności statystycznych.

Rozwiń kod źródłowy

library(ggplot2)
library(knitr)
set.seed(616)

1. Wnioskowanie dokładne dla tablic 2x2 (Test Fishera)

Rozważamy sytuację małej próby, gdzie klasyczny test niezależności $\chi^2$ z powodów asymptotycznych nie ma zastosowania. Do tego celu używamy Dokładnego Testu Fishera.

Rozwiń kod źródłowy

library(knitr)
x <- c(1,1,1,1, 0,0,0,0)
y <- c(1,1,1,1, 0,0,0,0)

tab_obs <- table(x, y)
dimnames(tab_obs) <- list(X = c("0", "1"), Y = c("0", "1"))
kable(tab_obs, caption = "Tabela obserwowana dla X i Y")

Tabela obserwowana dla X i Y
	0	1
0	4	0
1	0	4

W tabeli na przekątnej mamy wartości równe 4 – pełna zależność, jedynkom w wektorze $x$ odpowiadają jedynki w wektorze $y$.

Test dokładny bada, jak prawdopodobne (z użyciem rozkładu hipergeometrycznego) jest wystąpienie takiej tabeli przy ustalonej liczbie jedynek i zer w rzędach i kolumnach (czyli przy ustalonych sumach brzegowych).

Rozwiń kod źródłowy

test_dokladny <- fisher.test(tab_obs, alternative = "greater")
test_dokladny


    Fisher's Exact Test for Count Data

data:  tab_obs
p-value = 0.01429
alternative hypothesis: true odds ratio is greater than 1
95 percent confidence interval:
 2.003768      Inf
sample estimates:
odds ratio 
       Inf

2. Wnioskowanie symulacyjne (permutacyjne) na poziomie wektorów danych

Aby zasymulować rozkład danej statystyki z tabeli (np. różnicy ilości wystąpień między klasami, czy samego licznika w komórce [1,1]) pod hipotezą o całkowitym braku zależności, wystarczy w pętli losowo mieszać wektor y za pomocą sample(y).

Ten zabieg rozrywa naturalne powiązanie pomiędzy parą obserwacji $(X_i, Y_i)$, narzucając układ w którym sumy brzegowe pozostają identyczne (bo liczby jedynek i zer w $X$ i wymieszanym $Y$ nie ulegają zmianie), a jedyne co się zmienia to przyporządkowanie wewnątrz tabeli.

2.1 Symulacja dla licznika z komórki [1,1]

Budujemy rozkład empiryczny statystyki testowej opartej bezpośrednio o jedną wybraną komórkę tabeli. Jako statystykę wiodącą definiujemy ilość wspólnych wystąpień zer ($X=0$ oraz $Y=0$), co odpowiada komórce tab[1,1].

Rozwiń kod źródłowy

B <- 10000     
T0_komorka <- tab_obs[1,1] 
T_komorka <- numeric(B)

for (i in 1:B) {
  tab_sim <- table(x, sample(y))
  T_komorka[i] <- tab_sim[1,1]
}

table(T_komorka)

T_komorka
   0    1    2    3    4 
 146 2241 5135 2338  140

Aby móc odrzucić hipotezę zerową, musimy sprawdzić, jaka część rozkładu permutacyjnego osiąga lub przekracza wartość z naszej próby oryginalnej ($T_0$).

Rozwiń kod źródłowy

p_val_komorka <- sum(T_komorka >= T0_komorka) / B
cat("Empiryczne p-value dla komórki [1,1] wynosi:", p_val_komorka, "\n")

Empiryczne p-value dla komórki [1,1] wynosi: 0.014

Rozwiń kod źródłowy

df_komorka <- as.data.frame(table(T_komorka))
colnames(df_komorka) <- c("Statystyka", "Czestosc")
df_komorka$Statystyka <- as.numeric(as.character(df_komorka$Statystyka))

ggplot(df_komorka, aes(x = Statystyka, y = Czestosc)) +
  geom_bar(stat = "identity", fill = "#3498db", width = 0.5) +
  geom_vline(xintercept = T0_komorka, color = "darkred", 
             linetype = "dashed", linewidth = 1.2) +
  annotate("text", x = T0_komorka - 0.2, y = max(df_komorka$Czestosc)*0.9, 
           label = paste("Statystyka\nobserwowana T0 =", T0_komorka), 
           color = "darkred", hjust = 1) +
  scale_x_continuous(breaks = min(df_komorka$Statystyka):max(df_komorka$Statystyka)) +
  labs(title = "Rozkład empiryczny wartości komórki [1,1] po permutacji",
       x = "Wartość w komórce [1,1]", y = "Liczność (z symulacji)") +
  theme_minimal()

2.2 Symulacja różnicy zliczeń w kolumnie

W definiowaniu miar uchybień możemy, zdefiniować statystykę w postaci różnicy wartości pomiędzy pierwszą a drugą komórką pierwszej kolumny: $T_0 = tab[1,1] - tab[2,1]$.

Rozwiń kod źródłowy

T0_roznica <- tab_obs[1,1] - tab_obs[2,1]
T_roznica <- numeric(B)

for (i in 1:B) {
  tab_sim <- table(x, sample(y))
  T_roznica[i] <- tab_sim[1,1] - tab_sim[2,1]
}

p_val_roznica <- sum(T_roznica >= T0_roznica) / B
cat("Empiryczne p-value dla różnicy T0 =", T0_roznica, "wynosi:", p_val_roznica, "\n")

Empiryczne p-value dla różnicy T0 = 4 wynosi: 0.0148

Rozwiń kod źródłowy

df_roznica <- as.data.frame(table(T_roznica))
colnames(df_roznica) <- c("Statystyka", "Czestosc")
df_roznica$Statystyka <- as.numeric(as.character(df_roznica$Statystyka))

ggplot(df_roznica, aes(x = Statystyka, y = Czestosc)) +
  geom_bar(stat = "identity", fill = "#2ecc71", width = 0.5) +
  geom_vline(xintercept = T0_roznica, color = "darkred", 
             linetype = "dashed", linewidth = 1.2) +
  annotate("text", x = T0_roznica - 0.5, y = max(df_roznica$Czestosc)*0.9, 
           label = paste("Obserwowana\nróżnica T0 =", T0_roznica), 
           color = "darkred", hjust = 1) +
  scale_x_continuous(breaks = min(df_roznica$Statystyka):max(df_roznica$Statystyka)) +
  labs(title = "Rozkład empiryczny różnicy (tab[1,1] - tab[2,1]) po permutacji",
       x = "Różnica wartości", y = "Liczność (z symulacji)") +
  theme_minimal()

Zależność z prób oryginalnych była tak silna (1. przypadek na 70 możliwych do ułożenia kombinacji tablic!), że obie zaimplementowane metody symulacyjne odrzuciły hipotezę zerową z prawdopodobieństwem $\sim \frac{1}{70} \approx 0.014$.

3. Tablice wielodzielcze: Hipotezy złożone z warunkami logicznymi

Klasyczne testy odpowiadają na ogólne pytanie o zależność w tablicach $R \times C$. Często musimy odpowiedzieć na zapytania złożone. Badamy odsetek sukcesów pomiędzy trzema odrębnymi grupami.

Rozwiń kod źródłowy

dane_wielodzielcze <- matrix(c(30, 10,  
                               20, 20,  
                               15, 25), 
                             nrow = 3, byrow = TRUE)
rownames(dane_wielodzielcze) <- c("Grupa 1", "Grupa 2", "Grupa 3")
colnames(dane_wielodzielcze) <- c("Sukces", "Porazka")
kable(dane_wielodzielcze)

	Sukces	Porazka
Grupa 1	30	10
Grupa 2	20	20
Grupa 3	15	25

Za pomocą algorytmu symulacyjnego chcemy sprawdzić następujące hipotezy bazujące na odsetkach sukcesów ($p_1, p_2, p_3$):

Hipoteza złożona “LUB”: Grupa 1 ma lepsze wyniki od Grupy 2 lub Grupy 3. Wyznacznik $S_{lub} = \max(p_1 - p_2, p_1 - p_3)$.
Hipoteza złożona “I”: Grupa 1 ma lepsze wyniki od Grupy 2 oraz Grupy 3. Wyznacznik $S_{i} = \min(p_1 - p_2, p_1 - p_3)$.

Rozwiń kod źródłowy

get_prop <- function(tab) { tab[,1] / rowSums(tab) }
stat_lub <- function(p) { max(p[1] - p[2], p[1] - p[3]) }
stat_i <- function(p) { min(p[1] - p[2], p[1] - p[3]) }

p_obs <- get_prop(dane_wielodzielcze)
obs_S_lub <- stat_lub(p_obs)
obs_S_i <- stat_i(p_obs)

B_multi <- 10000
sim_tables_multi <- r2dtable(B_multi, rowSums(dane_wielodzielcze), colSums(dane_wielodzielcze))

sim_S_lub <- numeric(B_multi)
sim_S_i <- numeric(B_multi)

for(k in 1:B_multi) {
  p_sim <- get_prop(sim_tables_multi[[k]])
  sim_S_lub[k] <- stat_lub(p_sim)
  sim_S_i[k] <- stat_i(p_sim)
}

p_value_lub <- mean(sim_S_lub >= obs_S_lub)
p_value_i <- mean(sim_S_i >= obs_S_i)

wyniki_zlozone <- data.frame(
  Hipoteza = c("Złożona LUB (max z różnic)", "Złożona I (min z różnic)"),
  Stat_Obserwowana = c(obs_S_lub, obs_S_i),
  Empiryczne_P_value = c(p_value_lub, p_value_i)
)
kable(wyniki_zlozone)

Hipoteza	Stat_Obserwowana	Empiryczne_P_value
Złożona LUB (max z różnic)	0.375	0.0009
Złożona I (min z różnic)	0.250	0.0019

Wizualizacja dla hipotezy złożonej “I”

Rozkład z symulacji Monte Carlo statystyki $S_i$:

Rozwiń kod źródłowy

df_plot_i <- data.frame(S_i = sim_S_i)

ggplot(df_plot_i, aes(x = S_i)) +
  geom_histogram(aes(y = after_stat(density)), bins = 35, 
                 fill = "#9b59b6", color = "white") +
  geom_vline(xintercept = obs_S_i, color = "#e74c3c", 
             linetype = "solid", linewidth = 1.2) +
  annotate("text", x = obs_S_i + 0.05, y = 1.5, 
           label = paste("Obserwowane S_i =", round(obs_S_i, 3), 
                         "\np =", round(p_value_i, 4)), 
           color = "#e74c3c", fontface = "bold") +
  labs(title = 'Rozkład empiryczny różnic min(p1-p2, p1-p3) (Symulacja H0)',
       x = "Różnica proporcji najsłabszego ogniwa", 
       y = "Gęstość") +
  theme_minimal()

Podsumowanie

Prezentowane podejście pokazuje ewolucję wnioskowania analitycznego. Zaczynając od metod dokładnych (kombinatorycznych) dla tabel 2x2, przechodzimy przez metody mieszania wektorów $X$ i $Y$ dla małych wektorów by zasymulować rozkład pod H0 na prostych wskaźnikach własnych, by wreszcie stosując funkcję r2dtable testować elastyczne i wielowymiarowe hipotezy dla struktur o wyższych wymiarach niż tablica dwudzielcza.

Źródła

Majewski, J. Tablice wielodzielcze. Testy niezależności oraz test jednorodności $\chi^2$. Wykład z przedmiotu Statystyka w Medycynie, WFiIS AGH.
[Link do prezentacji (PDF)]
Agresti, A. (2012). Categorical Data Analysis (3rd ed.). John Wiley & Sons.[Link do książki (PDF)]

--- title: "Analiza tablic wielodzielczych i wnioskowanie symulacyjne" author: "Natalia Bąchor" format: html: theme: cosmo toc-title-document: "Spis treści" toc: true toc-location: left code-fold: show code-summary: "Rozwiń kod źródłowy" code-tools: true code-block-bg: "#f8f9fa" code-block-border-left: "#e83e8c" highlight-style: arrow embed-resources: true link-external-newwindow: true execute: echo: true warning: false message: false --- # Wstęp Projekt dotyczy analizy zmiennych jakościowych przy użyciu tablic kontyngencji. Skupia się na metodach obliczeniowych w środowisku R, które zastępują zawodne testy asymptotyczne (jak $\chi^2$) w przypadku małych prób. Takie podejście pozwala na automatyzację i elastyczne badanie nietypowych zależności statystycznych. ```{r} library(ggplot2) library(knitr) set.seed(616) ``` # 1. Wnioskowanie dokładne dla tablic 2x2 (Test Fishera) Rozważamy sytuację małej próby, gdzie klasyczny test niezależności $\chi^2$ z powodów asymptotycznych nie ma zastosowania. Do tego celu używamy **Dokładnego Testu Fishera**. ```{r} library(knitr) x <- c(1,1,1,1, 0,0,0,0) y <- c(1,1,1,1, 0,0,0,0) tab_obs <- table(x, y) dimnames(tab_obs) <- list(X = c("0", "1"), Y = c("0", "1")) kable(tab_obs, caption = "Tabela obserwowana dla X i Y") ``` W tabeli na przekątnej mamy wartości równe 4 – pełna zależność, jedynkom w wektorze $x$ odpowiadają jedynki w wektorze $y$. Test dokładny bada, jak prawdopodobne (z użyciem rozkładu hipergeometrycznego) jest wystąpienie takiej tabeli przy ustalonej liczbie jedynek i zer w rzędach i kolumnach (czyli przy ustalonych sumach brzegowych). ```{r} test_dokladny <- fisher.test(tab_obs, alternative = "greater") test_dokladny ``` # 2. Wnioskowanie symulacyjne (permutacyjne) na poziomie wektorów danych Aby zasymulować rozkład danej statystyki z tabeli (np. różnicy ilości wystąpień między klasami, czy samego licznika w komórce `[1,1]`) pod hipotezą o całkowitym braku zależności, **wystarczy w pętli losowo mieszać wektor y za pomocą `sample(y)`**. Ten zabieg rozrywa naturalne powiązanie pomiędzy parą obserwacji $(X_i, Y_i)$, narzucając układ w którym sumy brzegowe pozostają identyczne (bo liczby jedynek i zer w $X$ i wymieszanym $Y$ nie ulegają zmianie), a jedyne co się zmienia to przyporządkowanie wewnątrz tabeli. ## 2.1 Symulacja dla licznika z komórki \[1,1\] Budujemy rozkład empiryczny statystyki testowej opartej bezpośrednio o jedną wybraną komórkę tabeli. Jako statystykę wiodącą definiujemy ilość wspólnych wystąpień zer ($X=0$ oraz $Y=0$), co odpowiada komórce `tab[1,1]`. ```{r} B <- 10000 T0_komorka <- tab_obs[1,1] T_komorka <- numeric(B) for (i in 1:B) { tab_sim <- table(x, sample(y)) T_komorka[i] <- tab_sim[1,1] } table(T_komorka) ``` Aby móc odrzucić hipotezę zerową, musimy sprawdzić, jaka część rozkładu permutacyjnego osiąga lub przekracza wartość z naszej próby oryginalnej ($T_0$). ```{r} p_val_komorka <- sum(T_komorka >= T0_komorka) / B cat("Empiryczne p-value dla komórki [1,1] wynosi:", p_val_komorka, "\n") ``` ```{r} df_komorka <- as.data.frame(table(T_komorka)) colnames(df_komorka) <- c("Statystyka", "Czestosc") df_komorka$Statystyka <- as.numeric(as.character(df_komorka$Statystyka)) ggplot(df_komorka, aes(x = Statystyka, y = Czestosc)) + geom_bar(stat = "identity", fill = "#3498db", width = 0.5) + geom_vline(xintercept = T0_komorka, color = "darkred", linetype = "dashed", linewidth = 1.2) + annotate("text", x = T0_komorka - 0.2, y = max(df_komorka$Czestosc)*0.9, label = paste("Statystyka\nobserwowana T0 =", T0_komorka), color = "darkred", hjust = 1) + scale_x_continuous(breaks = min(df_komorka$Statystyka):max(df_komorka$Statystyka)) + labs(title = "Rozkład empiryczny wartości komórki [1,1] po permutacji", x = "Wartość w komórce [1,1]", y = "Liczność (z symulacji)") + theme_minimal() ``` ## 2.2 Symulacja różnicy zliczeń w kolumnie W definiowaniu miar uchybień możemy, zdefiniować statystykę w postaci różnicy wartości pomiędzy pierwszą a drugą komórką pierwszej kolumny: $T_0 = tab[1,1] - tab[2,1]$. ```{r} T0_roznica <- tab_obs[1,1] - tab_obs[2,1] T_roznica <- numeric(B) for (i in 1:B) { tab_sim <- table(x, sample(y)) T_roznica[i] <- tab_sim[1,1] - tab_sim[2,1] } p_val_roznica <- sum(T_roznica >= T0_roznica) / B cat("Empiryczne p-value dla różnicy T0 =", T0_roznica, "wynosi:", p_val_roznica, "\n") ``` ```{r} df_roznica <- as.data.frame(table(T_roznica)) colnames(df_roznica) <- c("Statystyka", "Czestosc") df_roznica$Statystyka <- as.numeric(as.character(df_roznica$Statystyka)) ggplot(df_roznica, aes(x = Statystyka, y = Czestosc)) + geom_bar(stat = "identity", fill = "#2ecc71", width = 0.5) + geom_vline(xintercept = T0_roznica, color = "darkred", linetype = "dashed", linewidth = 1.2) + annotate("text", x = T0_roznica - 0.5, y = max(df_roznica$Czestosc)*0.9, label = paste("Obserwowana\nróżnica T0 =", T0_roznica), color = "darkred", hjust = 1) + scale_x_continuous(breaks = min(df_roznica$Statystyka):max(df_roznica$Statystyka)) + labs(title = "Rozkład empiryczny różnicy (tab[1,1] - tab[2,1]) po permutacji", x = "Różnica wartości", y = "Liczność (z symulacji)") + theme_minimal() ``` Zależność z prób oryginalnych była tak silna (1. przypadek na 70 możliwych do ułożenia kombinacji tablic!), że obie zaimplementowane metody symulacyjne odrzuciły hipotezę zerową z prawdopodobieństwem $\sim \frac{1}{70} \approx 0.014$. # 3. Tablice wielodzielcze: Hipotezy złożone z warunkami logicznymi Klasyczne testy odpowiadają na ogólne pytanie o zależność w tablicach $R \times C$. Często musimy odpowiedzieć na zapytania złożone. Badamy odsetek sukcesów pomiędzy **trzema odrębnymi grupami**. ```{r} dane_wielodzielcze <- matrix(c(30, 10, 20, 20, 15, 25), nrow = 3, byrow = TRUE) rownames(dane_wielodzielcze) <- c("Grupa 1", "Grupa 2", "Grupa 3") colnames(dane_wielodzielcze) <- c("Sukces", "Porazka") kable(dane_wielodzielcze) ``` Za pomocą algorytmu symulacyjnego chcemy sprawdzić następujące hipotezy bazujące na odsetkach sukcesów ($p_1, p_2, p_3$): - **Hipoteza złożona "LUB":** Grupa 1 ma lepsze wyniki od Grupy 2 **lub** Grupy 3. Wyznacznik $S_{lub} = \max(p_1 - p_2, p_1 - p_3)$. - **Hipoteza złożona "I":** Grupa 1 ma lepsze wyniki od Grupy 2 **oraz** Grupy 3. Wyznacznik $S_{i} = \min(p_1 - p_2, p_1 - p_3)$. ```{r} get_prop <- function(tab) { tab[,1] / rowSums(tab) } stat_lub <- function(p) { max(p[1] - p[2], p[1] - p[3]) } stat_i <- function(p) { min(p[1] - p[2], p[1] - p[3]) } p_obs <- get_prop(dane_wielodzielcze) obs_S_lub <- stat_lub(p_obs) obs_S_i <- stat_i(p_obs) B_multi <- 10000 sim_tables_multi <- r2dtable(B_multi, rowSums(dane_wielodzielcze), colSums(dane_wielodzielcze)) sim_S_lub <- numeric(B_multi) sim_S_i <- numeric(B_multi) for(k in 1:B_multi) { p_sim <- get_prop(sim_tables_multi[[k]]) sim_S_lub[k] <- stat_lub(p_sim) sim_S_i[k] <- stat_i(p_sim) } p_value_lub <- mean(sim_S_lub >= obs_S_lub) p_value_i <- mean(sim_S_i >= obs_S_i) wyniki_zlozone <- data.frame( Hipoteza = c("Złożona LUB (max z różnic)", "Złożona I (min z różnic)"), Stat_Obserwowana = c(obs_S_lub, obs_S_i), Empiryczne_P_value = c(p_value_lub, p_value_i) ) kable(wyniki_zlozone) ``` ### Wizualizacja dla hipotezy złożonej "I" Rozkład z symulacji Monte Carlo statystyki $S_i$: ```{r} df_plot_i <- data.frame(S_i = sim_S_i) ggplot(df_plot_i, aes(x = S_i)) + geom_histogram(aes(y = after_stat(density)), bins = 35, fill = "#9b59b6", color = "white") + geom_vline(xintercept = obs_S_i, color = "#e74c3c", linetype = "solid", linewidth = 1.2) + annotate("text", x = obs_S_i + 0.05, y = 1.5, label = paste("Obserwowane S_i =", round(obs_S_i, 3), "\np =", round(p_value_i, 4)), color = "#e74c3c", fontface = "bold") + labs(title = 'Rozkład empiryczny różnic min(p1-p2, p1-p3) (Symulacja H0)', x = "Różnica proporcji najsłabszego ogniwa", y = "Gęstość") + theme_minimal() ``` # Podsumowanie Prezentowane podejście pokazuje ewolucję wnioskowania analitycznego. Zaczynając od metod dokładnych (kombinatorycznych) dla tabel 2x2, przechodzimy przez metody mieszania wektorów $X$ i $Y$ dla małych wektorów by zasymulować rozkład pod H0 na prostych wskaźnikach własnych, by wreszcie stosując funkcję `r2dtable` testować elastyczne i wielowymiarowe hipotezy dla struktur o wyższych wymiarach niż tablica dwudzielcza. # Źródła - **Majewski, J.** *Tablice wielodzielcze. Testy niezależności oraz test jednorodności* $\chi^2$. Wykład z przedmiotu Statystyka w Medycynie, WFiIS AGH. <br>\[<a href="https://kompilatory.agh.edu.pl/statwmed/wyklady/WFiIS-STATwMED-4-Tablice-kontyngencji.pdf" target="_blank">Link do prezentacji (PDF)</a>\] - **Agresti, A.** (2012). *Categorical Data Analysis* (3rd ed.). John Wiley & Sons.\[<a href="https://mybiostats.wordpress.com/wp-content/uploads/2015/03/3rd-ed-alan_agresti_categorical_data_analysis.pdf" target="_blank">Link do książki (PDF)</a>\]