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://www.karlin.mff.cuni.cz/~pesta/NMFM301/cv5.RData"))
load(url("http://www.karlin.mff.cuni.cz/~pesta/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”.
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.185423 0.049303
## Town or small city Country village
## 0.485888 0.277242
## Farm or home in countryside
## 0.002144Mohli bychom využít i funkce prop.table, která počítá totéž:
prop.table(tbl)##
## A big city Suburbs or outskirts of big city
## 0.185423 0.049303
## Town or small city Country village
## 0.485888 0.277242
## Farm or home in countryside
## 0.002144Otestujme 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<c3><ad> 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<c3><ad> 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.
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.08493 0.07671 0.08493 0.08219 0.08493 0.08219 0.08493 0.08493
## [9] 0.08219 0.08493 0.08219 0.08493Otestujeme 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, 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: invalid factor level, NA generated## skupina cetnost pst ocek rozdil statistika
## 1 1 21182 0.08493 21466 -2.836e+02 3.7466
## 2 2 19960 0.07671 19388 5.717e+02 16.8592
## 3 3 22787 0.08493 21466 1.321e+03 81.3454
## 4 4 22805 0.08219 20773 2.032e+03 198.7379
## 5 5 23120 0.08493 21466 1.654e+03 127.5099
## 6 6 21859 0.08219 20773 1.086e+03 56.7593
## 7 7 21367 0.08493 21466 -9.859e+01 0.4528
## 8 8 20357 0.08493 21466 -1.109e+03 57.2530
## 9 9 20946 0.08219 20773 1.728e+02 1.4382
## 10 10 20037 0.08493 21466 -1.429e+03 95.0762
## 11 11 18728 0.08219 20773 -2.045e+03 201.3484
## 12 12 19592 0.08493 21466 -1.874e+03 163.5332
## 13 <NA> 252740 1.00000 252740 1.819e-11 1004.0601V 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ě.
Nyní použijeme data kraje.
kraje## kraj obyv naroz posl zemr.15.19 obyv.15.19
## 1 Hlavn\355 m\354sto Praha 1246780 14233 24 9 47373
## 2 St\370edo\350esk\375 kraj 1291816 14483 25 27 61552
## 3 Jiho\350esk\375 kraj 636611 6672 12 22 31887
## 4 Plze\362sk\375 kraj 572687 5785 11 8 26926
## 5 Karlovarsk\375 kraj 301726 2831 5 13 15177
## 6 \332steck\375 kraj 826764 8246 14 15 42749
## 7 Libereck\375 kraj 438594 4609 8 7 22035
## 8 Kr\341lov\351hradeck\375 kraj 552946 5489 11 12 27779
## 9 Pardubick\375 kraj 516440 5405 10 8 26552
## 10 Kraj Vyso\350ina 511207 5166 11 12 27476
## 11 Jihomoravsk\375 kraj 1168650 12385 23 17 55909
## 12 Olomouck\375 kraj 637609 6319 12 9 31777
## 13 Zl\355nsk\375 kraj 587693 5512 12 10 29521
## 14 Moravskoslezsk\375 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.09284 0.12063 0.06249 0.05277 0.02974 0.08378 0.04318 0.05444
## [9] 0.05204 0.05385 0.10957 0.06228 0.05785 0.12455Nyní 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.38, df = 13, p-value = 0.01105kuchej.chisq(kraje$zemr.15.19,psti,skup=kraje$kraj)## Warning: invalid factor level, NA generated## skupina cetnost pst ocek rozdil
## 1 Hlavn\355 m\354sto Praha 9 0.09284 17.454 -8.454e+00
## 2 St\370edo\350esk\375 kraj 27 0.12063 22.678 4.322e+00
## 3 Jiho\350esk\375 kraj 22 0.06249 11.748 1.025e+01
## 4 Plze\362sk\375 kraj 8 0.05277 9.921 -1.921e+00
## 5 Karlovarsk\375 kraj 13 0.02974 5.592 7.408e+00
## 6 \332steck\375 kraj 15 0.08378 15.750 -7.503e-01
## 7 Libereck\375 kraj 7 0.04318 8.118 -1.118e+00
## 8 Kr\341lov\351hradeck\375 kraj 12 0.05444 10.235 1.765e+00
## 9 Pardubick\375 kraj 8 0.05204 9.783 -1.783e+00
## 10 Kraj Vyso\350ina 12 0.05385 10.123 1.877e+00
## 11 Jihomoravsk\375 kraj 17 0.10957 20.599 -3.599e+00
## 12 Olomouck\375 kraj 9 0.06228 11.708 -2.708e+00
## 13 Zl\355nsk\375 kraj 10 0.05785 10.877 -8.766e-01
## 14 Moravskoslezsk\375 kraj 19 0.12455 23.415 -4.415e+00
## 15 <NA> 188 1.00000 188.000 8.882e-16
## statistika
## 1 4.09471
## 2 0.82370
## 3 8.94570
## 4 0.37179
## 5 9.81483
## 6 0.03574
## 7 0.15409
## 8 0.30445
## 9 0.32487
## 10 0.34797
## 11 0.62877
## 12 0.62626
## 13 0.07065
## 14 0.83241
## 15 27.37596Testová 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.