Grammar of Graphics (ggplot2)

Knihovna ggplot2: implementace Grammar of Graphics (Leland Wilkinson), http://docs.ggplot2.org/

Krásný „tahák“ s přehledem základních možností: http://www.rstudio.com/resources/cheatsheets/ (Data Visualization Cheat Sheet)

Přehledné grafy od nejjednodušších po pokročilé: https://psyteachr.github.io/msc-data-skills/ggplot.html

Příklady grafů v R: http://www.cookbook-r.com/Graphs/

Grafy a obrázky zadáváme pomocí „podstatných jmen“, „přídavných jmen“ a „sloves“.

Vycházíme z dat v dobře uspořádaném dataframe (formát tidy data). Faktorové proměnné musejí být v dataframe skutečně označeny jako faktorové (viz kapitola 2. Faktorové proměnné), abychom znázornili, že se jedná o diskrétní kategorie, které nám mohou např. rozdělovat data do skupin. Pozor! Použití faktorových proměnných je obzvláště nutné u číselných hodnot, které mají být chápány jako diskrétní kategorie. Pokud máme např. proměnnou PořadíČtení s hodnotami 1 a 2, při jejím použití na x-ovou osu grafů může dojít k vytvoření reálné spojité osy, pokud je proměnná typu numeric. Jestliže ji převedeme na faktorovou proměnnou, budou správně vytvořeny dvě diskrétní hodnoty.

Vše by mělo být v datech hezky pojmenované, protože se ze jmen automaticky stávají popisky grafů. Tedy raději proměnná ‘Pohlaví’ než ‘pohl’, hodnoty spíše ‘muž’ a ‘žena’ než 0 a 1.

Dále přidáváme tzv. vzhled, estetiku (aesthetics), kde určujeme, co z dat se má mapovat na osy x a y, co má určovat velikost, tvar a barvu. Legendy (popisky) se tvoří podle těchto vlastností automaticky.

Pomocí facets je možné rozdělit obrázek do více panelů na základě jedné nebo více diskrétních kategorií (faktorové proměnné). Jedná se tedy opět o jistou formu aesthetics, kdy se automaticky zjistí, kolika různých hodnot daná proměnná nabývá, podle toho se data rozdělí do příslušných skupin a pro každou se vykreslí samostatný a přehledně popsaný graf.

Důležité je určit druh tvaru (geoms), jakým se mají data znázornit. Zda to mají být body, čáry, sloupce, histogramy, odhady hustoty pravděpodobnosti, krabicové grafy apod. V rámci jednoho zobrazení můžeme přidat současně klidně více různých geoms, pokud to dává smysl.

Aby toho nebylo málo, připojit můžeme statistické ukazatele (stats). Může se jednat například o regresi (automatické vykreslení pásu 95% pravděpodobnosti výskytu je samozřejmostí), zobrazení středu dat včetně případného 95% intervalu spolehlivosti, ale také třeba elipsy znázorňující 2D oblast výskytu 95 % nejčastějších hodnot. To vše pomocí jednoduchého příkazu, což ušetří mnoho náročné práce.

Dále je možné přidávat vlastní popisky (labels), měnit souřadnicový systém (coordinate systems) a volit různá témata (themes), protože základní vzhled je sice hezký, ale občas potřebujeme jednodušší zobrazení.

Instalace knihovny ggplot2, která je součástí package tidyverse

Pokud ještě není knihovna tidyverse nainstalována, učiníme tak příkazem (stačí provést jednou)

install.packages("tidyverse")

Před použitím knihovny je potřeba ji nejdříve zavést příkazem (nutné vždy po startu R)

library(tidyverse)

Příklad 1: korelační graf (scatter-plot) několika skupin

Stáhnout soubor s daty: modre_cervene.csv.

tabulka <- read.csv("modre_cervene.csv", encoding = "UTF-8")

ggplot(data = tabulka, aes(x = x, y = y)) +
    geom_point()

ggplot(data = tabulka, aes(x = x, y = y, color = group)) +
    geom_point()

ggplot(data = tabulka, aes(x = x, y = y, shape = group)) +
    geom_point()

ggplot(data = tabulka, aes(x = x, y = y, shape = group, color = group)) +
    geom_point()

ggplot(data = tabulka, aes(x = x, y = y, shape = group, color = group)) +
    geom_point(size = 3)

ggplot(data = tabulka, aes(x = x, y = y, shape = group, color = group)) +
    geom_point(size = 3) +
    stat_smooth(method = "lm")   # lineární regrese vč. 95% pásů

Library plotly

Pomocí knihovny plotly je možné z obrázků ggplot udělat interaktivní grafy s možnostmi kurzoru pro odečítání hodnot, zoomování a rolování.

library(plotly)

g <- ggplot(data = tabulka, aes(x = x, y = y, shape = group, color = group)) +
    geom_point(size = 3) +
    stat_smooth(method = "lm")
ggplotly(g)

Příklad 2: klasický čárový graf pro 2 proměnné - jak se postupně měnil tlak s teplotou.

ggplot(data = pressure, aes(x = temperature, y = pressure)) +
    geom_line()

Příklad 3: dvojice formantů vokálů rozdělených do skupin dle několika faktorů

Příprava dat

Stáhnout soubor s daty: vokalyhz.txt.

