Iterace

Jedním ze způsobů, jak redukovat duplicitní kód, jsou funkce. Dalším je iterace čili opakovaní stejné operace na různé vektory.

1 Balíčky

library(tidyverse)

2 For loops

Dejme tomu, že mámu tento jednoduchý dataframe s deseti řádky a čtyřmi sloupci.

df <- tibble(
  a = rnorm(10),
  b = rnorm(10),
  c = rnorm(10),
  d = rnorm(10)
)
df

># # A tibble: 10 × 4
>#          a      b        c      d
>#      <dbl>  <dbl>    <dbl>  <dbl>
>#  1 -0.240   1.74   0.989    1.12 
>#  2 -0.467  -0.459  1.82    -1.01 
>#  3 -0.264  -0.500 -1.06     1.06 
>#  4 -0.144   0.115 -1.94    -1.29 
>#  5  1.82    1.25  -1.38     0.350
>#  6  0.426  -0.473  0.163    0.281
>#  7  0.0374  0.870  0.105    2.11 
>#  8  0.0865  0.203 -1.35     0.841
>#  9  1.58    0.177  0.00523  0.843
># 10  0.712   1.86   1.65     1.03

Kdybyste chtěli spočítat např. medián každého sloupce, můžeme copy-pastovat kód.

median(df$a)

># [1] 0.06191261

median(df$b)

># [1] 0.1902308

median(df$c)

># [1] 0.05495897

median(df$d)

># [1] 0.841967

To ale není moc efektivní a navíc to zvyšuje pravděpodobnost, že někde uděláme chybu. Kromě toho jedno nepsané pravidlo říká, že bychom neměli copy-pastovat část kódu více než dvakrát.

Místo toho můžeme použít tzv. for loop (mohli bychom použít across(), ale i ta “na pozadí” využívá for loops). Nejdříve si připravíme prázdný číselný vektor output pomocí funkce double(), do kterého se bude ukládat výstup for loop a jehož délka se rovná počtu sloupců df (protože pro každý sloupec chceme vypočíst medián). (Jiné typy výsputních atomických vektorů bychom mohli připravit pomocí funkcí integer() nebo character().)

output <- double(length = length(df))
output 

># [1] 0 0 0 0

Takto pak vypadá samotná for loop

for (i in seq_along(df)) {            # 2. sekvence
  output[[i]] <- median(df[[i]])      # 3. tělo for loopu
}

output

># [1] 0.06191261 0.19023075 0.05495897 0.84196702

Pro vysvětlení for (i in seq_along(df)) vyvolává jednotlivé prvky ze sekvence seq_along(df), přičemž seq_along(df) vytvoří číselnou sekvenci s číselnými indexy sloupců v df.

seq_along(df)

># [1] 1 2 3 4

for loop vlastně vyvolává jednotlivé prvky i z seq_along(df)

for (i in seq_along(df)) {            
  print(i)    
}

># [1] 1
># [1] 2
># [1] 3
># [1] 4

Protože v těle funkce jsme definovali output[[i]] <- median(df[[i]]), vede to k tomu, že se do původně prázdného vektoru output ukládají popořadě mediány jednotlivých sloupců. Takže se vlastně děje toto:

1. output[[1]] <- median(df[[1]]).
2. output[[2]] <- median(df[[21]]).
3. output[[3]] <- median(df[[3]]).
4. output[[4]] <- median(df[[4]]).

Každá for loop má tři části.

1. Output (výstup). V našem případě jsme si jej vytvořili pomocí output <- double(length(x)). Musíme si připravit prázdný výstupní vektor, který má náležitý počet prvků a správný typ. Protože počítáme 4 mediány a máme spojité proměnné, připravili jsme si vektor typu double o délce = 4.
2. Sekvence. V našem případě i in seq_along(df). Ta určuje průběh for loop. Každá iterace for loop přiřadí objektu i jednu (a pouze jednu) z hodnot seq_along(df) a pak pokračuje k další, takže první iterace přidadí i první hodnotu ze seq_along(df), druhá iterace druhou hodnotu atd.
3. Tělo. V našem případě output[[i]] <- median(df[[i]]). To je vlastně hlavní část kódu, který dělá to, co chceme. Je spouštěň opakovaně, pokaždé s jinou hodnotou i. První iterace tedy vypadá takto output[[1]] <- median(df[[1]]), druhá takto output[[2]] <- median(df[[2]]) atd.

Cvičení

Zkuste napsat for loop, který vypočte průměr každého sloupce z datasetu mtcars. Základ pro kód máte níže. Ať se nám ve výstupu lépe orientuje, je vhodné si jeho prvky pojmenovat, např. pomocí set_names a to podle názvů sloupců z datasetu mtcars.

out <- double(length = ___) %>% 
  set_names(nm = ___)

for (i in ___) {
  out[[i]] <- ___
}

out

Zkuste napsat for loop, který určtí typ, typeof(), každého sloupce z datasetu starwars. Základ pro kód máte níže.

out <- double(length = ___) %>% 
  set_names(nm = ___)

for (i in ____) {
  out[[i]] <- typeof(___)
}

out

Zkuste dokončit for loop, který vypočte počet jedinečných hodnot z každé sloupce datasetu iris

out <- out <- integer(length = ___) %>% 
  set_names(nm = ___)

