Test permutacyji
(1930, R.A. Fisher)
Znany również jako test randomizacji, jest nieparametrycznym testem statystycznym używanym do określenia, czy istnieje znacząca różnica między dwiema lub więcej grupami.
Problem dwóch grup:
Rozważmy:
\(Z=(Z_1, \dots, Z_n)\) - i.d.d. z rozkładu \(F\),
\(Y=(Y_1, \dots, Y_m)\) - i.i.d. z rozkładu \(G\).
Zakładamy, że próby \(X\) i \(Y\) są niezależne.
Testujemy \[H_0: F\leq G\qquad \text{vs}\qquad H_1: F>G.\] Statystyka testowa: \(\hat{\theta}=\bar{Z}-\bar{Y}\).
Intuicyjne, duże wartości \(\hat{\theta}\) będą świadczyć przeciwko \(H_0\). Co to znaczy duże?
Zdefiniujmy achived significant level \(ACL=P_{H_0}\left(\hat{\theta}^{\ast}\geq \hat{\theta}\right)\), gdzie \(\hat{\theta}^{\ast}\) jest generowane przy prawdziwości hipotezy zerowej.
Małe wartości \(ACF\) świadczą przeciwko hipotezie zerowej.
Tradycyjne podjeście - test t-Studenta - wariant 3
Jeżeli założymy, że \(F=N(m_z, \sigma^2)\) oraz \(G=N(m_y, \sigma^2)\), to wtedy testujemy \[H_0: m_z\leq m_y\qquad \text{vs}\qquad H_1: m_z>m_y.\] Wtedy \(\hat{\theta}=\bar{Z}-\bar{Y}\sim N\left(0, \sigma^2\left(\frac{1}{n}+\frac{1}{m}\right)\right)\).
Przypomnijmy, że statystyka testowa w tym przypadku jest następującej postaci: \[U(Y,Z)=\frac{\sqrt{\frac{mn}{(m+n)}}\left(\bar{Z}-\bar{Y}\right)}{\sqrt{\frac{1}{m+n-2}\left(\sum_{i=1}^n(Z_i-\bar{Z})^2+\sum_{j=1}^m(Y_j-\bar{Y})^2\right)}}\sim t_{m+n-2}.\] Wobec tego, test postaci \[\phi(Y,Z)=1, \quad U(Y,Z)>c_{\alpha},\] gdzie \(c_{\alpha}=F^{-1}_{t_{m+n-2}}(1-\alpha)\) jest testem jednostajnie najmocniejszym wśród testów nieobciążonych.
P- wartość: \[p(Y,Z)=1-F_{t_{n+m-2}}\left(U(Y,Z)\right).\]
Test permutacyjny Fishera
W przepadku, gdy nie możemy założyć, że dane pochodzą z rozkładu normalnego lub mają niejednorodną wariancję możemy testu permutacyjnego Fishera.
Idea:
Skoro zmienny \(Z,Y\) mają ten sam rozkład (\(H_0: F=G\)), to możemy je ze sobą połączyć (stworzyć jedną próbę \(X\) o długości \(m+n\)), ignorując informację o przynależności do grup.
Następnie, losujemy bez powtórzeń \(m\) obserwacji z próby \(X\) i traktujemy je jako reprezentację pierwszej grupy, a pozostałe \(n\) obserwacji przypisujemy do drugiej grupy.
Obliczamy różnicę średnich \(\hat{\theta}^{\ast}\) między nowymi grupami utworzonymi w wyniku losowania.
Kroki 2 i 3 powtarzamy \(B\) razy, tworząc rozkład wynikowych różnic \(\hat{\theta}^{\ast}\). Następnie sprawdzamy, czy pierwotna różnica \(\hat{\theta}\) znajduje się poza \(95\%\) przedziałem ufności tego rozkładu. Jeśli tak, to odrzucamy hipotezę zerową.
Dla testu prawostronnego (\(H_0: F\leq G\), przeciwko \(H_1: F>G\)) \[\hat{ASL}_{perm}=\#\{\hat{\theta}^{\ast}(b)\geq \hat{\theta} \}/B.\] Dla testu lewostronnego (\(H_0: F\geq G\), przeciwko \(H_1: F<G\)) \[\hat{ASL}_{perm}=\#\{\hat{\theta}^{\ast}(b)\leq \hat{\theta} \}/B.\] Dla testu obustronnego (\(H_0: F= G\), przeciwko \(H_1: F\neq G\)) \[\hat{ASL}_{perm}=\#\left\{\left|\hat{\theta}^{\ast}(b)\right|\geq \left|\hat{\theta}\right| \right\}/B.\]
Zalety testu permutacyjnego:
Brak założenia normalności;
Brak założenia o jednorodności wariancji;
Użyteczny przy małych próbach, gdy ciężko sprawdzić założenie normalności;
Wady testu permutacyjnego:
Test permutacyjny może nie działać poprawnie, gdy struktura danych między grupami jest mocno niesymetryczna;
Wymagana jest duża moc obliczeniowa, zwłaszcza dla dużych próbek.
Przykład - Badanie skuteczności nowego suplementu
W małym miasteczku naukowcy postanowili zbadać skuteczność nowego suplementu diety, który miał zwiększać poziom energii. W badaniu wzięło udział \(55\) uczestników, którzy zostali podzieleni na dwie grupy: jedna grupa (\(30\) osób) przyjmowała suplement, a druga grupa (\(25\) osób) otrzymywała placebo.
Po \(4\) tygodniach badania, uczestnicy ocenili swój poziom energii w skali od 1 do 100. Naukowcy postanowili zastosować test t-Studenta, aby sprawdzić, czy średni poziom energii w grupie przyjmującej suplement nie jest większy od grupy placebo. Dodatkowo, chcieli przeprowadzić test permutacyjny, aby ocenić, czy różnice były statystycznie istotne, niezależnie od założeń o normalności rozkładu. W ten sposób zyskali pełniejszy obraz skuteczności suplementu.
Załaduj dane z pliku Baza_random i przeprowadź
testy.
library(readxl)
dane=read.csv("/home/karolina/Desktop/Dydaktyka/Testowanie hipotez statystycznych/Testy_randomizacji/Baza_random.csv")
Z=dane$lek
Y=na.omit(dane$placeo)
# Test t-studenta
m=length(Y)
n=length(Z)
test_statistic<-sqrt((m*n)/(n+m))*((mean(Z)-mean(Y))/sqrt((sum((Z-mean(Z))^2)+sum((Y-mean(Y))^2))/(n+m-2)))
p<-1-pt(test_statistic,n+m-2)
p## [1] 0.843722# Sprawdzamy za pomocą gotowej funkcji
t_test_result <- t.test(Z, Y, var.equal = TRUE, alternative="greater")
t_test_result##
## Two Sample t-test
##
## data: Z and Y
## t = -1.0196, df = 53, p-value = 0.8437
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
## -8.868161 Inf
## sample estimates:
## mean of x mean of y
## 50.68587 54.04253# Obliczamy różnicę średnich w danych oryginalnych
theta_obs <- mean(Z) - mean(Y)
# Łączymy dane w jedną próbkę
X <- c(Y, Z)
nx <- length(X)
B=999
theta_star <- numeric(B)
for (i in 1:B) {
permuted <- sample(X)
theta_star[i] <- mean(permuted[1:n]) - mean(permuted[(n+1):length(X)])
}
p_value <- mean(theta_star >= theta_obs)
p_value## [1] 0.8388388Zadanie
Wygeneruj \(20\) zmiennych losowych \(Z\sim N(0,5)\) i \(25\) zmiennych z rozkładu \(Y\sim N(d, 5)\), gdzie \(d\in[-3,3]\). Wyznacz moc empiryczną dla obustronnego testu t-Studenta i porównaj go z mocą empiryczną dla obustronnego testu permutacyjnego.
Następnie, zmień wariancję zmiennej losowej \(Z\sim N(0,1)\), jak to wpłynie na wyznaczone moce empiryczne testów?
Bootstrap
(1979, Bradley Efron)
I.i.d. bootstrap
Bootstrap jest metodą resamplingową, która służy do oszacowania rozkładu pewnej statystyki bez założeń o jej teoretycznym rozkładzie. W i.i.d. bootstrapie zakładamy, że dane są niezależne i mają taki sam rozkład, co pozwala na losowanie z powtórzeniami bezpośrednio z oryginalnej próbki.
Bootstrap jest często używamy do:
oszacowania wariancji statystyki,
tworzenia przedziałów ufności,
wyznaczenia rozkładu poszukiwanej statystyki.
Algorytm:
Załóżmy, że mamy próbę i.i.d. \(X=(X_1, \dots, X_n)\) pochodzącą z nieznanego rozkładu \(F\).
Tworzymy \(B\) prób bootstrapowych \(X^{\ast b}=(X_1^{\ast}, \dots, X_n^{\ast})\) poprzez losowanie z powtórzeniami obserwacji z wyjściowej próby. Prawdopodobieństwo wylosowania każdej obserwacji jest równe.
Dla każdej próby bootstrapowej obliczamy statystykę \(\hat{\theta}^{\ast b}\) np. średnią.
Wyznaczamy rozkład bootstrapowy dla statystyki na podstawie wartości \(\left\{\hat{\theta}^{\ast 1}, \dots, \hat{\theta}^{\ast B}\right\}\).
Uwaga: Bootstrap nie tworzy nowych obserwacji!
Przykład
Poniższy przykład demonstruje zastosowanie i.i.d. bootstrap do oszacowania empirycznego rozkładu średniej, wariancji bootstrapowej oraz przedziału ufności dla średniej na symulowanych danych.
set.seed(42)
n <- 50
B <- 1000
#original_sample <- rnorm(n, mean = 10, sd = 5)
original_sample <- rexp(n,5)
original_mean <- mean(original_sample)
bootstrap_means <- numeric(B)
for (b in 1:B) {
bootstrap_sample <- sample(original_sample, size = n, replace = TRUE)
bootstrap_means[b] <- mean(bootstrap_sample)
}
original_mean## [1] 0.2270545## [1] 0.2301901bootstrap_variance <- var(bootstrap_means)
confidence_interval <- quantile(bootstrap_means, c(0.025, 0.975))
bootstrap_variance## [1] 0.001634959## 2.5% 97.5%
## 0.1604564 0.3129517hist(bootstrap_means, breaks = 30, main = "Empiryczny rozkład średniej bootstrapowej",
xlab = "Średnia bootstrapowa", border = "blue", col = "lightblue")
abline(v = original_mean, col = "red", lwd = 2)
abline(v = confidence_interval, col = "darkgreen", lwd = 2, lty = 2)
Testy bootstrapowe
Problem dwóch prób
Niech
\(Z=(Z_1, \dots, Z_n)\) - i.d.d. z rozkładu \(F\),
\(Y=(Y_1, \dots, Y_m)\) - i.i.d. z rozkładu \(G\).
Zakładamy, że próby \(X\) i \(Y\) są niezależne.
Testujemy \[H_0: F\leq G\qquad \text{vs}\qquad H_1: F>G.\]
Podobnie jak wcześniej możemy użyć statystyki \(t(X)=\bar{Z}-\bar{Y}\) i \(ASL=P_{H_0}\left(t(X^{\ast}\geq t(X))\right)\).
Przy prawdziwości hipotezy zerowej \(X^{\ast}\) ma pewien rozkład, załóżmy \(F_0\). Jakie jest \(F_0\)?
W przypadku bootstrapu
Wykorzystujemy zasadę plug-in i zastępujemy nieznaną dystrybuantę \(F_0\) dystrybuantą empiryczną \(\hat{F}_0\), która przypisuje masę \(\frac{1}{m+n}\) każdej obserwacji.
Algorytm:
Podobnie jak w przypadku testów permutacyjnych łączymy próby \(Z\) i \(Y\) w jedną próbę o długości \(m+n\).
Tworzymy \(B\) prób bootstrapowych losując zmienne zgodnie z rozkładem \(\hat{F}_0\), czyli z powtórzeniami. Pierwszych \(n\) obserwacji uznajemy za pierwszą próbę \(Z^{\ast}\), pozostałe \(m\) za próbę \(Y^{\ast}\).
Dla każdej próby bootstrapowej obliczamy statystykę testową \(t(X^{\ast b}=\bar{Z}-\bar{Y}\).
Obliczamy odpowiednią p-wartość (ALS):
Dla testu prawostronnego (\(H_0: F\leq G\), przeciwko \(H_1: F>G\)) \[\hat{ASL}_{boot}=\#\left\{t \right(X^{\ast b}\left)\geq t(X \right\}/B.\]
Dla testu lewostronnego (\(H_0: F\geq G\), przeciwko \(H_1: F<G\)) \[\hat{ASL}_{perm}=\#\left\{t \right(X^{\ast b}\left)\leq t(X \right\}/B.\]
Dla testu obustronnego (\(H_0: F= G\), przeciwko \(H_1: F\neq G\)) \[\hat{ASL}_{perm}=\#\left\{\left|t\right(X^{\ast b}\left)\right|\geq \left|t(X)\right| \right\}/B.\] Zalety testów bootstrapowych:
prosta metoda,
może być zastosowany do wielu ogólnych problemów,
daje dobre, asymptotyczne wyniki.
Wady testów bootstrapowych:
warto podkreślić, że w tym laboratorium zakładamy, że dane są i.i.d.,
naiwne stosowanie bootstrapu nie zawsze daje asymptotycznie poprawne wyniki,
złożony obliczeniowo.
Przykład - Badanie skuteczności nowego suplementu
Zastosuj test bootstrapowy dla przykładu badania skuteczności suplementu (opisanego wcześniej).
library(readxl)
dane=read.csv("/home/karolina/Desktop/Dydaktyka/Testowanie hipotez statystycznych/Testy_randomizacji/Baza_random.csv")
Z=dane$lek
Y=na.omit(dane$placeo)
n=length(Z)
m=length(Y)
# Obliczamy różnicę średnich w danych oryginalnych
theta_obs <- mean(Z) - mean(Y)
# Łączymy dane w jedną próbkę
X <- c(Y, Z)
nx <- length(X)
B=999
theta_star <- numeric(B)
for (i in 1:B) {
permuted <- sample(X, size = nx, replace = TRUE)
theta_star[i] <- mean(permuted[1:n]) - mean(permuted[(n+1):length(X)])
}
p_value <- mean(theta_star >= theta_obs)
p_value## [1] 0.8478478Testy bootstrapowe dla testowania równości średnich
Załóżmy że mamy \(Z=(Z_1, \dots, Z_n)\sim N(m_z, \sigma_z^2)\) oraz \(Y=(Y_1, \dots, Y_n)\sim N(m_y, \sigma_y^2)\) i chcemy przeprowadzić następujący test
\[H_0: m_z= m_y\qquad \text{vs}\qquad H_1: m_z\neq m_y.\] Zauważmy, że zmienne pochodzą z rozkładu normalnego ale mają niejednorodna wariancję.
Wtedy możemy użyć następujący algorytm:
Niech \(\tilde{Z}_i=Z_i-\bar{Z}-\bar{X}\), \(X=(Z_1, \dots, Z_n, Y_1, \dots, Y_m)\) oraz \(\tilde{Y}_j=Y_j-\bar{Y}-\bar{X}\), gdzie \(i=1, \dots, n\), \(j=1, \dots, m\) oraz \(\bar{X}=\frac{1}{m+n}\sum_{i=1}^{m+n}X_i\).
Następnie tworzymy próby bootstrapowe \(Z^{\ast b}\) z \(\tilde{Z}\) oraz \(Y^{\ast b}\) z \(\tilde{Y}\).
Obliczamy statystykę testową \[T(X^{\ast b})=\frac{\bar{Z}^{\ast}-\bar{Y}^{\ast}}{\sqrt{\sigma_1^{2\ast}/n+\sigma_2^{2\ast}/m}},\]
Obliczamy odpowiednią p-wartość, tutaj: \[\hat{ASL}_{boot}=\#\left\{\left|T(X^{\ast b})\right|\leq \left|T(X)\right|\right\}/B,\] gdzie \(T(X)=\frac{\bar{Z}-\bar{Y}}{\sqrt{\hat{\sigma}^{2}_1/n+\hat{\sigma}^{2}_2/m}}\), \(\hat{\sigma}_1^{2}=\frac{1}{n-1}\sum_{i=1}^n\left(Z_i-\bar{Z}\right)^2\) oraz \(\hat{\sigma}_2^{2}=\frac{1}{m-1}\sum_{i=1}^m\left(Y_i-\bar{Y}\right)^2\).
Przykład - Badanie skuteczności nowego suplementu
Porównaj test bootstrapowy (opisany powyżej) i test t-Studenta o równości średniej pomiędzy grupami dla przykładu badania skuteczności suplementu (opisanego wcześniej).
library(readxl)
dane=read.csv("/home/karolina/Desktop/Dydaktyka/Testowanie hipotez statystycznych/Testy_randomizacji/Baza_random.csv")
Z=dane$lek
Y=na.omit(dane$placeo)
set.seed(43)
# Test t-studenta
m=length(Y)
n=length(Z)
test_statistic<-sqrt((m*n)/(n+m))*((mean(Z)-mean(Y))/sqrt((sum((Z-mean(Z))^2)+sum((Y-mean(Y))^2))/(n+m-2)))
p<-2*(1-pt(abs(test_statistic),n+m-2))
p## [1] 0.312556# Sprawdzamy za pomocą gotowej funkcji
t_test_result <- t.test(Z, Y, var.equal = TRUE, alternative="two.sided")
t_test_result##
## Two Sample t-test
##
## data: Z and Y
## t = -1.0196, df = 53, p-value = 0.3126
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -9.959946 3.246614
## sample estimates:
## mean of x mean of y
## 50.68587 54.04253# Obliczamy T(X)
sigma1 <- sum((Z - mean(Z))^2) / (n - 1)
sigma2 <- sum((Y - mean(Y))^2) / (m - 1)
theta_obs <- (mean(Z) - mean(Y)) / sqrt(sigma1 / n + sigma2 / m)
X <- c(Z, Y)
Z_tilde <- Z - mean(Z) + mean(X)
Y_tilde <- Y - mean(Y) + mean(X)
B=999
theta_star <- numeric(B)
for (i in 1:B) {
Z_star <- sample(Z_tilde, size = n, replace = TRUE)
Y_star <- sample(Y_tilde, size = m, replace = TRUE)
sigma1_star <- sum((Z_star - mean(Z_star))^2) / (n - 1)
sigma2_star <- sum((Y_star - mean(Y_star))^2) / (m - 1)
theta_star[i] <- (mean(Z_star) - mean(Y_star)) / sqrt(sigma1_star / n + sigma2_star / m)
}
p_value <- mean(abs(theta_star) >= abs(theta_obs))
p_value## [1] 0.3033033Zadanie
Wygeneruj \(20\) zmiennych losowych \(Z\sim Gamma(2,1)\) i \(25\) zmiennych z rozkładu \(Y\sim Gamma(d, 1)\), gdzie \(d\in[1,3]\). Wyznacz moc empiryczną dla lewostronnego testu bootstrapowego i porównaj go z mocą empiryczną dla testu permutacyjnego.