Labolatorium 7

Karolina Marek

2024-08-06

Test Wilcoxona

Jest to nieparametryczna wersja testu t-Studenta dla wariantu 1 i 2.

Przypomnijmy, że dla testu t-Studenta mieliśmy \(X_1, \dots, X_n\sim N(m, \sigma^2)\). W wariancie \(1\) testowaliśmy \(H_0: m\leq m_0\), przeciwko \(H_1: m> m_0\). Dla wariantu \(2\) natomiast \(H_0: m = m_0\), przeciwko \(H_1: m\neq m_0\).

Teraz mamy \(X_1, \dots, X_n\) próbę prostą o dystrybuancie \(F_X(t)=F(t)-\Delta(t)\), gdzie:

  • \(F\) jest nieznaną, ciągłą dystrybuantą, symetryczną wokół \(0\),

-\(\Delta\) jest funkcją stałego znaku.

Będą nas interesować dwa warianty:

  1. \(H_0: \Delta \leq 0\), przeciwko \(H_1: \Delta>0\),

  2. \(H_0: \Delta \leq 0\), przeciwko \(H_1: \Delta>0\),

Niech \((|R_1|, \dots, |R_n|)\) będzie wektorem rang dla \((|X_1|, \dots, |X_n|)\).

Oznaczmy również: \[\bar{R}_{|X_+|}=\frac{1}{n}\sum_{i=1}^n|R_i|\mathbb{I}(X_i>0)\quad\text{ oraz }\quad\bar{R}_{|X_-|}=\frac{1}{n}\sum_{i=1}^n|R_i|\mathbb{I}(X_i<0)\] Na przykład: Niech \(X=(-6, 7, -1, 4, -5, 3)\), wtedy \[Rang_{+}(X)=\] oraz \[Rang_{-}(X)=\]

Po co założenie symetrii?

Przy prawdziwości hipotezy zerowej, rozkład jest symetryczny wokół zera, więc gdy obłożymy wartością bezwzględną obserwacji, to ujemne obserwacje będą miały ten sam rozkład co dodatnie obserwacje.

Odwzorowanie \(\left(Rang_{+}(|X_1|, \dots, |X_n|), Rang_{-}(|X_1|, \dots, |X_n|)\right)\) jest maksymalnym niezmiennikiem względem grupy przekształceń \[G=\left\{g: \mathbb{R}^n\rightarrow \mathbb{R}: g(\mathbf{x})=\left(h(X_1), \dots, h(X_n)\right), \quad h\text{ - ciągła, ściśle rosnąca i niepatrzysta}\right\}.\] Test Wilcoxona odrzuci:

  1. \(H_0: \Delta\leq 0\), gdy \(\bar{R}_{|X_+|}-\bar{R}_{|X_-|}\) będzie duże,

  2. \(H_0: \Delta = 0\), gdy \(\left|\bar{R}_{|X_+|}-\bar{R}_{|X_-|}\right|\) będzie duże.

Zestandaryzowana statystyka testowa Wilcoxona przy prawdziwości hipotezy zerowej: \[W=\frac{\bar{R}_{|X_+|}-\bar{R}_{|X_-|}}{\sqrt{V_{\Delta=0}\left(\bar{R}_{|X_+|}-\bar{R}_{|X_-|}\right)}}= \frac{\bar{R}_{|X_+|}-\frac{n+1}{4}}{\sqrt{\frac{(n+1)(2n+1)}{24n}}}\xrightarrow[n\rightarrow \infty]{d} N(0,1).\] Zatem, dla \(H_0: \Delta \leq 0\), przeciwko \(H_1: \Delta>0\) mamy test postaci \(\phi(W)=1\) dla \(W>c_{\alpha}\), gdzie \(c_{\alpha}=F_{N(0,1)}^{-1}(1-\alpha)\) oraz \(p(W)=1-F_{N(0,1)}(W)\).

Natomiast dla \(H_0: \Delta = 0\), przeciwko \(H_1: \Delta\neq 0\) mamy test postaci \(\phi(W)=1\) dla \(|W|>c_{\alpha}\), gdzie \(c_{\alpha}=F_{N(0,1)}^{-1}(1-\alpha/2)\) oraz \(p(W)=2\left(1-F_{N(0,1)}(|W|)\right)\).

Zadanie 1

\(30\) samochodów użyto do zbadania dwóch rodzajów paliwa pod kątem spalania. Samochody zatankowały \(10\) litrów paliwa typu \(1\) i jeździły po torze do momentu aż paliwo się skończyło. Zanotowano liczę przejechanych kilometrów. Później te same samochody zatankowały paliwo typu \(2\) i powtórzono eksperyment. Na podstawie zebranych wyników przetestuj czy te dwa typy paliwa różnią się między sobą bez zakładania normalności różnic. Dane w pliku Baza 2.1.

