Vysoká škola ekonomická v Praze

Vysoká škola ekonomická v Praze Fakulta informatiky a statistiky Studijní program: Kvantitativní metody v ekonomice Studijní obor: Statistika a ekonometrie

Autor bakalářské práce: Petr Marek Vedoucí bakalářské práce: Mgr. Milan Bašta, Ph.D.

GRAFICKÉ ZOBRAZENÍ STATISTICKÝCH DAT V R

školní rok 2012/2013

Prohlášení Prohlašuji, že jsem bakalářskou práci zpracoval samostatně a že jsem uvedl všechny použité prameny a literaturu, ze kterých jsem čerpal.

V Praze dne 14. května 2013

......................................................... Jméno a příjmení studenta

Abstrakt Cílem práce je představit možnosti grafického zobrazení dat a ukázat implementaci těchto zobrazení v R. Grafická zobrazení jsou rozdělena dle počtu a charakteru proměnných. Jednotlivá zobrazení jsou popsána a je ukázán postup jejich tvorby. Zobrazení jsou porovnávána mezi sebou a je diskutován jejich přínos stejně tak jako výhody a nevýhody oproti ostatním. Je předvedena jejich implementace v R. Dále je rozvedena tvorba grafů v R obecně, včetně druhotného přizpůsobování a kombinování těchto grafů. Jsou prezentovány autorem vytvořené funkce umožňující tvorbu některých z méně tradičních grafů.

Klíčová slova R, Grafické zobrazení dat, Grafická analýza, Grafy

Abstract The aim of this thesis is to present ways of visualising data using R. Based on the number and types of variables suitable visualisation methods are presented. These methods are described and their creation is explained. They are further discussed and compared. Implementation of these methods in R is shown. Finally, the ways of customizing and combining graphs in R are presented, including some custom author-created functions.

Keywords R, Data visualisation, Graphical analysis, Graphs

Obsah

Obsah

Obsah 1 Předmluva

1

2 Úvod do R 2.1 O systému R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Práce s R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Grafy v R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2 2 2 4

3 Jednorozměrná data 3.1 Krabičkový graf . . . . . . . . . 3.2 Histogram . . . . . . . . . . . . 3.3 Jádrový odhad hustoty . . . . . 3.4 QQ graf . . . . . . . . . . . . . 3.5 Violin plot . . . . . . . . . . . . 3.6 Grafy pro kategoriální proměnné 3.6.1 Výsečový graf . . . . . 3.6.2 Wordcloud . . . . . . .

. . . . . . . .

5 5 9 13 15 18 20 20 21

. . . . .

23 23 24 25 28 30

. . . .

32 32 33 35 36

6 Univerzální nastavení a úpravy grafů 6.1 Kombinování grafů . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1.1 Více grafů v jednom obrázku . . . . . . . . . . . . . . . . . . . . . . 6.2 Přidání bodů do grafu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

38 38 38 39

4 Dvourozměrná data 4.1 XY graf . . . . . . . . . . . 4.2 Bagplot . . . . . . . . . . . 4.3 Vícenásobný krabičkový graf 4.4 Sectioned density plot . . . . 4.5 Histogramy . . . . . . . . . 5 Vícerozměrná data 5.1 Rainbow plot . . . . . . . . 5.2 Interaktivní 3D graf . . . . . 5.3 Mozaikový graf . . . . . . . 5.4 Vizualizace korelační matice

. . . . .

. . . .

. . . . .

. . . .

. . . . . . . .

. . . . .

. . . .

. . . . . . . .

. . . . .

. . . .

. . . . . . . .

. . . . .

. . . .

. . . . . . . .

. . . . .

. . . .

. . . . . . . .

. . . . .

. . . .

. . . . . . . .

. . . . .

. . . .

. . . . . . . .

. . . . .

. . . .

. . . . . . . .

. . . . .

. . . .

. . . . . . . .

. . . . .

. . . .

. . . . . . . .

. . . . .

. . . .

. . . . . . . .

. . . . .

. . . .

. . . . . . . .

. . . . .

. . . .

. . . . . . . .

. . . . .

. . . .

. . . . . . . .

. . . . .

. . . .

. . . . . . . .

. . . . .

. . . .

. . . . . . . .

. . . . .

. . . .

. . . . . . . .

. . . . .

. . . .

. . . . . . . .

. . . . .

. . . .

. . . . . . . .

. . . . .

. . . .

. . . . . . . .

. . . . .

. . . .

. . . . . . . .

. . . . .

. . . .

. . . . . . . .

. . . . .

. . . .

. . . . . . . .

. . . . .

. . . .

Obsah 6.3

Obsah Přidání čar do grafu . . . . . . . . . . . . 6.3.1 Přidání přímky . . . . . . . . . . 6.3.2 Přidání grafu dle předpisu funkce 6.3.3 Přidání specifických čar . . . . .

. . . .

39 39 41 41

7 Vytvořené funkce 7.1 Rainbow plot funkce . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Funkce pro tvorbu stromu života . . . . . . . . . . . . . . . . . . . . . . . . 7.3 Funkce popisující vznik sectioned density plotu . . . . . . . . . . . . . . . .

44 44 45 46

8 Závěr

47

9 Reference

48

A Příloha A.1 Použitá data . . . . . . . . . . . . A.1.1 Spotřeba automobilů 2012 A.1.2 Produkce plynu v Austrálii A.1.3 Návštěvnost ČR . . . . . A.1.4 Míry úmrtnosti v ČR . . . A.1.5 Střední stav populace ČR . A.1.6 Barva očí a vlasů . . . . . A.1.7 Údaje o automobilech . . A.2 Grafy . . . . . . . . . . . . . . .

51 51 51 51 52 52 53 53 54 54

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . . . . . . .

1

PŘEDMLUVA

1 Předmluva Představme si, že jsme do rukou dostali data, o kterých se chceme něco dozvědět. Než je začneme jakkoliv dále analyzovat, je vhodné si data nějakým způsobem prohlédnout. Spousta vlastností je totiž nejlépe vidět z grafů. Je možné, že u některých dat pomocí grafů zjistíme, že jejich rozdělení neumožňuje určitý druh analýzy, nebo je pro něj krajně nevhodné. Existuje celá řada grafů, které umožňují primárně poznat data a pomoci s (až druhotnou) analýzou. Z relevantních grafů můžeme získat veškeré potřebné informace o tom, jakou část statistického aparátu bychom měli nasadit. Grafy odhalují nesrovnalosti v datech a pomáhají pro ověření předpokladů pro další postup. Celá práce je psaná v duchu prvotního poznání dat grafickými nástroji. Pro tvorbu grafů je použito programové prostředí R [20], které v sekci (2) stručně představím. Znalost R však není pro čtení práce nutná. Po teoretickém úvodu by se měl uživatel zvládnout orientovat v komentovaném kódu1 . Při tvorbě grafu je kód vždy uveden hned vedle grafu tak, že i méně zkušený čtenář zvládne takový graf napodobit a dále přizpůsobit. Hlavním obsahem práce je představení grafů, a to jak pro jednu proměnnou v sekci (3), tak pro dvě proměnné v sekci (4) a pro tři a více proměnných v sekci (5). Graf představím, popíši jaké má výhody a nevýhody a ukáži jeho implementaci v R. Dále uvedu, k čemu je graf vhodné použít a co vše nám může o datech ukázat. Zmíním se také o případných nedostatcích. Informace o většině méně známých grafů, kterými se budu zabývat, jsou dostupné jen v anglické literatuře. Názvy těchto grafů se nesnažím otrocky překládat, ale držím se originálních anglických názvů, jak je standardem. Rozhodl jsem se do sekce (6) vypsat některé univerzální postupy a k nim náležící funkce používané při tvorbě grafů v R, které mají pomoci zejména nezkušenému čtenáři. V sekci (7) jsem shrnul zajímavé a užitečné funkce, které jsem při psaní této práce vytvořil. Pro ilustraci grafů jsem se vždy snažil vybrat co možná nejvhodnější data. Datových souborů je proto hned několik. V textu se na ně odkazuji kurzívou jako například Spotřeba automobilů 2012. Seznam datových souborů a informace o nich jsem zahrnul do přílohy.

1 Vždy

když poprvé použiji některou z pokročilejších či méně známých funkcí, příkládám komentář, popřípadě odkaz do sekce, kde se funkcí přímo zabývám.

1

2 ÚVOD DO R

2 Úvod do R 2.1 O systému R R je volně šířitelné programové prostředí2 podporující velkou škálu operačních systémů - lze jej spustit na *nix systémech (všechny Linux distribuce, MacOSX) i na operačním systému Windows. R je distribuováno jako základní jádro, které je možné rozšiřovat uživateli spravovanými balíčky, které významně rozšiřují jeho možné využití. Jedná se o velmi mocný nástroj s množstvím grafických výstupů, které umožňuje velmi detailně přizpůsobovat. Umí také generovat celou řadu dalších výstupů. R je také programovací jazyk, se kterým programové prostředí R pracuje. Syntaxí je podobný pythonu. Dle tabulky nejpoužívanějších programovacích jazyků (sestavuje společnost Tiobe Software3 ) je jazyk R za březen 2013 na 26. místě. Nutno podotknout, že se jedná o všechny programovací jazyky, nikoliv jen statistické (z těch se ještě na 23. místě umístil SAS, jinak žádný). Je tedy vidět, že R je velmi populární a široce používaný, a to především v zahraničí. Aktuální verze (v době psaní) je 3.0.

2.2 Práce s R V R je možné pracovat rovnou v konzoli, nebo psát R skripty. Konzoli je možné spustit příkazem R na příkazové řádce: petr@sova:~$ R R version 3.0.0 (2013-04-03) -- "Masked Marvel" Copyright (C) 2013 The R Foundation for Statistical Computing Platform: x86_64 -pc-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 2 Tvůrci R jej nevidí primárně jako statistický program, jak je občas prezentován. Preferují označení programové

prostředí, ve kterém jsou implementovány statistické postupy [21]. tabulka je k nahlédnutí na http://www.tiobe.com/index.php/content/paperinfo/tpci/index.html.

3 Celá

2

2.2 Práce s R

2 ÚVOD DO R

'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. >

Speciální znak > uvozuje řádek a je možné za ním psát příkazy4 . Jazyk R má velmi jednoduchou a intuitivní syntax. Lze ho vlastně použít i jako obyčejnou kalkulačku bez jakékoliv předchozí znalosti: > 1 + 2 + 3 [1] 6

> a > b > c > c [1]

> 3^3 [1] 27

= 8 = 3 = a/b 2.666667

Pro tvorbu grafů budeme používat funkce. Funkce voláme jejich jménem společně s argumenty v kulatých závorkách. Například c() pro tvorbu vektoru a sum() pro vypočtení sumy: > vektor = c(0,1,2) > sum(vektor) [1] 3

# toto je komentar

S R se velmi snadno začíná, jelikož obsahuje velmi kvalitní nápovědu. Pro nápovědu ohledně specifické funkce stačí napsat ?jmenoFunkce , např. tedy ?sum . Tato nápověda obsahuje detailní popis parametrů, včetně vysvětlení a příkladů. Kolem R je velmi silná a přátelská komunita5 , takže není problém se při případných potížích obrátit na zkušenější uživatele na internetu6 . Během práce v R se nevyhneme použití externích balíčků. Instalaci externího balíčku lze provést pomocí příkazu install.packages('jmenoBalicku') . To stačí provést jednou. Balíček se stáhne a bude nadále k dispozici. Pro využití funkcionality tohoto balíčku je potřeba jej ve skriptu zahrnout pomoci library('jmenoBalicku') . Balíčků je celá řada. Pro vyhledávání je vhodné použít Google a zadat relevantní klíčová slova (např. R violin plot package“) a nebo vyhledávat na stránce se seznamem všech evidovaných ” balíčků7 . 4 Stejně

