Labolatorium 6

Karolina Marek

2024-08-06

Rangi

Rozważmy nieparametryczne problemy testowania, które będą niezmiennicze względem grupy ciągłych i ściele rosnących (niekoniecznie liniowych) przekształceń.

Przypomnijmy z wykładu.

Niech \((X_1, \dots, X_n)\) będzie próbą prostą o ciągłej dystrybuancie \(F\) oraz niech \[G=\{g: \mathbb{R}^n\rightarrow \mathbb{R}^n: g(\mathbf{X})=\left(h(x_1), \dots, h(x_n)\right), \text{ h - ciągła i ściśle rosnąca funckja}\}\] będzie grupą przekształceń.

Rangą i-tej obserwacji wektora \((x_1, \dots, x_n)\) nazywamy liczbę współrzędnych tego wektora mniejszych lub równych \(x_i\), tzn: \[r_i(x_1, \dots, x_n)=\sum_{j=1}^n\mathbb{I}_{(-\infty, x_i]}(x_j).\]

Przykład: \(Rang(9,6,1,4,5)=\)

Odwzorowanie Rang jest maksymalnym niezmiennikiem względem grupy przekształceń G (pokazane na wykładzie).

Zadanie 1

Wygenerować \(20\) obserwacji \(X_1, \dots, X_{20}\) z rozkładu \(N(0,1)\). Narysować punkty \((X_i, R_i)\).

set.seed(321)
n <- 20
X <- rnorm(n)
R <- numeric(n)  

for (i in 1:n) {
  R[i] <- sum(X <= X[i])  
}
plot(X,R)

Zadanie 2

Wygenerować \(20\) obserwacji \(X_1, \dots, X_{20}\) z rozkładu Poissona \(P(5)\). Narysować punkty \((X_i, R_i)\).

set.seed(321)
n <- 20
X <- rpois(n,5)
R <- numeric(n)  

for (i in 1:n) {
  R[i] <- sum(X <= X[i])  
}
X
##  [1] 9 9 3 3 4 4 5 4 5 7 5 4 7 2 5 3 6 4 4 6
plot(X,R)

Zadanie 2 - cd

Narysować wykres jeszcze raz uwzględniając poprawki na powtarzające się rangi.

Dodajemy do rang pewną losowość aby uniknąć powtórzeń.

set.seed(321)
n <- 20
X <- rpois(n,5)
R <- numeric(n)

for (i in 1:n) {
  R[i] <- sum(X <= X[i])
}

U <- runif(n, 0, 1)
C <- R + U

RC <- numeric(n)
for (i in 1:n) {
  RC[i] <- sum(C <= C[i])
}
plot(X,RC)
points(X,R,pch=3,col="red")

Możemy również użyć gotowej funkcji rank().

R<-rank(X,ties.method = "max")
R
##  [1] 20 20  4  4 10 10 14 10 14 18 14 10 18  1 14  4 16 10 10 16
RC<-rank(X,ties.method = "random")
RC
##  [1] 20 19  2  4  9  8 11 10 12 18 13  6 17  1 14  3 15  5  7 16

Test Manna- Whitneya

Przypomnimy, test t-Studenta (wariant 3 i 4) mieliśmy: \(X_1, \dots, X_m\sim N(m_x, \sigma^2)\) oraz \((Y_1, \dots, Y_n)\sim(N(m_y, \sigma^2))\) i testowaliśmy:

(wariant 3) \(H_0: m_y\leq m_x\), przeciwko \(H_1: m_y> m_x\),

(wariant 4) \(H_0: m_y = m_x\), przeciwko \(H_1: m_y\neq m_x\).

Co gdy mamy sytuację, że nie możemy założyć normalności rozkładów, ani jednorodności wariancji?

Zatem mamy:

\((X_1, \dots, X_m)\), \((Y_1, \dots, Y_n)\) niezależne próby proste o dystrybuantach odpowiednio \(F_x\), \(F_Y\) takich, że \[\exists_{\Delta: \mathbb{R}\rightarrow \mathbb{R}}F_Y(t)=F_X\left(t\right)-\Delta(t),\] gdzie \(\Delta\) jest stałego znaku oraz \(F_X\) jest ciągła.

