x GEJZA WIMMER BIOŠTATISTIKA MODEL S NÁHODNÝMI EFEKTMI 1 ÚVOD Semestrálna prednáška ,,Bioštatistika obsahuje základy teórie odhadu v modeli s náhodnými efektmi a v zmiešanom lineárnom modeli. Oba modely vznikajú aj pri riešení problémov v biológii, medicíne, genetike, poľnohospodárstve atď. Cieľom tejto prednášky je ukázať na príkladoch genézu modelu s náhodnými efektmi ako aj zmiešaného lineárneho modelu, odhady parametrov, testy o parametroch, súvis s modelom analýzy rozptylu a kovariančnou analýzou. Jej súčasťou sú aj dôležité vety týkajúce sa kvadratických štatistík. 2 MODEL ANALÝZY ROZPTYLU MODEL S PEVNÝMI EFEKTMI Pri mnohých pokusoch (vyšetreniach) sa meria (sleduje, vyšetruje, určuje) hodnota kvantitatívnej náhodnej premennej (znaku) na štatistických jednotkách, na ktoré pôsobia rôzne faktory A, B, C, . . . , ktorých účinky nás zaujímajú. Každý faktor vystupuje na rôznych kvalitatívnych alebo kvantitatívnych úrovniach A1, A2, . . . , AI, B1, B2, . . . , BJ atď. Podmienky pri ktorých sa merania uskutočňujú, a ktorých vplyv na výsledok pokusu nás zaujíma, môžu byť vyjadrené v termínoch úrovní určitých faktorov. Faktormi môžu byť aj javy, ktoré z experimentu nemôžeme vylúčiť, ale môžeme od nich očakávať vplyv na výsledok experimentu, hoci tento vplyv na výsledok nie je predmetom nášho výskumu. Predpokladáme, že každá úroveň faktora spôsobuje na výsledok pokusu pevný (aj keď nie známy) efekt. Pokusy robíme pri rôznych známych úrovniach jednotlivých faktorov. Vznikajú otázky: Vplývajú faktory na výsledok pokusu? Ako ohodnotíme (odhadneme) vplyv jednotlivých faktorov na výsledok pokusu? Zaujímajú nás všetky úrovne faktorov, ktoré s meraním nejakým spôsobom súvisia. Jednotlivé merania sú nezávislé. Ďalej predpokladáme, že výsledky merania sú normálne rozdelené, pričom vplyv určitej úrovne niektorého faktora je vyjadrený pevným efektom ­ (pevným) posunom strednej hodnoty merania. Príklad 1 ([7], str. 316) Skúmajme vplyv rôznych druhov penicilínu na bacil B (substilis). Vyšetroval sa vplyv štyroch druhov penicilínu na rast bacilu B. Faktor A (druh penicilínu) má 4 3 úrovne. Na každej úrovni sa päťkrát nezávisle merala veľkosť bacila pod vplyvom úrovne faktora A. Výsledky merania sú v tabuľke 1. Tabuka 1 číslo druh penicilínu ­ i pokusu j 1 2 3 4 1 10,6 7,3 8,2 7,5 2 8,5 9,1 7,7 6,6 3 9,8 8,4 8,0 5,1 4 8,3 8,8 7,2 7,1 5 8,1 7,6 6,4 6,7 Pokus z príkladu 1 modelujeme nasledujúcim spôsobom. Výsledok j-teho pokusu (j = 1, 2, . . . , 5) pri pôsobení i-tej úrovne faktora A (i-teho druhu penicilínu) (i = 1, 2, 3, 4) je yij . Je to realizácia náhodnej premennej Yij , pričom (1) Yij = ij + ij , i = 1, 2, 3, 4 , j = 1, 2, 3, 4, 5 , kde ij N(0, 2 ); ij a 2 sú neznáme konštanty. Pretože predpokladáme, že každý druh penicilínu (úroveň faktora A) má na výsledok rovnaký (pevný) efekt, je pre i = 1, 2, 3, 4 ij = i, j = 1, 2, 3, 4, 5, teda Yij = i + ij, i = 1, 2, 3, 4 , j = 1, 2, 3, 4, 5 . Ak označíme 1 4 4 i=1 i = ; i - = i, i = 1, 2, 3, 4, môžeme (1) písať ako (2) Yij = + i + ij, i = 1, 2, 3, 4 , j = 1, 2, 3, 4, 5 , kde ij N(0, 2 ) sú nezávislé a modelujú náhodnú chybu i, j-teho výsledku (merania), i je (pevný) efekt i-tej úrovne Ai faktora A, je celková stredná hodnota. Model (2) nazývame jednofaktorový vyvážený model analýzy rozptylu (model s pevnými efektmi). Vyváženosť znamená, že pre každú úroveň faktora A bol vykonaný rovnaký počet meraní. Analýzou rozptylu testujeme hypotézu H0 : 1 = 2 = 3 = 4 4 týkajúcu sa zhody jednotlivých úrovní faktora A, alebo (čo je to isté) H0 : 1 = 2 = 3 = 4 = 0 týkajúcu sa zhody efektov, resp. významnosti efektov (t.j. či je niektorý z efektov štatisticky významne rôzny od nuly). V prípade zamietnutia nulovej hypotézy, metódami analýzy rozptylu odpovieme aj na otázku experimentátora: ,,Ktoré skupiny sa od seba líšia? (V príklade 1 ,,Ktoré druhy penicilínu sa od seba líšia vzhľadom na rast bacilu B? ). Všeobecné riešenie pozri napr. v [2] alebo [14]. Príklad 2 ([1], str. 49) Analyzujeme problém, či vplýva u detí vek a nezávisle od toho telesná výška na vitálnu kapacitu pľúc. Zároveň sa vyskytuje aj otázka, či vplyv výšky v spojitosti s vekom je vždy rovnaký, alebo sa mení. Označme yijk vitálnu kapacitu pľúc nameranú na k-tom náhodnom dieťati s telesnou výškou na i-tej úrovni (Ai) a vekom na j-tej úrovni (Bj ). V tabuľke 2 sú v okienku patriacom Ai, Bj (i, j = 1, 2, 3, 4) vždy dve čísla. Dolné je počet detí nij meraných pri daných úrovniach faktorov Ai a Bj a horné je vitálna kapacita nij k=1 yijk v dm3 . Ďalej vieme, že 4 i=1 4 j=1 nij k=1 y2 ijk = 2561, 4759. Tabuka 2 telesná vek v rokoch výška v cm 12­13 13­14 14­15 15­16 140­149 23,2 6,8 12,6 2,7 11 3 5 1 150­159 55,98 102,2 112,75 61,15 23 42 46 23 160­169 28,3 100,3 158,35 160,45 10 36 56 53 170­179 3,1 8,8 24,3 50,5 1 3 8 15 Nameraný výsledok yijk považujeme za realizáciu náhodnej premennej Yijk , pri- čom (3) Yijk = ij + ijk, i = 1, 2, 3, 4 , j = 1, 2, 3, 4 , k = 1, 2, . . . , nij , ijk N(0, 2 ) sú nezávislé a ij sú konštanty. Model (3) píšeme vo forme (4) Yijk = + i + j + ij + ijk, i = 1, 2, 3, 4 , j = 1, 2, 3, 4 , k = 1, 2, . . . , nij , 5 ktorý je preparametrizovaný, i je efekt i-tej úrovne Ai faktora A, j je efekt j-tej úrovne Bj faktora B a ij je tzv. interakcia t.j. v tomto prípade efekt spojitosti i-tej úrovne výšky s j-tou úrovňou veku. Model (4) je dvojfaktorový nevyvážený model analýzy rozptylu s interakciami. Metódami analýzy rozptylu testujeme v modeli (4) hypotézy o nulovosti interakcií, nulovosti i pre všetky i (vplyv faktora A na výsledok), nulovosti j pre všetky j (vplyv faktora B). Taktiež sa dá odpovedať na otázku, ktoré skupiny sa od seba líšia vzhľadom na faktor A, resp. vzhľadom na faktor B. 3 MODEL S NÁHODNÝMI EFEKTMI V medicíne, ale aj v iných vedeckých, ekonomických a technických disciplínach vzniká aj iný typ problému. Mnohé (napr. fyziologické) veličiny sú premenlivé (nie sú konštantné). Túto premenlivosť treba vziať do úvahy, keď porovnávame hodnoty takejto veličiny namerané na rôznych indivíduách (inter-individuálne kolísanie), alebo medzi skupinami (intra-individuálne kolísanie). Aby sme mohli vyjadriť toto kolísanie, opakujeme meranie na tom istom indivíduu, preto sú tieto merania závislé. Nejde nám o odhad pevných efektov, ale o zistenie ich variability. Príklad 3 ([1], str. 80) U desiatich potkaních samcov sa spočítal počet kapilár na 1000 svalových vláknach v troch rôznych rezoch na svale gastrocnemius. Zaujíma nás variabilita jednotlivých meraní u toho istého zvieraťa a variabilita hodnôt medzi zvieratami. Výsledky meraní sú v tabuľke 3. Tabuka 3 číslo číslo pokusného zvieraťa ­ i rezu j 1 2 3 4 5 6 7 8 9 10 1 1049 1152 1166 1534 1368 1300 1589 1223 1254 2 955 1208 1068 1252 1295 1443 1314 1377 1277 3 1242 1127 1532 1507 1424 1328 1231 1298 1588 Ak označíme yij výsledok merania počtu kapilár na 1000 svalových vláknach u i-teho potkana (i = 1, 2, . . . , 10) v j-tom reze (j = 1, 2, 3), tak túto hodnotu považujeme za realizáciu náhodnej premennej (5) Yij = + Ai + ij, i = 1, 2, . . . , 10 , j = 1, 2, 3 , kde je (konštantná) stredná hodnota (pre všetky merania), Ai je náhodná premenná vyjadrujúca premenlivosť hodnôt medzi zvieratami a ij vyjadruje variabilitu jednotlivých hodnôt u toho istého zvieraťa. Predpokladáme, že stredná hodnota (6.a) E(Ai) = 0, i = 1, 2, . . . , 10, 6 ďalej (6.b) E(ij) = 0, i = 1, 2, . . . , 10, j = 1, 2, 3, disperzia (6.c) D(Ai) = 2 A, i = 1, 2, . . . , 10, (6.d) D(ij) = 2 0 i = 1, 2, . . . , 10, j = 1, 2, 3, kovariancia (6.e) cov(Ak, Al) = 0 pre všetky k = l, (6.f) cov(Ai, jk) = 0 pre všetky i, j, k, (6.g) cov(ij, kl) = 0 pre všetky i = k, (6.h) cov(ij, il) = 0 pre všetky j = l. Platí (7) D(Yij) = 2 A + 2 0. Úlohou je odhadnúť 2 A a 2 0 , ďalej testovať, či 2 A = 0. Model (5) s podmienkami (6.a) až (6.h) sa nazýva model s jedným náhodným efektom. 2 A a 2 0 sú variančné komponenty jednotlivého merania (vidieť to z (7)). Príklad 4 ([1], str. 83) Na detskej klinike v Halle skúmali okamžitú adaptáciu detí nyktometrom (zisťuje sa počet identifikovaných plôšok pod určitou krivkou). Vykonalo sa niekoľko meraní v 4 rôznych obdobiach roka u 5 náhodne vybratých detí navštevujúcich druhú triedu. Výsledky sú v tabuľke 4. Zaujímajú nás variančné komponenty. Ak označme yijk hodnotu registrovanú nyktometrom u i-teho žiaka (i = 1, 2, 3, 4, 5) v j-tom termíne pri k-tom vyšetrení (k = 1, 2, . . . , nij ), tak v tabuľke 4 v (i, j)tom okienku sú dve čísla. Horné je počet bodov, ktoré spolu získal i-ty žiak v jtom termíne a dolné je nij (počet vyšetrení i-teho žiaka v j-tom termíne). Ďalej 5 i=1 4 j=1 nij k=1 y2 ijk = 143033, 5. 7 Tabuka 4 číslo číslo vyšetrenia ­ j žiaka i 1 2 3 4 1 177 131 255,5 246 4 4 5 5 2 119,5 109 173,5 167,5 4 4 5 5 3 167 127 215 166 4 4 5 5 4 142,5 153 256 171,5 4 4 6 4 5 155,5 131 269,5 178 4 4 6 4 Hodnotu yijk považujeme za realizáciu náhodnej premennej Yijk , pričom predpokladáme, že (8) Yijk = + Ai + Bj + Cij + ijk, i = 1, 2, 3, 4, 5 , j = 1, 2, 3, 4 , k = 1, 2, . . . , nij , je konštanta a (9.a) E(Ai) = E(Bj) = E(Cij) = E(ijk) = 0 pre všetky i, j, k, (9.b) D(Ai) = 2 A, i = 1, 2, 3, 4, 5, (9.c) D(Bj) = 2 B, j = 1, 2, 3, 4, (9.d) D(Cij) = 2 AB, pre všetky i, j, (9.e) D(ijk) = 2 0, pre všetky i, j, k. Náhodné premenné Ai, Bj, Cij, ijk sú medzi sebou nekorelované (pre všetky i, j, k). 8 Ai predstavujú (náhodný) efekt žiaka, Bj predstavujú (náhodný) efekt termínu a Cij predstavujú (náhodný) efekt interakcie žiaka s termínom. Platí D(Yijk) = 2 A + 2 B + 2 AB + 2 0, kde 2 A je miera variability medzi žiakmi, 2 B je miera variability medzi termínmi, 2 AB je miera variability interakcie (vedľajšieho účinku žiaka s termínom), 2 0 je miera variability nameraných hodnôt jednotlivých vyšetrení. Model (8) s podmienkami (9.a) až (9.e) nazývame modelom s dvoma náhodnými efektmi a interakciou. Poďme si ešte raz analyzovať rozdiel medzi modelmi s pevnými (analýza rozptylu) a náhodnými efektmi. Príklad 5 ([7], str. 307) Presnosť chromatografického určenia dietylénglykolu. V laboratóriu pracujú traja laboranti na troch chromatografoch a rutinne určujú obsah (v %) dietylénglykolu v etylénglykole. Chceme zistiť, či na výsledok analýzy má vplyv prístroj, či laborant a či vplyv súčinnosti laboranta a prístroja je rovnaká. Na každom prístroji B1, B2, B3 vykonal každý laborant A1, A2, A3 dve opakované merania (meral sa ten istý roztok). Namerané hodnoty percentuálneho obsahu dietylénglykolu sú v tabuľke 5. Tabuka 5 prístroj laborant B1 B2 B3 A1 0,110 0,116 0,101 0,102 0,108 0,109 A2 0,112 0,111 0,115 0,106 0,111 0,109 A3 0,114 0,112 0,107 0,109 0,113 0,110 Ak označíme yijk výsledok i-teho laboranta na j-tom prístroji pri k-tom určení (i, j = 1, 2, 3, k = 1, 2), tak yijk je realizácia náhodnej premennej Yijk , pričom Yijk = + i + j + ij + ijk, i, j = 1, 2, 3 , k = 1, 2 , i je (pevný) efekt i-teho laboranta, j je (pevný) efekt j-teho prístroja a ij je (pevný) efekt súčinnosti i-teho laboranta s j-tym prístrojom, ijk N(0, 2 0) považujeme za nezávislé (v tomto kroku sa vedome dopúšťame ,,deformácie modelu oproti skutočnosti). Prezentovaný model je modelom dvojfaktorovej analýzy rozptylu s interakciami (vyvážený model s pevnými efektmi). Príklad 6 ([7], str. 311) Vplyvy pôsobiace na výsledok chromatografického určenia dietylénglykolu . 9 V laboratóriu pracuje veľa laborantov a majú veľa prístrojov. Náhodne sa vybrali traja laboranti A1, A2, A3 a tri prístroje B1, B2, B3. Výsledky uvažujeme rovnaké ako v tabuľke 5. Máme určiť variabilitu laborantov, prístrojov a interakcií na výslednú variabilitu výsledkov. Opäť označíme yijk výsledok i-teho náhodne vybratého laboranta na j-tom náhodne vybratom prístroji pri k-tom určení (i, j = 1, 2, 3, k = 1, 2). Potom yijk považujeme za realizáciu náhodnej premennej Yijk , pričom Yijk = + Ai + Bj + Cij + ijk, i, j = 1, 2, 3 , k = 1, 2 , je konštanta, Ai, Bj, Cij, ijl sú nekorelované náhodné premenné s nulovými strednými hodnotami, D(Ai) = 2 A, D(Bj) = 2 B, D(Cij) = 2 AB, D(ijk) = 2 0, pre všetky i, j, k. Úlohu modelujeme dvojfaktorovým modelom s náhodnými efektmi a interakciou. Ešte na jednej dvojici príkladov si ukážme rozdiel medzi modelom analýzy rozptylu (model s pevnými efektmi) a modelom s náhodnými efektmi. Príklad 7 ([7], str. 313) Máme 5 laboratórií a 2 metódy stanovenia arzénika v potrave (metóda A­reakcia s molybdénovým roztokom, metóda B­reakcia s dietylditiokarbamátom strieborným podľa Vašáka a Šedivca). K vzorke potravy pridáme 15 g arzénika a vzorku nezávisle analyzujeme oboma metódami vo všetkých piatich laboratóriach. Vedú obe metódy vo všetkých laboratóriach k rovnakým výsledkom? Údaje sú v tabuľke 6. Tabuka 6 laboratórium metóda obsah arzénika v g A 12,9 13,2 12,9 12,9 13,1 13,0 1 B 13,4 13,0 13,0 17,0 16,6 A 14,6 16,2 14,0 15,0 15,5 13,7 2 B 14,8 15,2 14,6 15,0 A 13,4 13,0 13,2 13,2 13,1 13,2 3 B 14,8 14,8 15,0 14,9 14,8 A 13,3 13,8 12,5 13,5 13,6 12,8 4 B 14,8 14,8 15,0 14,5 15,4 15,2 A 15,9 14,8 15,3 15,6 14,9 15,2 5 B 13,8 14,1 13,8 13,9 14,0 10 Úlohu riešime modelom dvojfaktorovej analýzy rozptylu s interakciami Yijk = + i + j + ij + ijk, i = 1, 2, 3, 4, 5 , j = 1, 2 , k = 1, 2, . . . , nij . Hodnota yijk je nameraná hodnota v i-tom laboratóriu pri j-tej metóde k-ty krát. Je to realizácia náhodnej premennej Yijk ; i je (pevný) efekt i-teho laboratória, j je (pevný) efekt j-tej metódy, ij je (pevná) interakcia j-tej metódy v i-tom laboratóriu; ijk N(0, 2 0) sú nezávislé náhodné premenné (pre všetky i, j, k). Príklad 8 ([7], str. 313, modifikovaný) Máme veľa laboratórií a veľa metód na určenie arzénika v potrave. Náhodne vyberieme 5 laboratórií a 2 metódy. Meranie opakujeme niekoľkokrát. Výsledky sú v tabuľke 6. Zaujíma nás variabilita jednotlivých meraní (v tom istom laboratóriu), variabilita medzi laboratóriami, variabilita medzi metódami, variabilita náhodného efektu interakcie laboratória a metódy. Úlohu riešime dvojfaktorovým modelom s náhodnými efektmi Yijk = + Ai + Bj + Cij + ijk, i = 1, 2, 3, 4, 5 , j = 1, 2 , k = 1, 2, . . . , nij . Ai predstavuje (náhodný) efekt laboratória, Bj predstavuje (náhodný) efekt metódy a Cij predstavuje (náhodný) efekt interakcie laboratória a metódy. Platia vzťahy ako v modeli (8) s podmienkami (9.a) až (9.e). 4 MODEL SO ZMIEŠANÝMI EFEKTMI Príklad 9 ([7], str. 311, modifikovaný) Uvažujme rovnaké zadanie ako v príklade 5 s tým rozdielom, že z laborantov náhodne vyberieme troch. Faktor A (laboranti) je faktor s náhodnými efektmi a faktor B (prístroja) je faktor s pevnými efektmi. Príklad riešime modelom (10) Yijk = + j + Ai + Cij + ijk, i = 1, 2, 3 , j, k = 1, 2 , kde je (neznáma) konštanta, j je (pevný, neznámy) efekt j-teho prístroja, Ai je (náhodný) efekt i-teho laboranta, Cij je (náhodný) efekt interakcie (súčinnosti) iteho laboranta s j-tym prístrojom, ijk N(0, 2 0) sú nezávislé náhodné premenné. Ďalej predpokladáme, že Ai N(0, 2 A) sú medzi sebou navzájom nezávislé aj nezávislé so všetkými ijk . Cij N(0, 2 AB) sú medzi sebou nezávislé aj nezávislé 11 s Ai, ijk pre všetky i, j, k. Maticovo môžeme model (10) zapísať ako Y18,1 = Y111 Y112 Y121 Y122 Y131 Y132 Y211 Y212 Y221 Y222 Y231 Y232 Y311 Y312 Y321 Y322 Y331 Y332 = 1 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 1 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 1 0 0 0 0 0 0 1 1 0 0 0 1 0 0 0 0 1 0 0 0 0 0 1 1 0 0 0 1 0 0 0 0 1 0 0 0 0 0 1 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 1 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 1 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 1 0 0 1 0 1 0 0 0 0 0 0 1 0 0 0 1 1 0 0 0 0 1 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 1 0 0 0 0 0 0 1 0 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1 1 2 3 A1 A2 A3 C11 C12 C13 C21 C22 C23 C31 C32 C33 + 18,1, čiže (11) Y = X + Z + = (X ...Z ) + , kde = (, 1, 2, 3) je vektor pevných (neznámych) parametrov, = (A1, A2, A3, C11, C12, . . . , C33) je vektor náhodných parametrov (efektov). Matica plánu (X ...Z ) je známa. Niekedy je užitočné rozdeliť maticu plánu na submatice podľa jednotlivých ,,skupín efektov. V našom prípade dostávame model (12) Y = (X1 ...X2) + (Z1 ...Z2) + , kde matica X1 = 1 = (1, 1, . . . 1) je rozmeru 18 × 1 a prislúcha konštantnému (pevnému) parametru , matica X2 = 1 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 1 1 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 1 1 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 1 12 prislúcha pevnému vektoru (1, 2, 3) (efekt prístroja), matica Z1 = 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 prislúcha vektoru náhodných efektov (A1, A2, A3) (efekt laboranta) a matica Z2 = 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 prislúcha vektoru náhodných efektov (C11, C12, . . . , C33) (efekt interakcie). Vektor Y má strednú hodnotu (13) E(Y ) = X = (X1 ...X2) 1 2 3 a kovariančnú maticu (14) cov(Y ) = 2 AZ1Z1 + 2 ABZ2Z2 + 2 0I. Model (10) (resp. (11) alebo (12)), (pričom observovaný vektor má strednú hodnotu (13) a kovariančnú maticu (14)) patrí medzi modely so zmiešanými efektmi (alebo medzi zmiešané lineárne modely). 13 5 ODHADY V ZMIEŠANOM LINEÁRNOM MODELI PODĽA HENDERSONA V mnohých biologických, ale aj priemyselných, ekonomických a iných situáciach sa dostávame k modelu (15) Y = X + Z1b1 + Z2b2 + . . . + Zsbs + = (X ...Z ) b + , kde náhodný observačný vektor Y Rn , známa matica plánu (X ...Z ) = (X ...Z1 ...Z2 ... . . . ...Zs) je rozmerov n×k. Matica X rozmerov n×q ,,prislúcha pevným efektom (neznámym parametrom) Rq . Predpokladáme, že model (15) má s náhodných efektov, teda b = (b1, b2, . . . , bs) je vektor náhodných efektov (náhodných premenných), pričom bi Rmi , i = 1, 2, . . . , s, E(bi) = 0, i = 1, 2, . . . , s, cov(b) = cov( b1 b2 ... bs ) = 2 1Im1,m1 0 . . . 0 0 2 2Im2,m2 . . . 0 ... ... ... 0 0 . . . 2 s Ims,ms . Matice Z1, Z2, . . . , Zs (známe) ,,prislúchajúce týmto náhodným efektom sú rozmerov n×mi, i = 1, 2, . . . , s. Chybový vektor Rn má E() = 0 a var() = 2 0I. Predpokladáme, že nie je skorelovaný s náhodným vektorom b. Zrejme E(Y ) = X a cov(Y ) = =(Z1 ... . . . ...Zs ...I) 2 1Im1,m1 0 . . . 0 0 0 2 2Im2,m2 . . . 0 0 ... ... ... 0 0 2 s Ims,ms 0 0 0 . . . 0 2 0In,n Z1 ... Zs I = = 2 1Z1Z1 + 2 2Z2Z2 + . . . + 2 s ZsZs + 2 0I = s i=1 2 i Vi + 2 0I. V nasledujúcom budeme označovať b0 = , m0 = n, Z0 = In,n. Potom cov(Y ) = s i=0 2 i ZiZi. Model (15) s horeuvedenými podmienkami sa nazýva model s variančnými komponentmi alebo zmiešaný lineárny model. Je špeciálnym prípadom všeobecného lineárneho modelu s kovariančnými komponentmi v ktorom sa predpokladá cov(Y ) = s i=1 iVi, 1, . . . , s (otvorená v Rs ), Vi = Vi , i = 1, 2, . . . , s a matica s i-1 iVi je pozitívne semidefinitná. 14 Lema 1. V modeli (15) je p (p-daný vektor) lineárne nevychýlene odhadnuteľná funkcia (pevných) parametrov práve vtedy ak p (X ) = {X u: u Rn }. Dôkaz. p je lineárne nevychýlene odhadnuteľná l Rn : E(l Y ) = p Rq l Rn : l X = p Rq l X = p p (X ) = {X u: u Rn }. Henderson v [4] navrhol v modeli (15) na základe realizácie yn,1 náhodného (observačného) vektora Y 1. nevychýlený lineárny odhad lineárne odhadnuteľnej funkcie pevných parametrov (efektov) p , 2. lineárny prediktor vektora náhodných parametrov (b1, b2, . . . , bs) , 3. nevychýlené kvadratické odhady variančných komponentov 2 1, 2 2, . . . , 2 s a 2 0 . Definícia 2. Majme maticu A rozmerov m × n. Matica Arozmerov n × m pre ktorú platí (16) AAA = A je g - inverzia (pseudoinverzia) matice A. Pseudoinverzia existuje ku každej matici, nie je vzťahom (16) určená jednoznačne, bližšie pozri [2], str. 67, [10] a inde. Lema 3. BAA = B (B ) (A ) t.j. D, že B = DA. Dôkaz. Uvedomme si najprv, že (B ) (A ) D : B = A D D: B = DA. Ak (B ) (A ), tak B = A D = B = DA = BAA = DAAA = DA = B. Naopak, ak BAA = B, tak D(= BA): DA = B = (B ) (A ). Poznámka 4. Lema 3 platí aj v tvare B = AAB B = AD (B) (A). Lema 5. (A A ) = (A). Dôkaz. Vyplýva napr. z vety 1 v [2], str. 62. Lema 6. Ak Z = (Z1 ...Z2 ... . . . ...Zs) a S = X Z (X ...Z) = X X X Z Z X Z Z = A C C B , tak G= (X X)+(X X)X Z DZ X(X X)-(X X)X Z D- -DZ X(X X)D- , 15 je jedna g­inverzia matice S, pričom D = Z Z - Z X(X X)X Z. Dôkaz. Počítajme SGS = A C C B A+ A- CDC A- -A- CD- -DC A- DA C C B = = AA+ AA- CDC A- CDC A-... - AA- CD+ CDC A+ C A- CDC A- BDC A-... - C A- CD+ BDA C C B = = A + AA- CDC AA - CDC AA - AA- CDC + CD- C ... C AA + C A- CDC AA - BDC AA - C A- CDC + BD- C ... ... AAC + AA- CDC AC - CDC AC - AA- CDB + CD- B ... C AC + C A- CDC AC - BDC AC - C A- CDB + BD- B . Pre jednotlivé bloky matice SGS platí {SGS}11: A + AA- CDC AA - CDC AA - AA- CDC + CDC = = X X + X X(X X)X ZDC AA - X ZDC A- A-X X(X X)X Z DC + X ZDC = X X = A (lebo podľa lemy 3 je X X(X X)X = X , čiže AAC = C). {SGS}12: AA- C+AA- CDC A- C-CDC A- C-AA- CD- B+CDB =C =X Z. Podobne {SGS}21 = C = Z X. Konečne pre {SGS}22: C AC + C A- CDC AC - BDC AC - C A- CDB + BDB = = C AC + {C AC - B}{-C AC + B}C A- C+ + {B - C AC}{B - C A- C}B = = C AC + {-C AC + B}{-C AC + B}{-C AC + B} = = C AC - C AC + B = Z Z. 16 Poznámka 7. (i) Pre ľubovoľnú maticu A je [(A A)] tiež g-inverzia matice A A, preto 1 2 {(A A)+[(A A)] } je symetrická g-inverzia matice A A (bez ohľadu na výber (A A)- ). (ii) Pre ľubovoľnú maticu A matica A(A A)A nezávisí na voľbe (A A)a A(A A)A je vždy symetrická. Dôkazy (i) a (ii) nájdete v [2], str. 69. (iii) Na základe (i) a (ii) vždy existuje symetrická g-inverzia G v leme 6. Lema 8. Rovnica (17) X X X Z Z X Z Z ^ ^b1 ... ^bs = X Z Y je vždy riešiteľná (vzhľadom na ^, ^b1, . . . , ^bs). Jedno jej riešenie je (18) ^ ^b1 ... ^bs = G X Z Y , kde G je z lemy 6 a symetrická (existuje podľa poznámky 7). Pre toto riešenie platí: 1. p ^ je nevychýlený odhad p pre p (X ). 2. E ^b1 ... ^bs = 0. Dôkaz. Pretože X Z (X ...Z) = X Z , je (17) vždy riešiteľná. Jedno riešenie je (18) (pozri napr. [2], str. 70; všeobecné riešenie pozri tamtiež). 1. Počítajme pre p (X ) (teda p = X u): Rq E(p ^) = E(u X{(X X)X Y + + (X X)X Z DZ X(X X)X Y - (X X)X Z DZ Y }) = =u X(X X)X X+u X(X X)X Z DZ X(X X)X X- u X(X X)X Z DZ X = u X = p (využijúc lemu 5 a lemu 6). 2. Platí E ^b1 ... ^bs = -DZ X(X X)X X + DZ X = 0. 17 Poznámka 9. Odhad p ^ a predikciu ( ^b 1, ^b 2, . . . , ^b s) z lemy 8 navrhol Henderson v [4]. Lema 10. Ak pre náhodné vektory n,1 a m,1 platí E() = , E() = a P je matica n × m, tak E( P ) = P + tr P [cov(, )] . (tr P znamená stopu matice P ). Dôkaz. E( P ) = E([ - + ] P [ - + ]) = E([ - ] P [ - ]) + P = = P + E tr([ - ] P [ - ]) = P + E tr( - )( - ) P = = P+Etr P(-)(-) = P+tr PE(-)(-) = P+tr P[cov(, )] . Pri hľadaní odhadu variančných komponentov postupujeme zaužívanou cestou ako pri analýze rozptylu. Ukážeme si to pre prípad s = 3 (3 náhodné efekty plus ,,efekt chýb), teda (19) Y = X + Z1b1 + Z2b2 + Z3b3 + . Formálne predpokladáme, že aj b1, b2 a b3 sú pevné (nenáhodné) parametre a cov(Y ) = 2 0I. Uvažujme postupnosť lineárnych modelov M : E(Y ) = (X ...Z1 ...Z2 ...Z3) b1 b2 b3 , cov(Y ) = 2 0I; M1 : E(Y ) = (X ...Z1 ...Z2) b1 b2 , cov(Y ) = 2 0I; M2 : E(Y ) = (X ...Z1) b1 , cov(Y ) = 2 0I; M3 : E(Y ) = X, cov(Y ) = 2 0I. Je zrejmé, že M3 je submodelom M2 (lebo X = (X ...Z1) I 0 ), podobne M2 je submodelom M1 (lebo (X ...Z1) = (X ...Z1 ...Z2) I 0 0 I 0 0 ) a tiež M1 je submodelom M. Predpokladajme o hodnostiach (20) n > h(X ...Z1 ...Z2 ...Z3) > h(X ...Z1 ...Z2) > h(X ...Z1) > h(X) > 0. 18 Pri označení z (15) položme ^ = (X ...Z) X X X Z Z X Z Z - X Z Y , ^ = (X ...Z1 ...Z2) X X X Z1 X Z2 Z1X Z1Z1 Z1Z2 Z2X Z2Z1 Z2Z2 - X Z1 Z2 Y , ^ = (X ...Z1) X X X Z1 Z1X Z1Z1 - X Z1 Y , ^ = X(X X)X Y . Podľa [2], str. 138 ^ je BLUE (best linear unbiased estimator) (NNLO - najlepší nevychýlený lineárny odhad) E(Y ) v modeli M, podobne ^ je BLUE E(Y ) v modeli M1 atď. Se označme reziduálny súčet štvorcov v modeli M, t.j. Se = (Y - ^) (Y - ^). Veta 11. Veličiny ^, ^, ^ a ^ spĺňajú identitu (22) Se + (^ - ^) (^ - ^) + (^ - ^) (^ - ^) + (^ - ^) (^ - ^) + ^ ^ = Y Y . Dôkaz. Platí (23) Se = (Y - ^) (Y - ^) = Y Y - Y (X ...Z) X X X Z Z X Z Z - X Z Y = = Y Y - ^ ^. Ďalej (^ - ^) (^ - ^) = = ^ ^ - Y (X ...Z1 ...Z2) X X X Z1 X Z2 Z1X Z1Z1 Z1Z2 Z2X Z2Z1 Z2Z2 - X Z1 Z2 × × (X ...Z) X X X Z Z X Z Z - X Z Y - Y (X ...Z) X X X Z Z X Z Z - X Z × × (X ...Z1 ...Z2) X X X Z1 X Z2 Z1X Z1Z1 Z1Z2 Z2X Z2Z1 Z2Z2 - X Z1 Z2 Y + ^ ^ = 19 = ^ ^ + ^ ^ - Y (X ...Z1 ...Z2) X X X Z1 X Z2 Z1X Z1Z1 Z1Z2 Z2X Z2Z1 Z2Z2 - × × Im1,m1 0 0 0m1,m4 0 Im2,m2 0 0m2,m4 0 0 Im3,m3 0m3,m4 X Z × × (X ...Z) X X X Z Z X Z Z - X Z Y - (24) - Y (X ...Z) X X X Z Z X Z Z - X Z (X ...Z) Im1,m1 0 0 0 Im2,m2 0 0 0 Im3,m3 0m4,m1 0m4,m2 0m4,m3 × × X X X Z1 X Z2 Z1X Z1Z1 Z1Z2 Z2X Z2Z1 Z2Z2 - X Z1 Z2 Y = ^ ^ - ^ ^. Analogicky (25) (^ - ^) (^ - ^) = ^ ^ - ^ ^ a (26) (^ - ^) (^ - ^) = ^ ^ - ^ ^. Použitím (23)­(26) dostávame Se + (^ - ^) (^ - ^) + (^ - ^) (^ - ^) + (^ - ^) (^ - ^) + ^ ^ = Y Y . Označme (27) ^ ^ = Y (I - D1)Y = S1, kde D1 = I - X(X X)X , ďalej (^ - ^) (^ - ^) = ^ ^ - ^ ^ = = Y (X ...Z1) X X X Z1 Z1X Z1Z1 - X Z1 Y - Y (I - D1)Y = = Y (X ...Z1)× × (X X)+(X X)X Z1DZ1X(X X)-(X X)X Z1D- -DZ1X(X X)D- × × X Z1 Y - Y (I - D1)Y = 20 = Y X(X X)X Y + Y (I - X(X X)X )Z1[Z1(I - X(X X)X )Z1]- × × Z1(I - X(X X)X )Y - Y (I - D1)Y = (28) = Y (D1 - D12)Y = S2, kde D12 = D1 - D1Z1(Z1D1Z1)- Z1D1. (Využili sme tvar g-inverzie matice S z lemy 6, lebo podľa poznámky 7 nezáleží na tom, ktorú g-inverziu použijeme.) Počítajme ďalej analogicky ako v predchádzajúcom prípade (29) (^ - ^) (^ - ^) = ^ ^ - ^ ^ = = Y (X ...Z1 ...Z2) X X X Z1 X Z2 Z1X Z1Z1 Z1Z2 Z2X Z2Z1 Z2Z2 - X Z1 Z2 Y - Y (X ...Z1) X X X Z1 Z1X Z1Z1 - X Z1 Y = = Y (D12 - D123)Y = S3, kde D123 = D12 - D12Z2(Z2D12Z2)- Z2D12. Podobne dostaneme (30) (^ - ^) (^ - ^) = =Y (D123 - D1234)Y =S4, kde D1234 =D123 -D123Z3(Z3D123Z3)- Z3D123. Zo vzťahov (27)­(30) a (22) vyplýva (31) Se = Y Y - S1 - S2 - S3 - S4 = = Y Y-Y (I-D1)Y-Y (D1-D12)Y-Y (D12-D123)Y-Y (D123-D1234)Y = = Y D1234Y . Takto sme dostali (vhodné) kvadratické formy Se, S1, S2, S3, S4 náhodného vektora Y (toto bol ,,návod k ich dosiahnutiu). Vráťme sa k modelu (15) pre s = 3. Platí X D1 = X D12 = X D123 = X D1234 = 0, Z1D12 = Z1D123 = Z1D1234 = 0, Z2D123 = Z2D1234 = 0 a Z3D1234 = 0. Pretože E(Y ) = X 21 a cov(Y ) = 2 1Z1Z1 + 2 2Z2Z2 + 2 3Z3Z3 + 2 0I, dostávame pomocou lemy 10 (32) E(S2) = E(Y (D1 - D12)Y ) = = tr(D1 - D12)(2 1Z1Z1 + 2 2Z2Z2 + 2 3Z3Z3 + 2 0I) = = 2 1 tr Z1D1Z1 + 2 2{tr Z2D1Z2 - tr Z2D12Z2}+ +2 3{tr Z3D1Z3 - tr Z3D12Z3} + 2 0 tr(D1 - D12). (33) E(S3) = E(Y (D12 - D123)Y ) = = 2 2 tr Z2D12Z2 + 2 3{tr Z3D12Z3 - tr Z3D123Z3} + 2 0 tr(D12 - D123), (34) E(S4) = E(Y (D123 - D1234)Y ) = = 2 3 tr Z3D123Z3 + 2 0 tr(D123 - D1234), (35) E(Se) = 2 0 tr D1234. Z (32) až (35) dostávame (36) E S2 S3 S4 Se = = tr Z1D1Z1 tr Z2(D1-D12)Z2 tr Z3(D1-D12)Z3 tr(D1-D12) 0 tr Z2D12Z2 tr Z3(D12-D123)Z3 tr(D12-D123) 0 0 tr Z3D123Z3 tr(D123-D1234) 0 0 0 tr D1234 2 1 2 2 2 3 2 0 = = E 2 1 2 2 2 3 2 0 . Teraz sa venujme riešiteľnosti systému (36). Opreme sa o známe skutočnosti: (37) Am,n tr A A = m i=1 n j=1 {A}2 ij, (38) Am,n {(A)} = {Au: u Rn } = Ker A , kde je symbol pre ortogonálny doplnok a Ker A = {w Rm : A w = 0}. 22 Ak tr Z1D1Z1 = 0 = tr Z1D1D1Z1 = 0 = = Z1D1 = 0 = Z1(I - X(X X)X ) = 0 = = Ker X (= (I - X(X X)X )) Ker Z1 = = (Z1) (X), čo je spor s predpokladom (20) modelu. Ak tr Z2D12Z2 = 0 = tr Z2D12D12Z2 = 0 = = Z2D12 = 0 = Z2(D1 - D1Z1(Z1D1Z1)Z1D1) = 0 = = Z2(I - X(X X)X - D1Z1(Z1D1Z1)Z1D1) = 0 = = (I - X(X X)X - D1Z1(Z1D1Z1)Z1D1) Ker Z2 = = (Z2) Ker(I - X(X X)X - D1Z1(Z1D1Z1)Z1D1) = = (X(X X)X + D1Z1(Z1D1Z1)Z1D1) = = (X(X X)X + + Z1(Z1D1Z1)Z1D1 - X(X X)X Z1(Z1D1Z1)Z1D1) (X ...Z1), čo je spor s (20). Rovnakým spôsobom ukážeme, že ak tr Z3D123Z3 = 0 = (Z3) (X ...Z1 ...Z2 ...Z3), čo je opäť v spore s (20). Venujme sa ešte výrazu tr D1234 = 1 2 0 E(Se) = 1 2 0 E Y (I - (X ...Z) X X X Z Z X Z Z - X Z Y (pozri (23)). Lema 12. E(Se) = 2 0[n - h(X ...Z)]. Dôkaz. Z = (Z1 ...Z2 ...Z3). Z (23) dostávame E(Se) = E Y I - (X ...Z) X X X Z Z X Z Z - X Z Y = = X I - (X ...Z) X X X Z Z X Z Z - X Z X+ 23 + tr I - (X ...Z) X X X Z Z X Z Z - X Z × × Z 2 1Im1,m1 0 0 0 2 2Im2,m2 0 0 0 2 3Im3,m3 Z + 2 0I = = X X - (X X ...X Z) X X X Z Z X Z Z X X Z X + + tr Z 2 1Im1,m1 0 0 0 2 2Im2,m2 0 0 0 2 3Im3,m3 Z - tr(X ...Z) X X X Z Z X Z Z X Z Z Z 2 1Im1,m1 0 0 0 2 2Im2,m2 0 0 0 2 3Im3,m3 Z + +2 0 tr I - 2 0 tr(X ...Z) X X X Z Z X Z Z - X Z = = X X - X X + tr Z 2 1Im1,m1 0 0 0 2 2Im2,m2 0 0 0 2 3Im3,m3 Z - tr(Z X ...Z Z) X X X Z Z X Z Z X Z Z Z × × 2 1Im1,m1 0 0 0 2 2Im2,m2 0 0 0 2 3Im3,m3 + + n2 0 - 2 0 tr(X ...Z) X X X Z Z X Z Z - X Z = = tr Z Z 2 1Im1,m1 0 0 0 2 2Im2,m2 0 0 0 2 3Im3,m3 - tr Z Z 2 1Im1,m1 0 0 0 2 2Im2,m2 0 0 0 2 3Im3,m3 + + 2 0 n - tr X X X Z Z X Z Z X X X Z Z X Z Z = 2 0 n - h(X ...Z) , lebo pre ľubovolnú maticu H je HH(vždy) idempotentná (napr. pozri [2], str. 66) a preto tr HH= h(HH). Ďalej HHH = H, teda h(H) h(HH- ) h(HHH) = h(H), čiže h(HH) = h(H). 24 Ak je splnená podmienka (20), je tr Z1D1Z1 = 0, tr Z2D12Z2 = 0, tr Z3D123Z3 = 0 aj tr D1234 = n - h(X ...Z) = 0, teda (36) je riešiteľná (E v (36) je nesingulárna). Odhad (39) E-1 S2 S3 S4 Se = 2 1 2 2 2 3 2 0 je nevychýleným odhadom (2 1, 2 2, 2 3, 2 4) . Poznámka 13. Odhad (39) navrhol Henderson v [4], maticovo vyjadril Searle v [12], preformuloval Rohde [11]. Poznámka 14. Z predchádzajúceho je jasné, ako postupovať pri určení odhadu typu (39) pre model (15) s ľubovoľným počtom s náhodných efektov. Skôr ako špecifikujeme odhad (39) pre konkrétne modely, uvedieme bez dôkazov ešte niekoľko základných viet ohľadne rozdelenia kvadratických foriem. Dôkazy viet nájdeme v uvedenej literatúre. Definícia 15. Nech Y1, Y2, . . . , Yn sú nezávislé N(i, 1) rozdelené, = n i=1 2 i . Y = Y 2 1 + Y 2 2 + . . . + Y 2 n má necentrálne 2 (n, ) rozdelenie s n stupňami voľnosti a parametrom necentrality . Ak = 0, Y má centrálne 2 (n, 0) rozdelenie. Veta 16. Nech Y Np(, ). Náhodná premenná Y AY + 2b Y + c má 2 (k, ) rozdelenie práve vtedy ak (i) AA = A, alebo ekvivalentne (A)3 = (A)2 , (ii) (A + b) (A ), (iii) (A + b) (A + b) = A + 2b + c. V takomto prípade k = tr A, = (A + b) A (A + b). Dôkaz nájdeme v [10], str. 171, [6], str. 220, [5], str. 178. Veta 17. Nech Y Np(, ); Q1 = Y A Y , Q2 = Y BY dve kvadratické formy. Q1 a Q2 sú nezávislé práve vtedy ak (i) AB = 0, AB = BA = 0, AB = 0, ak A, B sú symetrické, nemusia byť p.s.d., môže byť singulárna. (ii) AB = 0, AB = 0, 25 ak A je p.s.d. (iii) AB = 0, ak A aj B sú p.s.d. (iv) AB = 0, ak je regulárna, A, B symetrické (nemusia byť p.s.d.). Dôkaz nájdeme v [10], str. 178. Veta 18. Nech Y N(, 2 H), A, B sú symetrické. Potom E(Y A Y Y BY ) = 4 [2 tr(AHBH) + tr(AH) tr(BH)]+ + 2 [ A tr(BH) + B tr(AH) + 4 AHB + A B], cov(a Y + Y AY , b Y + Y BY ) = 24 tr A HBH + 2 (2A + a) H(2B + b). Dôkaz vykonáme pomocou [6], str. 223, [5], str. 181. Urobte ho ako cvičenie. 6 ODHADY VARIANČNÝCH KOMPONENBTOV V MODELI S JEDNÝM NÁHODNÝM EFEKTOM Majme model (40) Yij = + bi + ij, i = 1, 2, . . . , I , j = 1, 2, . . . , ni , kde náhodný vektor b1 = b1 b2 ... bI je N(0I,1; 2 1I) rozdelený. Chybový vektor N(0, 2 0I) je nezávislý od b, je číslo. Model (40) sa dá písať Y = Y11 Y12 ... Y1n1 Y21 Y22 ... Y2n2 ... YI1 YI2 ... YInI = 1 1 0 0 . . . 0 1 1 0 0 . . . 0 ... ... ... ... 1 1 0 0 0 1 0 1 0 0 1 0 1 0 0 ... ... ... ... 1 0 1 0 0 ... ... ... ... 1 0 0 0 1 1 0 0 0 1 ... ... ... ... 1 0 0 0 1 b1 b2 ... bI + = (X ...Z1) b1 + . 26 Zrejme E(Y ) = X = 1, 1 = (1, 1, . . . 1) je n = I i=1 ni-rozmerný vektor a cov(Y ) = 2 1Z1Z1 + 2 0I. Hodnosť matice (X ...Z1) je I. Pri odhade variančných komponentov postupujeme podľa návodu na str. 17. Máme tentoraz postupnosť submodelov M : E(Y ) = (X ...Z1) b1 , cov(Y ) = 2 0I; M1 : E(Y ) = X, cov(Y ) = 2 0I. Platí ^ = (X ...Z1) X X X Z1 Z1X Z1Z1 - X Z1 Y , ^ = X(X X)X Y . Spočítajme S2. Podľa (28) S2 = Y (D1 - D12)Y , kde D1 = I - X(X X)X , D12 = D1 - D1Z1(Z1D1Z1)- Z1D1 (upozorňujeme len, že D1 a D12 nezávisia podľa poznámky 7 od voľby g-inverzie). Platí D1 = I - 1 n 11 , ďalej Z1D1Z1 = Z1Z1 - 1 n Z111 Z1 = = n1 0 . . . 0 0 n2 ... ... 0 0 . . . nI - 1 n n1 n2 ... nI (n1, n2, . . . , nI) = = n1 - n2 1 n -n1n2 n . . . -n1nI n -n2n1 n n2 - n2 2 n -n2nI n ... ... -nI n1 n nI - n2 I n . Jedna g-inverzia (Z1D1Z1)- = 1 n1 0 . . . 0 0 1 n2 ... ... 0 0 . . . 1 nI (dokážte ako cvičenie). 27 Dostávame (41) D1 - D12 = (I - 1 n 11 )Z1 1 n1 0 . . . 0 0 1 n2 ... ... 0 0 . . . 1 nI Z1(I - 1 n 11 ) = = Z1 1 n1 0 . . . 0 0 1 n2 ... ... 0 0 . . . 1 nI Z1 - 1 n 1(n1, n2, . . . , nI) 1 n1 0 . . . 0 0 1 n2 ... ... 0 0 . . . 1 nI Z1- - 1 n Z1 1 n1 0 . . . 0 0 1 n2 ... ... 0 0 . . . 1 nI n1 n2 ... nI 1 + 1 n 11 = = Z1 1 n1 0 . . . 0 0 1 n2 ... ... 0 0 . . . 1 nI Z1 - 1 n 1(1, 1, . . . , 1)Z1 - 1 n Z1 1 1 ... 1 1 + 1 n 11 . Platí tr(D1 - D12)=tr Z1Z1 1 n1 0 . . . 0 0 1 n2 ... ... 0 0 . . . 1 nI - 2 1 n tr 1 Z1 1 1 ... 1 + 1 n tr 1 1 = I - 1. Počítajme pomocou (41) (42) S2 = Y (D1 - D12)Y = Y Z1 1 n1 0 . . . 0 0 1 n2 ... ... 0 0 . . . 1 nI Z1Y - - 2 n I i=1 ni j=1 Yij (1, 1, . . . , 1) n1 j=1 Y1j n2 j=1 Y2j ... nI j=1 YIj + 1 n I i=1 ni j=1 Yij 2 = = n1 j=1 Y1j, n2 j=1 Y2j, . . . , nI j=1 YIj 1 n1 0 . . . 0 0 1 n2 ... ... 0 0 . . . 1 nI n1 j=1 Y1j n2 j=1 Y2j ... nI j=1 YIj - 28 - 1 n I i=1 ni j=1 Yij 2 = I i=1 1 ni ni j=1 Yij 2 - 1 n I i=1 ni j=1 Yij 2 . Pretože (43) ^ ^ = Y X(X X)X Y = 1 n I i=1 ni j=1 Yij 2 , použitím (31) dostávame (44) Se = I i=1 ni j=1 Y 2 ij - S2 - ^ ^ = I i=1 ni j=1 Y 2 ij - I i=1 1 ni ni j=1 Yij 2 . Pre riešenie rovníc (36) potrebujeme ešte spočítať tr Z1D1Z1 = tr Z1(I - X(X X)X )Z1 = tr Z1Z1 - 1 n tr Z1XX Z1 = = I i=1 ni - 1 n I i=1 n2 i = n - 1 n I i=1 n2 i . Rovnica (36) má tvar E S2 Se = n - 1 n I i=1 n2 i I - 1 0 n - I 2 1 2 0 . Konečne z rovnice (39) dostávame hľadané nevychýlené odhady variančných kom- ponentov (45) 2 1 = 1 n - 1 n I i=1 n2 i S2 I - 1 n - 1 n I i=1 n2 i (n - I) Se, (46) 2 0 = 1 n - I Se, kde S2 = I i=1 1 ni ni j=1 Yij 2 - 1 n I i=1 ni j=1 Yij 2 a Se = I i=1 ni j=1 Y 2 ij - I i=1 1 ni ni j=1 Yij 2 . 29 Pre vyvážený model t.j. v prípade, že n1 = n2 = . . . = nI = t dostávame 2 1 = 1 t S2 I - 1 - Se n - I = = It - 1 I(I - 1)t2(t - 1) I i=1 t j=1 Yij 2 - 1 It(t - 1) I i=1 t j=1 Y 2 ij(47) - 1 I(I - 1)t2 I i=1 t j=1 Yij 2 , 2 0 = Se n - I = (48) = 1 I(t - 1) I i=1 t j=1 Y 2 ij - 1 It(t - 1) I i=1 t j=1 Yij 2 . Lema 19. V prípade normality rozdelenia vektora Y má Se 2 0 vo vyváženom modeli (40) 2 (n - I, 0) rozdelenie. Dôkaz. Y N(X; 2 1Z1Z1 + 2 0I), pričom cov(Y ) je p.d. matica. Podľa (44) je Se = Y Y - 1 t Y Z1Z1Y = Y (I - Z1(Z1Z1)-1 Z1)Y , teda Se 2 0 = Y I - Z1(Z1Z1)-1 Z1 2 0 Y . Podľa vety 16 (i) I - Z1(Z1Z1)-1 Z1 2 0 (2 1Z1Z1 + 2 0I) I - Z1(Z1Z1)-1 Z1 2 0 = = 1 4 0 {[2 1Z1Z1 - 2 1Z1Z1 + 2 0I - 2 0Z1(Z1Z1)-1 Z1][I - Z1(Z1Z1)-1 Z1]} = = 1 4 0 {2 0I - 2 0Z1(Z1Z1)-1 Z1} = I - Z1(Z1Z1)-1 Z1 2 0 , (ii) (2 1Z1Z1 + 2 0I) I - Z1(Z1Z1)-1 Z1 2 0 X = = (I - Z1(Z1Z1)-1 Z1)X (I - Z1(Z1Z1)-1 Z1), 30 a podľa (ii) (iii) X I - Z1(Z1Z1)-1 Z1 2 0 (2 1Z1Z1 + 2 0I) I - Z1(Z1Z1)-1 Z1 2 0 X = = X I - Z1(Z1Z1)-1 Z1 2 0 X. Preto má Se 2 0 rozdelenie 2 s k = tr I - Z1(Z1Z1)-1 Z1 2 0 (2 1Z1Z1 + 2 0I) = n - I stupňami voľnosti a s parametrom necentrality = X I - Z1(Z1Z1)-1 Z1 2 0 (2 1Z1Z1 + 2 0I) I - Z1(Z1Z1)-1 Z1 2 0 × × (2 1Z1Z1 + 2 0I) I - Z1(Z1Z1)-1 Z1 2 0 X = 2 2 0 (n - n) = 0 (centrálne 2 (n - I, 0) rozdelenie). Lema 20. V prípade vyváženého modelu (40) (t.j. n1 = n2 = . . . = nI = t) má S2 t2 1 + 2 0 rozdelenie 2 (I - 1, 0). Dôkaz. Podľa (42) je S2 = Y Z1(Z1Z1)Z1Y - 1 n Y XX Y = Y 1 t Z1Z1 - 1 n XX Y . Postupujeme úplne analogicky ako v leme 19 a dokážeme (i), (ii) a (iii) podľa vety 16 pre maticu kvadratickej formy 1 t Z1Z1 - 1 n XX t2 1 + 2 0 . Vskutku, S2 t2 1 + 2 0 má 2 (I - 1, 0) rozdelenie, lebo pre stupne voľnosti k platí k = tr 1 t Z1Z1 - 1 n XX t2 1 + 2 0 (2 1Z1Z1 + 2 0I) = I - 1 a koeficient necentrality je = X 1 t Z1Z1 - 1 n XX t2 1 + 2 0 (2 1Z1Z1 + 2 0I) 1 t Z1Z1 - 1 n XX t2 1 + 2 0 × × (2 1Z1Z1 + 2 0I) 1 t Z1Z1 - 1 n XX t2 1 + 2 0 X = = 2 t2 1 + 2 0 X 1 t Z1Z1 - 1 n XX X = = 2 t2 1 + 2 0 1 t 1 Z1Z11 - 1 n 1 11 1 = 0. 31 Lema 21. V prípade vyváženého modelu (40) a normálne rozdeleného Y sú Se 2 0 a S2 t2 1 + 2 0 nezávislé. Dôkaz. Pretože Y má nesingulárnu kovariančnú maticu a matice kvadratických foriem Se 2 0 a S2 t2 1 + 2 0 sú symetrické, stačí podľa vety 17 ukázať, že 1 2 0(t2 1 + 2 0) (I - 1 t Z1Z1)(2 1Z1Z1 + 2 0I)( 1 t Z1Z1 - 1 n XX ) = 0. Toto dostaneme po krátkom výpočte. V prípade vyváženého modelu môžeme štatistikou (49) (n - I) (I - 1) S2 Se testovať hypotézu H0 : 2 1 = 0. Za platnosti H0 má (49) FI-1,n-I rozdelenie (podľa lemy 19, lemy 20 a lemy 21). Cvičenie 1. Ukážte, že vo vyváženom modeli s jedným náhodným efektom v prípade normality rozdelenia b a sa disperzie odhadov (47) a (48) rovnajú D(2 1) = 2(2 1t + 2 0)2 t2(I - 1) + 24 0 It2(t - 1) , D(2 0) = 24 0 I(t - 1) . 7 ODHADY VARIANČNÝCH KOMPONENTOV VO VYVÁŽENOM MODELI S DVOMI NÁHODNÝMI EFEKTMI A INTERAKCIOU Vyvážený model s dvomi náhodnými efektmi a interakciou sa dá písať ako (50) Yijk = + b1i + b2j + b3ij + ijk, i = 1, 2, . . . , I , j = 1, 2, . . . , J , k = 1, 2, . . . , K . 32 Náhodný observačný vektor Y N(1n, cov(Y )), 1 = (1, 1, . . . , 1) Rn , n = IJK; náhodné vektory b1 N(0I,1, 2 1I), b2 N(0J,1, 2 2I), b3 N(0IJ,1, 2 3I) aj N(0n,1, 2 0I) sú navzájom nezávislé, je skalár (číslo, neznáma konštanta). Model (50) sa dá písať aj v tvare (51) Yn,1 = (X ...Z1 ...Z2 ...Z3) b1 b2 b3 + = = X + Z1b1 + Z2b2 + Z3b3 + , kde X = 1n, Z1 je rozmerov n×I, Z2 je rozmerov n×J a Z3 je rozmerov n×IJ. E(Y ) = X = 1, cov(Y ) = 2 1Z1Z1 + 2 2Z2Z2 + 2 3Z3Z3 + 2 0I. Platí Z1 = II,I 1JK = II,I (1J 1K) ( znamená Kroneckerov súčin matíc, t.j. ak prvky matice A sú {A }ij , prvky matice B sú {B}ij , tak prvky matice A B sú {A }ijB), Z2 = IJ,J IJ,J ... IJ,J 1K = (1I IJ,J ) 1K a Z3 = IIJ,IJ 1K = II,I IJ,J 1K. Ľahko sa ukáže, že X X = IJK, X Z1 = JK(1, 1, . . . , 1)1,I = JK1I, X Z2 = IK(1, 1, . . . , 1)1,J = IK1J , X Z3 = K(1, 1, . . . , 1)1,IJ = K1IJ , Z1Z1 = JKII,I, Z1Z2 = K(1I 1J ), Z1Z3 = K (II,I 1J ), Z2Z2 = IKIJ,J , Z2Z3 = K (1I IJ,J ), Z3Z3 = K (II,I IJ,J ). V ďalšom budeme potrebovať ešte nasledujúce tvrdenie: 33 Lema 22. Ak m je prirodzené číslo, tak jedna g-inverzia matice 1 - 1 m - 1 m - 1 m . . . - 1 m - 1 m 1 - 1 m - 1 m . . . - 1 m - 1 m - 1 m 1 - 1 m . . . - 1 m ... ... ... ... - 1 m - 1 m - 1 m . . . 1 - 1 m je matica Im,m. Dôkaz je jednoduchý a urobte ho ako cvičenie. Pri odhade variančných komponentov budeme postupovať podľa návodu zo str. 17. Z postupnosti modelov (formálne predpokladáme pevné parametre bi, i = 1, 2, 3) M : E(Y ) = (X ...Z1 ...Z2 ...Z3) b1 b2 b3 , cov(Y ) = 2 0I; M1 : E(Y ) = (X ...Z1 ...Z2) b1 b2 , cov(Y ) = 2 0I; M2 : E(Y ) = (X ...Z1) b1 , cov(Y ) = 2 0I; M3 : E(Y ) = X, cov(Y ) = 2 0I. dostaneme štatistiky ^ = (X ...Z) X X X Z Z X Z Z - X Z Y , ^ = (X ...Z1 ...Z2) X X X Z1 X Z2 Z1X Z1Z1 Z1Z2 Z2X Z2Z1 Z2Z2 - X Z1 Z2 Y , ^ = (X ...Z1) X X X Z1 Z1X Z1Z1 - X Z1 Y a ^ = X(X X)X Y . Platí (podľa (31)) Se = Y Y - S1 - S2 - S3 - S4, pričom podľa (27) S1 = Y (I - D1)Y , kde D1 = I - X(X X)X , podľa (28) S2 = Y (D1 - D12)Y , kde D12 = D1 - D1Z1(Z1D1Z1)- Z1D1, 34 podľa (29) S3 = Y (D12 - D123)Y , kde D123 = D12 - D12Z2(Z2D12Z2)- Z2D12 a podľa (30) S4 =Y (D123 -D1234)Y =S4, kde D1234 =D123-D123Z3(Z3D123Z3)- Z3D123. Aby sme určili odhady 2 1, 2 2, 2 3 a 2 0 , potrebujeme spočítať matice D1, D12, D123, D1234 a prvky matice E z (36) (alebo matice E-1 z (39)). Lema 23. V modeli (50) majú matice D1 a D12 (pozri (27) a (28)) tvar D1 = IIJK,IJK - 1 IJK 1IJK1IJK, D12 = IIJK,IJK - 1 JK Z1Z1. Kvadratické formy S1 a S2 sú S1 = Y 1 IJK 11 Y = 1 IJK I i=1 J j=1 K k=1 Yijk 2 , S2 = Y 1 JK Z1Z1 - 1 IJK 11 Y = = 1 JK I i=1 J j=1 K k=1 Yijk 2 - 1 IJK I i=1 J j=1 K k=1 Yijk 2 . Dôkaz urobte ako cvičenie. Lema 24. V modeli (50) majú matice D123 a D123 - D1234 (pozri (29) a (30)) tvar D123 = IIJK,IJK - 1 JK Z1Z1 - 1 IK Z2Z2 + 1 IJK 11 . D123 - D1234 = 1 K Z3Z3 - 1 JK Z1Z1 - 1 IK Z2Z2 + 1 IJK 11 . Dôkaz urobte ako cvičenie. Lema 25. V modeli (50) kvadratické formy S3 a S4 sú S3 = Y 1 IK Z2Z2 - 1 IJK 11 Y = = 1 IK J j=1 I i=1 K k=1 Yijk 2 - 1 IJK I i=1 J j=1 K k=1 Yijk 2 , 35 S4 = Y 1 K Z3Z3 - 1 JK Z1Z1 - 1 IK Z2Z2 + 1 IJK 11 Y = = 1 K I i=1 J j=1 K k=1 Yijk 2 - 1 JK I i=1 J j=1 K k=1 Yijk 2 - - 1 IK J j=1 I i=1 K k=1 Yijk 2 + 1 IJK I i=1 J j=1 K k=1 Yijk 2 . Dôkaz urobte ako cvičenie. Lema 26. V modeli (50) je kvadratická forma Se = Y I - 1 K Z3Z3 Y = I i=1 J j=1 K k=1 Y 2 ijk - 1 K I i=1 J j=1 K k=1 Yijk 2 , Se 2 0 má 2 (IJ(K - 1), 0) rozdelenie. Dôkaz. Využite (31) a urobte ho ako cvičenie. Lema 27. V modeli (50) platí tr Z1D1Z1 = JK(I - 1), tr Z2(D1 - D12)Z2 = 0, tr Z3(D1 - D12)Z3 = K(I - 1), tr(D1 - D12) = I - 1, tr Z2D12Z2 = IK(J - 1), tr Z3(D12 - D123)Z3 = K(J - 1), tr(D12 - D123) = J - 1, tr Z3D123Z3 = K(I - 1)(J - 1), tr(D123 - D1234) = (I - 1)(J - 1), tr D1234 = IJ(K - 1). Dôkaz. (Pri dôkaze posledného vzťahu využite (31) a lemu 12.) Urobte ho ako cvičenie. Lema 28. V modeli (50) sú kvadratické nevychýlené odhady 2 1, 2 2, 2 3 a 2 0 z rovnice (36) 2 0 = Se IJ(K - 1) , 2 3 = 1 K S4 (I - 1)(J - 1) - Se IJ(K - 1) , 2 2 = 1 IK S3 J - 1 - S4 (I - 1)(J - 1) , 2 1 = 1 JK S2 I - 1 - S4 (I - 1)(J - 1) . 36 Dôkaz urobte ako cvičenie. Lema 29. V modeli (50) má S4 2 0 + K2 3 rozdelenie 2 ((I - 1)(J - 1), 0). Dôkaz urobte ako cvičenie. Lema 30. V modeli (50) má S3 2 0 + K2 3 + IK2 2 rozdelenie 2 (J - 1, 0). Dôkaz urobte ako cvičenie. Lema 31. V modeli (50) má S2 2 0 + K2 3 + JK2 1 rozdelenie 2 (I - 1, 0). Dôkaz urobte ako cvičenie. Lema 32. V modeli (50) sú kvadratické formy Se 2 0 a S4 2 0 + K2 3 nezávislé. Dôkaz urobte ako cvičenie. Lema 33. V modeli (50) sú kvadratické formy S2 2 0 + K2 3 + JK2 1 a S4 2 0 + K2 3 nezávislé. Dôkaz urobte ako cvičenie. Lema 34. V modeli (50) sú kvadratické formy S3 2 0 + K2 3 + IK2 2 a S4 2 0 + K2 3 nezávislé. Dôkaz urobte ako cvičenie. Testujme v modeli (50) hypotézu H0 : 2 3 = 0. Za platnosti H0 má (52) S4 (I-1)(J-1) Se IJ(K-1) = IJ(K - 1) (I - 1)(J - 1) S4 Se rozdelenie F(I-1)(J-1),IJ(K-1) (podľa lemy 26, lemy 29 a lemy 32; I, J, K > 1). Testujme v modeli (50) hypotézu H0 : 2 1 = 0. Za platnosti H0 má (53) S2 I-1 S4 (I-1)(J-1) = (J - 1) S2 S4 37 rozdelenie FI-1,(I-1)(J-1) (podľa lemy 29, lemy 31 a lemy 33; I, J > 1). Testujme v modeli (50) hypotézu H0 : 2 2 = 0. Za platnosti H0 má (54) S3 J-1 S4 (I-1)(J-1) = (I - 1) S3 S4 rozdelenie FJ-1,(I-1)(J-1) (podľa lemy 29, lemy 30 a lemy 34; I, J > 1). 8 MINQUE ODHADY VARIANČNÝCH KOMPONENTOV Uvažujme zmiešaný lineárny model (55) Y = X + Z1b1 + Z2b2 + . . . + Zsbs + Zs+1bs+1. (Ak Zs+1 označíme Z0, bs+1 označíme b0, 2 s+1 označíme 2 0 , dostávame model (15).) Kovariančná matica cov(Y ) = 2 1Z1Z1 + 2 2Z2Z2 + . . . + 2 s ZsZs + 2 s+1Zs+1Zs+1 = s+1 i=1 2 i Vi. Náhodné vektory b1, b2, . . . , bs+1 sa nazývajú hypotetické (alebo tiež nepozorovateľné) premenné. Uvažujme kvadratický odhad Y AY funkcie (56) p = (p1, p2, . . . , ps+1) 2 1 ... 2 s 2 s+1 , kde p je daný (s+1)-rozmerný vektor. Požadujeme symetriu matice A kvadratickej formy Y AY . Ďalej chceme, aby (i) odhad Y AY bol invariantný, t.j. (57) Rq Y AY = (Y - X) A(Y - X); 38 (ii) odhad Y AY bol nevychýleným odhadom (56), t.j. (58) {2 i (0, ) i = 1, 2, . . . , s + 1} E(Y AY ) = s+1 i=1 pi2 i ; (iii) Z AZ - 2 = tr(Z AZ - )(Z AZ - ) = min, kde Z = (Z1 ...Z2 ... . . . ...Zs+1) a = p1 m1 Im1,m1 0 . . . 0 0 p2 m2 Im2,m2 . . . 0 ... 0 0 . . . ps+1 ms+1 Ims+1,ms+1 . . Odhad splňujúci (i), (ii), (iii), pričom A = A sa nazýva MINQUE (minimum norm quadratic unbiased estimator). Poznámka 35. Podľa [2], str. 71 platí (59) P{Y (X ... cov(Y ))} = 1. Lema 36. Kvadratický odhad Y AY (A = A ) je invariantný s pravdepodobnosťou 1 práve vtedy ak AX = 0. Dôkaz. Odhad Y AY je invariantný Rq Y AY = (Y -X) A(Y -X) = Y AY -2 X AY + X A X (60) X AY = 0 & X AX = 0. Pretože podľa (59) P{Y = Xu+cov(Y )w} = 1, platí (60) s pravdepodobnosťou 1 práve ak X AX = 0 & X A cov(Y ) = 0. Je zrejmé, že ak AX = 0, je odhad Y AY (A = A ) invariantný (s pravdepodobnosťou 1). Naopak, ak je odhad Y AY (A = A ) invariantný, teda platí (60), tak s pravdepodobnosťou 1 je Y AY = Y [I - X(X X)X + X(X X)X ] A [I - X(X X)X + X(X X)X ]Y = (61) = Y [I - X(X X)X ] A[I - X(X X)X ]Y , lebo s pravdepodobnosťou 1 je Y [I - X(X X)X ] AX(X X)X Y = 39 = Y [X(X X)X ] A[I - X(X X)X ]Y = = Y [X(X X)X ] A[X(X X)X ]Y = 0. Zo vzťahu (61) vidíme, že s pravdepodobnosťou 1 sa Y AY rovná Y [I - X(X X)X ] A[I - X(X X)X ]Y , pričom [I - X(X X)X ] A[I - X(X X)X ] je symetrická a [I - X(X X)X ] A[I - X(X X)X ]X = 0. Lema 37. Kvadratický invariantný odhad Y AY (A = A ) je nevychýleným odhadom (56) práve vtedy ak (62) tr A ZiZi = pi, i = 1, 2, . . . , s + 1. Dôkaz. Invariantný odhad Y AY je nevychýleným odhadom (56) práve vtedy ak {2 i (0, ) i = 1, 2, . . . , s + 1} E(Y AY ) = X AX + tr A cov(Y ) = = X AX + s+1 i=1 2 i tr A ZiZi = s+1 i=1 pi2 i , čo je splnené práve vtedy ak tr A ZiZi = pi, i = 1, 2, . . . , s + 1. Poznámka 38. Keby sme mohli (izolovane) realizovať hypotetickú premennú bi i {1, 2, . . . , s + 1}, ktorej stredná hodnota je 0mi,1 a cov(bi) = 2 i Imi,mi , tak ,,prirodzený odhad 2 i by bol 2 i = bibi mi . (Podľa lemy 10 je to nevychýlený odhad 2 i .) ,,Prirodzený odhad (56) je p = (b1 ...b2 ... . . . ...bs+1) b1 b2 ... bs+1 . Rozdiel medzi ,,prirodzeným odhadom a invariantným nevychýleným odhadom Y AY (A = A ) funkcie (56) je Y AY - (b1 ...b2 ... . . . ...bs+1) b1 b2 ... bs+1 = 40 = (Y - X) A (Y - X) - (b1 ...b2 ... . . . ...bs+1) b1 b2 ... bs+1 = = (b1 ...b2 ... . . . ...bs+1)[Z AZ - ] b1 b2 ... bs+1 . Jedna z možností minimalizovať tento rozdiel je (iii) (minimalizácia euklidovskej normy). Veta 39. Nech v modeli (55) je matica V = ZZ = Z1Z1 +Z2Z2 +. . .+Zs+1Zs+1 nesingulárna. Kvadratická forma Y AY , kde (63) A = s+1 i=1 i[V-1 - V-1 X(X V-1 X)X V-1 ]ZiZi× × [V-1 - V-1 X(X V-1 X)X V-1 ] je MINQUE odhadom (56) ak čísla 1, 2, . . . , s+1 riešia systém tr A ZiZi = pi, i = 1, 2, . . . , s + 1. Dôkaz. A = A . Podľa lemy 36 je Y AY invariantný, lebo AX = 0. Podľa lemy 37 je odhad Y AY nevychýlený. Pretože Z AZ - 2 = tr(Z AZ - )(Z AZ - ) = = tr Z AZZ AZ - 2 tr Z AZ + tr = = tr AVAV - 2 tr AZZ + tr = = tr AVAV - 2 tr A s+1 i=1 pi mi ZiZi + tr = = tr AVAV - 2 s+1 i=1 p2 i mi + tr = = tr AVAV - tr , v (iii) musíme vlastne minimalizovať tr AVAV, pričom A = A , tr AZiZi = pi, i = 1, 2, . . . , s + 1 a AX = 0. Keby Y AY bol tiež MINQUE odhadom (56), tak A sa vždy dá napísať ako A = A + D, kde D = D , DX = 0 a tr DZiZi = 0, i = 1, 2, . . . , s + 1. Potom ale tr AVAV = tr(A + D)V(A + D)V = tr AVAV + tr DVDV tr AVAV, 41 lebo tr AVDV = (tr DVAV =) = tr s+1 i=1 i[V-1 - V-1 X(X V-1 X)X V-1 ]ZiZi× × [V-1 - V-1 X(X V-1 X)X V-1 ]VDV = = s+1 i=1 i tr[I - X(X V-1 X)X V-1 ]ZiZi[I - V-1 X(X V-1 X)X ]D = = s+1 i=1 i tr DZiZi = 0 a tr DVDV = tr DZZ DZZ = tr Z DZZ DZ = tr(Z DZ)(Z DZ) 0. Dôsledok 40. Nech (s + 1)-rozmerný vektor g je g = Y BZ1Z1BY Y BZ2Z2BY ... Y BZs+1Zs+1BY = Y BV1BY Y BV2BY ... Y BVs+1BY a (s + 1) × (s + 1) matica S je S = tr BV1BV1 tr BV1BV2 . . . tr BV1BVs+1 tr BV2BV1 tr BV2BV2 . . . tr BV2BVs+1 ... ... tr BVs+1BV1 tr BVs+1BV2 . . . tr BVs+1BVs+1 , kde B = V-1 - V-1 X(X V-1 X)X V-1 . MINQUE odhad funkcie p = s+1 i=1 pi2 i je p S- g. Dôkaz. Podľa vety 39 je Y AY MINQUE odhadom funkcie p , keď (65) Y AY = Y s+1 i=1 iBViBY = g, pričom = (1, 2, . . . , s+1) spĺňa rovnice tr AZiZi = pi, i = 1, 2, . . . , s + 1 (A = s+1 i=1 iBZiZiB je dané vzorcom (63)). Toto zapísané formou matíc znie (66) S = p, čiže p (S). Potom ale = Sp, kde Sje symetrická g-inverzia matice S, je jedno riešenie (66) (pozri napr. [2], str. 70) a podľa (65) je p Sg MINQUE odhadom funkcie p . 42 Dôsledok 41. Nech = (2 1, 2 2, . . . , 2 s+1). MINQUE odhad funkcie p je p ^, kde ^ je riešením sústavy (67) S ^ = g. Dôkaz. Podľa dôsledku 40 je p Sg MINQUE odhadom funkcie p , kde p (S). Podľa [2], str. 70 {^ = Sg + (SS - I)z : z Rs+1 , Sje jedna (ľubovoľná) g-inverzia matice S} sú všetky riešenia sústavy (67). Preto p ^ = p (Sg + (SS - I)z) = p S- g je MINQUE odhadom funkcie p . Poznámka 42. Podľa vety 18 a lemy 10 sa disperzia MINQUE odhadu Y AY v prípade Y N(X; V) rovná D(Y AY ) = E(Y AY - E(Y AY ))2 = E(Y AY Y AY ) - E2 (Y AY ) = = 2 tr AVAV + tr2 AV + 2X AX tr AV + 4X AVA X+ + (X A X)2 - (X A X + tr AV)2 = = 2 tr AVAV. Preto MINQUE odhad je v tomto prípade V-LMVUIE (V-lokálny nevychýlený invariantný odhad s minimálnou disperziou) odhad funkcie p . Cvičenie 1. Pomocou lemy 5 dokážte, že pre ľubovoľné matice A, B vhodných rozmerov je (AA + BB ) = (A ...B). 2. Ak vo vete 39 je Xs+1 = I, tak V je nesingulárna. Dokážte. 9 MINQUE ODHADY VARIANČNÝCH KOMPONENTOV V MODELI S JEDNÝM NÁHODNÝM EFEKTOM Hľadajme MINQUE odhady variančných komponentov v modeli s jedným náhodným efektom, teda v modeli (40), pričom pri označení podľa (55) sa s = 1, X = = 1n,1 Z1 = 1n1,1 e1 1n2,1 e2 ... 1nI ,1 eI (w-ta súradnica I-rozmerného vektora ew sa rovná 1, w {1, 2, . . . , I}, ostatné súradnice sa rovnajú 0), Z2 = In,n, n = I i=1 ni (2 2 je v modeli (40) označené ako 2 0 ). 43 Lema 43. V modeli (40) (pri označení matice B podľa dôsledku 40) platí V-1 = (Z1Z1 + I)-1 = I - Z1 1 n1+1 0 . . . 0 0 1 n2+1 ... 0 ... ... ... 0 0 0 . . . 1 nI +1 Z1, (X V-1 X)- = I i=1 ni ni + 1 -1 , V-1 X(X V-1 X)X V-1 = = I i=1 ni ni + 1 -1 1n,111,n - Z1 n1 n1+1 n2 n2+1 ... nI nI +1 11,n- -1n,1 n1 n1 + 1 , . . . , nI nI + 1 Z1 + Z1 n1 n1+1 n2 n2+1 ... nI nI +1 n1 n1 + 1 , . . . , nI nI + 1 Z1 , tr BV1BV1 = tr(V-1 - V-1 X(X V-1 X)X V-1 )Z1× × Z1(V-1 - V-1 X(X V-1 X)X V-1 )Z1Z1 = = I i=1 ni ni + 1 2 - 2 I i=1 ni ni + 1 -1 I i=1 ni ni + 1 3 + + I i=1 ni ni + 1 -2 I i=1 ni ni + 1 2 2 , tr BV1BV2 = tr BZ1Z1BI = tr Z1BBZ1 = = I i=1 ni (ni + 1)2 - 2 I i=1 ni ni + 1 -1 I i=1 ni ni + 1 2 + +2 I i=1 ni ni + 1 -1 I i=1 ni ni + 1 3 + n I i=1 ni ni + 1 -2 I i=1 ni ni + 1 2 - -2 I i=1 ni ni + 1 -2 I i=1 ni ni + 1 2 I j=1 n2 j nj + 1 + 44 + I i=1 ni ni + 1 -2 I i=1 ni ni + 1 2 I j=1 n3 j (nj + 1)2 , tr BV2BV2 = tr BB = = n - 2 I i=1 ni ni + 1 - 2n I i=1 ni ni + 1 -1 + +6 I i=1 ni ni + 1 -1 I i=1 n2 i ni + 1 - 6 I i=1 ni ni + 1 -1 I i=1 n3 i (ni + 1)2 + + I i=1 ni ni + 1 2 + 2 I i=1 ni ni + 1 -1 I i=1 n4 i (ni + 1)3 + +n2 I i=1 ni ni + 1 -2 - 4 I i=1 ni ni + 1 -2 n I i=1 n2 i ni + 1 + +4 I i=1 ni ni + 1 -2 I i=1 n2 i ni + 1 2 + 2 I i=1 ni ni + 1 -2 n I i=1 n3 i (ni + 1)2 - -4 I i=1 ni ni + 1 -2 I i=1 n2 i ni + 1 I j=1 n3 j (nj + 1)2 + I i=1 ni ni + 1 -2 I i=1 n3 i (ni + 1)2 2 , Y BV1BY = Y BZ1Z1BY = = I k=1 nk nk + 1 nk j=1 Ykj - nk I i=1 ni ni+1 I s=1 1 ns + 1 ns v=1 Ysv 2 , Y BV2BY = Y BBY = = I i=1 ni j=1 Y 2 ij - I i=1 1 ni + 1 + 1 (ni + 1)2 ni j=1 Yij 2 - -2 I i=1 ni ni + 1 -1 I i=1 1 ni + 1 ni j=1 Yij I i=1 1 (ni + 1)2 ni j=1 Yij + + I i=1 ni ni + 1 -2 I i=1 ni (ni + 1)2 -1 I i=1 1 ni + 1 ni j=1 Yij . Dôkaz. Vykonáme jednoduchým ale zdĺhavým postupným výpočtom. Urobte ho ako cvičenie. 45 Dôsledok 44. Vo vyváženom modeli s jedným náhodným efektom (označenie podľa (55), n1 = n2 = . . . = nI = t) platí: tr BV1BV1 = (I - 1) t t + 1 2 , tr BV1BV2 = (I - 1) t (t + 1)2 , tr BV2BV2 = It3 + It2 - It - 1 (t + 1)2 , Y BV1BY = 1 (t + 1)2 I i=1 t j=1 Yij 2 - 1 I(t + 1)2 I i=1 t j=1 Yij 2 , Y BV2BY = I i=1 t j=1 Y 2 ij t + 2 (t + 1)2 I i=1 t j=1 Yij 2 - 1 It(t + 1)2 I i=1 t j=1 Yij 2 . Dôkaz. Vykonáme dosadením n1 = n2 = . . . = nI = t do vzorcov z lemy 43. Urobte ho ako cvičenie. Lema 45. Vo vyváženom modeli s jedným náhodným efektom platí S= S-1 = 1 It(t - 1) It3 + It2 - It - 1 (I - 1)t -1 -1 t (označenie z dôsledku 40). Dôkaz. Vykonáme inverziou matice S, ktorej prvky sú dané v dôsledku 44. Urobte ho ako cvičenie. Veta 46. Vo vyváženom modeli s jedným náhodným efektom sú MINQUE odhady 2 1 a 2 2 totožné s odhadmi získanými metódou Hendersona (vzorce (47) a (48), pričom 2 2 je pri našom označení totožné so 2 0 vo vzorci (48)). Dôkaz. Podľa dôsledku 40 sa MINQUE odhad 2 1 rovná 2 1 = (1 , 0)S- g, pričom matica Sje daná v leme 45 a prvky vektora g (definovaného v dôsledku 40) sú dané v dôsledku 44. Preto MINQUE odhad 2 1 = 1 It(t - 1) It3 + It2 - It - 1 (I - 1)t , -1 g = 46 = 1 It(t - 1) It3 + It2 - It - 1 (I - 1)t 1 (t + 1)2 I i=1 t j=1 Yij 2 - - 1 I(t + 1)2 I i=1 t j=1 Yij 2 - - I i=1 t j=1 Y 2 ij + t + 2 (t + 1)2 I i=1 t j=1 Yij 2 + 1 It(t + 1)2 I i=1 t j=1 Yij 2 = = It - 1 I(I - 1)t2(t - 1) I i=1 t j=1 Yij 2 - 1 It(t - 1) I i=1 t j=1 Y 2 ij- - 1 I(I - 1)t2 I i=1 t j=1 Yij 2 , čo je ten istý odhad ako (47). Podobne MINQUE odhad 2 2 = (0 , 1)Sg = = 1 It(t - 1) (-1 , t)× × 1 I(t+1)2 I I i=1 t j=1 Yij 2 - I i=1 t j=1 Yij 2 I i=1 t j=1 Y 2 ij - t+2 (t+1)2 I i=1 t j=1 Yij 2 - 1 It(t+1)2 I i=1 t j=1 Yij 2 = = 1 It(t - 1) -t2 - 2t - 1 (t + 1)2 I i=1 t j=1 Yij 2 + t I i=1 t j=1 Y 2 ij = = 1 I(t - 1) I i=1 t j=1 Y 2 ij - 1 It(t - 1) I i=1 t j=1 Yij 2 , čo je ten istý odhad ako (48). 47 10 ODHADY VARIANČNÝCH KOMPONENTOV ZÍSKANÉ METÓDOU MAXIMÁLNEJ VIEROHODNOSTI (ML ­ ODHADY) Jednou zo štandardných metód získania odhadov je metóda maximálnej vierohodnosti. Vyžaduje si predpoklad o rozdelení pravdepodobnosti pre náhodné premenné, ktorých realizáciami sú získané (namerané) údaje. Prirodzený predpoklad je normalita rozdelenia (v neposlednom rade pre matematické zvládnutie riešenia). Majme zmiešaný lineárny model (55), teda Yn,1 = (X ...Z1 ... . . . Zs+1) b1 ... bs+1 = X + s+1 i=1 Zibi, cov(Y ) = s+1 i=1 2 i Vi = s+1 i=1 2 i ZiZi = V , bi Rmi , i = 1, 2, . . . , s + 1, E(bi) = 0, var(bi) = 2 i I, i = 1, 2, . . . , s + 1, cov(bi, bj) = 0 pre i = j (opäť pri označení Zs+1 = Z0, bs+1 = b0, 2 s+1 = 2 0 dostávame model (15)). Predpokladáme Y Nn(X, V ), teda hustota rozdelenia pravdepodobnosti observačného vektora je f(y) = 1 (2) n 2 det V e- 1 2 (y-X) V -1 (y-X) . Funkcia vierohodnosti (,,z jednej realizácie y) je (68) L(, 2 1, . . . , 2 s+1) = 1 (2) n 2 det V e- 1 2 (y-X) V -1 (y-X) . Realizácie odhadov ^, 2 1, . . . , 2 s+1 získané metódou maximálnej vierohodnosti (ML­odhady) sú také hodnoty ^real, {2 1}real, . . . , {2 s+1}real, pre ktoré (68) nadobúda maximum (pri realizácii y). Maximum (68) sa dosiahne aj maximalizáciou (69) l = ln L = - 1 2 n ln 2 - 1 2 ln det V - 1 2 (y - X) V -1 (y - X). Pri vypočítaní takéhoto maxima využijeme pomocné tvrdenie: Lema 47. Nech a, x sú n × 1 vektory, a x x je n × 1 vektor so súradnicami a x x i = a x xi , i = 1, 2, . . . , n. 48 Nech A je symetrická matica n × n, x Ax x je n × 1 vektor so súradnicami x Ax x i = x Ax xi , i = 1, 2, . . . , n. Nech B je symetrická nesingulárna n × n matica, ktorej prvky bij sú funkciami premennej t, čiže bij = bij(t), i, j = 1, 2, . . . , n, B t je n × n matica, ktorej prvky sú bij(t) t , i, j = 1, 2, . . . , n, det B B je n × n matica, ktorej prvky sú det B bij , i, j = 1, 2, . . . , n. Platí (i) a x x = x a x = a, (ii) x Ax x = 2Ax, (iii) B-1 t = -B-1 B t B-1 , (iv) ln det B t = tr B-1 B t . Dôkaz. Tvrdenia (i) a (ii) dokážeme ľahko rozpísaním skalárnych súčinov (pozri tiež [9], str. 98). Pri tvrdení (iii) si pomôžeme nasledujúcim spôsobom: Prvky matice B-1 označme b (-1) ij , i, j = 1, 2, . . . , n. Tiež sú funkciami premennej t, čiže b (-1) ij = b (-1) ij (t), i, j = 1, 2, . . . , n. Pre i, j = 1, 2, . . . , n je {BB-1 }i,j = n k=1 bik(t)b (-1) kj (t) = ij (ij je tzv. Kroneckerovo delta, čiže ij = 0 pre i = j a ii = 1). Preto 0 = t {BB-1 }ij = n k=1 t bik(t)b (-1) kj (t) = = n k=1 bik(t) t b (-1) kj (t) + n k=1 bik(t) b (-1) kj (t) t , i, j = 1, 2, . . . , n, 49 čo v maticovom zápise je B t B-1 + B B-1 t = 0, čiže B-1 t = -B-1 B t B-1 . Pri dôkaze (iv) vyjdeme z rovnosti det B B = (det B)(2B-1 - diag B-1 ), kde diag B-1 = b (-1) 11 0 . . . 0 0 b (-1) 22 ... ... 0 0 . . . b (-1) nn (pozri [9], str. 98). Preto ln det B t = 1 det B det B t = 1 det B n i=1 i j=1 det B bij bij t = = 1 det B n i=1 i j=1 det B(2b (-1) ij - b (-1) ij ij) bij t = n i=1 n j=1 b (-1) ij bij t = tr B-1 B t (lebo B aj B-1 sú symetrické matice). Pomocou (i), (ii), (iii), a (iv) nájdeme maximum (69) riešením rovníc: (70) l = X V -1 y - X V -1 X = 0, (71) l 2 i = - 1 2 tr V -1 ZiZi + 1 2 (y - X) V -1 ZiZiV -1 (y - X) = 0, i = 1, 2, . . . , s + 1 (lebo V 2 i = ZiZi , i = 1, 2, . . . , s + 1). Vo všeobecnosti ^real, {2 1}real, . . . {2 s+1}real, ktoré maximalizujú (69) (pri realizácii y) sú riešením (70) a (71). Každé riešenie (70) a (71) nemusí však byť realizáciou ML­odhadu ^, 2 1, . . . , 2 s+1. Musíme sa ešte presvedčiť o matici druhých 50 derivácií (72) 2 l 2 l (2 1, . . . , 2 s+1) 2 l 2 1 ... 2 s+1 2 l 2 1 ... 2 s+1 (2 1, . . . , 2 s+1) , že je pre získané riešenia (70) a (71) negatívne definitná. 2 l je q × q matica, ktorej (i, j)-ty prvok je 2 l ij , 2 l (2 1, . . . , 2 s+1) je q × (s + 1) matica, ktorej (i, j)-ty prvok je 2 l i2 j , podobne 2 l 2 1 ... 2 s+1 a 2 l 2 1 ... 2 s+1 (2 1, . . . , 2 s+1) . Rovnako je nutné, aby (73) {2 s+1}real > 0, {2 i }real 0 i = 1, 2, . . . , s. V takomto prípade ^real, {2 1}real, . . . {2 s+1}real sú realizáciami ML­odhadu ^, 2 1 . . . , 2 s+1. Ak pre riešenia (70) a (71) neplatí (73), čo nastane keď niektoré {2 i }real, i {1, 2, . . . , s + 1} sú záporné, postupujeme tak, že vynecháme v modeli tie faktory, pre ktoré {2 i }real sú záporné a počítame ML­odhady v redukovanom modeli. Dá sa ukázať, že matica (72) je negatívne definitná v každom lokálnom maxime funkcie (69) pokiaľ toto maximum leží vo vnútri parametrického priestoru, t.j. ak pre riešenia (70) a (71) platí {2 1}real > 0, i = 1, 2, . . . , s + 1. Pre riešenie , 2 1, . . . 2 s+1 rovníc (70) a (71) platí By = (V -1 - V -1 X(X V -1 X)X V -1 )y = (74) = V -1 y - V -1 X(X V -1 X)X V -1 y = = V -1 y - V -1 X(X V -1 X)X V -1 X = V -1 (y - X) 51 (využijeme lemu 3 pre maticu V - 1 2 X, kde V - 1 2 =U 1 1 0 . . . 0 0 1 2 ... ... 0 0 . . . 1 n U , V = U 1 0 . . . 0 0 2 ... ... 0 0 . . . n U , I = UU , V - 1 2 = (V 1 2 )-1 , pozri [2], str. 64). Pretože pre i = 1, 2, . . . , s + 1 je (75) tr V -1 ZiZi = tr V -1 ZiZiV -1 V = tr V -1 ZiZiV -1 s+1 j=1 2 j ZjZj = = s+1 j=1 [tr(V -1 ZiZiV -1 ZjZj)]2 j , môžeme pomocou (74) a (75) písať ML­rovnice (70), (71) v ekvivalentnom tvare (76) X V -1 X = X V -1 y, (77) tr V -1 Z1Z1V -1 Z1Z1 tr V -1 Z1Z1V -1 Z2Z2 ... tr V -1 Z1Z1V -1 Zs+1Zs+1 tr V -1 Z2Z2V -1 Z1Z1 tr V -1 Z2Z2V -1 Z2Z2 ... tr V -1 Z2Z2V -1 Zs+1Zs+1 ... ... tr V -1 Zs+1Zs+1V -1 Z1Z1 tr V -1 Zs+1Zs+1V -1 Z2Z2 ... tr V -1 Zs+1Zs+1V -1 Zs+1Zs+1 × × 2 1 2 2 ... 2 s+1 = y (V-1 -V-1 X(X V-1 X)X V-1 )Z1Z1(V-1 -V-1 X(X V-1 X)X V-1 )y y (V-1 -V-1 X(X V-1 X)X V-1 )Z2Z2(V-1 -V-1 X(X V-1 X)X V-1 )y ... y (V-1 -V-1 X(X V-1 X)X V-1 )Zs+1Zs+1(V-1 -V-1 X(X V-1 X)X V-1 )y . ML­rovnice (70), (71) alebo (76), (77) sú nelineárne v parametroch 2 1, 2 2, . . . , 2 s+1 a riešia sa numericky (iteračne). Poznámka 48. Hartley a Rao v [3] zaviedli v modeli (55) transformáciu parametrov 2 1, 2 2, . . . , 2 s+1 na nové parametre 2 s+1, 1 = 2 1 2 s+1 , . . . , s = 2 s 2 s+1 . Preto V = 2 s+1H, V -1 = 1 2 s+1 H-1 . Toto vedie k preformulovaniu ML­rovníc a k ML­odhadom 2 s+1, ^1, . . . , ^s. ML - rovnice pre , 1, . . . , s, 2 s+1 sú X V -1 X = X V -1 y 52 2 s+1 = (y - X) H-1 (y - X) n tr H-1 ZiZi = (y - X) H-1 ZiZiH-1 (y - X) 2 s+1 i = 1, 2, . . . , s. Iteračné riešenie rovníc Hartleya a Raa môže byť v mnohých prípadoch jednoduchšie ako iteračné riešenie ML­rovníc (76) a (77). Poznámka 49. Pre ML-odhady ^, 2 1, . . . , 2 s+1 vieme určiť asymptotickú kovariančnú maticu (pre n ): cov( ^ 2 1 ... 2 s+1 ) = -E 2 l 2 l (2 1, . . . , 2 s+1) 2 l 2 1 ... 2 s+1 2 l 2 1 ... 2 s+1 (2 1, . . . , 2 s+1) = (78) = I-1 ^ 2 1 ... 2 s+1 , pričom v (78) je l(, 2 1, . . . , 2 s+1, Y ) náhodnou premennou, I ^ 2 1 ... 2 s+1 je informačná matica ML­odhadov ^, 2 1, . . . , 2 s+1. Platí cov( ^ 2 1 ... 2 s+1 ) = (X V -1 X)-1 0 0 2C-1 , kde {C}ij = tr V -1 ZiZiV -1 ZjZj i, j = 1, 2, . . . , s + 1. Samozrejme tu vzniká otázka, čo to znamená ,,n a ,,= . Hartley a Rao v [3] a Miller v [8] sa zaoberajú týmito problémami. Podrobnejšie tiež pozrite napr. v [13]. 53 11 ML­ODHADY VO VYVÁŽENOM MODELI S JEDNÝM NÁHODNÝM EFEKTOM Hľadajme ML­odhady ^, 2 1, 2 2 vo vyváženom modeli s jedným náhodným efektom, čiže v modeli (79) Y = Y11 Y12 ... Y1t Y21 Y22 ... Y2t ... YI1 YI2 ... YIt = (X ...Z1 ...Z2) b1 , X = 1n,1, Z1 = 1t,1 e1 1t,1 e2 ... 1t,1 eI , Z2 = In,n, cov(Y ) = 2 1Z1Z1 + 2 2I = V , Y je normálne rozdelený. (Ak preznačíme 2 2 = 2 0 , dostávame vyvážený model (40)). Aj v tomto jednoduchom modeli sú rovnice (76) a (77) komplikované. Jednoduchšie je nájsť riešenie rovníc (70) a (71). Platí [cov(Y )]-1 = V -1 = 1 2 2 I - 2 1 2 2(t2 1 + 2 2) Z1Z1, X V -1 X = 1 1 2 2 I - 2 1 2 2(t2 1 + 2 2) Z1Z1 1 = It (t2 1 + 2 2) a X V -1 Y = 1 1 2 2 I - 2 1 2 2(t2 1 + 2 2) Z1Z1 Y . Pre ^, 2 1 a 2 2 , ktoré sú riešením (70) platí (80) ^ = (X V -1 X)-1 X V -1 Y = (t2 1 + 2 2) It 1 2 2 1 - 2 1 2 2(t2 1 + 2 2) 1 Z1Z1 Y = 54 = (t2 1 + 2 2) It 1 2 2 1 - 2 1 2 2(t2 1 + 2 2) 1 Y = 1 It 1 Y = 1 It I i=1 t j=1 Yij. Lema 50. V modeli (79) platí tr V -1 Z1Z1 = tr Z1V -1 Z1 = It (t2 1 + 2 2) , (Y - X ^) V -1 Z1Z1V -1 (Y - X ^) = Y 1 (t2 1 + 2 2)2 Z1Z1 - 1 I 11 Y , kde ^ je dané vzťahom (80). Dôkaz. Lemu dokážeme jednoduchým dosadením a úpravami, pričom využijeme rovnosti uvedené pred vzťahom (80). Dôkaz urobte ako cvičenie. Dôsledok 51. Nech ^, 2 1 a 2 2 sú riešenia ML­rovníc (70) a (71) pre i = 1. Platí (81) t2 1 + 2 2 = 1 It Y Z1Z1 - 1 I 11 Y . Dôkaz. Vzťah (81) ľahko dostaneme ak využijeme výsledok lemy 50 a ML­rovnicu (71) pre i = 1. Dôkaz urobte ako cvičenie. Lema 52. V modeli (79) nech ^, 2 1 a 2 2 sú riešeniami rovníc (70) a (71) pre i = 1. Potom platí tr V -1 Z2Z2 = tr V -1 = It 2 2 (t - 1)2 1 + 2 2 t2 1 + 2 2 , (Y - X ^) V -1 Z2Z2V -1 (Y - X ^) = Y 1 4 2 I - 1 t Z1Z1 Y + I2 2 t2 1 + 2 2 . Dôkaz. Lemu dokážeme podobne ako lemu 50. Užitočné je udržiavať pokope výraz t2 1 + 2 2 . Dôkaz urobte ako cvičenie. Dôsledok 53. Nech ^, 2 1 a 2 2 sú riešeniami ML­rovníc (70) a (71). Platí 2 2 = Y 1 I(t - 1) I - 1 t Z1Z1 Y = (82) = 1 I(t - 1) I i=1 t j=1 Y 2 ij - 1 It(t - 1) I i=1 t j=1 Yij 2 . Dôkaz. Vzťah (82) dostaneme ak využijeme (80), (81) a lemu 52 v rovnici (71) pre i = 2. Dôkaz urobte ako cvičenie. 55 Pretože ML­odhad danej funkcie parametrov je ten istý ako táto funkcia ML ­ odhadov parametrov, pre ML­odhad 2 1 dostávame z (81) a (82) 2 1 = (t2 1 + 2 2) - 2 2 t = = 1 t 1 It Y Z1Z1Y - 1 I Y 11 Y Y - 1 I(t - 1) Y Y + 1 It(t - 1) Y Z1Z1Y = = Y 1 It(t - 1) Z1Z1 - I t - 1 It 11 Y = = 1 It(t - 1) I i=1 t j=1 Yij 2 - 1 It(t - 1) I i=1 t j=1 Y 2 ij(83) - 1 I2t2 I i=1 t j=1 Yij 2 . Poznámka 54. Samozrejme {2 1}real a {2 2}real sú realizácie ML­odhadov len ak {2 1}real 0 a {2 2}real > 0. Poznámka 55. ML­odhad 2 2 je ten istý ako odhad (48) navrhnutý Hendersonom, či MINQUE odhad parametra 2 2 . ML­odhad 2 1 je rôzny od Hendersonovho odhadu parametra 2 1 (Hendersonov odhad je ten istý ako MINQUE odhad). Samozrejme ML-odhad 2 2 je preto nevychýleným odhadom 2 2 . Počítajme teraz strednú hodnotu ML­odhadu 2 1 . Podľa lemy 10 je E(2 1) = 1 It(t - 1) 2 1 Z1Z1 - I t - 1 It 11 1+ + tr Z1Z1 - I t - 1 It 11 (2 2I + 2 1Z1Z1) = (84) = I - 1 I 2 1 - 2 2 It . Poznámka 56. Z (84) vyplýva, že 2 1 nie je nevychýleným odhadom 2 1 . Lema 57. Pre kovariančnú maticu ML­odhadov (2 1 , 2 2) platí (85) cov( 2 1 2 2 ) = 24 2 It2(t - 1) + 2(t2 1 + 2 2)2 It2 I - 1 I - 4 2 It(t - 1) - 4 2 It(t - 1) 24 2 I(t - 1) . Dôkaz. Lemu dokážeme pomocou vety 18 (pozri tiež cvičenie 1 na str. 32). Opäť odporúčame udržiavať pokope výraz t2 1 + 2 2 . Dôkaz urobte ako cvičenie. ML­odhady v zložitejších modeloch ako aj iné typy odhadov nájdete napr. v [13]. 56 12 ODHADY ZÍSKANÉ METÓDOU NAJMENŠÍCH ŠTVORCOV Majme zmiešaný lineárny model Y = X + Z1b1 + Z2b2 + . . . + Zsbs + Zs+1bs+1 s rovnakým označením ako v (55) (resp. (15), pričom v (15) predpokladáme I = = Zs+1, = bs+1, 2 s+1 = 2 0 ). Samozrejme cov(Y ) = 2 1Z1Z1 + 2 2Z2Z2 + . . . + 2 s ZsZs + 2 s+1Zs+1Zs+1 = s+1 i=1 2 i Vi. Definícia 58. Povieme, že ^ je odhad pevných efektov získaný metódou najmenších štvorcov, ak X ^ - Y 2 = inf Rq X - Y 2 , pričom je Euklidovská norma vektora, t.j. z = z z. Veta 59. Ak ^ je odhad pevných efektov získaný metódou najmenších štvorcov, tak X ^ = X(X X)X Y a pre lineárne nevychýlene odhadnuteľnú funkciu p je disperzia D(p ^) = p (X X)- X s+1 i=1 2 i ViX(X X)- p. Dôkaz. Platí X - Y 2 = (X - Y ) (X - Y ) = X X - 2Y X + Y Y . Nech inf Rq X - Y 2 = X[(X X)X Y + ] - Y 2 = = (X(X X)X Y + X - Y ) (X(X X)X Y + X - Y ) = = Y X(X X)X X(X X)X Y + Y X(X X)X X - Y X(X X)X Y + + X X(X X)X Y + X X- X Y -Y X(X X)X Y -Y X+Y Y = = X(X X)X Y - Y 2 + X 2 , teda vskutku X ^ = X(X X)X Y . Podľa lemy 1 je p lineárne nevychýlene odhadnuteľná práve vtedy ak p (X ), teda w Rn , že p = X w. Preto D(p ^) = D(w X ^) = D(w X(X X)X Y ). Dôkaz už ľahko dokončíme. 57 Definícia 60. Odhad ^ = (2 1, . . . , 2 s+1) , pre ktorý platí s+1 i=1 2 i Vi - 2 = min c1,...,cs+1 s+1 i=1 ciVi - 2 , kde = (X ^ - Y )(X ^ - Y ) , ( ^ je odhad pevných efektov získaný metódou najmenších štvorcov; je tzv. predbežný odhad matice cov(Y ); . je Euklidovská norma matice, teda A = = tr AA ) sa nazýva odhad variančných komponentov získaný metódou najmenších štvorcov. Veta 61. Označme q = tr V1 tr V2 ... tr Vs+1 , Q = tr V1V1 tr V1V2 . . . tr V1Vs+1 tr V2V1 tr V2V2 . . . tr V2Vs+1 ... tr Vs+1V1 tr Vs+1V2 . . . tr Vs+1Vs+1 . ^ = (2 1, . . . , 2 s+1) je odhad variančných komponentov získaný metódou najmenších štvorcov práve vtedy ak (86) Q^ = q. Dôkaz. Hľadajme min c1,...,cs+1 s+1 i=1 ciVi - 2 = = min c1,...,cs+1 tr( s+1 i=1 ciVi - )( s+1 j=1 cjVj - ) = = min c1,...,cs+1 s+1 i=1 s+1 j=1 cicjtrViVj - 2 s+1 i=1 ci tr Vi + tr = = min cRs+1 {c Qc - 2c q + tr } = (c). Nutná podmienka existencie minima je c = 2Qc - 2q = 0, čiže ak ^ je odhadom variančných komponentov získaný metódou najmenších štvorcov, tak Q^ = q. 58 Naopak, majme u(0) Rs+1 pre ktorý platí Qu(0) = q a ľubovoľný z Rs+1 . Potom min c1,...,cs+1 s+1 i=1 ciVi - 2 = = min 1,...,s+1 s+1 i=1 (u (0) i + i)Vi - 2 = = min 1,...,s+1 tr( s+1 i=1 (u (0) i + i)Vi - )( s+1 j=1 (u (0) j + j)Vj - ) = = min Rs+1 {(u(0) + ) Q(u(0) + ) - 2(u(0) + ) q + tr } = = min Rs+1 {u(0) Qu(0) + 2u(0) Q + Q - 2q u(0) - 2q + tr } = = min Rs+1 { Q - u(0) Qu(0) + tr } = -u(0) Qu(0) + tr , z čoho vyplýva, že minimum je dosiahnuté pre vektor ^ = u(0) +, pričom Q = 0 Q = 0 (lebo Q je p.s.d. matica). Teda pre odhad ^ variančných komponentov získaný metódou najmenších štvorcov platí Q^ = Q(u(0) + ) = q. Dôsledok 62. Odhad variančných komponentov získaný metódou najmenších štvorcov existuje práve vtedy ak q (Q). Dôsledok 63. Ak V1, . . . , Vs+1 sú lineárne nezávislé, tak odhad 2 variančných komponentov získaný metódou najmenších štvorcov je ^ = Q-1 q, kde q = Y MXV1MXY Y MXV2MXY ... Y MXVs+1MXY , MX = I - X(X X)X . Dôkaz. Oprieme sa o fakt, že V1, . . . , Vs+1 sú lineárne nezávislé práve vtedy, ak Q je regulárna matica a využijeme jednoduchú úpravu qi = tr Vi = tr(Y - X ^)(Y - X ^) Vi = = (Y - X ^) Vi(Y - X ^) = Y (I - X(X X)X )Vi(I - X(X X)X )Y = = Y MXViMXY , i = 1, 2, . . . , s + 1. Dôkaz urobte ako cvičenie. 59 Poznámka 64. Ak p = (p1, p2, . . . , ps+1) Rs+1 , tak odhad s+1 i=1 pi2 i = p ^, kde Q^ = q. Veta 65. Ak V1, V2, . . . , Vs+1 sú lineárne nezávislé, tak (87) E(2 i ) = (2 1, 2 2, . . . , 2 s+1)RQ-1 ei, i = 1, 2, . . . , s + 1, kde R = tr MXV1MXV1 tr MXV1MXV2 . . . tr MXV1MXVs+1 tr MXV2MXV1 tr MXV2MXV2 . . . tr MXV2MXVs+1 ... ... tr MXVs+1MXV1 tr MXVs+1MXV2 . . . tr MXVs+1MXVs+1 . Dôkaz. Podľa dôsledku 63 je pre i = 1, 2, . . . , s + 1 2 i = eiQ-1 q = q Q-1 ei, teda pomocou lemy 10 E(2 i ) = (E(Y MXV1MXY ), . . . , E(Y MXVs+1MXY ))Q-1 ei = = tr MXV1MX s+1 i=1 2 i Vi, . . . , tr MXVs+1MX s+1 i=1 2 i Vi Q-1 ei = = s+1 i=1 2 i tr MXV1MXVi, . . . , s+1 i=1 2 i tr MXVsMXVi Q-1 ei = = (2 1, 2 2, . . . , 2 s+1)RQ-1 ei. Veta 66. Ak Y N X, s+1 i=1 2 i Vi a V1, . . . , Vs+1 sú lineárne nezávislé, je D(2 i ) = 2eiQ-1 W Q-1 ei, i = 1, 2, . . . , s + 1, kde {W }i,j = s+1 k=1 s+1 l=1 2 k2 l tr MXViMXVkMXVjMXVl, i, j = 1, 2, . . . , s + 1. Dôkaz. Počítajme pre i, j = 1, 2, . . . , s + 1 použitím vety 18 a lemy 10 cov(Y MXViMXY , Y MXVjMXY ) = 60 = E(Y MXViMXY Y MXVjMXY )-E(Y MXViMXY )E(Y MXVjMXY ) = = 2 s+1 k=1 s+1 l=1 2 k2 l tr MXViMXVkMXVjMXVl+ + tr MXViMX s+1 k=1 kVk tr MXVjMX s+1 l=1 lVl- tr MXViMX s+1 k=1 kVk tr MXVjMX s+1 l=1 lVl = = 2 s+1 k=1 s+1 l=1 2 k2 l tr MXViMXVkMXVjMXVl. Preto cov(q) = 2W . Teraz už ľahko dostávame (88) D(2 i ) = D(eiQ-1 q) = eiQ-1 cov(q)Q-1 ei = 2eiQ-1 W Q-1 ei. Iné vlastnosti odhadov získaných metódou najmenších štvorcov pozri napr. v [15]. 13 ODHADY ZÍSKANÉ METÓDOU NAJMENŠÍCH ŠTVORCOV V MODELI S JEDNÝM NÁHODNÝM EFEKTOM Majme model (40), teda Yij = + bi + ij, i = 1, 2, . . . , I , j = 1, 2, . . . , ni , kde náhodný vektor b1 = b1 b2 ... bI N(0I,1; 2 1I) a chybový vektor N(0, 2 2I) sú nezávislé. Model zapíšeme ako Yn,1 = (X ...Z1 ...Z2) b1 , 61 kde X = 1n,1, Z1 = 1n1,1 e1 1n2,1 e2 ... 1nI ,1 eI , Z2 = In,n, n = I i=1 ni (pozri aj kapitolu 9), E(Y ) = X = 1, cov(Y ) = 2 1Z1Z1 + 2 2Z2Z2 = 2 1V1 + 2 2V2 = 2 1V1 + 2 2I. Pomocou vety 61, dôsledku 63, poznámky 64, vzorcov (87) a (88) spočítame odhady získané metódou najmenších štvorcov a určíme ich stredné hodnoty a disperzie. Platí tr V1V1 = tr Z1Z1Z1Z1 = tr Z1Z1Z1Z1 = I i=1 n2 i , tr V1V2 = tr V1I = tr Z1Z1 = I i=1 ni = n, tr V2V2 = tr I I = I i=1 ni = n. Preto pre maticu Q z vety 61 a jej inverziu platí Q = I i=1 n2 i n n n , (89) Q-1 = 1 n( I i=1 n2 i - n) n -n -n I i=1 n2 i . Ďalej budeme potrebovať vyjadriť maticu MX (z dôsledku 63): MX = I - X(X X)X = I - 1 n 1n,111,n. Spočítame aj MXV1MX = (I - 1 n 11 )Z1Z1(I - 1 n 11 ) = = Z1Z1 - 1 n 11 Z1Z1 - 1 n Z1Z111 + 1 n2 11 Z1Z111 . Využijúc výpočty z kapitoly 6 dostávame (90) Y MXV1MXY = 62 = Y (Z1Z1 - 1 n 11 Z1Z1 - 1 n Z1Z111 + 1 n2 11 Z1Z111 )Y = = I i=1 ni j=1 Yij 2 - 2 n I i=1 ni j=1 Yij I k=1 nk nk l=1 Ykl + 1 n2 I i=1 n2 i I k=1 nk l=1 Ykl 2 , (91) Y MXV2MXY = Y MXY = Y (I - 1 n 11 )Y = = I i=1 ni j=1 Y 2 ij - 1 n I i=1 ni j=1 Yij 2 . Podľa dôsledku 63, poznámky 64 a vzorcov (89), (90), (91) sa odhady jednotlivých variančných komponentov získané metódou najmenších štvorcov rovnajú (92) 2 1 = e1Q-1 q = n n( I i=1 n2 i - n) (1, -1) Y MXV1MXY Y MXV2MXY = = 1 I i=1 n2 i - n I i=1 ni j=1 Yij 2 - 2 n I i=1 ni j=1 Yij I k=1 nk nk l=1 Ykl+ + I i=1 n2 i + n n2 I i=1 ni j=1 Yij 2 - I i=1 ni j=1 Y 2 ij , (93) 2 2 = - 1 I i=1 n2 i - n - I i=1 ni j=1 Yij 2 + 2 n I i=1 ni j=1 Yij I k=1 nk nk l=1 Ykl- - 2 I i=1 n2 i n2 I i=1 ni j=1 Yij 2 + I s=1 n2 s n I i=1 ni j=1 Y 2 ij . Pre vyvážený model, teda v prípade n1 = n2 = . . . = nI = t dostávame z (92) (94) 2 1 = 1 It(t - 1) I i=1 t j=1 Yij 2 t - 1 It I i=1 t j=1 Yij 2 - I i=1 t j=1 Y 2 ij a z (93) (95) 2 2 = 1 It(t - 1) - I i=1 t j=1 Yij 2 + t I i=1 t j=1 Y 2 ij = 63 = 1 I(t - 1) I i=1 t j=1 Y 2 ij - I i=1 1 t t j=1 Yij 2 . Z (94) je vidieť, že odhad 2 1 získaný metódou najmenších štvorcov je vo vyváženom modeli s jedným náhodným efektom ten istý ako ML­odhad a rôzny od Hendersonovho (ktorý je totožný v tomto prípade s MINQUE odhadom). Pozri tiež poznámku 55. Z (95) vidíme, že vo vyváženom modeli s jedným náhodným efektom sú Hendersonov, MINQUE, ML aj odhad získaný metódou najmenších štvorcov parametra 2 2 totožné (v prípade Hendersonovho odhadu je označený ako 2 0 ). Je preto zrejmé, že odhad 2 2 získaný metódou najmenších štvorcov vo vyváženom modeli s jedným náhodným efektom, je nevychýlený. Počítajme teraz strednú hodnotu odhadu 2 1 získaného metódou najmenších štvorcov. Podľa vety 65, vzťahu (87) je E(2 1) = (2 1, 2 2)RQ-1 e1, kde (96) {R}11 = tr MXV1MXV1 = tr MXZ1Z1MXZ1Z1 = = tr n1 0 . . . 0 0 n2 ... ... 0 0 . . . nI - 1 n n1 n2 ... nI (n1, n2, . . . , nI) × × n1 0 . . . 0 0 n2 ... ... 0 0 . . . nI - 1 n n1 n2 ... nI (n1, n2, . . . , nI) = = I i=1 n2 i - 2 n I i=1 n3 i + 1 n2 I i=1 n2 i 2 = = 1 n I i=1 ni I j=1 n2 j + 1 n I i=1 n2 i 2 - 2 I i=1 n3 i , (97) {R}12 = {R}21 = tr MXZ1Z1MXI = I i=1 ni - 1 n I i=1 n2 i a (98) {R}22 = tr MXV2MXV2 = tr MXMX = n - 1. 64 Pomocou (89), (96), (97) a (98) dostávame (99) E(2 1) = (2 1, 2 2) 1 I i=1 n2 i - n × × 1 n n I i=1 n2 + 1 n I i=1 n2 i 2 - 2 I i=1 n3 i n - 1 n I i=1 n2 i n - 1 n I i=1 n2 i n - 1 1 -1 = = 1 n( I i=1 n2 i - n) (2 1, 2 2)× × (n + 1) I i=1 n2 i + 1 n I i=1 n2 i 2 - 2 I i=1 n3 i - n2 n - I i=1 n2 i = = 2 1 n( I i=1 n2 i - n) (n + 1) I i=1 n2 i + 1 n I i=1 n2 i 2 - 2 I i=1 n3 i - n2 - 2 2 n . Spočítajme ešte rovnakým spôsobom (100) E(2 2) = (2 1, 2 2) 1 n( I i=1 n2 i - n) × × n I i=1 n2 i + 1 n I i=1 n2 i 2 - 2 I i=1 n3 i n2 - I i=1 n2 i n2 - I i=1 n2 i n(n - 1) -1 1 n I i=1 n2 i = = 1 n( I i=1 n2 i - n) (2 1, 2 2) 2 I i=1 n3 i + 2 n I i=1 n2 i 2 n( I i=1 n2 i - n) = = 22 1 n( I i=1 n2 i - n) I i=1 n3 i - 1 n I i=1 n2 i 2 + 2 2. Vo vyváženom modeli dostávame z (99) (101) E(2 1) = 1 It(It2 - Tt) (It + 1)It2 + I2 t4 It - 2It3 - I2 t2 - 2 2 It = = I - 1 I 2 1 - 1 It 2 2 (samozrejme vo vyváženom modeli je E(2 2) = 2 2 ). Teraz spočítajme disperzie odhadov variančných komponentov získaných metódou najmenších štvorcov v modeli s jedným náhodným efektom. Potrebujeme pritom ešte nasledujúcu lemu. 65 Lema 67. V modeli s jedným náhodným efektom pre maticu W z vety 66 platí {W }11 = 4 1 tr(Z1MXZ1)4 + 22 12 2 tr(Z1MXZ1)3 + 4 2 tr(Z1MXZ1)2 , {W }12 = {W }21 = 4 1 tr(Z1MXZ1)3 + 22 12 2 tr(Z1MXZ1)2 + 4 2 tr Z1MXZ1, {W }22 = 4 1 tr(Z1MXZ1)2 + 22 12 2 tr Z1MXZ1 + 4 2 tr MX, tr MX = n - 1, tr Z1MXZ1 = n - 1 n I i=1 n2 i , tr(Z1MXZ1)2 = 1 n n I i=1 n2 i + 1 n I i=1 n2 i 2 - 2 I i=1 n3 i , tr(Z1MXZ1)3 = I i=1 n3 i - 3 n I i=1 n4 i + 3 n2 I i=1 n2 i I j=1 n3 j - 1 n3 I i=1 n2 i 3 , tr(Z1MXZ1)4 = I i=1 n4 i - 4 n I i=1 n5 i + 4 n2 I i=1 n2 i I j=1 n4 j + 2 n2 I i=1 n3 i 2 - - 4 n3 I i=1 n2 i 2 I j=1 n3 j + 1 n4 I i=1 n2 i 4 . Dôkaz vykonáme jednoduchým, ale zdĺhavým dosadzovaním a spočítavaním. Urobte ho ako cvičenie. Počítajme teraz podľa vety 66 (využijúc lemu 67): (102) D(2 1) = 2e1Q-1 W Q-1 e1 = 2 ( I i=1 n2 i - n)2 ({W }11 - 2{W }12 + {W }22) = = 2 ( I i=1 n2 i - n)2 {4 1[tr(Z1MXZ1)4 - 2 tr(Z1MXZ1)3 + tr(Z1MXZ1)2 ]+ +22 12 2[tr(Z1MXZ1)3 - 2 tr(Z1MXZ1)2 + tr Z1MXZ1]+ +4 2[tr(Z1MXZ1)2 - 2 tr Z1MXZ1 + tr MX]}. Podobne (103) D(2 2) = 2e2Q-1 W Q-1 e2 = = 2 ( I i=1 n2 i - n)2 {W }11 - 2 n I i=1 n2 i {W }21 + 1 n I i=1 n2 i 2 {W }22 = 66 = 2 ( I i=1 n2 i - n)2 4 1 tr(Z1MXZ1)4 - 2 n I i=1 n2 i tr(Z1MXZ1)3 + + 1 n I i=1 n2 i 2 tr(Z1MXZ1)2 + 22 12 2 tr(Z1MXZ1)3 - - 1 n I i=1 n2 i tr(Z1MXZ1)2 + 1 n I i=1 n2 i 2 tr Z1MXZ1 + + 4 2 tr(Z1MXZ1)2 - 2 n I i=1 n2 i tr Z1MXZ1 + 1 n I i=1 n2 i 2 tr MX . Lema 68. Vo vyváženom modeli s jedným náhodným efektom je D(2 1) = 24 2 It2(t - 1) + 2(t2 1 + 2 2)2 It2 I - 1 I , D(2 2) = 4 2 I(t - 1) . Dôkaz. Lemu dokážeme využitím lemy 67 a dosadením n1 = n2 = . . . = nI = t do (102) a (103). Dôkaz urobte ako cvičenie. Výsledky lemy 68 porovnajte s lemou 57 a cvičením 1 na str. 32. Pre ďalšie štúdium problematiky odhadov variančných komponentov odporúčame napr. monografiu [13]. 67 Literatúra [1] ADAM, J. ­ SCHARF, J. H. ­ ENKE, H. : Methoden der Statistischen Analyse in Medizin und Biologie. VEB Verlag Volk und Gesundheit, Berlin 1971. [2] ANDĚL, J. : Matematická statistika. SNTL/ALFA, Praha 1978. [3] HARTLEY, H. O. ­ RAO, C. R. : Maximum likelihood estimation for the mixed analysis of variance model. Biometrika 54 (1967), 93­108. [4] HENDERSON, C. R. : Estimation of variance and covariance components. Biometrics 9 (1953), 226­252. [5] KUBÁČEK, L. : Základy teórie odhadu. VEDA, Bratislava 1983. [6] KUBÁČEK, L. : Foundations of estimation theory. Elsevier, Amsterdam 1988. [7] MELOUN, M. ­ MILITKÝ, J. : Statistické zpracování experimentálních dat. PLUS, Praha 1994. [8] MILLER, J. J. : Asymptotic properties of maximum likelihood estimates in the mixed model of the analysis of variance. Annals of Statistics 5 (1977), 746­762. [9] RAO, C. R. : Lineární metody statistické indukce a jejich aplikace. ACADEMIA, Praha 1978. [10] RAO, C. R. ­ MITRA, S. K. : Generalized inverse of matrices and its applications. Wiley, New York 1971. [11] ROHDE, C. R. : Special applications of the theory of generalized matrix inversion to statistics. Proc. Symp. Theo. Appl. Generalized Inverses of Matrices, Mathematics Ser. No. 4, Texas Technological College, Lubbock (1968). [12] SEARLE, S. R. : Another look at Henderson's method of estimating variance components. Biometrics 24 (1968), 749­778. [13] SEARLE, S. R. ­ CASELLA, G. ­ McCULLOCH, Ch. E. : Variance components. Wiley, New York 1992. [14] SCHEFFE, H. : The analysis of variance. Wiley, New York 1959. [15] ŠTULAJTER, F. : Consistency of linear and quadratic least squares estimators in regression models with covariance stationary errors. Applications of Mathematics 36 (1991), 149­155. 68 Obsah 1 Úvod. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2 Model analýzy rozptylu ­ model s pevnými efektmi . . . . . . . . . . . . . . . . . . . 3 3 Model s náhodnými efektmi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 4 Model so zmiešanými efektmi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 5 Odhady v zmiešanom lineárnom modeli podľa Hendersona . . . . . . . . . . . . 14 6 Odhady variančných komponentov v modeli s jedným náhodným efektom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 7 Odhady variančných komponentov vo vyváženom modeli s dvoma náhodnými efektmi a interakciou. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 8 MINQUE odhady variančných komponentov . . . . . . . . . . . . . . . . . . . . . . . . . . 38 9 MINQUE odhady variančných komponentov v modeli s jedným náhodný efektom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 10 Odhady variančných komponentov získané metódou maximálnej vierohodnosti (ML­odhady) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 11 ML-odhady vo vyváženom modeli s jedným náhodným efektom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 12 Odhady získané metódou najmenších štvorcov . . . . . . . . . . . . . . . . . . . . . . . . 57 13 Odhady získané metódou najmenších štvorcov v modeli s jedným náhodným efektom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 Literatúra ......................................................................... 68 69