Badanie założenia o rozkładzie normalnym błędów losowych w dwuczynnikowym układzie split-plot*

H^{en r y k} Br z e s k w in ie w ic z, W^{ie sła w} W^{a g n e r}

Poznań

Badanie założenia o rozkładzie normalnym błędów losowych w dwuczynnikowym układzie split-plot*

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

1. Wstęp. Jednym z najczęściej przyjmowanych założeń stochastycznych przy opracowywaniu wyników doświadczalnych pochodzących z dwuczynni- kowego układu split-plot jest rozkład normalny błędów doświadczalnych odpowiadających dużym i małym poletkom. Założenie to powinno być weryfikowane dla aktualnie analizowanych danych doświadczalnych. Weryfi- kacji takiej dokonujemy testami normalności dla próby prostej. Odpowiednią próbę prostą uzyskujemy z połączenia przekształconych obserwowalnych zmiennych losowych oraz wygenerowanych zmiennych losowych zgodnie z zasadą randomizacji Durbina. Normalność sprawdzamy testem D’Agostino, odpowiednim dla dużych prób i charakteryzującym się dobrymi własnościa- mi mocy.

W pracy podajemy metodę wyznaczania próby prostej dla badania normalności błędów doświadczalnych małych i dużych poletek w układzie split-plot, którą również ilustrujemy przykładem liczbowym.

2. Model analizy wariancji. Niech A i B będą dwoma badanymi czynnika- mi doświadczalnymi występującymi odpowiednio na a i b poziomach. Pozio- my czynnika A rozmieszczamy losowo w blokach zgodnie ze schematem bloków zrandomizowanych kompletnych. Jednostki doświadczalne, na któ- rych umieszczono poziomy czynnika A, zwane dużymi poletkami, dzielimy na równe części, zwane małymi poletkami. Małym poletkom przyporządkowuje- my losowo poziomy czynnika B. Poziomy obu czynników traktujemy jako ustalone, tzn. zajmujemy się jedynie modelem stałym analizy wariancji.

Zapiszmy model matematyczny obserwowanej zmiennej losowej yUj dotyczą-

PAN 501-73/254/86

czej /-tego bloku, i-tego poziomu czynnika A oraz ./-tego poziomu czynnika B w postaci

(1)' ^{y u j} = H + ^?i + Qi + dii + ^{5 j} + ^{W i j} + ^{e iij,}

dla i = 1, a; / = 1, r; j = 1, b, gdzie n jest parametrem ogólnym, i,- — efektem i-tego poziomu czynnika A, — efektem /-tego bloku, <51 — efektem j-tego poziomu czynnika B, (ió)ij — efektem interakcji i-tego poziomu czynnika A z 7-tym poziomem czynnika B. Symbole dit oraz eitj oznaczają nieobserwowalne zmienne losowe związane z błędami doświadczalnymi od- powiednio dużych i małych poletek. O zmiennych losowych du oraz zakładamy, że są niezależne o wartościach oczekiwanych równych zero oraz wariancjach a\ i o\. Przyjęte założenia o zmiennych da oraz eaj implikują następujące kowariancje:

(2)

Cov (dih dvl') = 0,

i = i', l = /', poza tym, Covieuj, evvy) =

0,

i = i', / = /', j = /, poza tym,

Cov{dih evvy) = 0 dla wszystkich i, l,j, i', /',/.

Z (1) i (2) wynikają dla yiU kowariancje U l +*2, (3) Co\(ynj, yn -f ) = <aj,

