Labolatorium 8

Karolina Marek

2024-08-06

Test Kolmogorowa-Smirnowa

Niech \((X_1, \dots, X_m)\), \((Y_1, \dots, Y_n)\) będzie próbą prostą o rozkładach zadanych przez ciągłe dystrybuanty, odpowiednio \(F_X,F_Y\).

Chcemy przeprowadzić test dominacji, tj:

  1. \(H_0: \forall_{t\in \mathbb{R}} F_Y(t)\leq F_X(t)\), przeciwko \(H_1: \exists_{t\in \mathbb{R}}F_Y(t)>F_X(t)\).

Zauważmy, że przy prawdziwości alternatywy dystrybuanty mogą się przeplatać. Dlatego jest to uogólnienie testu Wilcoxona i Manna- Whitneya.

Na wykładzie pokazane zostało, ż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(\mathbf{x,y})=\left(h(x_1), \dots, h(x_m), h(y_1), \dots, h(y_n))\right), \quad h \text{ - ciągła i ściśle rosnaca}\right\}.\]

Idea testu Kołmogorowa-Smirnowa opiera się na badaniu różnicy pomiędzy dystrybuantami empirycznymi. \[KS_{+}=\sup_{t\in \mathbb{R}}\left(\hat{F}_Y(t)-\hat{F}_X(t)\right),\]

gdzie \(\hat{F}_Y(t)= \frac{1}{n}\sum_{i=1}^n\mathbb{I}(Y_i\leq t)\) oraz \(\hat{F}_X(t)= \frac{1}{m}\sum_{i=1}^m\mathbb{I}(X_i\leq t)\).

Jak pokazano na wykładzie możemy zapisać \(KS_+\) za pomocą unormowanych do przedziału \((0,1)\) rang: \[KS_+(R)=\sup_{t\in (0,1)}\left(\frac{1}{n}\sum_{i=m+1}^{m+n}\mathbb{I}\left(\frac{R_i}{m+n+1}\leq t \right)-\frac{1}{m}\sum_{i=1}^{m}\mathbb{I}\left(\frac{R_i}{m+n+1}\leq t \right)\right).\]

Duże wartości \(KS_+\) będą świadczyć o odrzuceniu hipotezy zerowej.

W przypadku gdy będzie nas interesować test zgodności postaci:

  1. \(H_0: F_Y(t)= F_X(t)\), przeciwko \(H_1: F_Y(t)\neq F_X(t)\),

wystarczy rozważyć statystykę testową postaci: \[KS(R)=\sup_{t\in (0,1)}\left|\frac{1}{n}\sum_{i=m+1}^{m+n}\mathbb{I}\left(\frac{R_i}{m+n+1}\leq t \right)-\frac{1}{m}\sum_{i=1}^{m}\mathbb{I}\left(\frac{R_i}{m+n+1}\leq t \right)\right|.\] I ponownie duże wartości \(KS(R)\) będą świadczyć o odrzuceniu hipotezy zerowej.

Można pokazać, że przy prawdziwości hipotezy zerowej (\(F_x=F_Y\)) rozkład statystyk testowych jest następujący:

  1. \(\sqrt{\frac{mn}{m+n}}KS_+\xrightarrow[\min{m,n}\rightarrow \infty]{d}sup_{t\in(0,1)}B(t)\),

  2. \(\sqrt{\frac{mn}{m+n}}KS\xrightarrow[\min{m,n}\rightarrow \infty]{d}sup_{t\in(0,1)}|B(t)|\),

gdzie \(B(t)\) jest mostem Browna.

Szybkie obliczanie statystk \(KS_+\), \(KS\)

Będziemy do tego używać wzorów z wykładu:

  1. \(KS_+=\max_{k\in \{1, \dots, m+n\}}S_k\),

  2. \(KS=\max_{k\in \{1, \dots, m+n\}}|S_k|\),

gdzie \(S_k\) jest zadane w sposób rekurencyjny: \(S_1=T_{[1]}\), \(S_{k+1}=S_k+T_{[k+1]}\) oraz

Obliczanie p-wartości testu Kołmogorowa-Smirnowa metodą MC

Przypomnijmy, że przy prawdziwości hipotezy zerowej (\(F_x=F_Y\)) rozkład rang \((R_1, \dots, R_{m+n})\) jest jednostajnie rozłożony na zbiorze permutacji \((1, \dots, m+n)\).

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

  2. Dla każdej próby z kroku \(1\) wyznaczamy statystyki \(KS\) lub \(KS_+\), otrzymując probe prostą \(V_1, \dots, V_{mc}\).

  3. Dla \(H_0: \forall_{t\in \mathbb{R}} F_Y(t)\leq F_X(t)\), przeciwko \(H_1:\exists_{t\in \mathbb{R}} F_Y(t)>F_X(t)\) mamy \[\hat{p}=\frac{1}{mc}\sum_{i=1}^{mc} \mathbb{I}(V_i\geq KS_+).\] Dla \(H_0: F_X = F_Y\), przeciwko \(H_1: F_X\neq F_Y\) mamy \[\hat{p}=\frac{1}{mc}\sum_{i=1}^{mc} \mathbb{I}(V_i\geq KS).\]

