NMFM301, ZS 2022/23

Cvičenie 11 (praktické cvičenie 5)

(metóda analýzy rozptylu, ANOVA)

Zdrojový soubor pro Knit: Rmd soubor

Na úvod si vyskúšajte v programe RStudio otvoriť podkladový Rmd soubor súbor a skompilovať ho pomocou tlačítka Knit (pre sprvne fungovanie kódovania diakritiky je nutné subor otvoriť v kódovaní UTF8 – pomocou ponuky Reopen with encoding v menu programu).

Užitočné materiály pre prácu so štatistickým softwarom R

  • Bína, V., Komárek, A. a Komárková, L.: Jak na jazyk R. (PDF súbor)
  • Komárek, A.: Základy práce s R. (PDF súbor)
  • Kulich, M.: Velmi stručný úvod do R. (PDF súbor)
  • De Vries, A. a Meys, J.: R for Dummies. (ISBN-13: 978-1119055808)


Pre samotnú štatistickú analýzu využijeme data a niektoré dodatočné funkcie (príkazy), ktoré do programu R načítame pomocou nasledujúceho kódu:

load(url("http://msekce.karlin.mff.cuni.cz/~maciak/NMFM301/cv5.RData"))
load(url("http://msekce.karlin.mff.cuni.cz/~maciak/NMFM301/functions.RData"))

K dispozíci máme tri datové soubory zam, deti a kraje. Data zam pocházejí ze sociálního průzkumu provedeného v ČR v roce 2004. Data deti udávají počty dětí narozených v Československu v jednotlivých měsících roku 1970 (zdroj: Anděl, Matematická statistika, SNTL 1977). Datový súbor kraje pochádza z veřejných databází Českého statistického úřadu.

Zaujímať nás bude niekoľko konkrétnych otázok formulovaných v zmysle štatistického testu. Doležitý je správny pravdepodobnostný model a formulácia otázky pomocou vhodnej nulovej a alternatívnej hypotézy. Voľbe pravdepodobnostného modelu (konkrétneho štatistického testu) samozrejme predchádza exploratívna analýza uskutočnená pomocou vhodných popisných charakteristík a obrázkov.

Data si prohlédněte poklepáním na příslušný řádek v okénku “Workspace” nebo pomoci příkazů summary(), head(), str(), a pod. Pomocou niekoľkých obrázkov a grafov preskúmajte štruktúru v datach – vždy s ohľadom na kladenú teoretickú otázku.

I. \(\boldsymbol{\chi^2}\) testy nezávislosti a dobré zhody

V tejto časti budeme používať datovy súbor deti, který udáva počty dětí narozených v Československu v jednotlivých měsících roku 1970. Pomocou grafických nástrojov bežne dostupných v programe R sa na data podívajte a pokuste se data zobrazit vždy tak, aby výsledný graf poskytoval určité možnosti jak nahlédnout do struktúry, které se týka daný statistický test.

Rodí se děti rovnoměrně během roku?

Datový súbor obsahuje len 12 pozorovaní, ktoré predstavujú súhrnné počty narodených detí pre každý mesiac v roku 1970.

barplot(deti$poc.deti,names.arg = c("Leden", " Únor", "Březen", "Duben", "Květen", "Červen", "Červenec", "Srpen", "Žárí", "Říjen", "Listopad", "Prosinec"), col = "lightblue", cex.names = 0.8, ylim = c(0,25000))
abline(mean(deti$poc.deti), 0, col = "red", lwd = 2, lty =2)
text(13.5, 21800, "Celkový priemer", col = "red")

Ověříme, jestli se děti rodí rovnoměrně během roku. Pokud by tomu tak bylo, pravděpodonosti jednotlivých měsíců by byly přímo úměrné počtu dní v měsíci. Zapíšeme tedy počty dní v měsících do vektoru a ověříme, že součet dá 365 (rok 1970 nebyl přestupný). Vieme teda rozdelenie pravdepodobnosti, ktoré zodpovedá nulovej hypotéze. Z matematického/štatistického hľadiska preto využijeme tzv. \(\chi^2\) test dobrej zhody.

