NMFM301, ZS 2017/18

Cvičenie 8 (praktické cvičenie 2)

(popisné charakteristiky a štatistické testy)



Účelom druhého praktického cvičenia je práca s reálnymi datovými súbormi a aplikácia jednoduchých štatistických metód (počítanie základných výberových popisných charakteristík, konštrukcia intervalov spoľahlivosti pre neznámy parameter a test hypotéz o neznámom parametri) pomocou štatistického softwaru R.

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)




I. Kvantilový diagram

Pre ľubovoľný náhodný výběr \(X_1,\ldots,X_n\) z nejakého rozdělení \(F_X\) používame kvantilový diagram ako neformálnu grafickú metódu slúžiacu k porovnaniu neznámeho rozdělení \(F_X\) s hypotetickým rozdělením \(F_0\) (napr. v prípade potreby overiť predpoklad normality je \(F_0\) štandardne Gaussovo rozdelenie \(N(0,1)\)).

Na os \(y\) sa vykreslia pořádkové statistiky \(X_{(i)}\), pro \(i = 1, \dots, n\). Na os \(x\) se vynášejí aproximace střední hodnoty pořádkových statistik rozdělení \(F_0\), tj. \(F_0^{-1}\bigl(\frac{i}{n+1}\bigr)\) [proč zrovna tohle?].

Jestliže data pocházejí z rozdělení \(F_0\), vykreslené body leží přibližně na přímce \(y=x\). Jestliže data pocházejí z lineární transformace rozdělení \(F_0\), vykreslené body stále leží přibližně na přímce (ale jiné než \(y=x\)). Jestliže data pocházejí z rozdělení odlišného od \(F_0\), vykreslené body tvoří rostoucí křivku. Z tvaru této křivky dokážeme poznat, v čem se rozdělení dat liší od \(F_0\).

V R se dá porovnávat rozdělení dat s rozdělením N(0,1) pomocí funkce qqnorm(). Argumentem funkce je vektor obsahující náhodný výběr.

Vygenerujme 50 pozorování z rozdělení N(10,4) a namalujme normální kvantilový diagram:

x <- rnorm(50,10,2)
qqnorm(x, pch = 21, bg = "red")
qqline(x)

Funkce qqline() namaluje přímku procházející dolním a horním výběrovým kvartilem dat, která usnadnuje čtení grafu.

Nyní vygenerujme 50 pozorování z rozdělení Exp(0.4) a namalujme opět normální kvantilový diagram:

x <- rexp(50,0.4)
qqnorm(x, pch = 21, bg = "red")
qqline(x)

Poznali byste, že tato data nepocházejí z normálního rozdělení?

Úkol:


Vytvořte normální kvantilový diagram pro data z rozdělení \(\Gamma(a,p)\) a sledujte, co se děje, když se zvětšuje hodnota argumentu/parametru \(p\).
  • Vygenerujte 50 pozorování z rozdělení \(\Gamma(a,p)\), kde a=0.4 a p=0.5. Použijte help k funkci rgamma() pro správnou parametrizaci Gama rozdělení.
  • Nakreslete normální kvantilový diagram pro takto vygenerované data.
  • Zvětšujte hodnotu parametru p na 1, 2, 20, 50 a sledujte, co se děje s kvantilovým diagramem. Jaké to má vysvětlení?



II. Práce s daty

V nasledujúcom budeme pracovať so skutočnými datami, ktoré pochádzajú z reálneho experimentu. Datovy súbor je možne stiahnuť na adrese http://msekce.karlin.mff.cuni.cz/~maciak/NMFM301/cv2.RData, prípadne je možné použiť link na data a data priamo načítať do programu R pomocou príkazu.

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

Momentálne máme v programe R k dispozíci datový soubor s názvom cv2. Stručny náhľad na štruktúru datového súboru získame pomocou príkazu

head(cv2)
##         id   pohl              vzdelani                stav  vaha vyska
## 6126 68573   Male            Elementary             Married  57.1 165.0
## 8610 71157 Female    College unfinished             Married  94.5 162.7
## 2886 65182   Male            Elementary Living with partner  71.6 170.3
## 7500 69993 Female            Elementary              Single 118.4 156.6
## 2430 64705 Female    College unfinished             Married  76.7 148.9
## 3591 65919 Female Elementary unfinished             Widowed  49.4 145.8

Data mají celkovo 500 pozorování (řádků) a 6 veličin (sloupců). Túto informáciu zjistíme pomocou funkcie dim():

dim(cv2)
## [1] 500   6

Pozorování tvoří náhodný výběr 500 obyvatel USA. V jednotlivých veličinách o nich máme následující údaje:

Název Význam
id ID jedince
pohl Pohlaví
vzdelani Vzdělání
stav Rodinný stav
vaha Hmotnost [kg]
vyska Výška [cm]

