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:

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.

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

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.185423                         0.049303 
##               Town or small city                  Country village 
##                         0.485888                         0.277242 
##      Farm or home in countryside 
##                         0.002144

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.185423                         0.049303 
##               Town or small city                  Country village 
##                         0.485888                         0.277242 
##      Farm or home in countryside 
##                         0.002144

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)
## 
## 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.

\(\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    19592

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ý).

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.08493 0.07671 0.08493 0.08219 0.08493 0.08219 0.08493 0.08493
##  [9] 0.08219 0.08493 0.08219 0.08493

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.

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-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. 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.0601

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.

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      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.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.12455

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.38, df = 13, p-value = 0.01105
kuchej.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.37596

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.