for (i in ____ ) {
  ____ <- n_distinct(___)
}

out

Dokončete for loop, který generuje 10 náhodných čísel z normální distribuce s různými průměry: -10, 0, 10, 100

mu <- c(-10, 0, 10, 100)
out <- vector("list", length = length(mu))

for (i in seq_along(mu)) {
  out[[i]] <- rnorm(n = ___, mean = ___)
}
out

3 Varianty for loops

Existují čtyři variace na základní for loops:

1. Modifikace stávajícího objektu, namísto tvorby nového
2. Loopování napříč jmény nebo hodnotami vektoru, namísto číselných indexů
3. Práce s výstupy neznámé délky.
4. Práce se sekvencemi neznémé délky

3.1 Iterace modifikující stávající objekt

Dejme tomu, že máme tento cvičný tibble.

df <- tibble(
  a = rnorm(10),
  b = rnorm(10),
  c = rnorm(10),
  d = rnorm(10)
)

A chceme každý sloupce tranformovat na hodnoty s minimem = 0 a maximem = 1. Vytvoříme si k tomu příslušnou funkci.

rescale01 <- function(x) {
  rng <- range(x, na.rm = TRUE)
  (x - rng[1]) / (rng[2] - rng[1])
}

Mohli bychom kód copy-pastovat, ale víme, že se to dělat nemá

df$a <- rescale01(df$a)
df$b <- rescale01(df$b)
df$c <- rescale01(df$c)
df$d <- rescale01(df$d)

Proto si vytvoříme for loop, kdy input bude zároveň outputem. A protože dataframe je vlastně vektor typu list s prvky (sloupci) stejné délky, můžeme iterační sekvenci vyvořit pomocí seq_along(df). A “těle” této for loop pak použijeme funkci rescale01().

for (i in seq_along(df)) {
  df[[i]] <- rescale01(df[[i]])
}
df

># # A tibble: 10 × 4
>#         a     b     c      d
>#     <dbl> <dbl> <dbl>  <dbl>
>#  1 0.520  0.476 0.256 0.225 
>#  2 0.0925 0     0.150 0.454 
>#  3 0.159  0.783 0     0.919 
>#  4 0.479  0.250 0.823 0.486 
>#  5 0.705  0.932 0.544 0.424 
>#  6 0.681  0.258 1     1     
>#  7 1      0.862 0.765 0.0237
>#  8 0      0.282 0.758 0.467 
>#  9 0.797  1     0.401 0.515 
># 10 0.855  0.295 0.558 0

3.2 Iterace napříč prvky samotnými

For loop pro samotné prvky for(x in xs) vovolává samotné prvky x z objektu xs (nikoli jejich číselné indexy nebo jména). Hodí se hlavně tehdy, když nechceme output ukládat, např. při tvorbě více grafů zároveň. Takto bychom vytvořili jednoduché histogramy každého prvku (sloupce) z datasetu mtcars

for (x in mtcars) {
  hist(x)
}

3.3 Iterace pomocí názvů prvků

for (nm in names(x)) vyvolává jednotlivá jména prvků. Dostaneme tak postupně jednotlivé názvy prvků, které můžeme využit # k subsettingu, např. x[[nm]]. Jména prvků pak můžeme vyuít k pojmenování výsledků, grafů nebo vytvořených souborů. Nemusíme si pak dělat starosti s tím, že bychom prvky vektoru připraveného pro output museli dopředu nebo dodatečně pojmenovat.

out <- double(length = length(mtcars))
out

>#  [1] 0 0 0 0 0 0 0 0 0 0 0

for (nm in colnames(mtcars)) {
  out[[nm]] <- mean(mtcars[[nm]])
}
out

>#                                                                              
>#   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000 
>#                                                    mpg        cyl       disp 
>#   0.000000   0.000000   0.000000   0.000000  20.090625   6.187500 230.721875 
>#         hp       drat         wt       qsec         vs         am       gear 
># 146.687500   3.596563   3.217250  17.848750   0.437500   0.406250   3.687500 
>#       carb 
>#   2.812500

3.4 Iterace pomocí číselných indexů

Iterace pomocí číselných indexů je základní formou iterace, protože na základě číselné pozice lze vyvolat jak názvy, tak hodnoty. Kódem níže nejprve přidělíme jména prvkům vektorou připraveného pro output names(out)[[i]] <- colnames(mtcars)[[i]] a potom “dosadíme” prvkům hodnoty, v tomto případě průměry sloupců mtcars(), pomocí out[[i]] <- mean(mtcars[[i]]).

out <- double(length = length(mtcars))
out

>#  [1] 0 0 0 0 0 0 0 0 0 0 0

for (i in seq_along(mtcars)) {
  names(out)[[i]] <- colnames(mtcars)[[i]]
  out[[i]] <- mean(mtcars[[i]])
}

out

>#        mpg        cyl       disp         hp       drat         wt       qsec 
>#  20.090625   6.187500 230.721875 146.687500   3.596563   3.217250  17.848750 
>#         vs         am       gear       carb 
>#   0.437500   0.406250   3.687500   2.812500

3.5 Iterace s neznámou délkou outputu

