Odporność pewnych testów dwupróbkowych na naruszenie założenia o niezależności prób*

ROCZNIKI POLSKIEGO TOWARZYSTWA MATEMATYCZNEGO Seria III: MATEMATYKA STOSOWANA XXXII (1989)

Ry s z a r d Zi e l i ń s k i

Warszawa

Odporność pewnych testów dwupróbkowych

na naruszenie założenia o niezależności prób*

(Praca wpłynęła do Redakcji 1988.03.20)

1. Sformułowanie zagadnienia. Próba losowa X l , X 2, . . . , X m pochodzi z populacji Fx oraz próba losowa Yl ,Y2, . . Yn pochodzi z populacji FY. Zakładamy, że Fx{x) = F(x — jix) oraz FY{x) = F{x — jiY), gdzie F jest pewną dystrybuantą, a interesują nas testy hipotezy H: /lx = fiY przy jednostronnej

alternatywie K + : jix > jiY. W całej pracy liczby m oraz n traktujemy jako ustalone.

Znane są liczne testy tej hipotezy, np. test t (test Studenta), gdy F jest dystrybuantą rozkładu normalnego, lub test WMW (test Wilcoxona-Man- na-Whitneya), gdy F jest dystrybuantą ciągłą. Wszystkie znane i powszechnie

stosowane testy są konstruowane przy dodatkowym założeniu, że próby

(Xl , X 2, ..., X m) oraz (Y1? Y2, ..., Y„) są niezależne. Trudno jest zresztą

wyob-razić sobie możliwość skonstruowania testu rozważanej hipotezy w przypadku prób zależnych bez precyzyjnego określenia tej zależności. Z drugiej strony, gdyby okazało się, że niektóre ze znanych testów są mało wrażliwe na naruszenie założenia o niezależności prób, konstrukcje specjalnych testów dla prób zależnych stałyby się zbędne.

