Rのこと。

R言語のことをまとめる備忘録

【和訳】Programming with dplyr

はじめに

この記事は、「Programming with dplyr 」というVignetteを和訳したものです。VignetteはMITライセンスで公開されているため、この記事もそれに準じます。

dplyrを使ったプログラミング

ほとんどのdplyr関数は非標準評価(NSE)を使用します。これは包括的な用語(catch-all term)で、通常のRの評価規則に従わないことを意味します。代わりに、入力した表現式を捕捉し、それを独自の方法で評価します。これはdplyrコードにとって2つの大きな利点があります。

  • データフレームの名前を繰り返す必要がないため、データフレームに対する操作を簡潔に表現できます。例えば、df[df$x == 1 & df$y ==2 & df$z == 3, ]の代わりにfilter(df, x == 1, y == 2, z == 3)と書くことができます。
  • dplyrは、baseのRとは異なる方法で結果を計算することを選択できます。これは、dplyr自体は機能しませんが、データベースに何をすべきかを指示するSQLを生成するデータベースのバックエンドにとって重要です。

残念ながら、これらの利点はそう簡単には手に入りません。主な欠点が2つあります。

dplyrのほとんどの引数は、参照透明性(referentially transparent)がありません。つまり、他の場所で定義した、見かけ上、同等のオブジェクトで値を置き換えることはできません。言い換えれば、このコードです。

df <- tibble(x = 1:3, y = 3:1)
filter(df, x == 1)
#> # A tibble: 1 x 2
#>       x     y
#>   <int> <int>
#> 1     1     3

このコードと同等ではありません:

my_var <- x
#> Error in eval(expr, envir, enclos): object 'x' not found
filter(df, my_var == 1)
#> Error: object 'my_var' not found

このコードも。

my_var <- "x"
filter(df, my_var == 1)
#> # A tibble: 0 x 2
#> # … with 2 variables: x <int>, y <int>

これは、dplyrの動詞(verbs)の計算方法を変更する引数を持つ関数を作成することを困難にします。

  • dplyrのコードは曖昧です。どの変数が、どこで定義されているかに応じてfilter(df, x == y)は、以下のいずれかと同等になる可能性があります。
df[df$x == df$y, ]
df[df$x == y, ]
df[x == df$y, ]
df[x == y, ]

これはインタラクティブに作業するときには便利ですが(入力を節約してすぐに問題を発見できるため)、あなたが考えるよりも、関数が予測不可能になります。

幸いにも、dplyrにはこれらの課題を克服するためのツールが提供されています。それらはもう少しタイピングを必要としますが、それらは長期的に時間を節約できるので、多少の先行作業も価値があります。

このビネットには2つの目標があります。

  • dplyrのpronouns(代名詞)quasiquotationの使い方を示すことで、データ分析におけるコードの重複を減らし、信頼性の高い関数を記述すること。

  • 表現式と環境の両方を格納するデータ構造であるquosure tidyevalを含む基本的な理論を教え、そして基本的なツールキットを提供します。

ウォームアップから始め、この問題をより慣れ親しんだものに結び付けてから、いくつかの実用的なツールに進み、そこからより深い理論に飛び込みます。

ウォームアップ

あなたはそれを理解していなかったかもしれませんが、あなたはもう一つの分野(ストリング)でこのタイプの問題を解決することを既に達成しています。この関数があなたが望むように動かないことは明らかです。

greet <- function(name) {
  "How do you do, name?"
}
greet("Hadley")
#> [1] "How do you do, name?"

"は、入力をクオートするためです。入力した内容を解釈するのではなく、単に文字列に格納するだけです。関数を望むように動かす1つの方法はpaste()で文字列を少しずつ組み立てることです。

greet <- function(name) {
  paste0("How do you do, ", name, "?")
}
greet("Hadley")
#> [1] "How do you do, Hadley?"

別のアプローチはglueパッケージを使うことです。Rの表現式の値を伴い、文字列を置き換えて、文字列の構成要素をアンクオートすることを可能にします。name引数が、{name}で置き換えられるので、関数の洗練された実装を可能にします。