data <- read.table("vokalyhz.txt", header = TRUE, sep = ",")
str(data) # Je vidět, že Vowel a Sex jsou faktorové, což je dobré. SpeakerID by ale také měl být faktor.

## 'data.frame':    3000 obs. of  7 variables:
##  $ Vowel    : chr  "a" "a" "a" "a" ...
##  $ Sex      : chr  "f" "f" "f" "f" ...
##  $ F0.mean  : num  212 219 252 244 231 ...
##  $ SpeakerID: int  1 2 3 4 7 8 10 11 12 13 ...
##  $ F1_Hz    : num  819 713 748 791 788 ...
##  $ F2_Hz    : num  1354 1393 1495 1754 1622 ...
##  $ F3_Hz    : num  2595 2634 2497 2443 2961 ...

data$SpeakerID <- factor(data$SpeakerID)

# bylo by hezké přejmenovat hodnoty faktorové proměnné Sex
data$Sex <- factor(data$Sex)
levels(data$Sex)     # původní názvy

## [1] "f" "m"

data$Sex <- case_match(data$Sex, "f" ~ "žena", "m" ~ "muž")  # starší funkce: recode()

### --- vsuvka pro "chytré hlavy": Pokročilejší přejmenování levels (v úplně jiné hypotetické tabulce) ---
### Hypotetická situace: V tabulce "tab" mám šikovně zakódované francouzské mluvčí (sloupec speaker)
###                      kódy Fr1 až Fr8. Ty bych chtěl nechat, jak jsou.
###                      Také tam mám ale ne příliš vhodně pojmenované české mluvčí pod různými
###                      číselnými kódy začínajícími písmeny HC, např. HC108, HC215, HC232 apod.
###                      Ty bych rád automaticky přejmenoval na styl Cz1 atd., nehledě na původní čísla.
###                      Jak na to?
### indHC <- grep("HC", levels(tab$speaker), fixed = TRUE)      # Naleznu indexy jen těch, co obsahují HC
### levels(tab$speaker)[indHC] <- paste0("C", 1: length(indHC)) # Zjistím, kolik jich je, pomocí paste0()
###                                                             # slepím "C" s postupně rostoucími indexy,
###                                      # jména na pozici nalezených indexů nahradím těmito novými názvy.
### --- konec vsuvky ---

# a ještě přejmenovat proměnné
data <- data %>% rename(Vokál = Vowel, Pohlaví = Sex)

str(data)  # výrazně hezčí forma

## 'data.frame':    3000 obs. of  7 variables:
##  $ Vokál    : chr  "a" "a" "a" "a" ...
##  $ Pohlaví  : chr  "žena" "žena" "žena" "žena" ...
##  $ F0.mean  : num  212 219 252 244 231 ...
##  $ SpeakerID: Factor w/ 75 levels "1","2","3","4",..: 1 2 3 4 7 8 10 11 12 13 ...
##  $ F1_Hz    : num  819 713 748 791 788 ...
##  $ F2_Hz    : num  1354 1393 1495 1754 1622 ...
##  $ F3_Hz    : num  2595 2634 2497 2443 2961 ...

Obrázky

ggplot(data = data, aes(x = F1_Hz)) +  # volba dat a nejjednodušší aesthetics: na ose x jsou hodnoty F1_Hz
    geom_histogram()                   # přidání geom: histogram

Pro začátek to není špatné, ale máme všechny vokály v jednom histogramu. Co je zkusit odlišit barvou? Přidáme tedy pokyn do aesthetics, ať barvu nastaví podle Vokál.

ggplot(data = data, aes(x = F1_Hz, color = Vokál)) +  # ať nastaví barvu podle vokálu Vokál
    geom_histogram()                   # histogram

Aha, tím se nastavila ale jen barva okraje. Chtělo by to výplň. Podíváme se do “Cheat Sheet” a vidíme, že lze u histogramu nastavit také fill.

ggplot(data = data, aes(x = F1_Hz, fill = Vokál)) +  # raději tedy barvu výplně podle vokálu Vokál
    geom_histogram()                   # histogram

Nebylo by hezší zobrazit spíše odhad hustoty pravděpodobnosti (tzv. kernel density function)?

ggplot(data = data, aes(x = F1_Hz, color = Vokál)) +
    geom_density()                   # odhad hustoty pravděpodobnosti

ggplot(data = data, aes(x = F1_Hz, fill = Vokál)) + geom_density(alpha = 0.5)

Hezké, ale pořád máme muže a ženy dohromady. Udělal bych raději pro každou skupinu obrázek zvlášť, na řadu přicházejí facets.

ggplot(data = data, aes(x = F1_Hz, fill = Vokál)) +
    geom_density(alpha = 0.5) +
    facet_grid(. ~ Pohlaví)   # formát: řádky ~ sloupce. Tečka znamená "nerozdělovat". Zde tedy říkáme, že chceme 1 řádek a sloupce podle Pohlaví.

Vidíme, že ženy jsou výše. Když to jde tak snadno, co se vrátit k histogramům a sloupce udělat podle vokálu a řádky podle pohlaví?

ggplot(data = data, aes(x = F1_Hz, fill = Vokál)) +
    geom_histogram() +
    facet_grid(Pohlaví ~ Vokál)   # formát: řádky ~ sloupce.

ggplot(data = data, aes(x = F1_Hz, fill = Vokál)) + geom_density(alpha = 0.5) +
    facet_grid(Pohlaví ~ Vokál)