dni <- c(31,28,31,30,31,30,31,31,30,31,30,31)
sum(dni)
## [1] 365

Nyní napočítáme pravděpodobnosti měsíců za nulové hypotézy.

psti <- dni/sum(dni)
psti
##  [1] 0.08493151 0.07671233 0.08493151 0.08219178 0.08493151 0.08219178
##  [7] 0.08493151 0.08493151 0.08219178 0.08493151 0.08219178 0.08493151

Popripade tie ist pravdepodobnosti s presnosťou na štyri desatinné miesta:

round(psti, 4)
##  [1] 0.0849 0.0767 0.0849 0.0822 0.0849 0.0822 0.0849 0.0849 0.0822 0.0849
## [11] 0.0822 0.0849

Odhadnuté pravdepodobnosti z datového súboru (opäť na štyri desatinne miesta):

round(deti$poc.deti / sum(deti$poc.deti), 4)
##  [1] 0.0838 0.0790 0.0902 0.0902 0.0915 0.0865 0.0845 0.0805 0.0829 0.0793
## [11] 0.0741 0.0775



Budeme testovať nulovô hypotézu, že vektor pravděpodobností multinomického rozdělení, které vygenerovalo počty narozených dětí, je roven vektoru pravdepodobnosti, ktorý máme uložený v R objekte psti. Formálne (matematicky) zapísané, zaujíma nas test nulovej hypotézy \[ H_0:~\boldsymbol{p} = \boldsymbol{p}_0, \] kde \(\boldsymbol{p} = (p_1, \dots, p_{12})^\top \in [0,1]^{12}\) je vektor pravdepodobnosti, že sa náhodne vybrané dieťa narodí v niektorom z 12tich mesiacov. Samozrejme platí, že \(\sum_{i = 1}^{12} p_i = 1\), inak by sa totiž nejednalo o pravdepodobnostné rozdelenie. Pripomeňte si teoretické základy \(\chi^2\) testu dobrej zhody.

V programe R oužijeme funkci chisq.test. Její první argument je vektor pozorovaných četností a druhý, argument p, obsahuje hypotetické pravděpodobnosti - vektor pravdepodobnosti v multinomickom rozdelení za platnosti nulovej hypotézy.

chisq.test(deti$poc.deti,p=psti)
## 
##  Chi-squared test for given probabilities
## 
## data:  deti$poc.deti
## X-squared = 1004.1, df = 11, p-value < 2.2e-16

Testová statistika je velmi velká (referenční rozdělení \(\chi^2_{11}\) má střední hodnotu 11), p-hodnota je takřka nulová. Hypotézu tedy s velkou rezervou zamítáme.

Podívejme se blíže na způsob výpočtu testové statistiky. Testová štatistika \(\chi^2\) testu dobrej zhody, je definovaná predpisom \[ \chi^2=\sum_{k=1}^K\frac{(X_k-np^0_k)^2}{np^0_k}. \]

Pre lepšie pochopenie pomůže nám funkce kuchej.chisq (není součástí R).

kuchej.chisq(deti$poc.deti,psti)
##    skupina cetnost        pst      ocek        rozdil   statistika
## 1        1   21182 0.08493151  21465.59 -2.835890e+02    3.7465892
## 2        2   19960 0.07671233  19388.27  5.717260e+02   16.8591929
## 3        3   22787 0.08493151  21465.59  1.321411e+03   81.3453998
## 4        4   22805 0.08219178  20773.15  2.031849e+03  198.7378661
## 5        5   23120 0.08493151  21465.59  1.654411e+03  127.5099237
## 6        6   21859 0.08219178  20773.15  1.085849e+03   56.7592636
## 7        7   21367 0.08493151  21465.59 -9.858904e+01    0.4528084
## 8        8   20357 0.08493151  21465.59 -1.108589e+03   57.2530136
## 9        9   20946 0.08219178  20773.15  1.728493e+02    1.4382453
## 10      10   20037 0.08493151  21465.59 -1.428589e+03   95.0762005
## 11      11   18728 0.08219178  20773.15 -2.045151e+03  201.3484323
## 12      12   19592 0.08493151  21465.59 -1.873589e+03  163.5331734
## 13  Celkem  252740 1.00000000 252740.00  1.818989e-11 1004.0601088

