Čas od času se mi stane, že na mě někdo obrátí s prosbou o statistické vyhodnocení dat, které zahrnuje větší množství opakovaných operací. Typicky taková zadání vypadá nějak takhle: Tady jsou data, kde je pět různých kategorií a u proměnné x je potřeba zjistit, jestli se hodnota mezi kategoriemi liší. Tohle zadání sice není žádná hrůza a dá se zpracovat poměrně rychle i ručně. Ale s ohledem na moji lenost a obecnou neochotu dělat něco, co po praktické ukázce zvládne v podstatě každý, jsem tomuhle řešení nebyl nikdy extra nakloněn. Navíc poslední dataset, co jsem dostal do ruky obsahoval kategorií fakt hodně (stovky) a ve výsledku bylo zajímavých jen několik (nízké desítky) kombinací, kde daný test vyšel podle očkávání.

Dřív jsem tyhle úlohy řešil v R nějakou kombinací vnořených cyklů a tvorbou výstupních datových rámců v průběhu těchto cyklů. Zásadní nevýhodou je, že tohle řešení bylo poměrně specifické pro danou situaci a bylo potřeba ho dost upravovat při přenesení na nový problém. Částečně asi i proto, že jsem si nikdy nevytvořil pořádný a ucelený koncept, jak tyhle problémy řešit. Takže jsem pokaždé strávil dost času nad úpravami kódu.

Zároveň mi vždycky v hlavě ležela představa, že by mělo existovat jednoduché řešení pomocí nested data, jak ukazuje nápověda tidyr. Existence výše popsáneho komplexního problému a trochy volného času mě nakopla k tomu, dát těmhle problémům konečně strukturu a najít jednoduché řešení, které bude dobře škálovat.

Na začátku vyřešme potřebné balíky. dplyr a tidyr zajišťují základní manipulaci s daty. purrr umožní rozumné použití funkcinálního programování nad daty, takže odpadá nutnost psát cykly. glue a DT jsou jen kosmetické balíčky jednak na tvorbu textu v datech a výstupy tabulek do html.

library(dplyr)
library(tidyr)
library(purrr)
library(glue)
library(DT)

Začneme u dat, zde si vezmeme “provařenou” datovou sadu mtcars. Kategorie na datech nám určí sloupec cyl (počet válců) a zajímat nás bude spotřeba (mpg). Vypočteme průměr spotřeby za jednotlivé kategorie.

data(mtcars)

mtcars %>% 
  group_by(cyl) %>% 
  summarise(mean = mean(mpg)) %>% 
  datatable(options = list(scrollX = TRUE), rownames = FALSE)

Úkol nyní spočívá v otestování, zda-li je pro kategorii x průměrná spotřeba vyšší než pro kategorii y. Data mají 3 kategorie, což znamená 6 testů (netřeba testovat x proti x). Při tomhle zadání, není test x $\times$ y adekvátní testu y $\times$ x, protože testujeme ne pouze na rozdílnost, ale přímo na variantu vyšší než.

Základem pro testy je vytvořit si tzv. grid kombinací kategorické proměnné. Tenhle problém řeší funkce expand_grid z balíku tidyr. Jednotlivé sloupce vytvoříme jako unikátní hodnoty ze sloupce cyl a odfiltrujeme situace, kdy jsou obě proměnné shodné.

cyl_grid <- expand_grid(cyl1 = unique(mtcars$cyl), cyl2 = unique(mtcars$cyl)) %>%
  filter(cyl1 != cyl2)

Následující krok je klíčový. Vytvoříme si vnořené (nested) datové rámce, což znamená, že jeden sloupce v našich datech bude obsahovat celé datové rámce. Tahle vlastnost a možnost R je naprosto kouzelná. Tento sloupce vytvoříme tak, že pomocí funkce map2() (balík purrr) budeme iterovat současně přes dvě proměnné datového ramce vytvořeného v předešlém kroku. Na každý řádek aplikujeme jednoduchou funkci, která bere původní dataset mtcars a filtruje z něj pouze hodnoty určené proměnnými cyl1 a cyl2. Na zápisu funkce jsou důležité dvě věci, začíná znakem ~, což je v tomto stylu zápisu povinné, a na proměnné se odkazujeme jako .x (první parametr funkce map2()) a .y (druhý parametr funkce map2()).

Pomocí funkcí mutate() a factor() převedeme původně číselné hodnoty cyl na kategorie. Při převodu na kategorie, což většinou vyžadují statistické testy, je důležité určit i tzv. levels, které učují pořadí kategorií, což je opět důležité pro testy. Zde je nastavíme tak, aby bylo dodrženo pořadí z datového rámce. Test pak bude odpovídat na otázku cyl1 $>$ cyl2.