Zkusme krabicové grafy. Osa x (jednotlivé krabice) budou podle Vokál, y (hodnoty v krabici) zase F1_Hz.

ggplot(data = data, aes(x = Vokál, y = F1_Hz)) +
    geom_boxplot()

A opět rozdělíme podle pohlaví

ggplot(data = data, aes(x = Vokál, y = F1_Hz)) +
    geom_boxplot() +
    facet_grid(Pohlaví ~ .)

A co ušetřit místo a usnadnit porovnání tím, že krabice rozdělíme jen barvami?

ggplot(data = data, aes(x = Vokál, y = F1_Hz, color = Pohlaví)) +
    geom_boxplot()

Nebo zase raději fill?

ggplot(data = data, aes(x = Vokál, y = F1_Hz, fill = Pohlaví)) +
    geom_boxplot()

Občas narazíme také na „profesionální“ krabicové grafy (tzv. notchech box plots) znázorňující „vykousnutím“ interval spolehlivosti mediánu. Tím se dá krásně porovnávat, zda jsou středy skupin významně odlišné. Následující obrázek má ovšem jednu zásadní skrytou vadu.

ggplot(data = data, aes(x = Vokál, y = F1_Hz, fill = Pohlaví)) +
    geom_boxplot(notch = TRUE)

V čem je problém? Podíváme-li se na počty vokálů od jednotlivých mluvčích, zjistíme, že jsme se v grafu dopustili problému pseudoreplikace, jelikož každý mluvčí namluvil daný vokál vždy osmkrát. S tím ovšem interval spolehlivosti nepočítá, bere totiž každý řádek tabulky jako nezávislý náhodný výběr ze základního souboru. Zde je ovšem vždy 8 řádků závislých, takže interval vyjde přehnaně optimisticky úzký.

table(data$Vokál, data$SpeakerID)

##    
##     1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
##   a 8 8 8 8 8 8 8 8 8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8
##   e 8 8 8 8 8 8 8 8 8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8
##   i 8 8 8 8 8 8 8 8 8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8
##   o 8 8 8 8 8 8 8 8 8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8
##   u 8 8 8 8 8 8 8 8 8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8
##    
##     29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
##   a  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8
##   e  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8
##   i  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8
##   o  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8
##   u  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8
##    
##     54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
##   a  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8
##   e  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8
##   i  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8
##   o  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8
##   u  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8

Řešením je pro každou osobu vypočítat průměr nebo medián, takže nová tabulka bude ve formátu, že co řádek, to jiná osoba. Interval spolehlivosti uvidí, že pozorování je osmkrát méně a nebude ani zmaten přehnaně nízkým rozptylem položek (když vokál vysloví osmkrát jedna osoba, rozptyl bude jistě menší, než kdyby se jednalo o osm různých osob).

dataMedian <- data %>% group_by(SpeakerID, Pohlaví, Vokál) %>% summarise(F1_Hz = median(F1_Hz))
ggplot(data = dataMedian, aes(x = Vokál, y = F1_Hz, fill = Pohlaví)) +
    geom_boxplot(notch = TRUE)

Ovšem, kdo jednou zkusí houslové grafy, už žádné jiné nechce! Kromě úspornosti známé z krabicových grafů poskytují detailnější představu o tvaru rozdělení (zobrazují tzv. kernel density function). Že nevypadají jako houslové? To záleží na tvaru rozdělení dat.

ggplot(data = dataMedian, aes(x = Vokál, y = F1_Hz, fill = Pohlaví)) +
    geom_violin()

Jedná se však o hodně moderní záležitost a než si na ně společnost zvykne, bude opatrnější zobrazovat houslové grafy dohromady s krabicovými.

ggplot(data = dataMedian, aes(x = Vokál, y = F1_Hz, fill = Pohlaví)) +
    geom_violin() +
    geom_boxplot(notch=TRUE, position = position_dodge(width=.9), width=.15)

Zkusme jiné zobrazení, F1 vs. F2 pomocí teček (scatter plot).

ggplot(data = data, aes(x = F2_Hz, y = F1_Hz, color = Vokál)) +
    geom_point()

A co nastavit tvar bodu podle vokálu, ať to není rozlišno jen barvou? A obrátíme osy, aby to bylo jak z učebnice i-e-a-o-u zleva doprava a osa y odpovídala vertikální poloze jazyka.

ggplot(data = data, aes(x = F2_Hz, y = F1_Hz, color = Vokál, shape = Vokál)) +
    geom_point() +
    scale_x_reverse() +
    scale_y_reverse()

Krásné. Teď už by to snad jen chtělo přidat elipsy znázorňující 95 % dat. Přidáme tedy stats.

ggplot(data = data, aes(x = F2_Hz, y = F1_Hz, color = Vokál, shape = Vokál)) +
    geom_point() +
    scale_x_reverse() +
    scale_y_reverse() +
    stat_ellipse()

A nešly by udělat vyplněné?

ggplot(data = data, aes(x = F2_Hz, y = F1_Hz, color = Vokál, shape = Vokál)) +
    geom_point() +
    scale_x_reverse() +
    scale_y_reverse() +
    stat_ellipse(geom = "polygon", aes(fill = Vokál))

A průhledně? :-) Paramerem alpha to zařídíme.