Zadanie 1

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\). Przetestować czy terapia \(1\) jest lepsza bez zakładania typu rozkładów oraz przesunięcia dystrybuant. (Dane są w pliku Baza 6.1).

Przeprowadzimy dwa testy:

a)\[H_0: \forall_{t\in \mathbb{R}} F_Y(t)\leq F_X(t)\qquad H_1:\exists_{t\in \mathbb{R}} F_Y(t)>F_X(t),\]

(postępowanie analogiczne, jak wyżej)

  1. \[H_0: \forall_{t\in \mathbb{R}} F_X(t)\leq F_Y(t)\qquad H_1:\exists_{t\in \mathbb{R}} F_X(t)>F_Y(t).\]

Zauważmy, że tutaj statystyka testowa będzie postaci:

\[KS'_+(R)=-KS_+(R)=\sup_{t\in (0,1)}\left(-\frac{1}{n}\sum_{i=m+1}^{m+n}\mathbb{I}\left(\frac{R_i}{m+n+1}\leq t \right)+\frac{1}{m}\sum_{i=1}^{m}\mathbb{I}\left(\frac{R_i}{m+n+1}\leq t \right)\right).\] Zatem \(S'_{k+1}=-S_{k+1}=-S_k-T_{[k+1]}\) i \(KS'_+=\max_{k\in \{1, \dots, m+n\}}S'_k\).

Najpierw policzmy statystykę \(KS_+\) oraz \(KS'_+\).

library(readxl)
Baza<- read_excel("Baza 6.1.xlsx")
Z<-Baza$"Czas"
G<- Baza$"Leczenie"
# Podział danych na grupy X (G=1) i Y (G=2)
X <- Z[G == 1]
Y <- Z[G == 2]

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

plot(ecdf(X), col="red")
lines(ecdf(Y),col="blue")
legend("bottomright", legend = c("Grupa 1 - X", "Grupa 2 - Y"), col = c("red", "blue"), lty = 1)

T <- ifelse(G == 1, -1/m, 1/n)
V <- cbind(Z, T)
V <- V[order(V[,1]),]

SP <- cumsum(V[, 2])
SN <- -cumsum(V[, 2])

KSP <- max(SP)
KSN <- max(SN)
KSP
## [1] 0.2632626

Następnie policzmy p-wartość.

set.seed(321)
MC <- 10000
KSPMC <- numeric(MC)
KSNMC <- numeric(MC)
N <- length(Baza$"Czas")
for(j in 1:MC){
  ZMC <- runif(N, 0, 1)
  
  TMC <- c(rep(-1/m, m), rep(1/n, n))
  VMC <- cbind(ZMC, TMC)
  VMC <- VMC[order(VMC[,1]), ]
  
  SPMC <- cumsum(VMC[, 2])
  KSPMC[j] <- max(SPMC)
  SNMC <- -SPMC
  KSNMC[j] <- max(SNMC)
}
pKSP <- mean(KSPMC > KSP)
pKSP
## [1] 0.0176
pKSN <- mean(KSNMC > KSN)
pKSN
## [1] 0.9032

Odrzucamy hipotezę zerową dla przypadku a), ale przyjmujemy hipotezę zerową z przypadku b). Zatem, pacjenci leczeni terapią \(1\) mają dłuższy czas życia, czyli jest ona lepsza.

Możemy również użyć gotowej funkcji ks.test():

ks.test(Y, X, alternative="greater")
## 
##  Exact two-sample Kolmogorov-Smirnov test
## 
## data:  Y and X
## D^+ = 0.26326, p-value = 0.01759
## alternative hypothesis: the CDF of x lies above that of y
ks.test(Y, X, alternative="less")
## 
##  Exact two-sample Kolmogorov-Smirnov test
## 
## data:  Y and X
## D^- = 0.034483, p-value = 0.911
## alternative hypothesis: the CDF of x lies below that of y

Test Andersona - Darlinga

Test ten opiera się na wykorzystaniu zbieżności statystyk do mostu Browna, przy prawdziwości hipotezy zerowej (\(F_X=F_Y\)).

  1. \(\sqrt{\frac{mn}{m+n}}KS_+\xrightarrow[\min{m,n}\rightarrow \infty]{d}sup_{t\in(0,1)}B(t)\),

  2. \(\sqrt{\frac{mn}{m+n}}KS\xrightarrow[\min{m,n}\rightarrow \infty]{d}sup_{t\in(0,1)}|B(t)|\).

Zatem, z definicji mostu Browna mamy:

\[\forall_{t\in (0,1)}\sqrt{\frac{mn}{m+n}}KS_+\xrightarrow[\min{m,n}\rightarrow \infty]{d}N(0,t(1-t)).\] Standaryzując mamy:

\[\forall_{t\in (0,1)}\sqrt{\frac{mn}{m+n}}\frac{KS_+}{\sqrt{t(1-t)}}\xrightarrow[\min{m,n}\rightarrow \infty]{d}N(0,1).\] Miarą odchyłki w teście Andersona- Darlinga będzie norma \(L_2\):

  1. dla testu stochastycznego uporządkowania mamy

\(AD_+=\int_0^1\frac{\left[\left(\frac{1}{n}\sum_{i=m+1}^{m+n}\mathbb{I}\left(\frac{R_i}{m+n+1}\leq t\right)-\frac{1}{m}\sum_{i=1}^{m}\mathbb{I}\left(\frac{R_i}{m+n+1}\leq t\right)\right)_+\right]^2}{t(1-t)}dt\), gdzie \(h(t)_+=h(t)\mathbb{I}(h(t)>0)\).

  1. dla testu zgodności mamy

\(AD=\int_0^1\frac{\left[\frac{1}{n}\sum_{i=m+1}^{m+n}\mathbb{I}\left(\frac{R_i}{m+n+1}\leq t\right)-\frac{1}{m}\sum_{i=1}^{m}\mathbb{I}\left(\frac{R_i}{m+n+1}\leq t\right)\right]^2}{t(1-t)}dt\).

Statystyki te można zapisać jako:

  1. \(AD_+=\sum_{k=1}^{m+n+1}\left[(S_k)_+\right]^2\log\left(\frac{(k-1)(m+n-k+1)}{k(m+n-k)}\right)\),

  2. \(AD=\sum_{k=1}^{m+n+1}(S_k)^2\log\left(\frac{(k-1)(m+n-k+1)}{k(m+n-k)}\right)\),

gdzie \(S_k\) jest obliczane tak jak w teście Kołmogorowa-Smirnowa.

P- wartości obliczane są metodą MC analogicznie, jak w przypadku testu Kołmogorowa-Smirnowa.

Zadanie 1 - cd.

Wykonaj teraz test Andersona-Darling i porównaj do wcześniej otrzymanych wyników. Podobnie, tutaj testujemy:

a)\[H_0: \forall_{t\in \mathbb{R}} F_Y(t)\leq F_X(t)\qquad H_1:\exists_{t\in \mathbb{R}} F_Y(t)>F_X(t),\]

  1. \[H_0: \forall_{t\in \mathbb{R}} F_X(t)\leq F_Y(t)\qquad H_1:\exists_{t\in \mathbb{R}} F_X(t)>F_Y(t).\]
set.seed(321)
ADP <- sum(ifelse(SP > 0, (SP^2) * log(((1:N) + 1) * (m + n - (1:N) + 1) / ((1:N) * (m + n - (1:N)))), 0))
ADN <- sum(ifelse(SN > 0, (SN^2) * log(((1:N) + 1) * (m + n - (1:N) + 1) / ((1:N) * (m + n - (1:N)))), 0))

ADPMC <- numeric(MC)
ADNMC <- numeric(MC)
for(j in 1:MC){
  ZMC <- runif(N, 0, 1)
  TMC <- c(rep(-1/m, m), rep(1/n, n))
  VMC <- cbind(ZMC, TMC)
  VMC <- VMC[order(VMC[,1]), ]
  SPMC <- cumsum(VMC[, 2])
  SNMC <- -SPMC
  indices <- 1:(N-1)
  log_terms <- log((indices + 1) * (m + n - indices + 1) / (indices * (m + n - indices)))
  ADPMC[j] <- sum((SPMC[indices] > 0) * (SPMC[indices]^2) * log_terms)
  ADNMC[j] <- sum((SNMC[indices] > 0) * (SNMC[indices]^2) * log_terms)
}
pADP <- mean(ADPMC > ADP)
pADN <- mean(ADNMC > ADN)
pADP
## [1] 0
pADN
## [1] 0.8578

Zadanie

Wyznaczyć \(6\) trajektorii mostów Browna, wykorzystując fakt, że statystyka \(\sqrt{\frac{mn}{m+n}}\left(\frac{1}{n}\sum_{i=m+1}^{m+n}\mathbb{I}\left(\frac{R_i}{m+n+1}\leq t \right)-\frac{1}{m}\sum_{i=1}^{m}\mathbb{I}\left(\frac{R_i}{m+n+1}\leq t \right)\right)\rightarrow B(t)\). W tym celu, należy wygenerować próbę prostą \(X_1, \dots, X_m, Y_1, \dots, Y_n\sim U[0,1]\), gdzie \(n=1000\), \(m=10000\), obliczyć statystyki \(S_k\), a następnie wygenerować wykresy dla których na osi \(y\) będą się znajdować statystyki \(\sqrt{\frac{mn}{m+n}} S_k\), a na osi \(x\) rangi wektora \((X,Y)\) zawężone do przedziału \((0,1)\), tzn powinniśmy mieć wektor \(\left(\frac{1}{m+n+1},\frac{2}{m+n+1}, \dots,\frac{m+n}{m+n+1}\right)\). Całą procedurę należy powtórzyć dla \(6\) różnych losowych prób.