28

´ ANALYZA ˇ ´ VÍCEROZMERN YCH DAT ˇ ˇ AV ´ AN ´ Í V AKREDITOVANYCH ´ URCENO PRO VZDEL STUDIJNÍCH PROGRAMECH

JOSEF TVRDÍK

ˇ ÍSLO OPERACN ˇ ÍHO PROGRAMU: CZ.1.07 C ´ ˇ ÍHO PROGRAMU: NAZEV OPERACN ˇ AV ´ AN ´ Í PRO KONKURENCESCHOPNOST VZDEL ˇ Í: 7.2 OPATREN ˇ ÍSLO OBLASTI PODPORY: 7.2.2 C ´ ´ ˇ ˇ U ˚ VE INOVACE VYUKY INFORMATICKYCH PREDM ET ´ UNIVERZITY STUDIJNÍCH PROGRAMECH OSTRAVSKE ˇ Í C ˇ ÍSLO PROJEKTU: CZ.1.07/2.2.00/28.0245 REGISTRACN

OSTRAVA 2013

Tento projekt je spolufinancován Evropsk´ ym sociáln´ım fondem a státn´ım rozpoˇctem ˇ e republiky Cesk´

Recenzent: Prof. RNDr. Ing. Ivan Kˇriv´ y, CSc.

Název: Autor: Vydán´ı: Poˇcet stran:

Anal´ yza v´ıcerozmˇern´ ych dat Josef Tvrd´ık tˇret´ı, 2013 128

Jazyková korektura nebyla provedena, za jazykovou stránku odpov´ıdá autor.

c Josef Tvrd´ık ⃝ c Ostravská univerzita v Ostravˇe ⃝

OBSAH

i

Obsah ´ 1 Uvod

3

2 Vektory a matice

5

2.1

Základn´ı pojmy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

2.2

Vlastn´ı ˇc´ısla a vlastn´ı vektory matice . . . . . . . . . . . . . . . . . .

8

2.3

Dalˇs´ı d˚ uleˇzité vlastnosti matic . . . . . . . . . . . . . . . . . . . . . .

9

2.4

Derivace skalárn´ıho v´ yrazu podle vektoru . . . . . . . . . . . . . . . .

9

3 N´ ahodn´ y vektor a mnohorozmˇ ern´ e rozdˇ elen´ı

11

3.1

Sdruˇzené rozdˇelen´ı . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

3.2

Margináln´ı rozdˇelen´ı . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

3.3

Podm´ınˇené rozdˇelen´ı . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

3.4

Nezávislost veliˇcin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3.5

Charakteristiky náhodného vektoru . . . . . . . . . . . . . . . . . . . 13

3.6

V´ıcerozmˇerná rozdˇelen´ı . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.6.1

V´ıcerozmˇerné normáln´ı rozdˇelen´ı . . . . . . . . . . . . . . . . 15

4 Mnohorozmˇ ern´ a data

21

4.1

V´ ybˇerové charakteristiky . . . . . . . . . . . . . . . . . . . . . . . . . 23

4.2

Lineárn´ı transformace promˇenn´ ych . . . . . . . . . . . . . . . . . . . 24

4.3

Vzdálenost dvou objekt˚ u . . . . . . . . . . . . . . . . . . . . . . . . . 24

4.4

Chybˇej´ıc´ı hodnoty v datech . . . . . . . . . . . . . . . . . . . . . . . 25

4.5

Ovˇeˇrován´ı normality . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

4.6

Grafické metody ovˇeˇrován´ı normality . . . . . . . . . . . . . . . . . . 28

4.7

Transformace dat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

5 Line´ arn´ı regrese

31

5.1

Klasick´ y lineárn´ı model, metoda nejmenˇs´ıch ˇctverc˚ u . . . . . . . . . . 31

5.2

Odhad parametr˚ u metodou maximáln´ı vˇerohodnosti . . . . . . . . . . 33

2

OBSAH

6 Geometrie metody nejmenˇ s´ıch ˇ ctverc˚ u a regresn´ı diagnostika

37

6.1

Geometrie metody nejmenˇs´ıch ˇctverc˚ u . . . . . . . . . . . . . . . . . 37

6.2

Rozklad souˇctu ˇctverc˚ u . . . . . . . . . . . . . . . . . . . . . . . . . . 38

6.3

Regresn´ı diagnostika . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

6.4

Autokorelace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

7 Parci´ aln´ı a mnohon´ asobn´ a korelace

52

8 V´ ybˇ er regresor˚ u v mnohorozmˇ ern´ e regresi

56

8.1

Kroková regrese . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

8.2

Hledán´ı nejlepˇs´ı mnoˇziny regresor˚ u . . . . . . . . . . . . . . . . . . . 59

9 Zobecnˇ en´ı klasick´ eho line´ arn´ıho modelu

66

9.1

Transformace p˚ uvodn´ıch regresor˚ u . . . . . . . . . . . . . . . . . . . . 66

9.2

Aitken˚ uv odhad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

9.3

Heteroskedascita . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

9.4

Stochastické regresory . . . . . . . . . . . . . . . . . . . . . . . . . . 69

9.5

Diskrétn´ı regresory, umˇelé promˇenné . . . . . . . . . . . . . . . . . . 70

10 Zobecnˇ en´ y line´ arn´ı model (GLM)

73

10.1 Logistická regrese . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 11 Neline´ arn´ı regresn´ı model

82

12 Mnohorozmˇ ern´ e metody

96

12.1 Test shody vektoru stˇredn´ıch hodnot . . . . . . . . . . . . . . . . . . 97 12.2 Diskriminaˇcn´ı anal´ yza . . . . . . . . . . . . . . . . . . . . . . . . . . 99 12.3 Shluková anal´ yza . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 12.3.1 Hierarchické metody . . . . . . . . . . . . . . . . . . . . . . . 108 12.3.2 Nehierarchické metody . . . . . . . . . . . . . . . . . . . . . . 113 12.4 Anal´ yza hlavn´ıch komponent . . . . . . . . . . . . . . . . . . . . . . . 116 12.5 Faktorová anal´ yza . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 13 Literatura

127

3

1

´ Uvod

Tento text je urˇcen jako studijn´ı opora pˇredmˇetu Anal´ yza v´ıcerozmˇern´ ych dat ve vˇsech formách studia na Pˇr´ırodovˇedecké fakultˇe Ostravské university ve studijn´ım oboru Informaˇcn´ı systémy. C´ılem textu je poskytnout nezbytné základy pro statistickou anal´ yzu v´ıcerozmˇern´ ych dat, zejména pro regresn´ı anal´ yzu a vybrané mnohorozmˇerné metody. V´ yklad je zamˇeˇren sp´ıˇse na porozumˇen´ı základn´ım pojm˚ um, které je nutné pro správnou aplikaci metod v´ıcerozmˇerné anal´ yzy neˇz na matematické d˚ ukazy. Pˇresto pˇredkládán´ y text nen´ı lehké ˇcten´ı do postele pˇred span´ım. Pros´ım, poˇc´ıtejte s t´ım, ˇze budete ˇcasto nuceni usilovnˇe pˇrem´ yˇslet, vykládanou látku si postupnˇe vyjasˇ novat a k mnoha témat˚ um se opakovanˇe vracet. Nˇekdy vám m˚ uˇze pomoci i studium citovan´ ych uˇcebnic a publikac´ı nebo zdroj˚ u z internetu. Snad k pochopen´ı pomohou i obsáhlé ˇreˇsené pˇr´ıklady, které jsou pˇripojeny k vˇetˇsinˇe d˚ uleˇzit´ ych témat. Data pro ˇreˇsené pˇr´ıklady jsou pˇriloˇzena k uˇcebn´ımu textu v samostatn´ ych souborech ve formátu tabulkového procesoru Excel, aby byla ˇcitelná na vˇetˇsinˇe bˇeˇzn´ ych poˇc´ıtaˇc˚ u s r˚ uzn´ ym softwarov´ ym vybaven´ım. Kaˇzdá kapitola zaˇc´ıná pokyny pro jej´ı studium. Tato ˇca´st je vˇzdy oznaˇcena jako Pr˚ uvodce studiem s ikonou na okraji stránky. Pojmy a d˚ uleˇzité souvislosti k zapamatován´ı jsou vyznaˇceny na okraji stránky textu ikonou, kterou vid´ıte u tohoto odstavce V závˇeru kaˇzdé kapitoly je rekapitulace nejd˚ uleˇzitˇejˇs´ıch pojm˚ u. Tato rekapitulace je oznaˇcena textem Shrnut´ı a ikonou zobrazenou na okraji. Odd´ıl Kontroln´ı ot´ azky oznaˇceny ikonou jako u tohoto odstavce by vám mˇel pomoci zjistit, zda jste prostudovanou kapitolu pochopili a snad vyprovokuje i vaˇse dalˇs´ı otázky, na které budete hledat odpovˇed’. U nˇekter´ ych kapitol je pˇripomenuta Korespondeˇ cn´ı u ´ loha. Pro kombinované studium budou korespondenˇcn´ı u ´lohy zadávány v rámci kurzu daného semestru. ´ Uspˇeˇsné vyˇreˇsen´ı korespondenˇcn´ıch u ´loh je souˇcást´ı podm´ınek pro ukonˇcen´ı pˇredmˇetu v kombinovaném studiu.

5

2

Vektory a matice

Pr˚ uvodce studiem Kapitola je opakován´ım látky o vektorech a matic´ıch, coˇz byste mˇeli uˇz znát z kurz˚ u lineárn´ı algebry. Na tuto kapitolu poˇc´ıtejte se dvˇema hodinami studia, pokud jste toho moc z lineárn´ı algebry nezapomnˇeli. Jinak bude potˇrebný ˇcas delˇs´ı. Vˇsechny uvedené operace s maticemi jsou pak uˇz´ıv´ any v dalˇs´ım textu, tak je nutné, abyste je bezpeˇcnˇe zvládli.

2.1

Z´ akladn´ı pojmy

Vektory budeme oznaˇcovat mal´ ymi tuˇcn´ ymi p´ısmeny. Sloupcov´ y vektor – pˇr´ıklady:    x=  

x1 x2 .. .





    

  y=  

xp



y1 y2 .. .

    

yp

ˇ adkov´ R´ y vektor dostaneme transpozic´ı sloupcového vektoru, napˇr. xT = [x1 , · · · , xp ] Skalárn´ı souˇcin vektor˚ u je T

T

x y=y x=

p ∑

xi yi

i=1

Specilelnˇe xT x =

∑p i=1

x2i

Norma vektoru je ∥ x ∥=

√

v u p u∑ T x x=t x2 i

i=1

Kosinus smˇerového u ´hlu dvou vektor˚ u je pak cos α =

xT y ∥ x ∥∥ y ∥

Je-li cos α = 0, tj. xT y = 0, pak ˇr´ıkáme, ˇze vektory jsou ortogon´ aln´ı (jsou na sebe kolmé, cos α = 0 ). Vid´ıme, ˇze kosinus smˇerového u ´hlu vektor˚ u je vlastnˇe v´ ybˇerov´ y korelaˇcn´ı koeficient veliˇcin x, y, tedy jsou-li vektory x, y ortogonáln´ı, znamená to, ˇze veliˇciny x, y jsou nekorelované.

6

2 VEKTORY A MATICE

Matice typu (n × p) (matice budeme oznaˇcovat velk´ ymi tuˇcn´ ymi p´ısmeny) je    X=  

x11 x12 · · · x21 x22 · · · .. .. .. . . . xn1 xn2 · · ·

x1p x2p .. .

     

xnp

Matice A, B lze sˇc´ıtat (a odˇc´ıtat), pokud jsou stejného typu (n × p). A+B=C Matice C je opˇet typu (n × p) a pro jej´ı prvky plat´ı cij = aij + bij Matice A, B lze násobit, pokud jsou typu (n × p) a (p × m). AB = C Matice C je typu (n × m) a pro jej´ı prvky plat´ı cij =

p ∑

aik bkj

k=1

Jelikoˇz vektor je specieln´ı pˇr´ıpad matice maj´ıc´ı jen jeden sloupec nebo ˇra´dek, lze stejné pravidlo uˇz´ıt i pro násoben´ı matice vektorem. Pak napˇr. soustavu lineárn´ıch rovnic m˚ uˇzeme struˇcnˇe zapsat jako Ay = b. Pˇresvˇedˇcte se, ˇze opravdu je to rovnost dvou vektor˚ u, kaˇzd´ y o délce rovné poˇctu ˇra´dk˚ u matice A. Vektory jsou si rovny, kdyˇz jsou si rovny jejich vzájemnˇe si odpov´ıdaj´ıc´ı prvky, tzn. máme soustavu lineárn´ıch rovnic. Transponovaná matice XT vznikne z matice tzn.  x11 x21   x12 x22 XT =  ..  .. .  .

X tak, ˇze zamˇen´ıme ˇrádky a sloupce, ··· ··· .. .

xn1 xn2 .. .

x1p x2p · · ·

xnp

      T

Transponovaná matice XT je typu (p × n). Je zˇrejmé, ˇze plat´ı (XT ) = X. Pro transponován´ı plat´ı následuj´ıc´ı pravidla:

2.1 Základn´ı pojmy

7 (A + B)T = AT + BT (AB)T = BT AT

Hodnost matice typu (m × n) je pˇrirozené ˇc´ıslo h(C) ≤ min(m, n). Je-li matice typu (n × n), ˇr´ıkáme, ˇze je to ˇctvercov´ a matice ˇra´du n. Symetrická matice je ˇctvercová matice, pro kterou plat´ı AT = A, tzn. je symetrická podle hlavn´ı diagonály. Diagon´ aln´ı matice je ˇctvercová matice, která má vˇsechny prvky mimo hlavn´ı diagonálu rovny nule. Jednotková matice je diagonáln´ı matice s jedniˇckami na hlavn´ı diagonále. Oznaˇcujeme ji I nebo In , je-li nutno zm´ınit jej´ı rozmˇer. ∑ Stopa matice je souˇcet diagonáln´ıch prvk˚ u T r(A) = ni=1 aii . Determinant matice oznaˇcujeme |A| nebo det(A). Je to skalár (ˇc´ıselná hodnota), kterou m˚ uˇzeme chápat jako m´ıru nevyváˇzenosti matice. Kdyˇz A je typu 2 × 2, pak |A| = a11 a22 − a12 a21 . Matice A ˇra´du n je regulárn´ı, kdyˇz |A| ̸= 0. Pak existuje inverzn´ı matice A−1 , pro kterou plat´ı A−1 A = AA−1 = I Kdyˇz matice A ˇra´du n je regulárn´ı (|A| ̸= 0), pak hodnost matice je h(A) = n. Dále plat´ı, ˇze (AT )−1 = (A−1 )T (AB)−1 = B−1 A−1 Je-li matice A ˇrádu 2 regulárn´ı, pak inverzn´ı matice A−1 je [ A−1 =

a22 /∆ −a12 /∆ −a21 /∆ a11 /∆

] ,

kde ∆ = a11 a22 − a12 a21 , tj. determinant matice A. Kvadratická forma matice je skalár T

x Ax =

n ∑ n ∑

aij xi xj

i=1 j=1

Kvadratická forma je urˇcena matic´ı A. Matice B, pro kterou plat´ı bii = aii a souˇcasnˇe bij + bji = aij + aji urˇcuje tutéˇz kvadratickou formu. Existuje vˇsak jedin´ a symetrick´ a matice dané kvadratické formy.

8

2.2

2 VEKTORY A MATICE

Vlastn´ı ˇ c´ısla a vlastn´ı vektory matice

Necht’ A je ˇctvercová matice ˇrádu n. Pak vlastn´ı ˇc´ıslo (charakteristické ˇc´ıslo, eigenvalue) je takov´ y skalár λ, aby pro nenulov´ y vektor u platilo: Au = λu u je vlastn´ı (charakteristick´ y) vektor. Vid´ıme, ˇze v´ yˇse uvedenou rovnost´ı nen´ı definován jednoznaˇcnˇe, nebot’ rovnost plat´ı pro kaˇzd´ y vektor cu, c ̸= 0. Nadále budeme tedy uvaˇzovat jen vektory normované (s normou rovnou jedné), tzn. v = cu, kde c = 1/ ∥ u ∥. Rovnost pak m˚ uˇzeme pˇrepsat na tvar (A − Iλ)v = 0 Protoˇze v ̸= 0, mus´ı platit, ˇze determinant matice v závorkách na levé stranˇe rovnice |A − Iλ| = 0 Tento determinant je polynom n-tého stupnˇe, ˇreˇsen´ı je λ1 , λ2 , · · · , λn a kaˇzdému vlastn´ımu ˇc´ıslu odpov´ıdá vlastn´ı vektor vi . Kdyˇz A je symetrická matice, pak vˇsechna vlastn´ı ˇc´ısla jsou re´ aln´ a a vlastn´ı vektory jsou ortogonáln´ı, takˇze plat´ı i = 1, 2, · · · , n i ̸= j

viT vi = 1 viT vi = 0

Kdyˇz vlastn´ı vektory uspoˇrádáme do matice V = [v1 , v2 , · · · , vn ], pak V je ortogonáln´ı matice, tj. plat´ı VT = V−1 a VT V = I. Matici A m˚ uˇzeme diagonalizovat:    V AV = L =    T

λ1 0 · · · 0 λ2 · · · .. .. . . . . . 0 0 ···

0 0 .. .

     

λn

∑ Pak stopa matice A je rovna souˇctu jej´ıch vlastn´ıch ˇc´ısel, T r(A) = ni=1 λi a determinant matice A je roven souˇcinu jej´ıch vlastn´ıch ˇc´ısel, |A| = λ1 λ2 · · · λn .

2.3 Dalˇs´ı d˚ uleˇzité vlastnosti matic

9

Spektráln´ı rozklad matice A je definován jako A=

n ∑

λi vi viT

i=1

2.3

Dalˇ s´ı d˚ uleˇ zit´ e vlastnosti matic

Jestliˇze C je ortogonáln´ı matice a y = Cx, pak yT y = xT x Symetrická matice A je pozitivnˇe definitn´ı, jestliˇze kvadratická forma xT Ax > 0 pro kaˇzd´ y vektor x ̸= 0. Kdyˇz xT Ax ≥ 0, pak A je pozitivnˇe semidefinitn´ı. Pozitivnˇe definitn´ı matice má vˇsechna vlastn´ı ˇc´ısla kladná. Kdyˇz B je matice typu (n × m), s hodnost´ı m, pak BT B je pozitivnˇe definitn´ı. Jestliˇze A je pozitivnˇe definitn´ı, pak existuje regulárn´ı matice P taková, ˇze PT AP = I

a

PT P = A−1

Pseudoinverzn´ı matice A− : Kdyˇz A je matice typu (m × n), pak A− je typu (n × m) a plat´ı AA− A = A A− vˇzdy existuje, ale nen´ı jednoznaˇcnˇe urˇcena.

2.4

Derivace skal´ arn´ıho v´ yrazu podle vektoru

Je-li y = f (x1 , x2 , . . . , xn ), potom   ∂y  = ∂x  

∂y/∂x1 ∂y/∂x2 .. .

     

∂y/∂xn   ∂aT x  =  ∂x 

∂aT x/∂x1 ∂aT x/∂x2 .. .





    =    

∂aT x/∂xn Vid´ıme, ˇze ∂xT a ∂aT x = = a. ∂x ∂x

a1 a2 .. . an

   =a  

10

2 VEKTORY A MATICE

Shrnut´ı

• vektor, matice, transponován´ı vektor˚ u a matic • determinant matice, hodnost matice, inverzn´ı matice, jednotkov´ a matice, symetrická matice, diagonáln´ı matice, stopa matice • kvadratická forma, pozitivnˇe definitn´ı matice • vlastn´ı ˇc´ısla a vlastn´ı vektory matice • derivace funkce podle vektoru Kontroln´ı ot´ azky

ˇ 1. Necht’ xT = [x1 , · · · , xn ], 1 je vektor n × 1, jehoˇz prvky jsou rovny 1. Cemu jsou rovny výrazy 1T x, 1xT ? Jsou si tyto výrazy rovny? 2. Necht’ x = [x1 , x2 , , x3 ]T , a = [1, 2, 3]T , B = axT . Spoˇc´ıtejte determinant matice B. ˇ 3. Necht’ A je ˇctvercová matice ˇr´ adu n, y je vektor n×1. Cemu je rovna derivace T y Ay podle vektoru y?

11

3

N´ ahodn´ y vektor a mnohorozmˇ ern´ e rozdˇ elen´ı

Pr˚ uvodce studiem Na tuto kapitolu poˇc´ıtejte nejménˇe se tˇremi hodinami studia s t´ım, ˇze se k prob´ırané látce budete jeˇstˇe vracet po pochopen´ı dalˇs´ıch souvislost´ı. Zejména se zamˇeˇrte na pochopen´ı rozd´ıl˚ u mezi sdruˇzeným, margin´ aln´ım a podm´ınˇeným rozdˇelen´ım. Náhodn´ y vektor X = [X1 , X2 , · · · , Xp ]T je vektor, jehoˇz sloˇzky X1 , X2 , · · · , Xp jsou náhodné veliˇciny. U náhodného vektoru mus´ıme rozliˇsovat rozdˇelen´ı sdruˇzené, margin´ aln´ı a podm´ınˇené.

3.1

Sdruˇ zen´ e rozdˇ elen´ı

Pro dvousloˇzkov´ y náhodn´ y vektor (p = 2) je sdruˇzená distribuˇcn´ı funkce definována F (x1 , x2 ) = P (X1 < x1 , X2 < x2 ) Jsou-li X1 , X2 diskrétn´ı veliˇciny, pak sdruˇzená pravdˇepodobnostn´ı funkce je P (x1 , x2 ) = P (X1 = x1 , X2 = x2 ) Existuje-li nezáporná funkce f (x1 , x2 ) taková, ˇze ∫ F (x1 , x2 ) =

x1

∫

x2

f (u, v)dudv, −∞

−∞

pak náhodn´ y vektor [X1 , X2 ]T má rozdˇelen´ı spojitého typu. Funkce f (·, ·) se naz´ yvá sdruˇzená hustota. Pravdˇepodobnost, ˇze náhodné veliˇciny X1 , X2 nab´ yvaj´ı hodnot z interval˚ u [a1 , b1 ), [a2 , b2 ) je urˇcena vztahem ∫ P (a1 ≤ X1 < b1 , a2 ≤ X2 < b2 ) =

b1

∫

b2

f (x1 , x2 )dx1 dx2 a1

a2

Sdruˇzená hustota je ∂ 2 F (x1 , x2 ) f (x1 , x2 ) = ∂x1 ∂x2 Pro p > 2 plat´ı analogické vztahy, mimo jiné sdruˇzená hustota je derivac´ı distribuˇcn´ı funkce: ∂ p F (x1 , x2 , · · · , xp ) f (x) = f (x1 , x2 , · · · , xp ) = ∂x1 ∂x2 · · · ∂xp

12

3.2

´ ´ VEKTOR A MNOHOROZMERN ˇ ´ ROZDELEN ˇ Í 3 NAHODN Y E

Margin´ aln´ı rozdˇ elen´ı

Kdyˇz F (x1 , x2 ) je sdruˇzená distribuˇcn´ı funkce, pak margin´ aln´ı distribuˇcn´ı funkce veliˇcin X1 a X2 jsou F1 (x1 ) = P (X1 < x1 , X2 < ∞) = F (x1 , ∞) F2 (x2 ) = P (X1 < ∞, X2 < x2 ) = F (∞, x2 ) Pro diskrétn´ı rozdˇelen´ı margináln´ı pravdˇepodobnostn´ı funkce jsou definovány takto: ∑

P1 (x1 ) =

P (x1 , x2 )

M2

∑

P2 (x2 ) =

P (x1 , x2 )

M1

kde Mi je mnoˇzina hodnot diskrétn´ı náhodné veliˇciny Xi . Pro spojité rozdˇelen´ı margináln´ı hustoty jsou ∫ f1 (x1 ) = f (x1 , x2 )dx2 M2

∫ f2 (x2 ) =

f (x1 , x2 )dx1 M1

kde Mi je obor hodnot spojité náhodné veliˇciny Xi . Kdyˇz p > 2, pak sdruˇzené rozdˇelen´ı kaˇzdé neprázdné podmnoˇziny {X1 , X2 , · · · , Xp } je margináln´ı rozdˇelen´ı. Napˇr. pro p = 3 sdruˇzená rozdˇelen´ı náhodn´ ych vektor˚ u [X1 , X2 ]T , [X1 , X3 ]T , [X2 , X3 ]T , X1 , X2 , X3 jsou margináln´ımi rozdˇelen´ımi (máme tedy 6 margináln´ıch rozdˇelen´ı).

3.3

Podm´ınˇ en´ e rozdˇ elen´ı

V podm´ınˇeném rozdˇelen´ı jedna nebo v´ıce sloˇzek (r < p) náhodného vektoru je konstantn´ı. Pro p = 2 je podm´ınˇená distribuˇcn´ı funkce definována jako F (x1 | x2 ) = lim P (X1 < x1 | x2 < X2 ≤ x2 + ∆x2 ) ∆x2 →0

3.4 Nezávislost veliˇcin

13

Pro diskrétn´ı rozdˇelen´ı je ∑ F (x1 | x2 ) =

P (t, x2 ) P2 (x2 )

pro P2 (x2 ) ̸= 0

t<x1

a podm´ınˇená pravdˇepodobnostn´ı funkce je P (x1 | x2 ) =

P (x1 , x2 ) P2 (x2 )

pro P2 (x2 ) ̸= 0

Pro spojitá rozdˇelen´ı je podm´ınˇená distribuˇcn´ı funkce ∫ x1

f (t, x2 )dt f2 (x2 )

−∞

F (x1 | x2 ) =

pro f2 (x2 ) ̸= 0

a podm´ınˇená hustota f (x1 | x2 ) =

3.4

f (x1 , x2 ) f2 (x2 )

pro f2 (x2 ) ̸= 0

Nez´ avislost veliˇ cin

Pro nezávislé veliˇciny plat´ı F (x1 , x2 ) = F1 (x1 )F2 (x2 ) P (x1 , x2 ) = P1 (x1 )P2 (x2 ) f (x1 , x2 ) = f1 (x1 )f2 (x2 ) Jsou-li veliˇciny nezávislé, pak podm´ınˇen´ a rozdˇelen´ı jsou rovna margin´ aln´ım. Pro p > 2 plat´ı podobné vztahy pro sdruˇzenˇe (vzájemnˇe) nezávislé veliˇciny.

3.5

Charakteristiky n´ ahodn´ eho vektoru

Máme náhodn´ y vektor o délce p, X = [X1 , X2 , · · · , Xp ]T . Margináln´ı charakteristiky vektoru diskrétn´ıch náhodn´ ych veliˇcin jsou: Stˇredn´ı hodnoty E(Xj ) =

∑

xj Pj (xj ),

j = 1, 2, . . . , p

Mj

Rozptyly var(Xj ) =

∑ Mj

[xj − E(Xj )]2 Pj (xj ),

j = 1, 2, . . . , p


14

Pro vektor spojit´ ych náhodn´ ych veliˇcin jsou stˇredn´ı hodnoty ∫ E(Xj ) = xj fj (xj )dxj , j = 1, 2, . . . , p Mj

a rozptyly ∫ [xj − E(Xj )]2 fj (xj )dxj ,

var(Xj ) =

j = 1, 2, . . . , p

Mj

Podm´ınˇené charakteristiky vektoru diskrétn´ıch náhodn´ ych veliˇcin pro p = 2 jsou definovány následovnˇe. Stˇredn´ı hodnoty

∑

E(X1 | x2 ) =

x1 P1 (x1 | x2 )

M1

Rozptyly

∑

var(X1 | x2 ) =

[x1 − E(X1 | x2 )]2 P (x1 | x2 )

M1

Podm´ınˇené charakteristiky vektoru spojit´ ych náhodn´ ych veliˇcin pak jsou: Stˇredn´ı hodnoty

∫ E(X1 | x2 ) =

x1 f (x1 | x2 )dx1 M1

Rozptyly

∫ var(X1 | x2 ) =

[x1 − E(X1 | x2 )]2 f (x1 | x2 )dx1 M1

Podm´ınˇená stˇredn´ı hodnota E(X1 | x2 ) se naz´ yvá regresn´ı funkce (závislost X1 na X2). Podm´ınˇen´ y rozptyl var(X1 | x2 ) se naz´ yvá skedastická funkce. Je-li podm´ınˇen´ y rozptyl var(X1 | x2 ) konstantn´ı pro vˇsechna x2 , pak o rozdˇelen´ı ˇr´ıkáme, ˇze je homoskedastické, nen´ı-li konstantn´ı, mluv´ıme o heteroskedastickém rozdˇelen´ı náhodného vektoru [X1 , X2 ]T .

3.6

V´ıcerozmˇ ern´ a rozdˇ elen´ı

Základn´ımi charakteristikami v´ıcerozmˇerného rozdˇelen´ı je vektor stˇredn´ıch hodnot    E(X) =   

E(X1 ) E(X2 ) .. . E(Xp )

     

3.6 V´ıcerozmˇerná rozdˇelen´ı

15

a kovarianˇcn´ı (varianˇcn´ı) matice [ ] Σ = var(X) = cov(X) = E (X − E(X))(X − E(X)T )

(1)

coˇz znamená, ˇze    Σ=  

σ12 σ12 · · · σ21 σ22 · · · .. .. ... . . σp1 σp2 · · ·

σ1p σ2p .. .

   ,  

σp2

kde σij je kovariance dvou náhodn´ ych veliˇcin, tj. σij = cov(Xi , Xj ) = E [(Xi − EXi )(Xj − EXj )] a σii = σi2 je rozptyl var(Xi ). Vid´ıme, ˇze kovarianˇcn´ı matice Σ je symetrická, nebot’ σij = σji .

3.6.1

V´ıcerozmˇ ern´ e norm´ aln´ı rozdˇ elen´ı

Náhodn´ y vektor X má v´ıcerozmˇerné normáln´ı rozdˇelen´ı, jestliˇze jeho hustota je dána vztahem ( ) (x − µ)T Σ−1 (x − µ) −1/2 p/2 f (x) = (2π) |Σ| exp − , (2) 2 kde µ je vektor stˇredn´ıch hodnot a Σ je kovarianˇcn´ı matice. V´ıcerozmˇerné normáln´ı rozdˇelen´ı má tyto vlastnosti: • lineárn´ı kombinace prvk˚ u z X maj´ı normáln´ı rozdˇelen´ı • vˇsechny podmnoˇziny X maj´ı normáln´ı rozdˇelen´ı • nekorelovanost veliˇcin z X (sloˇzek vektoru X) znamená i jejich nezávislost • vˇsechna podm´ınˇená rozdˇelen´ı jsou normáln´ı Pro jednorozmˇerné normáln´ı rozdˇelen´ı z rov. (2) dostaneme 1

f (x) = √ 2πσ 2

( ) (x − µ)2 exp − 2σ 2

V exponentu je ˇctverec vzdálenosti ( 2

u =

x−µ σ

)2 ,

16


tedy vzdálenosti x od stˇredn´ı hodnoty µ, kde jednotkou vzdálenosti je σ. I pro v´ıcerozmˇerné normáln´ı rozdˇelen´ı je moˇzno chápat kvadratickou formu v exponentu jako ˇctverec vzdálenosti vektoru x od vektoru µ, ve kterém je obsaˇzena i informace z kovarinaˇcn´ı matice C 2 = (x − µ)T Σ−1 (x − µ). C je Mahalanobisova vzdálenost, pro zvolenou hodnotu je jej´ı ˇctverec geometricky √ plocha elipsoidu se stˇredem µ a osami c λj vj pro j = 1, 2, · · · , p, kde λj jsou vlastn´ı ˇc´ısla matice Σ a vj jsou vlastn´ı vektory matice Σ. Pˇ r´ıklad 3.1 Pro dvourozmˇerné normáln´ı rozdˇelen´ı s parametry EX1 = µ1 , EX2 = µ2 , varX1 = σ12 , varX2 = σ22 a kovarianc´ı σ12 je korelaˇcn´ı matice [ Σ=

σ12 σ12 σ12 σ22

]

[ =

σ12 σ1 σ2 ρ σ1 σ2 ρ σ22

]

nebot’ korelaˇcn´ı koeficient ρ = σ12 /(σ1 σ2 ). Determinant kovarianˇcn´ı matice je pak 2 |Σ| = σ12 σ22 − σ12 = σ12 σ22 (1 − ρ2 ).

Vid´ıme, ˇze tento determinant je roven nule, kdyˇz ρ2 = 1.