jako za $ na *nixu. což R, stejně jako za svůj velký rozmach, vděčí především Linuxu [8]. 6 Asi nejlepší volbou je stránka http://stackoverflow.com/. 7 http://cran.r-project.org/web/packages/available_packages_by_name.html 5 Za

3

2.3 Grafy v R

2 ÚVOD DO R

2.3 Grafy v R Nejpřesnější popis tvorby grafů si vypůjčím od pánů Foxe s Weisbergem [8]. Tvorbu grafů v R přirovnávají k reálnému kreslení grafů na papír. V reálu většinou totiž začneme tak, že si načrtneme osy a určíme měřítko. Poté přidáváme jednotlivé prvky grafu, jako jsou křivky, čáry, body, nadpisy či legenda. Přesně tak se to dělá v R - jen místo ořezané tužky používáme R příkazy.

4

3

JEDNOROZMĚRNÁ DATA

3 Jednorozměrná data Jako jednorozměrná data označujeme taková data, která zachycují hodnoty jedné proměnné. Pro potřeby grafického zobrazování rozdělíme proměnné na kvantitativní a kategoriální. Podle charakteru proměnné použijeme různé typy grafů. Pro kvantitavní proměnné budeme zkoumat především míry polohy (aritmetický průměr, medián, modus a kvantily), míry variability (rozptyl, respektive směrodatná odchylka a mezikvartilové rozpětí) a míry asymetrie (šikmost a špičatost). Tyto míry se pokusíme graficky zobrazit. Budeme zkoumat, zda-li se v datech vyskytují nějaká odlehlá pozorování8 a zda-li jsou rozdělení těchto proměnných jednovrcholová. Využijeme při tom krabičkový graf, histogram a QQ graf. V případě spojitých kvantitativních proměnných se také blíže podíváme na jejich rozdělení, pro která vykreslíme jádrový odhad. U kategoriálních proměnných budeme zkoumat četnost, kterou znázorníme výsečovým grafem a wordcloudem. Především pro základní popis proměnných platí otřepané obrázek vydá za tisíc slov“. Grafy ” pro jednu proměnnou jsou lehce pochopitelné i úplným laikem. Grafy mu umožní na data nahlédnout a bez znalosti výpočtu statistik je popsat.

3.1 Krabičkový graf Krabičkový graf je jedním z nejjednodušších způsobů, jak zobrazit základní informace o kvantitativní proměnné. Sestává z krabičky, vousů“ a jednotlivých bodů. Hranice krabičky jsou ” vymezeny dolním a horním kvartilem. Krabička je protnuta silnou linkou, jejíž poloha je určena hodnotou mediánu. Z krabičky vystupují oběma směry vousy“ znázorňující úrovně nevychý” leného9 maxima a minima [29]. Pomocí bodů jsou znázorněny odlehlé hodnoty10 . Krabičkový graf vykresluje hodnoty základních měr polohy. Vyčteme z něj také informace o variabilitě díky jasně viditelné hodnotě mezikvartilového rozpětí [10]. Z krabičkového grafu se ovšem nedozvíme, jakých hodnot nabývají míry asymetrie. Asymetrii však můžeme odhalit na základě vykreslených odlehlých pozorování. Krabičkový graf je vhodné doplnit dokreslením průměru, abychom měli přesnější představu o jeho vzdálenosti od mediánu. Výraznější vzdálenost totiž může poukazovat na kladné či záporné zešikmení [6]. 8 Tedy

taková pozorování, jejichž hodnoty se zdají být v numerickém nesouladu se zbytkem dat [26]. maxima a minima z pozorování, která nejsou označena jako odlehlá. 10 Nejčastěji (v R pak ve výchozím nastavení) se za odlehlou hodnotu považuje taková hodnota, která je vzdá” lena“ od hrany krabičky (tj. od dolního, respektive horního kvartilu) o více než jedenapůlnásobek mezikvartilového rozpětí [20]. 9 Tj.

5

3.1 Krabičkový graf

3


> > > > > >

data = read.csv('data/car_fuel.csv') boxplot(data$LPH) # na.rm smaze NA (prazdne) zaznamy prumer = mean(data$LPH , na.rm = TRUE) # pro symbol krizku volime pch = 3 points(prumer , pch = 3)

5

10

15

20

V R je možné krabičkový graf vykreslit funkcí boxplot() . Pro dokreslení průměru můžeme použít funkci points() 11 . Ukažme si krabičkový graf ve výchozím nastavení (s přidaným křížkem pro průměr) pro proměnnou spotřeba na 100 kilometrů z dat Spotřeba automobilů 2012 na obrázku 3.1.1.

Obrázek 3.1.1: Krabičkový graf s křížkem v místě průměru R nabízí mnoho možností, jak krabičkový graf přizpůsobit. Do tabulky 3.1.1 jsem vypsal užitečné parametry, kterými můžeme některé z přizpůsobení provést. Tyto parametry se nám budou hodit především později v sekci (4.3) při práci s vícenásobným krabičkovým grafem. Kombinací těchto parametrů můžeme pro stejnou proměnnou získat například krabičkový graf na obrázku 3.1.2. Při grafické analýze je vhodné zkoumat data i pomocí jiných grafů. Jelikož je krabičkový graf jakýmsi kompromisem mezi jednoduchostí zobrazení a reflektováním maxima informací o proměnné12 , je potřeba dávat pozor při porovnávání různých rozdělení. Krabičkový graf někdy zobrazí velmi odlišná rozdělení velmi podobně [3]. Příkladem mohou být vygenerovaná data na obrázku 3.1.3, kde jsem pro porovnání přidal histogramy13 .

11 Více

o této funkci píši v sekci (6.2). se dá zobecnit víceméně pro všechna grafická zobrazení. Většinou totiž graficky zobrazujeme v méně rozměrech (v drtivé většině ve dvou) než je počet proměných [3]. 13 V kódu je použita funkce layout() - více o této funkci v (6.1.1). Funkce pro tvorbu histogramu hist() je společně s grafem jako takovým popsána v (3.2). 12 To

6


3


Tabulka 3.1.1: Specifické parametry pro tvorbu krabičkového grafu Parametr

Výchozí hodnota

Popis

horizontal

FALSE

pokud je TRUE , tak se krabičkový graf vykreslí na šířku, ne na výšku

range

1.5

nastavení délky“ krabičkového grafu; určuje, ” kolika-násobek mezikvartilového rozpětí je hranicí pro odlehlá pozorování

varwidth

FALSE

pokud je TRUE , tak jsou krabičky“ v grafu široké jako ” druhá odmocnina počtu pozorování v této skupině

> data = read.csv('data/car_fuel.csv') > boxplot(data$LPH , + horizontal = TRUE, + varwidth = TRUE, + range = 2 + )

5

10

15

20

Obrázek 3.1.2: Přizpůsobený krabičkový graf

7


3


Krabičkový graf je velmi užitečným nástrojem při základním popisu zkoumanné proměnné. Není to však jeho jediné využití. Jak uvidíme v sekci (4.3), lze krabičkový graf snadno použít také pro zkoumání závislostí mezi kvantitativní a kategoriální proměnnou14 . Krabičkový graf se pro množství informací, které obsahuje, stal základem dalších (pokročilejších) grafů, jako jsou například violin plot či bagplot15 . Při porovnávání různých rozdělení je potřeba dávat pozor, jelikož použití krabičkového grafu může být poněkud zavádějící, jak je vidět na obrázku 3.1.3. > > > > + + + > > + + + > > > + + + + + + + + + + > >

# pro reprodukovatelny vysledek set.seed(151252) # vykreslime 4 radky layout( mat = matrix(c(1,2,3,4)), heights = c(4,1,4,1)/10 ) # dve ruzna rozdeleni data = list( c(rnorm(200), rnorm(200, 5)), runif(1000) ) # funkce pro vykresleni histogramu # a krabickoveho grafu drawBoxAndHist = function(x) { lim = c(min(x), max(x)) hist(x, freq = FALSE, main = NULL, xlim = lim, axes = FALSE ) boxplot(x, horizontal = TRUE, frame = FALSE, ylim = lim, axes = FALSE ) } # vykoname funkci pro obe rozdeleni r = lapply(data, FUN = drawBoxAndHist)

Obrázek 3.1.3: Pomocí krabičkového grafu mnohdy nedokážeme rozeznat různé tvary rozdělení

14 Dokonce 15 Více

je to jeden z nejlepších nástrojů [29]. v sekci (3.5), resp. (4.2).

8

3.2 Histogram

3


3.2 Histogram Histogram patří mezi nejužívanější grafy pro popis jedné proměnné, jejíž rozdělení graficky znázorňuje pomocí sloupcového diagramu se stejně širokými sloupci [4]. Šířka sloupců vyjadřuje šířku intervalů, do kterých jsou hodnoty proměnné rozděleny. Výška sloupce pak znázorňuje četnosti těchto hodnot v jednotlivých intervalech. V R se pro tvorbu histogramu používá funkce hist() . Výchozí nastavení histogramu pro náhodný výběru z normálního rozdělení vidíme na grafu 3.2.1.

15

> data = rnorm(150) > h = hist(data)

0

5

10

Frequency

20

25

30

Histogram of data

−3

−2

−1

0

1

2

data

Obrázek 3.2.1: Histogram ve výchozím nastavení Při tvorbě histogramu je potřeba zvážit, do kolika intervalů budeme hodnoty proměnné dělit, respektive jaká bude šířka těchto intervalů. Volba šířky je naprosto zásadní a její důležitost roste s rozsahem souboru. Při špatně zvolené šířce intervalů nám histogram dokáže výrazně zkreslit představu o datech. Možností, jak volit šířku intervalu je hned několik. Výchozím nastavením pro funkci hist() je Sturgesovo pravidlo, které je jedním z prvních pravidel pro volbu šířky16 intervalů vůbec [28]. Sturges navrhuje volit počet intervalů k jako: k = 1 + log2 n

(1)

kde n je rozsah výběru. Šířku intervalu z (1) můžeme vyjádřit jako: hst = 16 Ve

R 1 + log2 n

skutečnosti se jedná spíše o pravidlo pro počet skupin, jak píše Scott [23].

9

(2)

3.2 Histogram

3


kde: n . . . . . . . . . . rozsah výběru, R . . . . . . . . . variační rozpětí, tj. xmax − xmin . Sturgesovo pravidlo není vhodné pro velké výběry, které má tendenci příliš vyhlazovat použitím nedostatečného počtu intervalů [23], jak uvidíme dále. Alternativou může být pravidlo Scottovo: 3.5 · σˆ hsc = 1/3 (3) n nebo Freedman-Diaconisovo: IQR( x ) hfd = 2 · , (4) n1/3 kde: n . . . . . . . . . . rozsah výběru, σˆ . . . . . . . . . .výběrová směrodatná odchylka, IQR( x ) . . . mezikvartilové rozpětí. Ukažme si na jednoduchém příkladu všechna tři pravidla. Vezměme data Spotřeba automobilů 2012 a vytvořme histogram pro spotřebu v litrech na sto kilometrů. Vykreslíme čtyři histo√ gramy - jeden pro každé pravidlo a jeden pro počet intervalů17 roven k = n, tedy odmocnině z počtu pozorování18 . scott 150 100 0

50

Frequency

200 100 0

Frequency

300

sturges

5

10

15

20

5

10

fd

20

100 60 0 20

60

Frequency

100

n

0 20

Frequency

15

5

10

15

20

5

10

15

> > + + + > > + + + + > + + + +

data = read.csv('data/car_fuel.csv') par( mfrow = c(2,2), cex = 1.2, mar = c(3,3,3,1) ) rules = c('sturges', 'scott', 'fd') for (rule in rules) hist( data$LPH , breaks = rule, main = rule ) hist( data$LPH , breaks = sqrt(length(data$LPH)), main = bquote(sqrt('n')) )

20

Obrázek 3.2.2: Důležitost správné šířky intervalu pro histogram 17 Je