greet <- function(name) {
  glue::glue("How do you do, {name}?")
}
greet("Hadley")
#> How do you do, Hadley?

プログラミングレシピ

次のレシピでは、dplyrコードの重複を減らすという目的のもと、tidyevalの基本について説明します。これらの例は、理解しやすいように単純化しているため、いくらか本物ではありません。

それらはとても単純なので、なぜ関数を書くのを面倒にするのか、疑問に思うかもしれません。しかし、簡単な例でアイデアを習得することをお勧めします。そうすれば、自分自身のコードで見る複雑な状況に、それらを適用する準備が整うからです。

異なるデータセット

あなたはすでにdplyr verbsの最初の引数、すなわちデータを扱う関数を書く方法を知っています。これは、dplyrがその引数に対して特別なことをしないため、参照的に透過的だからです。たとえば、次のようなコードが繰り返し表示されたとします。

mutate(df1, y = a + x)
mutate(df2, y = a + x)
mutate(df3, y = a + x)
mutate(df4, y = a + x)

あなたは、すでにその重複を捉えるための関数を書くことができます:

mutate_y <- function(df) {
  mutate(df, y = a + x)
}

残念ながら、この単純なアプローチには欠点があります。変数の1つがデータフレームには存在せず、グローバル環境に存在する場合、暗黙的に失敗する可能性があります。

df1 <- tibble(x = 1:3)
a <- 10
mutate_y(df1)
#> # A tibble: 3 x 2
#>       x     y
#>   <int> <dbl>
#> 1     1    11
#> 2     2    12
#> 3     3    13

その曖昧さは、明示的に.data代名詞を使用することによって解決できます。変数が存在しない場合、これは有益なエラーを投げます。

mutate_y <- function(df) {
  mutate(df, y = .data$a + .data$x)
}

mutate_y(df1)
#> Column `a` not found in `.data`

この関数がパッケージに含まれている場合、.dataを使用するとR CMD checkが未定義のグローバル変数についてのNOTEも表示されなくなります(rlang::.dataを一緒にインポートした場合:@importFrom rlang .data)。

異なる表現

引数の1つを、変数名(like x)または表現式(like x + y)にしたい場合、関数を書くのは困難です。これは、dplyrがこれらの入力を自動的にクオートするため、参照として透過的ではないためです。簡単な場合から始めましょう。データを要約するために、グループ化変数を変えたいとします。

df <- tibble(
  g1 = c(1, 1, 2, 2, 2),
  g2 = c(1, 2, 1, 2, 1),
  a = sample(5),
  b = sample(5)
)

df %>%
  group_by(g1) %>%
  summarise(a = mean(a))
#> # A tibble: 2 x 2
#>      g1     a
#>   <dbl> <dbl>
#> 1     1  4   
#> 2     2  2.33

df %>%
  group_by(g2) %>%
  summarise(a = mean(a))
#> # A tibble: 2 x 2
#>      g2     a
#>   <dbl> <dbl>
#> 1     1  3.67
#> 2     2  2

これがうまくいくことを願います。

my_summarise <- function(df, group_var) {
  df %>%
    group_by(group_var) %>%
    summarise(a = mean(a))
}

my_summarise(df, g1)
#> Error: Column `group_var` is unknown

しかし、うまくいきません。おそらく、変数名を文字列として提供することで、問題は解決するかもしれません。

my_summarise(df, "g2")
#> Error: Column `group_var` is unknown

そうはいきません。

エラーメッセージをよく見ると、どちらの場合も同じです。group_by()"のように動作します。その入力を評価せす、クオートします。

この機能を動作させるには、2つのことをする必要があります。私たちは自分自身で入力をクオートする必要があり、group_by()が、その入力をクオートしないように指示する必要があります(なぜならクオートしたので)。

入力をどのようにクオートすれば良いでしょうか。""は入力を引用符で囲むため、使用することはできません。文字列が得られるからです。 代わりに、表現式とその環境を捕捉する関数が必要です(後でこれが重要になる理由を説明します)。baseのRで使うことができる2つの可能なオプション、quote()~があります。これらはどちらも私たちが望むようには動作しないので、quo()という新しい関数が必要です。