ggplot(data = data, aes(x = F2_Hz, y = F1_Hz, color = Vokál, shape = Vokál)) +
    geom_point() +
    scale_x_reverse() +
    scale_y_reverse() +
    stat_ellipse(geom = "polygon", aes(fill = Vokál), alpha = 0.5)

To je ale překryv! Taková je realita. Můžeme samozřejmě znázornit menší procento dat, třeba jen 68 %.

ggplot(data = data, aes(x = F2_Hz, y = F1_Hz, color = Vokál, shape = Vokál)) +
    geom_point() +
    scale_x_reverse() +
    scale_y_reverse() +
    stat_ellipse(geom = "polygon", aes(fill = Vokál), alpha = 0.5, level = 0.68)

Teď by bylo zajímavé nechat jen elipsy a podívat se na rozdíl mezi muži a ženami. Pomocí size nastavíme tlustší čáru.

ggplot(data = data, aes(x = F2_Hz, y = F1_Hz, color = Pohlaví, linetype = Pohlaví)) +
    scale_x_reverse() +
    scale_y_reverse() +
    stat_ellipse(geom = "polygon", aes(fill = Vokál), alpha = 0.5, level = 0.68, size = 1)

Ještě by bylo pěkné, kdyby byla legenda seřazena i-e-a-o-u. Tak to ale musíme upravit data! Faktorové proměnné Vokál vynutíme pořadí hodnot.

head(data$Vokál)  # skutečně, ggplot2 bere pořadí tak, jak jsou uvedeny Levels

## [1] "a" "a" "a" "a" "a" "a"

data$Vokál <- factor(data$Vokál, ordered = TRUE, levels = c("i", "e", "a", "o", "u"))
head(data$Vokál)  # teď už to bude pěkné

## [1] a a a a a a
## Levels: i < e < a < o < u

ggplot(data = data, aes(x = F2_Hz, y = F1_Hz, color = Pohlaví, linetype = Pohlaví)) +
    scale_x_reverse() +
    scale_y_reverse() +
    stat_ellipse(geom = "polygon", aes(fill = Vokál), alpha = 0.5, level = 0.68, size = 1)

Themes

Jednou funkcí můžeme také změnit celkový styl obrázků

theme_minimal()

ggplot(data = data, aes(x = F2_Hz, y = F1_Hz, color = Pohlaví, linetype = Pohlaví)) +
    scale_x_reverse() +
    scale_y_reverse() +
    stat_ellipse(geom = "polygon", aes(fill = Vokál), alpha = 0.5, level = 0.68, size = 1) +
    theme_minimal()

theme_classic()

ggplot(data = data, aes(x = F2_Hz, y = F1_Hz, color = Pohlaví, linetype = Pohlaví)) +
    scale_x_reverse() +
    scale_y_reverse() +
    stat_ellipse(geom = "polygon", aes(fill = Vokál), alpha = 0.5, level = 0.68, size = 1) +
    theme_classic()

Příklad 4: rovnoměrné měřítko os

Pokud nám jde o sklon v nějaké závislosti, může být defaultní nastavení měřítek os matoucí. Zde je na x-ové ose rozsah 10 mnohem širší než na y-ové ose.

d <- data.frame(x = seq(from = 0, to = 10, length.out = 10), y = seq(from = 0, to = 20, length.out = 10))
ggplot(data = d, aes(x, y)) +
    geom_line()

Pokud nastavíme stejné měřítko 1:1 na obou osách, situace bude zřejmá: rozpětí 10 je nyní na výšku i na šířku stejně velké.

d <- data.frame(x = seq(from = 0, to = 10, length.out = 10), y = seq(from = 0, to = 20, length.out = 10))
ggplot(data = d, aes(x, y)) +
    geom_line() +
    coord_equal()

Příklad 5: sloupcový graf z již vysčítaných dat

Data obsahují sečtené počty. Jedná se o záznam počtu hodů kostkou u dvou hráčů v člověče.

Zkusíme hned čtyři způsoby zobrazení, můžeme si vybrat.

kostka <- data.frame(hrac = factor(rep(c(1,2), each = 6)), hodnota = factor(rep(1:6, times = 2)),
                     pocet = c(15, 9, 11, 16, 18, 11,   12, 11, 9, 14, 21, 13))

ggplot(data = kostka, aes(x = hodnota, y = pocet, fill = hrac)) +
    geom_bar(stat = "identity")

ggplot(data = kostka, aes(x = hodnota, y = pocet, fill = hrac)) +
    geom_bar(stat = "identity", position = "dodge")

ggplot(data = kostka, aes(x = hrac, y = pocet, fill = hodnota)) +
    geom_bar(stat = "identity")

ggplot(data = kostka, aes(x = hrac, y = pocet, fill = hodnota)) +
    geom_bar(stat = "identity", position = "dodge")

Příklad 6: sloupcové grafy a konfidenční intervaly binomického rozdělení

Četnost slova A v korpusu je 4114 z celkové velikosti 26617523, četnost slova B je 9624 z celkové velikosti korpusu 26551540. Porovnejme je na hladině významnosti \(\alpha = 0.01\).

V korpusech je zvykem místo procent či promilí používat i.p.m. (instances-per-milion), 1 i.p.m. je tedy tisícina promile. Naše poměry pro převod na i.p.m. vynásobíme milionem (1e6).