nutné podotknout, že jak pojmenovaná“ pravidla, tak zadání čísla bere funkce hist() pouze jako orien” tační hodnotu. Pokud chceme přesný počet intervalů je potřeba zadat vektor s krajními body těchto intervalů. 18 Taková hodnota je výchozím nastavením například v MS Excel.

10

3.2 Histogram

3


Obrázek 3.2.2 nám ukazuje, že menší (Sturgesův) počet intervalů pro velký výběr skutečně histogram příliš vyhladil“ a zatajil tak důležitou informaci o dvouvrcholovosti rozdělení. ” Tabulka 3.2.1: Specifické parametry pro tvorbu histogramu Parametr breaks

Výchozí hodnota Sturges

Popis

počet intervalů, do kterých jsou data rozdělena je hned několik možností nastavení: a) předdefinovaná pravidla: Sturges , Scott , Freedman-Diaconis b) přirozené číslo, udávající celkový počet intervalů c) vektor s body, které určují, kde se budou intervaly lámat - do těchto intervalů budou data skutečně rozdělena, předchozi dva bere funkce hist() jako přibližné

freq

TRUE

- zobrazí na ose y absolutní četnosti FALSE - zobrazí se hodnoty jádrového odhadu (více o jádrových odhadech dále v sekci (3.3)) TRUE

Správně zkonstruovaný histogram (tedy histogram se správnou šířkou intervalů) podává informace o tvaru rozdělení, z kterého můžeme přibližně vyčíst jak míry polohy, tak míry variability. Pomocí histogramu můžeme také odhalit přítomnost vybočujících pozorování. Pokud se v histogramu vyskytuje skupina prázdných intervalů, poukazuje to na odlehlost některých z hodnot. Takový příklad můžeme vidět na obrázku 3.2.3, do kterého jsem pro názornost přidal krabičkový graf. Zejména u vícenásobných grafů se nám může hodit vykreslit otočený histogram. Histogram v R bohužel přímo takové nastavení nemá. Histogramu je potřeba přidat atribut plot s hodnotou FALSE , čímž zamezíme vykreslení. Jeho obsah poté vykreslíme pomocí sloupcového diagramu jako na obrázku 3.2.4. Jak uvidíme dále v (4.5), je otočený histogram vhodný pro porovnávání dvou skupin dat mezi sebou, pro což se běžně používá v demografii.

11

3.2 Histogram

3


60

Histogram of data

30 20 10 0

Frequency

40

50

> set.seed(1) > data = c(rnorm(150), runif(5, 6, 8)) > data.max = max(hist(data, plot = → → FALSE)$counts) > hist(data, ylim = c(-4, data.max)) > boxplot(data, varwidth = TRUE, + horizontal = TRUE, at = -3, + add = TRUE + )

−2

0

2

4

6

8

data

Obrázek 3.2.3: Detekce odlehlých pozorování pomocí histogramu