Někdy neznéme délku outputu. Např, kdybychom chtěli simulovat vektory náhodné délky. Řešením je uložit si výsledky do listu, protože ten může obsahovat prvky rozdílné délky, a pak případně prvky listu sloučit do jednoho atomického vektoru. Kódem níže si nejprve připravíme tři průměry, které chceme použít, a prázdný tříprvkový list pro output. V každé ze tří iterací pak for loop vynegeruje jedno náhodné celé číslo n od 1 do 100 a do outputu uloží n hodnot vylosovaných z normální distribuce s různými průměry (0, 1 a 3) a směrodatnou odchylkou 1 (tu jsme nenastavovali, ale je to defaultní hodnota ve funkci rnorm()). Jak vidíme níže, vylosovalo se nám 31 hodnot z distribuce s průměrem 0, 63 hodnot z distribuce s průměrem 1, a 48 hodnot z distribuce s průměrem 2.

set.seed(123)
means <- c(0, 1, 2)
out <- vector("list", length(means))

for (i in seq_along(means)) {
  n <- sample(1:100, size = 1) 
  out[[i]] <- rnorm(n, means[[i]])
}
str(out)

># List of 3
>#  $ : num [1:31] 0.801 1.19 -1.69 1.239 -0.109 ...
>#  $ : num [1:63] 1.9 1.88 1.82 1.69 1.55 ...
>#  $ : num [1:48] 1.54 2.31 1.92 2.41 2.18 ...

unlist(out) %>% head()

># [1]  0.8005543  1.1902066 -1.6895557  1.2394959 -0.1089660 -0.1172420

3.6 Iterace s neznámou délkou sekvence

Někdy neznáme ani délku vstupní sekvence hodnot, např. když děláme nějaké simulace. Přestavme si například, že chceme simulovat hody mincí a iterování ukončit, až když padne hlava třirkát za sebou. To nelze udělat pomocí for loop, ale pomocí while loop ano.

While loop vyžaduje stanovení podmínky, po jejíž naplnění se iterace ukončí

  while (podmínka) {
     hlavní kód
  }

While loop je ve skutečnosti obecnější než for loop, protože každý for loop lze předělat na while loop, ale naopak to jít nemusí. Kdybychom měli tento for loop pro výpočet průměrů všech sloupců mtcars.

out <- double(ncol(mtcars))
out

>#  [1] 0 0 0 0 0 0 0 0 0 0 0

for (i in seq_along(mtcars)) {
  out[[i]] <- mean(mtcars[[i]]) # Výpočet průměrů
  names(out)[[i]] <- colnames(mtcars)[[i]] # Nastavení jmen prvků outputu
}

out

>#        mpg        cyl       disp         hp       drat         wt       qsec 
>#  20.090625   6.187500 230.721875 146.687500   3.596563   3.217250  17.848750 
>#         vs         am       gear       carb 
>#   0.437500   0.406250   3.687500   2.812500

Mohli bychom jej poměrně snadno předělat na while loop. To, co přibude, je i <- 1 jako počáteční hodnota pro první sloupec. V těle while loop je pak stanoveno, že na konci každé iterace se má i zvýšit o 1. A v podmínce i <= ncol(mtcars) je stanoveno, že se while loop má běžet, dokud i je menší nebo rovno počtu sloupců mtcars(); pak se iterace ukončí (kdyby se neukončila, počítali bychom průměry pro další, neexistující sloupce, což by pochopitelně nedávalo smysl).

out <- double(ncol(mtcars))
out

>#  [1] 0 0 0 0 0 0 0 0 0 0 0

i <- 1

while (i <= ncol(mtcars)) {
  out[[i]] <- mean(mtcars[[i]])
  names(out)[[i]] <- colnames(mtcars)[[i]]
  i <- i + 1 
}

out

>#        mpg        cyl       disp         hp       drat         wt       qsec 
>#  20.090625   6.187500 230.721875 146.687500   3.596563   3.217250  17.848750 
>#         vs         am       gear       carb 
>#   0.437500   0.406250   3.687500   2.812500

Takhle bychom mohli pomocí while loop simulovat hody mincí a zjistit, kolik hodů je nutných, aby padla třikrát za sebou hlava (panna). Nejprve si vytvoříme funkci flip(), která náhodně vylosuje z vektoru c("H", "T") jeden prvek.

flip <- function() {
  sample(c("H", "T"), 1)
}
flip() # funkce náhodně vybírá "H" (head) nebo "T" (tail)

># [1] "T"

Pak si nastavíme iniciální hodnoty pro počet hodů (flips) a počet toho, kolikrát doposud padla panna/hlava (heads). Pomocí while loop pak opakovaně spouštíme funkci flip(), po každé iteraci se hodnota flips zváší o 1. Pokud padne hlava, hodnota heads se zvýší o 1, ale pokud padne orel, hodnota heads se vynuluje. Iterace budou pokračovat, dokud nepadne třikrát za sebou panna, tj. dokud platí podmínka heads < 3.

flips <- 0
heads <- 0

while (heads < 3) {
  if (flip() == "H") { # když padne hlava
    heads <- heads + 1 # připočti 1 k objektu heads
  } else {
    heads <- 0 # jinak heads "vynuluj"
  } 
  flips <- flips + 1 # počítej poček pokusů (hodů mincí celkem)
}
flips

># [1] 5

Cvičení