V prvním sloupečku je číslo \(k\) nebo název skupiny (kategorie). Ve druhém je pozorovaná četnost \(X_k\). Dále následuje hypotetická pravděpodobnost \(p^0_k\) a očekávaná četnost (kdyby hypotéza platila) \(np^0_k\). V posledních dvou sloupečcích je rozdíl \(X_k-np^0_k\) mezi pozorovanou a očekávanou četností a příspěvek \((X_k-np^0_k)^2/(np^0_k)\) do testové statistiky. Vidíme, které kategorie (měsíce) nejvíce přispěly k hodnotě 1004.06, která vedla k zamítnutí hypotézy. Byly to duben a květen, kdy se narodilo mnohem více dětí, než kolik by mělo, kdyby se děti rodily rovnoměrně, a dále listopad a prosinec, kdy se naopak narodilo dětí méně.

Je poměr nezamestnaných stejný u mužů a u žen?

Funkciu chisq.test je možné použíť aj k testovaniu nulovej hypotézy o nezávislosti dvoch premenných – tzv. \(\chi^2\) test nezávislosti. Ak budeme pre jednoduchosť uvažovať pouze dve binárne premenné – napr. náhodnú veličinu \(X_1\) (ktorá nabýva hodnot 1 a 2 s pravdepodobnostou \(p_11\) a \(p_12\)) a druhú náhodnú veličinu \(X_2\) (ktorá nabýva hodnot 1 a 2 s pravdepodobnostou \(p_21\) a \(p_22\)), tak nulovú hypotézu nezávislosti lze formálne (matematicky) zapísať takto:

\[ H_0: p_{i j} = p_{i \bullet} p_{\bullet j} \qquad \textrm{pre všetky $i, j \in \{1,2\}$.} \] Samozrejme platí, že \(p_{11} + p_{12} = 1\) a \(p_{21} + p_{22} = 1\) (v opačnom prípade by sa nejednalo o rozdelenie pravdepodobnosti) a tiež

\[ p_{i \bullet} = p_{i 1} + p_{i 2} \qquad \textrm{a} p_{\bullet j} = p_{1 j} + p_{2 j} \qquad \textrm{opäť pre všetky $i,j \in \{1,2\}$.} \] Pozorovaná kontingenčná tabuľka (pre nezamestnanosť mužov a žien) je

tbl <- table(zam$unempl, zam$gndr)
tbl
##    
##     Male Female
##   0 1238   1407
##   1   57     97

Pomocou funkcie barplot() možeme na data nahliadnúť pomocou obrázku

barplot(tbl, cex.names = 0.8, horiz = T, col = c("darkblue","red"), xlim = c(0,1600))

Keďže nás ale zaujíma pomer, lepšie bude pozrieť sa na data prostredníctvom relatívnych čísel:

(tbl_rep<- 100 * cbind(tbl[,1]/sum(tbl[,1]), tbl[,2]/sum(tbl[,2])))
##        [,1]      [,2]
## 0 95.598456 93.550532
## 1  4.401544  6.449468
barplot(tbl_rep, cex.names = 0.8, horiz = T, col = c("darkblue","red"), xlim = c(0,100))

Spočítajte relatívne četnosti a pomocou funkcie chisq.test aplikujte test nezávislosti medzi pohlavím a mierou nezamestnanosti.

chisq.test(tbl)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tbl
## X-squared = 5.2261, df = 1, p-value = 0.02225

Aká je výsledná interpretácia testu?

Je úmrtnost ve věku 15-19 let stejná ve všech krajích?

Nyní použijeme data kraje.