quo()"のように動作します。評価するのではなく、その入力をクオートします。

quo(g1)
#> <quosure>
#> expr: ^g1
#> env:  global
quo(a + b + c)
#> <quosure>
#> expr: ^a + b + c
#> env:  global
quo("a")
#> <quosure>
#> expr: ^"a"
#> env:  empty

quo()quosureという特殊なオブジェクトを返します。後でquosureについて学びます。この表現式を捉えたので、それをどのようにgroup_by()で使用すればよいでしょうか。単純な方法でそれを推し進めるだけではうまくいきません。

my_summarise(df, quo(g1))
#> Error: Column `group_var` is unknown

これでは前と同じエラーを得ます。それは、group_by()にクオートをどう扱うのか、まだ指示をしていないからです。言い換えると、my_summarise()によって、プレクオート(pre-quoted)されたので、その入力をクオートしないようにgroup_by()に指示する必要があります。 同じことを言い表すもう一つの方法は、私たちがgroup_varをアンクオートしたいということです。

dplyr(そして一般的にはtidyeval)では、!!を使うことで、評価されるように入力をアンクオートすることが可能です。これにより、実際に必要な動作を実行する関数が得られます。

my_summarise <- function(df, group_var) {
  df %>%
    group_by(!! group_var) %>%
    summarise(a = mean(a))
}

my_summarise(df, quo(g1))
#> # A tibble: 2 x 2
#>      g1     a
#>   <dbl> <dbl>
#> 1     1  4   
#> 2     2  2.33

Huzzah!

まだ1つのステップが残っています。group_by()を呼ぶように、この関数を呼びたいということです。

my_summarise(df, g1)

g1は、呼び出されるオブジェクトがないため、これは機能しません。何を入力したのかを把握し、それをクオートする必要があります。quo()を、それをするために使ってみるとよいかもしれません。

my_summarise <- function(df, group_var) {
  quo_group_var <- quo(group_var)
  print(quo_group_var)

  df %>%
    group_by(!! quo_group_var) %>%
    summarise(a = mean(a))
}

my_summarise(df, g1)
#> <quosure>
#> expr: ^group_var
#> env:  0x70b4c60
#> Error: Column `group_var` is unknown

何が問題になっているのかを明確にするため、print()の呼び出しを追加しました。quo(group_var)~group_varを常に返します。文字通りすぎます。ユーザーが指定した値である~g1を返すように値に置き換えます。

文字列と同様に、""は必要ありません。代わりに、引数を文字列に変換する関数が必要です。それがenquo()です。enquo()は、ダークマジックを使用して、引数を調べ、ユーザーが入力した内容を確認して、その値をクオートして返します。(技術的には、関数引数はpromiseと呼ばれる特別なデータ構造を使用して遅延評価されるので、これはうまくいきます。)

my_summarise <- function(df, group_var) {
  group_var <- enquo(group_var)
  print(group_var)

  df %>%
    group_by(!! group_var) %>%
    summarise(a = mean(a))
}

my_summarise(df, g1)
#> <quosure>
#> expr: ^g1
#> env:  global
#> # A tibble: 2 x 2
#>      g1     a
#>   <dbl> <dbl>
#> 1     1  4   
#> 2     2  2.33

(あなたがRのbaseにおけるquote()substitute()に精通している場合、quo()quote()は等価であり、enquo()substitute()は等価であるとわかるはずです。)

複数のグループ化変数を処理するために、これを拡張する方法について、疑問に思っているかもしれませんが、少し後でそれは説明します。

異なる入力変数

それではもう少し複雑な問題に取り組みましょう。以下のコードsummarise()は、入力変数を変えて、3つの要約を計算する重複ステートメントを示しています。

summarise(df, mean = mean(a), sum = sum(a), n = n())
#> # A tibble: 1 x 3
#>    mean   sum     n
#>   <dbl> <int> <int>
#> 1     3    15     5
summarise(df, mean = mean(a * b), sum = sum(a * b), n = n())
#> # A tibble: 1 x 3
#>    mean   sum     n
#>   <dbl> <int> <int>
#> 1   8.6    43     5