> histHoriz = function(data, breaks = → → 'Sturges', ...) { + # histogram nevykreslime + h = hist(data, + plot = FALSE, + breaks = breaks + ) + # ale vykreslime sloupcovy graf + barplot(h$density , + space = 0, horiz = TRUE, + col = NA, axes = FALSE + ) + } > > histHoriz(rnorm(500))

Obrázek 3.2.4: Otočený histogram

12

3.3 Jádrový odhad hustoty

3


3.3 Jádrový odhad hustoty Jádrový odhad je neparametrickým typem odhadu hustoty rozdělení proměnné. Je to funkce proměnné x o předpisu: ( ) n 1 x − X i fˆh ( x ) = K (5) nh i∑ h =1 kde: K . . . . . . . . . jádrová funkce, symerická kolem nuly, pro kterou platí h . . . . . . . . . . šířka vyhlazovacího okna, pro kterou platí h > 0.

∫

K ( x )dx = 1,

Jedná se vlastně o vážený klouzavý průměr [22]. Jádrový odhad je vizuálně přesnější než histogram, jelikož je spojitý a umožňuje pozorovat i ty nejmenší změny v průběhu rozdělení [5]. Stejně jako u histogramu a volby šířky intervalů je i zda naprosto zásadní volba šířky vyhlazovacího okna a volba jádrové funkce. Standardní sada jádrových funkcí je v R definována jako: > eval(formals(density.default)$kernel) [1] "gaussian" "epanechnikov" "rectangular" [6] "cosine" "optcosine"

"triangular"

"biweight"

Průběh těchto funkcí jsem pro názornost vykreslil v příloze A.2.1. V R provedeme jádrový odhad pomocí funkce density() , jejíž výstup použijeme jako zdroj dat pro univerzální funkci plot() . Jádrový odhad hustoty výběru z normálního rozdělení o rozsahu n = 150 ve výchozím nastavení (tedy Gaussovskou jádrovou funkcí a šířkou vyhlazovacího okna nrd0 - zjednodušeného Scottova pravidla19 [24].) vidíme na obrázku 3.3.1. Jak správně zvolit šířku vyhlazovacího okna stejně tak jako volbu vhodné jádrové funkce velmi detailně rozebírá Sheather [25]. V této práci se omezíme na jejich ilustraci. Do obrázku 3.3.2 jsem vykreslil20 jádrový odhad pro výběr z normovaného normálního rozdělení o rozsahu n = 150 při použití Gaussovské jádrové funkce a různých šířek vyhlazovacího okna. Užití jádrového odhadu je velmi podobné jako u histogramu (3.2), a proto se tyto dva odhady často vykreslují společně (za předpokladu, že histogram vykreslujeme také jako odhad hustoty - tj. s freq = FALSE ), jako na obrázku 3.3.3.

19 Více

v ?bw.nrd0 . kódu jsem použil funkci rug() , která do spodní části grafu vykreslí jakousi rohož“, která znázorňuje ” reálné umístění hodnot proměnné. Pro přidání čar do grafu jsem použil funkci lines() , více o ní se rozepíši v sekci (6.3).

20 V

13

3.3 Jádrový odhad hustoty

3


0.2

> > > >

set.seed(1) data = rnorm(150) odhad_hustoty = density(data) plot(odhad_hustoty)

0.0

0.1

Density

0.3

0.4

density.default(x = data)

−3

−2

−1

0

1

2

3

N = 150 Bandwidth = 0.2907

Obrázek 3.3.1: Jádrový odhad hustoty ve výchozím nastavení

0.3 0.2 0.1 0.0

dnorm(x)

0.4

0.5

dnorm 1.5 1 0.5 0.25 0.125

−3

−2

−1

0

1

2

3

> > > > > + > + + + + > + + + + >

set.seed(1) data = rnorm(150) x = seq(-3.1, 3.1, 0.001) bws = c(1.5, 1, 0.5, 0.25, 0.125) plot(x, dnorm(x), ylim=c(0, 0.55), type = 'l', lwd = 2) for (i in seq_along(bws)) lines( density(data, bw = bws[i]), col = i+1, lty = 2 ) legend('topright', lty = c(1, rep(2, length(bws))), col = 1:(length(bws)+1), legend = c('dnorm', bws) ) rug(data)

x

Obrázek 3.3.2: Různé šířky vyhlazovacího okna pro stejnou jádrovou funkci

14

3.4 QQ graf

3


0.4

Histogram of data

Density

0.0

0.1

0.2

0.3

> data = rnorm(150) > > hist(data, freq = FALSE) > lines(density(data), col = 'red')

−2

−1

0

1

2

data

Obrázek 3.3.3: Jádrový odhad hustoty vkreslený do histogramu Stejně jako histogram umožňuje jádrový odhad (při korektně zvolených parametrech) nahlédnout na rozdělení zkoumané proměnné. Z tvaru rozdělení poznáme polohu a variabilitu proměnné. Z grafu můžeme vyčíst, zda se v datech objevují odlehlá pozorování21 . Tvar rozdělení ukazuje na případné zešikmení a na špičatost.

3.4 QQ graf QQ graf22 je graf porovnávající dvě rozdělení na základě hodnot kvantilů [15]. Jednotlivé sobě-odpovídající kvantily (tj. například medián jednoho a medián druhého rozdělení) jsou postupně vykreslovány do grafu jako body, jejichž x, resp. y souřadnice odpovidají hodnotě kvantilů prvního, resp. druhého rozdělení. Pomocí QQ grafu tedy narozdíl od ostatních grafů v této sekci neprovádíme popis proměnné, ale porovnání jejího rozdělení s rozdělením jiným (buď rozdělením jiné proměnné, nebo častěji s rozdělením teoretickým). Z QQ grafu nepoznáme, jaké jsou míry polohy, variability či symetrie. Poznáme však, jestli se liší od měr druhého rozdělení. Poznáme také, jestli se mezi hodnotami vyskytují odlehlá pozorování. Nejčastěji se QQ graf používá pro porovnání s normálním rozdělením23 (v R 21 Stejně

qqnorm()

), jako

jako u histogramu - pokud se mezi dvěma shluky dat objevuje prázdná“ část grafu, poukazuje to na ” odlehlost jednoho ze shluků. 22 Také označován jako kvantil-kvantil graf. 23 Zejména pak zkoumáme-li rozdělení náhodné složky v regresní analýze.

15

3.4 QQ graf

3


na obrázku 3.4.1, kde porovnáváme náhodný výběr z normálního rozdělení se samotným teoretickým rozdělením.

0

> data = rnorm(200) > qqnorm(data)

−3

−2

−1

Sample Quantiles

1

2

3

Normal Q−Q Plot

−3

−2

−1

0

1

2

3

Theoretical Quantiles

Obrázek 3.4.1: QQ graf pro náhodný výběr z normálního rozdělení V případě, že jsou hodnoty kvantilu přibližně stejně umístěné, leží body v QQ grafu na přímce y = x. Pokud hodnoty leží na jiné přímce, stále hovoříme o jejich lineární závislosti a tím pádem podobnosti rozdělení. Pokud body na přímce neleží, poukazuje to na neshodu rozdělení. Jednotlivými druhy odchýlení od přímky se budeme zabývat dále. Pokud mezi sebou chceme proměnnou porovnat s jiným rozdělením, než je rozdělení normální, použijeme funkci qqplot() . Porovnejme například rozdělení spotřeby automobilů z dat Spotřeba automobilů 2012 s náhodným výběrem z rovnoměrného rozdělení o stejném rozsahu. Na obrázku 3.4.2 pozorujeme zcela odlišný průběh rozdělení. Jedná se tedy o ukázkový případ QQ grafu různých rozdělení. Hodnoty kvantilů se rovnají jen v extrémech a mediánu. Rozdělení mají tedy naprosto odlišný průběh. Podle tvaru QQ grafu pro normální rozdělení můžeme rozpoznat zásadní informace o rozdělení zkoumané proměnné. Pro usnadnění můžeme použít funkci qqline() , která vykreslí přímku propojující dolní a horní kvartil teoretického rozdělení. Ta nám umožní pozorovat celkovou vzdálenost hodnot od jejich teoretických protějšků. Také nám usnadní identifikaci odlehlých pozorování. Některé z běžných tvarů QQ grafu jsem včetně interpretace shrnul do obrázku 3.4.3.

16

3.4 QQ graf

3

data = read.csv('data/car_fuel.csv') spotreba = data$LPH[complete.cases(data$LPH)] random = runif( length(spotreba), min = min(spotreba), max = max(spotreba) ) qqplot(spotreba , random)

5

10

random

15

20

> > + > + + + + > >


5

10

15

20

spotreba

Obrázek 3.4.2: QQ graf pro spotřebu automobilů a výběr z rovnoměrného rozdělení

−2

−1

0

1

2

3 1

3

−3

−1

0

1

2

Zaporne zesikmeni

−1

0

1

2

−1 −3

15

−2

3

−3

1

−1

0

1

2

3

2

3

5 10 0

Sample Quantiles

10 6 2

0

−2

Odlehla pozorovani

−2 −3 −2 −1

3

1

Kladne zesikmeni Sample Quantiles


Bimodalita Sample Quantiles

−2


0 5 −3

−3 −1

Sample Quantiles

1 −1 −3

Sample Quantiles

Rovnomerne rozdeleni

−3

Sample Quantiles

Normalni rozdeleni

−3

−2

−1

0

1

2

3

> + + > + > > + > + > + > + > + > +

qqn = function(x, y) { qqnorm(x, main = y); qqline(x) } par(mfrow = c(3,2), mar = c(3,3,3,3), cex = 1) set.seed(1) qqn(rnorm(300), 'Normalni rozdeleni') qqn(runif(300, -3, 3), 'Rovnomerne rozdeleni') qqn(rchisq(300, 5), 'Kladne zesikmeni') qqn(c(rnorm(100), runif(200, 0, 2)), 'Zaporne zesikmeni') qqn(c(rnorm(200), rnorm(200, 8)), 'Bimodalita') qqn(c(rnorm(290), runif(10, 15, 16)), 'Odlehla pozorovani')

Obrázek 3.4.3: Bežné tvary QQ grafu společně s jejich významem

17

3.5 Violin plot

3


3.5 Violin plot Violin plot24 je kombinací krabičkového grafu a jádrového odhadu hustoty [10]. Spojuje pozitivní vlastnosti obou grafů a umožňuje nám pracovat jen s jedním grafem namísto dvou. Odstraňuje dříve zmíněné nevýhody samostatného použití těchto grafů.

2

Pro R existuje vioplot balíček25 [1] s funkcí violinplot() , vytvářející tento typ grafu. Na obrázku 3.5.1 vidíme violin plot ve výchozím nastavení. Je samozřejmě velmi podobný krabičkovému grafu. Tečka udává polohu mediánu, délka krabičky“ pak mezikvartilové rozpětí. ” Barevná plocha vyjadřuje křivku jádrového odhadu hustoty rozdělení.

set.seed(1) # pro tvorbu violin plotu # je potreba balicek library('vioplot') vioplot(rnorm(150))

−2

−1

0

1

> > > > >

1

Obrázek 3.5.1: Violin plot Rozíření o jádrový odhad je naprosto zásadní a řeší problém, který znázorňuji na obrázku 3.1.3. Udělejme nyní obdobný obrázek - jen namísto krabičkového grafu použijeme violin plot. Výsledek na obrázku 3.5.2 nám ukazuje, že violin plot umí spolehlivě rozlišit různé tvary rozdělení.

24 Šlo

by přeložit jako houslový graf“, ale nepřekládá se. ” informací o balíčku lze získat na http://cran.r-project.org/web/packages/vioplot/index.html.

25 Více

18

3.5 Violin plot

3

1

x

−2

0

2

4

6

8

1

x

0.0

0.2

0.4

0.6

0.8

> > > + + + > > + + + > > > + + + + + + + > > +


set.seed(151252) library('vioplot') layout( mat = matrix(c(1,2,3,4)), heights = c(3,1,3,1)/10 ) # dve ruzna rozdeleni data = list( c(rnorm(200), rnorm(200, 5)), runif(1000) ) # funkce pro vykresleni histogramu # a violin plotu drawViolinAndHist = function(x) { lim = c(min(x), max(x)) hist(x, freq = FALSE, main = NULL, xlim = lim, axes = FALSE ) vioplot(x, horizontal = TRUE) } # vykoname funkci pro obe rozdeleni r = lapply(data, FUN = drawViolinAndHist)

1.0

Obrázek 3.5.2: Pomocí violin plotu dokážeme na rozdíl od krabičkového grafu rozeznat různá rozdělení

19

3.6 Grafy pro kategoriální proměnné

3


−5

0

5

10

Na obrázku 3.5.3 je vidět, kolik informací navíc nám o zkoumané proměnné violin plot oproti krabičkovému grafu dává.

rovnomerne

normalni

bimodalni

rovnomerne

normalni

−5

0

5

10

bimodalni

> > > > > + + + + + > > + + + + + + +

set.seed(151252) library('vioplot') par(mfrow = c(2,1), mar = c(3,3,1,3)) # tri ruzna rozdeleni data = data.frame( bimodalni = c(rnorm(200), rnorm(200, 5)), rovnomerne = runif(400, -2.7, 7.3), normalni = rnorm(400, 2.6, 3.4) ) boxplot(data) with(data, { vioplot( bimodalni , rovnomerne , normalni , names = colnames(data) ) })

Obrázek 3.5.3: Violin plot nám dává daleko více informací o rozdělení proměnné než krabičkový graf Z violin plotu vyčteme veškeré podstatné informace o proměnné. Přimo vidíme hodnoty měr polohy, z tvaru rozdělení a velikosti mezikvartilového rozpětí vidíme variabilitu a případnou asymetrii. Violin plot je proto dle mého názoru nejlepším nástrojem pro popis spojité kvantitativní proměnné.

3.6 Grafy pro kategoriální proměnné 3.6.1 Výsečový graf V případě, že zkoumáme kategoriální proměnnou, můžeme pro základní popis využít výsečový graf. Výsečový graf nám poskytuje rychlý přehled o relativních četnostech kategorií. Právě proto, že ukazuje relativní zastoupení, může být výsečový graf velmi zavádějící. Vždy bychom k němu měli dodávat, z jak velkého vzorku je tvořen a na jak velkou populaci se snažíme zobecnit výsledek. Ve výsečovém grafu nejsou rozsah výběru ani velikost populace vidět.

20


3


Na to se mnohdy zapomíná a dochází k nesmyslným závěrům nad nedostatečným počtem pozorováním. Podívejme se například na zastoupení značek automobilů v datech Spotřeba automobilů 2012. Pro výsečový graf musíme vypočítat absolutní četnosti, ke kterým přiřadíme náležící názvy. Pokud je jako v případě těchto dat mnoho kategorií, je vhodné (například porovnáním četností kategorií s mediánem četností) zobrazit popisky jen u nejčetnějších názvů. Výsečový graf vidíme na obrázku 3.6.1.

Ford Motor Company

Chrysler Group LLC

BMW

General Motors

Audi Honda Volkswagen

Hyundai Kia

Toyota

MAZDA Subaru

Mercedes−Benz Nissan

Porsche

> > > > > > > + + + + > >

data = read.csv('data/car_fuel.csv') # ziskame absolutni cetnosti cetnosti = table(data$Mfr.Name) # pomoci ifelse vymazeme ta jmena, # ktera maji cetnosti nizsi # nez median znacky = ifelse( cetnosti > median(cetnosti), row.names(cetnosti), '' ) pie(cetnosti , labels = znacky)

Obrázek 3.6.1: Výsečový graf pro četnosti výrobců automobilů Výsečový graf se pro výše zmíněný problém používá spíše v nestatistické literatuře. 3.6.2 Wordcloud Spíše pro zajímavost přidávám zejména na internetu oblíbené zobrazení četnosti slov - wordcloud26 . Nalezená slova jsou tisknuta proporčně k jejich četnosti v textu a uspořádána do jakéhosi oblaku“. To nám umožňuje již na první pohled vidět, které hodnoty jsou nejčetnější. ” Využitím a významem je to vlastně obdoba výsečového grafu, jen s jiným grafickým znázorněním. Nejjednodušší cestou k vytvoření tohoto zobrazení v R je použití balíčku wordcloud [7]. Ten obsahuje stejnojmennou funkci. Zápis je obdobný jako u výsečového grafu.

26 Též

tag cloud; ani jeden termín se do češtiny nepřekládá.

21

Land Rover Maserati

Hyundai


Mitsubishi Motors NA

Jaguar Cars

aston martin

Mercedes−Benz Ferrari

Kia Audi Nissan Lotus

Ford Motor Company

General Motors BMW Suzuki Volvo

SubaruRolls−Royce

Toyota

MAZDA Bentley

Porsche Honda

Chrysler Group LLC Volkswagen

Mitsubishi Motors Co

Saab Cars North America

3

> > > > > > > + > > + + + +


set.seed(1) library('wordcloud') data = read.csv('data/car_fuel.csv') counts = table(data$Mfr.Name) labels = row.names(counts) colors = rainbow(length(labels), v = 0.7) wordcloud(labels, counts , scale = c(5, 1), colors = colors , random.order = FALSE )

Obrázek 3.6.2: Wordcloud Na příkladu wordcloudu na obrázku 3.6.2 jsem použil zajímavou funkci rainbow , umožňující tvorbu vektoru s barvami duhy. Více o této funkci v nápovědě ?rainbow .

22

4

DVOUROZMĚRNÁ DATA

4 Dvourozměrná data U dvourozměrných dat můžeme kromě základních měr a rozdělení proměnných (jako v předchozí sekci) zkoumat také závislost jedné proměnné na druhé. Pokud u proměnných objevíme závislost, bude nás zajímat, jaký je směr této závislosti a jestli tuto závislost splňují všechna pozorování, respektive jestli se jí některá vymykají (tj. budeme hledat obdobu odlehlých pozorování pro jednu proměnnou). Grafická zobrazení rozdělíme podle charakteru zkoumaných proměnných. Nejdříve se budeme zabývat dvěmi kvantitativními proměnnými, pro které vykreslíme XY graf a bagplot. Poté se přesuneme ke vztahu kvantitativní a kategoriální proměnné, pro které využijeme vícenásobný krabičkový graf a sectioned density plot. Vztah dvou kategoriálních proměnných popíšeme až v sekci (5) Vícerozměrná data pomocí mozaikového grafu, jelikož užití a interpretace pro 2 rozměry je stejná jako pro n > 2 rozměrů. Na závěr této sekce porovnáme rozdělení dvou proměnných pomocí dvou histogramů. Opět se v duchu této práce vyhneme číselné analýze (ať už analýze rozptylu či regresní/korelační analýze). Jak uvidíme dále, není číselná analýza pro rozpoznání závislosti nutná.

4.1 XY graf Jedním ze základních typů grafů je XY graf. Znázorňuje průběh jedné z proměnných na základě hodnot té druhé. Na osu x zpravidla řadíme vysvětlující proměnnou a na osu y proměnnou vysvětlovanou. V R určujeme směr závislosti pořadím proměnných při vykreslování grafu, např. tedy plot(y ∼ x) - čteme jako y vysvětlované pomocí x [8]. Tento typ grafu se nejčastěji používá jako první nástroj pro hrubou představu o podobě dat. Podívejme se například na obrázku 4.1.1 na hodnoty spotřeby na 100 kilometrů v závislosti na objemu motoru z dat Spotřeba automobilů 2012. Dle tvaru, kterého nabývá shluk bodů, můžeme uvažovat závislost proměnných. U předchozího příkladu 4.1.1 se zřejmě jedná o přímou závislost, jelikož body leží okolo hypotetické přímky y = ax + b; a > 0, jejíž předpis bychom mohli získat použitím metody nejmenších čtverců. Z XY grafu lze dobře vyčíst, zda-li z této závislosti některé body vybočují. V obrázku 4.1.1 žádné takové body nejsou. Byly by to takové body, které by měly při nízkém objemu motoru velkou spotřebu a naopak. Vystoupily by tak z hlavního shluku bodů.

23

4

DVOUROZMĚRNÁ DATA

15

20

4.2 Bagplot

5

10

LPH

> data = read.csv('data/car_fuel.csv') > plot(LPH ~ Engine, data = data)

0

2

4

6

8

Engine

Obrázek 4.1.1: XY graf vykreslující hodnoty spotřeby na 100 kilometrů v závislosti na objemu motoru

4.2 Bagplot Bagplot27 je zobecněním jednorozměrného krabičkového grafu do dvou rozměrů [19]. Pro pochopení jeho podstaty je potřeba definovat, co je to hloubka bodu v poloprostoru. Hloubku můžeme z dvourozměrného28 grafu jednoduše vyčíst jako minimum bodů, které jsme schopni zahrnout uzavřenou polorovinou s hraniční přímkou procházející zkoumaným bodem. Čím blíže jsme pomyslnému středu všech bodů, tím těžší“ je zahrnout málo bodů a hloubka tím ” pádem roste. Díky hloubce můžeme definovat dvourozměrný medián jako bod s největší hloubkou [19]. Pro ilustraci jsem do přílohy A.2.2 vykreslil jednoduchý příklad zabývající se hloubkou bodů. Podrobněji se touto tématikou zabývá Tukey [27]. Bagplot vychází ze základního XY grafu, který dále rozšiřuje. Hlavní částí bagplotu je bag obsahující 50% prostředních“ pozorování (tj. pozorování s nejvyšší hloubkou), který je na ” grafu vyznačen nejtmavší barvou. V bagu je i dvourozměrný medián, znázorněný bodem (v R hvězdičkou). Další částí je fence, který odděluje odlehlá pozorování od zbytku souboru. Jeho rozměry jsou oproti bagu trojnásobné (resp. vzdálnost vrcholů je od mediánu oproti vrcholům bagu trojnásobná). Vybarven je světlejší barvou. V grafu jsou navíc vyznačena odlehlá pozorování. Pro R existuje balíček aplpack [30], obsahující funkci 27 Nepřekládá 28 Více

se. rozměry se ve spojeni s bagplotem zabývat nebudeme.

24

bagplot()

. Vykresleme tedy bagplot

4.3 Vícenásobný krabičkový graf

4

DVOUROZMĚRNÁ DATA

pro spotřebu na 100 kilometrů a objem motoru z dat Spotřeba automobilů 201229 .

15 5

10

Spotreba na 100 km

20

> library('aplpack') > data = read.csv('data/car_fuel.csv') > > bagplot( + data$Engine , + data$LPH , + xlab = 'Objem motoru', + ylab = 'Spotreba na 100 km', + na.rm = TRUE, + show.whiskers = FALSE + ) [1] "Warning: NA elements have been → → removed!!" 0

2

4

6

8

Objem motoru

Obrázek 4.2.1: Bagplot pro spotřebu automobilů a objem motoru

Tabulka 4.2.1: Specifické parametry pro tvorbu bagplotu Parametr

Výchozí hodnota

factor

3

na.rm

FALSE

Popis

udává zvětšení“ fence oproti bagu, viz. v textu výše ” pokud je TRUE tak se chybějící pozorování ignorují pokud je FALSE jsou hodnoty chybějících pozorování nahrazeny hodnotou mediánu

Bagplot nám ukazuje, kde jsou pozorování dvou proměnných soustředěna. Ukazuje jak míry polohy (medián), tak variabilitu (rozměr bagu) a šikmost (tvar bagu a fence). Stejně jako XY graf pak zobrazuje korelaci (orientace, nebo směr“ bagu). Oproti XY grafu navíc jednoznačně ” zobrazuje vybočující pozorování a dá se z něj rozeznat, zda-li je rozdělení tzv. heavy-tailed (pokud je mnoho vybočujících pozorování na hranici fence). V tabulce 4.2.1 shrnuji užitečná nastavení pro tvorbu bagplotu v R.

4.3 Vícenásobný krabičkový graf Krabičkový graf představený v (3.1) lze velmi výhodně použít i pro zkoumání dvourozměrných dat. Představme si například, že máme data rozdělená do několika skupin. V takovém případě 29 Záměrně

volím stejný datový soubor jako v případě XY grafu.

25


4

DVOUROZMĚRNÁ DATA

5

10

15

20

25

se nabízí provést analýzu rozptylu a zkoumat, zda-li se tyto skupiny z hlediska polohy liší. Stejně tak to ale můžeme zjistit pomocí krabičkových grafů30 . Podívejme se tedy na příklad na obrázek 4.3.1 na kterém jsem vykreslil krabičkové grafy spotřeby automobilů pro jednotlivé počty válců motoru pro data Spotřeba automobilů 2012.

3

4

5

6

8

10

12

16

> > > > + > > > + + + + > + + + > + + +

data = read.csv('data/car_fuel.csv') uq = unique(data$Cylinders) uq = sort(uq[!is.na(uq)]) cc.lph = data$LPH[complete.cases(data$LPH)] # prazdny graf plot.new() plot.window( xlim = c(2, max(uq)+1), ylim = c(0.9*min(cc.lph), 1.1*max(cc.lph)) ); axis(2) boxplot( LPH ~ Cylinders , data = data, at = uq, add = TRUE ) abline( lm(LPH ~ Cylinders , data = data), col = 'red' )

Obrázek 4.3.1: Krabičkové grafy pro spotřebu automobilů pro jednotlivé počty válců motoru Na obrázku 4.3.1 můžeme naprosto zřetelně pozorovat přímou závislost spotřeby na počtu válců. Podíváme-li se na variační rozpětí jednotlivých krabičkových grafů, vidíme, že zejména pro počty válců 4 a 8 je toto rozpětí velké. To by v závislosti na počtu pozorování v této kategorii ovlivnilo případnou regresní analýzu. Pro ilustraci jsem přidal regresní přímku získanou metodou nejmenších čtverců. Je vidět, že nízká spotřeba u počtu válců rovnému 4 a vysoká spotřeba u počtu válců rovnému 8 zvětšují směrnici regresní přímky. V obrázku však stále schází informace o tom, jak jsou jednotlivé skupiny četné. To má velký vliv na platnost našich úsudků o odlišnosti skupin, pro validaci modelu přímky obvzlášť. U krabičkového grafu v R můžeme použít atribut varwidth představený v tabulce 3.1.1, díky kterému se krabičkové grafy vykreslí šířkou úměrnou relativní četnosti dané skupiny. Výsledek na obrázku 4.3.2 je již o mnoho lepší. Porovnávání skupin pomocí krabičkového grafu je vhodné použít i při práci s časovými řadami. 30 Nemusíme

navíc zkoumat platnost předpokladů analýzy rozptylu.

26

4

DVOUROZMĚRNÁ DATA

20


5

10

15

> data = read.csv('data/car_fuel.csv') > > boxplot( + LPH ~ Cylinders , + data = data, + varwidth = TRUE + )

3

4

5

6

8

10

12

16

Obrázek 4.3.2: Krabičkové grafy pro spotřebu automobilů pro jednotlivé počty válců motoru šíroké dle relativní četnosti skupin Vezměme si například čtvrtletní časovou řadu Návštěvnost ČR. Hodnoty časové řady můžeme rozdělit do 4 skupin = 4 čtvrtletí a zkoumat, zda se v datech projevuje sezónnost. Udělejme tak na obrázku 4.3.3. Nejen, že z obrázku vidíme, že se v datech sezónnost vyskytuje, ale vidíme hned i v jakém měřitku. Pokud proměnnou jako v úkazce vycentrujeme odečtením průměru, získáme dokonce i velmi čitelné měřítko pro porovnávání jednotlivých čtvrtletí. Jak jsme se přesvědčili, je krabičkový graf velmi užitečným nástrojem i při práci s dvourozměnými daty.

27

4

1500000

4.4 Sectioned density plot

−1000000 −500000

0

500000

1000000

> > > + + + + + + >

1

2

3

DVOUROZMĚRNÁ DATA

load('data/visits.RData') # prevedme casovou radu na array ar = tapply( visits, cycle(visits), function(x){ x - mean(visits) } ) boxplot(ar)

4

Obrázek 4.3.3: Zkoumání sezónnosti časové řady pomocí krabičkového grafu

4.4 Sectioned density plot Alternativou k violin plotu představeném v (3.5) je sectioned density plot vytvořený Cohenem a Cohenem [3]. Je to vlastně pokus o vtěsnání histogramu, popř. jádrového odhadu na co nejmenší prostor bez ztráty informace. Ačkoliv se jedná o graf znázorňující jednu proměnnou, rozhodl jsem se ho zařadit až do této sekce, protože je určen pro porovnání dvou a více proměnných. Pro vysvětlení konstrukce grafu jsem napsal funkci31 , která na čtyřech grafech na obrázku 4.4.1 ukazuje konstrukci grafu, obdobně jako to udělal autor [3]. Graf vychází z jádrového odhadu hustoty. V prvním kroku se grafem hustoty v hodnotě n kvantilů vykreslí n různobarevných (obvykle se užívají odstíny šedé) horizontálních úseček, spojující body křivky rozdělení. Tyto úsečky se pak sloučí do jedné, která se stane základem stejnobarevných obdélníků. Obdélníky se poté postupně naposouvají, jako na obrázku 4.4.1. Section density plot nabízí vlastně jakýsi pohled shora“, jako kdybychom koukali z výšky na ” dům. V R je v balíčku denstrip [14] obsažena funkce sectioned.density() , která umožňuje snadnou tvorbu tohoto typu grafu. Ukažme si jej pro ilustraci s histogramem na obrázku 4.4.2. Po otočení (tak se používá nejčastěji) je sectioned density plot velice kompaktní. Přímo se pak nabízí k porovnání skupin či výběrů. Ilustrujme jeho použití porovnáním rozdělení šesti výběrů z dat Spotřeba automobilů 2012 o různém rozsahu. Výsledek vidíme na obrázku 4.4.3. 31 Funkce

je k nahlédnutí v (7.3).

28

4

DVOUROZMĚRNÁ DATA

0.0

0.3

0.6

4.4 Sectioned density plot

−3

−2

−1

0

1

2

3

−3

−2

−1

0

1

2

3

−3

−2

−1

0

1

2

3

−3

−2

−1

0

1

2

3

> load('R/plot.explSecDens.R') > plot.explSecDens(n = 6)

Obrázek 4.4.1: Tvorba sectioned density plotu

30

35

Histogram of data

20 15 5

10

set.seed(1) data = rnorm(150) library('denstrip') hist(data) sectioned.density(data, at = 0.1)

0

Frequency

25

> > > > >

−2

−1

0

1

2

data

Obrázek 4.4.2: Sectioned density plot s histogramem na pozadí

29

4.5 Histogramy

n = 50

n = 100

n = 300

n = 500

n = 1108

15 5

10

Spotreba

20

n = 30

4

0

1

2

3 Index

4

5

6

> > > > > + > + > > + > > + + + > > + + > >

DVOUROZMĚRNÁ DATA

set.seed(2) library('denstrip') # priprava dat data = read.csv('data/car_fuel.csv') ccd = data$LPH[complete.cases(data$LPH)] sizes = c(30, 50, 100, 300, 500, length(ccd)) d = list() for (s in 1:6) d[[s]] = sample(ccd, sizes[s]) # prazdny graf plot(NA, ylab = 'Spotreba', xlim = c(0, 6), ylim = c(min(ccd), max(ccd)) ) # grafy for (s in 1:6) sectioned.density(d[[s]], at = s, horiz = FALSE, offset = 0.07 ) labels = paste('n = ', sizes) text(1:6 - 0.5, 22, labels)

Obrázek 4.4.3: Porovnání rozdělení výběrů o různém rozsahu ze stejné populace Sectioned density plot je především nástrojem k porovnání většího počtu skupin. Umožňuje zobrazit veškeré informace z histogramu, respektive jádrového odhadu (viz. sekce (3.2), respektive (3.3)) na velmi malé ploše. Při porovnání různých skupin tedy nepřicházíme o množství informací (jako například bimodalita rozdělení) jako při porovnání použitím krabičkových grafů.

4.5 Histogramy Pro porovnání dvou rozdělení je vhodné použít dva back-to-back“ histogramy. Jak jsme se ” přesvědčili v (3.2), je histogram při správně zvolené délce intervalu velmi informativní a pro porovnání dvou rozdělení naprosto dostačuje. Se dvěma histogramy se velmi často setkáváme v demografii - konkrétně jde o četnosti ženské a mužské populace, ze kterých se vytváří tzv. strom života. Sestavme si jej z dat Střední stav populace ČR do obrázku 4.5.1. Pro tvorbu tohoto typu grafu jsem vytvořil funkci plot.lifeTree , kterou je možné prohlédnout v (7.2). Tím, že máme histogramy hned u sebe, můžeme pozorovat rozdíly v průběhu rozdělení. Například u demografického stromu z obrázku 4.5.1 je patrné, že se s věkem mění poměr počtu

30

4.5 Histogramy

4

DVOUROZMĚRNÁ DATA

100

žijících mužů a žen. Je více mladých mužů a naopak více starých žen.

muzi zeny

load('R/plot.lifeTree.R') data = read.table( 'data/cz_pop.data', header = TRUE ) plot.lifeTree( data$muzi , data$zeny , data$vek )

0

20

40

Vek

60

80

> > + + + > + + + +

−0.010

−0.005

0.000

0.005

0.010

Podil na populaci

Obrázek 4.5.1: Ukázka demografického stromu života

31

5

VÍCEROZMĚRNÁ DATA

5 Vícerozměrná data Při práci s troj- a vícerozměrnými daty nastává při grafickém zobrazování problém. V R sice umíme vykreslit 3D graf, jak uvidíme dále, ale pro více rozměrů už řešení není tak snadné. V zásadě se pak nabízejí dvě možnosti, obě hojně užívané: a) vykreslení matice 2D grafů pro jednotlivé dvojice, b) redukce rozměrů a následné vykreslení do 2D. První postup se hodí pro relativně nízký počet proměnných. Pokud máme více než řekněme 5 proměnných, stává se matice nepřehlednou a grafy v ní příliš malými. Pak už nám nezbývá nic jiného, než přistoupit k vypuštění rozměrů s nejméně informacemi. U troj- a vícerozměrných dat jsou možnosti zobrazení podobné jako u dvourozměrných dat. Kromě základního popisu jednotlivých proměnných se zaměřujeme především na vztahy mezi nimi. Pro výše uvedený problém však není snadné odhalit vybočující pozorování pro více než trojrozměrná data. První z předvedených grafů bude rainbow plot, kterým prozkoumáme průběh dvourozměrných dat v čase. Pro libovolná data ukážeme interaktivní 3D graf. Pro vyjádření četností kategoriálních proměnných použijeme mozaikový graf. Na závěr sekce použijeme R pro grafické znázornění korelační matice.