kraje
##                             kraj    obyv naroz posl zemr.15.19 obyv.15.19
## 1       Hlavn\xed m\xecsto Praha 1246780 14233   24          9      47373
## 2      St\xf8edo\xe8esk\xfd kraj 1291816 14483   25         27      61552
## 3           Jiho\xe8esk\xfd kraj  636611  6672   12         22      31887
## 4            Plze\xf2sk\xfd kraj  572687  5785   11          8      26926
## 5            Karlovarsk\xfd kraj  301726  2831    5         13      15177
## 6             \xdasteck\xfd kraj  826764  8246   14         15      42749
## 7              Libereck\xfd kraj  438594  4609    8          7      22035
## 8  Kr\xe1lov\xe9hradeck\xfd kraj  552946  5489   11         12      27779
## 9             Pardubick\xfd kraj  516440  5405   10          8      26552
## 10              Kraj Vyso\xe8ina  511207  5166   11         12      27476
## 11          Jihomoravsk\xfd kraj 1168650 12385   23         17      55909
## 12             Olomouck\xfd kraj  637609  6319   12          9      31777
## 13            Zl\xednsk\xfd kraj  587693  5512   12         10      29521
## 14       Moravskoslezsk\xfd kraj 1226602 11820   22         19      63552

Budeme testovat, zdali všechny kraje měly stejné riziko úmrtí mladých lidí ve věku 15-19 let. Ve sloupečku zemr.15.19 vidíme počty zemřelých v této věkové kategorii v jednotlivých krajích (v roce 2012). Nesmíme však zapomenout, že kraje jsou různě velké a úmrtnost musíme posuzovat úměrně k počtu obyvatel v dané věkové kategorii. Tyto počty jsou uvedeny ve sloupečku obyv.15.19.

Napočítáme nejprve pravděpodobnosti, které by odpovídaly stejné úmrtnosti ve všech krajích.

psti <- kraje$obyv.15.19/sum(kraje$obyv.15.19)
psti
##  [1] 0.09283999 0.12062752 0.06249106 0.05276866 0.02974337 0.08377804
##  [7] 0.04318344 0.05444034 0.05203571 0.05384653 0.10956856 0.06227548
## [13] 0.05785425 0.12454705

Nyní provedeme test hypotézy, že rozdělení počtu zemřelých do krajů se řídí právě těmito pravděpodobnostmi a vykucháme výpočet testové statistiky.

chisq.test(kraje$zemr.15.19,p=psti)
## 
##  Chi-squared test for given probabilities
## 
## data:  kraje$zemr.15.19
## X-squared = 27.376, df = 13, p-value = 0.01105
kuchej.chisq(kraje$zemr.15.19,psti,skup=kraje$kraj)
##                          skupina cetnost        pst       ocek        rozdil
## 1       Hlavn\xed m\xecsto Praha       9 0.09283999  17.453919 -8.453919e+00
## 2      St\xf8edo\xe8esk\xfd kraj      27 0.12062752  22.677973  4.322027e+00
## 3           Jiho\xe8esk\xfd kraj      22 0.06249106  11.748319  1.025168e+01
## 4            Plze\xf2sk\xfd kraj       8 0.05276866   9.920508 -1.920508e+00
## 5            Karlovarsk\xfd kraj      13 0.02974337   5.591753  7.408247e+00
## 6             \xdasteck\xfd kraj      15 0.08377804  15.750271 -7.502709e-01
## 7              Libereck\xfd kraj       7 0.04318344   8.118487 -1.118487e+00
## 8  Kr\xe1lov\xe9hradeck\xfd kraj      12 0.05444034  10.234784  1.765216e+00
## 9             Pardubick\xfd kraj       8 0.05203571   9.782713 -1.782713e+00
## 10              Kraj Vyso\xe8ina      12 0.05384653  10.123148  1.876852e+00
## 11          Jihomoravsk\xfd kraj      17 0.10956856  20.598889 -3.598889e+00
## 12             Olomouck\xfd kraj       9 0.06227548  11.707791 -2.707791e+00
## 13            Zl\xednsk\xfd kraj      10 0.05785425  10.876599 -8.765994e-01
## 14       Moravskoslezsk\xfd kraj      19 0.12454705  23.414845 -4.414845e+00
## 15                        Celkem     188 1.00000000 188.000000  8.881784e-16
##     statistika
## 1   4.09471059
## 2   0.82370304
## 3   8.94570219
## 4   0.37179053
## 5   9.81483197
## 6   0.03573948
## 7   0.15409449
## 8   0.30445078
## 9   0.32486544
## 10  0.34797223
## 11  0.62877181
## 12  0.62626095
## 13  0.07064952
## 14  0.83241457
## 15 27.37595759

