Nieklasyczne Metody Analizy Danych

Tablice wielodzielcze: Wnioskowanie dokładne, symulacyjne oraz testy z hipotezami złożonymi

Author

Filip Wawrzacz

1. Wstęp

Analiza tablic wielodzielczych (contingency tables) jest kluczowym elementem statystyki w badaniu zależności między zmiennymi jakościowymi. W klasycznym podejściu często stosuje się test \(\chi^2\) Persona, jednak posiada on ograniczenia przy małych liczebnościach oczekiwanych.

W niniejszym projekcie skupię się na alternatywnych metodach wnioskowania: od testu dokładnego Fishera, poprzez metody symulacyjne (Monte Carlo), aż po konstruowanie specyficznych testów permutacyjnych dla hipotez kierunkowych oraz złożonych (logika “i” / “lub”).

1.1 Biblioteki

Code 
library(tidyverse)
library(stats)
library(broom)

2. Wnioskowanie dla tablic 2x2

Tablica 2x2 to najprostszy przypadek badania zależności, gdzie sprawdzamy np. czy wystąpienie efektu zależy od przynależności do grupy.

2.1 Wnioskowanie dokładne (Test Fishera)

Test dokładny Fishera opiera się na rozkładzie hipergeometrycznym. W przeciwieństwie do testu \(\chi^2\), nie opiera się on na założeniach asymptotycznych, co czyni go idealnym dla małych prób. Zakładając ustalone sumy brzegowe, prawdopodobieństwo otrzymania danej konfiguracji liczb w tablicy \(2 \times 2\) wyraża się wzorem:\[P = \frac{\binom{a+b}{a} \binom{c+d}{c}}{\binom{n}{a+c}}\]

Code 
#| code-fold: show
tablica_2x2 <- matrix(c(8, 2, 3, 9), nrow = 2, 
                      dimnames = list(Grupa = c("Lek", "Placebo"),
                                     Wynik = c("Poprawa", "Brak")))

tidy(fisher.test(tablica_2x2))
# A tibble: 1 × 6
  estimate p.value conf.low conf.high method                         alternative
     <dbl>   <dbl>    <dbl>     <dbl> <chr>                          <chr>      
1     10.4  0.0300     1.19      158. Fisher's Exact Test for Count… two.sided
Code 
x_wykres <- 0:11
gęstość_hyper <- dhyper(x_wykres, 11, 11, 10) 

df_fisher <- data.frame(Sukcesy = x_wykres, P = gęstość_hyper)

ggplot(df_fisher, aes(x = Sukcesy, y = P)) +
  geom_col(fill = kolor_niebieski, alpha = 0.6) +
  geom_col(data = filter(df_fisher, Sukcesy == 8), fill = kolor_czerwony) +
  labs(title = "Rozkład prawdopodobieństwa w teście Fishera",
       subtitle = "Czerwony słupek oznacza naszą zaobserwowaną wartość (8 sukcesów)",
       x = "Liczba sukcesów w grupie 'Lek'", y = "Prawdopodobieństwo") +
  theme(
    panel.background = element_rect(fill = tlo_wykresu),
    plot.background = element_rect(fill = tlo_wykresu),
    text = element_text(color = kolor_tekstu_i_osi),
    axis.text = element_text(color = kolor_tekstu_i_osi)
  )

2.2 Wnioskowanie symulacyjne (Monte Carlo)

Gdy tablica jest większa lub chcemy uniknąć złożoności obliczeniowej rozkładu hipergeometrycznego dla dużych danych, stosujemy symulacje. Generujemy tysiące losowych tablic o tych samych sumach brzegowych i sprawdzamy, jaki odsetek z nich wykazuje silniejszą zależność niż nasza tablica zaobserwowana.

Code 
#| code-fold: show
test_sim <- chisq.test(tablica_2x2, simulate.p.value = TRUE, B = 10000)
test_sim

    Pearson's Chi-squared test with simulated p-value (based on 10000
    replicates)

data:  tablica_2x2
X-squared = 6.6, df = NA, p-value = 0.0296

3. Hipotezy kierunkowe w tablicach wielodzielczych

Standardowe testy (jak \(\chi^2\)) badają ogólną niezależność (hipoteza dwustronna). Często jednak badacz stawia hipotezę kierunkową, np.: “Prawdopodobieństwo sukcesu w Grupie A jest większe niż w Grupie B”.

3.1 Symulacyjny test kierunkowy

W przypadku tablic \(2 \times 2\) możemy badać iloraz szans (Odds Ratio). Jeśli \(OR > 1\), sugeruje to pozytywną zależność. Poniższy kod obrazuje, jak za pomocą symulacji ocenić istotność jednostronną.

Code 
#| code-fold: show
or_obs <- (tablica_2x2[1,1] * tablica_2x2[2,2]) / (tablica_2x2[1,2] * tablica_2x2[2,1])

sim_kierunkowy <- function(tab, n_sim = 5000) {
  successes <- sum(tab[,1])
  total_a <- sum(tab[1,])
  total_b <- sum(tab[2,])
  n_total <- sum(tab)
  
  wyniki_or <- replicate(n_sim, {
    draw_a <- rhyper(1, successes, n_total - successes, total_a)
    draw_b <- successes - draw_a
    (draw_a * (total_b - draw_b)) / ((total_a - draw_a) * draw_b)
  })
  
  p_value <- mean(wyniki_or >= or_obs, na.rm = TRUE)
  return(list(OR_obs = or_obs, p_val = p_value))
}

