Rのこと。

記事は引っ越し作業中。2023年中までに引っ越しを完了させてブログは削除予定

Advanced R: Functionals

はじめに

Advanced R 2nd EditionのPaperBook版が届いたので、Rについてのおさらい備忘録。また、書籍を読んでいるにも関わらず誤った解釈や用語の誤用もあるかもしれませんので、参考にする際は自己責任でお願いします。

この記事のライセンスはAdvanced R 2nd Editionと同様で下記の通りです。

Functionals

Functionalsでは、汎関数の一般的な使い方、特にmap()への理解について説明されています。

map()

purrr::map()はベクトルと関数を取り、ベクトルの各要素に対して関数を1回呼び出し、結果をリストに返す関数。つまり、map(1:3, f)とlist(f(1), f(2), f(3))は同等。

triple <- function(x) x * 3
map(1:3, triple)
[[1]]
[1] 3

[[2]]
[1] 6

[[3]]
[1] 9

forループで書き直すと下記の通りです。

x <- 1:3
out <- vector("list", length(x))
for (i in seq_along(x)) {
  out[[i]] <- 3 * x[[i]]
}
out

[[1]]
[1] 3

[[2]]
[1] 6

[[3]]
[1] 9

つまり、簡易版map()のコードは下記の通り。入力と同じ長さのリストを割り当ててから、forループでリストを埋める。実際にはC言語で実装されているみたいなので、下記より高速。

simple_map <- function(x, f, ...) {
  out <- vector("list", length(x))
  for (i in seq_along(x)) {
    out[[i]] <- f(x[[i]], ...)
  }
  out
}

アトミックベクトルとmap()

map()にはリストではなくベクトルを返すmap_lgl()、map_int()、map_dbl()、map_chr()が用意される。つまり、指定された型のアトミックベクトルを返す。

map_chr(mtcars, typeof)
      mpg       cyl      disp        hp      drat        wt 
 "double"  "double"  "double"  "double"  "double"  "double" 
     qsec        vs        am      gear      carb 
 "double"  "double" "integer"  "double"  "double" 

map_lgl(mtcars, is.double)
  mpg   cyl  disp    hp  drat    wt  qsec    vs    am  gear  carb 
 TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE 

n_unique <- function(x) length(unique(x))
map_int(mtcars, n_unique)
 mpg  cyl disp   hp drat   wt qsec   vs   am gear carb 
  25    3   27   22   22   29   30    2    2    3    6

つまり、入力と同じ長さの指定された型のベクトルを割り当ててから、forループでベクトルを埋める。

out <- vector("character", ncol(mtcars))
names(out) <- names(mtcars)
for (i in names(mtcars)) {
  out[[i]] <- typeof(mtcars[[i]])
}
out
      mpg       cyl      disp        hp      drat        wt 
 "double"  "double"  "double"  "double"  "double"  "double" 
     qsec        vs        am      gear      carb 
 "double"  "double" "integer"  "double"  "double" 

out <- vector("logical", ncol(mtcars))
names(out) <- names(mtcars)
for (i in names(mtcars)) {
  out[[i]] <- is.double(mtcars[[i]])
}
out
  mpg   cyl  disp    hp  drat    wt  qsec    vs    am  gear  carb 
 TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE 

out <- vector("integer", ncol(mtcars))
names(out) <- names(mtcars)
for (i in names(mtcars)) {
  out[[i]] <- length(unique((mtcars[[i]])))
}
out
 mpg  cyl disp   hp drat   wt qsec   vs   am gear carb 
  25    3   27   22   22   29   30    2    2    3    6 

# "am" is class factor 
mtcars <- mtcars %>% mutate(am = as.numeric(am)-1)
out <- vector("double", ncol(mtcars))
names(out) <- names(mtcars)
for (i in names(mtcars)) {
  out[[i]] <- mean(mtcars[[i]])
}
out
      mpg        cyl       disp         hp       drat         wt 
 20.090625   6.187500   3.780862 146.687500   3.596563   3.217250 
      qsec         vs         am       gear       carb 
 17.848750   0.437500   0.406250   3.687500   2.812500

無名関数とショートカット

map()は既存の関数を使用する代わりに、インラインの無名関数を作成できる。#1が基本的な記述で、#2が無名関数を使った例、#3はラムダ式と呼ばれるショートカットうぃ使った記述。

# 01
n_unique <- function(x) length(unique(x))
map_dbl(mtcars, n_unique)

# 02
map_dbl(mtcars, function(x) length(unique(x)))

# 03
map_dbl(mtcars, ~ length(unique(.x)))

ショートカットの裏側は下記の通り。