Testová statistika vyšla 27.38. Porovnává se s kvantilem rozdělení \(\chi^2_{13}\) (14 krajů - 1), které má střední hodnotu 13. Její hodnota stačí na zamítnutí nulové hypotézy (p-hodnota 0.011). Tedy zamítáme hypotézu, že všechny kraje měly stejné riziko úmrtí mladých lidí ve věku 15-19 let.

Které kraje se nejvíce liší? Prakticky celou pozorovanou hodnotu testové statistiky vytvořily pouze tři kraje: Praha, kde byla úmrtnost nižší než jinde, a Jihočeský a Karlovarský kraj, kde byla úmrtnost naopak velmi vysoká. Máte-li v těchto dvou krajích mladší kamarády, radši jim poraďte, aby se raději přestěhovali buď do Prahy anebo alespoň do kraje Ústeckého nebo Moravskoslezského.

Samostatne

Otázka, zda je úmrtnosť lidí ve věku 15-19 let stejná ve všech krajích je vpodstatě ekvivalentní s otázkou, zda lze považovat úmrtnost a kraj za nezávisle veličiny. Formálne to vypadá na použití statistického testu pro testování nezávislosti - \(\chi^2\) testu nezávislosti. Je možné tento test aplikovat na data z datového souboru zam? Pokud ano, udělejte to, pokud ne, vysvětlete proč.



II. Inference pre lineárne kombinácie pravdepodobnosti

Budeme pracovat s daty zam a zkoumat velikosti městské vs. venkovské populace v ČR v roce 2004. Velikost bydliště respondenta je uvedena ve veličině zam$domicil. Funkce table nám spočítá, kolik respondentů se nachází v existujících kategoriích této veličiny.

tbl <- table(zam$domicil)
tbl
## 
##                       A big city Suburbs or outskirts of big city 
##                              519                              138 
##               Town or small city                  Country village 
##                             1360                              776 
##      Farm or home in countryside 
##                                6

Celkový počet respondentů je

sum(tbl)
## [1] 2799

Vizuálne môžeme data reprezentovať napr. pomocou koláčového grafu (tzv. piechart):

library("colorspace") 
pie(table(zam$domicil), col = diverging_hcl(5, c = 100, l = c(50, 90)))

Vhodným modelem pro tuto tabulku je multinomické rozdělení s \(K=5\) kategoriemi, celkovým počtem pozorování \(n=2799\) a neznámým vektorem pravděpodobností \(p\). Ten můžeme odhadnout vektorem relativních četností \(\widehat{p}_n\):

tbl/sum(tbl)
## 
##                       A big city Suburbs or outskirts of big city 
##                      0.185423365                      0.049303323 
##               Town or small city                  Country village 
##                      0.485887817                      0.277241872 
##      Farm or home in countryside 
##                      0.002143623

Mohli bychom využít i funkce prop.table, která počítá totéž:

prop.table(tbl)
## 
##                       A big city Suburbs or outskirts of big city 
##                      0.185423365                      0.049303323 
##               Town or small city                  Country village 
##                      0.485887817                      0.277241872 
##      Farm or home in countryside 
##                      0.002143623

Otestujme nejprve, zdali je pravda, že procento lidí žijících na venkově je o 10 větší než procento lidí žijících ve velkoměstech. Venkovská populace odpovídá posledním dvěma kategoriím tabulky. Velkoměstská populace je v první kategorii tabulky. Jde tedy o to, zdali součet posledních dvou pravděpodobností, tj. \(p_4+p_5\), je o 0.1 větší než pravděpodobnost \(p_1\).