library(Hmisc)  # rychlejší výpočet binomických konfidenčních intervalů

fA <- 4114
sA <- 26617523

fB <- 9624
sB <- 26551540

alpha <- 0.01

ciA <- binconf(fA, sA, alpha = alpha, method = "exact") * 1e6
ciB <- binconf(fB, sB, alpha = alpha, method = "exact") * 1e6

tab <- data.frame(Word = c("Slovo A", "Slovo B"), i.p.m. = c(ciA[1], ciB[1]),
                  ciL = c(ciA[2], ciB[2]), ciH = c(ciA[3], ciB[3]))

ggplot(tab, aes(x = Word, y = i.p.m., ymin = ciL, ymax = ciH, fill = Word)) +
  geom_bar(stat = "identity", orientation = "x") + geom_errorbar(width=.5) +
  scale_fill_grey(start = 0.35, end = 0.65)

Příklad 7: konfidenční intervaly k percepčnímu experimentu 2AFC (two-alternative forced choice)

25 hodnotících (ID) mělo u každé položky (dvojice originální a manipulovaný zvuk) hodnotit, který ze dvou zvuků jim připadá kratší. Na výběr byla tlačítka První a Druhý. Každá položka byla v experimentu zařazena 2x s prohozeným pořadím zvuků vzhledem k tlačítkům (order), přičemž samozřejmě hodnotitel nevěděl, co je originál a co je manipulace (ani jsme jim nezdůrazňovali, že položky jsou nějak manipulované). Do tabulky jsme si pro přehlednost uložili odpověď hodnotícího (response) už rovnou jako shorter / longer (nehledě na pořadí tlačítek říká, zda jim manipulovaný zvuk přišel jako kratší / delší vzhledem k originálu). Celkové pořadí testovaných položek měl každý jednotlivý hodnotitel zcela náhodné.

Spočítejme průměrné hodnocení každé položky. Rozdělíme do skupin podle hodnoceného slova (word) a typu manipulace (manipulation).

Odpovědi “shorter” jsme přiřadili hodnotu -1, odpovědi “longer” +1 a vypočetli průměr z každé dvojice (v rámci hodnotícího). Výsledné možnosti jsou tedy následující: -1 = v obou případech označena manipulovaná položka jako kratší, 0 = přestože se jednalo o stejnou položku, jednou byla hodnocena tak a podruhé naopak (hodnotící měl evidentně problém slyšet rozdíl), +1 = v obou případech označena jako delší. Pochopitelně je jasné, že ani odpověď -1 nedává jistotu, že hodnotící slyšel jasně a netipoval. Mohl se jen náhodou 2x trefit na shodnou odpověď. Pokud je ale skutečně položka problematická (nerozhodná), ukáže se to v celkové statistice průměru ze všech hodnotících, kde konfidenční interval střední hodnoty bude zahrnovat hodnotu 0, což znamená, že nelze vyvrátit nulovou hypotézu o tom, že rozdíl není významný.

Přestože se jedná o diskrétní škálu (3 možnosti), průměrováním od většího množství hodnotitelů začíná veličina nabývat kontinuálního rázu.

Výpočet konfidenčních intervalů je proveden metodou bootstrap ze skupiny matematických metod využívajících tzv. resampling dat, tyto metody jsou nezávislé na statistickém rozdělení hodnot. Následující skript by tedy bylo možné použít např. i pro kontinuální škálu response, nejen pro náš případ stupnice tří hodnot.

library(readxl)
library(Hmisc)
tab <- read_excel("results.xlsx")
tab # náhled tabulky

## # A tibble: 2,050 × 5
##    ID    word   manipulation order response
##    <chr> <chr>  <chr>        <chr> <chr>   
##  1 0     provas 5a           1     longer  
##  2 0     provas 5a           2     shorter 
##  3 0     provas 5b           1     shorter 
##  4 0     provas 5b           2     shorter 
##  5 0     provas 1a           1     shorter 
##  6 0     provas 1a           2     shorter 
##  7 0     provas 1b           1     longer  
##  8 0     provas 1b           2     shorter 
##  9 0     provas 3a           1     longer  
## 10 0     provas 3a           2     shorter 
## # ℹ 2,040 more rows

tab2 <- tab %>%
    mutate(response = if_else(response == "shorter", -1, 1)) %>%
    group_by(ID, word, manipulation) %>% summarise(response = mean(response))
# odpovědi "shorter" jsme přiřadili hodnotu -1, odpovědi "longer" +1 a vypočetli průměr z každé dvojice
tab2

## # A tibble: 1,025 × 4
## # Groups:   ID, word [100]
##    ID    word   manipulation response
##    <chr> <chr>  <chr>           <dbl>
##  1 0     provas 1a                 -1
##  2 0     provas 1b                  0
##  3 0     provas 2a                  0
##  4 0     provas 2b                 -1
##  5 0     provas 3a                  0
##  6 0     provas 3b                  0
##  7 0     provas 3b_entire          -1
##  8 0     provas 4a                  0
##  9 0     provas 4b                 -1
## 10 0     provas 5a                  0
## # ℹ 1,015 more rows

alpha <- 0.05 / 41   # Bonferroniho korekce pro opakované testy, protože chceme, aby alpha = 0.05
                     # platila pro celý obrázek, nikoliv pro každý interval, kterých je celkem 41.

