# NMFM301 # 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: ```{r} 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". ## 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. ```{r} tbl <- table(zam$domicil) tbl ``` Celkový počet respondentů je ```{r} sum(tbl) ``` 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$: ```{r} tbl/sum(tbl) ``` Mohli bychom využít i funkce `prop.table`, která počítá totéž: ```{r} prop.table(tbl) ``` 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. ```{r} multitest(tbl,comb=c(-1,0,0,1,1),c0=0.1) ``` 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í. ```{r} multitest(tbl,comb=c(1,1,-1/2,0,0),c0=0) ``` 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. ```{r} deti ``` 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ý). ```{r} dni <- c(31,28,31,30,31,30,31,31,30,31,30,31) sum(dni) ``` Nyní napočítáme pravděpodobnosti měsíců za nulové hypotézy. ```{r} psti <- dni/sum(dni) psti ``` Otestujeme 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. ```{r} chisq.test(deti$poc.deti,p=psti) ``` 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. 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`). ```{r} kuchej.chisq(deti$poc.deti,psti) ``` V 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ě. ### Je úmrtnost ve věku 15-19 let stejná ve všech krajích? Nyní použijeme data `kraje`. ```{r} kraje ``` 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. ```{r} psti <- kraje$obyv.15.19/sum(kraje$obyv.15.19) psti ``` 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. ```{r} chisq.test(kraje$zemr.15.19,p=psti) kuchej.chisq(kraje$zemr.15.19,psti,skup=kraje$kraj) ``` 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.