Zagadnienia z Wykładu

Antoni Dałek

1. Tablica wielodzielcza

Tablica wielodzielcza służy do przedstawiania łącznego rozkładu dwóch, lub wiekszej ilości zmiennych. Przyjmuje ona zazwyczaj postać macierzy, w której poszczególne komórki określają liczebność lub odsetek badanej populacji dla danej kombinacji wariantów tych zmiennych.

Tablice wielodzielcze są często używane ponieważ:

  • Są łatwe do zrozumienia, także dla osób nierozumiejących bardziej wyszukanych miar
  • Mogą być używane w przypadku zmiennych mierzonych na dowolnym poziomie: nominalnym, porządkowym, interwałowym czy ilorazowym
  • Łatwiej zauważyć związki między zmiennymi, analizując taką tablicę zamiast oddzielnych statystyk
  • Rozwiązują problem braków danych

2. Testy na tablicy 2x2

2.1 “Lady tasting tea”

W latach 20. XX wieku podczas popołudniowej herbaty w Cambridge, pewna pani stwierdziła, że potrafi bezbłędnie rozpoznać, co zostało wlane do filiżanki jako pierwsze – mleko czy herbata. Obecny tam wybitny statystyk, Ronald A. Fisher, postanowił sprawdzić jej twierdzenie w sposób naukowy, co doprowadziło do narodzin nowoczesnego testowania hipotez statystycznych.

Na czym polegał eksperyment?

  • Fisher przygotował 8 filiżanek herbaty z mlekiem
  • W 4 filiżankach najpierw wlał mleko, a potem herbatę
  • W pozostałych 4 filiżankach kolejność była odwrotna
  • Filiżanki podano pani w losowej kolejności. Wiedziała ona, że są dokładnie 4 filiżanki każdego rodzaju, a jej zadaniem było wskazanie tych 4, w których mleko wlano jako pierwsze

2.2 Przygotowanie danych

W eksperymencie mieliśmy 8 filiżanek (4 z mlekiem na dole, 4 z herbatą na dole). Jeśli dama odgadła wszystkie bezbłędnie, tablica na przekątnej będzie miała same czwórki, a poza nią zera.

Dane te zapisujemy w postaci tablicy 2x2:

Code
dane = matrix(c(4, 0, 0, 4), nrow = 2, byrow = TRUE)

colnames(dane) = c("Prawda: Mleko", "Prawda: Herbata")
rownames(dane) = c("Wybór: Mleko", "Wybór: Herbata")
tablica1 = as.table(dane)

kable(tablica1, caption = "Lady tasting tea (Tablica 2x2)")
Lady tasting tea (Tablica 2x2)
Prawda: Mleko Prawda: Herbata
Wybór: Mleko 4 0
Wybór: Herbata 0 4

W ramach wnioskowania dokładnego i symulacyjnego przeprowadzimy testy niezależności. Zbadamy, czy istnieje statystycznie istotny związek między faktyczną zawartością filiżanki a deklaracją badanej osoby.

2.3 Wnioskowanie dokładne (Test Fishera)

Code
# 1. Funkcja obliczająca p-value metodą Fishera 
fisher_p <- function(tab) {
  a <- tab[1,1]; b <- tab[1,2]
  c <- tab[2,1]; d <- tab[2,2]
  n <- a + b + c + d
  
  r1 <- a + b; r2 <- c + d
  c1 <- a + c; c2 <- b + d
  
  prob_obs <- (factorial(r1) * factorial(r2) * factorial(c1) * factorial(c2)) / 
    (factorial(n) * factorial(a) * factorial(b) * factorial(c) * factorial(d))
  
  p_value <- prob_obs + ((factorial(r1) * factorial(r2) * factorial(c1) * factorial(c2)) / 
                           (factorial(n) * factorial(0) * factorial(4) * factorial(4) * factorial(0)))
  
  return(list(P_tablicy = prob_obs, P_value = p_value))
}