library(readxl)
Baza<- read_excel("Baza 2.1.xlsx")
X<-Baza$"Typ 1"
Y<-Baza$"Typ 2"
Z<-X-Y
n<-length(Z)
R<-rank(abs(Z),ties.method = "random")
M <- sum(R[Z > 0])
M=M/n
W<-(M-(n+1)/4)/sqrt((n+1)*(2*n+1)/(24*n))
p<-2*(1-pnorm(abs(W),0,1))
p
## [1] 0.5170479

Można również użyć gotowej funkcji wilcox.test() (funkcja podrasowana - z poprawką na ciągłość).

wilcox.test(Z, alternative = "two.sided")
## 
##  Wilcoxon signed rank exact test
## 
## data:  Z
## V = 201, p-value = 0.5291
## alternative hypothesis: true location is not equal to 0

Obliczenie p-wartości testu Wilcoxona metodą MC

  1. Generujemy mc (np. \(mc=1000\)) prób prostych \(U_1, \dots, U_n\sim U[-1,1]\) (może być dowolny rozkład symetryczny wokół \(0\)).

  2. Dla każdej próby z kroku \(1\) wyznaczamy wektor \(\left(|R_1|\mathbb{I}(X_1>0), \dots, |R_n|\mathbb{I}(X_n>0) \right)\).

  3. Dla każdego wektora wyznaczonego w korku \(2\) obliczamy wartość statystyki testowej Wilcoxona \(W_1, \dots, W_{mc}\).

  4. Dla \(H_0: \Delta \leq 0\), przeciwko \(H_1: \Delta>0\) mamy \[\hat{p}=\frac{1}{mc}\sum_{i=1}^{mc} \mathbb{I}(W_i\geq W).\] Dla \(H_0: \Delta = 0\), przeciwko \(H_1: \Delta\neq 0\) mamy \[\hat{p}=\frac{1}{mc}\sum_{i=1}^{mc} \mathbb{I}(|W_i|\geq |W|).\]

Zadanie 1 - cd.

Obliczyć dokładną p-wartość używając metody MC.

set.seed(123)  
MC <- 100000  
WMC <- numeric(MC)

for(i in 1:MC){
  ZMC <- runif(n, -1, 1)
  RMC <- rank(abs(ZMC))
  MMC <- sum(RMC[ZMC > 0]) / n
  WMC[i] <- (MMC - (n + 1) / 4) / sqrt((n + 1) * (2 * n + 1) / (24 * n))
}
pMC <- mean(abs(WMC) >= abs(W))

pMC
## [1] 0.52786

Test znków

Nieparametryczny test dla kwantyla dowolnego rzędu.

Niech \(X_1, \dots, X_n\) będzie próbą prostą z rozkładu o nieznanej ciągłej dystrybuancie \(F\).

Niech \(x_p\) będzie kwantylem rzędu \(p\in (0,1)\) z powyższego rozkładu.

Ponieważ zakładamy, że \(F\) jest ciągła, to \(P(X_i<x_p)=p\) oraz \(P(X_i>x_p)=1-p\).

Interesują nas dwa testy:

  1. \(H_0: x_p\leq m\), przeciwko \(H_1: x_p> m\),

  2. \(H_0: x_p = m\), przeciwko \(H_1: x_p\neq m\).

Wartość \(m\) jest znana i ustalona. Wobec tego, możemy rozważyć przesunięte dane
\(\mathbf{Z}=(Z_1, \dots, Z_n)=(X_1-m,\dots, X_n-m)\). Wtedy \(z_p\) będzie oznaczać kwantyl rzędu \(p\) dla próby \(\mathbf{Z}\). Wówczas, problem testowania sprowadza się do:

  1. \(H_0: z_p\leq 0\), przeciwko \(H_1: z_p> 0\),

  2. \(H_0: z_p = 0\), przeciwko \(H_1: z_p\neq 0\).

Na wykładzie pokazano, że odwzorowania:

\(Rang_+\left(|Z_1|, \dots, |Z_n|\right)=\left(|r_1|\mathbb{I}(Z_1>0), \dots, |r_n|\mathbb{I}(Z_n>0)\right),\)

\(Rang_-\left(|Z_1|, \dots, |Z_n|\right)=\left(|r_1|\mathbb{I}(Z_1<0), \dots, |r_n|\mathbb{I}(Z_n<0)\right).\)

są maksymalnym niezmiennikiem względem grupy przekształceń

\(G=\left\{g: \mathbb{R}^n\rightarrow \mathbb{R}^n:g(\mathbf{x})=\left(h(X_1), \dots, h(X_n)\right), \quad h \text{ - ciagła, ściśle rosnąca i nieparzysta.}\right\}\)

Oznaczmy odpowiednio frakcje dodatnich i ujemnych obserwacji: \[P_+=\frac{1}{n}\sum_{i=1}^n\mathbb{I}(Z_i>0)\quad \text{oraz}\quad P_-=\frac{1}{n}\sum_{i=1}^n\mathbb{I}(Z_i<0).\]

Test znaków odrzuci:

  1. \(H_0: z_p\leq 0\), gdy \(P_+-P_-\) jest duże,

  2. \(H_0: z_p = 0\), gdy \(|P_+-P_-|\) jest duże.