5.1 Rainbow plot Při práci s rozsáhlými datovými maticemi se často musejí hledat nové metody pro jednoduchou vizualizaci co největšího množství informací. Jednou z novějších metod je tzv. rainbow plot32 představený Hyndmanem a Shangem [13]. Data jsou rozdělena a relevantně seřazena (například časová řada je rozdělena na jednotlivé roky a ty jsou seřazeny) a postupně vykreslována do jednoho grafu jako jednoduché XY grafy různých barev, které jsou postupně vybírány z barev duhy. Než budu pokračovat, ukažme si příklad na panelových datech Míry úmrtnosti v ČR zachycujících úmrtnosti dle věku mezi lety 1950 a 2011. Při použití balíčku demography [12] vykreslíme rainbow plot do obrázku 5.1.1. Vidíme, že křivky pro jednotlivé roky nabývají postupně červené, žluté, zelené, modré a nakonec fialové barvy. Můžeme tedy pozorovat klesající trend i bez toho, že bychom se jej snažili 32 Šlo

by přeložit jako duhový graf, ale nepoužívá se.

32

5.2 Interaktivní 3D graf

5


jakkoliv modelovat.

0

2

Czech Republic: female death rates (1950−2011)

−4 −6 −10

−8

Log death rate

−2

> library('demography') > load('data/mort.RData') > # objekt typu demogdata > class(mort) [1] "demogdata" > plot(mort)

