NMFM301, ZS 2015/16
Cvičenie 11 (praktické cvičení 5 – Instruktáž)
Příprava k práci
Přesvědčte se, že pracujete se síťovou verzí R (na disku L:) a váš pracovní direktorář je umístěn na disku H: (nikoli C:).
Načtěte data pro dnešní cvičení a potřebné funkce:
load(url("http://msekce.karlin.mff.cuni.cz/~maciak/NMFM301/cv5.RData"))
load(url("http://msekce.karlin.mff.cuni.cz/~maciak/NMFM301/functions.RData"))Nyní máte k dispozici 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). Data kraje pocházejí z veřejných databází Českého statistického úřadu.
Data si prohlédněte poklepáním na příslušný řádek v okénku “Workspace” nebo pomoci příkazu summary(), head(), str(), a pod.
Odhadování a testování lineárních kombinací pravděpodobností
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
## 6Celkový počet respondentů je
sum(tbl)## [1] 2799Vhodný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.002143623Mohli 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.002143623Otestujme 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)##
## 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)##
## 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.
\(\chi^2\) test dobré shody
Rodí se děti rovnoměrně během roku?
Data deti udávají počty dětí narozených v Československu v jednotlivých měsících roku 1970.
deti## mesic poc.deti
## 1 1 21182
## 2 2 19960
## 3 3 22787
## 4 4 22805
## 5 5 23120
## 6 6 21859
## 7 7 21367
## 8 8 20357
## 9 9 20946
## 10 10 20037
## 11 11 18728
## 12 12 19592Ověří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ý).
dni <- c(31,28,31,30,31,30,31,31,30,31,30,31)
sum(dni)## [1] 365Nyní 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.08493151Otestujeme hypotézu, že vektor pravděpodobností multinomického rozdělení, které vygenerovalo počty narozených dětí, je roven vektoru psti. Použijeme funkci chisq.test. Její první argument je vektor pozorovaných četností. Argument p obsahuje hypotetické pravděpodobnosti.
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-16Testová 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. Připomeneme si její vzoreček: \[
\chi^2=\sum_{k=1}^K\frac{(X_k-np^0_k)^2}{np^0_k}
\] Pomůže nám funkce kuchej.chisq (není součástí R).
kuchej.chisq(deti$poc.deti,psti)## Warning in `[<-.factor`(`*tmp*`, iseq, value = "Celkem"): invalid factor
## level, NA generated## 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 <NA> 252740 1.00000000 252740.00 1.818989e-11 1004.0601088V první 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ě.
Metóda analýzy rozptylu (ANOVA)
Je úroveň vzdelania rovnaká v mestách a na vekově?
Ú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.
boxplot(zam$eduyrs ~ zam$domicil, col = 1:5, 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 sa zdá vykazujú nižiu mieru vzdelanosti u svojich obyvateľov (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 kombináciu príkazov anova() a lm() štandardne dostupných v programe R.
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 ' ' 1Aka nulová hypotéza sa testuje? Aký je výsledok testu a aká je interpretácia vzhľadom k pôvodnej otázke. Možeme z výsledkov testu určiť, kde a ako sa jednotlivé 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")
hist(zam$eduyrs[zam$gndr == "Female"], col = "pink", xlab = "Roky stravene na skole", main = "Females")
a pomocou štatistického testu sa pozrieme, či je rozdiel štatistický významný, alebo nie.
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 ' ' 1Keďž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 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.99468Ktoré 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?
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 63552Budeme 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.12454705Nyní 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.01105kuchej.chisq(kraje$zemr.15.19,psti,skup=kraje$kraj)## Warning in `[<-.factor`(`*tmp*`, iseq, value = "Celkem"): invalid factor
## level, NA generated## skupina cetnost pst ocek
## 1 Hlavn\xed m\xecsto Praha 9 0.09283999 17.453919
## 2 St\xf8edo\xe8esk\xfd kraj 27 0.12062752 22.677973
## 3 Jiho\xe8esk\xfd kraj 22 0.06249106 11.748319
## 4 Plze\xf2sk\xfd kraj 8 0.05276866 9.920508
## 5 Karlovarsk\xfd kraj 13 0.02974337 5.591753
## 6 \xdasteck\xfd kraj 15 0.08377804 15.750271
## 7 Libereck\xfd kraj 7 0.04318344 8.118487
## 8 Kr\xe1lov\xe9hradeck\xfd kraj 12 0.05444034 10.234784
## 9 Pardubick\xfd kraj 8 0.05203571 9.782713
## 10 Kraj Vyso\xe8ina 12 0.05384653 10.123148
## 11 Jihomoravsk\xfd kraj 17 0.10956856 20.598889
## 12 Olomouck\xfd kraj 9 0.06227548 11.707791
## 13 Zl\xednsk\xfd kraj 10 0.05785425 10.876599
## 14 Moravskoslezsk\xfd kraj 19 0.12454705 23.414845
## 15 <NA> 188 1.00000000 188.000000
## rozdil statistika
## 1 -8.453919e+00 4.09471059
## 2 4.322027e+00 0.82370304
## 3 1.025168e+01 8.94570219
## 4 -1.920508e+00 0.37179053
## 5 7.408247e+00 9.81483197
## 6 -7.502709e-01 0.03573948
## 7 -1.118487e+00 0.15409449
## 8 1.765216e+00 0.30445078
## 9 -1.782713e+00 0.32486544
## 10 1.876852e+00 0.34797223
## 11 -3.598889e+00 0.62877181
## 12 -2.707791e+00 0.62626095
## 13 -8.765994e-01 0.07064952
## 14 -4.414845e+00 0.83241457
## 15 8.881784e-16 27.37595759Testová 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.