# 2. Wyświetlenie wyników
wynik <- fisher_p(tablica1)
cat("Prawdopodobieństwo tej konkretnej tablicy (a=4):", wynik$P_tablicy, "\n")
## Prawdopodobieństwo tej konkretnej tablicy (a=4): 0.01428571
cat("Empiryczny poziom istotności (p-value):", wynik$P_value, "\n")
## Empiryczny poziom istotności (p-value): 0.02857143

Wnioskowanie dokładne opiera się na bezpośrednim wyliczeniu prawdopodobieństwa uzyskania danej konfiguracji tablicy przy ustalonych sumach brzegowych. W przypadku eksperymentu “Lady Tasting Tea” mamy do czynienia z \(n=8\) filiżankami, gdzie sumy wierszy i kolumn są stałe (\(r_1=4, r_2=4, c_1=4, c_2=4\)).

Liczba wszystkich możliwych układów tablicy o takich brzegach wynosi:\[K = \binom{n}{r_1} = \binom{8}{4} = 70\] Prawdopodobieństwo otrzymania konkretnego układu komórek \(a, b, c, d\) wyraża się wzorem:\[p = \frac{r_1! \cdot r_2! \cdot c_1! \cdot c_2!}{n! \cdot a! \cdot b! \cdot c! \cdot d!}\] Hipotezy prezentują się następująco:

  • \(H_0\): Zmienne klasyfikujące są niezależne (Dama zgaduje losowo)
  • \(H_1\): Zmienne klasyfikujące są zależne (Dama posiada zdolność rozróżniania naparu)

Wartość \(p \approx 0.0286\) jest mniejsza od poziomu istotności (\(\alpha = 0.05\)), co daje podstawy do odrzucenia hipotezy zerowej o niezależności zmiennych. Można więc stwierdzić, że deklaracje badanej nie są dziełem przypadku.

2.4 Wnioskowanie symulacyjne

Jako miarę rozbieżności między danymi zaobserwowanymi a oczekiwanymi wykorzystamy statystykę testową \(\chi^2\). Algorytm wylosuje 10 000 tablic z zachowaniem sum brzegowych przy założeniu niezależności zmiennych, obliczy dla każdej z nich wartość tej statystyki, a następnie sprawdzi, jak często pojawia się wynik tak skrajny jak nasz.

Wzór na statystykę testową: \[\chi^2 = \sum_{i=1}^{R} \sum_{j=1}^{C} \frac{(O_{ij} - E_{ij})^2}{E_{ij}}\]

  • \(R\) – liczba wierszy
  • \(C\) – liczba kolumn
  • \(O_{ij}\) – liczebności zaobserwowane w badaniu
  • \(E_{ij}\) – liczebności oczekiwane przy założeniu hipotezy zerowej \(H_0\)
Code
set.seed(123)
test = chisq.test(tablica1, simulate.p.value = TRUE, B = 10000)
print(test)
## 
##  Pearson's Chi-squared test with simulated p-value (based on 10000
##  replicates)
## 
## data:  tablica1
## X-squared = 8, df = NA, p-value = 0.0268

Hipotezy prezentują się następująco:

  • \(H_0\): Zmienne klasyfikujące są niezależne
  • \(H_1\): Zmienne klasyfikujące są zależne

Wartość \(p \approx 0.0268\) jest mniejsza od poziomu istotności (\(\alpha = 0.05\)), co daje podstawy do odrzucenia hipotezy zerowej o niezależności zmiennych. Można więc stwierdzić, że deklaracje badanej nie są dziełem przypadku.

3. Testy na tablicy 3x3

3.1 Przygotowanie danych

Załóżmy, że przeprowadziliśmy badanie kliniczne: zrekrutowaliśmy dokładnie po 50 pacjentów do trzech grup (Lek A, Lek B, Placebo). Sprawdzamy ich stan zdrowia po przyjęciu leku (Pogorszenie, Bez zmian, Poprawa).

Wyniki zapisujemy w postaci tablicy 3x3:

Code
dane_3x3 <- matrix(c( 5, 10, 35, 15, 20, 15, 30, 15,  5), nrow = 3, byrow = TRUE)

rownames(dane_3x3) <- c("Lek A", "Lek B", "Placebo")
colnames(dane_3x3) <- c("Pogorszenie", "Bez zmian", "Poprawa")
tablica_3x3 <- as.table(dane_3x3)

kable(tablica_3x3, caption = "Wyniki badania klinicznego (Tablica 3x3)")
Wyniki badania klinicznego (Tablica 3x3)
Pogorszenie Bez zmian Poprawa
Lek A 5 10 35
Lek B 15 20 15
Placebo 30 15 5

3.2 Test jednorodności

W teście jednorodności weryfikujemy ogólną hipotezę mówiącą o tym, czy rozkłady badanej cechy są identyczne we wszystkich analizowanych populacjach.

Hipotezy prezentują się następująco:

  • \(H_0\): Rozkład zmian stanu zdrowia jest jednorodny (taki sam) w każdej z trzech grup
  • \(H_1\): Rozkłady nie są jednorodne (przynajmniej jedna grupa różni się od pozostałych)

Wykorzystujemy wnioskowanie symulacyjne, losując tablice 10 000 razy.

Code
set.seed(42)
test3 <- chisq.test(tablica_3x3, simulate.p.value = TRUE, B = 10000)
test3
## 
##  Pearson's Chi-squared test with simulated p-value (based on 10000
##  replicates)
## 
## data:  tablica_3x3
## X-squared = 47.788, df = NA, p-value = 9.999e-05
cat("Wartość p :", test3$p.value, "\n")
## Wartość p : 9.999e-05

Wygenerowana wartość \(p = 9.999 × 10^{-5} = 0.00009999\) jest bliska 0. Oznacza to, że przy założeniu prawdziwości \(H_0\), szansa na wylosowanie tak zróżnicowanej tablicy w 10 000 próbach symulacji wynosi praktycznie zero.

W związku z tym, na poziomie istotności \(\alpha = 0.05\), odrzucamy hipotezę zerową \(H_0\) na rzecz hipotezy alternatywnej \(H_1\). Można stwierdzić, że rozkłady stanu zdrowia pacjentów nie są jednorodne, a przynależność do konkretnej grupy (Lek A, Lek B, Placebo) ma statystycznie istotny wpływ na szansę wystąpienia poprawy bądź pogorszenia.

Aby dowiedzieć się, która dokładnie grupa odpowiada za ten wynik, na następnych 2 slajdach przedstawię test symulacyjny weryfikujący hipotezę kierunkową dla konkretnych wycinków tablicy.

3.3 Test symulacyjny

Test ten ocenia istotność statystyczną poprzez fizyczne, wielokrotne symulowanie rozkładu danych przy założeniu prawdziwości hipotezy zerowej (\(H_0\)). Podstawą konstrukcji statystyki testowej jest estymator proporcji: \[\hat{p}_{ij} = \frac{n_{ij}}{n_{i\cdot}}\]

  • \(n_{ij}\) to liczba obserwacji w \(i\)-tej grupie (wierszu) i \(j\)-tej kategorii (kolumnie)
  • \(n_{i\cdot}\) to całkowita liczebność \(i\)-tej grupy (suma \(i\)-tego wiersza)

Różnica zaobserwowanych proporcji stanowi naszą statystykę empiryczną \(T_0\) (np. \(T_0 = \hat{p}_{13} - \hat{p}_{33}\)). W procesie symulacji algorytm wykonuje \(N\) losowych przydziałów obserwacji do grup i dla każdej permutacji wylicza nową statystykę \(T_i\).