as_mapper(~ length(unique(.x)))
<lambda>
function (..., .x = ..1, .y = ..2, . = ..1) 
length(unique(.x))
attr(,"class")
[1] "rlang_lambda_function"

as_mapper()は与えられた文字、数字またはリストなどの入力に対する抽出関数を作成する。

as_mapper(c(1, 2))
function (x, ...) 
pluck(x, 1, 2, .default = NULL)
<environment: 0x7fe1843125e8>

as_mapper(c("a", "b"))
function (x, ...) 
pluck(x, "a", "b", .default = NULL)
<environment: 0x7fe18174ebd8>

as_mapper(list(1, "b"))
function (x, ...) 
pluck(x, 1, "b", .default = NULL)
<environment: 0x7fe18431f158>

purrr::pluck()というベクトルから要素を抽出するためのショートカットもある。JSONなどを触るときに便利。

x <- list(
  list(-1, x = 1, y = c(2), z = "a"),
  list(-2, x = 4, y = c(5, 6), z = "b"),
  list(-3, x = 8, y = c(9, 10, 11))
)

map_dbl(x, "x")
[1] 1 4 8

map_dbl(x, 1)
[1] -1 -2 -3

map_dbl(x, list("y", 1))
[1] 2 5 9
 
map_chr(x, "z")
 エラー: Result 3 must be a single string, not NULL of length 0
Call `rlang::last_error()` to see a backtrace

map_chr(x, "z", .default = NA)
[1] "a" "b" NA 
>

...(dot-dot-dot)

map()が呼び出している関数に追加の引数を渡すことももちろん可能。

x <- list(1:5, c(1:10, NA))
map_dbl(x, ~ mean(.x, na.rm = TRUE))
[1] 3.0 5.5

map_dbl(x, mean, na.rm = TRUE)
[1] 3.0 5.5

map()において、「引数を渡すこと」と「無名関数の中に余分な引数を置くこと」を比較すると、微妙な違いがある。無名関数f()に引数を入れるというのは、呼び出されるときに1度だけではなく、実行のたびに評価されることを意味する。

plus <- function(x, y) x + y
x <- c(0, 0, 0, 0)

map_dbl(x, plus, runif(1))
[1] 0.7508426 0.7508426 0.7508426 0.7508426

map_dbl(x, ~ plus(.x, runif(1)))
[1] 0.3837969 0.4666181 0.9763213 0.1613520

この仕組を利用すればシュミレーションを簡単にできる。

trials <- map(1:10000, ~ t.test(rpois(10, 10), rpois(10, 7)))
trials_df <- tibble(p_value = map_dbl(trials, "p.value"))
head(trials_df,5)
# A tibble: 5 x 1
   p_value
     <dbl>
1 0.619   
2 0.0259  
3 0.124   
4 0.000418
5 0.000883

# trials_df %>% 
#   ggplot(aes(x = p_value, fill = p_value < 0.05)) + 
#   geom_histogram(binwidth = .01)

パースタイル

現実的な問題を解決するために複数の{purrr}を使用する方法についてみていく。下記では、mtcarsのcylごとに3つのデータフレームのリストを作成し、各データフレームに線形モデルを当てはめ、回帰係数を抽出している。

by_cyl <- split(mtcars, mtcars$cyl)
by_cyl %>% 
  map(~ lm(mpg ~ wt, data = .x)) %>% 
  map(coef) %>% 
  map_dbl(2)

        4         6         8 
-5.647025 -2.780106 -2.192438

紹介していなかったがmap()はlapply()というapplyファミリーの互換関数なので、書き直すことももちろんできる。

by_cyl %>% 
  lapply(function(data) lm(mpg ~ wt, data = data)) %>% 
  lapply(coef) %>% 
  vapply(function(x) x[[2]], double(1))
        4         6         8 
-5.647025 -2.780106 -2.192438

もちろんforループでも書き直すことは可能。

by_cyl <- split(mtcars, mtcars$cyl)
slope <- double(length(by_cyl))
for (i in seq_along(by_cyl)) {
  model <- lm(mpg ~ wt, data = by_cyl[[i]])
  slope[[i]] <- coef(model)[[2]]
}
slope
[1] -5.647025 -2.780106 -2.192438

map()を使えば、異なる組み合わせのモデル式を一度にモデリングできる。

formulas <- list(
  mpg ~ disp,
  mpg ~ I(1 / disp),
  mpg ~ disp + wt,
  mpg ~ I(1 / disp) + wt
)

models <- map(formulas, lm, data = mtcars)
map(models, coef)
[[1]]
(Intercept)        disp 
  29.599855   -2.515095 

[[2]]
(Intercept)   I(1/disp) 
   10.75202    25.52576 