Zauważmy, że nieznanymi parametrami są \(\theta=(\Delta,F_X)\) - zbiór takich funkcji nie należy do skończenie wymiarowej przestrzeni Euklidesowej - nieparametryczny problem testowania.

Przypomnijmy, że odwzorowanie \[Rang(x_1, \dots, x_m, y_1, \dots, y_n)=(r_1, \dots, r_{m+n})\] jest maksymalnym niezmiennikiem względem grupy przekształceń \[G=\left\{g: \mathbb{R}^{m+n}\rightarrow \mathbb{R}^{m+n}: g\left(\mathbf{x,y}\right)=\left(h(x_1), \dots, h(x_m), h(y_1), \dots, h(y_n)\right), \text{ h- ciągła i ściśle rosnąca}\right\}.\] Zatem, test będzie oparty na wektorze rang \((R_1, \dots, R_{m+n})\) odpowiadającym łącznej próbie \((X_1, \dots, X_m, Y_1, \dots, Y_n)\).

Idea testu Manna-Whitney opiera się na porównaniu sum rang dla \(X\) i \(Y\):

\(\bar{R}_X=\frac{1}{m}\sum_{i=1}^m R_i\) - rangi odpowiadające \(X\),

\(\bar{R}_Y=\frac{1}{n}\sum_{i=m+1}^{m+n} R_i\) - rangi odpowiadające \(Y\).

Zauważmy, że jeżeli:

  • \(\Delta\leq 0\), wtedy \(\mathbf{E}(\bar{R}_Y)\leq \mathbf{E}(\bar{R}_X)\),

  • \(\Delta> 0\), wtedy \(\mathbf{E}(\bar{R}_Y) >\mathbf{E}(\bar{R}_X)\).

Wobec tego,

  • (dla wariantu 3) tj. \(H_0: \Delta\leq 0\), przeciwko \(H_1: \Delta>0\), będziemy odrzucać hipotezę zerową, gdy \(\bar{R}_Y-\bar{R}_X\) będzie duże,

  • (dla wariantu 4) tj. \(H_0: \Delta = 0\), przeciwko \(H_1: \Delta\neq0\), będziemy odrzucać hipotezę zerową, gdy \(\left\|\bar{R}_Y-\bar{R}_X\right\|\) będzie duże,

Wykonując standaryzację, otrzymujemy statystykę testową Manna-Whitney: \[W=\frac{\bar{R}_Y-\bar{R}_X}{\sqrt{V_{\Delta=0}\left(\bar{R}_Y-\bar{R}_X\right)}}=\frac{\bar{R}_Y-\frac{(m+n+1)}{2}}{\sqrt{\frac{m(m+n+1)}{12m}}}.\] Gdy \(F_X=F_Y\) zachodzi \(W\rightarrow N(0,1)\) dla \(n\rightarrow \infty\) oraz \(\min\{m,n\}\rightarrow \infty\).

Ostatecznie, dla $H_0: $, przeciwko alternatywie \(\Delta>0\) mamy \(\phi(W)=1\) dla \(W>c_{\alpha}\), gdzie \(c_{\alpha}=F^{-1}_{N(0,1)}(1-\alpha)\) i jest testem asymptotycznie zgodnym.

Asymptotyczna p-wartość: \(p(W)=1-F_{N(0,1)}(W)\).

Natomiast dla $H_0: = 0 $, przeciwko alternatywie \(\Delta \neq 0\) mamy \(\phi(W)=1\) dla \(|W|>c_{\alpha}\), gdzie \(c_{\alpha}=F^{-1}_{N(0,1)}(1-\frac{\alpha}{2})\) i jest testem asymptotycznie zgodnym.

Asymptotyczna p-wartość: \(p(W)=2\left(1-F_{N(0,1)}(W)\right)\).

Zadanie 3