Podm´ınˇené rozdˇelen´ı X1 |x2 je normáln´ı se stˇredn´ı hodnotou β0 + β1 x2 a rozptylem σ12 (1 − ρ2 ) σ12 β0 = µ1 − β1 µ2 β1 = 2 σ2 Podm´ınˇená stˇredn´ı hodnota E(X1 | x2 ) závis´ı lineárnˇe na x2 . Rozptyl veliˇciny X1 nezávis´ı na x2 . Pro dvourozmˇerné normáln´ı rozdˇelen´ı m˚ uˇzeme elipsy konstantn´ı hustoty (f (x1 , x2 ) = const) znázornit graficky. Kdyˇz veliˇciny X1 , X2 jsou nekorelované, tj. v pˇr´ıpadˇe v´ıcerozmˇerného normáln´ıho rozdˇelen´ı i nezávislé, osy elipsy konstantn´ı hustoty jsou rovnobˇeˇzné s x1 , x2 , jinak jsou pootoˇceny. Názornˇe to ukazuj´ı následuj´ıc´ı obrázky.


17

Obrázek 1: Hustota dvourozmˇerného normáln´ıho rozdˇelen´ı a elipsy konstantn´ı hustoty, µ1 = µ2 = 0, σ1 = σ2 = 1, ρ = 0

Obrázek 2: Hustota dvourozmˇerného normáln´ıho rozdˇelen´ı a elipsy konstantn´ı hustoty, µ1 = µ2 = 0, σ1 = 1, σ2 = 2, ρ = 0

18


Obrázek 3: Hustota dvourozmˇerného normáln´ıho rozdˇelen´ı a elipsy konstantn´ı hustoty, µ1 = µ2 = 0, σ1 = σ2 = 1, ρ = −0.6

Obrázek 4: Hustota dvourozmˇerného normáln´ıho rozdˇelen´ı a elipsy konstantn´ı hustoty, µ1 = µ2 = 0, σ1 = σ2 = 1, ρ = 0.8


19

Obrázek 5: Hustota dvourozmˇerného normáln´ıho rozdˇelen´ı a elipsy konstantn´ı hustoty, µ1 = µ2 = 0, σ1 = 1, σ2 = 2, ρ = 0.8

Obrázek 6: Hustota dvourozmˇerného normáln´ıho rozdˇelen´ı a elipsy konstantn´ı hustoty, µ1 = µ2 = 0, σ1 = 1, σ2 = 2, ρ = −0.8


20

Shrnut´ı

• n´ ahodný vektor, sdruˇzené, margin´ aln´ı a podm´ınˇené rozdˇelen´ı • nezávislost veliˇcin • vektor stˇredn´ıch hodnot, kovarianˇcn´ı (varianˇcn´ı) matice • v´ıcerozmˇerné normáln´ı rozdˇelen´ı Kontroln´ı ot´ azky

1. Sdruˇzená pravdˇepodobnostn´ı funkce náhodného vektoru [X, Y ]T je zadána tabulkou X = x1 X = x2 X = x3 Y = y1 p11 p12 p13 Y = y2 p21 p22 p23 Kolik margináln´ıch a kolik podm´ınˇených rozdˇelen´ı má tento vektor? Vyj´ adˇrete vˇsechny margináln´ı pravdˇepodobnostn´ı funkce a jednu zvolenou podm´ınˇenou pravdˇepodobnostn´ı funkci. 2. Necht’ A je ˇctvercová matice ˇr´ adu n, jej´ıˇz prvky jsou re´ aln´ a ˇc´ısla (konstanty), x je náhodný vektor typu n × 1, y = Ax. Vyj´ adˇrete vektor stˇredn´ıch hodnot n´ ahodného vektoru y a kovarianˇcn´ı matici náhodného vektoru y. Korespondeˇ cn´ı u ´ loha Korespondenˇcn´ı u ´lohy budou zadáv´ any ke kaˇzdému kursu samostatnˇe.

21

4

Mnohorozmˇ ern´ a data

Pr˚ uvodce studiem Tato kapitola má poslouˇzit pro orientaci v problematice statistické analýzy v´ıcerozmˇerných dat. Jsou zde uvedeny d˚ uleˇzité výbˇerové charakteristiky, které lze z dat vyhodnotit. Také jsou zm´ınˇeny techniky ovˇeˇrov´ an´ı pˇredpoklad˚ u o rozdˇelen´ı, zejména o norm´ aln´ım rozdˇelen´ı a nˇekteré transformace, které lze uˇz´ıt v analýze dat. Poˇc´ıtejte se tˇremi hodinami studia s t´ım, ˇze se k prob´ırané látce budete jeˇstˇe podle potˇreby vracet. Prozat´ım jsme se zab´ yvali otázkami abstraktn´ıho popisu vztah˚ u mezi náhodn´ ymi veliˇcinami, pˇredevˇs´ım náhodn´ ym vektorem. Nyn´ı obrát´ıme pozornost k praktiˇctˇejˇs´ım problém˚ um mnohorozmˇern´ ych dat. Pˇripomeˇ nme, ˇze veliˇciny, které zjiˇst’ujeme na sledovan´ ych objektech, jsou r˚ uzného typu. Podle oboru jejich hodnot rozliˇsujeme: • veliˇciny spojité - mohou nab´ yvat nespoˇcetnˇe mnoho hodnot, napˇr. ˇcas, délka atd. • veliˇciny nespojité (diskrétn´ı) - nab´ yvaj´ı jen spoˇcetnˇe mnoho hodnot, v praxi jen koneˇcného poˇctu a ˇcasto nˇekolika málo moˇzn´ ych hodnot, napˇr. kategoriáln´ı veliˇciny, vyjadˇruj´ıc´ı pˇr´ısluˇsnost k nˇejaké skupinˇe (kategorii) objekt˚ u • alternativn´ı (dichotomické, binárn´ı) veliˇciny, patˇr´ı mezi nespojité, ale t´ım, ˇze mohou nab´ yvat jen dvou moˇzn´ ych hodnot, ˇcasto interpretovan´ ych jako ANO/NE, TRUE/FALSE nebo 1/0, b´ yvá nˇekdy uˇziteˇcné nahl´ıˇzet na nˇe jako na zvláˇstn´ı typ veliˇcin Dále veliˇciny m˚ uˇzeme rozliˇsovat podle ˇskály, ve které mˇeˇr´ıme: • nomináln´ı (kategoriáln´ı) • ordináln´ı (poˇradové) • rozd´ılové (intervalové) kvantitativn´ı (metrické, je definována vzdálenost dvou hodnot) • pomˇerové kvantitativn´ı (metrické, je definována vzdálenost dvou hodnot) Kromˇe tˇechto hledisek, která klasifikuj´ı veliˇciny, je d˚ uleˇzité m´ıt i na pamˇeti, jak vlastnˇe data vznikla, co zobrazuj´ı a jaké jsou vztahy mezi pozorovan´ ymi objekty. Podstatné je, zda objekty m˚ uˇzeme povaˇzovat za nezávislé nebo zda vznikla jako ˇrada pozorován´ı téhoˇz objektu v r˚ uzn´ ych obdob´ıch. R˚ uzné situace rozliˇsuje následuj´ıc´ı tabulka.

ˇ ´ DATA 4 MNOHOROZMERN A

22

´ Ulohy ˇ reˇ sen´ e anal´ yzou dat, ˇ casov´ y prvek v datech: poˇcet poˇcet poˇcet typ objekt˚ u veliˇcin obdob´ı u ´lohy 1 p 1 kasuistika, pˇr´ıpadová studie n 1 1 jednorozmˇerná anal´ yza 1 1 T jednorozmˇerná ˇcasová ˇrada n 1 T T = 2 párové srovnán´ı, T > 2 opakovaná mˇeˇren´ı n p 1 v´ıcerozmˇerná anal´ yza dat 1 p T v´ıcerozmˇerná ˇcasová ˇrada n p T longitudináln´ı studie Poznámka: n, p, T > 1

Anal´ yza v´ıcerozmˇern´ ych dat se vˇetˇsinou zab´ yvá daty, kdy máme p veliˇcin pozorovan´ ych na n objektech. Rozliˇsit m˚ uˇzeme následuj´ıc´ı situace:

• vˇsechny veliˇciny metrické

– n bod˚ u v Rp

• vˇsechny veliˇciny kategoriáln´ı

– mnohorozmˇerné kontingenˇcn´ı tabulky

• sm´ıˇsená data – datová matice se rozdˇel´ı na podsoubory – analogie s dvou/v´ıcev´ ybˇerov´ ymi u ´lohami

Pro charakterizaci v´ıcerozmˇern´ ych dat potˇrebujeme odhadnout charakteristiky prozmˇerného vektoru náhodn´ ych veliˇcin, tj. vektor stˇredn´ıch hodnot a kovarianˇcn´ı matici (je symetrická), tedy urˇcit následuj´ıc´ı poˇcet v´ ybˇerov´ ych charakteristik: 2p +

p(p − 1) p2 − p p2 + 3p = 2p + = 2 2 2

Poˇcet odhadovan´ ych parametr˚ u tedy roste kvadraticky s poˇctem veliˇcin, napˇr. pro p = 10 je to 65 charakteristik, pro p = 40 je to uˇz 860 charakteristik. Pokud jsou veliˇciny kategoriáln´ı, potˇrebujeme odhadovat i sdruˇzenou pravdˇepodobnostn´ı funkci. Poˇcet pol´ıˇcek v tabulce roste exponenciálnˇe s poˇctem veliˇcin. Tud´ıˇz pokud maj´ı b´ yt tyto odhady d˚ uvˇeryhodné, potˇrebujeme, aby v´ıcerozmˇerná data byla mˇela dostateˇcn´ y rozsah, tzn. n bylo velké.

4.1 V´ ybˇerové charakteristiky

4.1

23

V´ ybˇ erov´ e charakteristiky

V´ıcerozmˇerná dat jsou reprezentována datovou matic´ı (n × p)    X=  

x11 x12 · · · x21 x22 · · · .. .. .. . . . xn1 xn2 · · ·



x1p x2p .. .



    =    

xT1 xT2 .. .

     

xTn

xnp

tzn., ˇze i−t´ y ˇrádek datové matice je ˇra´dkov´ y vektor pozorován´ı p veliˇcin na i−tém objektu. Vektor pr˚ umˇer˚ u jednotliv´ ych veliˇcin 1∑ xi x= n i=1 n

tj.

1∑ xj = xij , n i=1 n

j = 1, 2, · · · , p

Wishartova matice typu (p × p) má prvky wj,j ′ =

n ∑

j, j ′ = 1, 2, · · · , p

(xij − x¯j )(xij ′ − x¯j ′ ),

i=1

a m˚ uˇzeme ji zapsat také jako W=

n ∑

(xi − x)(xi − x)T

i=1

V´ ybˇerová kovarianˇcn´ı matice je opˇet typu p × p S=

1 W, n−1

jej´ı prvky jsou v´ ybˇerové kovariance tj. 1 ∑ = (xij − x¯j )(xij ′ − x¯j ′ ) n − 1 i=1 n

sjj ′

Korelaˇcn´ı matice má prvky rjj ′ = sjj ′ /(sj sj ′ ) tj. v´ ybˇerové korelaˇcn´ı koeficienty a má tvar   1 r12 · · · r1p    r21 1 · · · r2p   . .. . . .   . . ..  .   . rp1 rp2 · · ·

1


24

4.2

Line´ arn´ı transformace promˇ enn´ ych

Pˇri anal´ yze v´ıcerozmˇern´ ych dat je ˇcasto v´ yhodné pracovat s odvozen´ ymi veliˇcinami, které vzniknou z p˚ uvodn´ıch lineárn´ı transformac´ı, napˇr. s centrovan´ ymi promˇenn´ ymi vij = xij − x¯j , pro které vektor pr˚ umˇer˚ u je nulov´ y, v ¯ = 0, a kovarianˇcn´ı a korelaˇcn´ı matice se nezmˇen´ı, tzn. Sv = Sx , Rv = Rx . Dalˇs´ı ˇcasto uˇz´ıvanou transformac´ı je normován´ı. Normovaná hodnoty promˇenn´ ych vzniknou vycentrován´ım a vydˇelen´ım smˇerodatnou odchylkou p˚ uvodn´ıch promˇenn´ ych, tj. xij − x¯j i = 1, 2, · · · , n j = 1, 2, · · · , p. zij = sj Pak vektor pr˚ umˇer˚ u je nulov´ y z = 0, vˇsechny rozptyly a smˇerodatné odchylky normovan´ ych promˇenn´ ych jsou rovny jedné a kovarianˇcn´ı i korelaˇcn´ı matice normovan´ ych promˇenn´ ych jsou si rovny, tj. Sz = Rz = Rx a jsou rovny korelaˇcn´ı matici p˚ uvodn´ıch netrasformovan´ ych promˇenn´ ych. Pro veliˇcinu, která vznikne lineárn´ı kombinac´ı p˚ uvodn´ıch promˇenn´ ych T

ui = c xi =

p ∑

cj xij

j=1

plat´ı ¯, u¯ = cT x

4.3

s2u = cT Sx c.

Vzd´ alenost dvou objekt˚ u

ˇ adek datové matice, tj. vektor xT m˚ R´ uˇzeme povaˇzovat za souˇradnice bodu v prozmˇerném prostoru. Pak je uˇziteˇcné zab´ yvat se vzdálenost´ı dvou objekt˚ u. Jedno′ rozmˇerná vzdálenost (kdyˇz p = 1) dvou objekt˚ u i, i je absolutn´ı hodnota rozd´ılu pozorovan´ ych hodnot, d(i, i′ ) = |xi − x′i |. Pro v´ıcerozmˇerná data m˚ uˇzeme definovat r˚ uzné vzdálenosti. Eukleidovská vzdále′ nost dvou objekt˚ u i, i je v v u∑ u∑ p u u p 2 ′ 2 t DE (i, i ) = (xij − xi′ j ) = t dj (i, i′ ) j=1

j=1

4.4 Chybˇej´ıc´ı hodnoty v datech

25

Normovaná vzdálenost v v u∑ u∑ p u u p d2j (i, i′ ) ′ 2 t DN (i, i ) = (zij − zi′ j ) = t s2j j=1 j=1 Mahalanobisova vzdálenost respektuje jak rozd´ılnost ve variabilitˇe, tak korelaˇcn´ı strukturu. Jej´ı ˇctverec je pak 2 DM (i, i′ ) = dT S−1 d = (xi − xi′ )T S−1 (xi − xi′ ),