Wynikiem testu jest empiryczny poziom istotności, tzw. ASL (Achieved Significance Level), wyznaczany ze wzoru: \[ASL = \frac{\text{card}\{i: T_i \ge T_0\}}{N}\]

  • \(\text{card}\) (kardynalność zbioru) oznacza liczbę permutacji, w których wylosowana statystyka \(T_i\) była równa lub bardziej ekstremalna (w kierunku \(H_1\)) od zaobserwowanej \(T_0\).
  • \(N\) to całkowita liczba przeprowadzonych permutacji

Jeżeli \(ASL < \alpha\) (poziom istotności), oznacza to, że uzyskanie takiej różnicy proporcji przez czysty przypadek jest skrajnie nieprawdopodobne, co stanowi podstawę do odrzucenia hipotezy zerowej na rzecz hipotezy alternatywnej.

3.4 Wynik testu symulacyjnego

Code
# 1. OBLICZENIE STATYSTYKI T_0 (Wiersz 1 vs 3, Kolumna 3)
n_1_dot <- sum(dane_3x3[1, ]) 
n_3_dot <- sum(dane_3x3[3, ])

n_13 <- dane_3x3[1, 3] 
n_33 <- dane_3x3[3, 3] 

p_hat_13 <- n_13 / n_1_dot 
p_hat_33 <- n_33 / n_3_dot 

T_0 <- p_hat_13 - p_hat_33

# 2. SYMULACJA
N <- 10000

sukcesy_suma <- n_13 + n_33
porazki_suma <- (n_1_dot - n_13) + (n_3_dot - n_33)
wspolna_pula <- c(rep(1, sukcesy_suma), rep(0, porazki_suma))

set.seed(123)
T_symulowane <- replicate(N, {
  pula_zmieszana <- sample(wspolna_pula)
  sym_p1 <- mean(pula_zmieszana[1:n_1_dot])
  sym_p3 <- mean(pula_zmieszana[(n_1_dot + 1):(n_1_dot + n_3_dot)])
  return(sym_p1 - sym_p3)
})

# 3. OBLICZENIE ASL
kardynalnosc <- sum(T_symulowane >= T_0)
ASL <- kardynalnosc / N

# 4. WYŚWIETLENIE WYNIKÓW
cat("Wartość poziomu ASL (p-value):", ASL, "\n")
## Wartość poziomu ASL (p-value): 0
cat("Obserwowana różnica proporcji (T_0):", T_0, "\n")
## Obserwowana różnica proporcji (T_0): 0.6

Korzystając z danych, chcemy zweryfikować szczegółową hipotezę kierunkową. Podejrzewamy, że w grupie przyjmującej Lek A odsetek pacjentów wykazujących “Poprawę” jest istotnie wyższy niż w grupie przyjmującej Placebo.

Hipotezy prezentują się następująco:

  • \[H_0: p_{13} = p_{33}\]
  • \[H_1: p_{13} > p_{33}\] Różnica między odsetkiem poprawy w grupie Leku A (70%) a w grupie Placebo (10%) wynosi aż 60 punktów procentowych (\(T_0 = 0.6\)). Podczas symulacji ani razu nie wylosowano tak drastycznej różnicy. Poziom \(ASL \approx 0\), w związku z czym odrzucamy \(H_0\) na rzecz hipotezy alternatywnej \(H_1\). Można więc stwierdzić, że lek A jest istotnie skuteczniejszy od Placebo w indukowaniu poprawy stanu zdrowia.

4. Bibliografia

  • https://pl.wikipedia.org/wiki/Tabela_krzyżowa

  • https://en.wikipedia.org/wiki/Monte_Carlo_method

  • Digital Education Strategies at The G. Raymond Chang School of Continuing Education, “Lady Tasting Tea - Inferential Statistics and Experimental Design” Ryerson University, https://www.youtube.com/watch?v=lgs7d5saFFc

  • https://en.wikipedia.org/wiki/Lady_tasting_tea

  • Kończak Grzegorz, Chmielińska Magdalena. “Zastosowanie metod symulacyjnych w analizie wielowymiarowych tablic wielodzielczych” Uniwersytet Ekonomiczny w Katowicach, (108, 112-117)