Zkuste místo hodů mincí simulovat hody šestistěnnou kostkou a pomocí while loop zjistit, jak dlouho potrvá (kolik hodů bude potřeba), dokud dvakrát za sebou nepadne šestka.

3.7 Hromadný import souborů

Dejme tomu, že chceme importovat do R více datasetů z SPSS (s příponou .sav)a nechceme to dělat po jednom, ale místo toho využít iteraci.

Nejprve můžeme pomocí funkce dir() vyhledat jména všech souborů s koncovkou .sav v podsložce data/ nebo, přesněji řečeno, uložit si cesty (file paths) ke všem souborům ze složky data, které končí příponou .sav.

files <- dir("data/", pattern = "\\.sav$", full.names = TRUE)
files

># [1] "data/eukids.sav"            "data/international.sav"    
># [3] "data/international_new.sav" "data/lem.sav"              
># [5] "data/rses_data.sav"         "data/titanic.sav"

Dále bychom chtěli, aby jména importovaných objektů byla podobná jako jména importovaných souborů, ale bez uvedení podsložky data a bez přípony .sav; a kdyby se v nějakém názvu vyskytl space (" "), chceme ho nahradit podtržítkem.

nms <- files %>% 
  str_remove_all(pattern = "(data/)|(.sav)") %>% 
  str_replace_all(pattern = " ", replacement = "_")
nms

># [1] "eukids"            "international"     "international_new"
># [4] "lem"               "rses_data"         "titanic"

Dále si jako obvykle připravíme output pro for loop, a to prázný list s tolika prvky, jako je počet importovaných souborů.

data_list <- vector("list", length = length(files)) %>% 
  set_names(nms)

Pak můžeme vytvořit samotný for loop

for (i in seq_along(data_list)) {
  data_list[[i]] <- haven::read_sav(files[[i]])
}

str(data_list, max.level = 1)

># List of 6
>#  $ eukids           : tibble [768 × 236] (S3: tbl_df/tbl/data.frame)
>#  $ international    : tibble [20 × 7] (S3: tbl_df/tbl/data.frame)
>#  $ international_new: tibble [20 × 7] (S3: tbl_df/tbl/data.frame)
>#  $ lem              : tibble [380 × 21] (S3: tbl_df/tbl/data.frame)
>#  $ rses_data        : tibble [11,089 × 13] (S3: tbl_df/tbl/data.frame)
>#  $ titanic          : tibble [891 × 12] (S3: tbl_df/tbl/data.frame)

A nakonec převést prvky listu jako samostatné objekty do global environment pomocí funkce list2env()

list2env(data_list, envir = .GlobalEnv)

># <environment: R_GlobalEnv>

Cvičení

Zkuste hromadně importovat soubory z jiného formátu, např. .csv. Klidně k tomu využijte předchozí kód, pouze ho vhodně pozměňte.

4 Funkce map()

Protože často potřebujeme provést určitou operaci pro několik vektorů, a protože tvorba for loops může být zbytečně pracná, nabízí balíček purrr soubor funkcí, které iteraci značně usnadňují. Tyto funkce se liší v závislosti na očekávaném outputu:

• map() vytvoří list.
• map_lgl() vytvoří logical vector. *map_int() vytvoří integer vector.
• map_dbl() vytvoří double vector.
• map_chr() vytvoří character vector.
• map_df() vytvoří tibble,

Nejdříve si opět vytvoříce cvičný dataframe.

df <- tibble(
  a = rnorm(10),
  b = rnorm(10),
  c = rnorm(10),
  d = rnorm(10)
)

Použití těchto základních funkcí balíčku purrr je poměrně jednoduché, prvním argumentem je vstupní vektor a druhým nejčastěji nějaká funkce.

map_dbl(df, mean)

>#          a          b          c          d 
># -0.2636540 -0.3225047  0.2936220 -0.1147862

map_dbl(df, median)

>#           a           b           c           d 
># -0.03359335 -0.55025607  0.14775946 -0.30473854

map_dbl(df, sd)

>#        a        b        c        d 
># 1.310293 1.215060 1.049160 1.195070

Funkce může být i anonymní, vytvořená adhoc.

# Standardní zápis anonymní funkce
map_dbl(df, function(x) {mean(x)/sd(x)} ) 

>#           a           b           c           d 
># -0.20121757 -0.26542286  0.27986395 -0.09604978

# Zkrácený zápis anonymní funkce, tzv. lambda funkce, uvozená vlnovkou
map_dbl(df, ~mean(.x) /  sd(.x)) 

>#           a           b           c           d 
># -0.20121757 -0.26542286  0.27986395 -0.09604978

Jako druhý argument lze uvést i character vektor či integer vektor. Kódem dole vlastně říkáme, že z každého prvku objektu df (prvkem je každý sloupec matice), chceme vybrat pátý “podprvek” (neboli hodnotu na pátém řádku)

map_dbl(df, 5) # Výběr pátého řádku

>#          a          b          c          d 
># -2.3482903 -0.5362426 -1.1779675 -0.8675661

map_dbl(df, ~.x[[5]]) # Totéž, akorát podrobněji rozepsáno

>#          a          b          c          d 
># -2.3482903 -0.5362426 -1.1779675 -0.8675661