これを関数に変えるには、まず基本的なアプローチをインタラクティブにテストすることから始めます。変数をquo()でクオートしてから、dplyrの!!を呼び出し、クオートを解除します。複雑な式の中のどこでも、クオートを外すことができることに注意してください。

my_var <- quo(a)
summarise(df, mean = mean(!! my_var), sum = sum(!! my_var), n = n())
#> # A tibble: 1 x 3
#>    mean   sum     n
#>   <dbl> <int> <int>
#> 1     3    15     5

またquo()は、dplyrの呼び出しをラップして、dplyrの観点から何が起こるのかを確認することもできます。これはデバッグに非常に便利なツールです。

quo(summarise(df,
  mean = mean(!! my_var),
  sum = sum(!! my_var),
  n = n()
))
#> <quosure>
#> expr: ^summarise(df, mean = mean(^a), sum = sum(^a), n = n())
#> env:  global

これで、コードを関数に置き換えて(quo()enquo()に置き換えることを忘れないでください)、関数が機能することを確認できます。

my_summarise2 <- function(df, expr) {
  expr <- enquo(expr)

  summarise(df,
    mean = mean(!! expr),
    sum = sum(!! expr),
    n = n()
  )
}
my_summarise2(df, a)
#> # A tibble: 1 x 3
#>    mean   sum     n
#>   <dbl> <int> <int>
#> 1     3    15     5
my_summarise2(df, a * b)
#> # A tibble: 1 x 3
#>    mean   sum     n
#>   <dbl> <int> <int>
#> 1   8.6    43     5

異なる入出力変数

次の課題は、出力変数の名前を変えることです。

mutate(df, mean_a = mean(a), sum_a = sum(a))
#> # A tibble: 5 x 6
#>      g1    g2     a     b mean_a sum_a
#>   <dbl> <dbl> <int> <int>  <dbl> <int>
#> 1     1     1     5     4      3    15
#> 2     1     2     3     2      3    15
#> 3     2     1     4     1      3    15
#> 4     2     2     1     3      3    15
#> # … with 1 more row
mutate(df, mean_b = mean(b), sum_b = sum(b))
#> # A tibble: 5 x 6
#>      g1    g2     a     b mean_b sum_b
#>   <dbl> <dbl> <int> <int>  <dbl> <int>
#> 1     1     1     5     4      3    15
#> 2     1     2     3     2      3    15
#> 3     2     1     4     1      3    15
#> 4     2     2     1     3      3    15
#> # … with 1 more row

このコードは前の例と似ていますが、2つの新しいコツがあります。

  • 文字列を貼り付けることで、新しい名前を作成し、入力式を文字列に変換するquo_name()が必要です。
  • !! mean_name = mean(!! expr)は有効なコードではないので、rlangが提供するヘルパー演算子:=を使用する必要があります。
my_mutate <- function(df, expr) {
  expr <- enquo(expr)
  mean_name <- paste0("mean_", quo_name(expr))
  sum_name <- paste0("sum_", quo_name(expr))

  mutate(df,
    !! mean_name := mean(!! expr),
    !! sum_name := sum(!! expr)
  )
}

my_mutate(df, a)
#> # A tibble: 5 x 6
#>      g1    g2     a     b mean_a sum_a
#>   <dbl> <dbl> <int> <int>  <dbl> <int>
#> 1     1     1     5     4      3    15
#> 2     1     2     3     2      3    15
#> 3     2     1     4     1      3    15
#> 4     2     2     1     3      3    15
#> # … with 1 more row

複数の変数を取り込む

任意の数のグループ化変数を受け入れるように、my_summarise()を拡張します。3つの変更を加える必要があります。

  • 関数定義で...を使用して、関数が任意の数の引数を受け付けるようにします。
  • enquos()を使い、...ですべての式をリストとして捕捉します。
  • !!の代わりに!!!を使用して、group_by()へ引数をつなぎ合わせます。