Stanovíme vektor \(c=(-1,0,0,1,1)^T\) a budeme testovat hypotézu \(H_0: c^T p=0.1\). Použijeme funkci multitest, která není součástí R, ale implementuje vzorec (7.4) z větníku (str. 68-69). Prvním argumentem této funkce je realizace vektoru s multinomickým rozdělením, druhý argument specifikuje vektor \(c\) pro lineární kombinaci parametrů. Třetí argument uvádí hodnotu lineární kombinace \(c^T p\) za nulové hypotézy.

multitest(tbl,comb=c(-1,0,0,1,1),c0=0.1)
## Warning in c(-1, 1) * d: Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.
## 
## Testovani linearnich kombinací pravdepodobnosti
## 
## Multinomicky vektor: 519 138 1360 776 6 
## 
## Odhadnute pravdepodobnosti: 0.18500 0.04930 0.48600 0.27700 0.00214 
## H0: - p1 + p4 + p5 = 0.1 
## 
## Odhad - p1 + p4 + p5 = 0.09396 
## 
## Testova statistika = -0.4731, p-hodnota = 0.636
## 
## Priblizny 95%-ni interval spolehlivosti:
## ( 0.0689 , 0.119 )

Odhad požadované lineární kombinace je 0.09396, což je blízko 0.1. Testová statistika (7.4) je -0.4731. Test má p-hodnotu 0.636, takže nulovou hypotézu nelze zamítnout. Interval spolehlivosti říká, že rozdíl mezi procentem lidí žijících na venkově a procentem lidí žijících ve velkoměstech je s pravděpodobností 0.95 někde mezi 6.9 a 11.9.

Podobně můžeme otestovat hypotézu, že počet lidí žijících v menších městech je dvakrát větší než počet lidí žijících ve velkoměstech a jejich bezprostředním okolí.

multitest(tbl,comb=c(1,1,-1/2,0,0),c0=0)
## Warning in c(-1, 1) * d: Recycling array of length 1 in vector-array arithmetic is deprecated.
##   Use c() or as.vector() instead.
## 
## Testovani linearnich kombinací pravdepodobnosti
## 
## Multinomicky vektor: 519 138 1360 776 6 
## 
## Odhadnute pravdepodobnosti: 0.18500 0.04930 0.48600 0.27700 0.00214 
## H0:  p1 + p2 -0.5 p3 = 0 
## 
## Odhad  p1 + p2 -0.5 p3 = -0.008217 
## 
## Testova statistika = -0.7285, p-hodnota = 0.466
## 
## Priblizny 95%-ni interval spolehlivosti:
## ( -0.0303 , 0.0139 )

Sami si promyslete, proč je vektor \(c\) nastaven takto a co znamenají výsledky testu.

Samostatne

  • Pomocou základných popisných charakteristík a pomocou rôznych grafických nástrojov sa podívajte na štruktúru dat a navrhnite vlastnú lineárnu kombináciu pre vektor neznámich pravdepodobnosti v multinomickom rozdelení a test pre nulovosť vami zvolenej lineárnej kombinácie pomocou programu R spočítajte.
  • Interpretujte výsledky testu vrámci kontextu celého problému a dát, ktoré máte k dispozícii.



III. Metóda analýzy rozptylu (ANOVA)

Metóda analýzy rozptylu (z anglického The Analysis of Variance - ANOVA) nám umožňuje ověřit, zda na nějaku (spojitou) hodnotu náhodné veličiny má statisticky významný vliv hodnota jinej náhodnej veličiny (diskrétnej). Táto diskrétní veličina (resp. charakteristika, nebo znak) musí nabývat pouze konečného počtu možných hodnot (nejméně dvou). Slouží teda k rozdělení jedinců do skupin, které pak pomoci ANOVA vzájemně porovnávame.

Je úroveň vzdelania rovnaká v mestách a na venkově?