Funkce map_lgl(), map_int(), map_dbl() i map_chr() kontrolují, zda je output požadovaného typu. Očekávají, že výsledek každé operace bude jen jeden prvke daného typu. Proto by např. tyto příkazy nefungovaly a ústily v chybové hlášení:

map_lgl(df, mean) # Vypočtené průměry jsou typu doble, ne logical.

># Error: Can't coerce element 1 from a double to a logical

map_dbl(df, range) # Funkce range nevrací jednu hodnotu, ale dvě (Min a Max)

># Error in `stop_bad_type()`:
># ! Result 1 must be a single double, not a double vector of length 2

Sada funkcí map_ používá dot-dot-dot (...), takže lze přidat jakékoli další argumenty pro funkci použitou v druhém argumentu map_. Níže jsme upravili argumenty trim a na.rm pro funkci mean()

map_dbl(df, mean, trim = 0.5, na.rm = TRUE)

>#           a           b           c           d 
># -0.03359335 -0.55025607  0.14775946 -0.30473854

Funkce map()také usnadňuje práci s více modely. Dejme tomu, že chcete fitovat lineárně regresní model zvlášť pro jednotlivé skupiny. Nejdříve pomocí funkce split() rozdělíme dataframe mtcars na list dataframů podle proměnné cyl (počet válců)

mtcars_split <- mtcars %>% 
  as_tibble() %>% # Není nutné, ale tibbles vypadají v konzoli lépe
  split(.$cyl)
mtcars_split

># $`4`
># # A tibble: 11 × 11
>#      mpg   cyl  disp    hp  drat    wt  qsec    vs    am  gear  carb
>#    <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
>#  1  22.8     4 108      93  3.85  2.32  18.6     1     1     4     1
>#  2  24.4     4 147.     62  3.69  3.19  20       1     0     4     2
>#  3  22.8     4 141.     95  3.92  3.15  22.9     1     0     4     2
>#  4  32.4     4  78.7    66  4.08  2.2   19.5     1     1     4     1
>#  5  30.4     4  75.7    52  4.93  1.62  18.5     1     1     4     2
>#  6  33.9     4  71.1    65  4.22  1.84  19.9     1     1     4     1
>#  7  21.5     4 120.     97  3.7   2.46  20.0     1     0     3     1
>#  8  27.3     4  79      66  4.08  1.94  18.9     1     1     4     1
>#  9  26       4 120.     91  4.43  2.14  16.7     0     1     5     2
># 10  30.4     4  95.1   113  3.77  1.51  16.9     1     1     5     2
># 11  21.4     4 121     109  4.11  2.78  18.6     1     1     4     2
># 
># $`6`
># # A tibble: 7 × 11
>#     mpg   cyl  disp    hp  drat    wt  qsec    vs    am  gear  carb
>#   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
># 1  21       6  160    110  3.9   2.62  16.5     0     1     4     4
># 2  21       6  160    110  3.9   2.88  17.0     0     1     4     4
># 3  21.4     6  258    110  3.08  3.22  19.4     1     0     3     1
># 4  18.1     6  225    105  2.76  3.46  20.2     1     0     3     1
># 5  19.2     6  168.   123  3.92  3.44  18.3     1     0     4     4
># 6  17.8     6  168.   123  3.92  3.44  18.9     1     0     4     4
># 7  19.7     6  145    175  3.62  2.77  15.5     0     1     5     6
># 
># $`8`
># # A tibble: 14 × 11
>#      mpg   cyl  disp    hp  drat    wt  qsec    vs    am  gear  carb
>#    <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
>#  1  18.7     8  360    175  3.15  3.44  17.0     0     0     3     2
>#  2  14.3     8  360    245  3.21  3.57  15.8     0     0     3     4
>#  3  16.4     8  276.   180  3.07  4.07  17.4     0     0     3     3
>#  4  17.3     8  276.   180  3.07  3.73  17.6     0     0     3     3
>#  5  15.2     8  276.   180  3.07  3.78  18       0     0     3     3
>#  6  10.4     8  472    205  2.93  5.25  18.0     0     0     3     4
>#  7  10.4     8  460    215  3     5.42  17.8     0     0     3     4
>#  8  14.7     8  440    230  3.23  5.34  17.4     0     0     3     4
>#  9  15.5     8  318    150  2.76  3.52  16.9     0     0     3     2
># 10  15.2     8  304    150  3.15  3.44  17.3     0     0     3     2
># 11  13.3     8  350    245  3.73  3.84  15.4     0     0     3     4
># 12  19.2     8  400    175  3.08  3.84  17.0     0     0     3     2
># 13  15.8     8  351    264  4.22  3.17  14.5     0     1     5     4
># 14  15       8  301    335  3.54  3.57  14.6     0     1     5     8

Pak pro každý prvek listu (v tomto případě pro každý dataframe) zvlášť fitujeme stejný model. Všimněte si, že .x je placeholder pro jednotlivé prvky mtcars_split() (jednotlivé dataframy)

models <- mtcars_split %>% 
  map(~lm(mpg ~ wt, data = .x))

Dosteneme list s jednotlivými odhadnutými modely

models