my_summarise <- function(df, ...) {
  group_var <- enquos(...)

  df %>%
    group_by(!!! group_var) %>%
    summarise(a = mean(a))
}

my_summarise(df, g1, g2)
#> # A tibble: 4 x 3
#> # Groups:   g1 [2]
#>      g1    g2     a
#>   <dbl> <dbl> <dbl>
#> 1     1     1     5
#> 2     1     2     3
#> 3     2     1     3
#> 4     2     2     1

!!!は要素のリストを受け取り、それらを現在の呼び出しにつなぎ合わせます。!!!の下部を見て、...を考えてください。

args <- list(na.rm = TRUE, trim = 0.25)
quo(mean(x, !!! args))
#> <quosure>
#> expr: ^mean(x, na.rm = TRUE, trim = 0.25)
#> env:  global

args <- list(quo(x), na.rm = TRUE, trim = 0.25)
quo(mean(!!! args))
#> <quosure>
#> expr: ^mean(^x, na.rm = TRUE, trim = 0.25)
#> env:  global

いくつかの実用的な例を通してtidyevalの基本を学んだ今、理論に飛び込みましょう。これはここで学んだことを新しい状況に一般化するのに役立ちます。

クオート

クオートは、式を評価するのではなく、式を捕捉することです。式ベースの関数は、それらの引数をクオートし、そのコードを評価した結果ではなく、式としてRコードを取得します。あなたがRユーザーであれば、あなたはたぶん定期的に式をクオートします。Rで最も重要な引用演算子の1つは式です。統計モデルの指定によく使われています。

disp ~ cyl + drat
#> disp ~ cyl + drat

Rのbaseのもう1つの引用演算子quote()です。式ではなく生の表現式を返します。

# Computing the value of the expression:
toupper(letters[1:5])
#> [1] "A" "B" "C" "D" "E"

# Capturing the expression:
quote(toupper(letters[1:5]))
#> toupper(letters[1:5])