(0,

Na parametry stałe modelu (1) nakładamy walne restrykcje

i = i', ł = r,j = /, i = i , / = /»j j ? poza tym.

bez straty ogólności nieestymo-

(4) _{i= 1}£ T«-= Z _{/= 1}Qi = Z ^ = Z (T<% = Z (T<% = °-_{j= 1 i= 1 /=1}

Przy restrykcjach (4) najlepszymi nieobciążonymi estymatorami liniowymi parametrów stałych modelu (1) z uogólnionej metody najmniejszych kwadra- tów (uogólnionej MNK) są (patrz np. Mikos i Niedokos, 1976):

A = ^{y . . . ,}

*i = * = h ..., Q,

(5) Qi = y.i-y.. 1 = 1, ..., r,

Sj = y.j-y.. j = U ..., b,

W u = Jij-y.-. 1 + i = 1, • ••, a; j = 1, ..., b,

gdzie , y,.., ... oznaczają odpowiednie średnie, przy czym uśrednienie

przebiega względem wskaźnika (lub wskaźników) zastąpionych kropką (lub kropkami).

W analizie wariancji rozpatrywanego układu doświadczalnego dla wyzna- czenia sum kwadratów dla błędów dużych i małych poletek wykorzystuje się reszty z uogólnionej MNK dla:

(a) dużych poletek

du = y n - f i - h - Q i = yn.-yi.-y.i.+y..., przy i = 1, a; l = 1, r;

(b) małych poletek

tu j = y u j - fi - i i - Q i - Sj - ( i 3)ij - du =

= y u j- y n - y u + y i.., przy i = 1, a; 1=1, ..., r; j = 1, b.

Podane reszty wykorzystujemy dalej do badania normalności. Nasze postępowanie opieramy na następującym lemacie Cramera, wyrażającym własność charakteryzacji rozkładu normalnego (patrz np. Rao 1982, s. 525).

Le m a t. Jeżeli U x i U 2 są dwiema niezależnymi zmiennymi losowymi takimi, że U = U X + U2 ma rozkład normalny, to zarówno U x, jak i U2 mają rozkład normalny.

Przyjmując, że błędy dn oraz eUj odpowiadają zmiennym Ux i U2 w podanym lemacie, wyznaczamy nowe reszty jako sumę reszt dit oraz eiljt które nazwiemy resztami połączonymi. Te ostatnie będą podstawą do badania normalności zmiennych d^ + euj, co na mocy podanego lematu będzie świad- czyć o normalności każdej ze zmiennych oddzielnie. Jest to przyjmowane zwykle założenie stochastyczne o błędach losowych w układach split-plot.

3. Wektor reszt połączonych. Wykorzystując dit oraz eiU podane w punk- cie 2, wyznaczamy reszty połączone, które oznaczymy przez

(6) fuj = i , + eitj = yaj - y Lj - y.,. f y....

Reszty f aj są zmiennymi losowymi o wartościach oczekiwanych równych zero dla wszystkich wskaźników i, l, j, spełniają z prawdopodobieństwem 1

a r b

warunek X Z Z fuj = a ich wariancje i kowariancje wyznaczamy z

i = i /= 1 ^{j =} 1 równości

(7) Cov (jw. Jrvf) = Cov {yiU, yiTf) - 2 Cov {ynj, yr j ) -

— 2 Co\(yuj, y.r.) + 2 Cov y...) + Cov(.Vij, yr ./ ) + + 2Cov(>’;.,, ^{y .,')-}2Cov^(yLjJ y j + Cov(>’.,., ^yr.r)~

- 2 Cov y...) + Cov( v >...).

4 — Matematyka Stosowana t. 31

Przy ustalaniu kowariancji (tabela 1) występujących po prawej stronie w (7) korzystamy z oznaczeń (2) i (3), a ponadto wprowadzamy y =

Co = ba2 + 02, v = (To/n, Ł, = a%jab i n = arb.

Tabela 1

Cov(v ) i = i’

j= f

i = i' 1 = 1' i

i = i' ///' j= f

i = i' /#/' i * f

i * i' 1 = 1' i =f

i * i' 1 = 1' j * f

i * i' i* r j = f

i * i' 1*1' )*)'

yrrr y 0 0 0 0 0 0

k u, k j y/r y/r a\/r a\/r 0 0 0 0

Kij, y.r. 0 z, 0 0 0 0 0

Ku, y... ^{V V V V V V V V}

K.j, K’.j■ y/r y/r a\/r 0 0 0 0

yi.j, y.v. ^{V V V V V V V V}

k.j, y... ^{V V V V V V V V}

y.i., K’.j- £ 0 £ 0 Ć 0 Ć 0

y.i., y... ^{V V V V V V V V}

y..., y... ^{V V V V V V V V}

Wielkości podane w tabeli 1 wyznaczono, korzystając z własności kowarian- cji, co ilustruje poniższy przykład:

Cov(j/ , , y j = Cov ( ^ £ £ t t ,1 ) - 1

fl2rb2 1 a2rb2

1 a2rb2

1 a2rb2

f £ Cov(yiu> £ £ £ - i=lj=l i' = 1 I' = 1 / = 1

Z Z Cov(yiy, X Z JVy') =

i = 1 j=-l i'=l/=l

[Z z

f Z Covly^,^)] =

\ab((rl + o\) + ab{b- \)o2) = - (o’?-I-<r|-ł-feo-? — o"f) =n

= {bo\ + ol)ln = al/n = v.

Otrzymany wynik w 9 wierszu będzie identyczny dla każdej kolumny w tabeli 1, gdyż występuje tutaj tylko jeden wskaźnik /. Wyniki zawarte w tabeli 1 pozwalają wyznaczyć kowariancje (7), przy różnych układach wskaźników i, /, j, i’, /', / :

' (1-1 /r)(y -0 ,

— y/r + v,

(8) Co v {fij, Iwy) -1 v,

Reszty f Uj nie spełniają warunków próby prostej, tzn. ciągu nieskorelowanych zmiennych losowych o identycznych rozkładach, gdyż istnieją różne od zera kowariancje między tymi resztami. Uniemożliwia to zastosowanie znanych testów normalności bezpośrednio do tych reszt. Reszty f Uj przekształcamy w myśl metody randomizacji Durbina (Durbin, 1961) do nowych reszt (co pokazujemy w punkcie czwartym), które nazywamy poprawionymi resztami połączonymi. Te ostatnie będą spełniały warunki próby prostej.

4. Poprawione reszty połączone. Wyznaczanie poprawionych reszt połą- czonych przeprowadzamy w dwóch etapach. W pierwszym zapisujemy kowa- riancje podane wzorem (8) w postaci macierzowej, wykorzystując w tym celu macierze ortonormalne. W drugim budujemy pewną macierz będącą ich kombinacją liniową ze współczynnikami zależnymi jedynie od a, r, b, i cr|.

Tak uzyskana macierz jest przemnożona przez wygenerowany wektor losowy o składowych będących realizacjami zmiennych losowych o rozkładzie N(0,1). Utworzony wektor jest składnikiem korygującym wektor reszt połą- czonych, co w efekcie prowadzi do poprawionego wektora reszt połączonych.

Zapiszmy reszty f nj leksykograficznie w postaci n-wymiarowego wektora f = (/m ,/ii2 , ...,farb)'• Dla zapisania macierzy wariancji-kowariancji wekto- ra f wprowadzamy macierze ortonormalne (Brzeskwiniewicz, 1980)

gdzie (g) oznacza iloczyn kroneckerowski, Xjf = 1 T/s, X(!S) = 1 — 1 V/s są ma- cierzami s-tego stopnia (s = a, r, b), I oznacza macierz jednostkową, nato- miast IT jest macierzą jedynek. Macierz wariancji-kowariancji wektora f wyraża się wówczas w postaci

(9) X,kcl c 2 c 3

(10) D(f) = X000 + g2 X001 + 03 Xolo + 04 X0u +9s X100 + + 06 Xioi + 07 Xno + 08 Xui,

gdzie współczynniki gt, i = 1, . . . , 8, wyznaczamy z równości g = Zł.

W równości tej g = (0X, ..., 08)', ł = ((1 — 1/r)(y- ol/ab), ..., ol/n) jest wekto- rem złożonym z elementów kowariancyjnych występujących w (8), natomiast

elementy ztt>, t, t' = 1, 8, macierzy Z wyznaczamy wzorem

(U )

ztt' = [(1 - a j ) ( a - lJ-flrT 1 C(1 — ci2)(r — 1) — a2]C2 [(1 - a 3)(b - l ) - a 3]c3, gdzie c!, c2, c3 = 0, 1 oraz t = 4at + 2a2 + a3 + 1, t' = 4c! + 2c2 + c3 + 1. Ma- cierz Z można nazwać macierzą przejścia między macierzami partnerów (Yamoto, Fujii, Hamada, 1956, Brzesk winiewicz 1980) a macierzami ortonor- malnymi (9). Dla macierzy Z suma elementów pierwszego wiersza jest równa n = arb, natomiast pozostałe wiersze sumują się do zera, co łatwo można sprawdzić, zapisując elementy (11) w postaci jawnej:

b - \ r — 1 ( r - l ) ( b - \ ) ,a - \ (a -l)(fe -l) ( o - l ) ( r - 1) (a—l)(r—l)(6 —

-1 ^{r — 1} 1—r a - 1 1 - a (fl-D (r -l) — («— l)(r - I) b - 1 -1 1^{- b} a - 1 ( a - l ) ( b - l ) 1 - a - ( a - l ) ( b - l )

-1 -1 1 ^{a - 1 1 - a} 1 ^{—a a — 1}

b - 1 r — 1 ( r - l ) ( b - \ ) -1 1 - b \ - r - ( r - l ) ( b - l )

-1 r — 1 1 ^{—r -1} 1 1 —r r - 1

b - 1 -1 1 - b -1 1 - b 1 b - 1

-1 -1 1 -1 1 1 -1

Postać macierzy wariancji-kowariancji (10) wektora f ^skazuje na sposób wyznaczania wektora poprawiającego. Niech h = max (grj, d( = yfh — gi, i

= 1,..., 8, oraz

(12) X — di X00o + ^X001 +d3 X010-M4X011 4-^5 X100 + + d6X 101+d7X 110 + d8X l l l ,

gdzie macierze ortonormalne podano w (9). Oznaczmy jeszcze przez w n- wymiarowy wektor o rozkładzie N(O, I) i niezależny od wektora f. Tworzy- my wektor poprawionych reszt połączonych,

(13) v = f+X w.

Wektor losowy v ma wektor wartości oczekiwanych E(\) = £(f) + X£(w) = 0 oraz macierz wariancji-kowarancji

D(\) = D(f) + XD (w) X' = Z) (f) + X X' =

= 9\ X o o o + ••• + 0 8 X i n +(h — gi ) X 000+ . . . +(h — gs)X m =

= MX0oo+ + X m ) = hl.

Pokazaliśmy więc, że składowe wektora v, przy znanych o\ i c\ spełniają warunki próby prostej. Badanie normalności tej próby przeprowadzamy testem normalności D’Agostino podanym w punkcie 5. W przypadku niezna-

nych parametrów o\ i a\, zastępujemy je ich nieobciążonymi estymatorami z uogólnionej MNK, odpowiednio:

8\ = Ż Z Z e y a ( r - l ) ( b - l) , i= 11=1 7=1

S l= {8 l-S j)/b ,

°l = b t t d f i / ( a - l)(r-l),

* = i i=i

gdzie reszty dit oraz eiy podano w punkcie 2.

5. Wnioskowanie statystyczne o rozkładzie normalnym błędów losowych.

Proponujemy dla sprawdzenia normalności próby losowej złożonej ze skła- dowych wektora v test DA D’Agostino (1971). Test ten charakteryzuje się właściwościami testu omnibus, tzn. jest wrażliwy na asymetrię i kurtozę rozkładu (patrz np. Wagner, 1982). Oznacza to, że hipoteza zerowa o normalności jest odrzucona, gdy rozkład z próby wykazuje asymetrię i kurtozę. Test DA ma wysoką moc dla szerokiej klasy alternatywnych (róż- nych od rozkładu normalnego) rozkładów, jest niezmienniczy ze względu na przekształcenia liniowe obserwacji próbkowych, a ponadto dostępne są dla niego tablice wartości krytycznych (podają je m.in. Wagner i Błażczak, 1986).

Niech ciąg vt , ..., vn stanowi próbę n-elementową składowych wektora v, v — średnią arytmetyczną z tej próby, a s2 = £ (u, —u)2 — sumę kwadratówn

odchyleń od średniej v. Próbę v1, . . . , v H uporządkowaną według wielkości i= 1 niemalejących wyraża ciąg nierówności u(1) ^ i?(2) ^ ... < v(n), który traktuje- my jako próbkową realizację statystyk pozycyjnych Kl) < V(2) < ... <

Statystyka testowa testu D’Agostino ma postać

(141 Y= V"(P^ - ° ’282095)

' ’ 0,029986

gdzie

Test D’Agostino jest dwustronny. Jeżeli zachodzi nierówność Y*,2 , ^

^ Ti-a/2;n»to na poziomie istotności a nie mamy podstaw do odrzucenia hipotezy o normalności rozkładu próby vt , ..., vn.

Wartości krytyczne Ya:n odczytuje się z odpowiednich tablic. Brak pod- staw do odrzucenia wspomnianej wyżej hipotezy oznacza, że nie mamy także podstaw do odrzucenia hipotezy o rozkładzie normalnym wektora f, co wynika z (13). Tym samym w myśl podanego lematu brak podstaw do odrzucenia hipotezy o rozkładzie normalnym dla dużych i małych poletek.

6. Przykład liczbowy. Dane doświadczalne pochodzące z COBORU Słu- pia Wielka z doświadczenia dwuczynnikowego split-plot przedstawiają plony ziarna zebrane z poletek o powierzchni 25 m2. Duże poletka (czynnik A) stanowiły 4 poziomy nawożenia mineralnego, natomiast małe poletka (czyn- nik B) stanowiło 6 odmihn jęczmienia ozimego. Doświadczenie było przepro- wadzone w SDO Rarwino. Plony z 3 powtórzeń doświadczenia założonego w układzie split-plot były następujące:

Powtórzenie 1 Powtórzenie 2 Powtórzenie 3

A\ a2 A3 A 4 A! A 2 «3 Aą A1 ^a2 A 3 A^ą.

Bi 6,8 3,8 3,7 4,4 8,5 11,5 3,7 6,5 8,5 ^10,8 8,1 8,6

B2 5,7 4,5 2,8 2,8 5,4 8,5 3,2 8,1 7,5 7,7 5,7 5,7

B3 5,7 5,5 5,5 3,4 5,6 10,6 4,0 5,5 8,1 10,4 5,3 6,2

S4 6,5 5,4 5,8 1,9 6,9 9,4 4,8 4,8 6,4 8,6 5,8 7,9

B5 6,1 3,9 4,6 3,4 6,5 11,0 3,8 5,7 6,0 8,9 6,9 7,9

B6 3,2 3,2 2,0 1,6 5,4 8,0 2,5 3,6 5,9 7,2 3,5 3,7

Powyżej przez A u .. . ., A^ą oznaczono kolejne dawki nawożenia (poziomy czynnika A), natomiast Bi , B6 oznaczają użyte w doświadczeniu odmia- ny jęczmienne (poziomy czynnika B). Zamierzamy przeanalizować wyniki doświadczalne od strony ewentualnego naruszenia założenia o rozkładzie normalnym błędów losowych dla dużych i małych poletek. Dla podanych obserwacji doświadczalnych mamy: a — 4, r = 3, b = 6 i n = 72. Oceny nieobciążone wariancji wynoszą: 8\ = 1,60324, 8\ = 0,98311 i 8% = 10,60255.

Wyznaczamy stałe gt, i = 1, ..., 8, jako składowe wektora g = Zł = (7.46035, 8.44347, -0.98311, 0.98311, -8.44347, 3.14219, -0.98311, 0.98311)', a na- stępnie reszty połączone f Uj według (6):

0.55278 1.18611 0.91944 1.58611 1.58611 0.05278 0.10694 -1.25972 -1.32639 -0.15972 -0.15972 0.10694 -0.65972 0.07361 0.40694 -1.42639 -1.42639 -0.15972 -3.21389 -0.71389 -1.64722 -0.71389 -2.34722 -1.24722 2.34028 1.14028 1.30694 1.14028 2.60694 1.40694 0.87361 -0.42639 0.34028 -0.42639 -0.25972 -0.15972 0.21944 0.25278 2.25278 2.08611 1.18611 1.01944 -1.92639 -1.49306 -1.39306 -1.05972 -1.75972 -0.62639 1.70694 1.24028 -0.85972 -1.02639 0.57361 -0.39306 -0.41389 -1.04722 0.05278 -1.28056 -0.58056 0.25278 -0.45972 2.10694 0.00694 -0.52639 -0.42639 0.30694 0.87361 -1.05972 -0.05972 1.80694 1.00694 -0.55972

Pomijamy szczegółowe obliczenia związane z wektorem poprawiającym Xw w (13), gdyż wymagałoby to podawania dużej ilości wyników pośrednich obliczeń. Liczby losowe normalne wygenerowane zostały według metody podanej m.in. przez Zielińskiego (1972, s. 83), przy wykorzystaniu liczb losowych z rozkładu jednostajnego, wyznaczonych algorytmem cytowanym w

monografii Hellwiga (1975, s. 405). Poprawione reszty połączone obliczone

w e d ł u g (13) s t a n o w i ą c i ą g :

3.63808 - 4 .0 5 0 7 2 4 .6 4 4 3 9 - 4 .1 2 5 9 1 - 3 .1 1 2 1 3 - 4 .1 4 5 0 6 - 5 .0 1 0 3 1 - 5 .6 6 1 6 7 - 5 .8 2 7 4 5 - 2 .6 5 7 2 8 0.79301 1.75282 - 4 .8 7 7 7 1 1.32968 1.75681 - 1 .1 1 0 5 5 0 .1 6454 - 3 .0 1 8 4 9 - 2 .7 6 4 0 8 - 0 .1 4 3 5 6 - 4 .0 2 3 4 3 3.22964 - 2 .0 7 7 1 9 - 3 .0 8 7 1 9

1.07492 0.25307 3.1 3 9 1 4 - 0 .6 2 3 7 5 1.98073 1.96475

1.23710 4.83451 0 .0 6 4 0 6 - 0 .4 9 2 1 6 - 1 .7 9 5 6 5 1.61398 0.52894 - 2.46962 4 .2 7 2 2 6 - 1 .9 4 3 8 9 - 3 .3 7 2 3 6 - 1 .5 2 7 0 1 - 5 .5 5 8 7 3 - 3.87278 - 3 .2 3 8 0 1 - 0 .6 3 8 0 6 - 3 .5 7 5 7 5 1.74871

3.64608 1.09093 - 2 .7 9 9 9 8 2.18967 1.37057 - 1 .8 6 7 1 2

- 0 .3 1 4 1 3 0 .1 7454 - 2 .4 6 4 2 7 1.99926 - 2 .5 9 7 1 2 5.44857

4.49225 5.14699 2.6 5 6 2 6 1.17734 2.39643 - 0 .4 6 7 0 1

2.48584 1.94690 0.39848 1.57394 8.40315 0 .1 6 1 3 0

Wartość statystyki testowej (14) wynosi 7 = 0,58185. Przyjmując poziom istotności a = 0,05, mamy dwie wartości krytyczne: 70j025:72 = —2,644 i

*0,975 72 = 1,186. Ponieważ spełniona jest nierówność — 2,644 ^ 7 =

= 0,58185 < 1,186, więc na poziomie istotności a = 0,05 nie mamy podstaw do odrzucenia hipotezy o rozkładzie normalnym błędów losowych dla du- żych i małych poletek.

Literatura

H. Brzeskw iniew icz, Częściowo zrównoważone układy bloków niekompletnych dla doświadczeń czynnikowych, Dziesiąte Colloquium Metodologiczne z Agro-Biometrii (1980), 81-97.

R. B. D’A gostino, An omnibus test of normality for moderate and large size samples, Biometrika 58 (1971), 341-348.

J. Durbin, Some methods of constructing exact test, Biometrika 48 (1961), 41-55.

Z. H ellwig, Maszyny cyfrowe i ich zastosowanie, PWE, Warszawa 1975.

H. M ikos, E. N iedokos, Estymacja parametrów w układach eksperymentalnych rozszczepionych jednostek, Szóste Colloquium Metodologiczne z Agro-Biometrii (1976), 64-92.

C. R. Rao, Modele liniowe stystystyki matematycznej, PWN, Warszawa 1982.

W. Wagner, Testy zgodności z rozkładem normalnym dla próby prostej, Listy Biometryczne 77 (1982), 6-37.

W. Wagner, P. Błażczak, Statystyka matematyczna z elementami doświadczalnictwa, Skrypty AR w Poznaniu, 1986.

S. Yam am oto, Y. Fuj ii, N. Ha mad a, Composition of some series of association algebras, J.

Sci. Hiroshima Univ. Ser. A-I (1965), 181-215.

R. Z ieliński, Generatory liczb losowych. Programowanie i testowanie na maszynach cyfrowych, WNT, Warszawa 1972.