># $`4`
># 
># Call:
># lm(formula = mpg ~ wt, data = .x)
># 
># Coefficients:
># (Intercept)           wt  
>#      39.571       -5.647  
># 
># 
># $`6`
># 
># Call:
># lm(formula = mpg ~ wt, data = .x)
># 
># Coefficients:
># (Intercept)           wt  
>#       28.41        -2.78  
># 
># 
># $`8`
># 
># Call:
># lm(formula = mpg ~ wt, data = .x)
># 
># Coefficients:
># (Intercept)           wt  
>#      23.868       -2.192

Pak si např. můžeme extrahovat různé souhrnné statistiky sami.

models %>% 
  map(summary) %>% 
  map_dbl("r.squared") # pomocí character vectoru

>#         4         6         8 
># 0.5086326 0.4645102 0.4229655

Snadnější je ale použít funkce z balíčku broom. Pomocí funkce tidy() tak můžeme extrahovat koeficienty jednotlivých modelů a pomocí funkce glance() informace o fitu modelu. Arugmentem .id = "cyl" jen říkáme, ať se nám do outputu přidá sloupec s názvem "cyl", který bude obsahovat názvy prků vstupního vektoru models(), ať máme jistotu, ke kterému datasetu statistiky náležejí.

names(models)

># [1] "4" "6" "8"

models %>% 
  map_df(broom::tidy, .id = "cyl")

># # A tibble: 6 × 6
>#   cyl   term        estimate std.error statistic    p.value
>#   <chr> <chr>          <dbl>     <dbl>     <dbl>      <dbl>
># 1 4     (Intercept)    39.6      4.35       9.10 0.00000777
># 2 4     wt             -5.65     1.85      -3.05 0.0137    
># 3 6     (Intercept)    28.4      4.18       6.79 0.00105   
># 4 6     wt             -2.78     1.33      -2.08 0.0918    
># 5 8     (Intercept)    23.9      3.01       7.94 0.00000405
># 6 8     wt             -2.19     0.739     -2.97 0.0118

models %>% 
  map_df(broom::glance, .id = "cyl")

># # A tibble: 3 × 13
>#   cyl   r.squared adj.r…¹ sigma stati…² p.value    df logLik   AIC   BIC devia…³
>#   <chr>     <dbl>   <dbl> <dbl>   <dbl>   <dbl> <dbl>  <dbl> <dbl> <dbl>   <dbl>
># 1 4         0.509   0.454  3.33    9.32  0.0137     1 -27.7   61.5  62.7   99.9 
># 2 6         0.465   0.357  1.17    4.34  0.0918     1  -9.83  25.7  25.5    6.79
># 3 8         0.423   0.375  2.02    8.80  0.0118     1 -28.7   63.3  65.2   49.2 
># # … with 2 more variables: df.residual <int>, nobs <int>, and abbreviated
># #   variable names ¹​adj.r.squared, ²​statistic, ³​deviance

Cvičení

Použite jednu z map_ funkcí podle očekávaného typu outputu pro vyřešení následujících úloh.

• Spočítejte průměr každého ze sloupců mtcars.
• Zjistěte typ (typeof()) každého sloupce v datasetu starwars
• Zjistěte počet jedinečných hodnot n_distinct() v každém sloupci iris
• Generujte 10 náhodných hodnot z normálního rozdělení s různými průměry, a to -10, 0, 10 a 100
• vytvořte logický vektor, který pro každý sloupec starwars indikuje, zda se jedná of factor (is.factor())

Co dělají tyto příkazy?

# Příkaz 1
map(1:5, runif)
  
# Příkaz 2
map(-2:2, rnorm, n = 5)

# Příkaz 3
map(c(1, 10, 100), rnorm, n = 8, mean = 0)

Proč tento příkaz nefunguje? (Zkuste n nastavit na 1)

map_dbl(-2:2, rnorm, n = 5) 

># Error in `stop_bad_type()`:
># ! Result 1 must be a single double, not a double vector of length 5

4.1 Funkce map() s více argumenty

Doposud jsme prováděli iterace pouze s jedním vstupním vektorem, ale často máme vstupních vektorů více, např. pokud bychom chtěli generovat data z normální distribuce a měnit jak populační průměr, tak směrodatnou odchylku. Pro podobné případy, kdy máme více vstupních vektorů, slouží funkce map2() a pmap()

Dejme tomu, že chceme měnit jak průměr, tak směrodatnou odchylku normální distribuce, ze které se budou náhodně generovat hodnoty. Připravíme si proto dva vstupní vektory.

mu <- list(5, 10, -3)
sigma <- list(1, 5, 10)

A oba je vložíme do funkce map2 v pořadí, ve kterém je funce rnorm očekává, a zároveň pomocí argumentu n = 5 nastavíme, že chceme vygenerovat právě pět hodnot.

map2(mu, sigma, rnorm, n = 5)

># [[1]]
># [1] 5.736195 6.365777 4.423736 4.195268 4.464934
># 
># [[2]]
># [1] 13.957740  6.461964  3.622283 21.877322  4.532051
># 
># [[3]]
># [1]  -1.075584  -4.261604 -16.884321   1.699022   6.603824

Případně použijeme tzv. lambda funkci ~funkce(), kde se pomocí .x nebo ..1 odkazujeme na první input (mu), a pomocí .y nebo ...2 na druhý input (sigma).

map2(mu, sigma, ~rnorm(5, mean = .x, sd = .y))