(2重引用符と呼ばれているにもかかわらず、"は、式ではなく文字列を生成します。したがって、この文脈では使われているクオート演算子ではありません。)

実際には、この式はコードとその実行環境を取り込むため、2つの良いオプションをもちます。単純な式でも、環境が異なると値が異なる可能性があるため、これは重要です。たとえばxは、次の2つの式では異なる値を参照しています。

f <- function(x) {
  quo(x)
}

x1 <- f(10)
x2 <- f(100)

表示すれば、式は同じに見えるかもしれません。

x1
#> <quosure>
#> expr: ^x
#> env:  0x6bf61e8
x2
#> <quosure>
#> expr: ^x
#> env:  0x6c4db30

しかし、rlang::get_env()で使用している環境を調べると、それらは異なっています。

library(rlang)

get_env(x1)
#> <environment: 0x6bf61e8>
get_env(x2)
#> <environment: 0x6c4db30>

さらに、rlang::eval_tidy()を使用して、これらの式を評価すると、異なる値になることがわかります。

eval_tidy(x1)
#> [1] 10
eval_tidy(x2)
#> [1] 100

これはRの重要な特性です。1つの名前は異なる環境では異なる値を参照することができます。これはdplyrにとっても重要です。なぜなら、あなたが呼び出した中で、変数とオブジェクトを組み合わせることを可能にするからです。

user_var <- 1000
mtcars %>% summarise(cyl = mean(cyl) * user_var)
#>      cyl
#> 1 6187.5

オブジェクトが環境を追跡するとき、エンクロージャを持っていると言われます。これが、Rの関数がクロージャと呼ばれることがある理由です。

typeof(mean)
#> [1] "closure"

このような理由から、quosuresという片側式(one-sided formulas)を参照するために特別な名前を使用します。片側式(one-sided formulas)は、環境とクオート(式を含む)です。

Quosuresは通常のRオブジェクトです。それらは変数に格納して調べることができます。

var <- ~toupper(letters[1:5])
var
#> ~toupper(letters[1:5])

# You can extract its expression:
get_expr(var)
#> toupper(letters[1:5])

# Or inspect its enclosure:
get_env(var)
#> <environment: R_GlobalEnv>

Quasiquotation

簡単に言えば、準引用符(Quasiquotation)は、与えられた事例において言語表現を表す記号を導入することを可能にし、異なる事例においてその言語表現として使用される。- Willard van Orman Quine

自動的なクオートはインタラクティブな使用のために、dplyrを非常に便利にします。しかし、もしあなたがdplyrでプログラムしたいのであれば、間接的に変数を参照するための何らかの方法が必要です。この問題の解決策はquasiquotationです。これにより、引用符でクオートされている式の中で直接評価できます。

quasiquotationは、1940年代にWillard van Orman Quineによって考案され、1970年代にLISPコミュニティによってプログラミングに採用されました。tidyevalフレームワーク内のすべての式ベースの関数はquasiquotationをサポートします。アンクオートは、式の一部の引用符が取り消されます。アンクオートには3つのタイプがあります。

  • 基本的(basic)
  • 引用符のないスプライシング(unquote splicing)
  • 引用符なしの名前(unquoting names)

アンクオート

最初の重要な操作は、基本的なアンクオートであり、関数式がUQ()で、!!と等価です。

# Here we capture `letters[1:5]` as an expression:
quo(toupper(letters[1:5]))
#> <quosure>
#> expr: ^toupper(letters[1:5])
#> env:  global

# Here we capture the value of `letters[1:5]`
quo(toupper(!! letters[1:5]))
#> <quosure>
#> expr: ^toupper(<chr: "a", "b", "c", "d", "e">)
#> env:  global
quo(toupper(UQ(letters[1:5])))
#> <quosure>
#> expr: ^toupper(<chr: "a", "b", "c", "d", "e">)
#> env:  global

他の引用符で囲まれた式の引用符を外すこともできます。そのような象徴的なオブジェクトであるアンクオートは、表現式を処理するための強力な方法を提供します。

var1 <- quo(letters[1:5])
quo(toupper(!! var1))
#> <quosure>
#> expr: ^toupper(^letters[1:5])
#> env:  global

quosuresは環境を追跡しており、tidyeval関数はそれらを評価する方法を知っているので、quosuresのクオートを安全に解除できます。これにより、クオートとアンクオートをより親密に使い分けれます。

my_mutate <- function(x) {
  mtcars %>%
    select(cyl) %>%
    slice(1:4) %>%
    mutate(cyl2 = cyl + (!! x))
}

f <- function(x) quo(x)
expr1 <- f(100)
expr2 <- f(10)

my_mutate(expr1)
#>   cyl cyl2
#> 1   6  106
#> 2   6  106
#> 3   4  104
#> 4   6  106
my_mutate(expr2)
#>   cyl cyl2
#> 1   6   16
#> 2   6   16
#> 3   4   14
#> 4   6   16

関数式は、!の優先順位が問題を引き起こす場合に役立ちます。

my_fun <- quo(fun)
quo(!! my_fun(x, y, z))
#> Error in my_fun(x, y, z): could not find function "my_fun"
quo(UQ(my_fun)(x, y, z))
#> <quosure>
#> expr: ^`~`(fun)(x, y, z)
#> env:  global

my_var <- quo(x)
quo(filter(df, !! my_var == 1))
#> <quosure>
#> expr: ^filter(df, (^x) == 1)
#> env:  global
quo(filter(df, UQ(my_var) == 1))
#> <quosure>
#> expr: ^filter(df, (^x) == 1)
#> env:  global

引用符なしスプライシング

2番目のアンクオートの操作は、アンクオートスプライシング(unquote-splicing)です。その関数式UQS()であり、構文上の近道は!!!です。これはベクトルを受け取り、ベクトルの各要素を周囲の関数呼び出しに挿入します。

quo(list(!!! letters[1:5]))
#> <quosure>
#> expr: ^list("a", "b", "c", "d", "e")
#> env:  global

アンクオートスプライシングの非常に便利な機能は、ベクトル名が引数名になることです。

x <- list(foo = 1L, bar = quo(baz))
quo(list(!!! x))
#> <quosure>
#> expr: ^list(foo = 1L, bar = ^baz)
#> env:  global

これにより、名前付きドットを取るdplyr verbを使ったプログラミングが簡単になります。

args <- list(mean = quo(mean(cyl)), count = quo(n()))
mtcars %>%
  group_by(am) %>%
  summarise(!!! args)
#> # A tibble: 2 x 3
#>      am  mean count
#>   <dbl> <dbl> <int>
#> 1     0  6.95    19
#> 2     1  5.08    13

変数名を設定する

最後のアンクオート操作は、引数名を設定することです。変数名を設定する方法の1つを説明しましたが、=の代わりに:=演算子を使用することもできます。 :=は、LHSとRHSの両方のアンクオートをサポートします。

LHSの規則は少し異なります。引用符で囲まれていないオペランドまたは文字列は、シンボルとして評価されます。

mean_nm <- "mean"
count_nm <- "count"

mtcars %>%
  group_by(am) %>%
  summarise(
    !! mean_nm := mean(cyl),
    !! count_nm := n()
  )
#> # A tibble: 2 x 3
#>      am  mean count
#>   <dbl> <dbl> <int>
#> 1     0  6.95    19
#> 2     1  5.08    13

ブートストラップ信頼区間(パーセンタイル法)

はじめに

ブートストラップ信頼区間(パーセンタイル法)を書くことがあったので、その備忘録。今回の例では、ワイブル分布 (Weibull distribution) を例にする。この分布は、サバイバル分析などで出てくるように、人の死亡確率をモデリングする際に役立つ。人の死亡確率というと角が立ちそうだが、ECサイトの会員の離脱率なんかをイメージすると良いかもしれない。つまり、観測している期間において死亡確率は一定というよりも、死亡確率は時間と共に高くなり、変化するという仮定を起きたい時に便利。この反対もしかり。その事象が発生するまでの時間を確率変数と考えると、その確率変数が従う分布はワイブル分布となる。

ブートストラップ

ブートストラップサンプリングとは、得られたサンプルデータから、リサンプリングを行うこと。そこから、リサンプルデータから統計量を求め、推測を行う方法のことをブートストラップ〇〇とか言う。ブートストラップ信頼区間(パーセンタイル法)ではリサンプリングされた値を使ってパーセンタイルを計算すると、その値が取りうる信頼区間が得られる。これを各値で計算していくことで得られる。

{fitdistrplus}パッケージを使っている。{fitdistrplus}パッケージの詳細はこちら

library(fitdistrplus)
library(tidyverse)

set.seed(123)
n <- 100
df_prob <- rweibull(n = n, shape = 2, scale = 10)
x <- seq(0.1, 30, len = n)


loops <- 100
booted_df <- map_dfc(1:loops, function(i) {
  xi <- sample(df_prob, size = length(df_prob), replace = TRUE)
  MLE_est <- suppressWarnings(fitdist(xi, distr = "weibull",method = "mle"))  
  dweibull(x, shape = MLE_est$estimate["shape"],  scale = MLE_est$estimate["scale"])
})

dat1 <- booted_df %>% 
  bind_cols(x %>% as_tibble() %>% rename(x = value)) %>% 
  gather(key = boot_ind, val = y, -x)

dat2 <- booted_df %>% 
  array_tree(., 1) %>% 
  map_dfc(., quantile, c(0.025, 0.5, 0.975)) %>% 
  t() %>% as_tibble() %>%
  bind_cols(x %>% as_tibble()) %>%
  setNames(c("2.5%", "50%", "97.5%", "x")) %>%
  gather(key = quantile, val = y, -x)

ggplot() +
  geom_line(data = dat1, aes(x, y, group = boot_ind), col = "#006E4F", alpha = 0.2, size = 0.25) +
  geom_line(data = dat2, aes(x, y, group = quantile), col = "#006E4F", size = 0.5, linetype = "longdash") +
  theme_classic() + xlab("x") + ylab("Probability density") + 
  ggtitle("Weibull Probability Distribution")

f:id:AZUMINO:20190601044905p:plain