Testy permutacyjne stanowią fundament nowoczesnego podejścia do wnioskowania statystycznego, szczególnie gdy dane nie spełniają rygorystycznych założeń o rozkładzie normalnym. Cała logika opiera się na prostym mechanizmie: jeśli grupy nie różnią się od siebie, to przypisanie konkretnej obserwacji do danej grupy jest czysto przypadkowe. Poprzez wielokrotne, losowe przemieszanie etykiet, budujemy empiryczny rozkład statystyki, który pozwala precyzyjnie wyznaczyć poziom istotności nawet dla bardzo małych prób.
Proces testowania można najłatwiej prześledzić, budując algorytm Monte Carlo od podstaw. Wykorzystując dane eksperymentalne, sprawdzamy, jak prawdopodobne jest uzyskanie zaobserwowanej różnicy średnich w świecie, w którym rządzi wyłącznie przypadek.
x1 <- c(4, 6, 8, 9, 11, 11, 13, 16)
x2 <- c(3, 4, 4, 5, 5, 6, 10)
n1 <- length(x1)
combined <- c(x1, x2)
srx <- mean(combined)
T0 <- mean(x1) - mean(x2)
set.seed(123)
N_sim <- 10000
T_perm <- numeric(N_sim)
for (i in 1:N_sim) {
srs1 <- mean(sample(combined, n1))
srs2 <- (srx * length(combined) - srs1 * n1) / (length(combined) - n1)
T_perm[i] <- srs1 - srs2
}
# Poziom istotności (p-value)
ASL <- sum(T_perm >= T0) / N_sim
cat("Uzyskane p-value:", ASL)Uzyskane p-value: 0.0124Interpretacja wyniku: Uzyskana wartość \(p\)-value wynosząca 0.0124 (zależnie od konkretnego przebiegu symulacji) sugeruje, że przy tradycyjnym progu istotności \(\alpha = 0,05\) nie mamy podstaw do odrzucenia hipotezy zerowej. Oznacza to, że zaobserwowana różnica między grupami \(x1\) i \(x2\) może być dziełem przypadku i nie jest wystarczająco silna, by uznać ją za statystycznie znaczącą w tej wielkości próby.
ggplot(data.frame(T = T_perm), aes(x = T)) +
geom_histogram(bins = 30, fill = "steelblue", color = "white", alpha = 0.8) +
geom_vline(xintercept = T0, color = "red", size = 1.2, linetype = "dashed") +
labs(title = "Rozkład permutacyjny", x = "Różnica średnich", y = "Częstość") +
theme_minimal()
Wizualne przedstawienie wyników pomaga zrozumieć, gdzie na tle tysięcy losowych permutacji znajduje się nasza rzeczywista różnica (zaznaczona czerwoną linią).
W metodach Monte Carlo liczba permutacji bezpośrednio wpływa na precyzję wyniku. Poniżej sprawdzamy, jak wynik zmienia się przy małej liczbie losowań w porównaniu do dużej próby symulacyjnej.
set.seed(123)
p_small <- sum(replicate(100, mean(sample(combined, n1))) >= mean(x1)) / 100
p_large <- sum(replicate(10000, mean(sample(combined, n1))) >= mean(x1)) / 10000
data.frame(
"Liczba_symulacji" = c(100, 10000),
"Wynik_p_value" = c(p_small, p_large)
) %>% kable(caption = "Porównanie precyzji p-value w zależności od liczby losowań")| Liczba_symulacji | Wynik_p_value |
|---|---|
| 100 | 0.0200 |
| 10000 | 0.0123 |
Dodatkowo możemy zaobserwować proces stabilizacji wyniku na wykresie zbieżności. Pokazuje on, że po osiągnięciu odpowiedniej liczby permutacji, wartość \(p\)-value przestaje gwałtownie oscylować.
p_val_convergence <- cumsum(T_perm >= T0) / (1:N_sim)
ggplot(data.frame(iter = 1:N_sim, p_val = p_val_convergence), aes(x = iter, y = p_val)) +
geom_line(color = "darkgreen") +
geom_hline(yintercept = ASL, linetype = "dashed", color = "red") +
labs(title = "Stabilizacja wyniku p-value w czasie symulacji",
x = "Liczba wykonanych permutacji", y = "Skumulowane p-value") +
theme_minimal()
Zasada działania permutacji sprawia, że wybór konkretnej miary (np. różnica średnich vs. ich iloraz) często prowadzi do identycznych wniosków statystycznych. Wynika to z faktu, że permutacje zachowują monotoniczne relacje między danymi, co czyni te metody niezwykle odpornymi na specyficzne transformacje statystyk.
x1_sim <- rnorm(11, 8, 3)
x2_sim <- rnorm(16, 6, 3)
x_comb <- c(x1_sim, x2_sim)
T0_diff <- mean(x1_sim) - mean(x2_sim)
T0_ratio <- mean(x1_sim) / mean(x2_sim)
res <- data.frame(diff = numeric(N_sim), ratio = numeric(N_sim))
for (i in 1:N_sim) {
xs <- sample(x_comb)
res$diff[i] <- mean(xs[1:11]) - mean(xs[12:27])
res$ratio[i] <- mean(xs[1:11]) / mean(xs[12:27])
}
cat("p-value dla różnicy średnich:", sum(res$diff >= T0_diff)/N_sim,
"\np-value dla ilorazu średnich:", sum(res$ratio >= T0_ratio)/N_sim)p-value dla różnicy średnich: 0.0095
p-value dla ilorazu średnich: 0.0095Powyższe porównanie pokazuje kluczową cechę testów permutacyjnych – niezależnie od tego, czy naszą statystyką testową jest różnica średnich, czy ich iloraz, otrzymujemy identyczny (lub niemal identyczny) wynik \(p\)-value. Potwierdza to teorię, że testy te bazują na wzajemnym położeniu obserwacji względem siebie (permutacjach), a nie na parametrach rozkładu, co czyni je niezwykle odpornymi na wybór konkretnej miary matematycznej.
#Testowanie alternatywnych statystyk: Różnica median Wielką zaletą metod permutacyjnych jest możliwość testowania dowolnej statystyki, np. mediany, która jest bardziej odporna na wartości odstające niż średnia arytmetyczna.
T0_med <- median(x1) - median(x2)
T_perm_med <- numeric(N_sim)
for (i in 1:N_sim) {
shuffled <- sample(combined)
T_perm_med[i] <- median(shuffled[1:n1]) - median(shuffled[(n1+1):length(combined)])
}
ASL_med <- sum(T_perm_med >= T0_med) / N_sim
cat("Zaobserwowana różnica median:", T0_med, "\nUzyskane p-value dla mediany:", ASL_med)Zaobserwowana różnica median: 5
Uzyskane p-value dla mediany: 0.0355Nowoczesne biblioteki w środowisku R oferują zoptymalizowane algorytmy, które automatyzują proces testowania i pozwalają na analizę bardziej złożonych zależności.
Gdy celem analizy jest sprawdzenie, czy wyniki pojedynczej grupy odbiegają od przyjętego standardu (\(\mu\)), testy permutacyjne pozwalają uniknąć błędów typowych dla klasycznego testu t-Studenta, szczególnie przy małych zbiorach danych. Poniżej sprawdzamy hipotezę dla \(\mu = 8\)
test_1p <- oneSamplePermutationTest(x1, mu = 8, alternative = "greater")
summary_table <- data.frame(
"Parametr" = c("Statystyka (Suma)", "Hipoteza zerowa (Mu)", "Wartość p-value", "Zastosowana metoda"),
"Wartość" = c(test_1p$statistic, test_1p$null.value, round(test_1p$p.value, 4), test_1p$method)
)
kable(summary_table, align = "lc")| Parametr | Wartość |
|---|---|
| Statystyka (Suma) | 14 |
| Hipoteza zerowa (Mu) | 8 |
| Wartość p-value | 0.139 |
| Zastosowana metoda | One-Sample Permutation Test |
(Based on Sampling
Permutation Distribution
5000 Times) |Wykorzystując pakiet EnvStats, sprawdziliśmy, czy średnia grupy \(x1\) jest istotnie większa od założonej wartości \(8\). Tabela wyników pokazuje, że mimo iż średnia z próby jest wyższa, test permutacyjny pozwala ocenić, czy ta nadwyżka jest trwała. Jest to podejście znacznie bezpieczniejsze niż klasyczny test \(t\), ponieważ nie zakładamy tutaj, że dane pochodzą z rozkładu normalnego, co przy tak małej liczbie obserwacji (8 osób) byłoby ryzykowne.
W sytuacjach, gdy porównujemy więcej niż dwie grupy, stosujemy podejście oparte na ogólnych testach permutacyjnych. Pozwala to na wykrycie istotnych różnic między grupami bez rygorystycznych założeń parametrycznych.
oneway_test(mpg ~ factor(cyl), data = mtcars, distribution = approximate(nresample = 9999))
Approximative K-Sample Fisher-Pitman Permutation Test
data: mpg by factor(cyl) (4, 6, 8)
chi-squared = 22.706, p-value < 1e-04kruskal_wynik <- kruskal.test(mpg ~ factor(cyl), data = mtcars)
data.frame(
"Miernik" = c("Statystyka Chi-kwadrat", "Stopnie swobody", "Wartość p-value"),
"Wynik" = c(round(kruskal_wynik$statistic, 2), kruskal_wynik$parameter, round(kruskal_wynik$p.value, 5))
) %>% kable(caption = "Weryfikacja różnic między wieloma grupami")| Miernik | Wynik | |
|---|---|---|
| Kruskal-Wallis chi-squared | Statystyka Chi-kwadrat | 25.75 |
| df | Stopnie swobody | 2.00 |
| Wartość p-value | 0.00 |
Podsumowanie analizy wielu grup: Zastosowanie ogólnego testu permutacyjnego (oneway_test) oraz testu Kruskala-Wallisa dla danych mtcars pozwoliło na weryfikację wpływu liczby cylindrów na spalanie. Bardzo niskie \(p\)-value jednoznacznie wskazuje, że liczba cylindrów ma statystycznie istotny wpływ na osiągi, co potwierdza skuteczność metod nieparametrycznych w wykrywaniu zależności w strukturach wielogrupowych.
Metody permutacyjne są nieocenione wszędzie tam, gdzie tradycyjne testy zawodzą z powodu małej próby lub nietypowego rozkładu. Wykorzystanie mocy obliczeniowej do symulacji rozkładów daje badaczowi pewność wyników, która jest odporna na wartości odstające.