Interesujące nas zagadnienie własności testów w przypadku prób zależnych pojawiało się w literaturze statystycznej tylko sporadycznie i fragmenta-rycznie. Na przykład Walsh (1947) rozważa tylko przypadek testu t, gdy (m + w)-wymiarowa zmienna losowa (Xl, X 2, ..., X m, Y1,Y2, ..., YJ ma roz-kład normalny o macierzy kowariancji różnej od macierzy diagonalnej. Do studiowania podobnych odstępstw od niezależności zmiennych gaussowskich powraca się później, w nieco innym kontekście, w pracy Schadera i Schmida (1986). Asymptotyczne własności testu WMW przy różnych rodzajach zależ-ności pomiędzy obserwacjami są przedmiotem studiów Serflinga (1968), Hollandem, Pledgera i Lina (1974) oraz Gastwirtha i Rubina (1974).

* Praca wykonana w ramach problemu CPBP 01.02.

W niniejszej pracy rozważamy przedstawione wyżej zagadnienie w duchu koncepcji odporności procedur statystycznych przedstawionej pracy Zieliń-skiego (1983). Ponieważ chcemy przyjrzeć się jednocześnie kilku testom, w tym testowi t, ograniczamy się do sytuacji adekwatnej dla wszystkich tych testów, tzn. do przypadku gaussowskich zmiennych losowych X l , X 2, ..., X m,

Y„Y2,..., Yn.

2. Próby zależne. Walsh (1947) rozważał następującą strukturę kowarian- cyjną zmiennych losowych X 1, X 2, ..., X m, Y1, Y2, ..., Yn:

VarXi — g2,

Var Yj = g'2, Cov(Xj,Xj) = QG2, CoviY^Yj) = q'g'2, Co\{Xt, Yj) = q"g g'.

Odrzucamy to podejście z następującym uzasadnieniem: jeżeli q < 0, to istnieje

tylko skończona liczba zmiennych losowych X l , X 2, X m spełniających

warunek Co\(Xi, X j) = Qa2 dla i# ;' (np. dla q = — j musi być m< 3);

podobne ograniczenie zachodzi dla zmiennych Yl ,Y2,...,Y n i dla całego układu X ltX 2, ..., X m, YX,Y2, ..., Yn. Z punktu widzenia praktycznych za-stosowań takie ograniczenie jest zupełnie nieuzasadnione.

W dalszym ciągu będziemy rozważali próby losowe (X \,X 2, ..., X m) i (yi5y2, ..., y„), których zależność jest opisana w następujący sposbb. Dla danej liczby naturalnej k ^ min{m,n} rozważamy następujący ciąg niezależ-nych elementów losowych

(X, , y j, (X2, Y2),...,(X k,Yk),Xk + l ,..., X m, Ykk + 15

i zakładamy, że (Xt, Yj) ma dwuwymiarowy rozkład normalny

N(h x,iiy,o2,a2,q) dla 1 i ^ k, X { ma rozkład normalny N(h x,g2) dla k-1-1 ^ i ^ m oraz Yt ma rozkład normalny N((j,y,g2) dla k + 1 ^ i ^ n. Liczbę k, podobnie jak liczby m oraz n, uważamy w dalszym ciągu za ustaloną.

Interesujące przykłady tzw. „problemów z życia”, prowadzących do takiego schematu zależności prób, można znaleźć w pracy Hollandem, Pledgera i Lina (1974).

Oto te przykłady.

1. Obserwuje się efekty X t , X 2, ..., X m zabiegu A na m obiektach oraz, w pewnym późniejszym terminie, efekty Y1, Y2, ..., Yn zabiegu B na n obiek-tach. Ze względu na oszczędności lub po prostu niedostatek obiektów, pewne obiekty mogą występować zarówno w pierwszym, jak i w drugim badaniu.

Odporność pewnych testów dwupróbkowych . . . 7 zabieg A), badacz, świadomie lub nieświadomie, dobiera pary, w których jednego szczura poddaje zabiegowi A, a drugiego zabiegowi B. Takie grupowa-nie może wprowadzać pewne zależności między obserwowanymi wynikami.

3. Pewne badania mogą być prowadzone korespondencyjnie i może się zdarzyć, że niektórzy respondenci pojawiają się zarówno w próbce

X x, X 2, ..., X m, jak i w próbce YX,Y2, Yn.

3. Zagadnienie odporności testu. Niech tp będzie testem hipotezy H: px = pY wobec hipotezy alternatywnej K +: px > pY. Rozważamy tylko takie testy, w których

1) rozkład statystyki testowej nie zależy od parametru zakłócającego a2 (przy hipotezie H),

2) moc zależy od px, ptY oraz a2 tylko przez A = (fix — pY)/cr.

Niech (Xx, Yx), (X2, Y2), ..., (Xk, Yk), X k + U ..., X m, Yx, ..., Y„ będą nieza-leżnymi elementami losowymi takimi, że (Xt, Yt), 1 ^ i < k, ma dwuwymiarowy . rozkład normalny N(0,1,0,1, g), X {, k-1-1 ^ i ^ m, oraz Yt, k + 1 ^ i ^ n, mają

rozkład normalny N(0,1).

Wprowadzamy następujące skrótowe oznaczenia:

2Ć + A = (Xx + A,X 2 + A, ..., X m + A), AeR1,

y = ( y 1,y2, y „ ) ,

(X + A,Y) = ((Xl + A,Yl),(X2 + A,Y2),...,(X k + A, Yk), X k + 1 + A, ... . . . , X m + A,Yk +1,..., 7„), A e R 1.

Niech PQ będzie (łącznym) rozkładem (X, 7) zdefiniowanym jako produkt

k gaussowskich rozkładów dwuwymiarowych iV(0,1,0, l,p) oraz m + n — 2k

gaussowskich rozkładów jednowymiarowych N(0,1) i niech EQ będzie operato-rem wartości oczekiwanej odpowiadającym rozkładowi P . Przez pe{x,y) będziemy oznaczali gęstość rozkładu PQ (względem miary Lebesgue’a).

Przy tych oznaczeniach, EQ (p(X_ + A,Y) jest mocą testu ę na alternatywie

A przy współczynniku korelacji w parach równym g oraz Eeę(X^,Y) jest

rozmiarem tego testu. Zakładamy, że test ę jest skonstruowany jako test na z góry ustalonym poziomie istotności a, tzn. że

E0ę{X,Y) = a.

Jeżeli założenie o niezależności prób jest naruszone w ten sposób, że współczynnik korelacji g faktycznie przebiega pewien przedział {gL,gv) zawie-rający zero, to rozmiar cp przyjmuje wartości z przedziału

( inf E„ę{X,Y), sup Eeę(X,Y)).

Odporność testu ę ze względu na jego rozmiar przy opisanym wyżej zaburzeniu założeń o niezależności opisujemy za pomocą liczby

r0(QL,Qu)= sup Eeę{X,Y)~ inf EQę{X,Y).

Oczywiście r^g^Qu) ^ 0, a mniejsze wartości tej wielkości świadczą o większej odporności testu.

Analogicznie będziemy się interesowali tym, jak zmiana współczynnika korelacji q wpływa na moc testu q> przy alternatywach A > 0: odporność testu

(p ze względu na jego moc na alternatywie A przy opisanym wyżej zaburzeniu założeń o niezależności opisujemy za pomocą liczby

^AiQ^Qv)= SUP EQę(X + A,Y)~ inf Ee y(X + A, 7).

Q L ^ e ^ e u e L ^ e^ eu

Zauważmy, że główną rolę w analizie odporności testu ę odgrywa funkcja M „(M ) = Ev (p{X + A,Y)t ql <q ^ qv, A ^ O .

W dalszym ciągu koncentrujemy się na badaniu tej funkcji. Ponieważ badania symulacyjne wykazują, że przy pewnych wartościach A ^ 0 funkcja Mę(-,A) jest monotoniczna (przynajmniej dla rozważanych dalej testów ę), to dla tych wartości A tablice funkcji Mv(q,A) pozwalają na natychmiastową ocenę liczb rA(QL,Qv) dla różnych wartości ql i qv.

4. Analiza wybranych testów. Przedstawimy explicite funkcje Mę(Q,A) dla

wybranych testów dwupróbkowych, a mianowicie dla testu t (Studenta), testu WMW (Wilcoxona-Manna-Whitneya) oraz testu P (testu znaków różnic w parach). Wybór tych właśnie testów uzasadniamy w następujący sposób.

Jednopróbkowe i dwupróbkowe testy t (testy jednostronne i dwustronne na różnych poziomach ufności) okazały się odporne (por. np. Posten (1978)) na wiele różnych odstępstw od założenia o normalności rozkładu; faktycznie wahania rozmiaru tych testów na całej rodzinie rozkładów Pearsona okazały się nieistotne. Z drugiej strony testy Studenta, w tym testy dwupróbkowe, są standardowo używane w wielu praktycznych zastosowaniach technicznych, przyrodniczych i ekonomicznych. W tej sytuacji wydaje się celowe zbadanie odporności tych testów również na różne odstępstwa od założeń o niezależ-ności prób.

Odporność pewnych testów dwupróbkowych ... 9

przy pewnych typach zależności między obserwacjami. Zachowanie się tego testu przy małych próbach nie było jednak, jak nam wiadomo, badane.

Test znaków różnic w parach, z samej swojej konstrukcji powinien być absolutnie odporny na rozważaną przez nas zależność między próbami, przynajmniej jeżeli chodzi o rozmiar tego testu. Pozostaje jednak do wyjaś-nienia, jak odporna jest moc tego testu na takie zaburzenia założeń i jak bardzo moc tego testu różni się w takich sytuacjach od mocy dwóch poprzednich testów. To ostatnie pytanie jest związane z ewentualnym wyko-rzystaniem testu znaków różnic w parach jako przypuszczalnie odpornego konkurenta testu t i testu WMW.

Będziemy używali następujących oznaczeń: f 1, gdy t > 0,

Xy(t) = j y, gdy t = 0,

l 0, gdy t < 0,

gdzie y e [0,1], oraz oznaczenia x(t) dla przypadku y = 0. Odnotujmy, że przy rozważanych założeniach gęstość pe(x,y) wyraża się wzorem

Pa(x,z) = (<J2 ri)~n~m(sJ\ —Q2)~k f i exP | ~ g'2) ~ 2Q W + yfi j

m n

x FI exP { -2X?} II exp { - i yf}.

i = k + 1 i = k + 1

Jeżeli k = 0, to pe{x,y) nie zależy od q.

4.1. Test t (Studenta). Niech t = t(X_,Y) będzie statystyką t Studenta

< ' V V ' - lnm(n + m — 2)'

yj(m — l)S^ + (n— l)Sy

V

gdzie, jak zwykle, X = Zr=i XJm, Y = Z?=1 YJn, S2X = (X,.-X)2/(m - 1) oraz Sy = £"=1 Y)2/(n — 1). Test t ma postać

(p(X,Y) = x(t(X ,Y )-ta,m + n- 2),

gdzie 2 jest kwantylem rzędu 1 -a rozkładu t Studenta o (m + n - 2) stopniach swobody.

Funkcja M ^q^A) tego testu, którą będziemy oznaczali przez Mt(g,A),

przyjmuje postać

Mt(Q,A) = $x{t(x,y)-t0l'm + n- 2)pe(x-A,y)dxdy

zadanym przedziale (QL,QV)• Chociaż formalne różniczkowanie tej funkcji

względem q oraz względem A nie nastręcza tu trudności (można różniczkować

pod znakiem całki), to jednak uzyskane na tej drodze wyrażenia są zbyt skomplikowane na to, żeby można było zbadać zachowanie się tej funkcji na interesującym nas obszarze parametrów. Wartości tej funkcji, uzyskane metodą Monte Carlo na drodze symulacji 1000 prób (X, Y) przy różnych wartościach m, n oraz k, podajemy w załączonej tablicy.

Przede wszystkim odnotujmy następujące

Sp o s t r z e ż e n i e 1. Dla małych wartości A ^ 0, funkcja Mt(-,A) argumentu

{>e(— 1, 1) jest funkcją malejącą; dla dużych wartości A > 0, funkcja M((-, A) jest rosnąca.

Ilustruje to rysunek 1.

Bliższa analiza tablicy i rysunku 1 sugeruje, że rozważany test t nie jest bardzo wrażliwy na małe zmiany współczynnika korelacji q i że ze wzrostem q (dla q > 0) test staje się konserwatywny. Na przykład w małych próbach

Odporność pewnych testów dwupróbkowych ... 11 q do wartości —0.4, oszacowanie rozmiaru testu na poziomie istotności 0.05

wynosi 0.062 („dwa sigma” estymatora dla wszystkich liczb podanych w załą-czonej tablicy nie przekracza ^/O.OOl a 0.032).

„Praktyczny wniosek” z tej analizy sugeruje, żeby — dla zapewnienia rozmiaru testu na poziomie nie przekraczającym założonego poziomu istotnoś-ci a — wybrać wartość krytyczną testu tak, by Mt(gL,0) = ot, gdy ql < 0.

4.2. Test WMW (Wilcoxona-Manna-Whitneya). Rozważamy statystykę

Wilcoxona („suma rang X-ów w połączonej próbie X 1, X 2 . . . , X m,

Y1,Y2,..., Y ”) w następującej postaci:

W(X, Y) = m-n +m(m+ 1) ₂

E E

i - 1 j~ 1

Wtedy test WMW ma postać

<PwMw(X,Y) = Xy{ W ( X ,Y ) - w a),

gdzie wartość krytyczna wa i prawdopodobieństwo randomizujące y są dobrane w taki sposób, że jest to test na poziomie istotności ot.

Podobnie jak w przypadku testu t (§ 4.1), funkcja

Mw mw(Q>A) = $xy(W (x,y)-w 0)pe{x-A,y)dxdy

może być łatwo różniczkowana względem q oraz A, ale wyznaczenie kresu

dolnego i kresu górnego tej funkcji na danym przedziale (ql,q v), przy

ustalonym A, nastręcza duże trudności techniczne. W załączonej tablicy podajemy wartości funkcji Mw mw{q, A) oszacowane metodą Monte Carlo.

Ogólne spostrzeżenie sugeruje, że funkcje Mt(g, A) oraz M WMW(g, A) zachowują się bardzo podobnie jakościowo (monotoniczność) i różnią się nieznacznie swoimi wartościami. Z rezultatów symulacji wynika, że jakościowa ocena wrażliwości testu WMW na rozważane odstępstwa od założenia o niezależności prób jest taka sama jak ocena wrażliwości testu t. Podobnie jak w przypadku testu t, można sformułować

Spo s t r z e ż e n ie 2. Dla małych wartości A ^0, funkcja M wmw(•, A)

ar-gumentu £>e(— 1, 1) jest malejąca; dla dużych wartości A > 0, funkcja

M WMw(‘,A) jest rosnąca.

Podobnie jak w przypadku testu t, na podstawie eksperymentu symulacyj-nego wnioskujemy, że test WMW jest mało wrażliwy na małą korelację i że w przypadku korelacji dodatniej staje się testem konserwatywnymi Rysunek 2 jest odpowiednikiem rysunku 1 dla testu WMW. Podobnie jak w przypadku testu t, odpowiednia korekta wartości krytycznej testu (taka, żeby

Mw mw(ql,0) = ot) może zagwarantować, że rozmiar testu nie przekroczy

4.3. Test P (test znaków różnic w parach). Rozważamy przypadek m = n,

statystykę

i X(X,-Yd

/=1 i test oparty na tej statystyce

(pP{X,Y) = xy( Z x(Xi~

i = i

gdzie wartość krytyczna 7ra i prawdopodobieństwo randomizujące y dobrane są w taki sposób, że jest to test na poziomie istotności a. Funkcja MP(g,A) tego testu ma postać

Mp(q,A) = j z y( Z X(xi- y i) - n <t)pQ{ x - A ti)dxdi.

i = 1

Odporność pewnych testów dwupróbkowych ... 13

i rozkład statystyki 277= i Z (^£ — nie zależy od q. Wynika stąd, że — ze

względu na rozmiar — test P jest testem absolutnie odpornym na rozważane odstępstwa od niezależności prób. Wyniki symulacji (por. załączona tablica) ilustrują ten fakt. Jest natomiast raczej zaskakujące, że moc tego testu reaguje na zmiany współczynnika korelacji. Na przykład dla m = n = k = 4 mamy

MP( — 0.9; 1) =0.180 oraz MP(0.9; 1) = 0.734. Co więcej, moc tego testu

(absolutnie odpornego ze względu na rozmiar) jest wyraźnie bardziej wrażliwa na zmiany współczynnika korelacji niż moc testów t i WMW. Ilustruje to rysunek 3, będący odpowiednikiem rysunku 1 dla testu t i rysunku 2 dla testu

WMW. Wyniki obliczeń metodą Monte Carlo znowu sugerują

Sp o s t r z e ż e n i e 3. Dla A > 0 funkcja MP(-,A) argumentu ge ( — 1,1) jest

rosnąca.

5. Uwagi końcowe. W pracy udzielono odpowiedzi na kilka problemów związanych z odpornością testu t i testu WMW (porównawczo również dla testu P) na pewne odstępstwa od założeń o niezależności prób.

Z praktycznego punktu widzenia zagadnienie można uznać za rozwiązane i zamknięte: praktyk nie tylko otrzymuje informacje o tym, co dzieje się z rozmiarem i mocą testów t i WMW pod wpływem opisanej wyżej korelacji prób, ale otrzymuje również praktyczną wskazówkę jak modyfikować wartości krytyczne tych testów, aby ich rozmiar nie przekraczał założonego poziomu istotności.

Z teoretycznego punktu widzenia praca sugeruje kilka drobnych nowych interesujących problemów (np. dlaczego M^i^A) jest dla pewnych d ^ 0 rosnąca, a dla innych malejąca? czy istnieje A takie, że M(-,d) = const i „co miało by to znaczyć?”), pozostawiając jednak otwartym główny problem teoretyczny: jak skonstruować dobry (= odpowiednio mocny) test odporny na odstępstwa od założeń o niezależności prób. W pracy badaliśmy tylko stosunkowo łatwy przypadek „odstępstw parametrycznych” — zależność między próbami była opisana za pomocą jednowymiarowego parametru q.

Wydaje się, że bardziej naturalne byłoby rozważanie tej problematyki na „zaburzeniach nieparametrycznych”, ale to wymaga dopiero precyzyjnego sformułowania.

Prace cytowane

J. L. G astw irth, H. Rubin (1975), The behavior of robust estimators on dependent data, Ann. Statist. 3, 1070-1100.

M. H ollander, G. Pledger, P. E. Lin (1974), Robustness of the Wilcoxon test to a certain

dependency between samples, Ann. Statist. 2, 177-181.

H. O. P osten (1978), The robustness of the two-sample t test over the Pearson system, J. Stat. Comp. Simul. 6, 295-311.

M. Schader, F. Schmid (1986), How robust is one way Anova with respect to within group

correlation, Compstat 1986. Physica —Verlag, Heidelberg, Wien.

R. J. Serfling (1968), The Wilcoxon two-sample statistic on strongly mixing processes, Ann. Math. Statist. 39, 1202-1209.

J. E. W alsh (1947), Concerning the effect of interclass correlation on certain significance tests, Ann.

Math. Statist. 18, 88-96. *

R. Z ieliń sk i (1983), Robust statistical procedures: a general approach, In: Lecture Notes in

Odporność pewnych testów dwupróbkowych . . . 15 Tablice wartości funkcji Mę(g, A) dla wybranych testów:

<p = t — test t

ę = WMW — test Wilcoxona-Manna-Whitneya

ę — P — test par

Oznaczenia w tablicy: delta = A, RO = q, LT = 1000 x Mt(g,A), LW = 1000 x Mw mw{q, A),

LP = 1000 x Mp(q, A).

Wszystkie wartości funkcji M^g,/.I) oszacowano metodą Monte Carlo na podstawie 1000 niezależnych prób.

Wartości Mv(q,A) dla m = 4, n = 4, II

Odporność pewnych testów dwupróbkowych ... 17

Wartości Mv(q,A) dla m = 5, n = 5, k = 5

delta = 0.0 0.2 0.4 0.6 0.8 1.0 1.5 2.0 3.0 RO = -0.90 LT = 126 165 217 274 351 433 645 829 974 LW = 116 165 206 270 348 420 629 813 973 LP = 61 69 103 143 152 196 336 499 751 RO = -0.60 LT = 93 143 192 259 327 407 671 851 980 LW — 96 140 195 252 313 388 632 838 980 LP = 55 80 99 135 181 209 381 536 810 RO = -0.40 LT = 81 126 176 243 320 394 677 865 988 LW = 77 114 173 231 304 390 632 850 988 LP - 67 87 103 134 193 241 399 596 846 RO = -0.20 LT = 68 103 157 219 297 392 686 879 990 LW = 62 103 147 215 296 383 644 865 987 LP = 52 83 115 160 204 260 447 631 891 RO = -0.10 LT - 56 99 156 217 295 393 694 888 994 LW = 57 93 146 201 282 372 660 862 992 LP = 57 81 110 155 220 287 460 666 924 RO = 0.00 LT = 46 93 140 209 284 395 703 894 995 LW = 48 80 133 203 273 376 665 865 994 LP = 47 80 118 174 233 291 502 699 931 RO = 0.10 LT = 40 81 135 198 273 387 710 905 997 LW = 39 73 126 190 265 368 663 882 993 LP = 54 86 132 168 245 306 511 735 943 RO = 0.20 LT = 35 67 127 195 267 383 717 912 998 LW = 36 63 120 181 262 379» 677 887 995 LP = 42 89 134 208 257 346 554 767 965 RO = 0.40 LT = 23 52 99 169 270 385 727 933 999 LW = 27 50 101 169 257 367 689 909 997 LP = 40 99 165 232 317 425 672 864 985 RO = 0.60 LT = 15 34 78 142 258 395 758 957 1000 LW = 14 41 78 152 241 383 728 932 999 LP = 49 105 172 289 406 518 799 945 999 RO = 0.90 LT = 2 11 38 118 240 403 810 970 1000 LW = 4 17 43 115 250 421 779 949 999 LP = 52 180 408 632 836 947 998 1000 1000