kde d = xi − xi′ . Pokud S = σ 2 I (vˇsechny rozptyly jsou shodné, veliˇciny nekorelované), pak DN = DM . Pro vyhledáván´ı odlehl´ ych pozorován´ı je uˇziteˇcná v´ ybˇerová Mahalanobisova vzdálenost od tˇeˇziˇstˇe naˇsich pozorován´ı (od vektoru pr˚ umˇer˚ u) ¯ i′ )T S−1 (xi − x ¯ i′ ). Di2 = (xi − x Pro p-rozmˇerné normáln´ı rozdˇelen´ı populace je (n − p)n 2 D ∼ F (p, n − p), (n2 − 1)p i podle toho tedy m˚ uˇzeme posudit, zda je pozorován´ı odlehlé.

4.4

Chybˇ ej´ıc´ı hodnoty v datech

Zpracován´ı mnohorozmˇern´ ych dat v reáln´ ych u ´lohách je nˇekdy komplikováno t´ım, ˇze data nejsou u ´plná, hodnoty nˇekter´ ych prvk˚ u datové nejsou k dispozici (slangovˇe oznaˇcovány jako missingy, missings). Obvykl´ y jednoduch´ y postup je vypustit veliˇciny s mnoha missingy a vypustit pˇr´ıpady (objekty) s mnoha missingy a u ´sudku o tˇechto veliˇcinách a objektech se zˇr´ıci. Nejbˇeˇznˇejˇs´ı postup uˇz´ıvan´ y ve vˇetˇsinˇe statistick´ ych paket˚ u je automatické vypouˇstˇen´ı pˇr´ıpad˚ u (objekt˚ u) s jedn´ım nebo v´ıce missingy, tzv. CASEWISE strategie. Tato strategie je z hlediska dalˇs´ıho zpracován´ı nejbezpeˇcnˇejˇs´ı, ale m˚ uˇze nˇekdy vést k pˇriliˇs velké ztrátˇe informace, kdy mnoho objekt˚ u je vyˇrazeno jen kv˚ uli jedné chybˇej´ıc´ı hodnotˇe. Ale pˇri této strategii v´ ybˇerová kovarianˇcn´ı matice z˚ ustane pozitivnˇe definitivn´ı T (kdyˇz hodnost (X) = p), nebot’ hodnost (X X) = p. Pro odhad kovarianˇcn´ı matice lze uˇz´ıt i strategii PAIRWISE, kdy kaˇzdá kovariance se poˇc´ıtá ze vˇsech moˇzn´ ych dvojic a pr˚ umˇery v n´ı jen z hodnot uˇzit´ ych pro v´ ypo-


26

ˇcet kovariance nebo strategii ALLVALUE, kdy pro pr˚ umˇery se uˇzij´ı vˇsechny moˇzné (dostupné) hodnoty. Pˇri tomto postupu se vyuˇzije dat d˚ ukladnˇeji, ale v´ ybˇerová kovarianˇcn´ı matice nemus´ı b´ yt pozitivnˇe definitivn´ı. Jin´ y pˇr´ıstup k práci s missingy je tak zvaná imputace, ˇcili doplnˇen´ı chybˇej´ıc´ıch hodnot nˇejak´ ymi vhodn´ ymi, obvykle náhodn´ ymi hodnotami z rozdˇelen´ı, které je shodné nebo podobné s rozdˇelen´ım pozorovan´ ych hodnot v datové matici. C´ılem imputace je zabránit ztrátˇe informace nevyuˇzit´ım vˇsech pozorovan´ ych hodnot v datové matici za cenu rizika, ˇze informaci obsaˇzenou v datech trochu zkresl´ıme. Bˇeˇzn´ ymi metodami imputace, které jsou implementovány ve standardn´ım statistickém software, jsou následuj´ıc´ı postupy doplnˇen´ı chybˇej´ıc´ıch hodnot: • pr˚ umˇerem • náhodnˇe z pˇredpokládaného rozdˇelen´ı s vyuˇzit´ım odhadu jeho parametr˚ u, nej2 ˇcastˇeji je chybˇej´ıc´ı prvek xij nahrazen hodnotami z N (µ, σ ), kde hodnoty parametr˚ u odhadneme z dostupn´ ych dat, µ ˆ = x¯j , σ b2 = s2j • regresn´ım modelem, jehoˇz parametry odhadneme ze zb´ yvaj´ıc´ıch (n − 1) objekt˚ u, chybˇej´ıc´ı hodnota se nahrad´ı hodnotou predikovanou modelem, pˇr´ıpadnˇe jeˇstˇe modifikovanou náhodn´ ym kol´ısán´ım.

4.5

Ovˇ eˇ rov´ an´ı normality

Mnoho metod v´ıcerozmˇerné anal´ yzy dat vycház´ı z pˇrepokladu, ˇze náhodná sloˇzka v datech má normáln´ı rozdˇelen´ı. Nˇekdy je vyˇzadována v´ıcerozmˇerná normalita, tj. p-rozmˇerné normáln´ı rozdˇelen´ı, nˇekdy staˇc´ı jen normalita nˇekter´ ych veliˇcin. Jak v´ıme, p-rozmˇerné sdruˇzené normáln´ı rozdˇelen´ı má margináln´ı rozdˇelen´ı normáln´ı, avˇsak neplat´ı to naopak, tzn. margináln´ı normalita nezaruˇcuje sdruˇzenou normalitu. Uvedeme struˇcnˇe nˇekteré metody, kter´ ymi lze normalitu (vˇetˇsinou jen margináln´ı) testovat. Jednorozmˇ ern´ a normalita Testy dobré shody, ve kter´ ych se empirické rozdˇelen´ı porovnává s normáln´ım rozdˇelen´ım. Rozpˇet´ı pozorovan´ ych hodnot se rozdˇel´ı na r interval˚ u a porovnaj´ı se ˇcetnosti pozorovan´ ych hodnot v jednotliv´ ych intervalech s teoretick´ ymi ˇcetnostmi, které bychom oˇcekávali pˇri v´ ybˇeru stejného rozsahu z normáln´ıho rozdˇelen´ı. Hranice interval˚ u se vol´ı • bud’ ekvidistantnˇe (intervaly jsou stejnˇe ˇsiroké) tak, aby teoretické (oˇcekávané) ˇcetnosti Ψi > 1 byly pro vˇsechny intervaly a Ψi > 5 pro 80% hodnot (tzv. Cochranovo pravidlo). • nebo hranice r interval˚ u se vol´ı tak aby teoretické ˇcetnosti Ψi byly konstantn´ı, ” nejˇcastˇeji se vol´ı Ψ1 = . . . = Ψr = 5. Pokud se rozhodneme pro tento zp˚ usob

4.5 Ovˇeˇrován´ı normality

27

volby hranic interval˚ u, volbou poˇzadované hodnoty teoretické ˇcetnosti je urˇcen pˇri daném rozsahu v´ ybˇeru i poˇcet interval˚ u r.

Testová statistika je pak

r ∑ (ni − Ψi )2

∼ χ2(r−1)

Ψ2i

i=1

Kolmogorov˚ uv test – porovnává se v´ ybˇerová (Fn a teoretická (F ) distribuˇcn´ı funkce. V´ ybˇerová distribuˇcn´ı funkce definována jako Fn (0) = 0,

Fn (i) =

i n

a testovou statistikou je maximum absolutn´ı hodnoty rozd´ılu porovnávan´ ych distribuˇcn´ıch funkc´ı k2 = max(|F − Fn |). Pro malé rozsahy v´ ybˇeru jsou kritické hodnoty této statistiky tabelovány, pro n > 50 se m˚ uˇze uˇz´ıt asymptotická aproximace . 1, 36 k2 (0, 95) = √ n

. 1, 63 k2 (0, 99) = √ n

Testy ˇsikmosti a ˇspiˇcatosti Pˇri v´ ypoˇctu ˇsikomosti a ˇspiˇcatosti se uˇz´ıvaj´ı centráln´ı v´ ybˇerové momenty Mk . 1∑ (xi − x¯)k . n i=1 n

Mk =

ˇ Sikmost je definována jako g1 =

M3 3/2

M2

a normáln´ı rozdˇelen´ı má ˇsikmost rovnou nule (je symetrické). ˇ catost je rovna Spiˇ g2 =

M4 −3 M22

a i ta je pro normáln´ı rozdˇelen´ı nulová, nebot’ u normáln´ıho rozdˇelen´ı je pomˇer M4 /M22 = 3. K testován´ı normality lze uˇz´ıt statistiky √ K (3) =

g12 (n + 1)(n + 3) 6(n − 2)


28 √

a K

(4)

=

(n + 1)2 (n + 3)(n + 5) 24n(n − 2)(n − 3)

( g2 +

6 n+1

) ,

které obˇe maji pˇribliˇznˇe normované normáln´ı rozdˇelen´ı N (0, 1).

4.6

Grafick´ e metody ovˇ eˇ rov´ an´ı normality

Tyto grafické metody umoˇzn ˇuj´ı rychlé vizuáln´ı posouzen´ı shody empirického rozdˇelen´ı s normáln´ım (pˇr´ıpadnˇe i jin´ ym) rozdˇelen´ım a proto jsou velmi ˇcasto vyuˇz´ıvány a také jsou implementovány v bˇeˇznˇe uˇz´ıvaném statistickém software. Kromˇe histogram˚ u, do kter´ ych se proloˇz´ı i teoretické rozdˇelen´ı a graf˚ u porovnávaj´ıc´ıch empirickou distribuˇcn´ı funkci s teoretickou se uˇz´ıvá tzv. QQ-graf, kvantilov´ y graf. QQ je zkratka pro kvantil (angl. quantile). Pro sestrojen´ı kvantilového grafu nejdˇr´ıve uspoˇrádáme v´ ybˇer, tj. x(1) ≤ x(2) ≤ . . . ≤ x(n) Hodnoty v´ ybˇerové distribuˇcn´ı funkce se spoˇctou jako V DF (x(i) ) =

i− n

1 2

nebo V DF (x(i) ) =

i n+1

a kvantily x(i) se vynesou do grafu proti odpov´ıdaj´ıc´ım kvantil˚ um normovaného normáln´ıho rozdˇelen´ı, (tj.napˇr. proti hodnotám u(i/(n + 1)), kde u(p) je p-kvantil normovaného normáln´ıho rozdˇelen´ı). Pokud je v´ ybˇerové rozdˇelen´ı normáln´ı, grafem je pˇribliˇznˇe pˇr´ımka.

4.7

Transformace dat

Pokud zjist´ıme, ˇze namˇeˇrená data nejsou z normáln´ıho rozdˇelen´ı, je nˇekdy uˇziteˇcné pouˇz´ıt vhodnou transformaci, aby transformac´ı p˚ uvodn´ı veliˇciny vznikla odvozená veliˇcina, která normáln´ı rozdˇelen´ı má. Potom lze aplikovat metody vyˇzaduj´ıc´ı normáln´ı rozdˇelen´ı na transformované veliˇciny. Transformac´ı rozum´ıme takov´ y pˇrepoˇcet, yi = f (xi ), aby se yi , i = 1, . . . , n pˇribl´ıˇzilo v´ ybˇeru z normáln´ımu rozdˇelen´ı. Ze zkuˇsenosti lze doporuˇcit tyto transformace: • odmocninová transformace y =

√ x, kdyˇz x jsou ˇcetnosti

x ), kdyˇz x jsou relativn´ı ˇcetnosti • logitová transformace y = ln( 1−x

4.7 Transformace dat

29

• logaritmická transformace y = ln x pokud mˇeˇrená veliˇcina má log-normáln´ı rozdˇelen´ı - napˇr. v´ ydaje na domácnost, náklady na v´ yrobek atd. Jinou moˇznost´ı je transformace odvozen´ a z dat, kterou navrhli Box a Cox v roce 1964. Tato transformace poskytne hodnoty odvozené veliˇciny y , které se nejv´ıce pˇribliˇzuj´ı normáln´ımu rozdˇelen´ı. { y=

x2 −1 λ

ln x

pro λ ̸= 0 pro λ = 0

λ se odhadne jako λ = λmax , které maximalizuje vˇerohodnostn´ı funkci ] [ n n ∑ n 1∑ 2 ln L(λ) = − ln (yi − y¯) + (λ − 1) ln xi 2 n i=1 i=1 Asymptotick´ y interval 100(1 − α)% – n´ı spolehlivosti pro λ: 2 [ln L(λmax ) − ln L(λ)] ≤ χ21−α (1), ˇcili v tomto intervalu jsou vˇsechna x, pro která plat´ı: 1 ln L(x) ≥ ln L(λmax ) − χ21−α (1) 2 Charakteristiky pro data, kter´ a nejsou z norm´ aln´ıho rozdˇ elen´ı Takov´ ymi charakteristikami jsou ty, jejichˇz hodnoty nejsou ovlivnˇeny odhlehl´ ymi hodnotami v datech. Jako pˇr´ıklad uvedeme uˇrezávan´ y pr˚ umˇer (trimmed mean): n−m ∑ 1 x¯(α) = x(i) n − 2m i=m+1

m = int

( αn )

α je % uˇr´ıznut´ ych poˇra´dkov´ ych statistik na kaˇzdém konci. Podobnˇe lze zavést i uˇrezávan´ y odhad rozptylu atd.

100


30

Shrnut´ı

• vektor výbˇerových pr˚ umˇer˚ u, výbˇerov´ a kovarianˇcn´ı matice, výbˇerov´ a korelaˇcn´ı matice • ovˇeˇrován´ı normality, QQ graf Kontroln´ı ot´ azky

1. Vygenerujte si (napˇr. v Excelu) náhodný výbˇer z rovnomˇerného rozdˇelen´ı o rozsahu 100 a zkonstruujte QQ graf 2. Vygenerujte náhodný výbˇer stejného rozsahu z norm´ aln´ıho rozdˇelen´ı, zkonstruujte QQ graf a porovnejte s QQ grafem z pˇredchoz´ı otázky

31

5

Line´ arn´ı regrese

Pr˚ uvodce studiem Tato kapitola je pro pochopen´ı moˇznost´ı a technik analýzy v´ıcerozmˇerných dat kl´ıˇcov´ a. Proto na tuto kapitolu poˇc´ıtejte nejménˇe se tˇremi hodinami usilovného studia s t´ım, ˇze se k prob´ırané látce budete jeˇstˇe vracet po pochopen´ı dalˇs´ıch souvislost´ı. Regrese je jednou z velice ˇcasto aplikovan´ ych statistick´ ych metod, uvád´ı se, ˇze dokonce naprostá vˇetˇsina aplikac´ı (70 − 90%) je nˇejakou formou regresn´ıch metod. Poˇcátky metody nejmenˇs´ıch ˇctverc˚ u jsou dokumentovány jiˇz na zaˇca´tku 19. stolet´ı (Legendre, Gauss), minimalizace souˇctu absolutn´ıch hodnot odchylek je pˇripisována Galileovi jeˇstˇe o pár des´ıtek let dˇr´ıve.

5.1

Klasick´ y line´ arn´ı model, metoda nejmenˇ s´ıch ˇ ctverc˚ u

Klasick´ y model lineárn´ı regrese lze zapsat jako yi = β0 + β1 xi1 + · · · + βk xik + εi

(3)

kde yi je pozorovaná hodnota náhodné veliˇciny Y xi1 , · · · , xik jsou hodnoty vysvˇetluj´ıc´ıch promˇenn´ ych (regresor˚ u, prediktor˚ u) β0 , β1 , · · · , βk jsou parametry modelu (fixn´ı, leˇc neznámé hodnoty) εi je náhodná sloˇzka i = 1, 2, · · · , n je index pozorován´ı (objektu) Rovnici (3) m˚ uˇzeme zapsat maticovˇe y = Xβ + ε

(4)

yi = xTi β + εi

(5)

[pro i-t´ y ˇra´dek pak

Vektory jsou oznaˇceny pˇr´ısluˇsn´ ymi mal´ ymi tuˇcn´ ymi p´ısmeny, matice velk´ ymi tuˇcn´ ymi p´ısmeny. Matice X má k + 1 sloupc˚ u, v prvn´ım sloupci jsou jedniˇcky, dalˇs´ıch sloupc´ıch jsou hodnoty vysvˇetluj´ıc´ıch veliˇcin. Obvykl´ ymi pˇredpoklady v klasickém lineárn´ı modelu jsou: 1. E(ε) = 0, tj. vektor stˇredn´ıch hodnot náhodné sloˇzky je roven nulovému vektoru

´ Í REGRESE 5 LINEARN

32

2. cov(ε) = σ 2 In , σ 2 > 0, I je jednotková matice ˇra´du n, tj. náhodné sloˇzky jsou nekorelované a jejich rozptyl je konstantn´ı 3. X je nenáhodná matice typu n × (k + 1) 4. hodnost matice X, h(X) = k + 1 ≤ n, tj. sloupce matice X nejsou lineárnˇe závislé a poˇcet pozorován´ı je alespoˇ n roven poˇctu parametr˚ u Z rovnice (??) dostaneme pro vektor podm´ınˇen´ ych stˇredn´ıch hodnot E(y | X) = Xβ + E(ε | X) = Xβ

(6)

a pro kovarianˇcn´ı matici s vyuˇzit´ım rov. (1) { } cov(y | X) = E [y − E(y | X)][y − E(y | X)]T | X = { } { } = E [y − Xβ][y − Xβ]T | X = E (ε | X)(ε | X)T | X = σ 2 I

(7)

Neznámé parametry β lze odhadnout metodou nejmenˇs´ıch ˇctverc˚ u, tj. nalezen´ım takového vektoru b, pro kter´ y je nejmenˇs´ı tzv. residuáln´ı suma ˇctverc˚ u (RSS) odchylek pozorovan´ ych hodnot od jejich odhad˚ u z modelu RSS = (y − Xb)T (y − Xb)

(8)

Poloˇz´ıme-li derivaci v´ yrazu (8) podle vektoru b, tj. ∂ RSS ∂ = ∂b (yT y − bT XT y − yT X b + bT XT X b) = ∂b ∂ = ∂b (−2 bT XT y + bT XT X b) = −2 XT y + 2XT X b

rovnu nulovému vektoru, dostaneme soustavu normáln´ıch rovnic XT y = XT Xb

(9)

a vzhledem k platnosti pˇredpokladu (4) m˚ uˇzeme ˇreˇsen´ı této soustavy lineárn´ıch rovnic (vzhledem k b) vyjádˇrit explicitnˇe jako b = (XT X)−1 XT y

(10)

Odhady urˇcené podle (10) se naz´ yvaj´ı OLS – odhady (Ordinary Least Squares). Tyto odhady jsou nestranné, nebot’ E(b) = E(b|X) = (XT X)−1 XT E(y | X) = (XT X)−1 XT Xβ = β

(11)

Lze ukázat, ˇze tyto odhady maj´ı dalˇs´ı dobré vlastnosti, jsou to BLU-odhady (Best Linear Unbiased). Kovarianˇcn´ı matice tˇechto odhad˚ u je [ ]T cov(b) = (XT X)−1 XT σ 2 I (XT X)−1 XT = σ 2 (XT X)−1

(12)

5.2 Odhad parametr˚ u metodou maximáln´ı vˇerohodnosti

33

Na diagonále této matice jsou rozptyly odhadu parametr˚ u var(bi ). Nestrann´ y odhad parametru σ 2 je (d˚ ukaz viz napˇr. Andˇel 1978) s2 =

RSS n−k−1

(13)

a nestrann´ y odhad kovarianˇcn´ı matice odhad˚ u parametru b je Sbb = s2 (XT X)−1

(14)

Na diagonále matice Sbb jsou tedy nestranné odhady rozptyl˚ u odhadu parametr˚ u, jejich odmocninu (smˇerodatnou odchylku odhadu) oznaˇcme s(bi ). Pˇridáme-li k pˇredpoklad˚ um (1) aˇz (4) jeˇstˇe pˇredpoklad o tvaru rozdˇelen´ı náhodné sloˇzky modelu (3), resp. (??), a to (5)

εi ∼ N (0, σ 2 ),

i = 1, 2, · · · , n,

S vyuˇzit´ım pˇredpokladu (2) ε ∼ N (0, σ 2 I) pro následuj´ıc´ı statistiku plat´ı b − βi √i ∼ N (0, 1) var(bi ) a pro

i = 0, 1, · · · , k

b i − βi ∼ tn−k−1 . s(bi )

(15)

Tuto statistiku (15) pak m˚ uˇzeme uˇz´ıt ke stanoven´ı intervalu spolehlivosti pro parametr βi a testován´ı hypotéz o tomto parametru.

5.2

Odhad parametr˚ u metodou maxim´ aln´ı vˇ erohodnosti

Za pˇredpokladu (5) m˚ uˇzeme odvodit i maximálnˇe vˇerohodné (Maximum Likelihood) ML odhady pro klasick´ y model lineárn´ı regrese. ML-odhady odhaduj´ı hodnoty parametr˚ u tak, aby tyto odhady maximalizovaly tzv. vˇerohodnostn´ı funkci, tj. odhady jsou urˇceny jako nejpravdˇepodobnˇejˇs´ı hodnoty parametr˚ u pro pozorovaná data. MLodhady obecnˇe maj´ı ˇradu dobr´ ych vlastnost´ı:

– jsou asymptoticky nestranné (s rostouc´ım n jejich stˇredn´ı hodnota konverguje k odhadovanému parametru) – skoro vˇzdy jsou konsistentn´ı (s rostouc´ım n rozptyl odhadu konverguje k nule)

´ Í REGRESE 5 LINEARN

34

Vˇerohodnostn´ı funkce (souˇcin hustot jednotliv´ ych pozorován´ı) má pro klasick´ y model lineárn´ı regrese tvar (

LM L = 2πσ

) 2 −n/2

( T ) ε ε exp − 2 2σ

a jej´ı logaritmus je n Xβ) ln(LM L ) = L = (− ) ln (2πσ 2 ) − (y − Xβ)T (y − 2 2σ 2

(16)

Maximálnˇe vˇerohodn´ ymi odhady jsou takové odhady β M L , pro které vˇerohodnostn´ı funkce (16) nab´ yvá maxima. Pˇri hledán´ı maxima funkce (16) poloˇz´ıme derivace podle hledan´ ych promˇenn´ ych rovny nule, tedy ∂L = 0, ∂β

∂L =0 ∂σ 2

(17)

a dostaneme XT y = XT Xβ M L

(18)

Vid´ıme, ˇze ML-odhady regresn´ıch koeficient˚ u jsou v klasickém lineárn´ım modelu stejné jako odhady z´ıskané metodou nejmenˇs´ıch ˇctverc˚ u (OLS-odhady), β M L = b, srovnej rov.(18) s rov.(10). ML-odhad parametru σ 2 z rov.(17) je σM L =

1 RSS (y − Xβ)T (y − Xβ) = n n

tedy se liˇs´ı od OLS-odhadu v rov.(13), je pouze asymptoticky nestrann´ y.

5.2 Odhad parametr˚ u metodou maximáln´ı vˇerohodnosti

Shrnut´ı

• lineárn´ı regresn´ı model • pˇredpoklady v klasickém lineárn´ım regresn´ım modelu • metoda nejmenˇs´ıch ˇctverc˚ u • metoda maximáln´ı vˇerohodnosti • kovarianˇcn´ı matice odhad˚ u parametr˚ u Kontroln´ı ot´ azky

1. Jaký je rozd´ıl mezi parametry a jejich odhady? 2. Jsou odhady parametr˚ u náhodné veliˇciny? 3. Jaký je tvar kovarianˇcn´ı matice odhad˚ u parametr˚ u v klasickém modelu? 4. Jsou odhady z´ıskané metodou nejmeˇs´ıch ˇctverc˚ u nestranné? Dokaˇzte to. Korespondeˇ cn´ı u ´ loha Korespondenˇcn´ı u ´lohy budou zadáv´ any ke kaˇzdému kursu samostatnˇe.

35

37

6

Geometrie metody nejmenˇ s´ıch ˇ ctverc˚ u a regresn´ı diagnostika

Pr˚ uvodce studiem I tato kapitola je velmi d˚ uleˇzitá pro pochopen´ı princip˚ u statistické analýzy v´ıcerozmˇerných dat. Poˇc´ıtejte nejménˇe se ˇctyˇrmi hodinami usilovného studia s t´ım, ˇze se k prob´ırané látce budete podle potˇreby jeˇstˇe vracet po pochopen´ı dalˇs´ıch souvislost´ı.

6.1

Geometrie metody nejmenˇ s´ıch ˇ ctverc˚ u

Uvaˇzujeme klasick´ y lineárn´ı regresn´ı model y = Xβ + ε, ε ∼ N (0, σ 2 I) Jak v´ıme z pˇredchoz´ı kapitoly, odhad parametr˚ u m˚ uˇzeme vyjádˇrit explicitnˇe b = (XT X)−1 XT y. Vektor y ˆ = Xb je lineárn´ı kombinac´ı vektor˚ u regresor˚ u, tj. leˇz´ı v prostoru (pˇr´ımce, rovinˇe, nadrovinˇe), jehoˇz dimenze je rovna poˇctu regresor˚ u. Dosad´ıme-li za b, dostaneme y ˆ = Xb = X(XT X)−1 XT y = Hy Matice H = X(XT X)−1 XT je matice projekce vektoru y do prostoru urˇceného vektory regresor˚ u. Poˇzadavek formulovan´ y v metodˇe nejmenˇsich ˇctverc˚ u, tj. RSS = T (y − y ˆ) (y − y ˆ) vlastnˇe znamená, ˇze tato projekce je ortogonáln´ı. Pak tedy vektory y ˆ a e = y−y ˆ jsou ortogonáln´ı vektory, tzn. y ˆ T e = eT y ˆ = 0, o ˇcemˇz se velmi snadno m˚ uˇzeme pˇresvˇedˇcit: (Xb)T (y − Xb) = bT XT y − bT XT X)b = bT (XT y − XT Xb) = 0, nebot’ v´ yraz v posledn´ı závorce je nulov´ y vektor, viz normáln´ı rovnice (9). Vektoru e = y − y ˆ se ˇr´ıká vektor residu´ı, jeho sloˇzkám ei = yi − yî pak residua. Souˇcet a tedy i pr˚ umˇer residu´ı je roven nule: n ∑ i=1

ei =

n ∑ i=1

(yi − yî ) =

n ∑ i=1

yi −

n ∑ i=1

yî = 0,

38

ˇÍCH CTVERC ˇ ˚ 6 GEOMETRIE METODY NEJMENS U A REGRESNÍ DIAGNOSTIKA

¯ , kde x ¯ T = [1, x¯1 , x¯2 , . . . , x¯k ], tud´ıˇz nebot’ z prvn´ı normáln´ı rovnice plat´ı, ˇze y¯ = bT x n ∑

(ˆ yi − y¯) =

i=1

n ∑ i=1

yî −

n ∑

yi = bT

i=1

n ∑

¯ ) = 0, (xi − x

i=1

nebot’ souˇcet odchylek od pr˚ umˇeru je nulov´ y.

6.2

Rozklad souˇ ctu ˇ ctverc˚ u

Variabilitu vysvˇetlované veliˇciny m˚ uˇzeme vyjádˇrit jako souˇcet ˇctverc˚ u odchylek pozorovan´ ych hodnot od jejich pr˚ umˇeru. Tuto charakteristiku naz´ yváme celkov´ y souˇcet ˇctverc˚ u, TSS. TSS = (y − y ¯)T (y − y ¯) Lze ukázat, ˇze tuto celkovou sumu ˇctverc˚ u m˚ uˇzeme rozloˇzit na dvˇe sloˇzky MSS = (ˆ y−y ¯)T (ˆ y−y ¯) a uˇz dˇr´ıve definovanou RSS = (y − y ˆ)T (y − y ˆ ) = eT e Plat´ı tedy, ˇze TSS = MSS + RSS, MSS je ta ˇca´st z celkového souˇctu ˇctverc˚ u, která je vysvˇetlena závislost´ı vysvˇetlované veliˇciny na regresorech, zbylou ˇca´st (RSS) lineárn´ı závislost´ı vysvˇetlit nelze. Nyn´ı m˚ uˇzeme zavést d˚ uleˇzitou charakteristiku toho, jak u ´spˇeˇsnˇe regresn´ı model vysvˇetluje variabilitu vysvˇetlované veliˇciny. Této charakteristice se ˇr´ıká koeficient (index) determinace, R2 . R2 =

TSS − RSS RSS MSS = =1− TSS TSS TSS

Vid´ıme, ˇze 0 ≤ R2 ≤ 1. Hodnota indexu determinace R2 = 1, kdyˇz RSS = 0, tzn. regresn´ı model vysvˇetluje závislost vysvˇetlované veliˇciny na regresorech u ´plnˇe 2 (dokonalá lineárn´ı závislost). Naopak, R = 0, kdyˇz model nevysvˇetluje nic, tedy RSS = TSS, coˇz nastane jen tehdy, kdyˇz vˇsechny odhady b1 = b2 = . . . = bk = 0 a b0 = y¯, napˇr. pro k = 1 je regresn´ı pˇr´ımka rovnobˇeˇzná s osou x v u ´rovni b0 = y¯. Z rozkladu celkového souˇctu ˇctverc˚ u vycház´ı i anal´ yza rozptylu, která je obvyklou souˇcást´ı regresn´ıch program˚ u. Tabulka anal´ yzy rozptylu má vˇetˇsinou tento formát:

6.3 Regresn´ı diagnostika

zdroj variability model error total

39

stupnˇe volnosti k

souˇcet ˇctverc˚ u MSS

pr˚ umˇern´ y ˇctverec MSS/k

n−k−1 n−1

RSS TSS

RSS/(n − k − 1)

F MSS/k RSS/(n−k−1)

p−value 0. . . .

Za pˇredpokladu, ˇze pr˚ umˇern´ y ˇctverec s2 = RSS/(n − k − 1) je opravdu nestrann´ ym odhadem rozptylu náhodné sloˇzky (σ 2 ), tzn. v modelu jsou zaˇrazeny vˇsechny relevantn´ı regresory a RSS nen´ı zvˇetˇseno systematickou závislost´ı na nezaˇrazeném regresoru (podrobnˇeji viz napˇr. Draper a Smith, kap. 2 a 24) a náhodné kol´ısán´ı má normáln´ı rozdˇelen´ı, má statistika F rozdˇelen´ı F ∼ Fk,n−k−1 a m˚ uˇzeme ji uˇz´ıt k testu hypotézy H0 : β1 = β2 = . . . = βk = 0 proti H1 : aspoˇ n jeden parametr βj ̸= 0,

j = 1, 2, . . . , k

Povˇsimnˇeme si, ˇze d˚ uleˇzitou informaci o variabilitˇe residu´ı ei = yi − y¯i a t´ım i o shodˇe modelem predikovan´ ych hodnot y¯i s pozorovan´ ymi hodnotami yi nám poskytuje smˇerodatná odchylka residu´ı (square root mean error) √ s=

RSS n−k−1

Index determinace má tendenci nadhodnocovat pod´ıl modelu na vysvˇetlen´ı celkové variability veliˇciny y, mimo jiné i proto, ˇze kv˚ uli náhodnému kol´ısán´ı jsou odhady bj ̸= 0 i tehdy, kdyˇz βj = 0, j = 1, 2, . . . , k. Proto se zavád´ı tzv. adjustovan´ y 2 index determinace Radj , 2 Radj =1−

RSS/(n − k − 1) n−1 =1− (1 − R2 ) TSS/(n − 1) n−k−1

2 Vid´ıme, ˇze Radj < R2 , rozd´ıl je v´ yrazn´ y tehdy, kdyˇz poˇcet pozorován´ı je jen o málo 2 vˇetˇs´ı neˇz poˇcet regresor˚ u v modelu. Naopak hodnota Radj se pˇribliˇzuje R2 pro n ≫ k.

6.3

Regresn´ı diagnostika

Dalˇs´ı informace o vhodnosti modelu a o tom, zda jsou splnˇeny pˇredpoklady uˇcinˇené pro klasick´ y lineárn´ı model m˚ uˇzeme z´ıskat z anal´ yzy residu´ı. Vektor residu´ı m˚ uˇzeme vyjádˇrit pomoc´ı projekˇcn´ı matice H: e=y−y ˆ = Iy − Hy = (I − H)y

40


Pak kovarianˇcn´ı matice residu´ı je cov(e) = cov [(I − H)y] = (I − H)cov(y)(I − H)T = (I − H)σ 2 I(I − H)T = σ 2 I(I − H)T = σ 2 (I − H)(I − H)T = σ 2 (I − H − HT + HHT ) = σ 2 (I − H) nebot’ projekˇcn´ı matice H je symetrická (HT = H) a idempotentn´ı (H2 = H): HHT = H2 = X(XT X)−1 XT X(XT X)−1 XT = H Poznámka – vektor rezidu´ı e je náhodn´ y vektor (je to v´ yraz, ve kterém jsou náhodné vektory y a b). Matice H s prvky hij , i, j = 1, 2, . . . , n je symetrická, ale nemus´ı b´ yt diagonáln´ı. Jak bylo v pˇredchoz´ım odstavci ukázáno, kovarianˇcn´ı matice vektoru residu´ı je rovna cov(e) = σ 2 (I − H) Nestrann´ ym odhadem parametru σ 2 je reziduáln´ı rozptyl (tzn. rozptyl εi ): s2 =

1 eT e n−k−1

Dále uvedeme nˇekteré charakteristiky, které se uˇz´ıvaj´ı v tzv. regresn´ı diagnostice, tj. pˇri anal´ yze vhodnosti modelu. Klasick´ a residua Jsou to residua, kter´ y uˇz jsme se zab´ yvali, e = y − Xb. Jejich rozptyly var(ei ) = s2e (1 − hii ), nejsou konstantn´ı, i kdyˇz var(ϵi ) = σ 2 konstantn´ı je. Normovan´ a residua Jsou to klasická residua, vydˇelená reziduáln´ı smˇerodatnou odchylkou: eN i =

ei s

Jejich rozptyl je roven var(eN i ) = 1 − hii , tedy nemus´ı b´ yt roven jedné.

6.3 Regresn´ı diagnostika

41

Standardizovan´ a rezidua Nˇekdy se jim ˇr´ıká vnitˇrnˇe studentizovaná residua (internally studentized), jsou definována takto: ei eSi = √ s 1 − hii a jejich rozptyl je kostantn´ı, roven jedné. Plnˇ e studentizovan´ a rezidua Podle techniky uˇzité v jejich definici se jim ˇr´ıká také JACKKNIFE residua, jsou konstruována tak, ˇze vˇzdy pro i−t´ y bod se residuum poˇc´ıtá z modelu, jehoˇz parametry byly odhadnuty ze zb´ yvaj´ıc´ıch n − 1 bod˚ u, tedy vˇzdy i−t´ y bod se vypust´ı. eJi =

e(−i) √ . s(−i) 1 − hii

kde s(−i) je residuáln´ı smˇerodatná odchylka pˇri vynechán´ı i-tého bodu (ˇrádku datové matice). Tato residua maj´ı t−rozdˇelen´ı, eJi ∼ t(n − k − 2). Leverage Tyto charakteristiky ohodnocuj´ı vliv i-tého bodu na hodnoty odhad˚ u parametr˚ u. Jsou to diagonáln´ı prvky projekˇcn´ı matice, tedy hodnoty hii . Plat´ı, ˇze 0 < hii < 1

a

n ∑

hii = k + 1,

i=1

kde k je poˇcet regresor˚ u. Hodnota hii je u ´mˇerná vzdálenosti i-tého pozorován´ı od tˇeˇziˇstˇe (v k-rozmˇerném prostoru regresor˚ u), hii se povaˇzuje za velké, kdyˇz hii je vˇetˇs´ı neˇz dvojnásobek pr˚ umˇerné hodnoty, tj. hii > 2(k + 1)/n). Cookova vzd´ alenost Tato charakteristika slouˇz´ı také k posouzen´ı vlivu i-tého pozorován´ı na odhady parametr˚ u modelu, tj. na hodnoty b. Je to vlastnˇe relativn´ı zmˇena reziduáln´ıho souˇctu ˇctverc˚ u zp˚ usobená vypuˇstˇen´ım i-tého pozorován´ı. Cookova vzdálenost pro i-té pozorován´ı je definována Ci =

(y−ˆ y(−i) )T (y−ˆ y(−i) ) ps2

=

(b−b(−i) )T (X T X)(b−b(−i) ) ps2

=

hii e2 p(1−hii ) Si

kde b(−i) jsou jsou jackknife odhady (spoˇc´ıtané pˇri vypuˇstˇen´ı i-tého bodu) a p je poˇcet odhadovan´ ych parametr˚ u. Cookova vzdálenost ohodnocuje vliv i-tého pozorován´ı na odhad vektoru regresn´ıch parametr˚ u b. Je-li Cookova vzdálenost Ci ≥ 1, i-pozorován´ı velmi podstatnˇe ovlivˇ nuje odhady parametr˚ u.

42

6.4


Autokorelace

Pˇri posuzován´ı pˇredpokladu o nekorelovanosti residu´ı se obvykle vycház´ı modelu autokorelaˇcn´ıho procesu prvn´ıho ˇrádu – AR(1): εi = ρ1 εi−1 + ui

kde

ui ∼ N (0, σ 2 )

Autokorelaˇcn´ı koeficient prvn´ıho ˇrádu ρ1 odhadujeme jako ∑n i=2 ei ei−1 ρˆ1 = ∑ n 2 i=1 ei K testován´ı korelovanosti residu´ı se pak uˇz´ıvá Wald˚ uv test Wa =

nˆ ρ21 ∼ χ2 (1) 1 − ρˆ21

nebo ve statistickém software bˇeˇznˇe implementovaná Durbin – Watsonova statistika ∑n DW =

(ei − ei−1 )2 i=2∑ n 2 i=1 ei

≃ 2 (1 − ρˆ1 )

Pro tuto statistiku plat´ı 0 ≤ DW ≤ 4, E(DW ) = 2 pˇri ρ1 = 0. Kvantily této statistiky je obt´ıˇzné vyjádˇrit explicitnˇe, proto pro Durbin–Watson˚ uv test statistické programy neposkytuj´ı u jin´ ych test˚ u obvykl´ y komfort, totiˇz i dosaˇzenou v´ yznamnost (p). Pˇri rozhodován´ı je pro hodnoty statistiky velmi bl´ızké dvˇema spoléhat na intuici a povaˇzovat residua za nekorelované. Pro serióznˇejˇs´ı u ´sudek lze vyuˇz´ıt pˇribliˇzné kritické hodnoty, které jsou tabelovány, napˇr. [16].

6.4 Autokorelace

43

Pˇ r´ıklad 6.1 Data pro tento pˇr´ıklad jsou pˇrevzata z knihy [8], pˇr´ıklad 01A, a jsou uvedena i v následuj´ıc´ı tabulce. Veliˇcina y je mˇes´ıˇcn´ı spotˇreba páry ve firmˇe, veliˇcina x6 je poˇcet pracovn´ıch dn´ı v mˇes´ıci a veliˇcina x8 je venkovn´ı teplota ve stupn´ıch ´ Fahrenheita. Uloha, kterou máme ˇreˇsit, je odhadnout parametry lineárn´ıho regresn´ıho modelu a posoudit, zda je tento model vhodn´ y pro vysvˇetlen´ı závislosti y na x6 a x8. V ˇreˇsen´ı tohoto pˇr´ıkladu byl uˇzit statistick´ y programov´ y systém NCSS [14], zejména modul Multiple Regression, old version. V´ ystupy uvád´ıme bez vˇetˇs´ıch editaˇcn´ıch u ´prav v surovém stavu.

i 1 2 3 4 5 6 7 8 9 10 11 12 13

y 10.98 11.13 12.51 8.40 9.27 8.73 6.36 8.50 7.82 9.14 8.24 12.19 11.88

x6 20 20 23 20 21 22 11 23 21 20 20 21 21

x8 35.3 29.7 30.8 58.8 61.4 71.3 74.4 76.7 70.7 57.5 46.4 28.9 28.1

i 14 15 16 17 18 19 20 21 22 23 24 25

y 9.57 10.94 9.58 10.09 8.11 6.83 8.88 7.68 8.47 8.86 10.36 11.08

x6 19 23 20 22 22 11 23 20 21 20 20 22

x8 39.1 46.8 48.5 59.3 70.0 70.0 74.5 72.1 58.1 44.6 33.4 28.6

Vyˇsetˇrujeme-li nˇejakou závislost, vˇzdy je dobré data nejdˇr´ıve prohlédnout jednoduch´ ymi prostˇredky popisné statistiky. Ty nám ˇcasto pomohou odhalit zaj´ımavé vˇeci v datech, napˇr. odlehlé hodnoty. To m˚ uˇzeme vidˇet i na následuj´ıc´ım grafu, kde pozorován´ı na ˇra´dc´ıch 7 a 19 jsou zcela mimo hodnoty ostatn´ıch pozorován´ı (patrnˇe je to poˇcet pracovn´ıch dn˚ u v mˇes´ıc´ıch, kdy byly dovolené a ve firmˇe se nepracovalo). Tato pozorován´ı jsou apriori podezˇrelá a pˇri diagnostice modelu je nutné si na nˇe dát pozor.

44


Dalˇs´ı uˇziteˇcné nahlédnut´ı do analyzovan´ ych dat je v´ ybˇerová korelaˇcn´ı matice, která je volitelnou souˇca´st´ı v´ ystupu z modulu Multiple Regression: Multiple Regression Report Dependent y Correlation Matrix Section x6 x8 y x6 1.000000 -0.209761 0.536122 x8 -0.209761 1.000000 -0.845244 y 0.536122 -0.845244 1.000000 Vid´ıme, ˇze regresory x6 a x8 jsou jen slabˇe korelovány (korelaˇcn´ı koeficient je 0,21), takˇze matice regresor˚ u má plnou hodnost, nehroz´ı numerické pot´ıˇze spojené se ˇspatnou podm´ınˇenost´ı. Dalˇs´ı pro nás d˚ uleˇzitou souˇcást´ı v´ ystupu z modulu Multiple Regression je Regression Equation Section, kde jsou odhady parametr˚ u modelu. Regression Equation Section Independent Regression Variable Coefficient Intercept 9.12688 x6 0.2028154 x8 -7.239294E-02

Standard T-Value Error (Ho: B=0) 1.102801 8.2761 4.576761E-02 4.4314 7.999381E-03 9.0498

Prob Level 0.000000 0.000210 0.000000

6.4 Autokorelace

45

Vid´ıme, ˇze u vˇsech tˇr´ı odhadovan´ ych parametr˚ u zam´ıtáme nulovou hypotézu. Model tedy neobsahuje nadbyteˇcné parametry. Jak vid´ıme v následuj´ıc´ı tabulce ANOVA, model vysvˇetluje v´ yznamnou ˇcást z cel2 kové varibility veliˇciny y, podle hodnoty R = 0.8491 asi 85% z celkové varibility. Odhad residuáln´ı smˇerodatné odchylky je uveden jako Root Mean Square Error a je roven pˇribliˇznˇe 0,66. Druhá mocnina této charakteristiky je pak odhadem residuáln´ıho rozptylu σ 2 . Analysis of Variance Section Sum of Mean Source DF Squares Square F-Ratio Intercept 1 2220.294 2220.294 Model 2 54.1871 27.09355 61.9043 Error 22 9.628704 0.4376684 Root Mean Square Error Mean of Dependent R-Squared Adj R-Squared

Prob Level 0.000000

0.6615651 9.424 0.8491 0.8354

Pro posouzen´ı vhodnosti modelu jsou d˚ uleˇzité v´ ystupy z regresn´ı diagnostiky. Durbin-Watsonova statistika je velmi bl´ızká hodnotˇe 2, tak se nemus´ıme znepokojovat autokorelac´ı residu´ı. Durbin-Watson Value

2.1955

V následuj´ıc´ı tabulce jsou jackknife residua oznaˇcena jako Rstudent. Pozornost vˇenujme pˇredevˇs´ım ˇrádk˚ um 7 a 19, které byly podezˇrelé uˇz pˇri pˇredbˇeˇzné jednoduché anal´ yze a pak tˇem ˇrádk˚ um, kde studentizovaná residua jsou v absolutn´ı hodnotˇe velká (zhruba vˇetˇs´ı neˇz 2, coˇz je pˇribliˇzná hodnota kvantilu tn−k−1 (0.975)).

46


Regression Diagnostics Section Studentized Row Residual Rstudent 1 0.556825 0.547897 2 0.156002 0.152500 3 1.531855 1.583464 4 -0.814515 -0.808066 5 0.511834 0.503070 6 0.487210 0.478597 7 0.789332 0.782341 8 0.436764 0.428584 9 -0.711757 -0.703540 10 0.184529 0.180426 11 -2.452153 -2.810441 12 1.442820 1.481481 13 0.853098 0.847621 14 -0.913679 -0.910106 15 0.843302 0.837562 16 -0.142369 -0.139160 17 1.241672 1.258005 18 -0.658977 -0.650276 19 1.086621 1.091328 20 0.798049 0.791237 21 -0.450525 -0.442212 22 -1.100280 -1.105840 23 -1.697613 -1.779205 24 -0.644223 -0.635433 25 -0.708085 -0.699826

Hat Diagonal 0.085491 0.118877 0.124826 0.045374 0.056434 0.117502 0.447410 0.184719 0.095491 0.043373 0.046418 0.118566 0.123991 0.079880 0.075758 0.043079 0.065527 0.109838 0.436460 0.167792 0.094228 0.048654 0.050307 0.095790 0.124217

Cook’s D 0.009662 0.001094 0.111564 0.010511 0.005223 0.010535 0.168151 0.014407 0.017827 0.000515 0.097567 0.093341 0.034337 0.024158 0.019431 0.000304 0.036036 0.017861 0.304828 0.042803 0.007039 0.020638 0.050886 0.014656 0.023705

Nejvˇetˇs´ı residuum v absolutn´ı hodnotˇe má pozorován´ı 11, ale jak vid´ıme z ostatn´ıch statistik, neovlivˇ nuje nijak v´ yznamnˇe hodnoty odhad˚ u, je to odlehl´ y bod, kter´ y zvˇetˇsuje residuáln´ı rozptyl. Statistiky hii maj´ı nejvˇetˇs´ı pozorován´ı 7 a 19, jejich hodnoty jako jediné pˇresahuj´ı mezn´ı hodnotu 2p/n = 6/25, maj´ı i nejvˇetˇs´ı Cookovu vzdálenost, ale zdaleka nedosahuj´ıc´ı hodnotu 1. Tyto dva body jsou tzv. vlivné, ale nevyboˇcuj´ıc´ı. Naopak pˇrisp´ıvaj´ı ke sn´ıˇzen´ı residuáln´ıho rozptylu. Následuj´ıc´ı QQ graf residu´ı ukazuje, ˇze rozdˇelen´ı residu´ı m˚ uˇzeme povaˇzovat za normáln´ı, odchylky od pˇr´ımky nejsou velké. Tedy data i model vyhovuj´ı pˇredpokladu (5) klasického modelu a v´ ysledky test˚ u o parametrech modelu m˚ uˇzeme povaˇzovat za spolehlivé (nezavádˇej´ıc´ı).

6.4 Autokorelace

47

Graf residu´ı proti odhadovan´ ym hodnotám yˆ ukazuje, ˇze residuáln´ı rozptyl m˚ uˇzeme povaˇzovat za konstatn´ı, kaz´ı“ to jen bod 11 s residuem menˇs´ım neˇz -2. ”

48


V´ ysledky tedy m˚ uˇzeme shrnout takto: Hodnoty veliˇciny y (mˇes´ıˇcn´ı spotˇreby páry ve firmˇe) lze uspokojivˇe odhadovat lineárn´ı závislost´ı na veliˇcinˇe x6 (poˇcet pracovn´ıch dn´ı v mˇes´ıci) a veliˇciny x8 (vnˇejˇs´ı teplota). Oˇcekávanou spotˇrebu páry lze vyjádˇrit vztahem yˆ = 9.127 + 0.2028 x6 − 0.07239 x8. Smˇerodatná odchylka pˇredpovˇedi je 0.662, model vysvˇetluje 85% z celkové variability spotˇreby páry.

6.4 Autokorelace

49

Pˇ r´ıklad 6.2 V´ yznam statistik pro diagnostiku ukazuje jednoduch´ y pˇr´ıklad modelu s jedn´ım regresorem. Na obrázc´ıch jsou data a proloˇzené regresn´ı pˇr´ımky.

Odliˇsnost je jen v hodnotách vertikáln´ı souˇradnice bodu 4. V obou u ´lohách maj´ı body 4 velké hodnoty hii , ale liˇs´ı se v hodnotách ostan´ıch diagnostick´ ych statistik. V pˇr´ıpadˇe prvn´ım má bod 4 malé jackknife residuum i Cookovu vzdálenost, tzn. jeho vypuˇstˇen´ım ˇci pˇridán´ım se hodnoty odhad˚ u pˇr´ıliˇs nemˇen´ı: Row Residual Rstudent 4 -0.029151 -0.026991

Diagonal 0.605906

Cook’s D 0.000653

50


Ve druhém pˇr´ıpadˇe jsou jackknife residuum i Cookova vzdálenost extrémnˇe velké. Row Residual Rstudent 4 -2.382046 -5.067308

Diagonal 0.605906

Cook’s D 4.361897

Tento bod tedy má zásadn´ı vliv na hodnoty odhad˚ u, coˇz je ostatnˇe na prvn´ı pohled vidˇet v této u ´loze, kdy model má je jeden regresor, i na grafech proloˇzen´ ych regresn´ıch pˇr´ımek. V pˇr´ıpadˇe model˚ u s v´ıce regresory ovˇsem takové vizuáln´ı posouzen´ı nen´ı moˇzné a regresn´ı diagnostika je uˇziteˇcn´ ym nástrojem pro ovˇeˇren´ı pˇredpoklad˚ u a posouzen´ı vhodnosti modelu. D˚ usledky um´ıstˇen´ı bodu 4 vid´ıme v následuj´ıc´ı tabulce smˇerodatn´ ych odchylek:

intercept smˇ ernice sm. odch. residu´ ı

1 (vlivn´ y) 0.849 0.130 1.15

2 (odlehl´ y) 1.951 0.299 2.65

bez bodu 4 1.152 0.211 1.25

6.4 Autokorelace

51

Shrnut´ı

• projekˇcn´ı matice, ortogonáln´ı projekce • rozklad celkového souˇctu ˇctverc˚ u, index determinace R2 • residua, jackknife residua, diagonáln´ı prvky projekˇcn´ı matice, Cookova vzdálenost • autokorelace, Durbin–Watsonova statistika Kontroln´ı ot´ azky

1. Zkuste dokázat, ˇze opravdu plat´ı TSS = MSS + RSS. 2. Naˇcrtnˇete, co znamená ortogonalita vektor˚ u y ˆ a e. Pro graf zvolte model se dvˇema regresory a výbˇer o rozsahu 3. 3. Jaká hypotéza se testuje v analýze rozptylu, uvádˇené ve výstupu statistických program˚ u pro lineárn´ı regresi? 4. K ˇcemu slouˇz´ı regresn´ı diagnostika? Korespondenˇ cn´ı u ´ loha Korespondenˇcn´ı u ´lohy budou zadáv´ any ke kaˇzdému kursu samostatnˇe.

52

7

´ Í A MNOHONASOBN ´ ´ KORELACE 7 PARCIALN A

Parci´ aln´ı a mnohon´ asobn´ a korelace

Pr˚ uvodce studiem Kapitola uvád´ı nˇekteré charakteristiky vyjadˇruj´ıc´ı vztahy ve v´ıcerozmˇerných datech. Jejich pochopen´ı vám bude uˇziteˇcné pˇri studiu dalˇs´ıch metod analýzy dat. Na tuto kapitolu poˇc´ıtejte s jednou aˇz dvˇema hodinami studia. Jak v´ıme, korelaˇcn´ı koeficient cov(X, Y ) ρ(X, Y ) = √ varXvarY vyjadˇruje m´ıru lineárn´ı závislosti dvou náhodn´ ych veliˇcin X, Y . Plat´ı, ˇze −1 ≤ ρ(X, Y ) ≤ 1. Nyn´ı uvaˇzujme vektor náhodn´ ych veliˇcin [Y, X1 , X2 , · · · , Xk ]. Parci´ aln´ı korelaˇ cn´ı koeficient veliˇcin Y, X1 pak zap´ıˇseme jako ρ(Y, X1 |X2 , X3 , · · · , Xk ) = ρY X1 ·X2 X3 ···Xk a znamená vlastnˇe korelaˇcn´ı koeficient mezi náhodn´ ymi veliˇcinami, které zbudou“ ” z veliˇcin Y, X1 po odeˇcten´ı regresn´ı funkce tˇechto veliˇcin na podmiˇ nuj´ıc´ıch veliˇcinách, tj. Y = α1 + α2 X2 + α3 X3 + · · · + αk Xk + ε(1)

X1 = β1 + β2 X2 + β3 X3 + · · · + βk Xk + ε(2) pak parciáln´ı korelaˇcn´ı koeficient je (obyˇcejn´ y) korelaˇcn´ı koeficient residuáln´ıch náhodn´ ych sloˇzek: ( ) ρY X1 ·X2 X3 ···Xk = ρ ε(1) , ε(2) . Parciáln´ı korelaˇcn´ı koeficienty lze vyjádˇrit pomoc´ı korelaˇcn´ıch koeficient˚ u dvojic veliˇcin náhodného vektoru. Máme-li korelaˇcn´ı matici   1 ρ(Y, X1 ) ρ(Y, X2 ) · · · ρ(Y, Xk )   1 ρ(X1 , X2 ) · · · ρ(X1 , Xk )   ρ(X1 , Y )    ρ(X , Y ) ρ(X , X ) 1 · · · ρ(X , X ) ϱ= 2 2 1 1 k   .. .. .. ..   .. .   . . . . ρ(Xk , Y ) ρ(Xk , X1 ) ρ(Xk , X2 ) · · ·

1

53

pak

|ϱY X1 | ρY X1 ·X2 X3 ···Xk = √ |ϱY,Y ||ϱX1 X1 |

kde |ϱU V | je determinant matice, která vznikne z korelaˇcn´ı matice ϱ, kdyˇz je vypuˇstˇen ˇra´dek U a sloupec V . Napˇr. pro k = 2 má korelaˇcn´ı matice tvar 

 1 ρ(Y, X1 ) ρ(Y, X2 )   1 ρ(X1 , X2 )   ρ(X1 , Y ) ρ(X2 , Y ) ρ(X2 , X1 ) 1 a parciáln´ı korelaˇcn´ı koeficient ρY X1 ·X2 je ρY X1 ·X2 =

ρ(X1 , Y ) − ρ(X2 , Y )ρ(X1 , X2 ) . [(1 − ρ2 (X1 , X2 ))(1 − ρ2 (Y, X2 ))]1/2

Mnohon´ asobn´ y (celkov´ y) koeficient korelace naz´ yvan´ y také celkový koeficient korelace vyjadˇruje korelaci náhodné veliˇciny Y na veliˇcinách [X1 , X2 , · · · , Xk ]. Vyjádˇr´ıme-li náhodnou veliˇcinu Yˆ jako regresn´ı funkci veliˇcin [X1 , X2 , · · · , Xk ], Yˆ = β0 + β1 X1 + · · · + βk Xk (

) ˆ pak jednoduch´ y korelaˇcn´ı koeficient ρ Y, Y je koeficient mnohonásobné korelace ρY ·X1 X2 ···Xk , ( ) ρ Y, Yˆ = ρY ·X1 X2 ···Xk . Koeficient mnohonásobné korelace lze spoˇc´ıtat z korelaˇcn´ı matice, ρY ·X1 X2 ···Xk

( )1/2 |ϱ| = 1− |ϱY Y |

Pro hodnoty mnohonásobného korelaˇcn´ıho koeficientu plat´ı 0 ≤ ρY ·X1 X2 ···Xk ≤ 1, pˇri ˇcemˇz nule je roven, kdyˇz vˇsechny dvojice Y, Xi , i = 1, 2, · · · , k jsou nekorelované, tedy ρ(Y, Xi ) = 0 pro i = 1, 2, . . . , k. Jedniˇcce je mnohonásobn´ y koreˆ laˇcn´ı koeficient roven tehdy, kdyˇz Y = Y , tj. pro ˇcistou“ lineárn´ı závislost Y na ” [X1 , X2 , · · · , Xk ].

´ Í A MNOHONASOBN ´ ´ KORELACE 7 PARCIALN A

54

Dále plat´ı: • ρY ·X1 X2 ···Xk ≥ max [ρ(Xi , Y )] i=1,··· ,k

• ρ(Y, X1 ) ≤ ρY ·X1 X2 ≤ · · · ≤ ρY ·X1 X2 ···Xk Analogick´ ym zp˚ usobem jsou definovány i výbˇerové koeficienty parci´ aln´ı a mnohonásobné korelace, jen se poˇc´ıtaj´ı z v´ ybˇerov´ ych korelaˇcn´ıch koeficient˚ u (z v´ ybˇerové korelaˇcn´ı matice R). K pochopen´ı toho, co znamenaj´ı v´ ybˇerové koeficienty parciáln´ı a mnohonásobné korelace nám pom˚ uˇze jejich souvislost s lineárn´ım regresn´ım modelem (i kdyˇz, pˇr´ısnˇe vzato, v klasickém lineárn´ım modelu uvaˇzujeme, ˇze regresory jsou deterministické a korelace se t´ yká náhodn´ ych veliˇcin, ale chápejme tyto charakteristiky tak, ˇze jsou poˇc´ıtány podle formul´ı zaveden´ ych pro v´ ybˇerové korelaˇcn´ı koeficienty a v´ ybˇerové kovariance). Pro model s jedn´ım regresorem Yˆ = a0 + a1 x1 plat´ı, ˇze korelaˇcn´ı koeficient ryx1 = ssx1yx1sy smˇernice regresn´ı pˇr´ımky koeficient determinace

syx1 s2x1 2 2 R = ry,x1

a1 =

Pro v´ ybˇerov´ y koeficient mnohonásobné korelace v modelu s k regresory plat´ı, ˇze jeho druhá mocnina je rovna koeficientu determinace, tedy 2 ry·x = R2 = 1 − 1 x2 ···xk

RSS TSS

a pro v´ ybˇerov´ y koeficient parciáln´ı korelace plat´ı to, ˇze jeho druhá mocnina je rovna relativn´ı zmˇenˇe residuáln´ı sumy ˇctverc˚ u, napˇr. 2 ryx = 2 ·x1

∆RSS(x2 |x1 ) , RSS(x1 )

kde jmenovatel je residuáln´ı souˇcet ˇctverc˚ u modelu s jedn´ım regresorem x1 (Yˆ = a0 +a1 x1 ) a ˇcitatel zlomku je sn´ıˇzen´ı residuáln´ıho souˇctu ˇctverc˚ u zp˚ usobené pˇridán´ım regresoru x2 do modelu s jedn´ım regresorem, tj. ∆RSS(x2 |x1 ) = RSS(x1 ) − RSS(x1 , x2 ), kde RSS(x1 , x2 ) je )residuáln´ı souˇcet ˇctverc˚ u modelu se dvˇema regresory ( Yˆ = b0 + b1 x1 + b2 x2 .

55

Shrnut´ı

• koeficient parciáln´ı korelace, obor jeho hodnot • koeficient mnohonásobné korelace, obor jeho hodnot • vztah koeficientu mnohonásobné korelace a indexu determinace Kontroln´ı ot´ azky

1. Vyjádˇrete slovnˇe, co charakterizuje koeficient parci´ aln´ı korelace 2. Vyjádˇrete slovnˇe, co charakterizuje koeficient mnohonásobné korelace

56

8

´ ER ˇ REGRESOR˚ ˇ ´ REGRESI 8 VYB U V MNOHOROZMERN E

V´ ybˇ er regresor˚ u v mnohorozmˇ ern´ e regresi

Pr˚ uvodce studiem Kapitola se zabývá d˚ uleˇzitou a ˇcasto aplikovanou u ´lohou hled´ an´ı vhodného regresn´ıho modelu. Na tuto kapitolu poˇc´ıtejte nejménˇe se tˇremi aˇz ˇctyˇrmi hodinami studia. D˚ ukladnˇe prostudujte a promyslete i ˇreˇsený pˇr´ıklad na konci této kapitoly. Pomˇernˇe ˇcasto se v anal´ yze dat setkáváme s u ´lohami, které formálnˇe mohou b´ yt zapsány jako klasick´ y lineárn´ı regresn´ı model, y = Xβ + ε, ale v matici X typu n × (p + 1) je poˇcet regresor˚ u p velk´ y. Velk´ y poˇcet regresor˚ u má ˇcasto za následek, ˇze pak pro mnoho z parametr˚ u βi , i = 1, 2, . . . , p nem˚ uˇzeme zam´ıtnout H0 : βi = 0, tzn. i-t´ y regresor nevysvˇetluje zmˇeny hodnot veliˇciny y. Naˇs´ım c´ılem je naj´ıt jednoduˇsˇs´ı model s k regresory (k < p), obsahuj´ıc´ı jen takové regresory, které v´ yznamnou mˇerou vysvˇetluj´ı variabilitu hodnot y. ˇ sit takou u Reˇ ´lohu prozkoumán´ım vˇsech lineárn´ıch model˚ u k regresory, k = 1, 2, . . . , p, je pro vˇetˇs´ı hodnoty p ˇcasovˇe ne´ unosné. Znamenalo by to prozkoumat p ( ) ∑ p k=0

k

= 2p

model˚ u, tedy napˇr. jen pro p = 10 v´ ypoˇcet a interpretaci v´ ysledk˚ u odhad˚ u 1024 model˚ u, coˇz by pˇredstavovalo práci na nˇekolik t´ ydn˚ u. Nav´ıc pro velké hodnoty p b´ yvá ˇcasto matice X ˇspatnˇe podm´ınˇen´ a, tzn. determinant T . X X = 0 a odhady parametr˚ u jsou pak numericky nestabiln´ı a maj´ı velké rozptyly, takˇze v´ ysledky nejsou prakticky vyuˇzitelné.

8.1

Krokov´ a regrese

Jedn´ım ze zp˚ usob˚ u, jak nalézt podmnoˇzinu k regresor˚ u do vhodného lineárn´ıho regresn´ıho modelu, je kroková (stepwise) regrese. Jeˇstˇe dˇr´ıve, neˇz vysvˇetl´ıme tento postup, pˇripomeneme základn´ı pojmy a myˇslenky, na kter´ ych je zaloˇzen. V´ıme, ˇze celkov´ y souˇcet ˇctverc˚ u lze rozloˇzit: TSS = MSS + RSS

8.1 Kroková regrese

57

Dále, i modelovou sumu MSS m˚ uˇzeme rozloˇzit. Pˇredstavme si model s k regresory. Potom ˇca´st MSS pˇripadaj´ıc´ı i-tému regresoru, MSS(i · 1, 2, . . . , i − 1, i + 1, . . . , k),

oznaˇcme ji zkratkou MSSk(−i)

MSSk(−i) je tedy rozd´ıl modelové sumy ˇctverc˚ u MSS(k) pˇri zaˇrazen´ı vˇsech k regresor˚ u a modelové sumy ˇctverc˚ u MSS(k − 1) s (k − 1) regresory (i−t´ y regresor vynechán): MSSk(−i) = MSS(k) − MSS(k − 1). Pˇridán´ım i−tého regresoru se souˇcasnˇe zmˇen´ı odpov´ıdaj´ıc´ım zp˚ usobem i residuáln´ı souˇcet ˇctverc˚ u ∆RSSi = RSS(k − 1) − RSS(k) = MSS(k) − MSS(k − 1) = MSSk(−i) , nebot’ plat´ı, ˇze TSS = MSS(k) + RSS(k) = MSS(k − 1) + RSS(k − 1). Souˇcasnˇe v´ıme, ˇze ∆RSSi ≥ 0, tzn. pˇridán´ım regresoru se residuáln´ı souˇcet ˇctverc˚ u nezv´ yˇs´ı. Myˇslenka krokové regrese spoˇc´ıvá v tom, ˇze v kaˇzdém kroku (pˇredpokládejme, ˇze k − 1 regresor˚ u je uˇz zaˇrazeno v modelu) budeme z dosud nezaˇrazen´ ych vyb´ırat takov´ y regresor, kter´ y nejv´ıce sniˇzuje residuáln´ı sumu ˇctverc˚ u, tj. ten, jehoˇz ∆RSSi je nejvˇetˇs´ı. Pˇritom zaˇrad´ıme jen takov´ y regresor, kter´ y residuáln´ı sumu ˇctverc˚ u sniˇzuje v´ yznamnˇe. Kritériem (statistikou) pro posouzen´ı v´ yznamnosti je tzv. parci´ aln´ı F , coˇz je ∆MSSi ∼ F1,ν , s2 kde ∆MSSi oznaˇcuje zv´ yˇsen´ı modelové sumy ˇctverc˚ u odpov´ıdaj´ıc´ı zaˇrazen´ı i−tého 2 regresoru z dosud v modelu nezaˇrazen´ ych, s je nestrann´ y odhad parametru σ 2 a ν je jeho poˇcet stupˇ n˚ u volnosti. Implementace procedury krokové regrese se mohou liˇsit v tom, jak´ ym zp˚ usobem je poˇc´ıtán s2 . Jedna z moˇznost´ı je poˇc´ıtat s2 z residuáln´ı sumy ˇctverc˚ u v aktuáln´ım kroku, tj. RSS(k − 1) . s2 = n−k Pak parciáln´ı Fi i−tého nezaˇrazeného regresoru lze urˇcit z parciáln´ıho a celkového korelaˇcn´ıho koeficientu Fi =

2 [riY ·(k−1) ]/(n − k)

1−

rY2 ·(k−1)

=

∆RSSi , RSS(k − 1)/(n − k)

58


kde riY ·(k−1) je parciáln´ı korelaˇcn´ı koeficient Y a xi po odeˇcten´ı vlivu“ (k − 1) uˇz ” zaˇrazen´ ych regresor˚ u a rY ·(k−1) je mnohonásobn´ y (celkov´ y) koeficient korelace Y s (k − 1) uˇz zaˇrazen´ ymi regresory. V k−tém kroku tedy zaˇrad´ıme ten regresor, kter´ y má nejvˇetˇs´ı parciáln´ı Fi , a to jen tehdy, je-li Fi vˇetˇs´ı, neˇz zadaná hodnota F-to-entry, kterou obvykle vol´ıme jako takov´ y kvantil F rozdˇelen´ı, aby parciáln´ı Fi a tud´ıˇz i zmˇena v residuáln´ım souˇctu ˇctverc˚ u byly v´ yznamné, tedy F-to-entry = F1,(n−k−1) (1 − α1 ), kde α1 je zvolená hladina v´ yznamnosti pro zaˇrazen´ı regresoru do modelu. Po zaˇrazen´ı i−tého regresoru m˚ uˇze kv˚ uli korelaci mezi zaˇrazen´ ymi regresory nastat situace, ˇze parciáln´ı F nˇekterého ze zaˇrazen´ ych regresor˚ u pˇrestane b´ yt v´ yznamné. Jin´ ymi slovy, vypuˇstˇen´ı tohoto regresoru z modelu pak nezv´ yˇs´ı v´ yznamnˇe residuáln´ı sumu ˇctverc˚ u, tzn. regresor je v modelu nadbyteˇcn´ y. Proto se po zaˇrazen´ı regresoru spoˇc´ıtaj´ı parciáln´ı Fi vˇsech dosud zaˇrazen´ ych regresor˚ u Fi ==

∆RSS(i) , s2

kde ∆RSS(i) znamená zmˇenu (zv´ yˇsen´ı) residuáln´ı sumy ˇctverc˚ u pˇri vypuˇstˇen´ı i−tého regresoru z modelu. Najde se nejmenˇs´ı z tˇechto parciáln´ıch Fi a posuzuje se, zda bychom vypuˇstˇen´ım tohoto regresoru zv´ yˇsili RSS jen nepodstatnˇe. Kriterium pro toto rozhodován´ı je to, zda minimáln´ı Fi je menˇs´ı neˇz zadaná hodnota F-to-remove. Vˇetˇsinou vol´ıme F-to-remove = F1,(n−k−1) (1−α2 ), kde α2 je zvolená hladina v´ yznamnosti pro vypuˇstˇen´ı regresoru z modelu. Abychom pˇredeˇsli moˇznosti nekoneˇcného cyklu zaˇrazován´ı a vyˇrazován´ı regresor˚ u, obvykle se vol´ı F-to-remove < F-to-entry, tj. α2 > α1 . Po tomto vysvˇetlen´ı tedy m˚ uˇzeme algoritmus krokové regrese zapsat takto:

krok 0: zvol model se ˇzádn´ ym regresorem , tj. yˆ = y¯, a z p nezaˇrazen´ ych regresor˚ u zvol ten, kter´ y má nejvˇetˇs´ı absolutn´ı hodnotu korelaˇcn´ıho koeficientu s vysvˇetlovanou veliˇcinou y (pˇri jednom zaˇrazeném regresoru 2 , tedy nejv´ıce korelovan´ y regresor nejv´ıce sniˇzuje residuáln´ı je R2 = rxy sumu ˇctverc˚ u) if Fi < F-to-entry then konec else k = 1 krok k: mezi nezaˇrazen´ ymi regresory vyber ten s nejvˇetˇs´ım Fi . if Fi < F-to-entry then konec else zaˇrad’ i−t´ y regresor, k = k + 1 mezi zaˇrazen´ ymi regresory najdi ten s nejmenˇs´ım Fi , if min Fi < F-to-remove then vyˇrad’ i−t´ y regresor, k = k − 1 go to krok k

8.2 Hledán´ı nejlepˇs´ı mnoˇziny regresor˚ u

59

Analogick´ y krokov´ y (stepwise) postup se, jak uvid´ıme dále, uˇz´ıvá nejen v lineárn´ı regresi, ale i v dalˇs´ıch metodách anal´ yzy mnohorozmˇern´ ych dat. Existuj´ı i nˇekteré dalˇs´ı varianty, tzv. postupn´ ych procedur v´ ybˇeru veliˇcin do modelu, napˇr. zpˇetná backward procedura, která vycház´ı z modelu, ve kterém je zaˇrazeno vˇsech p veliˇcin a postupnˇe vyˇrazuje nev´ yznamné, nebo dopˇredná forward procedura, která je podobná v´ yˇse popsané krokové proceduˇre, avˇsak neumoˇzn ˇuje vypouˇstˇen´ı zaˇrazen´ ych veliˇcin, které se stanou nev´ yznamné. Obecnˇe lze ˇr´ıci, ˇze krokové procedury jsou uˇziteˇcn´ ym nástrojem pro hledán´ı vhodn´ ych model˚ u v mnohorozmˇern´ ych datech, ale negarantuj´ı nalezen´ı nejvhodnˇejˇs´ıho modelu, nebot’ ho mohou minout“. Pro hledán´ı vhodného lineárn´ıho regresn´ıho mo” delu je spolehlivˇejˇs´ı procedura popsaná v dalˇs´ı kapitole, ale ta je v´ ypoˇcetnˇe podstatnˇe nároˇcnˇejˇs´ı, takˇze pro velmi rozsáhlá data m˚ uˇze b´ yt jej´ı vyuˇzit´ı problematické.

8.2

Hled´ an´ı nejlepˇ s´ı mnoˇ ziny regresor˚ u

Systematiˇctˇeji neˇz kroková regrese pracuj´ı procedury oznaˇcované jako all possible ” regressions“ nebo best subset of regressors“. Pro kaˇzdé k = 1, . . . , p hledaj´ı takovou ” k-tici regresor˚ u, aby R2 = bylo pro dan´ y poˇcet regresor˚ u maximáln´ı. Jak bylo dˇr´ıve uvedeno, poˇcet model˚ u roste exponenciálnˇe s poˇctem potenciáln´ıch regresor˚ u p, procedury vyuˇz´ıvaj´ıc´ı jen hrubou s´ılu, tj. opravdu zkoumaj´ı vˇsechny modely, mohou b´ yt uˇzity jen pro pomˇernˇe mal´ y poˇcet potenciáln´ıch regresor˚ u p, napˇr. v NCSS 2000 je to p ≤ 15. V nˇekter´ ych statistick´ ych programech je implementována heuristika, která sice nezajiˇst’uje vyˇcerpávaj´ıc´ı prohledán´ı vˇsech model˚ u, ale zato dovoluje vˇetˇs´ı poˇcet regresor˚ u. V NCSS je to procedura Multivariate Selection v Regression – Variable selection routines. Jelikoˇz s rostouc´ım k index determinace R2 neklesá (obvykle roste), nen´ı vhodn´ ym kritériem pro optimalizaci modelu. Vhodnˇejˇs´ım kritériem je adjustovan´ y index determinace 2 =1− Radj

) ( ) n−1 ( k 1 − R2 = R2 − 1 − R2 n−k−1 n−k−1

nebo nejˇcastˇeji uˇz´ıvaná Mallowsova statistika Cp s2 Cp = [n − (k + 1)] k2 − [n − 2(k + 1)] = n s

(

) ( 2 ) s2k sk − 1 − (k + 1) 2 − 2 , s2 s

kde s2k je residuáln´ı rozptyl pˇri k zaˇrazen´ ych regresorech a s2 je residuáln´ı rozptyl pˇri vˇsech zaˇrazen´ ych regresorech. Stˇredn´ı hodnota této statistiky je E(Cp ) = k + 1.

60


Kdyˇz (s2k /s2 ) ≈ 1, tzn. model uˇz nelze podstatnˇe vylepˇsit, pak Cp ≈ k + 1 a zaˇrazen´ım zbyteˇcného dalˇs´ıho regresoru se Cp zvˇetˇs´ı o 1. Tedy vzhledem k Cp je nejlepˇs´ı ten model, kter´ y má Cp nejmenˇs´ı, pˇribliˇznˇe rovné poˇctu zaˇrazen´ ych regresor˚ u zvˇetˇsen´ ych o jedniˇcku. Obvykle je vˇsak Cp jen jedn´ım z kritéri´ı pˇri hledán´ı nejvhodnˇejˇs´ıho modelu, mus´ıme vz´ıt do u ´vahy i residuáln´ı rozptyl a dalˇs´ı, vˇetˇsinou nestatistická kriteria, jako poˇcet regresor˚ u (ˇc´ım ménˇe, t´ım obvykle lépe), cena jejich mˇeˇren´ı (levnˇejˇs´ı má pˇrednost), interpretaci vlivu regresoru na vysvˇetlovanou promˇennou atd. v závislosti na konkrétn´ı u ´loze.


61

Pˇ r´ıklad 8.1 Uˇzit´ı krokové regrese a prohledán´ı vˇsech podmnoˇzin regresor˚ u pˇri hledán´ı vhodného regresn´ıho modelu si ukáˇzeme na datech, která jsou v souboru STEPWISE.XLS. V datech máme 30 pozorován´ı vysvˇetlované veliˇciny y a deseti ´ potenciáln´ıch regresor˚ u, x1, x2, . . . , x10. Ukolem je naj´ıt vhodn´ y lineárn´ı regresn´ı model. Pro v´ ypoˇcty byly uˇzity programy stepwise regression a all subset z [14]. V´ ystupy jsou opˇet uvedeny v surovém stavu, jen s drobn´ ym zkrácen´ım. Stepwise Regression Report Dependent Y Iteration Detail Section Iter. No. 0 1 2 3 4 5 6 7 8 9 10 11 12

Action Unchanged Added Unchanged Added Unchanged Added Unchanged Added Unchanged Added Unchanged Added Unchanged

Variab R-Squared 0.000000 x1 0.765361 0.765361 x5 0.826587 0.826587 x6 0.985724 0.985724 x8 0.988579 0.988579 x2 0.990822 0.990822 x3 0.993080 0.993080

Sqrt(MSE) 3.010984 1.484322 1.484322 1.29947 1.29947 0.3799527 0.3799527 0.3465689 0.3465689 0.3170781 0.3170781 0.2812623 0.2812623

Max R-Sqrd Other X’s 0.000000 0.000000 0.000000 0.300558 0.300558 0.822143 0.822143 0.918500 0.918500 0.978202 0.978202 0.978713 0.978713

Pˇri implicitn´ım nastaven´ı kritéri´ı pro zaˇrazován´ı a vyˇrazován´ı regresor˚ u (α1 = 0.05, α2 = 0.10) byly postupnˇe zaˇrazovány regresory x1, x5, x6, x8, x2 a x3, ˇza´dn´ y nebyl 2 vyˇrazen. Pˇri pohledu na v´ ysledky vid´ıme, podstatná zmˇena v R a residuáln´ı smˇerodatné odchylky nastala po zaˇrazen´ı regresoru x6, tedy pro model se tˇremi regresory x1, x5, x6. Pˇridáván´ı dalˇs´ıch regresor˚ u uˇz index determinace R2 nijak v´ yznamnˇe nezv´ yˇsilo a ani zmenˇsen´ı residuáln´ı smˇerodatné odchylky nen´ı nikterak dramatické. Model se tˇremi regresory x1, x5, x6 je tedy nejnadˇejnˇeˇs´ım kandidátem na model, kter´ y vhodnˇe vysvˇetluje variabilitu veliˇciny y. Zda je to opravdu vhodn´ y model je nutno zkoumat podrobnˇeji s vyuˇzit´ım postup˚ u uveden´ ych v pˇr´ıkladu o odhadu regresn´ıch parametr˚ u a pak i posoudit vˇecné souvislosti s ˇreˇsen´ ym problémem.


62

All Possible Results Section Model Root Size R-Squared MSE

Cp

Model

1 1 1 1 1 1 1 1 1 1

0.765361 0.471363 0.395810 0.377170 0.356763 0.350940 0.347381 0.235375 0.126059 0.065507

1.484322 2.227958 2.381853 2.418317 2.457615 2.468714 2.475473 2.679494 2.864636 2.962214

848.389398 1943.984931 2225.534947 2295.000669 2371.046612 2392.745524 2406.008366 2823.404786 3230.773865 3456.422815

A E C G H I F J D B

(x1) (x5) (x3) (x7) (x8) (x9) (x6) (x10) (x4) (x2)

2 2 2

0.976958 0.826587 0.823812

0.4736817 1.29947 1.309826

61.867264 622.230164 632.570989

BE AE AC

3 3 3

0.985724 0.980175 0.978947

0.3799527 0.4477463 0.4614055

31.201415 51.880208 56.456614

AEF BCE BDE

4 4 4

0.988579 0.988347 0.987601

0.3465689 0.3500752 0.3611051

22.560822 23.426346 26.205962

AEFH ABEF ACEF

5 5 5

0.991966 0.990822 0.990702

0.2966695 0.3170781 0.3191541

11.939729 16.200645 16.649960

ACDEF ABEFH ACEFH

6 6 6

0.993525 0.993080 0.992527

0.2720509 0.2812623 0.29228

8.127863 9.789421 11.849453

ACDEFH ABCEFH ABDEFH

7 ... 8 ... 10

0.994075

0.2660938

8.079170

ABCDEFH

0.994333

0.266362

9.118089

ABCDEFHJ

0.994901

0.2656163

11.000000

ABCDEFGHIJ


63

Ve v´ ystupu z programu All possible si povˇsimnˇeme nejlepˇs´ıho modelu se dvˇema regresory x2 a x5, kter´ y má v´ yraznˇe vyˇsˇs´ı R2 neˇz ostatn´ı modely se dvˇema regresory a pˇribliˇzuje se hodnotám R2 s vˇetˇs´ım poˇctem regresor˚ u. Pˇritom kroková procedura tento model minula“. To je dosti názorná ilustrace nev´ yhod stepwise procedur, ” které jsou sice v´ ypoˇcetnˇe ménˇe nároˇcné neˇz u ´plné prohledáván´ı, ale za cenu rizika takového minut´ı vhodného modelu. Mallowsovo Cp má nejmenˇs´ı hodnotu pro model se sedmi regresory, jen o málo je Cp vˇetˇs´ı pro nejlepˇs´ı model se ˇsesti regresory. Vˇsimnˇeme si, ˇze tento nejlepˇs´ı model se ˇsesti regresory nen´ı shodn´ y s t´ım, kter´ y byl nalezen krokovou procedurou, na m´ısto regresoru x2 je zaˇrazen regresor x4. Vid´ıme, ˇze procedura All possible nám nab´ız´ı v´ıce kandidát˚ u na vhodn´ y model neˇz procedura Stepwise. Mezi tˇemito kandidáty je nutno peˇclivˇe vyb´ırat, rozhodnˇe nen´ı jedin´ y vhodn´ y model s minimáln´ım Cp . Pro v´ ybˇer vhodn´ ych model˚ u jsou uˇziteˇcná i grafická zobrazen´ı statistik pro nalezené modely proti poˇctu regesor˚ u. Jako pˇr´ıklad uvád´ıme grafy pro index determinace a Mallowsovo Cp , na kter´ ych je jasnˇe vidˇet v´ yrazn´ y skok v hodnotách statistik pro nejlepˇs´ı model se dvˇema regresory. Podobnˇe je uˇziteˇcn´ y i graf závislosti residuáln´ı smˇerodatné odchylky.

64


Opravdu vhodn´ y model je vˇsak moˇzno doporuˇcit aˇz po podrobnˇejˇs´ı anal´ yze a porovnán´ı jednotliv´ ych kandidát˚ u. Jak kroková procedura, tak All possible regressions nám jen generuj´ı návrhy, které je nutno podrobnˇeji analyzovat.


65

Shrnut´ı

• výbˇer regresor˚ u do modelu • kroková (stepwise) regrese • hledán´ı nejlepˇs´ı mnoˇziny regresor˚ u • kritéria pro posouzen´ı vhodnosti modelu Kontroln´ı ot´ azky

1. Vysvˇetlete principy a algoritmus krokové regrese. 2. Proˇc maximalizace R2 nen´ı dobrou strategi´ı pˇri hled´ an´ı vhodného modelu? 3. Proˇc minimalize Mallowsovy statistiky je pˇrijatelnou strategi´ı pˇri hled´ an´ı vhodného modelu? 4. Podle ˇceho se posuzuje, který model je vhodný pro danou u ´lohu? 5. Porovnejte výhody a nevýhody krokové regrese a prohled´ av´ an´ı vˇsech model˚ u. Korespondenˇ cn´ı u ´ loha Korespondenˇcn´ı u ´lohy budou zadáv´ any ke kaˇzdému kursu samostatnˇe.

ˇ Í KLASICKEHO ´ ´ ÍHO MODELU 9 ZOBECNEN LINEARN

66

9

Zobecnˇ en´ı klasick´ eho line´ arn´ıho modelu

Pr˚ uvodce studiem V této kapitole jsou ukázány nˇekteré postupy, které rozˇsiˇruj´ı oblast aplikace lineárn´ıho modelu, zejména za okolnost´ı, kdy pˇredpoklady klasického modelu nejsou splnˇeny. Prostudován´ı kapitoly a pochopen´ı souvislost´ı vyˇzaduje nejménˇe tˇri aˇz ˇctyˇri hodiny.

9.1

Transformace p˚ uvodn´ıch regresor˚ u

Prozat´ım jsme se zab´ yvali lineárn´ım modelem, kter´ y obsahoval pˇr´ımo hodnoty regresor˚ u x·,1 , x·,2 , . . . , x·,k , yi = β0 + β1 xi1 + · · · + βk xik + εi

(19)

Jedn´ım z mnoha moˇzn´ ych zobecnˇen´ı je model ve tvaru yi = β0 + β1 Zi,1 + · · · + βp−1 Zi,p−1 + εi ,

(20)

kde kaˇzdé Z·,j je nˇejakou funkc´ı p˚ uvodn´ıch regresor˚ u x·,1 , x·,2 , . . . , x·,k . Jako pˇr´ıklady takov´ ych model˚ u m˚ uˇzeme uvést (pro pˇrehlednost zápisu jsou ˇrádkové indexy vynechány): 1. k = 1, polynom stupnˇe p − 1 y = β0 + β1 x + β2 x2 + · · · + βp−1 xp−1 + ε 2. k = 2, p = 6, tzv. model 2. ˇra´du y = β0 + β1 x1 + β2 x2 + β11 x21 + β22 x22 + β12 x1 x2 + ε 3. k = 2, p = 10, tzv. model 3. ˇra´du y = β0 + β1 x1 + β2 x2 + +β11 x21 + β12 x1 x2 + β22 x22 + +β111 x31 + β112 x21 x2 + β122 x1 x22 + β222 x32 + ε Dalˇs´ı pouˇzitelné transformace jsou • reciproká transformace, tj. zj = 1/xj , kdyˇz ∀xj > 0 • logaritmická transformace, tj. zj = ln(xj ), kdyˇz ∀xj > 0 √ • odmocninová transformace, tj. zj = xj , kdyˇz ∀xj ≥ 0

9.2 Aitken˚ uv odhad

67

a mnoho podobn´ ych dalˇs´ıch transformac´ı a kombinace jejich uˇzit´ı v jednom modelu. D˚ uleˇzité vˇsak je, ˇze modely (20) jsou lineárn´ı v parametrech, tj. soustava normáln´ıch rovnic je soustava p lineárn´ıch rovnic pro p odhadovan´ ych parametr˚ u. Splˇ nuje-li model (20) pˇredpoklady klasického lineárn´ıho modelu, m˚ uˇzeme pro anal´ yzu a interpretaci takového modelu uˇz´ıt vˇsechny techniky, které jsme dosud uˇz´ıvali pro klasick´ y lineárn´ı model ve tvaru (19), tedy pro situaci, kdy jsme pracovali pˇr´ımo s regresory, nikoliv s jejich funkcemi. Na tvar lineárn´ıho modelu je moˇzno nˇekdy pˇrevést i modely, které na prvn´ı pohled lineárn´ı nejsou, napˇr. kdyˇz závislost vysvˇetlované veliˇciny na regresorech x1 , x2 , x3 m˚ uˇze b´ yt proloˇzena funkc´ı η = αxβ1 xγ2 xδ3 (21) Po zlogaritmován´ı dostaneme ln η = ln α + β ln x1 + γ ln x2 + δ ln x3

(22)

Pokud hodnoty vysvˇetlované náhodné veliˇciny y lze popsat modelem ln y = ln η + ε

(23)

a plat´ı, ˇze ε ∼ N (0, σ 2 I), pak k odhad˚ um parametr˚ u α, β, γ, δ a jejich interpretaci opˇet m˚ uˇzeme uˇz´ıt postup˚ u znám´ ych z klasického lineárn´ıho modelu. Pˇri linearizaci vztah˚ u typu (21) vˇsak mus´ıme b´ yt opatrn´ı v tom, jakou roli má náhodné kol´ısán´ı ε (tzv. chybov´ y ˇclen, error). Pˇredstava (23) znamená, ˇze náhodné kol´ısán´ı je multiplikativn´ı, nikoliv aditivn´ı, tj. hodnota náhodné veliˇciny y je vyjádˇrena modelem y = η exp(ε) = αxβ1 xγ2 xδ3 exp(ε), nikoliv modelem s aditivn´ı chybou ve tvaru y = η + error. To, zda je oprávnˇené uˇz´ıt multiplikativn´ı model, m˚ uˇze vyplynout z vˇecné anal´ yzy u ´lohy, ale ˇcasto i v situac´ıch, kdy model chyb je jin´ y, m˚ uˇze b´ yt v´ ysledek b´ yt v´ ysledek z´ıskan´ y linearizac´ı a aplikac´ı lineárn´ıho modelu uˇziteˇcn´ ym prvn´ım pˇribl´ıˇzen´ım k ˇreˇsen´ı problému.

9.2

Aitken˚ uv odhad

Tento odhad 5eˇs´ı problém, kdy nen´ı splnˇen pˇredpoklad (2) klasického model, ˇze náhodné sloˇzky maj´ı konstantn´ı rozptyl a jsou nekorelované. Pˇripust’me, ˇze náhodné sloˇzky mohou b´ yt korelované a nemus´ı m´ıt konstantn´ı rozptyl: cov(ε) = E(εεT ) = σ 2 Ω,

σ 2 > 0,

(24)

68


kde Ω je pozitivnˇe definitn´ı matice. Pak existuje regulárn´ı matice P, pro kterou plat´ı PΩPT = I a PT P = Ω−1

(25)

Vynásob´ıme-li rov.(??), tj. klasick´ y lineárn´ı model, matic´ı P zleva (transformujeme veliˇciny), dostaneme Py = PXβ + Pε (26) Oznaˇc´ıme-li y∗ = Py, X∗ = PX a ε∗ = Pε, pak rov.(26) m˚ uˇzeme pˇrepsat y∗ = X∗ β + ε∗

(27)

Kovarianˇcn´ı matice náhodn´ ych sloˇzek v rov.(27) je pak cov(ε∗ ) = E(ε∗ ε∗T ) = E(PεεT PT ) = σ 2 PΩPT = σ 2 I tzn., ˇze pro hvˇezdiˇckované veliˇciny je rov.(27) klasick´ y lineárn´ı model. Vyjádˇr´ıme rovnice pro odhady parametr˚ u v modelu (27) pomoc´ı p˚ uvodn´ıch netransformovan´ ych veliˇcin a dostaneme vztahy pro odhady parametr˚ u modelu b = (XT Ω−1 X)−1 XT Ω−1 y,

(28)

o kter´ ych v´ıme, ˇze to jsou BLU-odhady s kovarianˇcn´ı matic´ı cov(b) = σ 2 (XT Ω−1 X)−1

(29)

Nestrann´ y odhad parametru σ 2 je s2 = (y − Xb)T Ω−1 (y − Xb)

(30)

kter´ y pak m˚ uˇzeme uˇz´ıt k odhadu kovarianˇcn´ı matice a tedy i rozptyl˚ u odhad˚ u bi . Odhady z´ıskané t´ımto postupem jsou nestranné, nˇekteré jsou dokonce BLU-odhady, avˇsak k jejich v´ ypoˇctu potˇrebujeme znát matici Ω . Tu bohuˇzel v anal´ yze dat v naprosté vˇetˇsinˇe pˇr´ıpad˚ u neznáme. Nem˚ uˇzeme tuto matici z dat ani konsistentnˇe odhadnout, v datech máme n nezávisl´ ych pozorován´ı a potˇrebujeme odhadnout 2 (n + n)/2 jej´ıch prvk˚ u (diagonálu a polovinu nediagonáln´ıch prvk˚ u matice, matice Ω je symetrická). Vˇetˇsinou nezb´ yvá, neˇz na m´ısto pˇredpokladu (24) pˇrijmout nˇejaké vˇetˇs´ı omezen´ı.

9.3 Heteroskedascita

9.3

69

Heteroskedascita

Jedna z moˇznost´ı je ˇreˇsit tzv. heteroskedastickou regresi, tj. pˇripustit, ˇze rozptyly náhodné sloˇzky nejsou konstantn´ı, ale náhodné sloˇzky v lineárn´ım regresn´ım modelu (3) jsou nekorelované: cov(ε) = σ 2 diag(w12 , w22 , · · · , wn2 )

(31)

kde wi2 > 0 je váha rozptylu i-tého pozorován´ı. Pak matice Ω je také diagonáln´ı, Ω = diag(w12 , w22 , · · · , wn2 ), inverzn´ı matice je Ω−1 = diag(w1−2 , w2−2 , · · · , wn−2 ) a P = diag(1/w1 , 1/w2 , · · · , 1/wn ). Pak v modelu (3), tj. v datové matici vydˇel´ıme ˇrádek vahou rovnou smˇerodatné odchylce pozorován´ı yi /wi = β0 /wi + β1 xi1 /wi + · · · + βk xik /wi + εi /wi , m˚ uˇzeme uˇz´ıt OLS-odhady, které budou m´ıt dobré vlastnosti jako v klasickém modelu. Otázkou je, jak urˇcit váhu pozorován´ı, wi . Máme nˇekolik moˇznost´ı, záleˇz´ı na ˇreˇsené u ´loze: (1) V modelu máme je jeden regresor xi1 a pˇredpokládáme, ˇze pozorován´ı závislé veliˇciny yi maj´ı konstantn´ı relativn´ı chybu. Pak m˚ uˇzeme poloˇzit wi = xi1 . (2) Pozorován´ı závislé veliˇciny yi maj´ı konstantn´ı relativn´ı chybu a v modelu je v´ıce regresor˚ u. Pak nezb´ yvá, neˇz vybrat jeden podle subjektivn´ıho rozhodnut´ı, moˇzná ten, kter´ y nejv´ıce koreluje se závisle promˇennou y. (3) Nejdˇr´ıve spoˇc´ıtat yî jako OLS-odhad podle modelu (3) a pak ve druhém kroku poloˇzit wi = yî . (4) Postupovat jako ve variantˇe (9.3) a dále pokraˇcovat v iterac´ıch, dokud dva po sobˇe následuj´ıc´ı odhady nejsou dostateˇcnˇe bl´ızké. Tomuto postupu se ˇr´ıká metoda iterovaných váˇzených ˇctverc˚ u, iterated WLS (Weighted Least Squares).

9.4

Stochastick´ e regresory

Také pˇredpoklad v klasickém modelu, ˇze matice X obsahuje pevné hodnoty, u kter´ ych nen´ı tˇreba uvaˇzovat s jejich rozptylem a korelac´ı, je v mnoha aplikac´ıch nerealistick´ y. Pro tzv. nezávislou stochastickou regresi, kdy pˇredpokládáme, ˇze matice X je stochastická, tvoˇr´ı (k + 1) rozmˇern´ y náhodn´ y proces a náhodná sloˇzka ε nezávis´ı na X, uvedeme struˇcnˇe d˚ uleˇzité v´ ysledky, podrobnˇeji viz napˇr. [7].


70

OLS-odhady b, spoˇc´ıtané podle rov.(10), s2 podle rov.(13) a Sbb = s2 (XT X)−1 jsou nestranné. Ale nejsou to lineárn´ı odhady, nebot’ jsou stochastickou funkc´ı náhodného vektoru y a nejsou to BLU-odhady. Za pˇredpokladu, ˇze v pravdˇepodobnosti konverguje εT ε/n → σ 2 a XT X/n → ΣXX (kde ΣXX je kovarianˇcn´ı matice regresor˚ u) vˇsak plat´ı: • b = (XT X)−1 XT y je konzistentn´ım odhadem vektoru regresn´ıch koeficient˚ u β • s2 podle rov.(13) je konzistentn´ım odhadem parametru σ 2 • Sbb = s2 (XT X)−1 lze vz´ıt za konsistentn´ı odhad asymptotické kovarianˇcn´ı matice odhadovan´ ych parametr˚ u b. To znamená, ˇze OLS-odhady lze uˇz´ıt k bˇeˇzn´ ym test˚ um a urˇcen´ı interval˚ u spolehlivosti pro parametry. Pokud jsou náhodné sloˇzky modelu normálnˇe rozdˇeleny, jsou OLS-odhady podle rov.(10) také ML-odhady, takˇze maj´ı dobré asymptotické vlastnosti, jsou konsistentn´ı a asymptoticky eficientn´ı.

9.5

Diskr´ etn´ı regresory, umˇ el´ e promˇ enn´ e

Dosud jsme se zab´ yvali u ´lohami, ve kter´ ych vysvˇetluj´ıc´ı veliˇciny byly spojité. Docela ˇcasto se v anal´ yze dat stává, ˇze data pocházej´ı ze dvou nebo v´ıce populac´ı, vzpomeˇ nme napˇr. na dvouv´ ybˇerové testy ˇci anal´ yzu rozptylu. I na taková data m˚ uˇzeme aplikovat lineárn´ı regresi. Uvaˇzujme nyn´ı nejjednoduˇsˇs´ı pˇr´ıpad – lineárn´ı regresn´ı model s jedn´ım regresorem EYi = β0 + β1 xi . Parametr β1 je smˇernice pˇr´ımky, tzn. vyjadˇruje zmˇenu stˇredn´ı hodnoty náhodné veliˇciny Yi , zmˇen´ı-li se hodnota regresoru o jedniˇcku. Uvaˇzujme, ˇze regresor x je diskrétn´ı a má hodnoty {0, 1}, jin´ ymi slovy jen rozdˇeluje data do dvou skupin (v´ ybˇer˚ u) ze dvou populac´ı 0 a 1. Pak test hypotézy β1 = 0 znamená totéˇz jako test hypotézy µ0 = µ1 (shoda stˇredn´ıch hodnot obou populac´ı), tj. dvouv´ ybˇerov´ y t-test pˇri shodn´ ych rozptylech.

9.5 Diskrétn´ı regresory, umˇelé promˇenné

71

Pˇ r´ıklad 9.1 Máme otestovat hypotézu, ˇze stˇredn´ı hodnoty veliˇciny Y ze dvou populac´ı jsou shodné. V´ ybˇerové charakteristiky pro oba nezávislé v´ ybˇery jsou v následuj´ıc´ı tabulce a obrázku. v´ ybˇer n pr˚ umˇer sm. odchylka 0 30 5,06 1,04 1 20 5,96 2,09 K testu si m˚ uˇzeme vybrat nˇekolik metod, které nám daj´ı shodné v´ ysledky: metoda H0 t-test(shodné rozptyly) µ0 = µ1 lineárn´ı regrese β1 = 0 ANOVA µ0 = µ1

pˇredpoklad statistika p 2 2 σ0 = σ1 2,01 0,05 σ02 = σ12 2,01 0,05 2 2 σ0 = σ1 4,03 0,05

Jak vid´ıme, ve vˇsech tˇrech pˇr´ıpadech nám vyˇsla stejná hodnota p, pro t-test a lineárn´ı regresi i stejná hodnota statistiky, aˇckoliv testujeme r˚ uzné hypotézy, v anal´ yze rozptylu je hodnota F -statistiky rovna druhé mocninˇe t-statistiky u ostatn´ıch dvou metod. V pˇr´ıpadˇe lineárn´ı regrese je odhad b1 = y¯1 − y¯0 a b0 = y¯0 a testujeme, zda rozd´ıl pr˚ umˇer˚ u je dostateˇcnˇe velik´ y k zam´ıtnut´ı hypotézy µ0 = µ1 (pˇripomeˇ nme, ˇze smˇernice pˇr´ımky je zmˇena veliˇciny y pˇri zmˇenˇe veliˇciny x o jedniˇcku). To, ˇze uvedené statistiky vyˇsly stejnˇe, nen´ı ˇza´dné pˇrekvapen´ı, nebot’ se vyˇc´ısluj´ı ze stejn´ ych formul´ı. Také pˇredpoklady pro vˇsechny uvedené testy jsou shodné, normálnˇe rozdˇelená residua a shodné rozptyly v obou populac´ıch. Uveden´ y pˇr´ıklad ilustruje moˇznost podobného pohledu na anal´ yzu rozptylu a lineárn´ı regresi, ukazuje, ˇze diskrétn´ı regresory mohou b´ yt docela snadno interpretovány a naznaˇcuje smˇery dalˇs´ıho zobecnˇen´ı lineárn´ıho modelu. Pokud diskrétn´ı regresor nab´ yvá v´ıce neˇz dvou hodnot, lze k rozliˇsen´ı uˇz´ıt tzv. umˇelé promˇenné, dummy variables. Obvykle se jedna z r kategori´ı vybere jako referenˇcn´ı a r − 1 dummy promˇenn´ ych s hodnotami {0, 1} pak kóduje kategorie. Odhad smˇernice u konkrétn´ı dummy promˇenné znamená odhad zmˇeny stˇredn´ı hodnoty vysvˇetlované veliˇciny oproti referenˇcn´ı kategorii. Podrobnˇeji viz kapitola Logistická regrese. Pˇri v´ıce diskrétn´ıch regresorech pomocné promˇenné dovoluj´ı zkoumat regresn´ı anal´ yzou i velmi komplikované struktury závislosti, pˇr´ıpadnˇe i v kombinaci s dalˇs´ımi spojit´ ymi regresory tyto závislosti oˇciˇst’ovat“ od vlivu jin´ ych veliˇcin. Podrobnˇejˇs´ı ” v´ yklad takov´ ych postup˚ u pˇresahuje rozsah tohoto kursu, v pˇr´ıpadˇe potˇreby se obrat’te na literaturu, napˇr. Draper a Smith nebo Andˇel atd. Tam najdete i dalˇs´ı moˇznosti zaveden´ı pomocn´ ych promˇenn´ ych.


72

Shrnut´ı

• tranformace regresor˚ u, polynom, model druhého ˇr´ adu, linearizace • heteroskedascita, metoda váˇzených nejmenˇs´ıch ˇctverc˚ u • diskrétn´ı regresory, umˇelé promˇenné Kontroln´ı ot´ azky

1. Jaké zjednoduˇsen´ı pˇredstavuje metoda váˇzených nejmenˇs´ıch ˇctverc˚ u proti Aitkenovu odhadu? 2. Jak interpretovat smˇernici regresn´ı pˇr´ımky v pˇr´ıpadˇe spojitých regresor˚ u a jak v pˇr´ıpadˇe diskrétn´ıch regresor˚ u? 3. Co jsou umˇelé (dummy) promˇenné?

73

10

Zobecnˇ en´ y line´ arn´ı model (GLM)

Pr˚ uvodce studiem Zobecnˇený lineárn´ı model (GLM) projdˇete sp´ıˇse pro celkový pˇrehled, snaˇzte se vˇsak d˚ ukladnˇe pochopit ˇcást o logistické regresi vˇcetnˇe ˇreˇseného pˇr´ıkladu. Na tuto kapitolu poˇc´ıtejte nejménˇe se ˇctyˇrmi hodinami studia s t´ım, ˇze se k prob´ırané látce budete jeˇstˇe vracet po pochopen´ı dalˇs´ıch souvislost´ı. Zobecnˇen´ y lineárn´ı model ( Generalized Linear Model) oznaˇcovan´ y jako GLM nebo u GLIM, je podrobnˇeji popsán v knize McCullagh a Nelder [19] a do základn´ıch pojm˚ tohoto modelu nyn´ı nahlédneme. Zobecnˇen´ y lineárn´ı modle (GLM) zahrnuje: lineárn´ı regresi r˚ uzné modely anal´ yzy rozptylu (ANOVA) logistickou regresi probitov´ y model log-lineárn´ı model (multinomick´ y model pro ˇcetnosti v anal´ yze mnohorozmˇern´ ych kontingenˇcn´ıch tabulek) Oznaˇcen´ı datov´ ych struktur a v´ yznam symbol˚ u v GLM: Pozorován´ı závisle promˇenné (response) je sloupcov´ y vektor náhodn´ ych veliˇcin a je typu (n × 1), tedy y = [y1 , y2 , · · · , yn ]T . Pokud z kontextu je zˇrejmé, ˇze se jedná o libovoln´ y prvek vektoru y, bude oznaˇcován y (netuˇcná kursiva bez indexu) Matice X nezávisl´ ych promˇenn´ ych (regresor˚ u, covariates) je typu (n × p). Jej´ı j-t´ y sloupec oznaˇcujeme xj . Vektor parametr˚ u je β = [β1 , β2 , · · · , βp ]T . Náhodná sloˇzka modelu má vektor stˇredn´ıch hodnot E(Y) = µ typu (n × 1) a kovarianˇcn´ı matici cov(Y) Lineárn´ı prediktor η je systematická sloˇzka v lineárn´ım modelu, tedy η=

p ∑

xj βj

j=1

kde xj je j – t´ y sloupec matice X, tj. vektor (n × 1). Kaˇzdá sloˇzka vektoru Y má rozdˇelen´ı z exponenciáln´ı rodiny rozdˇelen´ı s hustotou fY (y, θ, ϕ) = exp{(yθ − b(θ))/a(θ) + c(y, θ)},

(32)

ˇ Y ´ LINEARN ´ Í MODEL (GLM) 10 ZOBECNEN

74

kde θ a ϕ jsou parametry rozdˇelen´ı, a(.), b(.), c(.) jsou funkce, jejichˇz tvar je dán konkrétn´ım rozdˇelen´ım z exponenciáln´ı rodiny. Pokud ϕ je známé, je rov.(32) hustota rozdˇelen´ı z exponenciáln´ı rodiny a má kanonický parametr θ. Pokud ϕ je neznámé, pak to m˚ uˇze, ale nemus´ı b´ yt dvouparametrické rozdˇelen´ı z exponenciáln´ı rodiny. Napˇr. pro normáln´ı rozdˇelen´ı, Y ∼ N (µ, σ 2 ) fY (y, θ, ϕ) = = exp{(yµ −

√ 1 exp{−(y − µ)2 /2σ 2 } = 2πσ 2 µ2 /2)/σ 2 − 21 (y 2 /σ 2 + ln(2πσ 2 ))},

takˇze v tomto pˇr´ıpadˇe θ = µ,

ϕ = σ2,

1 c(y, ϕ) = − (y 2 /σ 2 + ln(2πσ 2 )) 2 Logaritmus vˇerohodnostn´ı funkce (pˇri známém y funkce parametr˚ u θ, ϕ) je a(ϕ) = ϕ b(θ) = θ2 /2,

l(θ, ϕ, y) = ln fY (y, θ, ϕ) Stˇredn´ı hodnota a rozptyl mohou pak b´ yt urˇceny ze vztah˚ u znám´ ych pro vˇerohodnostn´ı funkci: ( ) ( 2 ) ( )2 ∂l ∂ l ∂l E =0 E + E =0 ∂θ ∂θ2 ∂θ Vˇerohodnostn´ı funkci pro jedno pozorován´ı z jakéhokoli rozdˇelen´ı z exponenciáln´ı rodiny lze zapsat l(θ, ϕ, y) = (yθ − b(θ))/a(ϕ) + c(y, ϕ) Derivace podle kanonického parametru jsou ∂l/∂θ = (y − b′ (θ))/a(ϕ) a ∂ 2 l/∂θ2 = b′′ (θ)/a(ϕ). Poloˇz´ıme-li je rovny nule, m˚ uˇzeme vyjádˇrit stˇredn´ı hodnotu a rozptyl vysvˇetlované náhodné veliˇciny: ( E

∂l ∂θ

) = (µ − b′ (θ))/a(ϕ) = 0 ⇒ E(Y ) = µ = b′ (θ)

a var(Y ) −

b′′ (θ) = 0 ⇒ var(Y ) = b′′ (θ)a(ϕ). a(ϕ)

Stˇredn´ı hodnota je funkc´ı pouze kanonického parametru θ. Rozptyl náhodné veliˇciny Y je souˇcinem dvou funkc´ı. Jedna, b′′ (θ) závis´ı pouze na kanonickém parametru rozdˇelen´ı (a tedy na stˇredn´ı hodnotˇe µ náhodné veliˇciny Y ). Naz´ yvá se varianˇcn´ı funkce (variance function) a m˚ uˇzeme ji zapsat jako funkci stˇredn´ı hodnoty, V (µ). Druhá funkce v souˇcinu je nezávislá na kanonickém parametru θ a závis´ı jen na ϕ.

75

Funkce a(ϕ) má obvykle tvar a(ϕ) = ϕ/w. Parametr ϕ se naz´ yvá disperzn´ı parametr a je konstantn´ı pro vˇsechna pozorován´ı, w je apriornˇe známá váha pozorován´ı, m˚ uˇze b´ yt r˚ uzná pro r˚ uzná pozorován´ı. GLM tedy dovoluje i jiná rozdˇelen´ı z exponenciáln´ı rodiny neˇz jen normáln´ı uˇzité v klasickém modelu. Dalˇs´ı zobecnˇen´ı je v tom, ˇze lineárn´ı prediktor nemus´ı vysvˇetlovat jen (podm´ınˇenou) stˇredn´ı hodnotu náhodné veliˇciny, ale i nˇejakou jej´ı funkci. Vztah mezi lineárn´ım prediktorem η a stˇredn´ı hodnotou µ vysvˇetlované náhodné veliˇciny Y vyjadˇruje spojovac´ı funkce (link): ηi = g(µi ) Spojovac´ı funkce g(.) m˚ uˇze b´ yt jakákoli monotónn´ı diferencovatelná funkce. V klasickém lineárn´ım modelu je spojovac´ı funkc´ı identita, tj. η = µ. U jin´ ych model˚ u se uˇz´ıvaj´ı zejména tyto spojovac´ı funkce: logit probit

η = ln{µ/(1 − µ)} 0<µ<1 −1 η = Φ (µ) 0<µ<1 Φ(.) je distribuˇcn´ı funkce rozdˇelen´ı N(0,1)

komplementárn´ı log-log

η = ln{− { ln(1 − µ)} µλ pro λ ̸= 0 mocninové funkce η = ln µ pro λ = 0

0<µ<1 µ>0

Jelikoˇz stˇredn´ı hodnota µ = b′ (θ), je tedy jen funkc´ı kanonického parametru θ a spojovac´ı funkce je monotónn´ı, existuje inverzn´ı funkce, kterou m˚ uˇzeme vyjádˇrit jako funkci stˇredn´ı hodnoty, θ(µ). Nˇekterá rozdˇelen´ı maj´ı zvláˇstn´ı spojovac´ı funkce, kdy kanonick´ y parametr rozdˇelen´ı je roven lineárn´ımu prediktoru, θ = η. Tyto spojovac´ı funkce se naz´ yvaj´ı kanonické (canonical link). Pro bˇeˇzná rozdˇelen´ı jsou kanonick´ ymi následuj´ıc´ı spojovac´ı funkce: normáln´ı rozdˇelen´ı,N (µ, σ 2 ) Poissonovo, P (µ) alternativn´ı, A(π) gamma, G(µ, v) inversn´ı Gaussovo, IG(µ, σ 2 )

η η η η η

=µ = ln µ = ln{π/(1 − π)} = µ−1 = µ−2


76

V klasickém modelu se v´ ystiˇznost modelu (tˇesnost proloˇzen´ı) vyjadˇruje obvykle pomoc´ı koeficientu determinace R2 , tj. jako pod´ıl variability závislé veliˇciny vysvˇetlené modelem na celkové variabilitˇe. R2 = 1 − ∑

RSS (yi − y¯)2

∑ Celková variabilita (yi − y¯)2 odpov´ıdá RSS lineárn´ıho modelu s jedn´ım parametrem, jehoˇz odhad je b0 = y¯. Pro takov´ y model je ∑ (yi − y¯)2 = 0. R =1− ∑ (yi − y¯)2 2

Pro model vysvˇetluj´ıc´ı variabilitu veliˇciny y u ´plnˇe je RSS = 0 a tedy je R2 = 1. Ve zobecnˇeném lineárn´ım modelu lze model (tˇesnost proloˇzen´ı) posuzovat analogicky. Uvaˇzujme tak zvan´ yu ´pln´ y model s n parametry, kter´ y by vysvˇetloval pozorované hodnoty y pˇresnˇe, tzn. yi = µi . Jelikoˇz m˚ uˇzeme kanonick´ y parametr vyjádˇrit jako funkci stˇredn´ı hodnoty, θ(µ), m˚ uˇzeme vˇerohodnostn´ı funkci zapsat jako l(y, ϕ, y). To je maximálnˇe dosaˇzitelná hodnota vˇerohodnostn´ı funkce. Pro model jen s jedn´ım parametrem, kdy µi = konst (nulov´ y model, obsahuje jen intercept), bychom dostali vˇerohodnostn´ı funkci minimáln´ı hodnoty pro daná data. Dvojnásobek rozd´ılu mezi tˇemito vˇerohodnostn´ımi funkcemi je analogi´ı k celkové variabilitˇe v klasickém modelu. Oznaˇc´ıme-li odhad stˇredn´ıch hodnot v modelu s p ˆ a odhad kanonického parametru pro tento model jako θˆ = θ(µ) a parametry jako µ θ˜ = θ(y) a pˇredpokládáme-li ai (ϕ) = ϕ/wi , pak dvojnásobek rozd´ılu vˇerohodnostn´ıch funkc´ı l(y, ϕ, y) a l(ˆ µ, ϕ, y) je ∑

2wi {yi (θ˜i − θî ) − b(θ˜i ) + b(θî )}/ϕ = D(y, µ ˆ)/ϕ.

D(y, µ ˆ)/ϕ, která je funkc´ı pozorovan´ ych dat, se naz´ yvá deviance a je to analogie residuáln´ı sumy ˇctverc˚ u, RSS. Klasick´ y lineárn´ı model je zvláˇstn´ım pˇr´ıpadem zobecnˇeného modelu, kdy spojovac´ı funkce je identita a pak pro normálnˇe rozdˇelenou náhodnou sloˇzku modelu je deviance rovna residuáln´ımu souˇctu ˇctverc˚ u, ∑ 2 D(y, µ ˆ)/ϕ = (yi − µî ) = RSS D∗ (y, µ ˆ) = D(y, µ ˆ)/ϕ je tzv. scaled deviance, je to deviance vyjádˇrená jako násobek disperzn´ıho parametru. Ve zobecnˇeném lineárn´ım modelu je tedy c´ılem nalézt model, kter´ y zmenˇsuje celkovou devianci (´ umˇernou rozd´ılu logaritm˚ u vˇerohodnostn´ıch funkc´ı mezi u ´pln´ ym modelem a nulov´ ym modelem s jedn´ım parametrem). Takov´ y model m˚ uˇze b´ yt vytváˇren i postupnˇe, mohou b´ yt do modelu zaˇrazovány ty regresory, které nejv´ıce sniˇzuj´ı devianci vzhledem k aktuáln´ımu modelu se zaˇrazen´ ymi k parametry, tedy m˚ uˇze b´ yt pouˇzit krokov´ y (stepwise) postup pro vyhledáván´ı regresn´ıho modelu. Regresory ve

10.1 Logistická regrese

77

zobecnˇeném lineárn´ım modelu mohou b´ yt i kvalitativn´ı (faktory) a regresory mohou b´ yt i interakce (souˇciny) p˚ uvodn´ıch regresor˚ u, takˇze pomoc´ı zobecnˇeného modelu je moˇzné odhadovat parametry i sloˇzit´ ych model˚ u anal´ yzy rozptylu.

10.1

Logistick´ a regrese

K logistickému regresn´ımu modelu dojdeme ze zobecnˇeného lineárn´ıho modelu (GLM) g[E(Y |x)] = xT β, (33) ve kterém nˇejaká funkce g podm´ınˇené stˇredn´ı hodnoty náhodné veliˇciny Y je vyjádˇrena jako lineárn´ı funkce vektoru regresor˚ u xT = (1, x1 , x2 , . . . , xs ) s regresn´ımi koeficienty β T = (β0 , β1 , . . . , βs ). Pokud má náhodná veliˇcina Y alternativn´ı rozdˇelen´ı, tedy Y ∼ A(p), které má, jak známo, stˇredn´ı hodnotu E(Y ) = p, a jako spojovac´ı (t. zv. link) funkci ve zobecnˇeném lineárn´ım modelu zvol´ıme logit, ( logit(p) = ln

p 1−p

) ,

(34)

dojdeme k logistickému regresn´ımu modelu ( ln

p 1−p

) = xT β,

(35)

ve kterém logit podm´ınˇené stˇredn´ı hodnoty je vyjádˇren jako lineárn´ı funkce regresor˚ u. Parametry β0 , β1 , . . . , βs regresn´ıho modelu (35) lze odhadovat metodou maximáln´ı vˇerohodnosti. Algoritmy pro nalezen´ı tˇechto odhad˚ u b0 , b1 , . . . , bs jsou jiˇz ˇradu let implementovány v dostupn´ ych statistick´ ych programech. Logistick´ y regresn´ı model má pomˇernˇe snadnou a pˇr´ımoˇcarou interpretaci. Pomˇer p/(1 − p), tedy pomˇer pravdˇepodobnosti u ´spˇechu“ ku pravdˇepodobnosti ne´ uspˇechu“, je v anglosaském svˇetˇe ” ” oznaˇcován jako odds a je zcela samozˇrejmˇe pouˇz´ıván i mimo statistiku, napˇr. pˇri sázˇ a terminologie nen´ı ustálená, uˇz´ıvá se pomˇer ˇsanc´ı nebo sázkové riziko. kách. Cesk´ Necht’ tedy odds0 =

p0 pˇri hodnotách regresor˚ u x = x0 1 − p0

odds1 =

p1 pˇri hodnotách regresor˚ u x = x1 1 − p1

Pomˇer dvou odds je oznaˇcován jako odds ratio, zkratkou OR.

OR =

odds1 p1 /(1 − p1 ) = odds0 p0 /(1 − p0 )

(36)


78

Odhad regresn´ıho koeficientu bi , i ∈ [1, s], znamená odhad zmˇeny logitu pˇri zmˇenˇe regresoru xi o jedniˇcku a pˇri konstantn´ıch hodnotách regresor˚ u ostatn´ıch, tedy d jestliˇze x1,i − x0,i = 1 a x1,j = x0,j , bi = ln(OR),

j ̸= i,

j = 1, 2, . . . , s

Odhad OR pˇri zmˇenˇe regresoru xi o jedniˇcku lze spoˇc´ıtat jednoduˇse jako d = ebi OR Interpretaci v´ ysledk˚ u logistické regrese ilustruje následuj´ıc´ı pˇr´ıklad nejjednoduˇsˇs´ıho logistického modelu s jedn´ım dichotomick´ ym regresorem. Pro vˇetˇs´ı názornost si pˇredstavme, ˇze regresor X znamená expozici (vystaven´ı riziku), vysvˇetlovaná promˇenná ˇ Y znamená pˇr´ıtomnost pˇr´ıznaku nemoci. Cetnosti pozorovan´ ych pˇr´ıpad˚ u pak m˚ uˇzeme zapsat do ˇctyˇrpoln´ı tabulky Expozice Nemoc X = 1 X = 0 Y =1 a b Y =0 c d Pak, je-li a, b, c, d > 0 odds1 =

a/(a + c) a = , c/(a + c) c d = ad , OR bc

odds0 = ( b1 = ln

b/(b + d) b = d/(b + d) d

ad bc

)

Pro rozptyl tohoto odhadu asymptoticky plat´ı - viz na pˇr. [5] ( ( )) ad var(b1 ) = var ln = 1/a + 1/b + 1/c + 1/d bc Je tedy zˇrejmé, ˇze logistickou regresi je moˇzno aplikovat i v pˇr´ıpadech, kdy regresor je diskrétn´ı dichotomická veliˇcina. Pokud je regresor nomináln´ı, lze takovou promˇennou transformovat na dichotomické veliˇciny s hodnotami {0, 1}, tzv. indikátory (dummy variables). Uvaˇzujme regresor xi , i ∈ [1, s], kter´ y je nomináln´ı s ki kategoriemi a má pozorované hodnoty xli , l = 1, 2, . . . , n, n je poˇcet pozorován´ı. Hodnoty kategori´ı m˚ uˇzeme oznaˇcit ˇc´ıseln´ ymi kódy {0, 1, . . . , ki − 1}. Kategorii s kódem 0 zvol´ıme jako referenˇcn´ı (t.zv. baseline category) a vytvoˇr´ıme ki − 1 indikátor˚ u s ohodnocen´ım podle následuj´ıc´ıho pravidla    1 kdyˇz xli = j (dij )l =

 

0 jinak

j = 1, 2, . . . , ki − 1,

l = 1, 2, . . . , s

(37)


79

Regresn´ı koeficienty koresponduj´ıc´ı s tˇemito indikátory m˚ uˇzeme oznaˇcit βij , j = 1, 2, . . . , ki − 1, i = 1, 2, . . . , s, jejich odhady pak oznaˇc´ıme bij . Odhady regresn´ıch koeficient˚ u u jednotliv´ ych indikátor˚ u jsou vlastnˇe logaritmem odhadovaného pomˇeru odds pˇr´ısluˇsné kategorie k odds kategorie referenˇcn´ı, tedy logaritmem pˇr´ısluˇsného odds ratio. Pro velké v´ ybˇery m˚ uˇzeme 100(1 − α)-procentn´ı oboustrann´ y interval spolehlivosti pro regresn´ı koeficient βij vyjádˇrit jako ⟨bij − u(1 − α/2)SE(bij ),

bij + u(1 − α/2)SE(bij )⟩

a interval spolehlivosti pro OR ⟨exp [bij − u(1 − α/2)SE(bij )] ,

exp [bij + u(1 − α/2)SE(bij )]⟩,

(38)

kde u(1−α/2) je kvantil normovaného normáln´ıho rozdˇelen´ı N (0, 1) a SE(bij ) je smˇerodatná odchylka odhadu regresn´ıho koeficientu. Neobsahuje-li interval spolehlivosti pro OR jedniˇcku, lze odds v této kategorii povaˇzovat za odliˇsn´ y od odds kategorie referenˇcn´ı, takˇze interpretace v´ ysledk˚ u regresn´ıho modelu je velice pˇr´ımoˇcará. Pokud máme regresn´ı model s v´ıce regresory, odhad regresn´ıho parametru vyjadˇruje lineárn´ı závislost predikované veliˇciny na daném regresoru po adjustován´ı vlivu ostatn´ıch regresor˚ u. Tedy v logistické regresi je odhad regresn´ıho koeficientu roven logaritmu odhadovaného odds ratio po adjustaci vlivu ostatn´ıch regresor˚ u. Pˇ r´ıklad 10.1 Data pro tuto u ´lohu jsou v souboru LOGREG2.XLS. Vysvˇetlovaná veliˇcina Y je dichotomická s hodnotami {0, 1}. Hodnota 1 znamená, ˇze pozorovaná osoba je nemocná, hodnotu 0 má osoba zdravá. Regresory jsou veliˇciny expozice (dichotomická, hodnota 1 znamená, ˇze osoba pracuje v rizikovém provozu, hodnota 0 znamená opak), vek (roky) a koureni (poˇcet cigaret za den) jsou spojité, resp. i poˇcet cigaret za spojit´ y m˚ uˇzeme povaˇzovat. Zkrácen´ y v´ ystup z modulu Logistic Regression [14] následuje: Logistic Regression Report Response

Y

Parameter Estimation Section Regression Variable Coefficient Intercept -25.35205 expozice 2.285141 vek 0.6799906 koureni 5.641818E-02

Standard Chi-Square Error Beta=0 5.291554 22.95 0.4990213 20.97 0.1548485 19.28 1.658742E-02 11.57

Prob Level 0.000002 0.000005 0.000011 0.000671

80


Model in Transformation Form -25.35205 + 2.285141*expozice + .6799906*vek + + 5.641818E-02*koureni Note that this is XB. Prob(Y=1) is 1/(1+Exp(-XB)). Odds Ratio Estimation Section Regression Odds Variable Coefficient Ratio Intercept -25.352054 expozice 2.285141 9.827076 vek 0.679991 koureni 0.056418 Model Summary Section Model Model R-Squared D.F. 0.345596 3

Model Chi-Square 77.63

Lower 95% Conf.Limit

Upper 95% Conf.Limit

3.695359

26.133163

Model Prob 0.000000

V prvn´ı ˇca´sti jsou odhady parametr˚ u logistického modelu a statistiky pro test hypotéz o nulovosti parametr˚ u. Vid´ıme, ˇze u vˇsech ˇctyˇr parametr˚ u zam´ıtáme nulovou hypotézu βj = 0, odhadované parametry u vˇsech tˇr´ı regresor˚ u jsou kladné, tzn. logit roste s hodnotou regresoru, je tedy vyˇsˇs´ı u exponovan´ ych, roste s vˇekem a poˇctem vykouˇren´ ych cigaret. V ˇca´sti Odds Ratio Estimation jsou znovu uvedeny odhady pad a 95%-n´ı interval spolehlivosti. rametr˚ u a pro dichotomick´ y regresor je uveden i OR Jelikoˇz tento interval neobsahuje hodnotu 1 (doln´ı hranice intervalu je 3,7), znamená d je v´ to, ˇze OR yznamnˇe vˇetˇs´ı neˇz 1 i po odeˇcten´ı (adjustaci) vlivu vˇeku a kouˇren´ı a ˇze expozice v´ yznamnˇe zvyˇsuje riziko onemocnˇen´ı. Model Summary Section je analogi´ı sekce ANOVA v lineárn´ı regresi a slouˇz´ı k testu hypotézy, ˇze vˇsechny parametry jsou nulové.


Shrnut´ı

• zobecnˇený lineárn´ı model (GLM) • spojovac´ı funkce, kanonické spojovac´ı funkce pro bˇeˇzn´ a rozdˇelen´ı • logistická regrese, odds ratio Kontroln´ı ot´ azky

1. Vysvˇetlete hlavn´ı myˇslenky zobecnˇeného modelu. 2. Co je lineárn´ı prediktor? 3. Jaké rozdˇelen´ı má vysvˇetlovan´ a veliˇcina v logistické regresi? 4. Co je to logit? Je funkc´ı stˇredn´ı hodnoty vysvˇetlované veliˇciny? Korespondenˇ cn´ı u ´ loha Korespondenˇcn´ı u ´lohy budou zadáv´ any ke kaˇzdému kursu samostatnˇe.

81

´ Í REGRESNÍ MODEL 11 NELINEARN

82

11

Neline´ arn´ı regresn´ı model

Pr˚ uvodce studiem Kapitola je vˇenována základ˚ um nelineárn´ı regrese. Na tuto kapitolu poˇc´ıtejte nejménˇe se tˇremi hodinami studia. Prostudujte d˚ ukladnˇe i ˇreˇsený pˇr´ıklad na konci kapitoly. Základn´ı pˇredstava pro nelineárn´ı regresn´ı model je, ˇze stˇredn´ı hodnoty sloˇzek náhodného vektoru y, tj. E(yi ), i = 1, 2, . . . , n, m˚ uˇzeme vyjádˇrit jako nˇejakou funkci regresor˚ u Eyi = f (xi , β), (39) y ˇrádek kde xi je k-ˇclenn´ y vektor nenáhodn´ ych vysvˇetluj´ıc´ıch promˇenn´ ych (xTi je i-t´ matice regresor˚ u X, matice X je typu n × k) a β je p-ˇclenn´ y vektor parametr˚ u. Obvyklou základn´ı u ´lohou nelineárn´ı regrese je pro daná data [y, X] a dan´ y tvar funkce f (xi , β) odhadnout hodnoty parametr˚ u β tak, aby model (39) co nejlépe vysvˇetloval pozorované hodnoty náhodného vektoru y. Zda je model (39) opravdu nelineárn´ı v parametrech poznáme podle parciáln´ıch derivac´ı ∂f (xi , β) gj = ∂βj Pokud gj = const pro vˇsechna j = 1, 2, . . . , p

(40)

tzn. parciáln´ı derivace nejsou závislé na βj , pak model (39) je lineárn´ı v parametrech, pokud alespoˇ n pro jeden z parametr˚ u βj podm´ınka (40) neplat´ı, je model nelineárn´ı. Pokud podm´ınka (40) nen´ı splnˇena pro ˇzádn´ y z parametr˚ u modelu, ˇr´ıkáme, ˇze model je neseparabiln´ı. To je napˇr. následuj´ıc´ı model s jedn´ım regresorem (k = 1) f (xi , β) = exp(β1 xi ) + exp(β2 xi ). Pokud pro nˇekteré parametry je podm´ınka (40) splnˇena, je model separabiln´ı, napˇr. f (xi , β) = β1 + β2 exp(β3 xi ). To je model, kter´ y je nelineárn´ı jen vzhledem k parametru β3 . Nˇekter´ y tvar nelineárn´ıch model˚ u (39) m˚ uˇzeme vhodnou transformac´ı linearizovat, napˇr. f (xi , β) = β1 exp(

β2 ) xi

po zlogaritmován´ı pˇrejde na tvar ln(Eyi ) = ln[f (xi , β)] = γ1 + γ2 zi ,

83

coˇz je lineárn´ı funkce promˇenné zi = 1/xi a parametry γ1 = ln β1 a γ2 = β2 . Hodnoty regresor˚ u z matice X] m˚ uˇzeme zobrazit jako n bod˚ u v (k)-rozmˇerném prostoru. Nelineárn´ı regresn´ı model (39) je plocha v k + 1-rozmˇerném prostoru (srovnej s lineárn´ım modelem s jedn´ım regresorem – hodnoty regresoru jsou zobrazeny na vodorovné ose, tj. v 1D, pˇr´ımka v rovinˇe, tj. v 2D). Na rozd´ıl od lineárn´ıho modelu tuto plochu nelze vyjádˇrit jako lineárn´ı kombinaci regresor˚ u ( nen´ı rovná“), ” má zakˇriven´ı. Pˇri dan´ ych datech a modelové funkci je tvar této plochy závisl´ y na ´ hodnotách parametr˚ u βj . Ulohou odhadu parametr˚ u nelineárn´ıho regresn´ıho modelu je tedy nalézt takové hodnoty parametr˚ u, pro které tato plocha dobˇre aproximuje pozorované hodnoty náhodného vektoru y. Podobnˇe jako u lineárn´ı regrese m˚ uˇzeme tuto u ´lohu ˇreˇsit metodou nejmenˇs´ıch ˇctverc˚ u. Uvaˇzujme tzv. aditivn´ı model pro náhodnou sloˇzku yi = f (xi , β) + εi

(41)

Pˇredpokládejme, ˇze pro i = 1, 2, . . . , n

• εi jsou vzájemnˇe nezávislé náhodné veliˇciny (nejsou korelovány), • E(εi ) = 0, hodnoty yi náhodnˇe kol´ısaj´ı okolo prokládané plochy, stˇredn´ı hodnota tohoto kol´ısán´ı je nulová, • var εi = σ 2 , tzn. maj´ı konstantn´ı rozptyl σ 2 .

Potom souˇcet ˇctverc˚ u rozd´ıl˚ u pozorovan´ ych a modelov´ ych hodnot vyjádˇren´ y jako funkce parametr˚ u je n ∑ Q(β) = [yi − f (xi , β)]2 (42) i=1

Odhad metodou nejmenˇs´ıch ˇctverc˚ u znamená nalézt takové odhady βˆ parametr˚ u β, aby Q(β) bylo minimáln´ı. Zderivujeme-li Q(β) podle βj , j = 1, 2, . . . , p a poloˇz´ıme-li derivace rovny nule, dostaneme soustavu p rovnic n [ ∑ i=1

] ˆ ∂f (xi , β) | ˆ, yi − f (xi , β) β=β ∂βj

(43)

Na rozd´ıl od lineárn´ı regresn´ı funkce, kdy je soustava normáln´ıch rovnic lineárn´ı (parciáln´ı derivace jsou konstantn´ı) a Q(β) eliptick´ y paraboloid s minimem v bodˇe βˆ = b, které lze za dosti obecn´ ych podm´ınek jednoznaˇcnˇe urˇcit, je u nelineárn´ıho modelu soustava normáln´ıch rovnic nelineárn´ı. Jej´ı ˇreˇsen´ı je nároˇcné, nebot’ nemus´ı vˇzdy existovat jednoznaˇcnˇe, nav´ıc je nutno uˇz´ıvat iterativn´ıch metod, které nezaruˇcuj´ı nalezen´ı globáln´ıho minima funkce (42).


84

U nelineárn´ıch model˚ u lze z´ıskat informace o tvaru funkce Q(β) v okol´ı bodu βk z jej´ıho Taylorova rozvoje druhého stupnˇe: 1 Q(β) ∼ = Q(β k ) + ∆β k gk + ∆β Tk Hk ∆β k , 2

(44)

kde ∆β k = β − β k , gk je gradient (vektor) s prvky gj =

∂Q(β) |β=βk ∂βj

a Hk je symetrická matrice ˇrádu p (Hessián) s prvky Hij =

∂ 2 Q(β) |β=βk ∂βi ∂βj

Gradient lze vyjádˇrit gk = −2JT d, kde d je vektor typu (n × 1) s prvky di = yi − f (xi , β k ), J je matice typu (n × p) (Jakobián) s prvky ∂f (xi , β k ) , ∂βj

Jij =

i = 1, 2, . . . , n, j = 1, 2, . . . , p.

Pak pro Hessián plat´ı: Hk = 2JT J + Bk , kde Bk (ˇrádu p) má prvky: Bjl = −2

n ∑

[yi − f (xi , βk )]

i=1

∂ 2 f (xi , β k ) , ∂βj ∂βl

j, l = 1, 2, . . . , p.

Regresn´ı parametry lze odhadnout jednoznaˇcnˇe jen tehdy, jsou-li ∂f (xi , βk ) ∂βj lineárnˇe nezávislé. To znamená, ˇze neexistuj´ı cj ̸= 0, aby p ∑ j=1

cj

∂f (xi , βk ) =0 ∂βj

(45)

Pokud plat´ı (45), tj. JT J je singulárn´ı, pak je model nevhodnˇe specifikován a nelze ˇ ıkáme, ˇze model je pˇreurˇcen. Jediná cesta k nápravˇe je odhadnout jednotlivá βj . R´ zmˇena modelu, která vˇetˇsinou vede pˇres sn´ıˇzen´ı poˇctu parametr˚ u, tj. zjednoduˇsen´ı funkce f (xi , β). Pokud rov. (45) plat´ı pˇribliˇznˇe (analogie s multikolinearitou v lineárn´ım modelu), pak jsou odhady βˆj silnˇe korelované a jejich odhad nen´ı spolehliv´ y.

85

Numerické metody odhadu regresn´ıch parametr˚ u, uˇz´ıvané ve statistickém software, jsou vˇetˇsinou algoritmy, které minimalizuj´ı Q(β) iterativnˇe. Vycházej´ı z poˇcáteˇcn´ıho, ˇ sen´ı prob´ıhá po kroc´ıch uˇzivatelem zadaného nástˇrelu“ hodnot parametr˚ u. Reˇ ” β (0) , β (1) , . . . , β (H) . Zmˇenu mezi jednotliv´ ymi iteraˇcn´ımi kroky m˚ uˇzeme vyjádˇrit jako pˇriˇcten´ı pˇr´ır˚ ustkového vektoru ∆h β (h+1) = β (h) + ∆h , pˇri ˇcemˇz iteraˇcn´ı proces by mˇel splˇ novat podm´ınku: Q(β (h+1) ) < Q(β h ) Pˇr´ır˚ ustkov´ y vektor m˚ uˇzeme vyjádˇrit jako souˇcin tzv. smˇerového vektoru vh a koeficientu αh , ∆h = αh vh Jednotlivé algoritmy se liˇs´ı ve volbˇe smˇerového vektoru vh a zp˚ usobu adaptace koeficientu αh bˇehem v´ ypoˇctu. Z rovnice (44) dostaneme smˇerov´ y vektor ve tvaru vh = −H−1 g = (JT J − B)−1 JT ˆ e,

(46)

kde ˆ e je vektor residu´ı. Optimáln´ı je α = 1. Tato metoda se naz´ yvá Newtonova. Jej´ı nev´ yhodou je, ˇze potˇrebuje v´ ypoˇcet matice druh´ ych derivac´ı kriteriáln´ı funkce Q(β). V metodˇe Gauss-Newtonovˇe se zanedbává matice B a smˇerov´ y vektor se urˇcuje ze vztahu vh = (JT J)−1 JT ˆ e. (47) Dalˇs´ımi bˇeˇznˇe uˇz´ıvan´ ymi metodami jsou r˚ uzné modifikace Marquardtovy metody, kdy smˇerov´ y vektor se urˇcuje podle vztahu e, vh = (JT J + λDTh Dh )−1 JT ˆ

(48)

kde Dh je diagonáln´ı matice eliminuj´ıc´ı vliv r˚ uzn´ ych velikost´ı sloˇzek matice J a λ je parametr, jehoˇz hodnota se adaptuje bˇehem iteraˇcn´ıho procesu. Podrobnˇeji viz napˇr. Meloun a Militk´ y [20]. Jedna z variant této metody je implementována i v NCSS 2000 [14]. Pro vˇsechny iterativn´ı algoritmy je vˇsak velmi podstatn´ y uˇzivatelem zadan´ y nástˇrel“ ” ˇ poˇcáteˇcn´ıch hodnot parametr˚ u a jejich minima a maxima. Spatná volba poˇca´teˇcn´ıch hodnot m˚ uˇze zp˚ usobit bud’ pomalou konvergenci, pˇr´ıpadnˇe ukonˇcen´ı v nˇejakém lokáln´ım minimu, nebo dokonce u ´plné selhán´ı algoritmu. Vˇzdy se vyplat´ı obor hod-

86


not parametr˚ u vˇecnˇe i numericky pˇredem d˚ ukladnˇe analyzovat a teprve pak spustit v´ ypoˇcet. Pokud standardn´ı iterativn´ı algoritmy selhávaj´ı, lze k odhadu parametr˚ u nelineárn´ıho regresn´ıho model uˇz´ıt nˇekter´ y ze stochastick´ ych algoritm˚ u globáln´ı optimalizace. Odhadneme-li parametry nelineárn´ıho regresn´ıho modelu, m˚ uˇzeme (a vˇetˇsinou i mus´ıme, nebot’ správnost odhadu, tj. nalezen´ı globáln´ıho minima funkce Q(β) nen´ı zaruˇcena) posoudit vhodnost modelu a správnost nalezen´ ych odhad˚ u. K prvn´ımu 2 hrubému posouzen´ı poslouˇz´ı index determinace R , R2 = 1 −

RSS , TSS

ˆ celková suma ˇctverc˚ kde RSS = Q(β), u TSS je definována stejnˇe jako u lineárn´ıho modelu. Dalˇs´ımi uˇziteˇcn´ ymi jednoduch´ ymi charakteristikami je odhad residuáln´ıho rozptylu RSS , σˆ2 = s2 = n−p pˇr´ıpadnˇe odhad smˇerodatné odchylky residu´ı, která je odmocninou residuáln´ıho rozptylu. Dále je vˇzdy potˇreba posoudit, zda regresn´ı funkce s nalezen´ ymi hodnotami parametr˚ u dobˇre vystihuje pozorované veliˇciny vysvˇetlované veliˇciny y. Pokud máme v modelu jen jeden regresor, vˇetˇsinou postaˇc´ı jen vizuáln´ı posouzen´ı grafu nalezené funkce a hodnot y proti hodnotám regresoru. Pokud je regresor˚ u v´ıce, m˚ uˇzeme uˇz´ıt grafy residu´ı podobnˇe jako je popsáno v kapitole o lineárn´ım regresn´ım modelu. Je také vhodné posoudit, jak silnˇe jsou odhady parametr˚ u korelovány. Kovarianˇcn´ı matice odhad˚ u se obvykle nejjednoduˇseji aproximuje Sb = s2 (JT J)−1 , po pˇr´ıpadˇe pˇresnˇejˇs´ımi aproximacemi s vyuˇzit´ım i druh´ ych derivac´ı kriteriáln´ı funkce ˆ Q(β) v bodˇe β, tj. s vyuˇzit´ım matice H. V matici Sb jsou na diagonále odhady rozptyl˚ u jednotliv´ ych odhad˚ u parametr˚ u, takˇze lze pak snadno z kovarianˇcn´ı matice Sb i korelaˇcn´ı matici. Pokud se absolutn´ı hodnota nˇekterého z korelaˇcn´ıch koeficient˚ u bl´ıˇz´ı jedné, je model bud’ pˇreurˇcen´ y nebo ˇspatnˇe podm´ınˇen´ y. Pak je potˇreba cel´ y problém znovu analyzovat a bud’ zjednoduˇsit model nebo domˇeˇrit dalˇs´ı data. Za pˇredpokladu, ˇze v modelu (41) maj´ı náhodné sloˇzky normáln´ı rozdˇelen´ı, tj. ε ∼ N (0, σ 2 I), lze pak spoˇc´ıtat i intervalové odhady pro parametry a testovat hypotézy H0: βj = 0, j = 1, 2, . . . , p. Jak intervaly spolehlivosti, tak testy je vˇsak nutno uˇz´ıvat s opatrnost´ı, nebot’ v pˇr´ıpadˇe nelineárn´ı regrese odhady parametr˚ u z´ıskané metodou nejmenˇs´ıch ˇctverc˚ u obecnˇe nemus´ı b´ yt nestranné a intervaly spolehlivosti i

87

testy hypotéz jsou zaloˇzeny jen na aproximaci tvaru funkce Q(β) v okol´ı nalezen´ ych ˆ hodnot odhad˚ u β. Pˇ r´ıklad 11.1 Data v souboru NLR1.XLS obsahuj´ı 44 pozorován´ı veliˇcin X a Y. Regresn´ı funkce má tvar Y = A + (0.49 − A) ∗ exp(−B ∗ X). Máme odhadnout parametry A, B a posoudit vhodnost modelu. V´ ystup z modulu nonlinear regression [14] je následuj´ıc´ı: Nonlinear Regression Report Dependent

Y

Minimization Phase Section Itn Error Sum No. Lambda Lambda A 0 3.58328 0.00004 0.1 1 2.387233E-02 0.000016 0.3856243 Stepsize reduced to 0.9189852 by bounds. Stepsize reduced to 0.9202212 by bounds. Stepsize reduced to 0.9325328 by bounds. 2 1.453589E-02 0.064 0.3846167 3 9.619339E-03 0.0256 0.3635308 4 9.364404E-03 0.01024 0.3440655 5 8.930465E-03 0.004096 0.3211968 6 8.764675E-03 0.0016384 0.3023451 7 8.722906E-03 6.5536E-04 0.2922916 8 8.720914E-03 2.62144E-04 0.2899626 9 8.720909E-03 1.048576E-04 0.2898918 10 8.720909E-03 4.194304E-05 0.2898947 Convergence criterion met.

B 0.13 7.893059E-02

3.805048E-02 3.545998E-02 2.731663E-02 2.272922E-02 2.005132E-02 0.0189681 1.874566E-02 1.874063E-02 0.018741

Ve v´ ystupu vid´ıme, jak z poˇca´teˇcn´ıch hodnot parametr˚ u (A = 0.1, B = 0.13) postupuje iterativn´ı proces hledán´ı minima residuáln´ı sumy ˇctverc˚ u.


88

Model Estimation Section Parameter Name A B

Parameter Estimate 0.2898947 0.018741

Asymptotic Standard Error 6.939709E-02 8.529059E-03

Lower 95% C.L. 0.1498457 1.528661E-03

Upper 95% C.L. 0.4299437 3.595333E-02

Z odhad˚ u parametr˚ u a jejich interval˚ u spolehlivosti vid´ıme, ˇze odhady jsou v´ yznamnˇe odliˇsné od nuly. Model Y = A+(0.49-A)*EXP(-B*X) R-Squared 0.779217 Iterations 10 Estimated Model (.2898947)+(0.49-(.2898947))*EXP(-(.018741)*(X)) Index determinace ukazuje, ˇze model vysvˇetluje asi tˇri ˇctvrtiny z celkové variability veliˇciny Y . Analysis of Variance Table Sum of Source DF Squares Mean 1 7.9475 Model 2 7.978279 Model (Adj) 1 3.077909E-02 Error 42 8.720909E-03 Total (Adj) 43 0.0395 Total 44 7.987

Mean Square 7.9475 3.98914 3.077909E-02 2.076407E-04

Tabulka anal´ yzy rozptylu má podobn´ yu ´ˇcel, jako u lineárn´ıho modelu, m˚ uˇzeme zam´ıtnout hypotézu, ˇze vektor parametr˚ u je nulov´ y. Asymptotic Correlation Matrix of Parameters

A B

A 1.000000 0.996008

B 0.996008 1.000000

Z korelaˇcn´ı matice odhad˚ u je zˇrejmé, ˇze odhady jsou silnˇe korelovány. To je u nelineárn´ıch model˚ u tohoto typu (jeden parametr je v souˇciniteli v´ yrazu, ve kterém je druh´ y parametr v exponentu) ˇcast´ y jev. Indikuje, ˇze minimalizovaná u ´ˇcelová funkce

89

(souˇcet residuáln´ıch ˇctverc˚ u) je v okol´ı nalezeného minima málo zakˇrivená a zmˇena v hodnotách odhad˚ u nezp˚ usobuje dramatickou zmˇenu v hodnotˇe minimalizované funkce. Z diagnostick´ ych graf˚ u vid´ıme, ˇze pˇredpoklady modelu jsou zhruba splnˇeny, rozptyl residu´ı m˚ uˇzeme povaˇzovat za konstantn´ı, residua jsou zhruba normálnˇe rozdˇelená.

90


Nalezen´ y model ukazuje následuj´ıc´ı graf. Pozorované hodnoty jsou vyznaˇceny krouˇzky, odhadované (modelové) hodnoty ˇctvereˇcky.

91

Z grafu je zˇrejmé, ˇze model dobˇre prokládá pozorované hodnoty. M˚ uˇzeme tedy uzavˇr´ıt, ˇze navrˇzen´ y model je pro tuto závislost vhodn´ y, index determinace je 0.78, odhadovaná residuáln´ı odchylka je 0.0144, odhady parametr˚ u jsou v´ yˇse.

Pˇ r´ıklad 11.2 Na tomto pˇr´ıkladu si ukáˇzeme postup pˇri ˇreˇsen´ı jedné praktické u ´lohy, ve které je potˇreba proloˇzit nˇejakou vhodnou funkci namˇeˇren´ ymi daty. Tato u ´loha je pˇrevzata od Ing. Morávky z firmy Tˇrineck´ y inˇzen´ yring, a.s. Data pro tento pˇr´ıklad jsou v souboru moravka prikl.xls. Mˇeˇrila se závislost koncentrace (conc) na ˇcase (t). Zjiˇstˇená empirická závislost je nakreslena na následuj´ıc´ım obrázku:

4.00

conc

3.00

2.00

1.00

0.00 0

50

100 t

150

200

Poˇzadavek byl proloˇzit tuto závislost funkc´ı, která bude m´ıt následuj´ıc´ı vlastnosti: • v ˇcase t = 0 má m´ıt koncentrace hodnotu rovnou nule f (0) = 0 • v ˇcase t → ∞ má m´ıt koncentrace hodnotu rovnou jedné • funkce má dobˇre prokládat empirickou závislost Z grafu závislosti vid´ıme, ˇze poˇca´teˇcn´ı (rostouc´ı) ˇca´st závislosti i za n´ı následuj´ıc´ı prudk´ y pokles bychom mohli vystihnout souˇcinem dvou funkc´ı: • nˇejaké závislosti procházej´ıc´ı poˇcátkem, napˇr. ve tvaru y = A t, kde A je nˇejak´ y neznám´ y parametr, jehoˇz hodnotu pak odhadneme z dat, nebo funkc´ı s dvˇema parametry y1 = A tB , která umoˇzn´ı rychlost r˚ ustu aproximovat pruˇznˇeji.


92

• exponenciáln´ı funkce ve tvaru y2 = exp(C t), která m˚ uˇze dobˇre prokládat klesaj´ıc´ı ˇcást závislosti koncentrace na ˇcase, C je opˇet parametr modelu. (Exponenciáln´ı závislost y2 klesá rychleji neˇz roste y1 ). Pak nám jeˇstˇe zb´ yvá vyˇreˇsit to, aby s rostouc´ım ˇcasem se funkce bl´ıˇzila hodnotˇe jedna. To znamená, ˇze k v´ yˇse uvedenému souˇcinu potˇrebujeme pˇriˇc´ıst nˇejakou funkci y3 , která má malé hodnoty pˇri mal´ ych hodnotách t a lim y3 (t) = 1 pro t → ∞. Tomuto poˇzadavku vyhovuje napˇr. funkce ve tvaru y3 = 1 −

1 exp(D t)

Pro lepˇs´ı pochopen´ı tˇechto u ´vah si zkuste nakreslit pr˚ ubˇehy funkc´ı y1 , y2 , y3 pro r˚ uzné hodnoty jejich parametr˚ u a 0 ≤ t ≤ 100 pomoc´ı nˇejakého software, m˚ uˇzete tˇreba uˇz´ıt moˇznost function plot“ v nab´ıdce Graphics“ v NCSS. ” ” Empirickou závislost tedy m˚ uˇzeme zkusit modelovat funkc´ı f (t) = y1 × y2 + y3 = A tB exp(C t) + 1 −

1 , exp(D t)

kde A, B, C, D jsou ˇctyˇri neznámé parametry, jejichˇz hodnoty odhadneme z dat metodou nejmenˇs´ıch ˇctverc˚ u. K tomu m˚ uˇzeme vyuˇz´ıt standardn´ı statistick´ y software, napˇr. NCSS. V´ıme uˇz, ˇze pro u ´spˇeˇsn´ y odhad parametr˚ u nelineárn´ıho modelu iterativn´ımi metodami je velmi d˚ uleˇzité vhodnˇe volit jejich poˇcáteˇcn´ı nástˇrely“ a dovolen´ y obor ” hodnot, tj. minimum a maximum. Pˇri tom mus´ıme vz´ıt v u ´vahu jak tvar funkce, tak i obor hodnot nezávisle promˇenné, tj. ˇcasu t. U parametru A je to celkem snadná u ´loha: jelikoˇz hodnoty koncentrace jsou nezáporné, i hodnota A mus´ı b´ yt nezáporná. Hodnoty koncentrace jsou menˇs´ı neˇz 4, takˇze maximáln´ı hodnota parametru A nem˚ uˇze b´ yt pˇr´ıliˇs velká, interval [0, 20] je postaˇcuj´ıc´ı, startovac´ı hodnota A = 1 m˚ uˇze b´ yt dobrá volba. U parametru B mus´ıme uváˇzit, jak rychle má funkce r˚ ust. Pokud by B = 1, pak by r˚ ust podle y2 byl lineárn´ı, pro B < 1 pomalejˇs´ı, pro B > 1 rychlejˇs´ı. M˚ uˇzeme zkusit poˇcáteˇcn´ı hodnotu B = 1 a dovolené hodnoty z intervalu [0, 2]. U parametr˚ u C a D mus´ıme b´ yt opatrnˇejˇs´ı. Oba parametry jsou v exponentu, nav´ıc v souˇcinu s veliˇcinou t, která má hodnoty z intervalu [0, 160], tzn. hroz´ı numerické problémy – pˇreteˇcen´ı nejvˇetˇs´ı v poˇc´ıtaˇci reprezentované hodnoty ˇc´ısla v pohyblivé ˇca´rce. Je zˇrejmé, ˇze hodnoty C a D mus´ı b´ yt kladné, ale malé. Tedy zkus´ıme nástˇrely C = 0.01 a D = 0.01 a interval pro oba parametry [0, 0.1] a zjist´ıme, ˇze v pr˚ ubˇehu iterac´ı doˇslo k pˇreteˇcen´ı (overflow). Zmˇen´ıme-li nástˇrel na D = 0.001, iteraˇcn´ı proces konverguje po 25 kroc´ıch a dostaneme následuj´ıc´ı v´ ysledky:

93

Database D:\UZIVATELE\Moravka\moravka_prikl.S0 Dependent conc Model Estimation Section Parameter Name A B C D

Parameter Estimate 1.685585 0.5552467 5.248904E-02 2.545075E-02

Asymptotic Standard Error 9.694253E-02 3.312053E-02 3.513403E-03 7.617775E-03

Lower Upper 95% C.L. 95% C.L. 1.493156 1.878015 0.489503 0.6209905 4.551499E-02 5.946309E-02 1.032959E-02 4.057191E-02

Model conc = (A*T^B)*EXP(-C*T)+1-1/EXP(D*T) R-Squared 0.963863 Iterations 25 Estimated Model ((1.685585)*(T)^(.5552467))*EXP(-(5.248904E-02)*(T)) +1-1/EXP((2.545075E-02)*(T)) Analysis of Variance Table Sum of Mean Source DF Squares Mean 1 378.8356 Model 4 469.3829 Model (Adjusted) 3 90.54733 Error 96 3.394817 Total (Adjusted) 99 93.94215 Total 100 472.7778

Square 378.8356 117.3457 30.18244 3.536267E-02

Asymptotic Correlation Matrix of Parameters

A B C D

A 1.000000 -0.913029 -0.572957 -0.339223

B -0.913029 1.000000 0.810758 0.531945

C -0.572957 0.810758 1.000000 0.885083

D -0.339223 0.531945 0.885083 1.000000

Vid´ıme, ˇze index determinace je pˇribliˇznˇe 0.964, tedy model dobˇre prokládá empirickou závislost. Graf rezidu´ı na následuj´ıc´ım obrázku ukazuje, ˇze vˇetˇs´ı odchylka modelov´ ych hodnot od pozorovan´ ych se vyskytuje pˇredevˇs´ım pro malé hodnoty t,


94

kde i z grafu empirické závislosti je patrné, ˇze v této oblasti byla pˇresnost mˇeˇren´ı nejmenˇs´ı.

Residuals vs Predictor 1.0

Residuals

0.5

0.0

-0.5

-1.0 0.0

50.0

100.0

150.0

200.0

t Vhodnost zvoleného modelu ukazuje i následuj´ıc´ı obrázek, kde jsou kromˇe empirické závislosti nakresleny i modelové hodnoty (body vyznaˇcené troj´ uheln´ıky).

4.00

conc

3.00

2.00

1.00

0.00 0

50

100 t

150

200

95

Shrnut´ı

• nelineárn´ı regresn´ı model, metoda nejmenˇs´ıch ˇctverc˚ u • aproximace tvaru kriteriáln´ı funkce v okol´ı nalezeného minima • metody odhadu parametr˚ u Kontroln´ı ot´ azky

1. Vysvˇetlete hlavn´ı rozd´ıly mezi lineárn´ım a nelineárn´ım regresn´ım modelem. 2. V ˇcem je nalezen´ı odhad˚ u parametr˚ u modelu obt´ıˇzné? 3. Vysvˇetlete principy algoritm˚ u pro odhad parametr˚ u nelineárn´ıch regresn´ıch model˚ u. Korespondenˇ cn´ı u ´ loha Korespondenˇcn´ı u ´lohy budou zadáv´ any ke kaˇzdému kursu samostatnˇe.

ˇ ´ METODY 12 MNOHOROZMERN E

96

12

Mnohorozmˇ ern´ e metody

Pr˚ uvodce studiem Na tuto rozsáhlou kapitolu poˇc´ıtejte nejménˇe s deseti hodinami usilovného studia s t´ım, ˇze se k prob´ırané látce budete jeˇstˇe vracet po pochopen´ı dalˇs´ıch souvislost´ı. Vˇenujte pozornost ˇreˇseným pˇr´ıklad˚ um. Dosud jsme se zab´ yvali regres´ı, kdy jedna náhodná veliˇcina je vysvˇetlována (nebo predikována) pomoc´ı jin´ ych veliˇcin. Hledá se závislost podm´ınˇené stˇredn´ı hodnoty náhodné veliˇciny na regresorech. Je to nejˇcastˇeji aplikovaná statistická metoda. Odhaduje se, ˇze v´ıce jak 90% aplikac´ı statistiky se op´ırá o regresi. Pochopen´ı princip˚ u regrese je velmi uˇziteˇcné pro pochopen´ı ostatn´ıch metod anal´ yzy mnohorozmˇern´ ych dat. Regresn´ı anal´ yza b´ yvá povaˇzována za zcela samostanou ˇca´st stoj´ıc´ı vedle mnohorozmˇern´ ych metod (methods of multivariate analysis). Ve vˇetˇsinˇe statistického software je regresn´ı a korelaˇcn´ı anal´ yza uvádˇena jako samostatná poloˇzka stoj´ıc´ı vedle mnohorozmˇern´ ych metod. Také uˇcebnice a monografie b´ yvaj´ı vˇenovány samostatnˇe regresi a samostatnˇe zb´ yvaj´ıc´ım metodám mnohorozmˇerné anal´ yzy dat. Mezi mnohorozmˇerné metody jsou zaˇrazovány pˇredevˇs´ım: • testy shody vektor˚ u stˇredn´ıch hodnot MANOVA (multivariate analysis of variance) - mnohorozmˇerná analogie anal´ yzy rozptylu • kanonické korelace, které m˚ uˇzeme povaˇzovat za jisté zobecnˇen´ı lineárn´ı regrese, kdy vysvˇetlujeme ne jednu náhodnou veliˇcinu, ale vektor náhodn´ ych veliˇcin • metody klasifikace, kdy pˇredpokládáme, ˇze data pocházej´ı z v´ıce populac´ı a – hledáme pravidlo umoˇzn ˇuj´ıc´ı zaˇradit (klasifikovat) objekt charakterizovan´ y vektorem hodnot do jedné z populac´ı (diskriminaˇcn´ı anal´ yza, logistická regrese, neuronové s´ıtˇe atd.) – pokouˇs´ıme se naj´ıt v datech podmnoˇziny podobn´ ych objekt˚ u (shluková anal´ yza – cluster analysis) • metody redukce dimenze u ´lohy, kdy promˇenlivost a závislosti v datech se pokouˇs´ıme vyjádˇrit pomoc´ı ménˇe veliˇcin. Anal´ yza hlavn´ıch komponent (principal components) vysvˇetluje rozptyl. Faktorová anal´ yza vysvˇetluje kovarianˇcn´ı (korelaˇcn´ı) strukturu.

12.1 Test shody vektoru stˇredn´ıch hodnot

12.1

97

Test shody vektoru stˇ redn´ıch hodnot

Pro test shody vektoru stˇredn´ıch hodnot se uˇz´ıvá Hotteling˚ uv T 2 (ˇcti té-kvadrát) test. Je to mnohorozmˇerná analogie t-test˚ u. Ve struˇcnosti uvedeme základn´ı myˇslenky. Jednovýbˇerový Hotteling˚ uv T 2 test: Testuje se hypotéza, ˇze p-rozmˇern´ y vektor stˇredn´ıch hodnot µ je roven nˇejakému danému konstantn´ımu vektoru. Pˇredpokládá se, ˇze v´ ybˇer je z mnohorozmˇerného normáln´ıho rozdˇelen´ı. Testovou statistikou je pak T 2 = n(¯ x − µ)T S−1 (¯ x − µ).

(49)

Tato statistika má Hottelingovo rozdˇelen´ı. Lze také uˇz´ıt statistiku T2 n − p ∼ Fp,n−p n−1 p

(50)

Intervaly spolehlivosti pro p-rozmˇern´ y vektor stˇredn´ıch hodnot odvod´ıme z (49) a (50). [ 2 ] T n−p P < F1−α (p, n − p) = 1 − α n−1 p Po u ´pravˇe dostaneme (x − µ)T S−1 (x − µ) <

n−1 p F1−α (p, n − p), n n−p

kde (x − µ)T S−1 (x − µ) = c znamená plochu elipsoidu se stˇredem x ¯, jehoˇz tvar a velikost závis´ı na v´ ybˇerové kovarianˇcn´ı matici S. Volbou α urˇc´ıme hodnotu c a m˚ uˇzeme urˇcit intervaly spolehlivosti pro vektor stˇredn´ıch hodnot µ. Dvouvýbˇerový Hotteling˚ uv T 2 test: Testujeme shodu dvou vektor˚ u stˇredn´ıch hodnot (mnohorozmˇerná analogie dvouv´ ybˇerového t-testu). Máme dva v´ ybˇery z p-rozmˇerného normáln´ıho rozdˇelen´ı o rozsaz´ıch n1 , n2 , n1 + n2 = n. Vektory v´ ybˇerov´ ych pr˚ umˇer˚ u jsou x ¯1 , x ¯2 . Za pˇredpokladu shody kovarianˇcn´ıch matic Σ1 = Σ2 = Σ m˚ uˇzeme z v´ ybˇerov´ ych kovarianˇcn´ıch matic S1 , S2 odhadnout spoleˇcnou v´ ybˇerovou kovarianˇcn´ı matic S=

(n1 − 1)S1 + (n2 − 1)S2 n1 + n2 − 2

Oznaˇc´ıme δ = µ1 − µ2 . Pak statistika T2 =

n1 n2 (¯ x1 − x ¯2 − δ)T S−1 (¯ x1 − x ¯2 − δ) n


98

má Hottelingovo rozdˇelen´ı a n − p − 1 T2 ∼ F (p, n − p − 1), p n−2 kterou m˚ uˇzeme uˇz´ıt k testu hypotézy H0 : µ1 = µ2 Pokud Σ1 ̸= Σ2 , jedná se o mnohorozmˇernou analogii dvouv´ ybˇerového t-testu s nestejn´ ymi rozptyly. Pak je test ponˇekud komplikovanˇejˇs´ı, viz napˇr. Hebák a kol.

12.2 Diskriminaˇcn´ı anal´ yza

12.2

99

Diskriminaˇ cn´ı anal´ yza

Diskriminaˇcn´ı anal´ yza je postup, kter´ y hledá vhodné pravidlo (rozhodovac´ı funkci) umoˇzn ˇuj´ıc´ı na základˇe zadan´ ych hodnot vektoru x zaˇradit objekt do nˇekteré, ˇreknˇeme h-té skupiny. Napˇr. velmi jednoduchou u ´lohou tohoto typu je rozhodnout podle zmˇeˇrené teploty osoby o tom, zda je zdravá ˇci nemocná. V tomto pˇr´ıpadˇe x je skalár (teplota) a rozhodovac´ı pravidlo velice jednoduché: je-li teplota vyˇsˇs´ı neˇz 370 C, pak zaˇrad’ do skupiny nemocn´ ych, jinak do skupiny zdrav´ ych. Tedy klasifikujeme osobu, u n´ıˇz pˇr´ısluˇsnost do skupiny neznáme a uˇz´ıváme rozhodovac´ıho pravidla z´ıskaného z dat popisuj´ıc´ıch vztah teploty pˇr´ısluˇsnosti ke skupinˇe. Je jasné, ˇze naˇs´ım zájmem je naj´ıt takové pravidlo, které by klasifikovalo pokud moˇzno správnˇe. V reálném svˇetˇe vˇetˇsinou nen´ı moˇzné naj´ıt pravidlo klasifikuj´ıc´ı správnˇe vˇzdy. Proto dobré pravidlo bude takové, které minimalizuje pravdˇepodobnost chybn´ ych rozhodnut´ı. Jak uvid´ıme za chvilku, za jist´ ych pˇredpoklad˚ u je takov´ ym pravidlem lineárn´ı diskriminaˇcn´ı funkce. Odvozen´ı jej´ıho tvaru si ukáˇzeme pro klasifikaci do dvou skupin. Zavedeme následuj´ıc´ı oznaˇcen´ı: h = 1, 2 – index skupiny Ah – jev pˇr´ısluˇsnost k h-té skupinˇe“ ” P (Ah ) = πh – apriorn´ı pravdˇepodobnost fh (x) – sdruˇzená hustota pro h-tou skupinu P (Ah |x) – aposteriorn´ı pravdˇepodobnost, tj. pravdˇepodobnost pˇr´ısluˇsnosti k hté skupinˇe za podm´ınky dan´ ych hodnot x Hustotu m˚ uˇzeme zapsat fh (x) = f (x|Ah ) pro h = 1, 2, tj. sdruˇzená hustota pro h-tou skupinu je hustota za podm´ınky, ˇze nastane jev Ah . Podle Bayesova vzorce vyjádˇr´ıme aposteriorn´ı pravdˇepodobnost: P (Ah |x) =

P (Ah )fh (x|Ah ) πh fh (x) = , P (A1 )f (x|A1 ) + P (A2 )f (x|A2 ) π1 f1 (x) + π2 f2 (x)

(51)

h = 1, 2. Klasifikovat budeme do skupiny s nejvˇetˇs´ı aposteriorn´ı pravdˇepodobnost´ı. Dále oznaˇcme S – v´ ybˇerov´ y prostor (mnoˇzinu vˇsech moˇzn´ ych v´ ysledk˚ u x). Naˇs´ım c´ılem je rozdˇelit tento v´ ybˇerov´ y prostor na dvˇe ˇcásti splˇ nuj´ıc´ı podm´ınky: S = S1 ∪ S 2 ,

S1 ∩ S2 = ⊘.

Pak kdyˇz x ∈ Sh , zaˇrad´ıme do h - té skupiny.


100

Pravdˇepodobnost chybného zaˇrazen´ı objektu z h-té skupiny do h′ -té skupiny je ∫ P (x ∈ Sh′ |Ah ) = fh (x)dx, h = 1, 2. Sh′

Podle vˇety o u ´plné pravdˇepodobnosti je celková pravdˇepodobnost chybné klasifikace ∫ ∫ ω = π1 f1 (x)dx + π2 f2 (x)dx. (52) S2

S1

Pokud obˇe chyby klasifikace maj´ı stejnou váhu, je optimáln´ı rozhodovac´ı pravidlo, které minimalizuje ω dané vztahem (52). Chceme-li chybám klasifikace dát r˚ uznou váhu, uˇzijeme ztrátovou matici: [ Z=

0 z(2|1) z(1|2) 0

]

Pak celková ztráta z chybné klasifikace je: ∫ ∫ τ = z(2|1)π1 f1 (x)dx + z(1|2)π2 f2 (x)dx S2

S1

a optimáln´ı je postup, kter´ y minimalizuje τ . Objekt ˇrad´ıme do skupiny s vyˇsˇs´ı aposteriorn´ı pravdˇepodobnost´ı, napˇr. z rov. (51) do skupiny 1 zaˇrad´ıme objekt, kdyˇz π1 f1 (x) > π2 f2 (x) (jmenovatel je shodn´ y pro obˇe skupiny). Klasifikaˇcn´ı pravidlo pro zaˇrazen´ı do skupiny 1 je tedy f1 (x) π2 > f2 (x) π1

(53)

Pˇredpokládáme-li p-rozmˇerné normáln´ı rozdˇelen´ı vektoru x, tj. Np (µ1 , Σ1 ) v 1. skupinˇe a Np (µ2 , Σ2 ) ve 2. skupinˇe, pak hustota je: [ ] 1 p fh (x) = (2π)− 2 |Σh |− 2 exp −(x − µh )T Σ−1 h (x − µh )/2 Po dosazen´ı do (53) a zlogaritmován´ı dostaneme xT Γx + η T x + ξ > 0, kde −1 Γ = 0, 5(Σ−1 2 − Σ1 ), −1 η T = µT1 Σ−1 1 − µ2 Σ2

ξ=

) π2 1 ( T −1 1 |Σ2 | µ ln − ln − µ1 Σ1 µ1 − µT2 Σ−1 2 2 2 |Σ1 | π1 2


101

Jsou-li kovarianˇcn´ı matice v obou skupinách shodné, tj. Σ1 = Σ2 , pak odpadne kvadratick´ y ˇclen a rozhodovac´ı pravidlo se podstatnˇe zjednoduˇs´ı: β T x + γ > 0, kde β T = (µ1 − µ2 )T Σ−1 a

1 1 π2 γ = − β T (µ1 + µ2 ) − ln 2 2 π1

Funkce L(x) = β T x

(54)

se naz´ yvá lineárn´ı diskriminaˇcn´ı funkce, zkratkou LDF. Pokud x má p-rozmˇerné normáln´ı rozdˇelen´ı, pak i L(x) má normáln´ı rozdˇelen´ı. ˇ Ctverec Mahalanobisovy vzdálenosti vektor˚ u stˇredn´ıch hodnot je △2 = (µ1 − µ2 )T Σ−1 (µ1 − µ2 ), stˇredn´ı hodnoty LDF jsou pak pro skupinu 1 a skupinu 2 jsou 1 E1 [L(x)] = △2 2

1 E2 [L(x)] = − △2 2

Oba rozptyly jsou shodné var1 [L(x)] = var2 [L(x)] = △2 Oba podprostory S1 a S2 v p-rozmˇerném podprostoru S oddˇeluje nadrovina urˇcená rovnic´ı β T x + γ = 0 ˇcili L(x) = −γ LDF (54) lze vyjádˇrit také jako 1 Lh (x) = µTh Σ−1 x − µTh Σ−1 µh 2 a klasifikovat do té skupiny, pro kterou je Lh (x) nejvˇetˇs´ı. Tak se postupuje, kdyˇz se klasifikuje do v´ıce neˇz dvou skupin. LDF je optimáln´ı rozhodovac´ı pravidlo pro klasifikaci do skupin, pokud náhodn´ y vektor x má normáln´ı rozdˇelen´ı a skupiny se liˇs´ı jen vektorem stˇredn´ıch hodnot, nikoliv kovarianˇcn´ı strukturou. Procedura diskriminaˇcn´ı anal´ yzy z dat, u kter´ ych je klasifikace známa, odhaduje hodnoty parametr˚ u lineárn´ı diskriminaˇcn´ı funkce β. Pak LDF ve tvaru (54) s hod-


102

notami odhad˚ u lze uˇz´ıt pro klasifikaci objekt˚ u, jejichˇz pˇr´ısluˇsnost do skupiny známa nen´ı. Pˇ r´ıklad 12.1 V souboru DISKRIM.XLS jsou na 30 objektech, které pocházej´ı ze dvou populac´ı (veliˇcina skup), zmˇeˇreny hodnoty 10 spojit´ ych veliˇcin (x1 aˇz x10). Naˇs´ım u ´kolem je nalézt pravidlo pro klasifikaci objekt˚ u. Pravidlo má b´ yt co nejjednoduˇsˇs´ı (ˇc´ım ménˇe veliˇcin, t´ım lépe). Pouˇzijeme jednak diskriminaˇcn´ı anal´ yzu (lineárn´ı diskriminaˇcn´ı funkci), jednak logistickou regresi [14]. Abychom ovˇeˇrili, zda v˚ ubec m˚ uˇzeme lineárn´ı diskriminaˇcn´ı funkci uˇz´ıt, je nutné, aby se populace liˇsily ve stˇredn´ıch hodnotách. To lze zjistit pomoc´ı dvouv´ ybˇerového Hotellinova testu: Two-Sample Hotelling’s T2 Report Group

skup

Descriptive Statistics Means Variable 0 1 x1 12.45333 17.25333 x2 14.996 13.638 x3 12.05333 17.23333 x4 183.5267 236.06 x5 180.28 232.2533 x6 187.74 233.14 x7 190.1267 240.2667 x8 188.4933 233.7267 x9 190.2333 238.74 x10 188.08 231.4333 Count 15 15

Standard Deviations 0 1 2.289375 2.207088 2.427947 2.424848 3.346612 2.916129 73.80299 62.96599 37.71711 47.72365 43.81628 48.71579 42.31836 47.44453 37.21266 52.34204 31.271 54.68576 50.21964 61.33117 15 15

Uˇz letm´ y pohled na popisné statistiky v´ yˇse naznaˇcuje, ˇze mezi pr˚ umˇery nˇekter´ ych veliˇcin ve skupinách jsou v´ yznamné rozd´ıly. To potvrzuje jak Hotelling˚ uv test, tak dvouv´ ybˇerové t-testy pro jednotlivé veliˇciny.

Hotelling’s T2 Test Section Covariance Assumption

T2

DF1

DF2

Prob Level


Equal Unequal

101.988 101.988

103

10 10

28 27

0.0002 0.0002

Student’s T-Test Section Prob Variable |Student’s T| Level All (T2) 101.988 0.0002 x1 5.846 0.0000 x2 1.533 0.1366 x3 4.520 0.0001 x4 2.097 0.0451 x5 3.309 0.0026 x6 2.684 0.0121 x7 3.055 0.0049 x8 2.728 0.0109 x9 2.982 0.0059 x10 2.118 0.0432 These individual t-test significance levels should only be used when the overall T2 value is significant.

Stepwise procedura diskriminaˇcn´ı anal´ yzy poskytne následuj´ıc´ı v´ ystup: Discriminant Analysis Report Dependent

skup

Variable-Selection Summary Section Action Independent Pct Chg In Prob Iteration This Step Variable Lambda F-Value Level 0 None 1 Entered x1 54.97 34.18 0.000003 2 Entered x3 25.24 9.12 0.005481 Variable-Selection Detail Section - Step 2 Independent Pct Chg In Prob Status Variable Lambda F-Value Level In x1 41.77 19.37 0.000152 In x3 25.24 9.12 0.005481

R-Squared Other X’s 0.226687 0.226687


104

Out x2 4.74 1.29 Out x4 6.85 1.91 Out x5 0.17 0.04 Out x6 6.30 1.75 Out x7 1.46 0.39 Out x8 6.35 1.76 Out x9 0.33 0.09 Out x10 2.25 0.60 Overall Wilks’ Lambda = 0.336670 Action this step: None

0.265589 0.178413 0.836049 0.197650 0.539842 0.195840 0.772237 0.445960

0.118719 0.749175 0.398169 0.505648 0.485859 0.505363 0.485821 0.320911

Linear Discriminant Functions skup Variable Constant x1 x3

0 -22.93688 2.481297 1.242258

1 -44.95984 3.438486 1.775299

Classification Count Table for skup Predicted Actual 0 1 Total 0 15 0 15 1 1 14 15 Total 16 14 30

Jako veliˇciny v´ yznamnˇe odliˇsuj´ıc´ı dvˇe skupiny byly do klasifikaˇcn´ıho pravidla krokovou procedurou vybrány x1 a x3. Ze 30 objekt˚ u je pak 29 klasifikováno správnˇe a jen jeden chybnˇe. Toto empirické ovˇeˇren´ı spolehlivosti klasifikace je vˇsak nutno brát s opatrnost´ı, nebot’ spolehlivost klasifikaˇcn´ıho pravidla je ovˇeˇrována na datech, ze kter´ ych byly koeficienty lineárn´ı diskriminaˇcn´ı funkce spoˇc´ıtány, nikoliv na nezávisl´ ych pozorován´ıch, kde mus´ıme poˇc´ıtat s niˇzˇs´ı u ´spˇeˇsnost´ı. Realistiˇctˇejˇs´ı odhad oˇcekávané u ´spˇeˇsnosti klasifikace lze z´ıskat bud’ tak, ˇze data rozdˇel´ıme náhodnˇe na skupinu uˇc´ıc´ı a testovac´ı, parametry lineárn´ı diskriminaˇcn´ı funkce se spoˇc´ıtaj´ı jen z uˇc´ıc´ı skupiny a u ´spˇeˇsnost klasifikace se odhadne ze zjiˇstˇené klasifikace testovac´ı skupiny. Tento postup má ovˇsem nev´ yhodu, ˇze parametry lineárn´ı diskriminaˇcn´ı funkce se poˇc´ıtaj´ı z podstatnˇe menˇs´ıho poˇctu pozorován´ı, tzn. nevyuˇzije se informace v datech. V nˇekter´ ych programech (v NCSS prozat´ım nikoliv) je proto jackknife procedura, která postupnˇe spoˇc´ıtá parametry lineárn´ı diskriminaˇcn´ı funkce z n − 1


105

objekt˚ uau ´spˇeˇsnost klasifikace se vˇzdy ovˇeˇruje na objektu, kter´ y byl vyjmut. Tak se z´ıská odhad spolehlivosti klasifikace na nezávisl´ ych pozorován´ıch. Na obrázku, kter´ y následuje, vid´ıme nesprávnˇe klasifikovan´ y objekt 15 (ze skupiny 1), kter´ y leˇz´ı uvnitˇr shluku objekt˚ u skupiny 0 ( kaz´ı“ lineárn´ı separabilitu). Pro ” ostatn´ı body v rovinˇe lze skupiny od sebe oddˇelit lineárn´ı funkc´ı (pˇr´ımkou).

Stejnou u ´lohu hledán´ı klasifikaˇcn´ıho pravidla pro klasifikaci do dvou skupin lze ˇreˇsit i logistickou regres´ı. Pˇr´ısluˇsnost do skupiny je nutno oznaˇcit { 0, 1 }. Klasifikace je pak zaloˇzena na odhadu pravdˇepodobnosti, ˇze pro dané hodnoty regresor˚ u má veliˇcina Y má hodnotu 1. Tvar klasifikaˇcn´ı funkce lze snadno vyjádˇrit z modelu logistické regrese (35) exp(xT β) p= 1 + exp(xT β) Je-li p vˇetˇs´ı neˇz zvolená hodnota (vˇetˇsinou 0,5), pak objekt klasifikujeme do skupiny 1, jinak do skupiny 0. Klasifikaˇcn´ı pravidlo na rozd´ıl od lineárn´ı diskriminaˇcn´ı funkce je sloˇzitˇejˇs´ı, pˇredevˇs´ım nen´ı lineárn´ı funkc´ı regresor˚ u. To v nˇekter´ ych pˇr´ıpadech je v´ yhodou, nebot’ lze naj´ıt vhodné klasifikaˇcn´ı pravidlo i pro skupiny, které nejsou lineárnˇe separabiln´ı nebo v situac´ıch, kdy nejsou splnˇeny pomˇernˇe pˇr´ısné pˇredpoklady pro aplikaci lineárn´ı diskriminaˇcn´ı funkce (mnohorozmˇerné normáln´ı rozdˇelen´ı, shoda kovarianˇcn´ıch matic ve skupinách). Nˇekdy je ovˇsem tato v´ yhoda problematická, jak ukazuje následuj´ıc´ı v´ ystup z postupné logistické regrese (metoda


106

forward, tj. postupné pˇridáván´ı v´ yznamn´ ych regresor˚ u bez vyluˇcován´ı nadbyteˇcn´ ych). Pˇ r´ıklad 12.2 Logistic Regression Report Response

skup

Forward Variable-Selection Action Variable Added x1 Added x3 Added x5 Added x2 Parameter Estimation Section Regression Standard Chi-Square Variable Coefficient Error Beta=0 Intercept -231.4874 68355.56 0.00 x1 4.703548 2773.934 0.00 x3 4.343548 1478.741 0.00 x5 2.484694 322.7473 0.00 x2 -27.89017 5084.675 0.00 Model Summary Section Model Model Model R-Squared D.F. Chi-Square 0.624562 4 41.59

Prob Level 0.997 0.998 0.997 0.993 0.995

Model Prob 0.000000

Classification Table Predicted Actual 0 1 Total 0 Count 15 0 15 1 Count 0 15 15 Total Count 15 15 30 Percent Correctly Classified=100 Logistickou regres´ı se podaˇrilo nalézt pravidlo se ˇctyˇrmi regresory x1, x3, x5 a x2, které má stoprocentn´ı u ´spˇeˇsnost klasifikace. Lineárn´ı diskrimininaˇcn´ı funkce pro tato data stoprocentn´ı u ´spˇeˇsnosti klasifikace nedosáhne ani pˇri zaˇrazen´ı vˇsech deseti veliˇcin. Ale pˇri podrobnˇejˇs´ım pohledu na odhady parametr˚ u logistického modelu a


107

statistiky s nimi spojené vid´ıme, ˇze smˇerodatné odchylky odhad˚ u parametr˚ u jsou velmi vysoké v porovnán´ı s hodnotami odhad˚ u, takˇze vlastnˇe ˇza´dn´ y z odhad˚ u parametr˚ u nem˚ uˇzeme povaˇzovat za v´ yznamnˇe odliˇsn´ y od nuly. Klasifikaˇcn´ı pravidlo je uˇsito na m´ıru“ dat˚ um. Pokud zpˇr´ısn´ıme kritérium pro zaˇrazován´ı regresor˚ u, dosta” neme následuj´ıc´ı v´ ystup: Logistic Regression Report Forward Variable-Selection Action Variable Added x1 Added x3 Parameter Estimation Section Regression Variable Coefficient Intercept -26.49609 x1 1.113222 x3 0.7096213 Model Summary Section Model Model R-Squared D.F. 0.540694 2

Standard Error 10.43108 0.5043868 0.3449301

Model Chi-Square 31.78

Chi-Square Prob Beta=0 Level 6.45 0.011082 4.87 0.027309 4.23 0.039658

Model Prob 0.000000

Classification Table Predicted Actual 0 1 Total 0 Count 15 0 15 1 Count 1 14 15 Total Count 16 14 30 Percent Correctly Classified=96.67 Toto klasifikaˇcn´ı pravidlo obsahuje dva regesory (stejné jako LDF), vˇsechny parametry tohoto logistického modelu m˚ uˇzeme povaˇzovat za nenulové a klasifikaˇcn´ı pravidlo má stejnou u ´spˇeˇsnost jakou mˇela LDF.


108

12.3

Shlukov´ a anal´ yza

C´ılem shlukové anal´ yzy je nalézt v datech podmnoˇziny podobn´ ych objekt˚ u. Mˇejme mnoˇzinu m objekt˚ u, tuto mnoˇzinu oznaˇcme M . Pro kaˇzdé dva objekty a, b ∈ M máme ˇc´ıslo σ(a, b), kterému ˇr´ıkáme numerická podobnost. σ :M ×M →R Poˇzadavky na vlastnosti numerické podobnosti: 1. 0 ≤ σ(a, b) ≤ 1 2. σ(a, a) = 1 3. σ(a, b) = σ(b, a) 4. min(σ(a, b), σ(b, c)) ≤ σ(a, c) – slabˇs´ı troj´ uheln´ıková nerovnost Charakteristiku (1 − σ(a, b)) m˚ uˇzeme chápat jako normovanou vzdálenost dvou objekt˚ u a, b. ´ Ulohou shlukové anal´ yzy je naj´ıt rozklad {Mi }ki=1 , mnoˇziny M tj. 1.

k ∪

Mi = M

i=1

2. Mi ∩ Mj = 0 pro i ̸= j 3. vágn´ı kritérium: objekty uvnitˇr Mi jsou si podobnˇejˇs´ı mezi sebou neˇz s objekty z mnoˇziny Mj , napˇr. kdyˇz a, b ∈ Mi , c ∈ Mj , pak σ(a, b) ≤ σ(a, c), σ(a, b) ≤ σ(b, c) Je mnoho moˇznost´ı, – jak definovat numerickou podobnost, – jak formulovat postup zaˇrazován´ı objekt˚ u do podmnoˇzin, tedy existuje mnoho metod shlukován´ı. 12.3.1

Hierarchick´ e metody

Vycház´ı se z matice podobnosti objekt˚ u (symetrická matice s jedniˇckami na diagonále), nejˇcastˇejˇs´ı je aglomerativn´ı procedura, zaˇcne od m shluk˚ u (kaˇzd´ y shluk je tvoˇren jedn´ım objektem a spojuje ty shluky, které jsou si nejpodobnˇejˇs´ı, aˇz skonˇc´ı jedn´ım shlukem, obsahuj´ıc´ım vˇsech m objekt˚ u. Pro takovou posloupnost rozklad˚ u {Mij }kj=1 , pro i1 < i2 plat´ı Mi1 ,j ⊆ Mi2 ,j , tj. rozklady jsou do sebe zasunuty, objekty

12.3 Shluková anal´ yza

109

jednou spojené do shluku z˚ ustávaj´ı spolu. Posloupnost spojován´ı m˚ uˇzeme je graficky znázornit dendrogramem. Podobnou u ´lohu ˇreˇs´ı taxonomie v biologii. Nejˇcastˇeji uˇz´ıvané strategie spojován´ı shluk˚ u jsou:

– single linkage (nejbliˇzˇs´ı soused, nearest neiborough) - shluk tvoˇr´ı souvisl´ y podgraf, tj. existuje aspoˇ n jedna cesta mezi dvˇema uzly podgrafu, nejménˇe pˇr´ısná metoda na podobnost uvnitˇr shluk˚ u, shluky maj´ı tvar souhvˇezd´ı“ ” – complete linkage (nejvzdálenˇejˇs´ı soused, furthest neiborough) shluk tvoˇr´ı u ´pln´ y podgraf, tj. kaˇzdé dva uzly podgrafu jsou spojeny hranou, nejpˇr´ısnˇejˇs´ı na podobnost uvnitˇr shluku – average linkage - spojuje shluky podle jejich pr˚ umˇerné vzdálenosti – centroidn´ı - spojuje shluky podle vzdálenost´ı jejich tˇeˇziˇstˇe

Rozd´ıly mezi strategiemi shlukován´ı ilustruj´ı následuj´ıc´ı pˇr´ıklady.

Pˇ r´ıklad 12.3 Data pro tento pˇr´ıklad jsou z knihy [18] a jsou uvedena i v souboru EMPLOY.XLS. Obsahuj´ı u ´daje o pod´ılu zamˇestnan´ ych v dev´ıti odvˇetv´ıch ve 26 ´ evropsk´ ych zem´ı. Udaje jsou z konce 70. let 20. stolet´ı, proto jsou v nich uvedeny i státy, které uˇz nyn´ı neexistuj´ı. Jednotlivé veliˇciny znamenaj´ı: AGR = agriculture (zemˇedˇelstv´ı), MIN = mining (tˇeˇzba), MAN = manufacturing (tˇeˇzk´ y pr˚ umysl), PS = power supplies (energetika), CON = construction (stavebnictv´ı), SER = service industries (lehk´ y pr˚ umysl), FIN = finance, SPS = social and personal services (sociáln´ı sluˇzby), TC = transport and communications (doprava a spoje). Ve vˇsech ukázkách je zvolena Eukleidovská vzdálenost mezi objekty, vˇsechny veliˇciny byly standardizovány, po standardizaci tedy maj´ı jednotkov´ y rozptyl. V´ ystupy se liˇs´ı jen podle pouˇzité strategie (metody) shlukován´ı.

Hierarchical Clustering Report Variables AGR to TC Clustering Method Single Linkage (Nearest Neighbor) Distance Type Euclidean Scale Type Standard Deviation

110


Na dendrogramu vid´ıme, jak postupnˇe byly vytváˇreny shluky. Jako prvn´ı se spojily ˇ edsko, pak Belgie s Franci´ı atd. Dendrogram nejpodobnˇejˇs´ı objekty, tj. Dánsko a Sv´ ukazuje typické chován´ı metody nejbliˇzˇs´ıho souseda, kdy k velk´ ym shluk˚ um jsou pˇripojovány jednotlivé objekty nebo shluky s mal´ ym poˇctem objekt˚ u. Na u ´rovni nepodobnosti (vzdálenosti) zhruba 0,85 jsou zemˇe rozdˇeleny do 6 shluk˚ u: ˇ 1. západoevropské zemˇe s v´ yjimkou Spanˇ elska a Lucemburska 2. Lucembursko 3. zemˇe b´ yvalého socialistického bloku ˇ 4. Spanˇ elsko 5. b´ yvalá Jugoslávie 6. Turecko

Pˇ r´ıklad 12.4 Hierarchical Clustering Report


Variables Clustering Method Distance Type Scale Type

111

AGR to TC Complete Linkage (Furthest Neighbor) Euclidean Standard Deviation

ˇ edsko, pak Belgie Jako prvn´ı se opˇet spojily nejpodobnˇejˇs´ı objekty, tj. Dánsko a Sv´ s Franci´ı atd. Metoda nejvzdálenˇejˇs´ıho souseda má tendenci vytváˇret kompaktnˇejˇs´ı shluky, ve kter´ ych je poˇcet objekt˚ u rovnomˇernˇejˇs´ı. Na u ´rovni nepodobnosti (vzdálenosti) zhruba 0,95 jsou zemˇe rozdˇeleny do 6 shluk˚ u: 1. západoevropské Belgie aˇz Finsko, seznam viz dendrogram ˇ ycarsko, Itálie, Lucembursko a Spanˇ ˇ 2. SRN, Sv´ elsko ˇ 3. Recko, Portugalsko a pˇr´ımoˇrské zemˇe b´ yvalého socialistického bloku s v´ yjimkou NDR ˇ 4. Ceskoslovensko, NDR a Mad’arsko 5. Turecko 6. b´ yvalá Jugoslávie

112


Pˇ r´ıklad 12.5 Hierarchical Clustering Report Variables AGR to TC Clustering Method Simple Average (Weighted Pair-Group) Distance Type Euclidean Scale Type Standard Deviation

Metoda pr˚ umˇerné vzdálenosti vedla k v´ ysledk˚ um, které jsou podobné metodˇe nejvzdálenˇejˇs´ıho souseda, rozd´ıly jsou pˇredevˇs´ım uprostˇred shlukovac´ı procedury, napˇr. ˇ Spanˇ elsko se ke shluku západoevropsk´ ych zem´ım pˇripojilo pozdˇeji neˇz SSSR ke shluku Rumunska, Polska atd.


12.3.2

113

Nehierarchick´ e metody

Mezi nejpopulárnˇejˇs´ı nehierarchické metody patˇr´ı metoda k-means . Poˇcet shluk˚ u k je pˇredem znám, objekty se rozdˇeluj´ı do shluk˚ u tak, aby rozptyl uvnitˇr shluk˚ u (within sum of squares) byl co nejmenˇs´ı. Jde tedy o to, abychom nalezli takové pˇriˇrazen´ı objekt˚ u do shluk˚ u tak, aby stopa matice W byla minimáln´ı. W=

k ∑

Wg ,

(55)

g=1

Wg je Wishartova matice pro shluk g, tj. Wg =

ng ∑

(g)

(g)

(xj − x ¯(g) )(xj − x ¯(g) )T ,

(56)

j=1 (g)

kde xj je vektor hodnot veliˇcin j-tého objektu v g-tém shluku, x ¯(g) = (∑ ) ng (g) /ng je vektor pr˚ umˇer˚ u (centroid) g-tého shluku. Kritériem, jeˇz má b´ yt j=1 xj minimalizovano, je pak TRW = tr(W). (57) Naj´ıt globáln´ı minimum je algoritmicky obt´ıˇzn´ y problém, kter´ y neum´ıme vyˇreˇsit v polynomiáln´ım ˇcase. Obvykle se uˇz´ıvá se Hartigan˚ uv algoritmus k-means, kter´ y um´ı naj´ıt pˇrijatelné lokáln´ı minimum pro vˇetˇsinu jednoduˇsˇs´ıch klasifikaˇcn´ıch u ´loh nebo se v posledn´ı dobˇe pro optimalizaci klasifikace vyuˇz´ıvaj´ı evoluˇcn´ı algoritmy. Algoritmus k-means je velmi jednoduch´ y: 1. Nejdˇr´ıve se k centroid˚ u (tˇeˇziˇst’ shluk˚ u) zvol´ı náhodnˇe, bud’ se vybere náhodnˇe k objekt˚ u ze zadan´ ych dat nebo se objekty náhodnˇe klasifikuj´ı do k shluk˚ ua spoˇc´ıtaj´ı jejich tˇeˇziˇstˇe (vektor pr˚ umˇer˚ u). 2. Objekty se zaˇrad´ı do shluku, jehoˇz tˇeˇziˇsti jsou nejbliˇzˇs´ı a spoˇc´ıtá se nové tˇeˇziˇstˇe kaˇzdého shluku. 3. Krok 2 se opakuje tak dlouho, dokud docház´ı ke zmˇenˇe klasifikace objekt˚ u.


114

Pˇ r´ıklad 12.6 Vyuˇzijeme opˇet data ze souboru EMPLOY.XLS. Struˇcné v´ ysledky pro ˇsest shluk˚ u následuj´ı. K-Means Cluster Analysis Report Iteration Section Iteration No. of No. Clusters 1 2 2 3 3 4 4 5 5 6

Percent of Variation 72.59 52.48 43.44 36.60 33.70

Bar Chart of Percent |||||||||||||||||||||| |||||||||||||||| |||||||||||||| ||||||||||| |||||||||||

Vid´ıme, jak s poˇctem shluk˚ u klesá pod´ıl variability uvnitˇr shluk˚ u (within sum of squares) na celkové variabilitˇe. Nejv´ yraznˇejˇs´ı skok je mezi 2 a 3 shluky,pak uˇz se rychlost sniˇzován´ı zmenˇsuje. Pr˚ umˇery (tj. souˇradnice tˇeˇziˇstˇe shluk˚ u) jsou v následuj´ıc´ı tabulce. Podobnou tabulku volitelnˇe obsahuje v´ ystup z procedury k-means [14] i pro smˇerodatné odchylky uvnitˇr shluk˚ u. Cluster Means Variab Clust1 AGR 31.7 MIN 0.95 MAN 25.17 PS 0.62 CON 9.17 SER 10.1 FIN 3.72 SPS 12.8 TC 5.72 Count 4

Clust2 12.9 0.97 26.75 1.35 7.7 16.3 4.72 22.55 6.77 4

Clust3 20.13 2.45 31.68 1.08 8.33 8.56 0.85 19.18 7.71 6

Clust4 9.76 1.20 31.9 0.78 8.98 17.06 4.50 19.92 5.88 5

Clust5 6.78 0.4 24.04 0.82 8.44 16.6 6.04 29.46 7.46 5

Clust6 57.75 1.1 12.35 0.6 3.85 5.8 6.2 8.6 3.6 2

Zemˇe byly do shluk˚ u zaˇrazeny takto: ˇ ˇ Portugalsko, Spanˇ elsko, Rumunsko 1. Recko, 2. Irsko, Británie, Rakousko, Finsko ˇ 3. Bulharsko, Ceskoslovensko, NDR, Mad’arsko, Polsko, SSSR ˇ ycarsko 4. Francie, SNR, Itálie, Lucembursko, Sv´


115

ˇ edsko 5. Belgie, Dánsko, Holandsko, Norsko, Sv´ 6. Turecko, Jugoslávie V´ ysledky se s t´ım, co poskytly hierarchické procedury, shoduj´ı jen ˇca´steˇcnˇe, ovˇsem podstatné rysy moˇzné klasifikace jsou spoleˇcné. Právˇe porovnán´ı v´ ysledk˚ u v´ıce shlukovac´ıch procedur a nalezen´ı jejich spoleˇcn´ ych rys˚ u je uˇziteˇcné pro u ´vahy o moˇzné klasifikaci.

K této u ´loze se jeˇstˇe vrát´ıme v kapitole o hlavn´ıch komponentách.


116

12.4

Anal´ yza hlavn´ıch komponent

Anal´ yza hlavn´ıch komponent (Principal components analysis, PCA) je jedna z metod redukce dimenze u ´lohy. Snaˇz´ı se vysvˇetlit celkov´ y rozptyl vektoru náhodn´ ych veliˇcin, resp. jeho podstatnou ˇcást pomoc´ı ménˇe veliˇcin. x = (X1 , X2 , . . . , Xp )T náhodn´ y vektor V = [cov(Xi , Xj )] = [σij ], i, j = 1, 2, . . . , p kovarianˇcn´ı (varianˇcn´ı) matice Kovarianˇcn´ı matice V typu (p×p) je symetrická (vlastn´ı ˇc´ısla jsou reálná) a pozitivnˇe semidefinitn´ı, tj. yT Vy ≥ 0 pro jak´ ykoliv vektor y ̸= 0. Tzn., ˇze vˇsechna vlastn´ı ˇc´ısla jsou nezáporná, – d˚ ukaz viz napˇr. Andˇel, str. 28 [2] Vlastn´ı ˇc´ısla m˚ uˇzeme uspoˇrádat: λ1 ≥ λ2 ≥ . . . ≥ λp ≥ 0 Matici V m˚ uˇzeme napsat jako V = UΛUT ,

UUT = I,

(58)

kde Λ = diag(λ1 , . . . , λp ) a k – t´ y sloupec matice U je vlastn´ı vektor vk matice V, kter´ y pˇr´ısluˇs´ı vlastn´ımu ˇc´ıslu λk a pro kter´ y plat´ı vkT vk = 1. Tedy vk jsou ortonormáln´ı vlastn´ı vektory kovarianˇcn´ı matice V, plat´ı

V

Vvk = λk vk T = λ1 v1 v1 + λ2 v2 v2T + . . . + λp vp vpT I = v1 v1T + v2 v2T + . . . + vp vpT

Vektory v1 , v2 , . . . , vp tvoˇr´ı bázi prostoru Rp , takˇze libovoln´ y vektor y ∈ Rp lze vyjádˇrit jako y = c1 v1 + c2 v2 + . . . + cp vp Plat´ı pro kaˇzd´ y vektor y ∈ Rp jednotkové délky (yT y = 1), ˇze kvadratická forma yT Vy ≤ λ1 , pˇri ˇcemˇz rovnost plat´ı pro c = v1 , tj. pro prvn´ı vlastn´ı vektor – viz Andˇel, str. 297 [2]. Celková variabilita náhodného vektoru x = (X1 , X2 , . . . , Xp )T je souˇcet diagonáln´ıch prvk˚ u kovarianˇcn´ı matice (rozptyl˚ u jednotliv´ ych náhodn´ ych veliˇcin): σ 2 = var(X1 ) + var(X2 ) + . . . + var(Xp ) Hledáme novou náhodnou veliˇcinu, která vznikne lineárn´ı transformac´ı vektoru x = (X1 , X2 , . . . , Xp )T , tj. vlastnˇe hledáme c ∈ Rp , cT c = 1, aby náhodná veliˇcina mˇela

12.4 Anal´ yza hlavn´ıch komponent

117

co nejvˇetˇs´ı rozptyl. T´ım tato nová veliˇcina vyˇcerpá co nejvˇetˇs´ı ˇcást celkové variability, tzn. jej´ı rozptyl je roven λ1 . Tedy c = v1 . Náhodnou veliˇcinu Y1 = v1T x naz´ yváme prvn´ı hlavn´ı komponentou. Pak hledáme dalˇs´ı náhodnou veliˇcinu, tj. dalˇs´ı c ∈ Rp , cT c = 1, aby náhodná veliˇcina cT x byla nekorelovaná s Y1 = v1T x. Tomu vyhovuje c = v2 , takˇze druhá hlavn´ı komponenta je Y2 = v2T x a jej´ı rozptyl je var(Y2 ) = var(v2T x) = λ2 . Podobnˇe i pro tˇret´ı a dalˇs´ı hlavn´ı komponenty. Celková variabilita σ 2 = var(X1 ) + var(X2 ) + . . . + var(Xp ) = = var(Y1 ) + var(Y2 ) + . . . + var(Yp ) =

p ∑

λi

i=1

Relativn´ı pod´ıl celkové variability vysvˇetlovan´ y i-tou hlavn´ı komponentou je λi /σ 2 . ∑k Prvn´ıch k hlavn´ıch komponent vysvˇetluje i=1 λi /σ 2 z celkové variability. V praktick´ ych aplikac´ıch se vycház´ı z v´ ybˇerové kovarianˇcn´ı matice nebo ˇcastˇeji z v´ ybˇerové korelaˇcn´ı matice (abychom se vyhnuli vlivu volby jednotek, ve kter´ ych mˇeˇr´ıme veliˇciny, na hodnoty v´ ybˇerov´ ych rozptyl˚ u). Korelaˇcn´ı matice je vlastnˇe kovarianˇcn´ı matice standardizovan´ ych veliˇcin. Pak celková variabilita je rovna poˇctu veliˇcin, σ2 =

p ∑ i=1

λi = p


118

Pˇ r´ıklad 12.7 Základn´ı moˇznosti anal´ yzy hlavn´ıch komponent ukáˇzeme na pˇr´ıkladu o struktuˇre zamˇestnanosti ve 26 evropsk´ ych zem´ıch, data jsou v souboru EMPLOY.XLS. Vycház´ıme z v´ ybˇerové korelaˇcn´ı matice. Principal Components Report Database

C:\avdat\employ.S0

Eigenvalues No. 1 2 3 4 5 6 7 8 9

Eigenvalue 3.487151 2.130173 1.098958 0.994483 0.543218 0.383428 0.225754 0.136790 0.000046

Individual Percent 38.75 23.67 12.21 11.05 6.04 4.26 2.51 1.52 0.00

Cumulative Percent Scree Plot 38.75 |||||||| 62.41 ||||| 74.63 ||| 85.68 ||| 91.71 || 95.97 | 98.48 | 100.00 | 100.00 |

Vid´ıme, ˇze tˇri vlastn´ı ˇc´ısla jsou vˇetˇs´ı neˇz 1, ˇctvrté témˇeˇr rovno jedné. Prvn´ı dvˇe hlavn´ı komponenty vysvˇetluj´ı 62 % z celkové variability, prvn´ı ˇctyˇri uˇz 85 % celkové variability. Prvn´ı ˇctyˇri vlastn´ı vektory jsou v následuj´ıc´ı tabulce: Eigenvectors Variables Factor1 AGR 0.523791 MIN 0.001323 MAN -0.347495 PS -0.255716 CON -0.325179 SER -0.378920 FIN -0.074374 SPS -0.387409 TC -0.366823

Factor2 0.053594 0.617807 0.355054 0.261096 0.051288 -0.350172 -0.453698 -0.221521 0.202592

Factor3 0.048674 -0.201100 -0.150463 -0.561083 0.153321 -0.115096 -0.587361 0.311904 0.375106

Factor4 0.028793 0.064085 -0.346088 0.393309 -0.668324 -0.050157 -0.051567 0.412230 0.314372

Na dvourozmˇern´ ych grafech v rovinách dvojic prvn´ıch tˇr´ı hlavn´ıch komponent vid´ıme, ˇze pˇri takto podstatném sn´ıˇzen´ı dimenze, ve kter´ ych objekty zobrazujeme, lze sledovat podobnosti a odliˇsnosti zobrazen´ ych objekt˚ u. Zejména na obrázku prvn´ıch

12.4 Anal´ yza hlavn´ıch komponent

119

dvou hlavn´ıch komponent jsou viditelné shluky objekt˚ u (zem´ı) v souladu s klasifikac´ı nalezenou shlukovou anal´ yzou.

120


12.5 Faktorová anal´ yza

12.5

121

Faktorov´ a anal´ yza

Faktorová anal´ yza je jedna z metod redukce dimenze u ´lohy. Snaˇz´ı se vysvˇetlit kovarianˇcn´ı strukturu (korelaˇcn´ı matici) vektoru náhodn´ ych veliˇcin, pomoc´ı ménˇe tzv. faktor˚ u, tj. jak´ ychsi skryt´ ych veliˇcin, které nem˚ uˇzeme nebo neum´ıme pˇr´ımo mˇeˇrit. Uvaˇzujme, ˇze namˇeˇrená data jsou matice X typu n × p (n objekt˚ u, p veliˇcin). Standardizovaná matice dat je matice Z opˇet typu n × p, kde zij =

xij − x¯j sj

x¯j , sj jsou v´ ybˇerov´ y pr˚ umˇer a smˇerodatná odchylka j-té veliˇciny. Model faktorové anal´ yzy m˚ uˇzeme pak zapsat zij = aj1 fi1 + aj2 fi2 + . . . + ajm fim + eij ,

m
tj. namˇeˇrenou hodnotu j-té veliˇciny na i-tém objektu vysvˇetlujeme jako váˇzen´ y souˇcet m faktor˚ u a nˇejaké sloˇzky, která tˇemito faktory vysvˇetlit nelze. Zavedeme následuj´ıc´ı matice: A je typu (p × m), má prvky ajk , oznaˇcuj´ı se jako faktorové zátˇeˇze (sycen´ı, saturace, loadings) F je typu (n × m), má prvky fik , ˇr´ıká se jim faktorové skóry E je typu (n × p), má prvky eij , jsou to rezidua, tj. to z hodnot zij , co nem˚ uˇzeme vysvˇetlit pomoc´ı faktor˚ u Maticovˇe pak m˚ uˇzeme model faktorové anal´ yzy zapsat takto: Z = FAT + E Pˇredpokládejme, ˇze faktory jsou ortonormáln´ı (nekorelované, jednotkové délky), tzn. FT F = I. Oznaˇcme U = ET E. U je diagonáln´ı matice typu p × p, tedy ∑ U = diag(u1 , u2 , . . . , up ), kde uj = ni=1 e2ij , tj. variabilita j-té veliˇciny, kterou nelze vysvˇetlit faktory, tzv. specificita. Pak hj = 1 − uj =

m ∑

a2jk

k=1

je tzv. komunalita j-té veliˇciny, tj. variabilita vysvˇetlitelná faktory. Korelaˇcn´ı matici m˚ uˇzeme vyjádˇrit jako R = AAT + U


122

Matice AAT vysvˇetluje korelaˇcn´ı matici aˇz na diagonáln´ı prvky, kde m´ısto jedniˇcek jsou komunality. Faktorová anal´ yza hledá model, aby – pˇr´ıspˇevek specifick´ ych faktor˚ u (uj ) byl co nejmenˇs´ı – faktorové zátˇeˇze byly v absolutn´ı hodnotˇe co nejbliˇzˇs´ı jedné nebo nule – poˇcet faktor˚ u byl co nejmenˇs´ı (podstatnˇe menˇs´ı neˇz poˇcet veliˇcin) Tyto poˇzadavky jsou ve vzájemném rozporu a umˇen´ı“ faktorové anal´ yzy spoˇc´ıvá ” v nalezen´ı vhodného a pˇrijatelného kompromisu: extrakce faktor˚ u - stanovit jejich poˇcet, napˇr. z vlastn´ıch ˇc´ısel korelaˇcn´ı matice probl´ em komunalit Rj2 ≤ hj ≤ 1, kde Rj2 je koeficient mnohonásobné korelace (j-tá veliˇcina na ostatn´ıch) rotace faktor˚ u tj. nalezen´ı jednoduché struktury“, aby faktorové zátˇeˇze v abso” lutn´ı hodnotˇe byly co nejbliˇzˇs´ı jedné nebo nule Ortogonáln´ı rotace znamená naj´ıt takovou transformaci matice A na B = AT, aby B splˇ novalo poˇzadavek jednoduché struktury, T je ortogonáln´ı matice, tj. TT T = I, takˇze AAT = BBT . Je mnoho moˇzn´ ych metod takové rotace faktor˚ u, nejbˇeˇznˇejˇs´ı je tzv. VARIMAX, zaloˇzená na maximalizaci v´ yrazu 1 ∑∑ 2 (a − a2·j )2 , p j=1 i=1 ij m

kde a2·j =

1 p

∑p i=1

p

a2ij , tj. pr˚ umˇer ˇctverc˚ u zátˇeˇz´ı pro j-t´ y faktor.

Po rotaci m˚ uˇzeme nahlédnout, která veliˇcina patˇr´ı“ kterému faktoru (faktorové ” zátˇeˇze v absolutn´ı hodnotˇe bl´ızké jedné), pˇr´ıpadnˇe faktorové zátˇeˇze jednotliv´ ych veliˇcin vynést do graf˚ u pro dvojice faktor˚ u. Lze pak spoˇc´ıtat i faktorové skóry a na grafech faktorov´ ych skór˚ u hledat, zda zobrazené objekty nevytváˇrej´ı nˇejaké shluky signalizuj´ıc´ı moˇzn´ y rozklad pozorovan´ ych objekt˚ u do dvou ˇci v´ıce podskupin. Pˇ r´ıklad 12.8 Data pro tento pˇr´ıklad jsou pˇrevzata z knihy [6]. P˚ uvodn´ı data byla v´ ysledky dosaˇzené ve v´ ybˇeru 220 chlapc˚ u v ˇsesti pˇredmˇetech - galˇstinˇe, angliˇctinˇe, dˇejepisu, aritmetice, algebˇre a geometrii. K dispozici pro anal´ yzu máme v´ ybˇerovou korelaˇcn´ı matici (doln´ı troj´ uheln´ık, pˇredmˇety jsou ve v´ yˇse uvedeném poˇrad´ı): 1.00 0.44

1.00


0.41 0.29 0.33 0.25

0.35 0.35 0.32 0.33

1.00 0.16 0.19 0.18

123

1.00 0.59 0.47

1.00 0.46

1.00

Abychom si uvˇedomili rozd´ıly mezi anal´ yzou hlavn´ıch komponent a faktorovou anal´ yzou, uvád´ıme nejdˇr´ıve struˇcné v´ ysledky anal´ yzy hlavn´ıch komponent, která vysvˇetluje celkov´ y rozptyl, tedy hlavn´ı diagonálu korelaˇcn´ı matice. Principal Components Report Database

C:\avdat\subject.S0

Eigenvalues No. 1 2 3 4 5 6

Eigenvalue 2.728683 1.128792 0.615291 0.602809 0.522514 0.401910

Individual Percent 45.48 18.81 10.25 10.05 8.71 6.70

Factor Loadings Variables Factor1 Gaelic -0.660803 English -0.688465 History -0.516356 Arithmetic -0.735620 Algebra -0.741868 Geometry -0.678168

Cumulative Percent 45.48 64.29 74.55 84.59 93.30 100.00

Scree Plot |||||||||| |||| ||| ||| || ||

Factor2 -0.444475 -0.289771 -0.639552 0.417018 0.372759 0.354100

Dvˇe vlastn´ı ˇc´ısla jsou vˇetˇs´ı neˇz jedna, pro faktorovou anal´ yzu tedy zvol´ıme poˇcet faktor˚ u 2, rotaci varimax a dostaneme následuj´ıc´ı v´ ysledky:


124

Factor Analysis Report Database

C:\avdat\subject.S0

Eigenvalues after Varimax Rotation Individual Cumulative No. Eigenvalue Percent Percent 1 1.596863 56.94 56.94 2 1.207981 43.08 100.02 3 0.050820 1.81 101.83 4 0.011910 0.42 102.26 5 -0.008657 -0.31 101.95 6 -0.054642 -1.95 100.00

Scree Plot |||||||||||| ||||||||| | | | |

Vlastn´ı ˇc´ısla se t´ ykaj´ı matice AAT = R − U po rotaci, c´ılem faktorové anal´ yzy je pˇredevˇs´ım vysvˇetlit korelaˇcn´ı strukturu, tj. mimodiagonáln´ı prvky korelaˇcn´ı matice. Celkov´ y vysvˇetlen´ y rozptyl je roven jen souˇctu komunalit. Po rotaci jsou faktorové zátˇeˇze následuj´ıc´ı: Factor Loadings after Varimax Rotation Variables Factor1 Factor2 Gaelic -0.233132 -0.659253 English -0.322810 -0.552071 History -0.084713 -0.589192 Arithmetic -0.765986 -0.170657 Algebra -0.718105 -0.214689 Geometry -0.573340 -0.214994

Absolutn´ı hodnoty faktorov´ ych zátˇeˇz´ı jsou znázornˇeny graficky: Bar Chart of Absolute Factor Loadings after Varimax Rotation Variables Gaelic English History Arithmetic Algebra

Factor1 ||||| ||||||| || |||||||||||||||| |||||||||||||||

Factor2 |||||||||||||| |||||||||||| |||||||||||| |||| |||||


Geometry

125

||||||||||||

|||||

Faktor 1 je sycen veliˇcinami aritmetika, algebra a geometrie, faktor 2 veliˇcinami galˇstina, angliˇctina a dˇejepis. Pod´ıly faktor˚ u na komunalitách ukazuje následuj´ıc´ı tabulka: Communalities Factors Variables Gaelic English History Arithmetic Algebra Geometry

after Varimax Rotation Factor1 0.054350 0.104206 0.007176 0.586735 0.515675 0.328719

Factor2 0.434614 0.304783 0.347147 0.029124 0.046091 0.046222

Communality 0.488965 0.408989 0.354324 0.615859 0.561766 0.374942

Vztah veliˇcin k faktor˚ um je graficky znázornˇen v grafu faktorov´ ych zátˇeˇz´ı, ve kterém vid´ıme shluk veliˇcin 1, 2 a 3 s n´ızk´ ymi absolutn´ımi hodnotami zátˇeˇz´ı prvn´ıho faktoru, které syt´ı pˇredevˇs´ım druh´ y faktor, a shluk veliˇcin 4, 5 a 6 s n´ızk´ ymi absolutn´ımi hodnotami zátˇeˇz´ı druhého faktoru, syt´ıc´ıch pˇredevˇs´ım prvn´ı faktor.

126


Korelaˇcn´ı strukturu pozorovan´ ych dat (studijn´ıch v´ ysledk˚ u v ˇsesti pˇredmˇetech) lze tedy uspokojivˇe vysvˇetlit dvˇe faktory. Prvn´ı faktor vyjadˇruje matematickou dispozici ˇza´ka, druh´ y dispozici jazykovˇe-humanitn´ı.

Shrnut´ı

• u ´lohy ˇreˇsené mnohorozmˇernými metodami • test shody vektor˚ u stˇredn´ıch hodnot • metody klasifikace objekt˚ u, diskriminaˇcn´ı analýza, logistick´ a regrese • metody shlukován´ı, numerická podobnost, hierarchické a nehierarchické metody shlukován´ı • redukce dimenze, analýza hlavn´ıch komponent, faktorov´ a analýza Kontroln´ı ot´ azky

1. V ˇcem jsou shodné a v ˇcem odliˇsné c´ıle diskriminaˇcn´ı analýzy a shlukové analýzy? 2. Proˇc se pˇri hledán´ı lineárn´ı diskriminaˇcn´ı funkce mus´ı liˇsit stˇredn´ı hodnoty skupin? 3. Jaké jsou výhody a nevýhody lineárn´ı diskriminaˇcn´ı funkce oproti jiným klasifikaˇcn´ım pravidl˚ um (napˇr. logistick´ a regrese, neuronové s´ıtˇe ap.)? 4. Jaké jsou rozd´ıly mezi hierarchickými a nehierarchickými metodami? 5. V ˇcem se liˇs´ı faktorová analýza od analýzy hlavn´ıch komponent? 6. Jak zjist´ıte souˇradnice jednotlivých objekt˚ u (ˇr´ adk˚ u datové matice) v rovinˇe prvn´ıch dvou hlavn´ıch komponent? 7. Jak urˇcit poˇcet faktor˚ u? 8. Co je to komunalita? 9. Co jsou faktorové zátˇeˇze? Co lze vyˇc´ıst z grafu faktorových zátˇeˇz´ı? 10. Jakou ˇcást celkové variability vysvˇetluje prvn´ı hlavn´ı komponenta? Korespondenˇ cn´ı u ´ loha Korespondenˇcn´ı u ´lohy budou zadáv´ any ke kaˇzdému kursu samostatnˇe.

127

13

Literatura

[1] Afifi A., Clark V. A., May S., Computer-Aided Multivariate Analysis, Chapman & Hall/CRC, 2004 [2] Andˇel J.: Matematická statistika, SNTL, Praha, 1978 [3] Andˇel J.: Statistické metody, Matfyzpress, Praha, 1993 [4] Antoch, J., Vorl´ıˇcková, D.: Vybrané metody statistické anal´ yzy dat. Academia, Praha, 1992. [5] Armitage P., Berry G.: Statistical Methods in Medical Research, Blackwell Sci Publ., 1994 [6] Bartholomew D. J., Steel F., Moustaki I., Galbraith J.I., The Analysis and Interpretation of Multivariate Data for Social Sciences, Chapman and Hall/CRC, 2002 [7] Cipra., T., Ekonometrie, MFF UK, SPN Praha, 1984 [8] Draper N. R., Smith H.: Applied Regression Analysis, John Wiley, 1998 [9] Hair J. F.jr, Black W. C., Babin B. J., Anderson R. E., Multivariate Data Analysis, Prentice Hall, 2010 [10] Havránek T.: Statistika pro biologické a lékaˇrské vˇedy, Academia, Praha, 1993 [11] Havránek T. a kol.: Matematika pro biologické a lékaˇrské vˇedy, Academia, Praha, 1981 [12] Hebák P. a kol. V´ıcerozmˇerné statistické metody (1), Informatorium, Praha, 2004 [13] Hebák P. a kol. V´ıcerozmˇerné statistické metody (2),(3), Informatorium, Praha, 2005 [14] Hintze, J., NCSS 8. NCSS, LLC. Kaysville, Utah, USA. www.ncss.com, 2012 [15] Kleinbaum D. G.: Logistic Regression: A Self-Learning Text, Springer-Verlag, New York Berlin Heidelberg, 1994 [16] Likeˇs J., Laga J., Základn´ı statistické tabulky, SNTL, Praha, 1978 ˇ [17] Lukasová A., Sarmanov´ a J., Metody shlukové anal´ yzy, SNTL, Praha, 1981

128

[18] Manly, B. F. J., Multivariate Statistical Methods - A primer, Chapman and Hall/CRC, 1994 [19] McCullagh, P., Nelder, J. A.: Generalized Linear Model, 2nd Edition, Chapman and Hall, 1989 [20] Meloun M., Militk´ y J.: Statistické zpracován´ı experimentáln´ıch dat, PLUS, 1994 [21] Tvrd´ık, J., Základy pravdˇepodobnosti a statistiky, OU, 2010 [22] Tvrd´ık, J., Anal´ yza dat, OU, 2008 [23] Zvára, K., Regrese, Academia, Praha, 1989

28

Recommend Documents