V tejto úlohe budeme používať datový súbor zam, které pocházejí ze sociálního průzkumu provedeného v ČR v roce 2004. Úroveň vzdelania v tejto úlohe budeme sledovať pomocou premennej eduyrs, ktorá zaznamenáva počet rokov strávených štúdiom u každého subjektu zvlášť. Na jednotlivé populácie definované miestom trvalého bydliska sa možeme pozrieť pomocou boxplotu.

levels(zam$domicil) <- c("Big City", "Suburbs", "Town/City", "Village", "Countryside")
cols <- rainbow(5, start = 0, end = 0.1)
boxplot(zam$eduyrs ~ zam$domicil, col = cols, main = "Dlzka vzdelania v rokoch")

Prvé štyri poplulácie sa zdajú byť rovnaké čo sa týka školskej dochádzky (v rokoch), ale vidiecké lokality vykazujú v priemere o niečo nižiu mieru vzdelanosti medzi svojími obyvateľmi (t.j. v primere menší počet rokov trávených na škole).

Pomocou metódy analýzy rozptylu sa pokúsime zistiť, či je tento rozdiel štatisticky významný, alebo nie. Využijeme k tomu najprv príkaz aov() (pre podrobnú implementáciu príkazu použijte help, t.j. ?aov). štandardne dostupných v programe R.

aov(zam$eduyrs ~ zam$domicil)
## Call:
##    aov(formula = zam$eduyrs ~ zam$domicil)
## 
## Terms:
##                 zam$domicil Residuals
## Sum of Squares       319.82  15219.46
## Deg. of Freedom           4      2794
## 
## Residual standard error: 2.333922
## Estimated effects may be unbalanced

Vyššie uvedený príkaz porovnajte ešte s výstupom nižšie, ktorý získame dvoj-kombináciou iných dvoch príkazov – anova() a lm():

anova(lm(zam$eduyrs ~ zam$domicil))
## Analysis of Variance Table
## 
## Response: zam$eduyrs
##               Df  Sum Sq Mean Sq F value    Pr(>F)    
## zam$domicil    4   319.8  79.955  14.678 7.183e-12 ***
## Residuals   2794 15219.5   5.447                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dokážete v jednotlivých výstupoch identifikovať jednotlivé súčty štvorcov a uvedené čísla interpretovať?

Ako je formulovaná nulová hypotéza, ktorú testujeme? Aký je výsledok testu a aká je interpretácia vzhľadom k pôvodnej otázke? Je možné z výsledkov testu určiť, kde a ako sa jednotlivé (sub-)populácie od seba líšia?

Je rozdíl ve vzdělání mužov a žien?

Odpoveď budeme opäť zisťovať pomocou metódy analýzy rozptylu. Najprv sa ale pozrieme na obe populácie (mužov a ženy) pomocou grafických nástrojov. Využijeme boxplot

boxplot(zam$eduyrs ~ zam$gndr, col = c("lightblue","pink"), main = "Dlzka vzdelania v rokoch")

prípadne histogram

par(mfrow = c(1,2))
hist(zam$eduyrs[zam$gndr == "Male"], col = "lightblue", xlab = "Roky stravene na skole", main = "Males", breaks = 10)
hist(zam$eduyrs[zam$gndr == "Female"], col = "pink", xlab = "Roky stravene na skole", main = "Females", breaks = 10)

Možeme skúsiť porovnať aj empirický odhad hustoty (neparametrický) s odhadom hustoty založenom na predpoklade, že rozdelenia sú oba normálne, a líšía sa iba vrámci parametru strednej hodnoty (homoskedasticita), prípadne aj rozptylu (heteroskedasticita).
K tomu je potrebný trochu detailnejší R kod…

par(mfrow = c(1,2))
hist(zam$eduyrs[zam$gndr == "Male"], col = "lightblue", xlab = "Roky stravene na skole", main = "Males", freq = F, ylim = c(0, 0.35), breaks = 20)
lines(density(zam$eduyrs[zam$gndr == "Male"], adjust = 3), col = "red", lwd = 2)
lines(density(rnorm(10000, mean(zam$eduyrs[zam$gndr == "Male"]), sd(zam$eduyrs))), col = "blue", lwd = 1, lty = 2)
lines(density(rnorm(10000, mean(zam$eduyrs[zam$gndr == "Male"]), sd(zam$eduyrs[zam$gndr == "Male"]))), col = "green", lwd = 1, lty = 2)
legend(15, 0.35, c("Empirical", "Homoscedastic", "Heteroscedastic"), col = c("red", "blue", "green"), lty = c(1,2,2), cex = 0.8)

