Poznámky z kurzu statistiky. Statistikem? Těžko, ale rychle

Poznámky z kurzu statistiky Jak to udělat v Rku

Statistikem? Těžko, ale rychle. Jáchym Čepický 16. března 2006

Tento text vznikl z poznámek ze statistického kurzu, vedeného panem Ing. Karlem Drápelou, CSc., který proběl v lednu roku 2005. Řešení jednotlivých úloh statistické analýzy jsem se pokusil ukázat na statistickém programu R. I když v textu přímo necituji, uvádím na koncii seznam použité literatury.

Tento text je uvolněn pod licencí GNU/FDL: Copyright (c) 2005 Jáchym Čepický Je povoleno kopírovat, šířit a/nebo upravovat tento dokument za podmínek GNU Free Documentation License, verze 1.2 nebo jakékoli další verze vydané nadací Free Software Foundation; bez neměnných oddílů, bez textů předních desek a bez textů zadních desek. Kopie této licence je zahrnuta v oddílu jménem ”GNU Free Documentation License”.

k

c

Jáchym Čepický 2005 [email protected] domovská stránka Les-ej (http://les-ejk.cz)

Obsah 1 Základní pojmy

1

2 Postup statistické analýzy

3

3 Metody průzkumové analýzy dat 3.1 Metody grafické . . . . . . . . . . . . . 3.1.1 Krabicový „box“ graf . . . . . 3.1.2 Kvantil-kvantilový graf „Q-Q“ 3.1.3 Histogram . . . . . . . . . . . . 3.1.4 Autokorelace . . . . . . . . . . 3.2 Metody početní - testy . . . . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

4 Odhad parametrů základního souboru 4.1 Bodový odhad . . . . . . . . . . . . . . . . . . . . . 4.1.1 Odhad střední hodnoty . . . . . . . . . . . . 4.2 Intervalový odhad . . . . . . . . . . . . . . . . . . . 4.2.1 Střední hodnota . . . . . . . . . . . . . . . . 4.2.2 Faktory ovlivňující interval spolehlivosti . . . 4.2.3 Interval spolehlivosti směrodatné odchylky σ

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . .

5 5 6 6 7 9 10

. . . . . .

13 13 13 13 13 14 14

5 Transformace

17

6 Práce s malými výběry

19

7 Výběry 7.1 Metody výběru . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Zjišťování velikosti výběrů . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21 21 21

8 Testy 8.1 Testové kritérium . . 8.2 Postup testu . . . . 8.3 P-hodnota . . . . . . 8.4 Neparametrické testy

. . . .

25 25 25 26 31

9 Síla testu 9.1 Analýza síly testu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9.1.1 Postup analýzy síly testu . . . . . . . . . . . . . . . . . . . . . . . . .

33 34 34

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

i

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

9.2 9.3

Pro statistické testy o střední hodnotě . . . . . . . . . . . . . . . . . . . . . . Testy – přehled . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10 Korelační a regresní analýza 10.1 Korelace . . . . . . . . . . . . . . . . . . . . . . . . 10.1.1 Míra korelační závislosti . . . . . . . . . . . 10.1.2 Korelační koeficient . . . . . . . . . . . . . 10.2 Regrese . . . . . . . . . . . . . . . . . . . . . . . . 10.2.1 Příklady regresních modelů . . . . . . . . . 10.2.2 Postup regresní analýzy . . . . . . . . . . . 10.2.3 Předpoklady metody . . . . . . . . . . . . . 10.2.4 Multikolinearita . . . . . . . . . . . . . . . 10.2.5 Řešení multikolinearity . . . . . . . . . . . . 10.2.6 Homoskedasticita a Heteroskedasticita . . . 10.2.7 Postup regresní analýzy (lineární model) . . 10.2.8 Postup regresní analýzy (nelineární model) 11 Anova 11.0.9 11.1 Model 11.2 Model 11.3 Model

Podmínky použití . . . . . . . . . jednofaktorové anovy . . . . . . . . vícefaktorové anovy bez opakování vícefaktorové anovy s opakováním .

ii

. . . .

. . . .

. . . .

. . . .

. . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

39 39

. . . . . . . . . . . .

41 41 41 42 45 46 46 47 47 48 48 48 59

. . . .

65 66 67 68 69

Seznam obrázků 3.1 3.2 3.3 3.4

Čtyři různé „boxploty“ prvních čtyř sloupců Různé průběhy Q-Q grafu. . . . . . . . . . . Různé histogramy . . . . . . . . . . . . . . Graf autokorelace . . . . . . . . . . . . . . .

ze souboru x. . . . . . . . . . . . . . . . . . . . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

7 8 9 10

8.1

Znázornění kritického bodu a chyby prvního řádu . . . . . . . . . . . . . . . .

26

10.1 10.2 10.3 10.4 10.5 10.6 10.7 10.8

Princip regresní analýzy . . . . . . . . . . . . . . Regresní model . . . . . . . . . . . . . . . . . . . Homo- a heteroskedasticita . . . . . . . . . . . . Graf regresní přímky a dat . . . . . . . . . . . . Několik důležitých diagnostických grafů . . . . . Různé tvary rozložení reziduí . . . . . . . . . . . Interval spolehlivosti pro lineární model a interval Graf Michajlovovy růstové funkce . . . . . . . . .

42 46 49 50 51 56 59 62

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . spolehlivosti predikce modelu. . . . . . . . . . . . . . . . .

11.1 Znázornění rozdílů mezi zásahy pomocí boxplotu . . . . . . . . . . . . . . . .

iii

67

iv

1 Základní pojmy Základní soubor • Skládá se ze všech prvků

• Je to to, co zkoumám, soubor, jehož charakteristiky mě zajímají. • Jeho popisem se zabývá popisná statistika

• Musí být jasně definován

Základní soubor může mít nekonečně mnoho prvků. Výběrový soubor • Jedná se podmnožinu základního souboru.

• Prvky jsou vybrány náhodným výběrem

• Výběrový soubor je pouze prostředek k poznání základního souboru (V S → ZS).

• Používá se statistická indukce.

• Jako nástroj se používá matematická statistika – pravděpodobností počet. Výběrové statistiky • Charakteristické vlastnosti výběrového souboru.

• Parametry základního souboru se odhadují na základě charakteristik souboru výběrového.

1. ZS → VZ → Výb. statistiky → OdhadparametrůZS Momentové charakteristiky (σ, σ 2 , µ, koef. špičatosti, koef. nesouměrnosti, . . .) Aby mělo smysl počítat momentové charakteristiky, musí výběrový soubor splňovat několik podmínek: 1. Musí se jednat o normální rozdělení 2. Nesmí obsahovat extrémní hodnoty 3. Data musí být navzájem nezávislá 4. Musí jich být nějaký minimální počet Při splnění těchto podmínek můžeme získat dobré popisné vlastnosti souboru. V případě, že tyto podmínky splněny nejsou, musíme použít kvantilové charakteristiky. 1

Kvantilové charakteristiky jsou závislé na pořadí prvků v souboru. Jedná se zejména o x ˆ, kvartily, kvantilové odchylky. • Výhodou je nezávislost na na podmínkách, které musí splňovat momentové charakteristiky, takže soubor dat může vypadat v podstatě jakkoliv. • Nevýhodou je, že nevychází ze všech prvků v souboru, ale pouze z jejich pořadí.

2

2 Postup statistické analýzy Každá statistická analýza se musí sestávat z následujících kroků: 1. Stanovení cíle analýzy 2. Získání dat 3. Jejich zpracování 4. Statistické vyhodnocení výsledků 5. Interpretace výsledků z hlediska statistické analýzy s ohledem na cíl analýzy. Abychom mohli rozhodnout, jaké charakteristiky použijeme, musíme provést následující kroky: 1. Průzkumovou analýzu dat, ze které se získají hodnoty jako stupeň špičatosti, stupeň symetrie, odlehlost dat, porovnání s rozděleními (hlavně s normálním). 2. Ověření předpokladů o datech, hlavně ověření jejich normality, nezávislosti prvků výběru, homogenity, potřebné velikosti výběru. 3. Odhad parametrů ZS, zejména výpočet charakteristiky, bodových odhadů parametrů a intervalových odhadů parametrů. 4. Testování hypotéz, což znamená především jejich formulaci, potvrzení nebo zamítnutí, a ověření síly testu.

3

4

3 Metody průzkumové analýzy dat 3.1

Metody grafické

• ⊕ výsledky jsou „dobře vidět“ • ⊖ ale neřeknou nic jednoznačně.

Načtení dat do programu R Mějme např. soubor data.txt, který obsahuje ve sloupcích dvoje různé hodnoty: tloustky 42.9 36.1 46.7 53.9 53.4 57.1 29.7 43.4 52.7 37.4 40.2 33.2 ...

studny 32 75 16 99 80 78 28 29 170 86 80 81

letokruhy 2.1 2.8 2.5 3.3 3.3 1.9 3.1 2.5 2.6 3.2 3.8 2.9

med 1.098 1.876 2.149 0.904 1.128 0.468 0.428 0.175 0.05 0.074 0.05 0.64

zinek 1.231 0.654 5.8 0.963 1.795 3.808 0.982 0.191 0.161 0.56 0.295 0.435

kadmium 0.08 0.018 0.015 0.052 0.08 0.02 0.036 0.013 0.032 0.06 0.07 0.019

ovzdusi 4.61 3.8 10.53 34.4 31.48 109.58 16.34 0.5 13.07 18.31 22.39 23.36

BSK 6.5 5.8 16.7 6.4 7 6.3 7 9.2 6.7 6.7

Tabulka 3.1: Soubor data.txt, obsahuje sloupce hodnot. Sloupce jsou od sebe odděleny znakem . Čísla jsou oddělena desetinnou tečkou (ne čárkou).

$ # spuštění programu R se děje příkazem $ R [...] R> # načtení hodnot do proměnné ’x’ R> # soubor se jmenuje ’data1.txt’, má hlavičku a jednotlivé hodnoty jsou R> # odděleny znakem tabulátoru R> x <- as.matrix(read.table("data1.txt",header=T,sep="\t")) R> R> # na data se můžeme podívat příkazem 5

R> x [...] R> R> # všechny prvky z prvního sloupce si vypíšeme příkazem R> x[,1] [...] R> R> # třetí prvek druhého sloupce si můžeme nechat zobrazit příkazem R> x[3,2] [1] 16 R> R> # prvky 1 až 10 z pátého sloupce získáme R> x[1:10,5] 1 2 3 4 5 6 7 8 9 10 1.231 0.654 5.800 0.963 1.795 3.808 0.982 0.191 0.161 0.560 R> R> # celkovou charakteristiku dat získáme příkazem R> summary(x) [...] R> summary(x[,1]) Min. 1st Qu. Median Mean 3rd Qu. Max. 27.00 39.25 44.40 44.72 50.08 61.60

3.1.1

Krabicový „box“ graf

Krabicový graf ukazuje základní charakteristiku sbíraného souboru dat. Krabicový graf ukazuje základní charakteristiku dat. Krabice („box“) obsahuje horní a dolní kvantil dat (tedy 25% a 75% všech dat). Střední osa je medián, „nožičky“ ukazují tzv. rozsah nevybočujících hodnot. „Kolečka“ jsou hodnoty odlehlé. Extrémní hodnoty bývají kandidáty na vyřazení. R> R> R> R> R> R> R> R> R> R> R> R>

# vytištění boxplotu z dat prvního sloupce: boxplot(x[,1],horizontal=T,main="Tloustky x[,1]") boxplot(x[,2],horizontal=T,main="Tloustky x[,1]") #... # všechny grafy do jednoho okna par(mfrow=c(2,2)) # nastavení okna na 2x2 buněk boxplot(x[,1],horizontal=T,main="Tloustky x[,1]") boxplot(x[,2],horizontal=T,main="Studny x[,2]") boxplot(x[,3],horizontal=T,main="Letokruhy x[,3]") boxplot(x[,4],horizontal=T,main="Med x[,4]") par(mfrow=c(1,1)) # nastavení okna zpět na jedno velké

3.1.2

Kvantil-kvantilový graf „Q-Q“

Ukazuje na ose seřazení hodnot podle velikosti jejich výskytu (vzestupně uspořádána). Tyto hodnoty bývají zobrazovány bodově. V grafu bývá linie teoretického ideálního rozdělení. 6

Tloustky x[,1]

30

40

50

Studny x[,2]

60

0

100

Letokruhy x[,3]

300

500

Med x[,4]

2.0 2.5 3.0 3.5 4.0 4.5

0

2

4

6

8

Obrázek 3.1: Čtyři různé „boxploty“ prvních čtyř sloupců ze souboru x.

R> R> R> R> R>

# vytištění qqplotu qqplot(x[,1],main="Normální Q-Q plot Tloustek") # zobrazení linie teoretického rozdělení qqline(x[,1])

Na Q-Q grafu je dobře vidět, o jaká data se vlastně jedná. Jsou-li data „normální“, následují skoro přesně teoretickou přímku. Druhý graf ukazuje průběh, jaký mají data pravostranná (viz část 3.1.3), třetí graf, s průběhem okolo teoretické linie ve tvaru písmena „S“ ukazuje data špičatá, pokud by se jednalo o tvar obráceného „S“, jednalo by se o data plochá.

3.1.3

Histogram

Histogram ukazuje rozdělení dat do jednotlivých tříd podle četnosti. Ideální počet tříd lze vypočítat podle vzorců L = 2, 46 · (n − 1)0,4 √ L=2· n

(3.1) (3.2)

V obou případech L – počet tříd, n – počet prvků. Je na čase, seznámit se s příkazem ’help’, se kterým se většinou člověk dobere toho, co vlastně potřebuje Jelikož chceme vytisknout ’histogram’, necháme si toto slovo najít 7

0

200

400

600

Normální Q−Q plot Studny

Sample Quantiles

60 50 40 30

Sample Quantiles

Normální Q−Q plot Tloustek

−1

0

1

2

−2

−1

0

1

2

Normální Q−Q plot Letokruhy

Normální Q−Q plot Med

6 4 0

2

4.0 3.0

8

Theoretical Quantiles

Sample Quantiles

Theoretical Quantiles

2.0

Sample Quantiles

−2

−2

−1

0

1

2

−2

Theoretical Quantiles

−1

0

1

2

Theoretical Quantiles

Obrázek 3.2: Různé průběhy Q-Q grafu.

R> help.search(’histogram’) [...] hist.scott(MASS) Plot a Histogram with Automatic Bin Width Selection ldahist(MASS) Histograms or Density Plots of Multiple Groups truehist(MASS) Plot a Histogram hist.POSIXt(graphics) Histogram of a Date or Date-Time Object hist(graphics) Histograms nclass.Sturges(graphics) Compute the Number of Classes for a Histogram plot.histogram(graphics) Plot Histograms [...] Ukazuje se, že pro vytisknutí histogramu existuje celá řada funkcí. Každá funkce má v závorce napsáno, v jakém balíku se vyskytuje. Použijeme základní funkci ’hist’ z balíku ’graphics’, která vytiskne pouze základní histogram. Funkce ’truehist’ z balíku ’MASS’ umí histogramy o poznání lepší. O instalaci balíku MASS viz Google, alespoň prozatím. . . Když už víme, jakou funkci chceme, zobrazíme si nápovědu, která nám napíše jak se daná funkce vlastně používá. Možností máme několik. Jednak praktický příklad: R> example(hist) 8

Histogram Studen

20

40

Cetnost

10

0

0

5

Cetnost

20

60

80

Histogram Tloustek

30

40

50

60

0

100

300

500

Studny

Histogram Letokruhu

Histogram Medu

15 5

10

Cetnost

4 3 2

0

0

1

Cetnost

5

20

6

Tloustky

1.5

2.5

3.5

4.5

0

Letokruhy

2

4

6

8

10

Med

Obrázek 3.3: Různé histogramy

A jednak vlastní popis funkce s vysvětlením parametrů: R> help(hist) R> # alternativně a rychleji R> ?hist R> Zobrazí krátký popis funkce (Description), použití (Usage), podrobný popis argumentů pro použití (Arguments) a další a další text. My si z toho vezmeme následující R> R> R> R>

hist(x[,1],main="Histogram hist(x[,2],main="Histogram hist(x[,3],main="Histogram hist(x[,4],main="Histogram

Tloustek",xlab="Tloustky",ylab="Cetnost") Studen",xlab="Tloustky",ylab="Cetnost") Letokruhu",xlab="Tloustky",ylab="Cetnost") Medu",xlab="Tloustky",ylab="Cetnost")

Na obrázcích histogramech (obrázek 3.3) je dobře patrno to, co jsme už mohli usoudit na základě Q-Q plotu. Jsou vidět rozdělení silně levostranná a plochá.

3.1.4

Autokorelace

Při autokorelaci data „sami se sebou“ korelují. Hovoříme o korelaci řádu n. Jedná-li se např. o autokorelaci prvního řádu, pak je prvek x + 1 závislý na prvku x. Autokorelaci můžeme 9

Series x[1:99, 2]

0.6 −0.2

5

10

15

20

0

5

10

15

Lag

Lag

Series x[1:38, 3]

Series x[1:34, 4]

−0.2

0.2 −0.2

0.2

ACF

0.6

0.6

1.0

1.0

0

ACF

0.2

ACF

0.2 −0.2

ACF

0.6

1.0

1.0

Series x[1:100, 1]

0

5

10

15

0

Lag

5

10

15

Lag

Obrázek 3.4: Graf autokorelace pro první čtyři sloupce dat.

odhalit např. pomocí grafu (obrázek 3.4). Na ose x je časová řád autokorelace a na ose y její stupeň. Přesáhne-li hodnota stupně hranice, vyjádřené modrou čárkovanou čárou, jedná se o autokorelaci řádu n. Takže na příklad u letokruhů dochází k autokorelaci druhého řádu. Ostatní data vypadají bez problémů. Grafy vznikly následujícím způsobem: R> acf(x[1:34,3] # je potřeba specifikovat rozsah dat, pro další parametry R> # viz ?acf Zvláštním příkladem je analýza časových řad, při které se předpokládá autokorelace.

3.2

Metody početní - testy

Stanovíme si hypotézu H0 : Výběr pochází z normálního rozdělení.. Provedeme některé testy, které nám ji buď potvrdí nebo vyvrátí. Mezi použitelné testy patří např. KolmogorovSmirnovův test, Shapiro-Wilkokův test a další. V Rku najdeme například právě Shapiro-Wilkokův test normality. Ověříme si jej na našich datech: R> help.search(’normality’) R> ?shapiro.test R> shapiro.test(x[,1]) 10

Shapiro-Wilk normality test data: x[, 1] W = 0.9916, p-value = 0.7886 R> shapiro.test(x[,2]) Shapiro-Wilk normality test data: x[, 2] W = 0.6361, p-value = 2.616e-14 R> Nyní máme ověřeno i početně, co již tušíme graficky (viz obrázky 3.3 a 3.2). Všechny testy musíme číst následně: Porovnáváme hodnotu p-value (pravděpodobnost, že jev nastane) s hraniční testovou hodnotou, tzv. hladinou významnosti α, která bývá obvykle volena jako 0.05. α vyjadřuje pravděpodobnost, že jev nenastane. V tomto případě tedy můžeme tvrdit, že s pravděpodobností 1 − α, tedy 95%, se jedná o normální rozdělení.

11

12

4 Odhad parametrů základního souboru Parametry základního souboru jsou odvozeny z parametrů výběrového souboru. Jedná se o bodový a intervalový odhad.

4.1 4.1.1

Bodový odhad Odhad střední hodnoty

¯ = µ. Tedy že odhad středí hodnoty se rovná průměru. Odhad středí hodnoty E(x) Odhad rozptylu se spočítá jako S2 ·

n = σ2 n−1

(4.1)

n je korekce vychýlení. kde S 2 je rozptyl výběru a n−1 Charakteristiky základního souboru σ, µ jsou nazývány parametry, charakteristiky výběrového souboru S, x ¯ jsou nazývány statistikami. Vlastnosti bodových souborů:

• nespornost • nevychýlenost • vydatnost

4.2 4.2.1

Intervalový odhad Střední hodnota

Pro interval spolehlivosti T1 , T2 pro parametr τ při hladině významnosti α ∈ (0, 1) platí P (T1 ≤ τ T2 ) = 1 − α

(4.2)

kde α je pravděpodobnost, že parametr v intervalu ležet nebude. Parametr bude ležet v daném intervalu s pravděpodobností P = 1 − α. Zatímco pro charakteristiku výběrového souboru stačí bodové odhady, pro jejich vztažení na základní soubor potřebujeme intervalové odhady. Pro symetrické rozdělení je intervalový odhad symetrický, pro nesymetrické nesymetrický. Jednostranné odhady se používají při šetření o překročení nějaké hraniční hodnoty. Existují v dva základní přístupy k určení intervalového intervalu: 13

1. Pro velký počet prvků, což se obvykle bere jako > 30. Odhad se provede z kvantilu normálního rozdělení. σ σ ¯ + zα/2 · √ x ¯ − zα/2 · √ ≤ µ ≤ x n n

(4.3)

2. Pro malý počet prvků, (< 30) se používá Studentovo T-rozdělení. S S ¯ + t(α/2,n−1) · √ x ¯ − t(α/2,n−1) · √ ≤ µ ≤ x n n Veličina

4.2.2

x ¯−µ √ S/ n

(4.4)

má T rozdělení s k = (n − 1) stupni volnosti.

Faktory ovlivňující interval spolehlivosti

Mezi faktory ovlivňující velikost výběru patří zejména • velikost výběru • hladina významnosti α • variabilita • použitý vzorec Rozdělení T funguje vždy.

4.2.3

Interval spolehlivosti směrodatné odchylky σ

1. Pro malé výběry < 30 prvků platí rozdělení χ2 , které je nesouměrné. v u u n · S2 t ≤σ≤ 2

χα/2

v u u n · S2 t

χ21−α/2

(4.5)

2. Pro malé výběry > 30 prvků platí S σ = S ± zα/2 · √ 2n

Jak na to v Rku Zjišťování základních charakteristik souboru: R> # standardní odchylka S sloupce ’studny’ R> sd(x[,1]) [1] 7.596336 R> R> # rozptyl S^2 R> var(x[,1]) [1] 57.70432 R> 14

(4.6)

R> # středí hodnota R> mean(x[,1]) [1] 44.715 R> R> # počet prvků R> length(x[,1]) [1] 100 R> Podle rovnice 4.3 se interval spolehlivosti pro µ spočítá R> # z-kvantil z(alpha/2) dostaneme příkazem R> qnorm(p=0.975) [1] 1.959964 R> R> # vlastní interval je R> qnorm(p=0.975) * sd(x[,1])/sqrt(length(x[,1])) [1] 1.488854 R> α je 0.05. Protože funkce qnorm() počítá s pravděpodobností a ne s αou, a protože pravděpodobnost P = 1 − α, tak 1 − α/2 = 0.975, za předpokladu, že α = 0.05. Stejného výsledku bychom dosáhli obráceně: R> qnorm(p=0.025) * sd(x[,1])/sqrt(length(x[,1])) [1] -1.488854 Takže interval spolehlivosti pro µ je 44.715 ± 1.488854. Pro případy, kdy pracujeme s méně, než 30 prvky, pracujeme se Studentovým T rozdělením (rovnice 4.4). R> # hodnota pro T rozdělení potřebuje znát počet stupňů volnosti R> # což je vždy n-1: R> qt(p=0.975,df=length(x[,1])) [1] 1.983972 R> R> # a nyní interval spolehlivosti: R> qt(p=0.975,df=length(x[,1]))*sd(x[,1])/sqrt(length(x[,1])) [1] 1.507091 R> qt(p=0.025,df=length(x[,1]))*sd(x[,1])/sqrt(length(x[,1])) [1] -1.507091 Při výpočtu tiše předpokládáme, že prvků je méně, než 30. Důvod je ten, že při 30-ti prvcích se hodnoty T a normálního rozdělení začínají vyrovnávat. Pro výpočet intervalu σ postupujeme stejně. Nejdříve pro méně, než 30 prvků (rovnice 4.5): R> # hodnota chi^2 pro n-1 stupňů volnosti a prst alfa/2: R> qchisq(p=0.05/2, df=length(x[,1])-1) [1] 73.36108 15

R> R> #hodnotu chi^2 můžeme uložit do proměnné ’chi.t’ R> chi.t <- qchisq(p=0.05/2, df=length(x[,1])-1) R> chi.t [1] 73.36108 R> R> # a proměnnou ’chi.t’ můžeme hned použít při výpočtu horní hranice: R> sqrt(length(x[,1])*var(x[,1])/(chi.t)) [1] 8.868931 R> R> # spodní hranice se vypočítá podobně: R> chi.t <- qchisq(p=(1-0.05/2), df=99) R> chi.t [1] 128.422 R> sqrt(length(x[,1])*var(x[,1])/(chi.t)) [1] 6.703235 Takže s pravděpodobností 95% je σ ∈ (6.703235, 8.868931). Pro velký počet prvků (30 <) se postupuje podobně podle rovnice 4.6: R> hranice <- qnorm(p=0.025)*sd(x[,1])/sqrt(2*length(x[,1])) R> hranice [1] -1.052779 R> R> # horní hranice R> sd(x[,1]) - hranice [1] 8.649115 R> R> # dolní hranice R> sd(x[,1]) + hranice [1] 6.543557 σ je tak s pravděpodobností 95% ∈ (6.543557, 8.649115). Všimněme si, že intervaly odvozené pomocí normálního i podle rozdělení χ2 se skoro rovnají.

16

5 Transformace Na transformaci je dobrá funkce box.cox(), ale zatím jsem nepřišel na to, jak na ní.

17

18

6 Práce s malými výběry Pro práci s malými výběry, tedy výběry do deseti prvků, se používají jiné funkce. Pro malé výběry lze těžko ověřit předpoklady normality a homogenity. Používá se Hornův test. Používají se tzv. pivoty. Mějme 10 prvků: 5.8, 6.3, 6.4, 6.5, 6.7, 6.7, 7, 7, 9.2, 16.7 H=

int( n+1 int 10+1 2 + 1) 2 +1 = =3 2 2

(6.1)

Určí se horní a dolní pivot, což jsou třetí prvky (podle 6.1) se shora a ze zdola, takže 6.4 a 7. V dalším kroku se počítá pivotová polosuma: PL = (6.4 + 7)/2 = 6.7

(6.2)

RL = (7 − 6.4) = 0.6

(6.3)

A pivotové rozpětí

Nakonec interval spolehlivosti PL − RL · tL, 1−α ,n ≤ µ ≤ PL − RL · tL, 1−α ,n 2

2

Používá se speciální testovací statistika t, což není žádné běžné rozdělení. Google mlčí. . .

19

(6.4)

20

7 Výběry Vlastnosti, které by měl každý representativní výběr mít, jsou zejména 1. homogenní výběr 2. normální rozdělení 3. prvky jsou vzájemně nezávislé 4. všechny prvky základního souboru mají stejnou pravděpodobnost, že budou vybrány

7.1

Metody výběru

Jednoduchý výběr je nejpoužívanější pro soubory se známým a relativně malým počtem prvků. Systematický zvolí se náhodný počátek a náhodná velikost kroku. Následně se systematicky vyberou prvky. Oblastí výběr pokud lze soubor rozlišit na oblasti. Tento výběr se ještě dělí na úměrný kde počet prvků je úměrný velikosti oblasti a variabilní kde se zohledňuje variabilita jednotlivých oblastí. Vícestupňový . . . Vícefázový kdy se používají ještě jiné metody.

7.2

Zjišťování velikosti výběrů

Snažíme se odpovědět na otázku, jakou minimální velikost výběru potřebujeme vzhledem k účelu analýzy a požadované vypovídací hodnoty základního souboru? 1. Jaká bude přesnost odhadu? 2. Jakou požaduji spolehlivost, ze skutečné vzdálenosti mezi odhadem a skutečným parametrem bude menší než D, tedy spolehlivost odhadu? 3. Jaká je variabilit základního souboru? Ta se získá pomocí rozptylu S 2 nebo variačního koeficientu S. 21

Příklad 1. 1. Přesnost odhadu: tloušťka s přesností 3 cm, tedy ±1 cm → D = ±1. 2. Spolehlivost odhadu 95% → α = 0.05. 3. S = 18.04 4. n =? Jakou velikost výběru potřebují, abych zjistil, že µ bude x ¯ ± D? S2 D2 S% n = t2α/2;n−1 ∗ 2 D n = t2α/2;n−1 ∗

(7.1)

kde t – kvantil Studentova rozdělení pro hladinu významnosti α. Pokud je výběr větší, než 30, lze použít normální rozdělení zα/2 . Postupuje se aproximačně: 1. V prvním kroku se odhadne velikost výběru n a z toho se určí testovací statistika pro Studentovo rozdělení tα/2;n−1 . 2. Spočítá se podle vzorce 7.1 přesnější velikost výběru. Celý postup se opakuje, dokud se n mění. Pokud se očekává více, než 30 prvků, postupuje se okamžitě bez iterací. R> n_odh <- 5 # odhad R> D <- 1.5 # přesnost R> S2 <- 18.04 # směr. odchylka R> n <- 0 # výsledný počet prvků R> R> # výpočet provedeme automaticky pomocí smyčky typu ’for’ R> # ’tiše’ předpokládáme, že požadovaná přesnost bude dosažena v pěti krocích R> R> for(i in c(1:5)) { + n <- qt(p=1-0.025, df=n_odh-1)^2*S2/D^2 + print(n); + n_odh<-n; +} [1] 61.80622 [1] 32.06317 [1] 33.34537 [1] 33.23870 [1] 33.24723 R> R> # stejná úloha při použití normálního rozdělení: R> qnorm(p=1-0.025)^2*S2/D^2 [1] 30.79996 R> 22

Příklad 2. S jakou přesností jsem schopen určit interval odhadu střední hodnoty, když mám změřeno 15 stromů? Počet prvků n získáme z rovnice 7.1: D=

S · tα/2 ; n − 1 √ n

R> alpha <- 0.05 R> # směrodatná odchylka je známa z předchozího příkladu R> sqrt(S2)*qt(p=1-alpha/2, df=15-1)/sqrt(15) [1] 2.352105 R> Výslednou střední hodnotu jsem schopen určit s přesností ±2.35cm.

23

(7.2)

24

8 Testy Na základně údajů zjištěných z náhodného výběru ověřujeme, má-li být hypotéza potvrzena nebo zamítnuta. Pro každý test tedy musí být formulovány dvě hypotézy: H0: Nulová hypotéza. Např. že střední hodnota je rovna 50. Předpokládáme že platí, dokud test neprokáže něco jiného. H1: Alternativní hypotéza. Jedná se o opak. Platí, pokud je H0 zamítnuta. Hypotézy jsou buď jedno- nebo dvou stranné. O jednostrannou hypotézu se jedná, jde-li o překonání nějaké hranice (≤, ≥). O oboustrannou hypotézu se jedná, jde-li o rozhodnutí, že se něco rovná něčemu (=).

8.1

Testové kritérium

Je to obor možných hodnot. Hranici pro obor přijatých čísel je kritický bod, který je určen na základě α.

8.2

Postup testu

1. Formulace nulové hypotézy a alternativní hypotézy 2. Zvolíme hladinu významnosti α (chyba prvního druhu) 3. Zvolíme druh testu, v závislosti na rozdělení. 4. Určíme kritický bod, na jeho základě rozhodneme o přijetí H0 nebo jejím zamítnutí (obr. 8.1.). Výsledek testu je signifikantní na hladině významnosti α, pokud jeho výsledek vede k zamítnutí hypotézy.

Příklad Srovnání měření dvou výškoměrů. Máme 15 měření jedné výšky. Hodnota výšky je (zjištěno pásmem) 20 m. Průměr všech výšek je 19.2 m a směrodatná odchylka S=1.1 m. Měří výškoměr správně? 1. H0: Pro výběrový soubor o parametrech x ¯, S 2 platí, že µ0 = µ. H1: Výšky neměří správně. . . 25

obsah = 1 obsah = alpha/2

Obor pøijetí

obor zamítnutí

krit

krit

Obrázek 8.1: Znázornění kritického bodu a chyby prvního řádu

2. Stanovení hladiny významnosti, že test rozhodne špatně α = 0.05. Pravděpodobnost zamítnutí hypotézy je 1 − α. 3. Výběr testu. Vybereme test střední hodnoty pro jeden výběr: jednovýběrový t-test. Testovací statistika má tvar √ n−1 tkrit = (¯ x − µ0 ) (8.1) S Což je v podstatě upravená rovnice pro intervalový odhad 6.4. √ n−1 = −2.72 t = (19.2 − 20) 1.1 4. Stanovení kritické hodnoty qt(1-0.05/2, 14) = 2.145 5. Posouzení: Je kritická hodnota <> než testovací statistika? 2.145 < 2.72 → t − krit < t → H0 se zamítá: Odchylka naměřených hodnot je statisticky významná. nebo

R> dt((19.2-20)*sqrt(14)/1.1,df=14) # pt(t,n-1) = p-value [1] 0.008277159

0.008 < 0.025 → p − value < α/2 → H0 se zamítá.

S jakou pravděpodobností by se nezamítla? 1-pt((19.2-20)*sqrt(14)/1.1,df=14), 99%.

8.3

S pravděpodobností 1-p-value:

P-hodnota

Nebo taky p-value je pravděpodobnost, že získáme stejné nebo vyšší kritérium než vypočítané za předpokladu, že platí H0. Viz také obrázek 8.1. P-hodnotu dostaneme v Rku funkcí dt(kritická hodnota, počet stupňů volnosti) pro Studentovo a dnorm(kritická hodnota) pro normální rozdělení. Z předchozího příkladu bereme vstupní hodnoty: 26

R> pt(2.72,df=14) [1] 0.9834754 (V sešitě se to počítalo poněkud jinak, nevyšla hodnota 0.98, ale 1-0.98) S pravděpodobností 98.34% zamítáme hypotézu, že výškoměr měří správně. Hodnota p závisí zejména na • rozdíl mezi testovanou a skutečnou hodnotou • na velikosti výběru (protože stupeň volnosti je roven n − 1) Srovnávat má smysl pouze stejné velikosti výběru. Statistická významnost (p) Praktická důležitost pozoro- nevýznamná významná vaného rozdílu H0 − H1 nedůležitý n je OK n je velké důležitý n je malé n je OK Např. rozdíl mezi průměrnou výškou dvou porostů je 2mm. Před pokusem je potřeba udělat analýzu síly testu. P-hodnota může jen v omezené míře sloužit ke srovnání dvou výsledků a musíme rozlišovat mezi statistickou a reálnou významností testu. Buď bychom dostali zbytečně přesný výsledek (nákladný sběr dat) nebo naopak neprůkazný.

Příklady Příklad 1. Při zkoušce správnosti měření novým typem výškoměru byla 30 x změřena kontrolní výška, jejíž správná hodnota byla 20 m. Posuďte, zda výškoměr měří správně. Jedná se o jednovýběrový t-test. Testujeme hypotézu H0: x ¯ = 20: R> R> R> R> R>

# načteme data vysky <- read.table("vyska.txt")

# střední hodnota mean(vysky) V1 20.22 R> R> # jednovýběrový t-test R> ?t.test R> t.test(vysky, alternativ="two.sided", mu=20) One Sample t-test data: vysky t = 1.7514, df = 29, p-value = 0.09045 alternative hypothesis: true mean is not equal to 20 95 percent confidence interval: 27

19.96309 20.47691 sample estimates: mean of x 20.22 R> Ve výpisu vidíme několik důležitých řádků: Testovací statistika t, počet stupňů volnosti df a p-hodnota p-value. t = 1.7514, df = 29, p-value = 0.09045 95-ti procentní interval spolehlivosti pro data vysky: 19.96309 20.47691 A průměr dat 20.22. Spočítali-li bychom si kritickou hodnotu tkrit (viz. např. strana ??): R> qt(1-0.05/2,29) [1] 2.045230 Můžeme nyní na základě všech těchto údajů říct: 1. p-hodnota ¿ α/2 (0.09 ¿ 0.025) 2. tkrit > t (2.054 ¿ 1.75) 3. interval spolehlivosti µ ∈ (19.96309; 20.47691), že H0 se nezamítá na hladině významnosti α = 0.05. Tedy že rozdíl není statisticky významný.

Příklad 2. Variabilita teploty vzduchu v laboratoři je stanovena směrodatnou odchylkou S0 = 3◦ C. Při kontrole laboratorních podmínek bylo provedeno 20 měření teploty a ze zjištěných údajů byla vypočtena výběrová směrodatná odchylka S = 3,27 ◦ C. Je třeba posoudit, zda tato hodnota nesignalizuje zvýšení variability teploty vzduchu. Směrodatná odchylka má rozdělení χ2 . Posuzujeme jednostranně, zda-li není směrodatná odchylka S0 větší, než určitá hodnota. Hodnota χ2 se vypočítá podle vzorce χ2 = tedy

(n − 1) · S 2 , σ02

(20−1)·3.272 32

R> (20-1)*3.27^2/3^2 [1] 22.5739 Stanovení kritické hodnoty pro rozdělení χ2 : 28

(8.2)

R> qchisq(1-0.05, 19) [1] 30.14353 R> R> # p-hodnota R> pchisq(30.14353, 19) [1] 0.95 (V sešitě to mám p-hodnotu jako 0.256) Závěr: p-hodnota ¿ α, t ¡ tk rit → H0 se nezamítá, kolísání směrodatné odchylky je statisticky nevýznamné.

Příklad 3. Ve dvou porostech byly měřeny výčetní tloušťky v cm. Posuďte, zda jsou oba porosty stejně tloušťkově vyspělé. Jedná se o dvouvýběrový t-test. Ještě před vlastním testem musíme rozhodnout, mají-li výběry stejné rozptyly či ne. K tomuto testu se používá tzv. F-test. Nulová hypotéza říká, že rozptyly jsou shodné: R> R> R> R> R>

# načteme data porosty <- read.table("porosty1.txt",header=T,sep="\t") # f-test var.test(vysky[,1],vysky[,2]) F test to compare two variances

data: vysky[, 1] and vysky[, 2] F = 2.4411, num df = 24, denom df = 19, p-value = 0.05114 alternative hypothesis: true ratio of variances is not equal to 1 95 percent confidence interval: 0.9954193 5.7247255 sample estimates: ratio of variances 2.441087 R> Testovací statistika F se rovná 2.4411, Fkrit má hodnotu 2.441087, p-hodnota má hodnotu 0.05112/2 = 0.02556. (Proč děleno 2?). Takže p-hodnota ¡ α, F ¿ Fk rit → hypotéza H0 se zamítá a rozdíl je statisticky významný. Nyní můžeme přikročit k vlastnímu t-testu. Jedná se tedy o oboustranný párový t-test s nestejnými rozptyly: R> t.test(vysky[,1],vysky[,2], alternative = c("two.sided"), paired = TRUE) Paired t-test data:

vysky[, 1] and vysky[, 2] 29

t = 5.2608, df = 19, p-value = 4.455e-05 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: 4.34148 10.07852 sample estimates: mean of the differences 7.21 R> T = 7.21, Tkrit = 5.2608 → T > Tkrit , p − hodnota = 0, α = 0.05, α > p − hodnota. Zamítáme hypotézu H0, rozdíl středních hodnot je statisticky významný.

Příklad 4. Při zkoušce přesnosti výškoměrů byly při 15 kontrolních měřeních získány následující hodnoty. Posuďte, zda oba výškoměry měří stejně. Jedná se o párový t-test s nezávislými výběry. Zatímco v předchozím případě byla data sbírána ve dvou různých porostech na zcela různých stromech, v tomto případě porovnáváme stejné výšky dvěma různými výškoměry. Nulová hypotéza říká, že rozdíl středních hodnot je nulový. R> vyskomery <- read.table("poznamky/data/vyskomery1.txt",header=T,sep="\t") R> t.test(vyskomery[,1],vyskomery[,2], alternative = c("two.sided"), paired = TRUE) Paired t-test data: vyskomery[, 1] and vyskomery[, 2] t = -1.5671, df = 14, p-value = 0.1394 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -0.33160704 0.05160704 sample estimates: mean of the differences -0.14 R> P-hodnota ¿ α, interval spolehlivosti pro rozdíl je (-0.33160704, 0.05160704): Rozdíl je statisticky nevýznamný.

Příklad 5. Testy shody testují hypotézu, že experimentální data jsou ve shodě s nějakým rozdělením. Jedná se o jedno a vícevýběrové testy. V porostu borovice bylo změřeno 98 výšek a roztříděno do 2 m tříd. Řešila se otázka, zda pro rozdělení výšek v porostu je vhodným teoretickým modelem normální rozdělení. Na to se používá Kolmogorov-Smirnovův test. V Rku je na to funkce ks.test(). 30

Ručně to lze spočítat tak, že se data rozdělí do tříd a následně se vypočítá testovací statistika a kritická hodnota 1 max|Ne1 i − Ne2 i | n r 1 1 α D1,α = √ · − ln n 2 2

(8.3)

D1 =

(8.4)

V souboru shoda1.txt jsou data normálního rozdělení a data z měření. Data nelze přímo použít pro ks.test() protože jsou zapsána v distribučním pořadí a ne jako jednotlivé hodnoty.

Příklad 6. Testování, zda-li se dvě rozdělení rovnají lze provést pomocí rovnic pro testovací statistiku a kritickou hodnotu D2 = max|We1 i − We2 i |D2,α,krit =

1 α − ln 2 2

r

s

n1 + n2 n1 · n2

(8.5)

Soubor s daty je shoda2.txt. A zadání úkolu zní Ve dvou porostech byly měřeny výčetní tloušťky ve 4-cm tloušťkových třídách. Zjištěné počty v jednotlivých tloušťkových třídách jsou uvedeny v tabulce. Posuďte, zda lze rozdělení tlouštěk v obou porostech považovat za shodné.

8.4

Neparametrické testy

Doposud jsme používaly parametrické testy, které testovali parametry σ, µ,. . . Tyto testy fungují pouze za předpokladu splnění podmínek pro normální rozdělení. Co dělat, když tyto podmínky splněny nejsou? Použijí se neparametrické testy. Ty testují shodu rozdělení, nejsou na ničem závislé, ale mají menší sílu testu. Spíš se H0 nezamítne. Každý test má svou neparametrickou obdobu: T-test – Mann-Whiteyův Párový t-test – Willcoxův. Pokud to lze, měly by se používat parametrické testy.

31

32

9 Chyby testování a síla testu Chyba I druhu je pravděpodobnost, že test zamítne hypotézu, která skutečně platí α Chyba II druhu je pravděpodobnost, že test nezamítne hypotézu, která ve skutečnosti neplatí β. Síla testu je pravděpodobnost, že test správně zamítne hypotézu, která skutečně neplatí 1 − β. skutečnost testy platí neplatí platí správně p = 1 − α chyba II druhu β neplatí chyba I druhu α správně s = 1 − β Pokud test nezamítne nulovou hypotézu, může se dopustit chyby druhého druhu. Chybu druhého druhu nelze určiti před testem, protože je určena již chyba druhého druhu. β = Prst (t < C|H0neplat) β je pravděpodobnost, že (t-test ¡ kritická hodnota, za předpokladu, že H0 neplatí). C je kritický bod a t je hodnota testovacího kritéria.

Příklad H0: µ = 60 H1: µ = 65 n = 100 σ = 20 Stanovení kritického bodu: σ C = µ0 + zα · √ = n 20 = 60 + 1.1645 · = 63.29 10

(9.1) (9.2)

R> 60 + qnorm(0.95)*20/sqrt(100) [1] 63.28971 > Výpočet pravděpodobnosti odpovídající hodnotě β: β = Prst



c−µ x ¯−n √ < √ 1 σ n σ n 33



(9.3)

Standardizací dostaneme z=

x−µ σ

(9.4)

což je převod na standardizované normální rozdělení. Prst



x ¯ − µ1 C −µ √ < √ 1 σ n σ n



= Prst

63, 29 − 65 Z< 20/10





= P (Z < −0.855) = 0.1963

(9.5)

(Rovnice 9.1 až 9.5 bych asi potřeboval vysvětlit ještě jednou.) Síla testu závisí na 1. Velikosti efektu – jak daleko je od sebe H0 a H1 2. Variabilitě ZS 3. α+ → β− α− → β+ 4. Typ testu (parametrický/neparametrický)

9.1

Analýza síly testu

1. Apriorní analýza – před vlastním provedením pokusu. Zjišťuje se, jak veliký musí být výběr n, když známe α a sílu testu (1 − β) a známe velikost efektu, který chceme testem detekovat (např. přesnost přístroje). 2. Aposteriorní analýza – po měření. Zjišťuje skutečnou sílu testu 1 − β, s jakou pravděpodobností bude test schopen zamítnout hypotézu, které neplatí. Napevno se určí vždycky ta chyba, jejíž důsledky jsou závažnější.

9.1.1

Postup analýzy síly testu

1. Formulace problému a určení testovaných hypotéz 2. Určení optimálního uspořádání experimentu 3. Stanovení variability ZS. Určení kritické velikosti efektu, stanovení požadovaného poměru chyb α : β 4. Upřesnění experimentu: stanovení přijatelných hodnot α β. 5. Stanovení potřebné velikosti výběru na základě údajů stanovených v bodech 3 a 4. 6. Je stanovená velikost výběru přijatelná? Pokud ne a nižší hodnoty α a β jsou nežádoucí, musíme se vrátit zpět k bodům 1 a 2. Pokud nižší hodnoty akceptovatelné jsou a nebo pokud je stanovená velikost výběru přijatelná, pokračujeme dalším krokem: 7. Naměření potřebných dat a provedení testu nulové hypotézy 34

8. Je-li nulová hypotéza zamítnuta máme výsledek testu. Je-li potvrzena, musíme provést post hoc (následnou) analýzu síly testu. 9. Pokud je síla testu přijatelná, máme výsledek. Pokud není, musíme si opět položit otázku, je-li možné akceptovat nižší hodnoty α a β (bod 6).

Příklad 1. Byl sledován obsah v pitné vodě. Povolený obsah chlóru je 0.3 mg.l−1 . Určete, zda norma není překročena. Dále chceme zjistit, kolik vzorků je potřeba odebrat, aby možná chyba testu nepřesáhla 5%. Předběžně bylo odebráno 23 vzorků (soubor chlor1.txt): 0.1 0.15 0.25 0.15 0.3 0.25 0.25 0.3 0.35 0.55 0.7 0.7 0.25 0.2 0.15 0.65 0.55 0.5 0.3 0.35 0.3 0.25 0.8 Tato hodnota je typický jednostranný t-test, protože nás zajímá pouze překročení normy. Formulace hypotéz je tedy následující: H0: Cl ≤ 0.3 H1: Cl > 0.3 Chyba I druhu znamená zamítnutí platné nulové hypotézy. V tomto případě bychom tedy řekli, že obsah chlóru je překročen, i když by byl v pořádku. To by znamenalo další náklady na čištění vody. Při chybě II druhu by ve vodě byl vyšší, než povolený obsah chlóru a přitom by test doporučil nulovou hypotézu nezamítnout. Zde jsou následky vážnější. Testovat budeme klasickým t-testem. Několik základních údajů o vzorcích: R> chlor <- read.table("poznamky/data/chlor1.txt",) R> summary(chlor) V1 Min. :0.1000 1st Qu.:0.2500 Median :0.3000 Mean :0.3630 3rd Qu.:0.5250 Max. :0.8000 R> sd(chlor) V1 0.2023821 R> Testujeme hypotézu, že obsah chlóru není překročen: R> t.test(chlor,alternative = c("less"),mu=0.3) One Sample t-test data: chlor t = 1.4939, df = 22, p-value = 0.9253 alternative hypothesis: true mean is less than 0.3 95 percent confidence interval: -Inf 0.4355062 sample estimates: 35

mean of x 0.3630435 R> Hodnota testovacího kritéria t je tedy 1.4939 a kritická hodnota je qt(0.95,df=22) [1] 1.717144 Test tedy neumožnil zamítnou nulovou hypotézu a tak bychom se přiklonili k názoru, že voda je v pořádku. Nyní ovšem zjistíme sílu testu. V Rku jsou to tradičně moduly s označením power před názvem testu. power.t.test potřebuje pro svou práci následující argumenty: • počet prvků n, v našem případě 23 • velikost efektu delta, 0.3-0.363 = 0.063 • standardní odchylku, sd 0.202 • hladina významnosti α sig.level, 0.05 • sílu testu, power • typ testu type: dvouvýběrový, jednovýběrový, párový • alternativa: dvoustranný, jednostranný • ... Vynechán musí být jeden z argumentů n, delta, sig.level nebo power, který je z ostatních argumentů dopočítán. Nyní nás zajímá síla testu:

R> power.t.test(n=23,sd=sd(chlor),type="one.sample", delta=0.063, alternative="one.sided") One-sample t test power calculation n delta sd sig.level power alternative

= = = = = =

23 0.063 0.2023821 0.05 0.4215251 one.sided

R> Nyní již víme, že síla testu je 0.412 – tedy chyba II druhu je 0.579. Jinými slovy, test nesprávně zamítne nulovou hypotézu téměř v 60% případů. Test by téměř vždy ohlásil bezproblémový obsah chlóru, zatímco ve skutečnosti by to nebyla pravda. Musíme tedy určit velikost výběru takovou, aby se tato chyba minimalizovala. Při 95% síle testu potřebujeme 36

R> power.t.test(power=0.95,n=NULL,sd=sd(chlor),type="one.sample", delta=0.063, alternative One-sample t test power calculation n delta sd sig.level power alternative

= = = = = =

113.0484 0.063 0.2023821 0.05 0.95 one.sided

R> tedy 114 vzorků. To by bylo příliš drahé a tak, protože chyba I druhu pro nás nemá takové důsledky, snížíme její možnou úroveň na hodnotu 0.15 (v 15% případů připustíme, že test nesprávně zamítne nulovou hypotézu a bud nám „vnucovat“, že voda obsahuje příliš vysoký obsah chloru):

R> power.t.test(power=0.95,sig.level=0.15,n=NULL,sd=sd(chlor),type="one.sample", delta=0.0 One-sample t test power calculation n delta sd sig.level power alternative

= = = = = =

74.7377 0.063 0.2023821 0.15 0.95 one.sided

R> Výsledkem jsou tedy 75 měření. Nejčastější akceptovatelná síla testu je 80%. Na této úrovni začíná počet potřebných prvků prudce stoupat. Můžeme ještě použít kompromisní volbu poměru chyb vzhledem k jejich závažnosti, tedy β : α = 1 : 2 a zvolíme počet prvků např. 30. Tady to ale s Rkem tak lehce nejde. . . asi to chce pokusem a omylem (výsledek má být 30 měření).

Příklad 2. Mějme 20 měření tlouštěk dvou porostů. Chceme odpovědět na otázku, mají-li oba porosty stejné střední tloušťky. (Soubor vseobec data.txt, sloupce Tloušťka1 a Tloušťka2. R> obecdata <- read.table("vseobec_data.txt",header=T,sep="\t") R> summary(obecdata[,7:8]) Tloušťka1 Tloušťka2 Min. :17.60 Min. :20.70 1st Qu.:22.82 1st Qu.:23.40 Median :24.20 Median :26.10 37

Mean :24.61 Mean :26.34 3rd Qu.:26.48 3rd Qu.:29.10 Max. :30.90 Max. :32.20 NA’s :13.00 NA’s :13.00 R> sd(obecdata[1:20,7:8]) Tloušťka1 Tloušťka2 3.529559 3.364720 > Testujme nulovou hypotézu, že střední hodnoty jsou si rovné: 1. F-test, mají-li oba porosty stejnou směrodatnou odchylku: R> > var.test(obecdata[1:20,7],obecdata[1:20,8]) F test to compare two variances data: obecdata[1:20, 7] and obecdata[1:20, 8] F = 1.1004, num df = 19, denom df = 19, p-value = 0.837 alternative hypothesis: true ratio of variances is not equal to 1 95 percent confidence interval: 0.4355442 2.7800585 sample estimates: ratio of variances 1.100381 R> Jejich rozptyly se rovnají 2. Dvouvýběrový párový t-test s rovností rozptylů R> t.test(obecdata[1:20,7],obecdata[1:20,8],paired =T) Paired t-test data: obecdata[1:20, 7] and obecdata[1:20, 8] t = -1.5208, df = 19, p-value = 0.1448 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -4.099058 0.649058 sample estimates: mean of the differences -1.725 R> (Jiné hodnoty, než v souboru testy 1 a 2 výběry příklady.xls). Rozdíl v průměrech je statisticky nevýznamný. 38

Nyní je čas na analýzu síly testu. K tomu si nejdříve potřebujeme spočítat průměrnou směrodatnou odchylku σ: σ=

q

2 + σ 2 − 2Rσ σ σA A B B

(9.6)

Kde R je korelační koeficient obou výběrů.

(Poje poznámky, strana 27/z druhé strany. Korelační koeficient mi vychází cor(obecdata[1:20,7],obecdata[1:20,8]) [1] -0.08219599 což je něco úplně jiného, než se uvádí v souboru testy 1 a 2 . . .

9.2

Pro statistické testy o střední hodnotě

Pro situace, kdy není žádný statistický software k dispozici, následuje několik vzorců: Pro jeden výběr HO : µ = µ0 n=

2 s2  t + t α,n−1 β(1),n−1 σ2 s

δ=

 s2  tα,v + tβ(1),v n

(9.7) (9.8)

Pro dva výběry HO : µ1 = µ2 2 2s2p  n = 2 tα,v + tβ(1),v sδ

δ=

 2s2p  tα,v + tβ(1),v n

(9.9) (9.10)

Kde: n – je velikost výběru s2 – bodový odhad rozptylu základního souboru δ – požadovaná přesnost → µ0 − µ1 α, β – chyba I a II druhu. Chyba II druhu β má vždy jednostranné rozdělení. (2 × β).

9.3

Testy – přehled

Kromě T testů t.test() a jeho neparametrické varianty wilcox.test() disponuje R ve svém základu následujícími testy: binom.test friedman.test prop.test

– test parametru binomické distribuce – neparametrický test pro „complete block design“ bez opakování – test hody poměrů, výpočet konfidenčního intervalu pro pravděpodobnost úspěchu 39

– test χ2 pro dvourozměrnou kontingenční tabulku s možností tzv. Yatesovy korelace kruskal.test – neparametrický test pro jednoduché třídění cor.test – testuje hypotézu o nulové hodnotě korelačního koeficientu – Pearsonovou, Kendallovou nebo Spearmanovou metodou mantelhaen.test – Mantel-Haenzselův χ2 test pro trojrozměrné kontingenční tabulky var.test – F test srovnávající variance dvou výběrů barlett.test – porovnávání variance pro více než dva výběry fisher.tets – Fischerův exaktní test pro dvourozměrnou kontingenční tabulku mcnemar.test – χ2 test hypotéz o symetrii pro dvourozměrné kontingenční tabulky

chisq.test

40

10 Korelační a regresní analýza – Statistické modely Vzorce udávají závislou („vysvětlovanou“) a nezávislou („vysvětlující“) proměnné a charakter vztahu mezi nezávislými proměnnými. Jedná se o vícerozměrné statistiky a souvisí s anovou (více rozměrnou analýzou). Statistická závislost znamená, že pokud znaku x přísluší libovolné hodnoty všech ostatních znaků, nazýváme znaky x1 , x2 , . . . , xn statisticky nezávislé. Data, která nekorelují, jsou bez trendu. Stochastická závislost znamená, že existuje jakýsi trend. Pro hodnoty x existuje přesnější určení hodnoty y, než je x ¯. Hodnoty y můžeme přiřadit na základě určité pravděpodobnosti. Funkční závislost každé hodnotě x odpovídá právě jedna hodnota y. Korelace popisuje vliv změny úrovně jednoho znaku na změnu úrovně jiných znaků kvantitativních (měřených) dat. Kontingence popisuje vliv změny kvalitativních znaků (např. „druh“) Asociace popisuje vliv změny alternativních znaků (boolean).

10.1

Korelace

Podle počtu korelovaných znaků rozlišujeme Jednoduchou korelaci, kde se shodují 2 znaky Mnohonásobnou korelaci, kde zkoumáme 3 a více znaků Parciální , kde zkoumáme dva znaky ve vícerozměrném souboru s vyloučením ostatních znaků. Dále může být korelace kladná, záporná, lineární, nelineární, . . . Při korelační analýze zjišťujeme těsnost závislosti, existuje-li jaká. Při regresní analýze se snažíme vytvořit vhodný model a stanovit jeho parametry.

10.1.1

Míra korelační závislosti

Celkovou variabilitu tvoří odchylky od průměru. Dělí se na část vysvětlenou modelem a na část reziduální (zbytkovou, viz obr 10.1). Těsnost závislosti se hodnotí poměrem části vysvětlené modelem a částí nevysvětlenou (jako ANOVA). ANOVA se používá pro test významnosti 41

Reziduum

Cást vysvìtlená modelem

Obrázek 10.1: Princip regresní analýzy

regresního modelu. Jako míra korelace se používá koeficient determinace (rovnice 10.1), což je poměr variance vysvětlené modelem a celkovému rozptylu, tedy podíl celkové variability vysvětlené regresním modelem. R2 =

s2x′

2

s2x2

=1−

s2x1 x2 s2x2

(10.1)

Dále používáme koeficient korelace R=

10.1.2

v u 2 u sx′ t 2

s2x2

v u u s2 = t1 − x21 x2

sx2

(10.2)

Korelační koeficient

Pro jednoduchou korelaci se používá tzv. párový korelační koeficient. Je to zvláštní případ vícenásobného korelačního koeficientu, vyjadřuje míru lineární stochastické závislosti mezi náhodnými veličinami xi a xj . Máme dva: Pearsonův , pokud je splněna podmínka dvojrozměrného normálního rozdělení a Spearmanův , pokud není. Pro vícenásobnou korelaci se používá koeficient vícenásobný , který definuje míru lineární stochastické závislosti mezi náhodnou veličinou x1 a nejlepší lineární kombinací složek x2 , x3 , . . . , xm náhodného vektoru x a parciální , který definuje míru lineární stochastické závislosti mezi náhodnými veličinami xi a xj při zkonstantnění dalších složek vektoru x. 42

Pearsonův korelační koeficient r Používá se, mají-li obě proměnné normální rozdělení. covx1,x2 rx1 ,x2 = rx2 ,x1 = sx1 · sx2

(10.3)

kde covx1 x2 je normovaná kovariance. Kovariance

• je míra intenzity vztahu mezi složkami vícerozměrného souboru • je mírou intenzity lineární závislosti, • je vždy nezáporná, • její limitou je součin směrodatných odchylek • je symetrickou funkcí svých argumentů covx1 x2 =

2 1X (x1i − x1 )(x2i − x2 ) n i=1

(10.4)

Mezi základní včasnosti Pearsonova korelačního koeficientu patří • že je to bezrozměrná míra lineární korelace • nabývá proto hodnot od 0 - 1 pro kladnou a 0 - -1 pro zápornou korelaci • hodnota korelačního koeficientu je stejná pro obě proměnné Příklad Mějme datové soubory x1 a x2. Zjistíme jejich korelační koeficient: R> x1<-c(5,2,4,5,6,2,4) R> x2<-c(5,2,5,1,5,4,1) R> ?var R> cor(x1, x2,method="pearson") [1] 0.2309401 R> Korelační koeficient je tedy 0.23. Spearmanův korelační koeficient r Je neparametrický, nevychází z hodnot, ale z jejich pořadí. Používá se, pokud není splněna podmínka normality. 6 · ni=1 d2i rs = 1 − n3 − n P

Kde d2i je diference mezi pořadími hodnot x a y v jednom řádu. R> cor(x1, x2,method="spearman") [1] 0.3077492 R> 43

(10.5)

Mnohonásobný korelační koeficient Mezi základní vlastnosti mnohonásobného korelačního koeficientu patří: • 0‘R‘1 • pokud R = 1 znamená to, že závislá proměnná x1 je přesně lineární kombinací veličin x2 , . . . , xm . • pokud je R = 0, potom jsou také všechny párové korelační koeficienty nulové • s růstem počtu vysvětlujících (nezávislých) proměnných hodnota vícenásobného korelačního koeficientu neklesá, R1(2) ‘R1(2,...,m) . R1(2,...,m) =

s

1−

det(R) det(R(11)

(10.6)

kde det(R) je determinant korelační matice a det(R(11) ) je determinant korelační matice s vypuštěným sloupcem a řádkem odpovídajícím té proměnné, jejíž závislost na zbytku matice se vypočítává. Jako příklad si zkusíme vypočítat mnohonásobný korelační koeficient koncentrace kadmia v různých částech rostli. Zjistíme, zda-li obsahy kadmia v různých částech spolu navzájem korelují: R> kadmium<-read.table("kadmium.txt",sep="\t",header=T) R> kadmium [...] R> # korelační matici získáme funkcí cor() R> kadmium.cor<-cor(kadmium) R> kadmium.cor zrno otruby stonek kořen zrno 1.0000000 0.9836940 0.9934748 0.9948218 otruby 0.9836940 1.0000000 0.9934409 0.9869345 stonek 0.9934748 0.9934409 1.0000000 0.9884743 kořen 0.9948218 0.9869345 0.9884743 1.0000000 R> R> # Mnohonásobný korelační koeficient pro zrno: R> sqrt(1-det(kadmium.cor)/det(kadmium.cor[2:4,2:4])) [1] 0.9985795 R> R> # Mnohonásobný korelační koeficient pro stonek: R> sqrt(1-det(kadmium.cor)/det(kadmium.cor[c(1,2,4),c(1,2,4)])) [1] 0.9985228 R> Parciální korelační koeficient Používá se k posouzení síly závislosti dvou veličin ve vícerozměrném souboru při vyloučení vlivu ostatních veličin. Podle počtu vyloučených proměnných se stanovují řády parciálního R. Parciální korelační koeficient III. řádu má tedy tři vyloučené proměnné. 44

Výpočet vychází z korelační matice, ze které se počítají parciální korelace I. řádu s jednou vyloučenou proměnnou, z nich II. řádu (dvě vyloučené proměnné) atd. až do potřebného řádu. 

−1j · det R(ij)



Rij(1,2,...,m) = r     det R(ii) · det R(jj)

(10.7)

Na našem příkladu s kadmiem v různých částech rostliny, chceme zjistit parciální korelační koeficient pro zrno a otruby (1,2). Protože chceme vyloučit stonek a kořen, jedná se parciální derivaci II. řádu:

R> # korelační matice R> kadium.cor<-cor(kadmium) R> R> R11 <- kadmium.cor[2:4,2:4] R> R22 <- kadmium.cor[c(1,3:4),c(1,3:4)] R> R12 <- kadmium.cor[c(2:4),c(1,3:4)] R> R> # vlastní parciální korelační koeficient R> (-1^2*det(R12))/sqrt(det(R11)*det(R22)) [1] 0.7181145 R> Můžeme tedy říct, že mezi zrnem a otrubami je korelace 0.71 s vyloučením stonku a kořene.

10.2

Regrese

Základní úlohou regresní analýzy je nalezení vhodného modelu studované závislosti. Snažíme se nahradit každou meřenou (empirickou) hodnotu závislé proměnné y hodnotou teoretickou (modelovou), t.j hodnotou ležící na spojité funkci nezávislé proměnné x. Regresní model předpokládá, že nezávislá proměnná je nenáhodná a závislá proměnná je náhodná. Tento předpoklad nebývá v praxi striktně naplněn. Často je náhodně měřená proměnná i ta nezávislá. Pak se hovoří o tzv. korelačním modelu. Regresní modely jsou lineární , jejichž parametry jsou v vzájemném lineárním postavení a nelineární , jejichž parametry jsou v nelineárním postavení. 45

Regresní parametr 1

Absolutni clen Obrázek 10.2: Regresní model

10.2.1

Příklady regresních modelů

Matematická rovnice Lineární modely y = a + bx y = a + bx + cx2 y = a + b/x Nelineární modely x = a · xb ...

10.2.2

kód pro modul lm() v Rku a poznámky Rovnice přímky. y~1+x nebo y~ x Absolutní člen a může být nastaven na hodnotu 0 pomocí y~0+x Parabola. y~1+x+I(x^2) nebo y~poly(x,2). V Rku lze provést i regresi matice X: y~X + poly(x,2). Hyperbola. y~1 + 1/x.

Postup regresní analýzy

• Najít nejvhodnější tvar regresního modelu (tedy určit příslušnou rovnici, která bude popisovat závislost y na x). • Stanovit jeho parametry • Stanovit statistickou významnost modelu • Výsledky dané modelem interpretovat z hlediska zadání Stanovení vhodného tvaru modelu 1. Najít množinu modelů, které svými vlastnostmi vyhovují řešenému problému (např. růstové funkce) 2. Teprve mezi nimi najít podle statistických kritérií ten model, který nejlépe vyhovuje měřeným datům. 46

Jedna z nejčastěji používaných metod je metoda nejmenších čtverců. Řeší se rovnice X

(yi + yˆi )2 = min

(10.8)

Obecný vztah pro výpočet regresních parametrů b = (xT · x)−1 · xT · y

(10.9)

kde xT – x – −1 – y –

transponovaná matice x normální matice inverze závislá proměnná

Obecný vztah pro výpočet predikovaných (modelových) hodnot lineárního modelu yˆ = x(xT · x) · xT · y

(10.10)

kde vztah x(xT · x) · xT je projekční matice H.

10.2.3

Předpoklady metody

• Regresní parametr b může nabývat jakékoliv hodnoty • Regresní model je lineární v parametrech • nezávislé proměnné jsou skutečně vzájemně nezávislé, dochází ke kolinearitě • podmíněný rozptyl D(y/x) = σ 2 je konstantní (tzv. podmínka homoskedasticity). • ...

10.2.4

Multikolinearita

Každý sloupec by měl být nezávislý na ostatních. Multikolinearita způsobuje početní problémy , hlavně špatnou podmíněnost matice a problém při inverzi matice statistické problémy , jako že nelze odděleně sledovat vliv jednotlivých nezávislých proměnných na závislou, nespolehlivé parametry, nestabilita odhadů regresních parametrů – malá změna hodnot způsobí velkou změnu parametrů Příčiny multikolinearity Mezi příčiny paří zejména • přeurčenost regresního modul (zbytečně moc nezávislých proměnných) • nevhodné rozmístění experimentálních bodů • povaha modelu (např. polynomu) 47

Testování multikolinearity VIF – variance inflation factor – diagonální prvky inverzní matice ke korelační matici nezávisle proměnných. Kritická hodnota je 10. Pokud je hodnota větší, než 10, dochází k multikolinearitě. Postup 1. Vytvoří se korelační matice 2. Provede se její inverze 3. Sledují se hodnoty na diagonále. Kde je VIF větší, než 10, dochází k multikolinearitě. R> R> R> R> R> R> R> R> R>

# Můžeme použít buď funkci solve() nebo ginv() z balíku MASS, # použijeme standardní funkci solve(). Viz nápověda k funkci solve() # Použijeme data kadmium, které již máme načtena # Nejdříve korelační matice R<-cor(kadmium)

# zbývá i invertovat solve(R) zrno otruby stonek kořen zrno 352.2457 176.7276 -293.4094 -234.8127 otruby 176.7276 171.9395 -212.1790 -135.7720 stonek -293.4094 -212.1790 338.7247 166.4761 kořen -234.8127 -135.7720 166.4761 204.0375 R> Na příkladu vidíme, že všechna data mezi sebou vykazují značnou multikolinearitu.

10.2.5

Řešení multikolinearity

K odstranění nebo snížení multikolinearity lze použít snížení nezávislých proměnných, použití jiného modelu nebo použití metody PCR – principal component regression.

10.2.6

Homoskedasticita a Heteroskedasticita

Homoskedasticita znamená, že hodnoty závislé proměnné y mají všechny hodnoty nezávislé proměnné x konstantní rozptyl. Heteroskedasticita se testuje buď pomocí trendu reziduí nebo pomocí Cookova-Weisbergova testu. Nejobvyklejším řešením je použití metody vážených nejmenších čtverců, kdy se podmínka sumy reziduí násobí vhodně zvolenými váhami.

10.2.7

Postup regresní analýzy (lineární model)

Mějme soubor listy1.txt, obsahující dva sloupce: šířku a délku listu. Pokusíme se ukázat si několik analýz a výstupů které by neměly chybět jejich interpretaci. 48

Závislá promìnná

Závislá promìnná

Nezávislá promìnná

Nezávislá promìnná

Homoskedasticita

Heteroskedasticita

Obrázek 10.3: Homo- a heteroskedasticita

1. Vlastní regresní model je jednoduchý. Budeme zkoumat závislost šířky listů na jejich délce. Chceme použít lineární model a na to je v Rku funkce lm(). Zkrátka: R> # načteme data R> listy<-read.table("poznamky/data/listy1.txt", header=T, sep="\t") R> listy délka šířka 1 24 14 2 29 19 3 30 14 4 33 22 [...] R> R> # lineární model uděláme pomocí funkce lm() R> lm(listy$šířka ~ listy$délka) Call: lm(formula = listy$šířka ~ listy$délka) Coefficients: (Intercept) listy$délka -2.7521 0.6763 R> Mohlibychom tedy jednoduše uzavřít, že datům odpovídá model šířka = −2.752 + 0.6763 × délka. Nyní několik grafů (obrázky 10.4, 10.5):

49

R> R> R> R> R> R> R> R> R>

# graf prměnných, tak jak jsou plot(listy$délka,listy$šířka) výroba hodnot X a Y pro model xfit <- c(min(listy$délka):max(listy$délka)) yfit <- -2.752+0.6763*xfit # vytiskneme linii do grafu lines(xfit,yfit)

30 15

20

25

listy$šířka

35

40

Délka / Šířka listů

30

40

50

60

listy$délka

Obrázek 10.4: Graf regresní přímky a dat

R> R> R> R> R> R> R> R>

# nastavení rozdělení okna pro grafy na 4 části: par(mfrow=c(2,2)) # několik důležitých diagnostických grafů pomocí funkce lm() vnořené # ve funkci plot() plot(lm(listy$šířka~listy$délka)) par(mfrow=c(1,1))

Obrázek 10.5 ukazuje několik zajímavých analytických grafů dat. Jako první je graf reziduí okolo regresní linie modelu. Jejich hodnoty můžeme získat pomocí funkce residuals(objekt()), v našem případě tedy 50

Normal Q−Q plot

6

Residuals vs Fitted 19

9

1 0 −1

2 0 −4

−2

Residuals

4

Standardized residuals

2

19

9

−6

12 12

15

20

25

30

35

40

−2

−1

Fitted values

1

2

Theoretical Quantiles

Cook’s distance plot 0.5

Scale−Location plot 1.5

0

19

19

0.4 0.3 0.2

3 12

0.0

0.1

0.5

1.0

Cook’s distance

9

0.0

Standardized residuals

12

15

20

25

30

35

40

5

Fitted values

10 Obs. number

Obrázek 10.5: Několik důležitých diagnostických grafů

51

15

20

R> residuals(lm(listy$šířka~listy$délka)) 1 2 3 4 5 6 7 0.5197185 2.1379734 -3.5383756 2.4325774 3.4325774 -1.2437716 -2.2728187 8 9 10 11 12 13 14 3.3744833 -3.9782147 0.3454363 2.3163892 -5.0363088 -1.7126578 2.2873422 15 16 17 18 19 20 -1.3890068 -1.0944029 -2.4471009 -0.1234499 5.8475031 0.1421070 R> Číslem jsou označeny tzv. vlivné body, body vírazně ovlivňující model. Jsou to horcí kandidáti na vymazání ze souboru a tím by se mělo dosáhnout přesnějšího modelu. O vlivných bodech (snad) viz níže. Druhý graf je již známý Kvantil-kvantilový graf ovšem reziduí. Zajímavý je také graf Cookovy vzdálenosti. Cook-Weisbergův test se používá pro test heteroskedasticity (část 10.2.6.). Kromě grafů potřebujeme znát také nějaké číselné statistiky o lineárním modelu. Celkový přehled získáme pomocí funkce summary() R> summary(lm(listy$šířka~listy$délka)) Call: lm(formula = listy$šířka ~ listy$délka) Residuals: Min 1Q -5.036309 -1.852698

Median 0.009329

3Q 2.294604

Max 5.847503

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -2.7521 2.7293 -1.008 0.327 listy$délka 0.6764 0.0613 11.034 1.93e-09 *** --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.9 on 18 degrees of freedom Multiple R-Squared: 0.8712, Adjusted R-squared: 0.864 F-statistic: 121.7 on 1 and 18 DF, p-value: 1.925e-09 R> Postupně jsme dostali přehled kvantilových charakteristik reziduí. Dále jsme dostali charakteristiku jednotlivých parametrů modelu. (Předpokládáme formu rovnice y = ax + b). Parametr b (Intercept) má p-hodnotu 0.33, která je větší, než hodnota α 0.05. Hypotéza H0 zní: parametr a je roven 0. Můžeme tedy na hladině významnosti α = 0.05 parametr b vynechat. 52

P-hodnota parametru a (1.93 · 10−9 ) je poměrně nízká a my můžeme říci, že zamítáme nulovou hypotézu, že parametru a = 0 (což se jaksi předpokládalo. . . ). Hodnota koeficientu determinace (R2 ) 0.87 říká, že 87% všech hodnot je vysvětleno modelem a 13% jsou rezidua. F-test testuje koeficient determinace – test významnosti R. P-hodnota celého testu (1.925 · 10−9 ) nám říká, že test je významný. Celý test ještě budeme testovat pomocí anovy. R> anova(lm(listy$šířka~listy$délka)) Analysis of Variance Table Response: listy$šířka Df Sum Sq Mean Sq F value Pr(>F) listy$délka 1 1023.65 1023.65 121.75 1.925e-09 *** Residuals 18 151.35 8.41 --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 R> Anova testuje významnost regresního modelu. Stejně jako u anovy testuje celkový rozptyl, tedy rozptyl vysvětlený faktorem + rezidua. Významnost testu je podle p-hodnoty anovy poměrně vysoká (0.001). Protože p-hodnota ¡ α, zamítáme hypotézu H0, že průměry jsou si rovné – tedy v případě modelu, že koeficienty rovnice a a b jsou si rovné. 2. Ok, rozhodli jsme, že koeficient b je zbytečný, protože na základě testu se ukazuje, že by s pravděpodobností 0.95% mohl být 0: R> lm(listy$šířka~ 0 + listy$délka) Call: lm(formula = listy$šířka ~ 0 + listy$délka) Coefficients: listy$délka 0.6163 R> Model vychází tedy y = 0.6163 · x. A nyní tesy: R> summary(lm(listy$šířka~ 0 + listy$délka)) Call: lm(formula = listy$šířka ~ 0 + listy$délka) Residuals: 53

Min 1Q Median -4.9665 -1.9113 -0.2746

3Q 1.8129

Max 6.6378

Coefficients: Estimate Std. Error t value Pr(>|t|) listy$délka 0.61631 0.01457 42.3 <2e-16 *** --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.901 on 19 degrees of freedom Multiple R-Squared: 0.9895, Adjusted R-squared: 0.9889 F-statistic: 1790 on 1 and 19 DF, p-value: < 2.2e-16 R> R> anova(lm(listy$šířka~ 0 + listy$délka)) Analysis of Variance Table Response: listy$šířka Df Sum Sq Mean Sq F value Pr(>F) listy$délka 1 15060.1 15060.1 1789.6 < 2.2e-16 *** Residuals 19 159.9 8.4 --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 R> Kvalita modelu poněkud poklesla (např. viz suma koeficient determinace R2 ), ale model se zjednodušil. Podle testu anovy vychází celková významnost modelu pořád poměrně dobře. 3. Celkovou charakteristiku lineárního modelu získáme pomocí funkce profile(). Místo funkce lm() (lineární modely) použijeme ale funkci glm() (obecné lineární modely): R> profile(glm(listy$šířka~listy$délka)) $"(Intercept)" tau par.vals.(Intercept) par.vals.listy$délka 1 -4.558368 -15.1930477 0.9477663 2 -3.798640 -13.1195556 0.9025301 3 -3.038912 -11.0460634 0.8572939 4 -2.279184 -8.9725712 0.8120577 5 -1.519456 -6.8990791 0.7668214 6 -0.759728 -4.8255869 0.7215852 7 0.000000 -2.7520947 0.6763490 8 0.759728 -0.6786026 0.6311128 9 1.519456 1.3948896 0.5858766 10 2.279184 3.4683818 0.5406404 11 3.038912 5.5418739 0.4954041 12 3.798640 7.6153661 0.4501679 54

$"listy$délka" tau par.vals.(Intercept) par.vals.listy$délka 1 -3.798640 7.3185527 0.4435017 2 -3.038912 5.3044232 0.4900711 3 -2.279184 3.2902937 0.5366406 4 -1.519456 1.2761643 0.5832101 5 -0.759728 -0.7379652 0.6297795 6 0.000000 -2.7520947 0.6763490 7 0.759728 -4.7662242 0.7229185 8 1.519456 -6.7803537 0.7694879 9 2.279184 -8.7944832 0.8160574 10 3.038912 -10.8086127 0.8626269 11 3.798640 -12.8227422 0.9091964 12 4.558368 -14.8368717 0.9557658 attr(,"original.fit") Call:

glm(formula = listy$šířka ~ listy$délka)

Coefficients: (Intercept) listy$délka -2.752095 0.676349 Degrees of Freedom: 19 Total (i.e. Null); 18 Residual Null Deviance: 1175 Residual Deviance: 151.3458 AIC: 103.2342 attr(,"summary") Call: glm(formula = listy$šířka ~ listy$délka) Deviance Residuals: Min 1Q -5.036308792 -1.852698023

Median 0.009328567

3Q 2.294603955

Max 5.847503072

Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -2.75209474 2.72925573 -1.00837 0.32664 listy$délka 0.67634901 0.06129755 11.03387 1.9251e-09 *** --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for gaussian family taken to be 8.408098) Null deviance: 1175.00000 Residual deviance: 151.34577

on 19 on 18 55

degrees of freedom degrees of freedom

AIC: 103.23424 Number of Fisher Scoring iterations: 2 attr(,"class") [1] "profile.glm" "profile" R> Hlášek je již podstatně více. Kromě koeficientů a a b (Coefficients) a charakteristik reziduí zde máme i podstatnou charakteristiku Akaikovo informační kritérium AIC = 103.23. Tato statistika se používá při porovnávání různých modelů mezi sebou. Platí, že čím vyšší hodnota AIC, tím lepší model. Můžeme také porovnat druhý model, který zkoušíme, kdy b = 0: R> profile(glm(listy$šířka~ 0 + listy$délka)) [...] AIC: 102.33327 [...] R> A skutečně, víme, že tento model by měl být poněkud horší, než předchozí. 4. Ještě k obrázku 10.5, konkrétně ke grafu reziduí. Rozložení reziduí – jejich trend – ukazuje na kvalitu modelu. Má několik charakteristických tvarů rozložení reziduí (obr. 10.6). První graf ukazuje situaci tak, jak to má být. Rezidua bez výrazného trendu.

Obrázek 10.6: Různé tvary rozložení reziduí Druhý graf ukazuje rezidua tak, jak by vypadaly při použití nevhodného modelu. Některé hodnoty „prostě nesedly“. Třetí graf ukazuje situaci, kdy závislá proměnná nemá konstatntí rozpty. Nastává heteroskedasticita – porušení jedné z podmínek pro metodu nejmenších čtverců. Jedním z řešní je použít váženou metodu nejmenších čverců a jako váhu lze použít např. hodnotu 1/y 2 . Tím se oslabí váha hodnot s velkou variabilitou. 5. Interval spolehlivosti parametrů lineárního modelu lze spočítat jednoduše pomocí funkce confint() z balíku MASS. To je asi poprvé, co použijeme nějaký jiný balík, než standardní funkce z Rka. Do R přispívá mnoho lidí svými vlastními funkcemi a tím je R obohacováno. Knihova MASS slouží k různým pokročilým statistickým analýzám. Nejdříve 56

ji musíme nainstalovat a následně natáhnout do R. Poněkud matoucí je, že knihovna MASS se skrývá v balíku „VR“: R> install.packages("VR") [...] # hlášky o stahování a instalaci R> # natažení knihovny MASS R> library(MASS) R> ?confint R> confint(lm(listy$šířka~listy$délka)) 2.5 % 97.5 % (Intercept) -8.4860483 2.9818588 listy$délka 0.5475676 0.8051304 R> Vidíme, že interval spolehlivosti pro parametr b (y = ax + b) obsahuje 0 a tak se bez něj s 95% pravděpodobností obejdeme (stejně tak jsme rozhodli i výše). 6. Ještě musíme udělat intervalový odhad všech prvků a modelu jako celku. Ten se počítá pomocí Fisherovy transformace, která má přibližně normální rozdělení: 1 z(R) ± z1−α/2 · √ n−3

(10.11)

kde z je kvantil normovaného normálního rozdělení. Musí se určit horní a dolní hranice a následně provést retransformace zpět. Interval spolehlivosti je nesymetrický. Nám ale postačí Fisherova testovací statistika z bodu 2. 7. Dále je potřeba udělat intervalový odhad spolehlivosti jednotlivých hodnot. K tomu slouží rovnice (pro model přímky) σ My′ = y1 ± tα/2;n−2 · √ · n−2 ′′

s

Kde

n(xi − x ¯)2 1 + Pn ¯)2 i=1 (xi − x

(10.12)

tα/2 – směrodatná odchylka yi′ – modelovaná hodnota σ – směrodatná odchylka reziduí 8. A neškodil by ani interval spolehlivosti odhadovaných hodnot yi′ ± tα/2;n−m σ

(10.13)

kde m je počet parametrů (např. pro rovnici y = ax + b platí, že m = 2. Oba poslední body získáme následovně: Nejdříve si uložíme model do nové proměnné lm.vs (lineární model výška-šířka). R>

lm.vs<-(glm(listy$délka ~ listy$šířka)) 57

Kromě toho, že nyní můžeme opět vytisknout charakteristické grafy a získávat souhrny o takto vytvořeném modelu R> plot(lm.vs) R> summary(lm.vs) [...] Můžeme nyní jednoduše spočítat intervalové odhady pro jednotlivé hodnoty modelu. To se dělá pomocí funkce predict() R> predict.lm(lm.vs,interval="confidence") fit lwr upr 1 27.14894 23.55272 30.74515 2 33.58936 30.95924 36.21948 3 27.14894 23.55272 30.74515 [...] 20 60.63915 56.83169 64.44661 R> pro interval spolehlivosti těchto hodnot a R> predict.lm(lm.vs,interval="prediction") fit lwr upr 1 27.14894 18.00498 36.29289 2 33.58936 24.78046 42.39826 3 27.14894 18.00498 36.29289 [...] 20 60.63915 51.41007 69.86823 R> pro interval spolehlivosti předpovězených hodnot. Tyto hodnoty si můžeme vytisknout do grafu (obr. 10.7): R> R> R> R> R> R> R> R> R>

# tisk závislosti délky na šířce: plot(listy$šířka, listy$délka, main="Graf závislosti délky listu na jeho šířce") # tisk funkce ’předpovědi’ jako linie: nejdříve interval spolehlivosti matlines(listy$šířka,predict.lm(lm.vs,interval="c"),lty=c(1,2,2),col=’blue’)

# a dále interval předpovědi matlines(listy$šířka,predict.lm(lm.vs,interval="p"), lty=c(1,9,9),col=c(’black’,’r

Více o analýze reziduí (Jack-Knife rezidua atd.) viz. (Faraway 02). 58

60

Graf závislosti délky listu na jeho šířce

40 30

listy$délka

50

délka = −2.752 + 0.676*šířka Interval spolehlivosti Interval predikce

15

20

25

30

35

40

listy$šířka

Obrázek 10.7: Interval spolehlivosti pro lineární model a interval spolehlivosti predikce modelu.

Postup regresní analýzy při mnohonásobné regresi 1. Nejdříve ověříme, nedochází-li k multikolinearitě. K tomu použijeme VIF (viz 10.2.4) transpozici korelační matice. Kritická hodnota je 10 na ose matice. R> # použijeme pouze sloupečky d13, h, koruna a V..1000.m3 R> solve(cor(stromy[,c(1,2,3,5)])) d13 h koruna V..1000.m3. d13 27.0933405 -2.675107 0.8315288 -24.829190 h -2.6751070 12.023025 -5.5362276 -3.577005 koruna 0.8315288 -5.536228 8.9740653 -3.863416 V..1000.m3. -24.8291902 -3.577005 -3.8634163 32.263736 R> A skutečně, nejméně koreluje s ostatními grafy koruna. 2. V první řadě bychom měli rozhodnout, která data mají být do modelu zahrnuta. . . (ale nevím jak)

10.2.8

Postup regresní analýzy (nelineární model)

Parametry spolu nejsou v lineárním postavení. Mějme data rust1.txt, obsahující sloupec s věkem a sloupec s výškou porostu: 59

R> rust <- read.table("rust1.txt",header=T,sep="\t") R> rust Věk výška 1 8 4.1 2 11 6.3 3 16 9.5 [...] R> Zvolíme si tvar růstové funkce podle Michajlova k

y = a · et

(10.14)

kde y je výška stromu, a je asymptota, ke které se výška stromu blíží, k je koeficient křivosti křivky. Na nelineární modely je v Rku funkce nlm() (nelineární modelování) a nls() (nelineární nejmenší čtverce). Použijeme nls(), protože mi přijde snazší. Funkce nls() potřebuje ke své spokojenosti jednak vzorec, podle kterého bude hledat řešení (ze vzorce 10.14) a*exp(k/rust$Věk), dále počáteční hodnotu parametrů a a k. Ta je předána funkci nls() parametrem start a funkcí list(). Jak nastavit počáteční hodnoty? Víme, že parametr a je asymptota – teoreticky největší možná výška, jaké může strom dosáhnout. Parametr k by měl být záporný, ale to „ jako“ nevíme a nastavíme jej na „nějakou“ hodnotu. R> nls(rust$výška ~ a*exp(k/rust$Věk),start=list(a=60,k=1),data=rust) Nonlinear regression model model: rust$výška ~ a * exp(k/rust$Věk) data: rust a k 46.0972 -37.1337 residual sum-of-squares: 257.7884 R> −37.1337

vyšlo nám tedy, že y = 46.0972 · e x . Nyní můžeme otestovat jednotlivé parametry a celou regresi tak, jak jsme to již několikrát udělali. Uložíme výsledek regrese do proměnné michajlov a dále budeme pracovat s ní. R> michajlov <- nls(rust$výška ~ a*exp(k/rust$Věk),start=list(a=60,k=1),data=rust) R> profile(michajlov) $a tau par.vals.a par.vals.k 1 -3.3805725 39.52522 -27.03273 2 -2.6915422 40.73843 -28.86380 3 -2.0014528 42.01336 -30.80996 4 -1.3103458 43.35514 -32.87633 5 -0.6201533 44.76539 -35.06197 6 0.0000000 46.09720 -37.13370 7 0.6254316 47.50745 -39.33076 8 1.2500278 48.98816 -41.63654 60

9 10 11 12

1.8737982 2.4967326 3.1187901 3.7399263

50.54493 52.18384 53.91152 55.73536

-44.05487 -46.59000 -49.24619 -52.02810

$k 1 2 3 4 5 6 7 8 9 10 11 12

tau par.vals.a par.vals.k -3.7260192 54.96943 -53.10715 -3.0972670 53.25111 -50.11365 -2.4685749 51.63231 -47.24772 -1.8400030 50.10522 -44.50447 -1.2116353 48.66296 -41.87949 -0.5824720 47.29703 -39.36449 0.0000000 46.09720 -37.13370 0.6872286 44.75615 -34.61871 1.3775041 43.48446 -32.21552 2.0684685 42.28162 -29.92929 2.7600913 41.14284 -27.75650 3.4524121 40.06361 -25.69349

attr(,"original.fit") Nonlinear regression model model: rust$výška ~ a * exp(k/rust$Věk) data: rust a k 46.0972 -37.1337 residual sum-of-squares: 257.7884 attr(,"summary") Formula: rust$výška ~ a * exp(k/rust$Věk) Parameters: Estimate Std. Error t value Pr(>|t|) a 46.097200 2.028100 22.72926 < 2.22e-16 *** k -37.133703 3.397076 -10.93108 5.5304e-12 *** --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.931373 on 30 degrees of freedom Correlation of Parameter Estimates: a k -0.8970736 attr(,"class") [1] "profile.nls" "profile" R> 61

Opět vidíme, že parametry koeficientů a a k vyšly 46.097 a -37.133. Suma čtverců vyšla 257.788. Chyba reziduí při 30 stupních volnosti je 2.932. Interval spolehlivosti parametrů získáme pomocí funkce confint() z balíku MASS: R> library(MASS) R> confint(michajlov) Waiting for profiling to be done... 2.5% 97.5% a 41.93627 50.97988 k -45.37536 -30.01519 R> Křivku můžeme do grafu vložit např. následujícím způsobem (obr. 10.8): R> yfit <- 46.0972 * exp(-37.133/rust$Věk) # uložení hodnot výšky R> plot(rust$výška, rust$Věk) # graf růstu a výšky R> lines(rust$Věk, yfit) # linie podle modelu

35

Graf závisloti výšky na věku

20 5

10

15

rust$výška

25

30

y = 46.097 exp(−37.133/x))

20

40

60

80

100

120

rust$Věk

Obrázek 10.8: Graf Michajlovovy růstové funkce Nyní si zkusíme to samé s Korfovou funkcí: k

y = a · e (1−n)·tn−1

(10.15)

R> korf <- nls(rust$výška ~ a*exp(k/((1-n)*rust$Věk^(n-1))),start=list(a=60,k=1,n=1),data= 62

Error in numericDeriv(form[[3]], names(ind), env) : Missing value or an Infinity produced when evaluating the model R> (Tak přes tohle jsem se nedostal . . .)

63

64

11 Anova Analýza rozptylu. Pokud pracujeme s třemi a více výběry zároveň, tak nelze použít opakované t-testy. Např.: H0: µ1 = µ2 = µ3 = . . . = µn H1: Alespoň mezi dvěma středními hodnotami existuje statistický rozdíl Zvyšuje se chyba I. druhu – pro k výběrů platí αβ = 1 − (1 − α)k

(11.1)

Např. pro 7 výběrů a α = 0.05 platí αβ = 0.0302.

Příklad 1. Zkoumáme, zda-li vliv hnojení semenáčků prokazatelně zvýší jejich růst. H0: hnojení nemá vliv (střední hodnoty jsou stejné), H1: hnojení má vliv (alespoň mezi dvěma je statistický rozdíl). Data jsou v souboru semenacky1.txt. Data nemohou být uložena ve formě „matice“, ale spíše ve „sloupečcích“:

R> > semenacky <- read.table("semenacky1.txt",header=T,sep="\t") davka vyska 1 davka1 16 2 davka1 18 3 davka1 17 [...] Data představují čtyři různé lokality L1-4 a způsoby hnojení A-C. Posuzuje se jednak variabilita mezi výběry – což jej rozdíl mezi µ1 , µ2 , . . . , µn a jednak variabilita uvnitř výběru tedy variabilita pro každý výběr zvlášť. Pokud je variabilita mezi výběry větší, než uvnitř výběru, map jsou jednotlivé výběry odlišitelné mezi sebou . 65

Anova

je

parametrická jednofaktorová zkoumá působnost jednoho faktoru s pevnými efekty – přesně se stanová dávka hnojiva. To lze jenom u laboratorních měření s náhodnými efekty – pouze se měří, často situace v terénu vícefaktorová zkoumá působnost dvou a více faktorů. Je důležité rozmyslet si kolik faktorů je důležitých. Pokud možno nezkoumat více, než 3 faktory najednou s pevnými efekty s náhodnými efekty se smýšenými efekty neparametrická . . .

11.0.9

Podmínky použití

1. Vždy musí být porovnávány nezávislé výběry 2. Výběry pochází se ZS s normálním rozdělením 3. Všechny výběry pochází ze souboru se shodnými rozptyly Anova je docela odolná proti porušení třetí podmínky. Nezávislost mezi výběry Lze ji otestovat buď graficky (plot(proměnná1, proměnná2), pak-li že je tam trend, jak např. u obrázku 10.4, nebo pomocí lineárního modelu summary(lm(proměnná1 proměnná2)). Shoda rozptylů Testuje se homoskedasticita, chochranův test, barttletův test. . . > bartlett.test(semenacky$vyska, semenacky$davka) Bartlett test for homogeneity of variances data: vyska and davka Bartlett’s K-squared = 2.0454, df = 3, p-value = 0.563 R> Vidíme, že nezamítneme nulovou hypotézu o rovnosti variancí v jednotlivých skupinách hnojiv. 66

11.1

Model jednofaktorové anovy yij = µ + αi + ǫij

(11.2)

kde yij µ αi ǫij

měřená hodnota průměrná hodnota změna měřené hodnoty způsobená faktorem chyba experimentu rozdíly mezi zásahy si můžeme nechat zobrazit pomocí boxpotu (obr. 11.1):

R> R> R> R> R>

# zpřístupníme si názvy sloupců v proměnné ’semenacky’ attach(semenacky) # nyní můžeme používat přímo názvy sloupců boxplot(split(vyska,davka))

10

12

14

16

18

20

Rozdíly mezi zásahy

davka1

davka2

davka3

davka4

Obrázek 11.1: Znázornění rozdílů mezi zásahy pomocí boxplotu Stejného grafu bychom dosáhli pomocí funkce plot(): R> plot(vyska ~ davka, data=semenacky) 67

Průměr jednotlivých skupin dostaneme např. pomocí funkce tapply(), která aplikuje operaci na každou ze skupin: R> tapply(vyska, davka, mean) davka1 davka2 davka3 davka4 17.1 12.6 14.3 14.5 Stejným způsobem získáme směrodatnou odchylku. Řekli jsme si, že nás zajímá, je-li variabilita uvnitř výběru větší, než variabilita mezi výběry (str. 65). R> tapply(vyska,davka,sd) davka1 davka2 davka3 davka4 1.791957 1.264911 2.002776 1.957890 R> bartlett.test(vyska, davka) Bartlett test for homogeneity of variances data: vyska and davka Bartlett’s K-squared = 2.0454, df = 3, p-value = 0.563 R> Nezamítneme nulovou hypotézu o rovnosti variancí. Vlastní analýzu varianci provedeme následovně: R> analyza<-aov(vyska ~ davka, data=semenacky) R> summary(analyza) Df Sum Sq Mean Sq F value Pr(>F) davka 3 103.475 34.492 10.902 3.073e-05 *** Residuals 36 113.900 3.164 --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 R> P-hodnota Pr je menší, než α (0.05). Můžeme tedy říci, že rozdíly v průměrech jsou statisticky významné a že zamítáme nulovou hypotézu o rovnosti středních hodnot.

11.2

Model vícefaktorové anovy bez opakování

Zkoumáme vliv více faktorů (lokalita, odrůda) a data se nám neopakují (soubor psenice.txt). R> psenice <- read.table("psenice.txt",header=T,sep="\t") R> attach(psenice) R> psenice lokalita odruda vyska 1 Lokalita1 A 27 2 Lokalita1 B 27 3 Lokalita1 C 30 4 Lokalita1 D 19 [....] 68

A vlastní analýza R> analyza<-aov(vyska ~ odruda + lokalita, data=psenice) R> summary(analyza) Df Sum Sq Mean Sq F value Pr(>F) odruda 11 440.72 40.07 12.9592 8.227e-07 *** lokalita 4 32.67 8.17 2.6415 0.06411 . Residuals 20 61.83 3.09 --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 R> Můžeme tedy říci, že lokalita nemá na hodnotu výšky vliv, ale zároveň můžeme zamítnout nulovou hypotézu, že odrůda nemá na výšku vliv.

11.3

Model vícefaktorové anovy s opakováním

Zkoumáme vliv více faktorů (lokalita, způsob hnojení). Data se opakují – máme několik měření pro každou kombinaci faktorů (semenacky2.txt): R> semenacky <- read.table("poznamky/data/semenacky2.txt",header=T,sep="\t") R> attach(semenacky) R> semenacky lokalita hnojeni vyska 1 L1 A 23 2 L1 A 15 3 L1 A 26 4 [...] R> R> analyza<-aov(vyska ~ hnojeni + lokalita, data=semenacky) R> summary(analyza) Df Sum Sq Mean Sq F value Pr(>F) hnojeni 2 344.93 172.47 9.4149 0.0003107 *** lokalita 3 46.05 15.35 0.8379 0.4789924 Residuals 54 989.20 18.32 --Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 R> Můžeme říci, že rozdíly mezi středními hodnotami různých způsobů hnojení jsou statisticky významné.

69

70

Literatura Bonk, R. 2002: Zabudnite na GNUPlot, prichádza ”R”, http://www.root.cz/clanky/zabudnite-na-gnuplot-prichadza-r/, 2005 Čepický, J. 2004: Tvorba grafů pomocí programu ”R”, http://www.root.cz/serialy/tvorba-grafu-pomoci-programu-r/, 2005 Faraway, J. 2002: Practical Regression and Anova using R, http://www.stat.lsa.umich.edu/~faraway/book/, 2005 Šmilauer, P. 2004: Úvod k S (program R 1.7), www.fle.czu.cz/~jachym/skripta/r/r-intro.pdf 2005 Venables W. N., Smith D. M. 2004: An Introduction to R, http://cran.r-project.org/manuals.html, 2005

71