# musíme vytvořit vlastní funkci pro zobrazovanou statistiku,
# protože chceme zobrazit přes sebe jak čáru s průměrem, tak "errorbar" (vousy),
# což by se v ggplot počítalo každé zvlášť, a protože bootstrap používá generátor náhodných čísel,
# každá věc by vyšla malinko jinak posunutá
mean_boot <- function(x) {
    set.seed(1)  # trik - reset generátoru náhodných čísel na stejný výchozí bod
    return(mean_cl_boot(x, conf.int = 1-alpha))
}

ggplot(tab2, aes(x = word, y = response)) +
    geom_abline(intercept = 0, slope = 0, color = "darkgray") +  # referenční vodorovná čára "0"
    stat_summary(fun.data = "mean_boot", fatten = 2) +           # rovná čára s tečkou
    stat_summary(fun.data = "mean_boot", geom = "errorbar", width = .2) +   # vousy
    facet_wrap(~ manipulation)

Příklad 8: přehledné znázornění párů

Graf vhodný pro znázorněný párů - čárami spojíme příslušné body z 1. a 2. skupiny. Do aes() je nutno v tomto případě přidat ještě parametr group, aby funkce věděla, podle čeho čáry tvořit.

ggplot(sleep, aes(x = group, y = extra, fill=ID, color=ID)) +  # 
    geom_point(size=10, pch=21, alpha=.5) +
    geom_line(aes(group = factor(ID)), size=1)

Možná by bylo přehlednější napsat čísla přímo k bodům. Na barvy totiž není dobré spoléhat.

ggplot(sleep, aes(x = group, y = extra, color=ID)) +
    geom_line(aes(group = factor(ID)), size=1) +
    geom_text(aes(label = ID), size = 7)

Výhodou těchto exploratorních grafů je, že přestože je mezi jednotlivci velká variabilita, takže při pouhém porovnání středů obou skupin bychom si moc rozdílností jistí nebyli, v párovém grafu je rostoucí trend krásně patrný. Dopředu tak lze dobře odhadnout, že párový t-test by rozdíly mezi skupinami našel. Pomocí této exploratorní analýzy bychom také velice rychle nalezli a identifikovali případné podezřelé „ustřelené “jedince (outliers).

Příklad 9: propojení středů skupin a konfidenční intervaly

A jak na náš oblíbený graf uvedený v kapitole 16. Testy středních hodnot?

library(Hmisc)  # pro výpočet konfidenčních intervalů

ToothGrowth$dose <- factor(ToothGrowth$dose)
ggplot(ToothGrowth,
    aes(y = len, x = supp, color = dose, shape = supp, group = dose)) +
    geom_point(size=3) +
    stat_summary(fun = mean, geom = "line") +
    stat_summary(fun.data = mean_cl_boot, geom = "errorbar", width = 0.1)

Příklad 10: chytré umístění popisků

Pokud obsahuje graf větší množství popisků, knihovna ggrepel umí doslova zázraky. Snaží se je umístit tak, aby se nepřekrývaly, grafy pak působí výrazně lépe.