0

20

40

60

80

100

Age

Obrázek 5.1.1: Rainbow plot pro míry úmrtnosti v Českých zemích Rainbow plot lze využít i pro obyčejné časové řady jedné proměnné. V R jsem pro něj vytvořil vlastní funkci plot.rainbow() 33 , jelikož balíček demography umí pracovat jen se specifickým formátem dat. Podívejme se na rainbow plot z dat Produkce plynu v Austrálii na obrázku 5.1.2, kde časovou řadu dělím podle let. Problémem rainbow plotu je, že nemáme příliš představu o tom, kde leží mediánová“ křivka ” a nejsme schopni identifikovat odlehlá pozorování. Řešením může být functional bagplot, který je stejně jako rainbow plot navrhován Hyndmanem a Shangem [13]. Základem tohoto grafu je jako u bagplotu bag, který obsahuje prostředních“ 50% křivek34 . Tvorba takového ” grafu je však nad rámec této práce.

5.2 Interaktivní 3D graf Pokud máme trojrozměrná data, nebo se na tři rozměry zvládneme statistickými metodami dostat, můžeme využít interaktivního 3D grafu. Funkce plot3d() obsažená v balíčku rgl [2] vykreslí krychli, se kterou můžeme otáčet a prohlížet data z jakéhokoliv úhlu. Graf je také možné přibližovat a oddalovat.

33 Funkci 34 Více

je možné prozkoumat v sekci (7.1). o bagplotu v sekci (4.2).

33

5.2 Interaktivní 3D graf

5


load('R/plot.rainbow.R') library('forecast') data(gas) plot.rainbow(gas)

10000

20000

30000

40000

> > > >

0

bounds$y

50000

60000

1956 1969 1982 1995

2

4

6

8

10

12

bounds$x

Obrázek 5.1.2: Rainbow plot poukazující na přítomnost trendové složky v časové řadě

> data = read.csv('data/car_fuel.csv') > library('rgl') > plot3d(data[,c(4,5,8)])

Obrázek 5.2.1: Interaktivní 3D graf

34

5.3 Mozaikový graf

5


Dal by se označit jako XYZ graf, jelikož nám poskytuje stejné informace jako XY graf pro dva rozměry. Umožňuje nám pozorovat směr závislosti a celkovou variabilitu. Tento graf se také hodí při použití metody nejmenších čtverců, kdy můžeme jednoduše identifikovat odlehlá pozorování, která vybočují z námi vypočítaného modelu roviny.

5.3 Mozaikový graf Předpokládejme, že jsme do rukou dostali data sestávající z kategoriálních proměnných. Zpracovali jsme je do formy kontingenční tabulky. Pro více proměnných se kontingenční tabulka stává poněkud nepřehledná. Mozaikový graf nám může pomoci. Je to zobrazení kontingenční tabulky pomocí dlaždic, které jsou vykreslovány proporčně k hodnotám v tabulce [9]. Šířka dlaždic je vykreslována proporčně k marginální četnosti n+ j . Výška je vykreslována podle podmíněných četností ni| j jednotlivých řádků. Obsah jednotlivých dlaždic pak znázorňuje četnosti v jednotlivých buňkách tabulky. Mozaikový graf nám umožňuje rozpoznat četnosti napříč kategoriemi na první pohled. V R budeme pracovat s balíčkem vcd [16] a to na ukázkových trojrozměrných datech Barva očí a vlasů. Ukázku výchozího (původního) mozaikového grafu vidíme na obrázku 5.3.1. Eye Hazel Green Female Male

Blue

> library('vcd') > mosaic(HairEyeColor)

Female Male Female Male

Blond

Red

Female Sex

Hair

Brown

Male

Black

Brown

Obrázek 5.3.1: Mozaikový graf ve výchozím nastavení Porovnáváním rozměrů jednotlivých dlaždic můžeme zkoumat závislost kategorií. V případě nezávislosti jsou všechny dlaždice na řádků stejně vysoké [9]. Pro zřetelnější výstup se často používá šrafování či výplně, která znázorňuje velikost odchylky od hodnot v případě nezávislosti. Pro vybarvení dlaždic přidáme parametr shade jako na obrázku 5.3.2.

35

5.4 Vizualizace korelační matice

5


Eye Hazel Green Female Male

Blue

Pearson residuals: 8.02

Male

Black

Brown

Female Sex

Hair

Brown

4.00

2.00

Female Male Female Male

Blond

Red

0.00

> library('vcd') > mosaic( + HairEyeColor , + shade = TRUE + )

−2.00

−4.19 p−value = < 2.22e−16

Obrázek 5.3.2: Mozaikový graf s barevným znázorněním velikosti reziduí Mozaikový graf je vhodné použít pro získání rychlého přehledu o zastoupení kategorií, včemž předčí samotnou kontingenční tabulku. Je to jakési zobecnění sloupcového diagramu do více rozměrů.

5.4 Vizualizace korelační matice Pokud máme velké množství spojitých kvantitativních proměnných a chceme zkoumat, jak se vzájemně ovlivňují, konstruujeme korelační matici. Problémem (stejně jako u kontingenční tabulky v (5.3)) je však nepřehlednost při větším počtu proměnných. Než abychom se snažili vyznat v rozsáhlé tabulce, použijeme raději grafické řešení. V R můžeme v balíčku ellipse [18] najít funkci plotcorr() , která vykreslí korelační matici pomocí elips udávajících směr a sílu závislosti. Použijme ji na datech Údaje o automobilech na obrázku 5.4.1. Kromě tvaru elipsy můžeme využít také barev. Chladné barvy (odstíny modré) by vyjadřovaly přímou závislost a barvy teplé (odstíny červené) závislost nepřímou. Sytost barvy by pak udávala intezitu závislosti. Ilustrujme takové nastavení na obrázku 5.4.2, s drobnou úpravou převzatého od Murdocha s Chowem [17]. S obrázkem 5.4.2 už můžeme být spokojeni. Na první pohled je vidět, která proměnná závisí na které, jaký je směr a síla závislosti.

36

5


carb

gear

qsec

wt

drat

hp

disp

cyl

mpg

5.4 Vizualizace korelační matice

mpg

> > + > >

cyl disp hp drat

library('ellipse') data = mtcars[, !(names(mtcars) %in% c('vs','am'))] # cex.lab pro vetsi pismo plotcorr(cor(data), cex.lab = 1.5)

wt qsec gear carb

wt cyl disp hp carb qsec gear drat mpg

mpg

drat

gear

qsec

carb

hp

disp

cyl

wt

Obrázek 5.4.1: Vizualizace korelační matice pomocí funkce

> > + > > > > > > + + + > > + +

plotcorr()

library('ellipse') data = mtcars[, !(names(mtcars) %in% c('vs','am'))] # ziskame korelacni matici corr.mtcars = cor(data) ord = order(corr.mtcars[1,]) xc = corr.mtcars[ord, ord] # odstiny modre a cervene colors = c("#A50F15","#DE2D26", "#FB6A4A","#FCAE91","#FEE5D9", "white", "#EFF3FF","#BDD7E7", "#6BAED6","#3182BD","#08519C") # vykreslime graf plotcorr(xc, col=colors[5*xc + 6], cex.lab = 1.5)

Obrázek 5.4.2: Vizualizace korelační matice pomocí funkce s využitím barev

37

6 UNIVERZÁLNÍ NASTAVENÍ A ÚPRAVY GRAFŮ

6 Univerzální nastavení a úpravy grafů 6.1 Kombinování grafů Pokud chceme grafy kombinovat, máme hned dvě možnosti: a) Vykreslíme průběh všech funkcí do jednoho grafu v jednom souřadnicovém systému. b) Vykreslíme jednotlivé grafy zvlášť, ale do jednoho obrázku. V prvním případě si vystačíme s parametrem add = TRUE u funkcí vytvořujících graf. V druhém případě musíme nastavit parametry grafu, jak popíši dále. 6.1.1 Více grafů v jednom obrázku Pro základní rozdělení obrázku na několik částí nám postačí nastavit parametr mfrow pomocí funkce par() 35 . Hodnotou parametru musí být dvoučlenný vektor, určující počet řádků a sloupců. Použitím můžeme například vytvořit čtveřici grafů a umístit je do mřížky jako na obrázku 6.1.1. Histogram of runif(50)

6 4

Frequency

8

10 8 6

0

0

2

2

4

Frequency

10

14

Histogram of rnorm(50)

−1

0

1

2

0.0