Data si prohlédněte poklepáním na příslušný řádek v okénku “Workspace”. Strukturu dat vypíšeme pomocí funkce str:

str(cv2)
## 'data.frame':    500 obs. of  6 variables:
##  $ id      : num  68573 71157 65182 69993 64705 ...
##  $ pohl    : Factor w/ 2 levels "Male","Female": 1 2 1 2 2 2 1 1 1 1 ...
##  $ vzdelani: Factor w/ 5 levels "Elementary unfinished",..: 2 4 2 2 4 1 3 4 2 4 ...
##  $ stav    : Factor w/ 6 levels "Married","Widowed",..: 1 1 6 5 1 2 3 6 1 1 ...
##  $ vaha    : num  57.1 94.5 71.6 118.4 76.7 ...
##  $ vyska   : num  165 163 170 157 149 ...
##  - attr(*, "na.action")=Class 'omit'  Named int [1:385] 13 32 48 75 83 93 103 105 119 145 ...
##   .. ..- attr(*, "names")= chr [1:385] "30" "62" "87" "130" ...

Základní popisné statistiky všech veličin získáme pomocí funkce summary():

summary(cv2)
##        id            pohl                      vzdelani  
##  Min.   :62172   Male  :259   Elementary unfinished: 40  
##  1st Qu.:64728   Female:241   Elementary           : 82  
##  Median :67260                High school          :100  
##  Mean   :67133                College unfinished   :155  
##  3rd Qu.:69597                College              :123  
##  Max.   :71897                                           
##                   stav          vaha            vyska      
##  Married            :268   Min.   : 35.30   Min.   :139.3  
##  Widowed            : 34   1st Qu.: 65.67   1st Qu.:160.4  
##  Divorced           : 45   Median : 78.60   Median :167.8  
##  Separated          : 18   Mean   : 81.81   Mean   :168.1  
##  Single             : 97   3rd Qu.: 94.08   3rd Qu.:175.3  
##  Living with partner: 38   Max.   :188.50   Max.   :204.5

S datovým souborem se zachází podobně jako s maticí. Můžeme vypsat jeden nebo více řádků pomocí hranaté závorky:

cv2[1:5,]
##         id   pohl           vzdelani                stav  vaha vyska
## 6126 68573   Male         Elementary             Married  57.1 165.0
## 8610 71157 Female College unfinished             Married  94.5 162.7
## 2886 65182   Male         Elementary Living with partner  71.6 170.3
## 7500 69993 Female         Elementary              Single 118.4 156.6
## 2430 64705 Female College unfinished             Married  76.7 148.9

K jednotlivým veličinám máme přístup buď pomocí čísla sloupce (cv2[,2] je vektor hodnot pohlaví) nebo lépe pomocí názvu veličiny napsané za názvem dat a oddělené znakem \(, např `cv2\)pohl`. Výšky prvních pěti obyvatel vypíšeme takto:

cv2$vyska[1:5]
## [1] 165.0 162.7 170.3 156.6 148.9

Hmotnosti lidí se základním vzděláním získáme pomocí podmínky v hranaté závorce (rovnítko musí být zdvojeno):

cv2$vaha[cv2$vzdelani=="Elementary"]
##  [1]  57.1  71.6 118.4  77.6  82.8  62.1  35.3  73.2 130.0  67.0  61.0
## [12]  83.0  78.7 109.0 104.1  66.6  87.1  84.5 137.4  91.2  62.1  63.0
## [23] 124.7  87.2  84.4  86.6  82.4  70.0  76.7  83.7 113.8  90.4  62.1
## [34]  90.6  84.0  64.6  81.0 167.8  71.3 129.2  73.4  83.1  69.8  51.9
## [45]  77.2  88.3 109.3  83.0  78.2  87.9  75.0  94.5  78.1 140.1 101.4
## [56]  51.8  87.4  81.1  70.0 109.7  69.1  73.3  64.5  57.7 110.2  76.1
## [67]  60.2 143.5  62.2  90.2  59.6  88.6  68.2  70.3 104.5  54.5  88.1
## [78]  86.4  65.6 110.1  75.9  64.5



Užitočné


Několik užitočných funkcii, ktoré uľahčia prácu v programe R (napr. funkcie plotCI,plotmeans,ci.asym,ci.t,sign.test).

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



Zkusme se nyní podívat na normální kvantilový diagram výšky.

qqnorm(cv2$vyska, pch = 21, bg = "red")
qqline(cv2$vyska)

Obecně se věří, že výška v populaci má normální rozdělení. Nasvědčuje tomu tento obrázek?

Podívejme se ještě na výšky mužů a žen zvlášť.

qqnorm(cv2$vyska[cv2$pohl=="Male"], pch = 21, bg = "red")
qqline(cv2$vyska[cv2$pohl=="Male"])