># [[1]]
># [1] 4.994859 4.129536 4.537121 3.088849 5.369816
># 
># [[2]]
># [1]  7.688560 11.609779  8.579993 15.459282 10.801082
># 
># [[3]]
># [1]  24.0816828  -0.1371007  -3.3251025 -12.8967802  -0.8819007

map2(mu, sigma, ~rnorm(5, mean = ..1, sd = ..2))

># [[1]]
># [1] 4.027982 5.549224 3.950978 5.449134 4.614144
># 
># [[2]]
># [1]  2.149322 16.134713 13.333472  4.241984  9.196058
># 
># [[3]]
># [1] -12.729138  -9.104617   1.992769   5.793641   6.231292

I když by mohly existovat funkce map3(), map4(), map5(), map6() atd., není to zapotřebí, protože existuje obecná funkce pmap(), která zvládá neomezené množství vstupních vektorů, pokud jsou součástí listu (i dataframes a tibbles jsou v základu lists, takže je můžeme použít také).

Postupujeme tak, že do jednoho listu si uložíme vektory všech vstupních argumentů Nejlépe uděláme, pokud prvky listu nazveme přesně podle názvů vstupních argumentů následně použité funkce, zde rnorm()

args(rnorm)

># function (n, mean = 0, sd = 1) 
># NULL

params <- list(
  n = c(1, 3, 5),
  mean = c(5, 10, -3),
  sd = c(1, 5, 10)
)

A pak list vložíme do funce pmap() a iterativně aplikujeme funkci rnorm(), takže se nám vygeneruje jedna hodnota z normálního rozdělení s průměrem 3 a SD = 1, tři hodnoty z normálního rozdělení s průměrem 10 a SD = 5 a pět hodnot z normálního rozdělení s průměrem -3 a SD = 10

params %>%
  pmap(rnorm)

># [[1]]
># [1] 5.693134
># 
># [[2]]
># [1]  9.302498 -6.553394 12.692329
># 
># [[3]]
># [1]   5.924859   1.894620 -14.963517   1.627379   6.717036

Pro ještě sofistikovanější iterace můžeme (kromě toho, že měníme vstupní argumenty funkce), měnit i funkci samotnou. Nejprve si vytvoříme vektor s názvy funkcí. Funkce runif() generuje hodnoty z uniformní distribuce, rnorm() z normální distribuce a rpois() z Poissonovy distribuce

fns <- c("runif", "rnorm", "rpois")

Pak si vytvoříme vektor se vstupními parametry

params <- list(
  list(min = -1, max = 1), 
  list(mean = 10, sd = 5), 
  list(lambda = 10)
)

A nakonec použijeme funkci invoke_map(), kde jako první argument použijeme názvy funkcí (zapsané jako textový vektor), jako druhý argument list s argumenty, které jsou specifické pro dané funkce a které budou do jednotlivýh funkcí zadány: min a max bude zadáno do první funkce runif(), protože uniformní rozdělení je definováno minimem a maximem, mean a sd do druh0 funkce rnorm, protože normální rozdělení je definováno průměrem a směrodatnou odchylkou a lambda definuje tvar Poissonova rozdělení. Jako další argumenty můžeme uvést argumenty, které mají všechny funkce stejné a chceme je nechat konstatní (zde n = 5 pro vygeneraování pěti hodnot z každého rozdělení)

invoke_map(fns, params, n = 5)

># [[1]]
># [1] -0.85902204 -0.08042961  0.40317024 -0.82611966  0.98588803
># 
># [[2]]
># [1]  6.676153 12.427300  8.121986  7.190618  8.280414
># 
># [[3]]
># [1] 10 15  9 13 11

Také si můžeme názvy funkcí uložit do jednoho tibble, Funkce tribble je podobná, jako funkce tibble(), akorát hodnoty zadáváme po řádcích a jména sloupců definujeme jako prnví pomocí tildy. Všimněte si, že sloupce v tibbles nemusejí být atomické vektory, např. níže je params typu list a obsahuje jako prvky další listy.

sim <- tribble(
  ~fns,      ~params,
  "runif", list(min = -1, max = 1),
  "rnorm", list(sd = 5),
  "rpois", list(lambda = 10)
)

sim

># # A tibble: 3 × 2
>#   fns   params          
>#   <chr> <list>          
># 1 runif <named list [2]>
># 2 rnorm <named list [1]>
># 3 rpois <named list [1]>

Abychom měli názvy požitých funkcí, argumenty vložené do použitých funkcí i vygenerované hodnoty pohromadě v jednom tibbles, můžeme si do sim doplnit pomocí mutate nový sloupes s vygenerovanými hodnotami.

sim <- sim %>% 
  mutate(simulated_values = invoke_map(fns, params, n = 5))
sim

># # A tibble: 3 × 3
>#   fns   params           simulated_values
>#   <chr> <list>           <list>          
># 1 runif <named list [2]> <dbl [5]>       
># 2 rnorm <named list [1]> <dbl [5]>       
># 3 rpois <named list [1]> <int [5]>

sim$simulated_values

># [[1]]
># [1] -0.09047674  0.54040952 -0.87470000  0.63016299 -0.39771495
># 
># [[2]]
># [1] -1.730004 -8.910285  2.324709 -9.755153 -2.580567
># 
># [[3]]
># [1] 11 12 12 13  8