Pro více příkladů viz [https://cran.r-project.org/web/packages/ggrepel/vignettes/ggrepel.html]

library(ggrepel)

dat <- subset(mtcars, wt > 2.75 & wt < 3.45)
dat$car <- rownames(dat)

p <- ggplot(dat, aes(wt, mpg, label = car)) +
  geom_point(color = "red")

p1 <- p + geom_text() + labs(title = "geom_text()")

p2 <- p + geom_text_repel() + labs(title = "geom_text_repel()")

gridExtra::grid.arrange(p1, p2, ncol = 2)

ggplot(mtcars, aes(x = wt, y = 1, label = rownames(mtcars))) +
  geom_point(color = "red") +
  geom_text_repel(
    nudge_y      = 0.05,
    direction    = "x",
    angle        = 90,
    vjust        = 0,
    segment.size = 0.2
  ) +
  xlim(1, 6) +
  ylim(1, 0.8) +
  theme(
    axis.line.y  = element_blank(),
    axis.ticks.y = element_blank(),
    axis.text.y  = element_blank(),
    axis.title.y = element_blank()
  )

Další příklady

Pěkná kuchařka na grafy v ggplot2: http://www.cookbook-r.com/Graphs/

A colorblind-friendly palette

Není dobré odlišovat skupiny jen barvami. Nemalé procento mužské populace má snížený barvocit a některé barevné kombinace od sebe nerozpozná (vypadají zcela stejně). Dalším problémem může být černobílý tisk.

Pro alespoň částečné usnadnění rozpoznání barev byla vytvořena následující barevná škála. Více informací zde: http://jfly.iam.u-tokyo.ac.jp/color/ http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/

# The palette with grey
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

# The palette with black
cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

Pro výplně: přidat + scale_fill_manual(values = cbPalette)

Pro čáry a body: přidat + scale_color_manual(values = cbPalette)

ggplot(data = dataMedian, aes(x = Vokál, y = F1_Hz, fill = Pohlaví)) +
    geom_violin() +
    geom_boxplot(notch=TRUE, position = position_dodge(width=.9), width=.15) +
    scale_fill_manual(values = cbPalette)

Ukládání obrázků pro publikace

Potřebujeme nastavit finální rozměr obrázku, velikosti fontů, typ fontu shodný s textem publikace, pohrát si s legendou, nastavit šedivou škálu, aby byly kategorie rozpoznatelné i při černobílém tisku. Ideální je uložit obrázek ve vektorovém formátu, aby nebyl při zvětšování zubatý.

Zde záměrně ukazuji pro inspiraci pouze fragmenty kódů, které řeší zajímavé problémy. Jsou inspirovány naším výzkumem, ale obsahují pozměněné hodnoty, nereflektují reálná data.

Př. 1 - Basic plotting system

Čtyři proměnné (pt1 až pt4) obsahující vždy složku s vektorem časových okamžiků $t a vektor příslušných hodnot frekvencí $f. Každou proměnnou zobrazíme v samostatném panelu. Vykreslení každého dílčího grafu zajišťuje funkce jedenPt, parametr index určující pořadí se uvnitř přetaví na popis ve stylu (a), (b) atd., který zapisuji na místo běžné legendy. Funkce zajišťuje jednotný rozsah osy y. Popis os x a y zde záměrně vypínám, z důvodu maximální úspory místa doplňuji společný popis os všech čtyř panelů v samotném závěru funkcemi mtext().

Rozměry obrázku jsou pevně dané: 9 * 5 cm (velmi limitující konferenční dvousloupcový formát), velikost fontu fixně 8 (automaticky nastavená velikost by mohla vést na nečitelně malé fonty). Funkcí expression() vytvářím hezký matematický formát základní frekvence f0. Parametrem mar nastavuji velikost okrajů obrázků, aby byly akorát natěsno. Výsledek ukládám do vektorového formátu svg (lze editovat v InkScape a Word ho též v pořádku zobrazí) a eps (lze editovat v CorelDRAW).

windowsFonts(Times=windowsFont("Times New Roman"))  # nastavení skupiny fontů z dostupných v OS Windows
par(family = "Times")
fontSize <- 8

jedenPt <- function(pt, index) {
    plot(pt$t, pt$f, xlab = "", ylab = "", pch = 16, type = "b", lty = 3,
         ylim = c(20, 33), xlim = c(pt$tmin, pt$tmax))
    m <- lm(pt$f ~ pt$t)  # lineární regrese
    abline(m, col = "black", lwd = 2, lty = 1) # vykreslení přímky
    
    legend("topright", legend = paste0("(", letters[index], ")"), bty = "n", cex = 1.2); # označení (a), (b) atd.
}

svg(filename = "basic1.svg", width = 9.0 / 2.54, height = 5 / 2.54, pointsize = fontSize)
par_mfrow <- par()$mfrow; par_mar <- par()$mar; par(mfrow = c(2,2), mar = c(3.1, 3, 0.1, 0.1), family = "Times")
jedenPt(pt1, 1); jedenPt(pt2, 2); jedenPt(pt3, 3); jedenPt(pt4, 4)
mtext(side = 2, expression(italic("f")[0]*" (ST re 50 Hz)"), line = -1.2, cex = 1, outer = TRUE)
mtext(side = 1, "Time (sec)", line = -1.2, cex = 1, outer = TRUE)
par(mfrow = par_mfrow, mar = par_mar)
dev.off()

cairo_ps(filename = "basic1.eps", width = 9.0 / 2.54, height = 5 / 2.54, pointsize = fontSize, family = "Times")
par_mfrow <- par()$mfrow; par_mar <- par()$mar; par(mfrow = c(2,2), mar = c(3.1, 3, 0.1, 0.1))
jedenPt(pt1, 1); jedenPt(pt2, 2); jedenPt(pt3, 3); jedenPt(pt4, 4)
mtext(side = 2, expression(italic("f")[0]*" (ST re 50 Hz)"), line = -1.2, cex = 1, outer = TRUE)
mtext(side = 1, "Time (sec)", line = -1.2, cex = 1, outer = TRUE)
par(mfrow = par_mfrow, mar = par_mar)
dev.off()

Př. 2 - Rozměry a velikost fontů v ggplot2

Nastavení fontů v ggplot2, víceřádkové popisky hodnot faktorových proměnných a uložení do svg a eps s rozměry 9 * 8 cm.

Návody na detailní nastavení legendy:

• http://www.sthda.com/english/wiki/ggplot2-legend-easy-steps-to-change-the-position-and-the-appearance-of-a-graph-legend-in-r-software
• https://www.datanovia.com/en/blog/ggplot-legend-title-position-and-labels/

tab <- read.table("data.txt", header = TRUE, sep = "\t")
tidy$Method <- factor(tidy$Method, ordered = TRUE, levels = c("M1", "M2", "M3", "M4")) # určuji pořadí zobrazení

# přejmenovávám hodnoty, protože chci, aby byl popisek zalomený na 2 řádky
tidy$Condition <- case_match(tidy$Condition, "SP-P1 (n = 376)" ~ "SP-P1\n(n = 376)",
                         "SP-P2 (n = 404)" ~ "SP-P2\n(n = 404)",
                         "RT-P1 (n = 729)" ~ "RT-P1\n(n = 729)",
                         "RT-P2 (n = 1310)" ~ "RT-P2\n(n = 1310)")

# nastavuji font a jeho velikost
windowsFonts(Times=windowsFont("Times New Roman"))
fontSize <- 8
g <- ggplot(tidy, aes(x = Condition, y = Gradient, fill = Method)) +
  geom_bar(stat = "identity", position = position_dodge()) +
  facet_grid(Parameter ~ ., scales = "free") + theme_bw(base_family = "Times") + scale_fill_grey() +
  theme(axis.title=element_text(size=fontSize), axis.text = element_text(size=fontSize),
    legend.text = element_text(size=fontSize), legend.title = element_text(size=fontSize),
    strip.text = element_text(size=fontSize),
    legend.position = "right", legend.key.width = unit(12, "points"), legend.key.height = unit(12, "points"),
    legend.margin = margin()) + labs(y = "Slope (dB/dec)", x = "")

# ukládám jako svg s rozměry 9 * 8 cm
svg(filename = "gg2.svg", width = 9.0 / 2.54, height = 8 / 2.54, pointsize = fontSize)
print(g)
dev.off()

# Ukládám do eps (také vektorový formát, který bez problémů otevře i CorelDRAW).
# Minimálně pod Windows je ale zásadní specifikovat device = cairo_ps, jinak Corel může zkazit barvy.
ggsave(filename = "gg2.eps", width = 9, height = 8, units = "cm", plot = g, device = cairo_ps)

Př. 3 - Více obrázků ggplot2 v jednom Figure

Zde jsem použil balíček gridExtra, u prvního obrázku g1 jsem vypnul legendu, která bude společná s g2. Kvůli tomu pak na závěr ve funkci grid.arrange volím ručně poměr šířek jednotlivých grafů.

Pro maximální úsporu místa vkládám popisky (a) a (b) “ručně”, což je poněkud krkolomné, jelikož se musí zadávat souřadnice.

Určení velikosti fontu v annotate() je v jiných jednotkách a je potřeba použít šikovný přepočet.

Automatické číslování by šlo řešit pomocí jiných balíčků, např. cowplot nebo multipanelfigure, viz například diskuze zde: https://stackoverflow.com/questions/1249548/side-by-side-plots-with-ggplot2

library(gridExtra)
windowsFonts(Times=windowsFont("Times New Roman"))
fontSize <- 8
legendKeySize <- 12

g1 <- ggplot(...) + geom_bar(stat = "identity", position = position_dodge()) +
    theme_bw(base_family = "Times") + scale_fill_grey() + 
    theme(axis.title=element_text(size=fontSize), axis.text = element_text(size=fontSize),
          legend.position = "none", strip.text = element_text(size=fontSize)) +
    labs(y = "Slope (dB/dec)", x = "LOR (months)") +
    annotate(geom = "text", x = "6+", y = -3.75, label = "(a)", size = 8/72.27 * 25.4)

g2 <- ggplot(...) + geom_bar(...) +
    theme_bw(base_family = "Times") + scale_fill_grey() + 
    theme(axis.title=element_text(size=fontSize), axis.text = element_text(size=fontSize),
          legend.text = element_text(size=fontSize), legend.title = element_text(size=fontSize),
          strip.text = element_text(size=fontSize), legend.position = "right",
          legend.key.width = unit(legendKeySize, "points"), legend.key.height = unit(legendKeySize, "points"),
          legend.margin = margin()) + labs(y = "Gradient (ST/sec)", x = "LOR (months)") +
    annotate(geom = "text", x = "6+", y = -1.85, label = "(b)", size = 8/72.27 * 25.4)

svg(filename = "gg3.svg", width = 9.0 / 2.54, height = 5 / 2.54, pointsize = fontSize)
w1 = 0.39; grid.arrange(g1, g2, ncol = 2, widths = c(w1, 1-w1))
dev.off()

Př. 4 - Seřazení skupin dle velikosti hodnot

Obrázek se může velice zpřehlednit, seřadíme-li skupiny dle velikosti zobrazovaných hodnot.

library(tidyverse)
tidy <- read.table("data.txt", header = TRUE, sep = "\t")

head(tidy)

##   Speaker Language  Mean  ymin  ymax
## 1    CZ01       CZ 5.694 5.301 6.087
## 2    CZ02       CZ 6.961 6.432 7.490
## 3    CZ03       CZ 5.758 5.480 6.036
## 4    CZ04       CZ 7.363 7.183 7.543
## 5    CZ05       CZ 6.766 6.427 7.105
## 6    CZ06       CZ 4.966 4.614 5.318

# Seřadíme zaprvé podle Language (2 skupiny), dále dle hodnoty Mean v rámci každé skupiny.
# V seřazené tabulce pak trikem definujeme faktorovou proměnnou s daným pořadím hodnot,
# takže následně bude ggplot toto pořadí brát v potaz.
tidy <- tidy %>% arrange(Language, Mean) %>% mutate(Speaker = factor(Speaker, levels = unique(Speaker)))

windowsFonts(Times=windowsFont("Times New Roman"))
fontSize <- 8

g <- ggplot(tidy, aes(x = Speaker, y = Mean, ymin = ymin, ymax = ymax)) +
    geom_point() +geom_errorbar(width=0.2) +
    theme_bw(base_family = "Times") + scale_fill_grey() +
    theme(text=element_text(size = fontSize),
          axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
    labs(y = "AR (syll/sec)")
    

svg(filename = "gg4.svg", width = 9.0 / 2.54, height = 5 / 2.54, pointsize = fontSize)
print(g)
dev.off()

---

© 19. 4. 2021 Tomáš Bořil, borilt@gmail.com