Pacjenci u których zdiagnozowano pewną chorobę są leczeni jedną z dwóch terapii. Zmierzono czas życia od momentu rozpoczęcia terapii dla \(58\) pacjentów leczonych terapią \(1\) oraz dla \(52\) pacjentów leczonych terapią \(2\). Czy terapie się różnią, a jeśli tak, to która jest lepsza? (Dane są w pliku Baza 6.1)

library(readxl)
Baza<- read_excel("Baza 6.1.xlsx")
Z<-Baza$"Czas"
G<- Baza$"Leczenie"

X <- Z[G == 1]
Y <- Z[G == 2]

m <- length(X)
n <- length(Y)

# Zwrócić uwagę na przesunięcie dystrybuant
par(mfrow=c(2,2))
hist(X,prob=TRUE)
hist(Y,prob=TRUE)

# Wykresy dystrybuant empirycznych
plot(ecdf(X))
lines(ecdf(Y),col="blue")
par(mfrow=c(1,1))

R<-rank(Z)

# Wybór rang dla grupy Y
RY <- R[G == 2]

W<-(mean(RY)-(m+n+1)/2)/sqrt(m*(m+n+1)/(12*n))
p<-2*(1-pnorm(abs(W),0,1))
p
## [1] 0.02400052
p1<-1-pnorm(W,0,1)
p1
## [1] 0.9879997

Zatem, terapie różnią się miedzy sobą rozkładem. Możemy również stwierdzić, że $ $, wobec tego terapia druga jest lepsza (krótsza) od terapii pierwszej.

Możemy również użyć gotowej funkcji wilcox.test.

wilcox.test(X,Y,alternative = "two.sided")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  X and Y
## W = 1885, p-value = 0.02419
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(X,Y,alternative = "less")
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  X and Y
## W = 1885, p-value = 0.9881
## alternative hypothesis: true location shift is less than 0

Obliczanie p-wartości testu Manna-Whitneya metodą MC

Przypomnijmy, że gdy \(\Delta=0\), to wektor rang ma rozkład jednostajny na zbiorze permutacji \((1, \dots, m+n)\). Zatem równoważnie możemy wygenerować zmienne \(\mathbf{U}=U_1, \dots, U_{m+n}\sim(U(0,1))\) i obliczyć \(Rang(U)\).

Kroki:

  1. Generujemy \(mc\) (np. \(mc=1000\)) prób prostych \(U_1, \dots, U_{m+n}\sim U[0,1]\).

  2. Dla każdej z prób z kroku \(1\) obliczamy rangi \(R_1, \dots, R_{m+n}\).

  3. Dla każdego wektora rang z korku \(2\) obliczamy statystykę Manna-Whitney \(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}\left(|W_i|\geq |W|\right).\] Statystyka \(W\) jest obliczona dla oryginalnych danych.

Zadanie 3 - ciąg dalszy

Obliczyć dokładną p wartość metodą MC.

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

for (i in 1:mc){
  UMC <- runif(m + n)
  RMC <- rank(UMC)
  RYMC <- RMC[(m + 1):(m + n)]  
  WMC[i] <- (mean(RYMC) - (m + n + 1) / 2) / sqrt(m * (m + n + 1) / (12 * n))
}
pMC <- mean(abs(WMC) >= abs(W))
pMC
## [1] 0.0229
p1MC <- mean(WMC >= W)
p1MC
## [1] 0.98855

Zadanie

Niech \((X_1, \dots, X_{20})\) - i.i.d. z rozkładu \(N(m_x=0, 1)\) oraz \((Y_1, \dots, Y_{30})\) - i.i.d z rozkładu \(N(m_y, 1)\).

  1. Sporządź wykres mocy empirycznej dla \(m_y\in [-1,2]\) dla testu t-Studenta. (wariant 4: \(H_0: m_X=m_y\), przeciwko \(H_1: m_X\neq m_y\)).

  2. Dla tych samych danych obliczyć p-wartości testu t-Studenta i testu Manna-Whitney. Narysuj wykres różnić p-wartości.

  3. Dodać do wykresu z podpunktu a) moc empiryczną testu Manna-Whitney.