[[3]]
(Intercept)        disp          wt 
  34.960554   -1.081628   -3.350825 

[[4]]
(Intercept)   I(1/disp)          wt 
  19.024206   18.723245   -1.797649

ブートストラップ回帰分析をシンプルに実装できる。

bootstrap <- function(df) {
  df[sample(nrow(df), replace = TRUE), , drop = FALSE]
}

bootstraps <- map(1:10, ~ bootstrap(mtcars))

bootstraps %>% 
    map(~ lm(mpg ~ disp, data = .x)) %>% 
    map(summary) %>% 
    map_dbl("r.squared")

 [1] 0.7976214 0.8573429 0.7347270 0.7233033
 [5] 0.6881601 0.7229928 0.7138933 0.8491383
 [9] 0.6885969 0.7887907

マップの変形

{map}ファミリは直交する入力と出力の関係を持つ。つまり、入力を行に、出力を列にして、すべてのファミリを行列にまとめることが可能。

リスト アトミックベクトル 同型 出力なし
1つの引数 map() map_lgl()… modify() walk()
2つの引数 map2() map2_lgl()… modify2() walk2()
1つの引数+インデックス imap() imap_lgl()… imodify() iwalk()
n個の引数 pmap() pmap_lgl()… pwalk()

modify()は常に入力と同じタイプの出力を返す。下記のようにデータフレームが入れば、データフレームをアウトする。

two <- function(x) 2*x
> map(df, two)
$x
[1] 2 4 6

$y
[1] 12 10  8

modify(df, two)
  x  y
1 2 12
2 4 10
3 6  8

modify()の簡易版実装は下記のようになる。

simple_modify <- function(x, f, ...) {
  for (i in seq_along(x)) {
    x[[i]] <- f(x[[i]], ...)
  }
  x
}

simple_modify(df, d)
  x  y
1 2 12
2 4 10
3 6  8

map2()は入力に2つのベクトルを用いることが可能。下記のように、加重平均は2つのベクトルの入力を必要とする。わかりやすいように引数名をつけた。

xs <- map(1:2, ~ runif(10))
xs
[[1]]
 [1] 0.53475331 0.36807929 0.27609387 0.69554501 0.09510626
 [6] 0.48716753 0.31880390 0.13628335 0.12398587 0.57005658

[[2]]
 [1] 0.03994327 0.19862998 0.71279430 0.69998526 0.75600462
 [6] 0.84803090 0.54106060 0.04433105 0.30520698 0.49661419

ws <- map(1:2, ~ rpois(10, 5) + 1)
ws
[[1]]
 [1]  7  5 10  6  7  8  4  8  8 11

[[2]]
 [1] 4 6 8 6 6 6 8 4 6 6

map2_dbl(.x = xs, .y = ws, .f = weighted.mean)
[1] 0.3609318 0.5032461

map2()の基本的な実装はmap()と似ている。1つのベクトルを反復する代わりに、2つのベクトルを並列に反復する。3つ以上を扱うpmap()はここでは扱わない。

simple_map2 <- function(x, y, f, ...) {
  out <- vector("list", length(x))
  for (i in seq_along(x)) {
    out[[i]] <- f(x[[i]], y[[i]], ...)
  }
  out
}

simple_map2(x = xs, y = ws, f = weighted.mean)
[[1]]
[1] 0.3609318

[[2]]
[1] 0.5032461

walk()は出力をしない関数。いくつかの関数は、主にその副作用のために呼ばれている。例えば、cat()とmap()を使った例はわかりよい。この場合、リストの部分は表示したくないので、map()はcat()を副作用(文字列を表示するだけ)として呼ばれている。

welcome <- function(x) {
  cat("Welcome ", x, "!\n", sep = "")
}
names <- c("Hadley", "Jenny")

map(names, welcome)
Welcome Hadley!
Welcome Jenny!
[[1]]
NULL

[[2]]
NULL

walk()はこのような場合に、戻り値を無視する。

walk(names, welcome)
Welcome Hadley!
Welcome Jenny!

walk2()はディスクにデータを保存する際に、オブジェクトと保存先のパスを指定できるので便利。

temp <- tempfile()
dir.create(temp)

cyls <- split(mtcars, mtcars$cyl)
paths <- file.path(temp, paste0("cyl-", names(cyls), ".csv"))
walk2(cyls, paths, write.csv)

dir(temp)
[1] "cyl-4.csv" "cyl-6.csv" "cyl-8.csv"

値とインデックスの繰り返し