hist(zam$eduyrs[zam$gndr == "Female"], col = "pink", xlab = "Roky stravene na skole", main = "Females", freq = F, ylim = c(0,0.35), breaks = 20)
lines(density(zam$eduyrs[zam$gndr == "Female"], adjust = 3), col = "red", lwd = 2)
lines(density(rnorm(10000, mean(zam$eduyrs[zam$gndr == "Female"]), sd(zam$eduyrs))), col = "blue", lwd = 1, lty = 2)
lines(density(rnorm(10000, mean(zam$eduyrs[zam$gndr == "Female"]), sd(zam$eduyrs[zam$gndr == "Female"]))), col = "green", lwd = 1, lty = 2)
legend(15, 0.35, c("Empirical", "Homoscedastic", "Heteroscedastic"), col = c("red", "blue", "green"), lty = c(1,2,2), cex = 0.8)

Pomocou vhodného štatistického testu sa pozrieme, či je rozdiel, ktorý prostredníctvom histogramov sledujeme, statistický významný, alebo je zanedbateľný. Test založíme na vyhodnotení rozptylu v jednotlivých skupinách - tzv. metoda analýzy rozptylu - ANOVA.

anova(lm(zam$eduyrs ~ zam$gndr))
## Analysis of Variance Table
## 
## Response: zam$eduyrs
##             Df Sum Sq Mean Sq F value    Pr(>F)    
## zam$gndr     1    120 120.049  21.777 3.207e-06 ***
## Residuals 2797  15419   5.513                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Keďže porovnávame iba dve populácie, môžeme použiť aj klasický dvojvýberový postup.

t.test(zam$eduyrs ~ zam$gndr, paired = F, var.equal = T)
## 
##  Two Sample t-test
## 
## data:  zam$eduyrs by zam$gndr
## t = 4.6665, df = 2797, p-value = 3.207e-06
## alternative hypothesis: true difference in means between group Male and group Female is not equal to 0
## 95 percent confidence interval:
##  0.2408303 0.5898852
## sample estimates:
##   mean in group Male mean in group Female 
##             12.41004             11.99468

Ktoré výsledky sú vo výstupe analýzy rozptylu a t.testu rovnaké a prečo? Aká je interpretácia výsledkov vzhľadom k pôvodnej otázke?

Samostatne

  • Metóda analýzy rozptylu pracuje s konkrétny štatistickým modelom. Ako tento model vyzerá a na akých predpokladoch je založeny?
  • Aké následky má nesplnenie jednotlivých predpokladov modelu ANOVA na výsledky a interpretáciu výsledkov získaných v teste pomocou metody analýzy rozptylu?



Neparametrické verze pro analýzu rozptylu

Neparametrická verzia metódy analýzy rozptylu - porovnávanie niekoľkých skupín, je implementována v programe R pomocou funkcie kruskal.test(). Pripomeňte si, z akým modelom Kruskal Wallisov test pracuje a aké sú predpokladty testu. Ako vyzerá nulová a ako vyzerá alternatívna hypotéza?

kruskal.test(zam$eduyrs ~ zam$domicil)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  zam$eduyrs by zam$domicil
## Kruskal-Wallis chi-squared = 59.495, df = 4, p-value = 3.703e-12



Samostatne

  • Aplikujte neparametrický postup - tzv. Kruskal-Wallisov test na oba príklady uvedené vyššie.
  • Líšia sa výsledky testov založených na metode ANOVA a na neparmetrickom postupe pomocou Kruskal-Wallisovho testu? Aká je interpretácia výsledkov vzhľadom k celkovému kontextu experimentu?