Statystyka testowa to zestandaryzowana różnica \(P_+-P_-\) przy prawdziwości hipotezy zerowej ma następujący rozkład: \[S=\sqrt{n}\frac{P_+-(1-p)}{\sqrt{p(1-p)}}\xrightarrow[n\rightarrow \infty]{d} N(0,1).\]

Zatem:

  1. Dla testu \(H_0: z_p\leq 0\), przeciwko \(H_1: z_p> 0\), mamy \(\phi(S)=1\) dla \(S>c_{\alpha}\), gdzie \(c_{\alpha}=F^{-1}_{N(0,1)}(1-\alpha)\). Asymptotyczna p-wartość wynosi \(p(S)=1-F_{N(0,1)}(S).\)

  2. Dla testu \(H_0: z_p= 0\), przeciwko \(H_1: z_p\neq 0\), mamy \(\phi(S)=1\) dla \(|S|>c_{\alpha}\), gdzie \(c_{\alpha}=F^{-1}_{N(0,1)}(1-\alpha/2)\). Asymptotyczna p-wartość wynosi \(p(S)=2\left(1-F_{N(0,1)}(|S|)\right).\)

Zadanie 2 W pewnej firmie produkującej żarówki wylosowano \(134\) egzemplarzy i sprawdzono ich czas świecenia (w miesiącach z dokładnością do \(1\) miejsca po przecinku). Przetestować hipotezę że mniej niż \(25\)% produkowanych żarówek świeci dłużej niż \(12\) miesięcy. Dane w pliku Baza 7.1

library(readxl)
Baza<- read_excel("Baza 7.1.xlsx")
Z<-Baza$"Czas"
plot(ecdf(Z))

Musimy zadbać o to, aby hipoteza zerowa była zbiorem domkniętym, zatem \(H_0: X_{0.75} \geq 12\), przeciwko \(H_1: X_{0.75} < 12\).

Po wycentrowaniu (\(Z' = X - 12\)) otrzymujemy:
\(H_0: Z'_{0.75} \geq 0\), przeciwko \(H_1: Z'_{0.75} < 0\).

Ponieważ test znaków ma w hipotezie zerowej znak mniejszości, to weźmy \(Z = -(X - 12)\), wówczas:
\(H_0: Z_{0.75} \leq 0\), przeciwko \(H_1: Z_{0.75} > 0\).

n<-length(Z)
rz<-0.75
K<-12

Z=-(Z-K)

M <- sum(Z > 0) / n

S<-sqrt(n)*(M-(1-rz))/sqrt(rz*(1-rz))

p<- (1-pnorm(S,0,1))
M
## [1] 0.7014925
p
## [1] 0

Możemy również użyć godowej funkcji binom.test, która domyślnie przeprowadza test dla mediany, ale można zmienić dla dowolnego kwantylu - argument p.

binom.test(sum(Z > 0), n, p=0.75, alternative = "greater" )
## 
##  Exact binomial test
## 
## data:  sum(Z > 0) and n
## number of successes = 94, number of trials = 134, p-value = 0.9166
## alternative hypothesis: true probability of success is greater than 0.75
## 95 percent confidence interval:
##  0.6297845 1.0000000
## sample estimates:
## probability of success 
##              0.7014925

Obliczanie p-wartości testu znaków za pomocą metody MC

  1. Generujemy mc (np. \(mc=1000\)) prób prostych \(U_1, \dots, U_n\sim U[-p,1-p]\).

  2. Dla każdej próby z kroku \(1\) wyznaczamy statystykę testu znaków otrzymując próbę prostą \(S_1, \dots, S_{mc}\).

  3. Dla \(H_0: z_p \leq 0\), przeciwko \(H_1: z_p>0\) mamy \[\hat{p}=\frac{1}{mc}\sum_{i=1}^{mc} \mathbb{I}(S_i\geq S).\] Dla \(H_0: z_p = 0\), przeciwko \(H_1: z_p\neq 0\) mamy \[\hat{p}=\frac{1}{mc}\sum_{i=1}^{mc} \mathbb{I}(|S_i|\geq |S|).\]

Zadanie 2 - cd.

Obliczmy dokładną p-wartość wykorzystując MC.

MC=1000
ZMC <- matrix(runif(MC * n, -rz, 1-rz), nrow = MC, ncol = n)
MMC <- rowMeans(ZMC > 0)
SMC <- sqrt(n) * (MMC - (1 - rz)) / sqrt(rz * (1 - rz))
pMC <- mean(SMC >= S)
pMC
## [1] 0

Zadanie

Wygeneruj \(n\in[10,50]\) obserwacji z rozkładu \(N(0,1)\). Przeprowadź test \(H_0: x_{0.5}=0.2\), przeciwko \(H_1: x_{0.5}\neq 0.2\). Porównaj na wykresie empiryczną moc testu, gdzie p wartość obliczana jest na dwa sposoby (asymptotycznie i metodą MC).