0.2

0.4

0.6

0.8

rnorm(50)

runif(50)

Histogram of rpois(50, 1.5)

Histogram of rt(50, 10)

1.0

# nastavime mrizku par(mfrow = c(2,2)) # vykreslime grafy hist(rnorm(50)) hist(runif(50)) hist(rpois(50, 1.5)) hist(rt(50, 10))

0

0

5

10

Frequency

10 5

Frequency

15

15

−2

> > > > > > >

0

1

2

3

rpois(50, 1.5)

4

−4

−3

−2

−1

0

1

2

3

rt(50, 10)

Obrázek 6.1.1: Rozdělení grafu do mřížky pomocí

mfrow

Pokud si nevystačíme s pravidelným rozložením grafů, musíme použít funkci layout() . Ta nám umožňuje nastavit graf úplně přesně. Atributem je matice, která určuje, kde se který 35 Parametr mfrow

je jedním z globálních nastavení grafu, podobně jako například mar .

38

6.2 Přidání bodů do grafu


graf vykreslí. Grafům přidelíme čísla, která použijeme při vyplnění těch řádků a sloupců, které chceme, aby zaplnil. Příklad vidíme na obrázku 6.1.2. Dalším příkladem může být obrázek 3.1.3 na straně 8. Histogram of rnorm(50)

10 8 6

Frequency

4

6 4

0

0

2

2

Frequency

8

12

Histogram of runif(50)

0.0

0.2

0.4

0.6

0.8

1.0

−2

−1

0

1

2

rnorm(50)

0.0 −1.0

−0.5

sin(x)

0.5

1.0

runif(50)

> # nastavime mrizku > (mat = matrix(c(1,3,2,3), ncol = 2)) [,1] [,2] [1,] 1 2 [2,] 3 3 > layout(mat) > # vykreslime grafy > hist(runif(50)) > hist(rnorm(50)) > curve(sin(x), from = 0, to = 4*pi)

0

2

4

6

8

10

12

x

Obrázek 6.1.2: Rozdělení grafu pomocí

layout()

6.2 Přidání bodů do grafu Pro přidání bodů do již vykresleného grafu můžeme použít funkci points() . Specifikujeme umístění bodů a je dobré také zvolit symbol (a případně barvu), kterým se tento bod vykreslí. Přidejme do obrázku 6.2.1 několik bodů s různými symboly a barvami. Seznam symbolů 6.2.2 můžeme snadno vygenerovat.

6.3 Přidání čar do grafu 6.3.1 Přidání přímky Pro přidání přímky je možno použít funkci abline() . Určime úrovňovou konstantu a směrnici přímky. Do obrázku 6.3.1 jsem vykreslil hned několik přímek. Funkce abline() je také možno použít při modelování regresní přímkou. Pokud do argumentu funkce vložíme výstup z lm() , vykreslí se nám regresní přímka.

39

6.3 Přidání čar do grafu


0 −2

−1

rnorm(100)

1

2

> > > > > > + + > > + + + + + 0

20

40

60

80

plot(rnorm(100), col = 'gray') souradnice_x = c(40, 60, 90) souradnice_y = c(2, -2, 1) kod_symbolu = c(2, 5, 8) barva_symbolu = c( 'red', 'blue', 'green' ) points( souradnice_x , souradnice_y , pch = kod_symbolu , col = barva_symbolu )

100

Index

Obrázek 6.2.1: Přidání bodu do grafu pomocí

0

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

> > > > + > > > > + + + + >

seq x = y = for y x = y =

points()

= 0:18 rep((1:4)*20, 5) c() (i in rev(1:5)) = c(y, rep(i, 4)) x[seq + 1] y[seq + 1]

plot(x, y, pch = seq, xaxt = 'n', yaxt = 'n', main = NULL, frame = FALSE, xlim = c(10, 90), ylim = c(0, 5.5) ) text(x, y - 0.45, seq, cex = 0.7)

Obrázek 6.2.2: Seznam základních znaků pro vykreslení bodu a jejich číselných protějšků

40



y

−40

−20

0

20

40

> # vytvoreni prazdneho grafu > plot(c(-50, 50), c(-50, 50), type = → → 'n', xlab = 'x', ylab = 'y') > > # y = 35 + 0.5x > abline(35, 0.5) > # y = 13 - 0.75x > abline(13, -0.75, col = 'red') > # y = 43 > abline(43, 0, col = 'orange')

−40

−20

0

20

40

x

Obrázek 6.3.1: Přidání přímek pomocí

abline()

6.3.2 Přidání grafu dle předpisu funkce Do grafu lze přidat i křivku dle předpisu. Pomocí funkce curve() můžeme přidat funkci libovolného zápisu. Do obrázku 6.3.2 jsem jich přidal hned několik. 6.3.3 Přidání speciﬁckých čar Pro spojení dvou bodů je nejsnazší použít univerzální funkci lines() . Jedná se o obdobu funkce points() . Rozdíl je pouze v tom, že se body neukazují samostatně, nýbrž jsou propojeny čárou. Ukázkou propojení hned několika dvojic bodů vidíme na obrázku 6.3.3. Čáru je možné si libovolně přizpůsobit, jak ukazuje tabulka (6.3.1) s parametry36 . Tabulka 6.3.1: Specifické parametry pro funkci lines() Parametr

Výchozí hodnota

Popis

lty

1

zkrácenina pro line-type, tedy druh čáry; základních druhů je 6

lwd

1

nastavení tloušťky čáry

Pro ilustraci možných druhů čar ( lty ) si můžeme vygenerovat seznam do obrázku 6.3.4. 36 Tato

nastavení je možno použít kdykoliv se v grafech setkáme s čarami.

41



−2

−1 −2

−1

0

1

# prazdny graf plot( x = c(-2, 2), y = c(-2, 2), type = 'n', xlab = 'x', ylab = 'y' ) # y = x^3 curve(x^3, add = TRUE) # y = cos(3*x) curve(cos(3*x), col = 'green', add = TRUE) # y ~ N(0,1) curve(dnorm(x), col = 'red', add = TRUE)

2

x

Obrázek 6.3.2: Přidání grafů funkcí do grafu

50

100

150

> + > + + + + + + + + + + >

0

y

0

1

2

> > + + + + + > > > > + > > +

0

20

40

60

plot(c(0, 90), c(0, 170), type = 'n', xlab = '', ylab = '') xy = rbind( c(20, 20), c(70, 20), c(70, 90), c(45, 150), c(20, 90), c(20, 20), c(70, 90), c(20, 90), c(70, 20) ) lines(xy)

80

Obrázek 6.3.3: Propojení bodů pomocí funkce

42

lines()


1 2

y

3 4 5


> > > > > > > + + > > + > > +

seq = rev(1:6) mar = 20 y = cbind(seq*mar, seq*mar) x = cbind(seq*0, seq*0 + 2) plot(x, y, type = 'n', main = NULL, xaxt = 'n', yaxt = 'n', ylim = c(0, 130), frame = FALSE) for (i in seq) lines(x[i,], y[i,], lty = i, lwd = 2) text(x = mean(x[1,]), y = y - 8, labels = rev(seq))

6

Obrázek 6.3.4: Seznam možných druhů čar

43

7

VYTVOŘENÉ FUNKCE

7 Vytvořené funkce V této sekci pro úplnost vypisuji mnou vytvořené funkce.

