NMFM301

Praktické cvičení 2 – Instruktáž


Kvantilový diagram

Mějme náhodný výběr \(X_1,\ldots,X_n\) z rozdělení \(F_X\). Kvantilový diagram je grafická metoda sloužící k neformálnímu porovnání rozdělení \(F_X\) s hypotetickým rozdělením \(F_0\).

Na ose \(y\) se kreslí pořádkové statistiky \(X_{(i)}\). Na ose \(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)
qqline(x)

plot of chunk unnamed-chunk-1

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)
qqline(x)

plot of chunk unnamed-chunk-2

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 argument \(p\).

  1. Vygenerujte 50 pozorování z rozdělení \(\Gamma(a,p)\), kde a=0.4 a p=0.5.
    Musíte dát pozor na správnou parametrizaci Gama rozdělení; použijte volání rgamma(50,rate=0.4,shape=0.5) (argument rate udává hodnotu parametru a, argument shape udává hodnotu parametru p, oba musíte pojmenovat).

  2. Nakreslete normální kvantilový diagram pro tato data.

  3. 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í?

Práce s daty

Nyní se naučíme pracovat s datovými soubory. Načtěte data z pracovního direktoráře nebo přímo z webu příkazem

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

Nyní máte k dispozici datový soubor nazvaný cv2.

Načteme ješte několik funkcí, s kterými dnes budeme pracovat (plotCI,plotmeans,ci.asym,ci.t,sign.test).

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

Data mají 500 pozorování (řádků) a 6 veličin (sloupců). Jejich velikost zjistíme funkcí 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.3   Min.   :139  
##  Widowed            : 34   1st Qu.: 65.7   1st Qu.:160  
##  Divorced           : 45   Median : 78.6   Median :168  
##  Separated          : 18   Mean   : 81.8   Mean   :168  
##  Single             : 97   3rd Qu.: 94.1   3rd Qu.:175  
##  Living with partner: 38   Max.   :188.5   Max.   :204

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

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

qqnorm(cv2$vyska)
qqline(cv2$vyska)

plot of chunk unnamed-chunk-11

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"])
qqline(cv2$vyska[cv2$pohl=="Male"])

plot of chunk unnamed-chunk-12

qqnorm(cv2$vyska[cv2$pohl=="Female"])
qqline(cv2$vyska[cv2$pohl=="Female"])

plot of chunk unnamed-chunk-13

Jsou výšky pro každé pohlaví normální?

Intervaly spolehlivosti pro střední hodnotu

Funkce ci.t počítá 95% interval spolehlivosti pro střední hodnotu založený na t-rozdělení: \[ \left(\overline{X}_{n}-\frac{S_n}{\sqrt{n}}t_{n-1}\Bigl(1-\frac{\alpha}{2}\Bigr), \ \overline{X}_{n}+\frac{S_n}{\sqrt{n}}t_{n-1}\Bigl(1-\frac{\alpha}{2}\Bigr)\right). \] (pro \(\alpha=0.05\)). Je to přesný interval pro normální rozdělení a přibližný pro jakékoli rozdělení s konečným druhým momentem.

smp <- rnorm(20,2,1)
prum <- mean(smp)
ci.t(smp)
## CI.low   Mean  CI.up 
##  1.827  2.247  2.667

Výsledkem funkce ci.t je vektor o třech složkách: dolní konec intervalu, průměr, horní konec intervalu.

Vygenerujme nyní 100 náhodných výběrů o rozsahu 20 se střední hodnotou 2. Nejprve nastavíme požadované vstupy:

nobs <- 20
nvyb <- 100
str.h <- 2

Nyní vygenerujeme 2000 veličin z rozdělení N(2,1) a uspořádáme je do matice 20 x 100. Ve sloupcích bude 100 výběrů o rozsahu 20 z rozdělení N(2,1).

vybery <- rnorm(nobs*nvyb,mean=str.h,sd=1)
data.mat <- matrix(vybery,nrow=nobs,ncol=nvyb)

Na každý sloupec (výběr) spočítáme interval spolehlivosti pro střední hodnotu funkcí ci.t:

vs.ci <- apply(data.mat,2,ci.t)

Zjistíme, kolik intervalů pokrylo skutečnou střední hodnotu (2) a jaká byla průměrná šířka intervalu. Kolik intervalů by ji teoreticky mělo pokrýt?

(pokryti <- sum(vs.ci[1,]<str.h & str.h<vs.ci[3,])/nvyb)
## [1] 0.97
(prum.sirka <- mean(vs.ci[3,]-vs.ci[1,]))
## [1] 0.9479

Nakonec nakreslíme všech 100 spočítaných intervalů do obrázku. Červená čára značí skutečnou střední hodnotu. Vidíte, které intervaly ji nepokryly?

plotCI(vs.ci[2,],uiw=vs.ci[3,]-vs.ci[2,],liw=vs.ci[2,]-vs.ci[1,],
       gap=0.15,sfrac=0.002,ylab="Int. spol. pro stredni hodnotu",xlab="Vyber")
abline(h=str.h,col="red")