Pozor na to, že pokud se ve funkci factor() parametr levels nenastaví ručně, tak se odvodí z dat, podle toho, v jakém pořadí se vyskytují hodnoty. Což u některých operací nemusí vadit, ale v zrovna pro testy to problém je, respektive může být. Pokud je takto specifikujeme, tak získáváme nad daty větší kontrolu.

mtcars_combinations <- cyl_grid %>%
  mutate(df_data = map2(cyl1, cyl2,
                        ~ mtcars %>% 
                          filter(cyl %in% c(.x, .y)) %>% 
                          mutate(cyl = factor(cyl, c(.x, .y))) 
  ))

Výsledek vypadá následovně. Tabulku nelze rozumně zobrazit, ale třeba RStudio ji umožňuje rozkliknout. To co na html výstupu vidíme jako [object Object] je ve skutečnosti vnořený datový rámec.

mtcars_combinations %>% 
  datatable(options = list(scrollX = TRUE), rownames = FALSE)

Důkazem může být vytažení třetího sloupce z prvního řádku a jeho zobrazení. Při kontrole opravdu obsahuje jen řádky s hodnotami cyl 4 a 6, což odpovídá předchozí tabulce a hodnotám cyl1 a cyl2 pro daný řádek.

mtcars_combinations[[3]][[1]] %>% 
  datatable(options = list(scrollX = TRUE), rownames = FALSE)

S nachystanými daty se můžeme pustit do testů. Jako příklad použijeme wilcox.test(), ale lze zaměnit za test jiný, či za jinou operaci (třeba výpočet regrese funkcí lm()). Tentokrát použijeme funkci map(), do které bude vstupovat sloupce s vnořenými datovými rámci. Funkce, která se bude iterativně aplikovat na jednotlivé prvky daného sloupce je samotný test (nezapomenout na uvození ~). Definice testu je zcela klasická, pouze u parametru data udáme .x, což je vnořený datový rámec, ostatní parametry testu můžeme zadat dle libosti. Výsledkem bude nový sloupec v datech (wilcox_test), kde bude uložen výsledek testu, v tomto případě objekt - seznam.

mtcars_tests <- mtcars_combinations %>% 
  mutate(wilcox_test = map(df_data, 
                           ~ wilcox.test(mpg ~ cyl,
                                         data = .x,
                                         alternative = "greater"))
  )

Výsledek testu v podobě seznamu není úplně užitečný, spíš nás budou zajímat specifické údaje. Ty můžeme extrahovat do nových sloupců za pomocí funkce map() a jejich odvozenin map_*(). Zde kontrétně map_chr() na textové hodnoty a map_dbl() na hodnoty číselné. První parametr vždy specifikuje sloupec s uloženými výsledky testů a druhý název prvku, který chceme z výsledku testu “vytáhnout” (tyto názvy jsou v nápovědě daného testu v položce Value).

mtcars_tests <- mtcars_tests %>% 
  mutate(data_name = map_chr(wilcox_test, "data.name"),
         alternative = map_chr(wilcox_test, "alternative"),
         text_alternative = glue("{cyl1} {alternative} {cyl2}"),
         p_value = map_dbl(wilcox_test, "p.value"))

Výsledek můžeme zobrazit, sloupec s vnořeným datovým rámcem pro zjednodušení vynecháme. Sloupec wilcox_test opět nelze rozumně zobrazit, jak bylo popsáno dříve.

mtcars_tests %>% 
  select(-df_data) %>% 
  datatable(options = list(scrollX = TRUE), rownames = FALSE)

Na závěr vyfiltrujeme pouze hodnoty, které nás zajímají.

mtcars_tests %>% 
  select(-df_data) %>% 
  filter(p_value < 0.05) %>% 
  datatable(options = list(scrollX = TRUE), rownames = FALSE)

Výsledný kód je jendoduše čitelný, snadno upravitelný a dobře škálovatelný. Na cca 50 řádcích máme obecnou definici řešení netriviálního problému. Tím, že nepoužíváme složité cykly a většinu práce přenecháme specializovaným funkcím limitujeme možnost, že se dopustíme “hloupé” chyby, která se bude složitě zjišťovat a dohledávat.

V případě, že je dataset opravdů velký, doporučil bych u kroku, se vytváří vnořené datové rámce, přidat kromě funkce filter() do pipeline ještě funkci select() a vybrat pouze sloupce nezbytné pro test (zde konrétně mpg a cyl), aby vnořené datové rámce neobsahovaly zbytečně velké objemy dat.