forループを使用してベクトルをループ処理する基本的な方法は3つある。最初の形式はmap()と似てる。2番目と3番目の形式は、ベクトルの値とインデックスを並列に反復処理できるimap()と同等。

  • 要素をループ：for (x in xs)
  • 数値インデックスをループ：for (i in seq_along(xs))
  • 名前をループ：for (nm in names(xs))

ベクトルが無名の場合、2番目の引数はインデックスになる。そうではない場合、カラム名を引っ張る。

imap_chr(iris, ~ paste0("The first value of ", .y, " is ", .x[[1]]))
                            Sepal.Length 
"The first value of Sepal.Length is 5.1" 
                             Sepal.Width 
 "The first value of Sepal.Width is 3.5" 
                            Petal.Length 
"The first value of Petal.Length is 1.4" 
                             Petal.Width 
 "The first value of Petal.Width is 0.2" 
                                 Species 
  "The first value of Species is setosa" 

x <- map(1:6, ~ sample(1000, 10))
imap_chr(x, ~ paste0("The highest value of ", .y, " is ", max(.x)))
[1] "The highest value of 1 is 761"
[2] "The highest value of 2 is 983"
[3] "The highest value of 3 is 993"
[4] "The highest value of 4 is 991"
[5] "The highest value of 5 is 834"
[6] "The highest value of 6 is 922"

reduce()は長さnのベクトルを取り、1度に1対の値で関数を呼び出すことにより、長さ1のベクトルを生成する。reduce(1:4, f)はf(f(f(1, 2), 3), 4)と等価。リストの中の各ベクトルで一致する値を抽出したい場合、reduce()は役立つ。考え方としては、非常に大きなデータセットの処理によく使われるmap-reduceフレームワーク。

l <- map(1:4, ~ sample(1:10, 15, replace = T))
l
[[1]]
 [1]  7  9  5  7  4  9  4  2 10  6  9  4  1  3  1

[[2]]
 [1]  6  5  1  7  1  3  1  1  4 10  8  4  9  6  1

[[3]]
 [1] 10  6  8  1  6  3  9 10 10  5  2  1  9  8  5

[[4]]
 [1]  8  5  1  3  3  4  1  1  5  8  3  9  1  8 10

reduce(l, intersect)
[1]  9  5 10  1  3

reduce(l, union)
 [1]  7  9  5  4  2 10  6  1  3  8

reduce(l, sum)
[1] 311

sum(unlist(l))
[1] 311

reduce()はforループを囲む単純なラッパーを減らすことが可能。

simple_reduce <- function(x, f) {
  out <- x[[1]]
  for (i in seq(2, length(x))) {
    out <- f(out, x[[i]])
  }
  out
}

挙動を確認してみる。

simple_reduce2 <- function(x, f) {
  out <- x[[1]]
  for (i in seq(2, length(x))) {
    out <- f(out, x[[i]])
    print(out)
  }
}

l <- map(1:100, ~ sample(1:10, 15, replace = T))
simple_reduce2(l, sum)
[1] 132
[1] 202
[略]
[1] 8092
[1] 8166

述語汎関数

  • every()：すべての要素が一致する場合、every(.x、.p)はTRUEを返す。
  • some()：最初のTRUEを見つけたときにTRUEを返し、every()は最初のFALSEを見つけたときにFALSEを返す。
  • detect(.x、.p)：最初に一致した値を返す。
  • detect_index(.x、.p)：最初に一致したインデックスを返す。
  • keep(.x、.p)：一致するすべての要素を保持する。
  • discard(.x、.p)：一致するすべての要素を削除する。
every(df, is.factor)
[1] FALSE

some(df, is.factor)
[1] TRUE

detect(df, is.factor)
[1] a b c
Levels: a b c

detect_index(df, is.factor)
[1] 2

keep(df, is.factor)
  y z
1 a z
2 b z
3 c z

discard(df, is.factor)
  x
1 1
2 2
3 3

mapを変形させる

map()に述語関数を取る変種を作成すると非常に強力です。

df <- data.frame(
  num1 = c(0, 10, 20),
  num2 = c(5, 6, 7),
  chr1 = c("a", "b", "c"),
  stringsAsFactors = FALSE
)

df
  num1 num2 chr1
1    0    5    a
2   10    6    b
3   20    7    c

map_if(df, is.numeric, mean)
$num1
[1] 10

$num2
[1] 6

$chr1
[1] "a" "b" "c"

modify_if(df, is.numeric, mean)
  num1 num2 chr1
1   10    6    a
2   10    6    b
3   10    6    c

map(keep(df, is.numeric), mean)
$num1
[1] 10

$num2
[1] 6

modify(keep(df, is.numeric), mean)
  num1 num2
1   10    6
2   10    6
3   10    6