7.1 Rainbow plot funkce Prototyp funkce pro tvorbu rainbow plotu37 z časové řady. Ukázkou je obrázek 5.1.2. > plot.rainbow = function(data, legend.count = 4, legend.pos = 'topright') { + # pomocne promenne + data.years = unique(floor(time(data))) + data.freq = frequency(data) + # rozdelime casovou radu dle let + i = seq(data) + y = split(data, ceiling(i/data.freq)) + y.len = length(y) + # prazdny graf + bounds = list( + x = c(1, data.freq), + y = c(min(unlist(y)), max(unlist(y))) + ) + plot(bounds , type = 'n') + # barvy duhy + cols = rainbow(y.len) + # vykresleni car + r = lapply(1:y.len, function(i) { + y_ = as.numeric(unlist(y[i])) + lines(1:length(y_), y_, col = cols[i]) + }) + ind = seq(1956, 1995, len = legend.count) + # legenda + legend( + x = legend.pos, + lty = c(1,1), + legend = ind, + col = rainbow(length(ind)) + ) + }

37 Více

o tomto grafu v (5.1).

44

7.2 Funkce pro tvorbu stromu života

7

VYTVOŘENÉ FUNKCE

7.2 Funkce pro tvorbu stromu života Prototyp funkce pro tvorbu demografického stromu života. Ukázka na obrázku 4.5.1. > plot.lifeTree = function(a, b, names = c()) { + # jen stejna struktura dat + if (length(a) != length(b)) + return(FALSE) + # relativni cisla + sum.ab = sum(c(a,b)) + a = a/sum.ab + b = b/sum.ab + # osa y + if (!length(names)) + names = 1:length(a) + # pomocne promenne + data.min = min(c(a,b)) + data.max = max(c(a,b)) + data.length = length(a) + data.seq = 1:length(a) + data.col = c('#4682AB', '#AB4650') + # okraje + # par(mar = c(2,3,2,0)) + # prazdny graf + plot( + NA, + xlim = c(-data.max, data.max), + ylim = c(names[1], names[length(names)]), + xlab = 'Podil na populaci', + ylab = 'Vek' + ) + # cary + rect(-a, data.seq, 0, data.seq + 1.5, col = data.col[1], border = 'NA') + rect(b, data.seq, 0, data.seq + 1.5, col = data.col[2], border = 'NA') + # legenda + legend('topright', + col = data.col, + legend = c('muzi', 'zeny'), + lwd = 3 + ) + }

45

7.3 Funkce popisující vznik sectioned density plotu

7

VYTVOŘENÉ FUNKCE

7.3 Funkce popisující vznik sectioned density plotu > plot.explSecDens = function(n = 6) { + # n = pocet car + newPlot = function() { + plot.new(); plot.window(xlim = c(-3,3), ylim = c(0,0.6)); axis(1) + } + shades = colorRampPalette(c('black', 'white')) + par(mfrow = c(4,1), lend = 2, mar = c(3,3,1,3)) + # graf rozdeleni + newPlot() + axis(2) + d = curve(dnorm(x), add = TRUE) + # vykresleni n car z kvantilu + q = quantile(d$y, probs = (1:n)*(0.85/n)) + coor = lapply(q, FUN = function(x) { + a = round(abs(d$y - x), 8) + ind = which(round(a, 8) == min(a)) + list(x = d$x[ind], y = d$y[ind]) + }) + colors = shades(n) + newPlot() + lns = lapply(coor, FUN = function(x) { + lines(x, col = colors , lwd = 3) + colors <<- colors[-1] + }) + # vykresleni obdelniku + coor = lapply(coor, FUN = function(x) { + x$y = c(0.3, 0.3) + x + }) + colors = shades(n) + newPlot() + lns = lapply(coor, FUN = function(x) { + lines(x, col = colors , lwd = 30) + colors <<- colors[-1] + }) + # posun -> finalni graf + colors = shades(n) + newPlot() + times = 0 + lns = lapply(coor, FUN = function(x) { + x$y = c(0.15, 0.15) + 0.07*times + lines(x, col = colors , lwd = 30) + times <<- times + 1 + colors <<- colors[-1] + }) + + }

46

8

ZÁVĚR

8 Závěr Cílem práce bylo představit možnosti grafického zobrazení dat a ukázat jejich implementaci v R. Grafická zobrazení byla roztříděna dle počtu proměnných, který jsou schopna zachytit. V každé kategorii byla podle charakteru proměnné postupně vybrána ta nejužitečnější zobrazení, která byla důsledně popsána. Ukázal jsem, pro co se jednotlivé typy zobrazení hodí, v čem vynikají a v čem naopak ztrácejí. Poukázal jsem na případné problémy použitelnosti a pokud to šlo, navrhl jsem řešení, jak tyto problémy řešit. Pro všechna vybraná zobrazení jsem ukázal implementaci v R. V případě, že v R nebyla implementace přímočará, vytvořil jsem vlastní funkce, které tvorbu takových zobrazení umožňují. Konkrétně jsem vytvořil funkce pro tvorbu demografického stromu života, rainbow plotu a funkci ilustrující tvorbu shaded density plotu. Ukázal jsem, že R je velmi silný nástroj, nejen co se grafického znázorňování dat týče. Kromě samotných grafických zobrazení a jejich významu jsem se zabýval také obecným přizpůsobením grafů v R. Po přečtení práce je i méně zkušený čtenář schopen zvolit vhodný typ grafu na základě předložených dat a zvládne jej v R realizovat. Vyčte z grafu veškeré dostupné informace a vyvodí z nich závěry. Je schopen graf dále přizpůsobovat a jednotlivé grafy kombinovat dohromady. Ví, na co si při práci s jednotlivými typy grafů musí dát pozor. V práci se povedlo ukázat, jak je grafické zobrazení dat důležité. Je vhodné, aby grafické zobrazení dat předcházelo každé číselné analýze. Jak jsme se několikrát přesvědčili, může číselnou analýzu dokonce i nahradit, popřípadě ukázat na nevhodnost jejího použití.

47

9 REFERENCE

9 Reference [1] Daniel Adler. vioplot: Violin plot, 2005. R package version 0.2. [2] Daniel Adler and Duncan Murdoch. rgl: 3D visualization device system (OpenGL), 2011. R package version 0.92.798. [3] Dale J. Cohen and Jon Cohen. The sectioned density plot. The American Statistician, 60:167–174, 2006. [4] P. Dalgaard. Introductory Statistics with R. Statistics and Computing. Springer, 2008. [5] A. Elgammal, R. Durai swami, D. Harwood, and L.S. Davis. Background and foreground modeling using nonparametric kernel density estimation for visual surveillance. Proceedings of the IEEE, 90(7):1151–1163, 2002. [6] Warren W. Esty and Jeff Banfield. The box-percentile plot. Journal of Statistical Software, 8(17):1–14, 10 2003. [7] Ian Fellows. wordcloud: Word Clouds, 2012. R package version 2.2. [8] J. Fox and S. Weisberg. An R Companion to Applied Regression. SAGE Publications, 2010. [9] Michael Friendly. Mosaic displays for multi-way contingency tables. Journal of the American Statistical Association, 89(425):pp. 190–200, 1994. [10] Jerry L. Hintze and Ray D. Nelson. Violin Plots: A Box Plot-Density Trace Synergism. The American Statistician, 52(2):181–184, 1998. [11] Rob J Hyndman, contributions from George Athanasopoulos, Slava Razbash, Drew Schmidt, Zhenyu Zhou, and Yousaf Khan. forecast: Forecasting functions for time series and linear models, 2013. R package version 4.03. [12] Rob J Hyndman, contributions from Heather Booth, Leonie Tickle, John Maindonald, Simon Wood, and R Core Team. demography: Forecasting mortality, fertility, migration and population data, 2012. R package version 1.14. [13] Rob J Hyndman and Han Lin Shang. Rainbow plots, bagplots and boxplots for functional data. Monash Econometrics and Business Statistics Working Papers 9/08, Monash University, Department of Econometrics and Business Statistics, 2008. [14] C. H. Jackson. Displaying uncertainty with shading. 62(4):340–347, 2008.

48

The American Statistician,

9 REFERENCE [15] John I. Marden. Positions and QQ Plots. Statistical Science, 19(4):606–614, 2004. [16] David Meyer, Achim Zeileis, and Kurt Hornik. The strucplot framework: Visualizing multi-way contingency tables with vcd. Journal of Statistical Software, 17(3):1–48, 10 2006. [17] D. J. Murdoch and E. D. Chow. A graphical display of large correlation matrices. The American Statistician, 50(2):178–180, 1996. [18] Duncan Murdoch and E. D. Chow. ellipse: Functions for drawing ellipses and ellipse-like confidence regions, 2012. R package version 0.3-7. [19] Ida Ruts & John W. Tukey Peter J. Rousseeuw. The Bagplot: A Bivariate Boxplot. The American Statistician, 1999. [20] R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, 2013. ISBN 3-900051-07-0. [21] r project.org. About. http://www.r-project.org/about.html, 2013. [Online; zobrazeno 3. března 2013]. [22] Martin ŘEZÁČ. Jádrové odhady hustoty [online]. Disertační práce, Masarykova univerzita, Přírodovědecká fakulta, 2007 [cit. 2013-04-04]. [23] David W. Scott. On optimal and data-based histograms. Biometrika, 66(3):pp. 605–610, 1979. [24] D.W. Scott. Multivariate Density Estimation: Theory, Practice, and Visualization. Wiley series in probability and mathematical statistics: Applied probability and statistics. Wiley, 1992. [25] Simon J. Sheather. Density estimation. Statistical Science, 19(4):pp. 588–597, 2004. [26] C. H. Sim, F. F. Gan, and T. C. Chang. Outlier labeling with boxplot procedures. Journal of the American Statistical Association, 100(470):pp. 642–652, 2005. [27] John W. Tukey. Mathematics and the Picturing of Data. In Ralph D. James, editor, International Congress of Mathematicians 1974, volume 2, pages 523–532, 1974. [28] M. P. Wand. Data-based choice of histogram bin width. The American Statistician, 51:59–64, 1996. [29] David F. Williamson, Robert A. Parker, and Juliette S. Kendrick. The box plot: A simple visual method to interpret data. Annals of Internal Medicine, 110(11):916–921, 1989.

49

9 REFERENCE [30] Peter Wolf and Uni Bielefeld. aplpack: Another Plot PACKage: stem.leaf, bagplot, faces, spin3R, and some slider functions, 2012. R package version 1.2.7.

50

A

PŘÍLOHA

A Příloha A.1 Použitá data A.1.1 Spotřeba automobilů 2012 Data jsou ze serveru fueleconomy.gov, spravovaného vládou USA. Ponechal jsem jen několik sloupců a přidal evropskou míru spotřeby - litry na 100 km38 - převodem dle (6).

L100km = LPH =

235.214 MPG

(6)

Data je možné volně stáhnout na http://www.fueleconomy.gov/feg/epadata/12data.zip. Vzorek dat: > data = read.csv('data/car_fuel.csv') > length(data[,1]) [1] 1142 > head(data) Model.Yr Mfr.Name Carline Engine Cylinders MPG Gears LPH 1 2012 aston martin V12 Vantage 5.9 12 13.1319 6 17.91165 2 2012 aston martin V8 Vantage 4.7 8 15.9819 6 14.71752 3 2012 aston martin V8 Vantage 4.7 8 16.0421 7 14.66229 4 2012 aston martin V8 Vantage 4.7 8 15.1971 6 15.47756 5 2012 aston martin V8 Vantage S 4.7 8 16.0421 7 14.66229 6 2012 Audi R8 4.2 8 15.9120 6 14.78218

A.1.2 Produkce plynu v Austrálii Měsíční časová řada pro roky 1956 až 1995 zachycující produkci plynu v Austrálii. Data jsou obsažena v balíčku forecast [11]. Data: > library('forecast') > str(gas) Time-Series [1:476] from 1956 to 1996: 1709 1646 1794 1878 2173 ... > frequency(gas) [1] 12

38 LPH

jako litres per hundred.

51

A.1 Použitá data

A

PŘÍLOHA

A.1.3 Návštěvnost ČR Tato čtvrtletní časová řada je z Českého statistického úřadu39 . Přesně se jedná o Návštěvnost v hromadných ubytovacích zařízeních v ČR. Data jsem stáhl jako xls soubor a převedl jej do R. Data: > load('data/visits.RData') > str(visits) Time-Series [1:52] from 2000 to 2013: 2082827 2915857 3903838 1961250 → → 2089561 ... > frequency(visits) [1] 4

A.1.4 Míry úmrtnosti v ČR Roční panelová data udávající míry úmrtnosti v Českých zemích mezi lety 1950 a 2011. Data pocházejí z Human Mortality Database a jsou (po registraci) k dispozici ke stažení40 . Pro jednoduchou práci s daty využívám balíčku demography [12], který velmi snadno zachází s daty strukturovaných dle Human Mortality Database. Data: > library('demography') > load('data/mort.RData') > mort Mortality data for Czech Republic Series: female male total Years: 1950 - 2011 Ages: 0 - 110

39 Lze

ji stáhnout z http://www.czso.cz/.

40 http://www.mortality.org/cgi-bin/hmd/country.php?cntr=CZE&level=1

52

A.1 Použitá data

A

PŘÍLOHA

A.1.5 Střední stav populace ČR Data jsou kombinací dvou xls souborů z Českého statistického úřadu41 , konkrétně Věkové složení mužů k 1.7.2011 a Věkové složení žen k 1.7.2011, udávajících střední stavy mužů a žen za rok 2011. Vzorek dat: > data = read.table('data/cz_pop.data', header = TRUE) > head(data) vek zeny muzi 1 0 57960 61150 2 1 59350 61777 3 2 59893 62882 4 3 57913 60335 5 4 52713 55963 6 5 50122 52726

A.1.6 Barva očí a vlasů Data obsahují jednoduchou křížovou tabulku s proměnnými pohlaví, barva očí a barva vlasů. Data jsou obsažena přímo v R [20] jako HairEyeColor . Data: > HairEyeColor , , Sex = Male Eye Hair Brown Blue Hazel Green Black 32 11 10 3 Brown 53 50 25 15 Red 10 10 7 7 Blond 3 30 5 8 , , Sex = Female Eye Hair Brown Blue Hazel Green Black 36 9 5 2 Brown 66 34 29 14 Red 16 7 7 7 Blond 4 64 5 8

41 Ke

stažení na http://www.czso.cz/csu/2012edicniplan.nsf/p/4003-12.

53

A.2 Grafy

A

PŘÍLOHA

A.1.7 Údaje o automobilech Data obsahují 11 proměnných o 32 automobilech z roku 1974. Jsou obsažena přímo v R [20]. Data: > head(mtcars) Mazda RX4 Mazda RX4 Wag Datsun 710 Hornet 4 Drive Hornet Sportabout Valiant

mpg cyl disp hp drat wt qsec vs am gear carb 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1

A.2 Grafy

0.5

Jadra v R pro sirku vyhlazovaciho okna = 1

0.0

0.1

0.2

Density

0.3

0.4

gaussian epanechnikov rectangular triangular biweight cosine optcosine

−3

−2

−1

0

1

2

3

> kernels = + eval(formals(density.default)$kernel) > plot(density(0, bw = 1), xlab = "", + main = "Jadra v R pro sirku → → vyhlazovaciho okna = 1", + ylim = c(0, 0.5)) > for(i in 2:length(kernels)) + lines( + density(0, bw = 1, + kernel = kernels[i]), + col = i + ) > legend('topright', legend = kernels , + col = seq_along(kernels), lty = 1, + cex = 1.2 + )

Obrázek A.2.1: Vykreslení jádrových funkcí do jednoho souřadnicového systému

54

A

20

A.2 Grafy

10

15

> > > > > > > > >

Hloubka 2

set.seed(1) pts = sample(1:25, 16, replace = TRUE) plot(pts) abline(2, 0, lty = 2, col = 2) text(10, 2 + 1, 'Hloubka 1', col = 2) abline(20, -5, lty = 2, col = 3) text(2, 10 + 1, 'Hloubka 2', col = 3) abline(0.5, 0.5, lty = 2, col = 4) text(11, 6 + 1, 'Hloubka 3', col = 4)

5

Hloubka 3

Hloubka 1

5

10

15

Obrázek A.2.2: Výpočet hloubky ve dvou rozměrech

55

PŘÍLOHA

Vysoká škola ekonomická v Praze

Recommend Documents