ROBUST’2004
c JČMF 2004
ROBUSTNOST V MODELU RŮSTOVÝCH KŘIVEK Ivan Žežula, Daniel Klein Klíčová slova: Mnohorozměrný lineární model, simulace. Abstrakt: Příspěvek se zabývá dvěmi variantami modelu růstových křivek – replikovaným modelem a modelem se speciálními variančními strukturami – a zkoumá chování odhadů parametrů v případě porušení předpokladů, zejména normality rozdělení.
1
Motivace a popis modelu
Model růstových křivek vznikl při analýze anatomických dat: na jisté zubní klinice byl sledován vývoj vzdálenosti středu hypofýzy od pterygomaxilární brázdy u chlapců a děvčat. Otázkou bylo, zda tato vzdálenost je u děvčat menší než u chlapců a zda rychlost růstu je stejná. Získaná data jsou na obrázku. Silnou čarou je vyznačen průměr; jeho průběh ukazuje, že růst má zhruba lineární trend.
Na první pohled to vypadá jako dvě nezávislé regrese. Potthoff a Roy si však uvědomili, že oba soubory – vzhledem k okolnostem vzniku – mají stejnou varianční matici. Proto navrhli společný model, který tuto skutečnost zohledňoval: E Y = XBZ
var(vec Y ) = Σ ⊗ I
Tento model přirozeným způsobem spojuje regresní analýzu s analýzou rozptylu – matice X je maticí analýzy rozptylu, Z maticí regresních konstant. Neznámé parametry jsou v maticích B (regresní koeficienty pro jednotlivé skupiny) a Σ (variance a kovariance pro jednotlivé časy pozorování). Později se tento model dočkal mnoha dalších aplikací i rozšíření. Z těchto aplikací vyplynuly různé struktury varianční matice Σ :
444
Ivan Žežula, Daniel Klein
• obecná p.d. matice • Σ=
P
θi Vi , Vi známé
• rovnoměrná varianční struktura, t.j. Σ = σ 2 ((1 − ρ) I + ρ11′ ) • seriální varianční struktura, t.j.
1 ρ .. .
Σ = σ2 ρp−1
ρ 1 .. .
... ... .. .
ρp−2
...
ρp−1 ρp−2 .. . 1
P Všimněme si, že rozklad Σ = θi Vi je zcela obecný a zahrnuje P všechny ′ ostatní případy. Např. obecnou matici Σ lze napsat ve tvaru Σ = i σii ei ei + P ′ ′ i<j σij ei ej + ej ei , kde ei jsou jednotkové vektory (t.j. obsahující 1 na i-tém místě a 0 jinde). Při praktickém použití modelu se však setkáváme s některými problémy. K odhadu B totiž potřebujeme znalost Σ: ˆ = (X ′ X)−1 X ′ Y Σ−1 Z ′ ZΣ−1 Z ′ B
−1
.
To přináší dva okruhy otázek: Jednak, jak odhadnout Σ resp. její komponenty? A následně, jak se změní jeho vlastnosti odhadu B, jestliže Σ odhadneme? Specielně, o kolik vzroste jeho rozptyl? V tomto směru již bylo dosaženo mnoha výsledků. Především je známo, že v obecné situaci, při rozkladu Σ na varianční komponenty, existují stejnoměrně nejlepší odhady jen v triviálních případech. Tedy ve většině případů můžeme použít buď jen maximálně věrohodné odhady (pro něž většinou neexistují explicitní vzorce, počítají se pouze numericky a známe jen asymptotické vlastnosti) anebo lokálně nejlepší odhady (u nichž sice známe explicitní vzorce, ale zase musíme dobře vědět, v kterém bodě parametrického prostoru je počítat). Stejnoměrně nejlepší odhad existuje v jednom důležitém případě: když Σ je zcela neznámá. Částečná znalost varianční struktury situaci komplikuje. Lokálně nejlepší odhady variančních komponent θ mohou být navíc závislé na hodnotě B – dochází pak ke kruhové závislosti. Proto důležitosti nabývají odhady invariantní, které na B nezávisí; i pro ně jsou známé explicitní vzorce, odhady však nemusí vždy existovat. V případě odhadování pouze σ 2 a ρ známe nestranné odhadové rovnice pro oba parametry. Ve všech těchto případech je odhad B asymptoticky nestranný. Podrobněji o těchto modelech viz např. [1], [3].
445
Robustnost v modelu růstových křivek
2
Replikovaný model
Tento model byl vytvořen k oslabení závislosti lokálně nejlepších odhadů na volbě počátečního bodu. Obsahuje nezávislá opakování měření ze základního modelu: Yj = XBZ + ej , j = 1, . . . , s Varianční struktura přitom zůstává stejná. Odhad B je také prakticky stejný, jako v základním modelu: ˆ = (X ′ X)−1 X ′ Y¯ Σ−1 Z ′ ZΣ−1 Z ′ B
−1
.
Odhady variančních komponent v tomto modelu vždy existují. Výsledky z tohoto modelu jsou shrnuty v [4]. Dřívější simulace (viz [5]) ukázaly rychlou konvergenci odhadů k známému asymptotickému rozdělení. Nás zajímal problém robustnosti odhadů v tomto modelu. Všechny známé výsledky jsou totiž založeny na předpokladu normality. Je tedy přirozené se ptát, co se stane, když rozdělení chyb není normální? Podobnou otázkou je: co se stane s odhady, když chyby nemají nulovou střední hodnotu, t.j. když pozorování obsahují systematickou chybu? K tomuto účelu jsme provedli rozsáhlou simulační studii. Uvažovali jsme přitom různě složité modely střední hodnoty (polynomy 1. až 3. stupně) a různě složité modely varianční struktury: 1. Σ obecná p.d. matice 2. Σ s konstantními diagonálami 3. Σ se dvěmi komponentami (Σ = θ1 V + θ2 I) Uvažovaná chybová rozdělení byla následující: 1. normální, ale se systematickou chybou (µ 6= 0) 2. směs normálního a exponenciálního
3. různé formy beta rozdělení
446
Ivan Žežula, Daniel Klein
4. Laplaceovo rozdělení
Všechna rozdělení s výjimkou 1) byla posunuta tak, aby měla nulovou střední hodnotu. Výsledky simulací: Následující vybrané grafy zobrazují empirické rozdělení statistiky χ2 (kteb b k jeho skutečné hodnotě) pro různé rá měří celkovou vzdálenost od odhadu B hodnoty s. Počet replikací roste zezadu dopředu, zcela vpředu je asymptotické rozdělení.
Kvadratická závislost, Σ obecná, normální rozdělení
Kvadratická závislost, Σ se 2 komponentami, normální rozdělení
447
Robustnost v modelu růstových křivek
Lineární závislost, Σ obecná, rozdělení B(0.5,6)
Kubická
závislost,
Σ
obecná, rozdělení B(2,2)
Lineární závislost, Σ se 2 komponentami, rozdělení B(0.5,6)
Kubická závislost, Σ s
Kubická závislost, Σ se 2
konstatními diagonálami, rozdělení B(2,2)
komponentami, rozdělení B(2,2)
V případě nenulovosti střední hodnoty chyby – t.j. existence systematické chyby – grafy neuvádíme, jelikož odhady ležely daleko od skutečných hodnot a s rostoucím počtem replikací divergovaly do nekonečna. Lze tedy udělat následující závěry: • Model je silně citlivý na přítomnost systematické chyby • Ve všech ostatních případech je konvergence k asymptotickému rozdělení velmi rychlá • Rychlost konvergence je nepřímo úměrná počtu odhadovaných parametrů
3
Speciální varianční struktury
Dále jsme zkoumali robustnost základního modelu s rovnoměrnou a seriální strukturou (t.j. bez replikací). Odhady σ ˆ 2 a ρˆ pro oba modely jsou uvedeny v [6]; porovnej také s [2]. V obou modelech nás především zajímalo, jak se projeví různá chybová rozdělení na střední kvadratické chybě (MSE) odhadu sledovaných parametrů σ 2 a ρ. V modelu s rovnoměrnou strukturou lze pro námi uvažované odhady odvodit, že platí:
448
Ivan Žežula, Daniel Klein
MSE σ ˆ2
2σ 4 1 + (p − 1)ρ2 · n − r(X) p 2 (1 − ρ)2 (1 + (p − 1)ρ)2 · + o n−1 n − r(X) p (p − 1)
= var σ ˆ2 =
MSE ρˆ =
Přirozenou otázkou bylo, jestli lze odhad MSE, který vznikne dosazením odhadnutých hodnot σ ˆ 2 a ρˆ do tohoto teoretického vzorce, je prakticky použitelný, t.j. jestli se příliš neliší od skutečnosti. V obou modelech simulace ukázaly, že odhady parametrů σ 2 i ρ jsou dostatečně robustní a téměř nezávisí na chybovém rozdělení. Největší vliv na MSE těchto odhadů měla hodnota ρ. Zjistili jsme, že 1. chyba odhadu σ ˆ 2 roste s |ρ|
Model s rovnoměrnou strukturou
Model se seriální strukturou
2. chyba odhadu ρˆ klesá s |ρ|
Model s rovnoměrnou strukturou
Model se seriální strukturou
Trochu jiná situace je u aproximace MSE pomocí odhadnutých hodnot parametrů. Zde jsme zjistili, že
449
Robustnost v modelu růstových křivek 1. aproximace MSE ρˆ je robustní
ρ malé
ρ velké
2. aproximace MSE σ ˆ 2 je citlivá na změnu rozdělení, hlavně pro malá ρ
ρ malé
ρ velké
Reference [1] von Rosen D. (1989). Maximum likelihood estimators in multivariate linear normal models. Journal of Multivariate Analysis 31, 187 – 200. [2] Lee J. C. (1988). Prediction and estimation of growth curves with special covariance structures. JASA 83, No. 402, 432 – 440. [3] Žežula I. (1993). Covariance components estimation in the growth curve model. Statistics 24, 321 – 330. [4] Žežula I. (1997). Asymptotic properties of the growth curve model with covariance components. Applications of Mathematics 42, No. 1, 57 – 69. [5] Žežula I. (1996). Simulation study in the growth curve model. Tatra Mountains Mathematical Publications 7, 183 – 188. [6] Žežula I. (2004). Special variance structures in the growth curve model. To appear. Adresa: I. Žežula, D. Klein, UMV PF UPJŠ, Jesenná 5, 041 54 Košice, SR E-mail :
[email protected],
[email protected]
450