sim_kierunkowy(tablica_2x2)
$OR_obs
[1] 12

$p_val
[1] 0.0148
Code 
#| code-fold: show
n_sim_plot <- 10000
successes <- sum(tablica_2x2[,1]); total_a <- sum(tablica_2x2[1,]); 
total_b <- sum(tablica_2x2[2,]); n_total <- sum(tablica_2x2)

sim_ors <- replicate(n_sim_plot, {
  d_a <- rhyper(1, successes, n_total - successes, total_a)
  d_b <- successes - d_a
  (d_a * (total_b - d_b)) / ((total_a - d_a) * d_b)
})

df_or <- data.frame(OR = sim_ors[is.finite(sim_ors)])

ggplot(df_or, aes(x = OR)) +
  geom_density(fill = kolor_niebieski, alpha = 0.4, color = kolor_niebieski) +
  geom_vline(xintercept = or_obs, color = kolor_czerwony, linetype = "dashed", size = 1) +
  labs(title = "Symulowany rozkład Ilorazu Szans (OR) pod H0",
       subtitle = "Linia przerywana to zaobserwowany efekt (OR > 1)",
       x = "Odds Ratio", y = "Gęstość") +
  theme(
    panel.background = element_rect(fill = tlo_wykresu),
    plot.background = element_rect(fill = tlo_wykresu),
    text = element_text(color = kolor_tekstu_i_osi),
    axis.text = element_text(color = kolor_tekstu_i_osi)
  )
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

4. Hipotezy złożone (ORAZ / LUB)

W tablicach wielowymiarowych (np. \(2 \times 2 \times 2\)) lub przy wielu tablicach \(2 \times 2\) możemy chcieć sprawdzić hipotezy złożone, np. “Zależność występuje w mieście X LUB w mieście Y”. ### 4.1 Symulacja dla hipotezy złożonej “LUB” Załóżmy badanie w dwóch oddzielnych centrach medycznych. Interesuje nas, czy w którymkolwiek z nich lek zadziałał istotnie.

Code 
centrum_1 <- matrix(c(10, 2, 4, 8), 2)
centrum_2 <- matrix(c(5, 5, 6, 4), 2)

n_iter <- 5000
istotnosc_lub <- replicate(n_iter, {
  s1 <- chisq.test(matrix(rbinom(4, 20, 0.5), 2), simulate.p.value = T)$p.value
  s2 <- chisq.test(matrix(rbinom(4, 20, 0.5), 2), simulate.p.value = T)$p.value
  
  s1 < 0.05 | s2 < 0.05
})

wynik_lub <- mean(istotnosc_lub)

Zjawisko to wiąże się z problemem porównań wielokrotnych – prawdopodobieństwo popełnienia błędu I rodzaju wzrasta, gdy testujemy alternatywę “LUB”. ## 5. Wizualizacja wyników symulacji Poniżej przedstawiam rozkład statystyki \(\chi^2\) uzyskany drogą symulacji dla tablicy o zadanych brzegach, w porównaniu z teoretycznym rozkładem gęstości.

Code 
n_sim <- 5000
statystyki <- replicate(n_sim, {
  m <- r2dtable(1, rowSums(tablica_2x2), colSums(tablica_2x2))[[1]]
  chisq.test(m, correct = FALSE)$statistic
})

df_sim <- data.frame(chi_val = statystyki)

ggplot(df_sim, aes(x = chi_val)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "steelblue", alpha = 0.6) +
  stat_function(fun = dchisq, args = list(df = 1), color = "red", linewidth = 1) +
  labs(title = "Symulowany vs Teoretyczny rozkład Chi-kwadrat (df=1)",
       x = "Wartość statystyki", y = "Gęstość") +
  theme(
    panel.background = element_rect(fill = tlo_wykresu),
    plot.background = element_rect(fill = tlo_wykresu),
    text = element_text(color = kolor_tekstu_i_osi),
    axis.text = element_text(color = kolor_tekstu_i_osi)
  )

6. Podsumowanie

Wnioskowanie dla tablic wielodzielczych przy użyciu metod nieklasycznych pozwala na:

  1. Precyzję przy małych próbach (Fisher).
  2. Elastyczność w definiowaniu hipotez (symulacje kierunkowe).
  3. Kontrolę nad złożonymi strukturami logicznymi badanych zjawisk.

Metody symulacyjne są szczególnie cenne w dobie dużej mocy obliczeniowej, gdyż pozwalają na “zobaczenie” rozkładu statystyki bez konieczności spełniania rygorystycznych założeń testów parametrycznych.

7. Bibliografia

  • Agresti, A. (2013). Categorical Data Analysis. Wiley.
  • Fisher, R. A. (1922). On the interpretation of \(\chi^2\) from contingency tables, and the calculation of P.
  • Wickham, H. et al. (2019). Welcome to the Tidyverse.
  • Dokumentacja R: help(fisher.test), help(r2dtable).