Cvičení

Vygenerujte různý počet hodnoty z normálního rozdělení podle tabulky:

n	M	SD
30	2.5	0.50
60	3.0	0.75
90	5.0	1.00

5 Funkce walk()

Funkce walk() je alternativou k funkci map(), pokud chceme spusti nějakou funkci kvůli jejím “vedlejším efektům”, jako je tvorba grafu nebo uložení souboru na disk, nikoli kvůli jejím výstupním hodnotám a tvorbě nového datového objektu.

Tady je jednoduchý příklad s funkcí print()

x <- list(1, "a", 3)

x %>% 
  walk(print)

># [1] 1
># [1] "a"
># [1] 3

Pomocí pwalk() např. můžeme uložit několik grafů zároveň. Nejdříve si vytvoříme list s bodovými grafy proměnných mpg a wt pro jednotlivé úrovně cyl

plots <- mtcars %>% 
  split(.$cyl) %>% 
  map(~ggplot(., aes(wt, mpg)) + geom_point())

Pak si připravíme jména souborů, které chceme uložit

plot_names <- str_c("cyl", names(plots), ".png")
plot_names

># [1] "cyl4.png" "cyl6.png" "cyl8.png"

Sloučíme obojí od jednoho listu, jehož prvky pojmenujeme podle názvů argumentů funkce ggsave(), Tento list pak vložíme do funkce pwalk(), čímt iterativně uplatníme funkci ggsave(). Dodatečným argumentem path pak můžeme specifikovat, do které podsložky se mají všechny grafy uložit.

list(
  filename = plot_names, # ggsave má tyto dva hlavní argumenty: filename a plot
  plot = plots
) %>% 
  pwalk(ggsave, path = "plots")

># Saving 7 x 4.32 in image
># Saving 7 x 4.32 in image
># Saving 7 x 4.32 in image

Cvičení

Balíček ggplot2 obsahuje také dataset msleep s informace i spánku různých savců. Zkuste pomocí nejprve pomocí funkce map() hromadně vytvořit boxploty všech kvantitativních proměnných (sleep_total až bodywt) v závislosti na druhu potravy (proměnná vore, na osu X) a pak se všechny grafy pokuste uložit pomocí funkce pwalk() nebo walk2().

6 Predikativní funkce

Dále balíček purrr nabízí soubor funkcí, které využívají jiné, predikativní funkce, která vracejí jako výstupní hodnoty buď TRUE, nebo FALSE. Funkce keep() zachová prvky, pro které je predikát TRUE, funkce discard() naopak zahodí prvky, pro které je predikát TRUE. Níže můžete vidět příklady použití těchto funkcí.

# Zachovat v datasetu iris pouze sloupce typu factor
iris %>% 
  keep(is.factor) %>% 
  glimpse()

># Rows: 150
># Columns: 1
># $ Species <fct> setosa, setosa, setosa, setosa, setosa, setosa, setosa, setosa…

# Zachovat v datasetu iris pouze sloupce, které nejsou typu factor
iris %>% 
  discard(is.factor) %>% 
  glimpse()

># Rows: 150
># Columns: 4
># $ Sepal.Length <dbl> 5.1, 4.9, 4.7, 4.6, 5.0, 5.4, 4.6, 5.0, 4.4, 4.9, 5.4, 4.…
># $ Sepal.Width  <dbl> 3.5, 3.0, 3.2, 3.1, 3.6, 3.9, 3.4, 3.4, 2.9, 3.1, 3.7, 3.…
># $ Petal.Length <dbl> 1.4, 1.4, 1.3, 1.5, 1.4, 1.7, 1.4, 1.5, 1.4, 1.5, 1.5, 1.…
># $ Petal.Width  <dbl> 0.2, 0.2, 0.2, 0.2, 0.2, 0.4, 0.3, 0.2, 0.2, 0.1, 0.2, 0.…

Funkce some() ověřuje, zda je nějaký predikát platný aspoň pro jeden prvek , a funkce every() ověřuje, zda je nějaký predikát platný pro všechny prvky.

x <- list(1:5, letters, list(10))
x

># [[1]]
># [1] 1 2 3 4 5
># 
># [[2]]
>#  [1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l" "m" "n" "o" "p" "q" "r" "s"
># [20] "t" "u" "v" "w" "x" "y" "z"
># 
># [[3]]
># [[3]][[1]]
># [1] 10

# Je aspoň jeden prvek typu character?
x %>% 
  some(is.character) 

># [1] TRUE

# Jsou všechny prvky numerické?
x %>% 
  every(is.numeric) 

># [1] FALSE

Cvičení

Zkuste ověřit, zda v datasetu starwars:

• je nějaký sloupec typu list;
• je nějaký sloupec typu factor;
• jsou všechny sloupce atomické vektory.

Vyberte z datasetu starwars:

• pouze sloupce typu character;
• pouze sloupce typu double;
• pouze sloupce typu integer;
• pouze sloupce typu list;
• pouze sloupce, které jsou atomickými vektory.

Vyřaďte z datasetu starwars:

• všechny sloupce typu character;
• všechny sloupce typu double;
• všechny sloupce typu integer;
• všechny sloupce typu list.