まだ厨二病

RとPyhtonを使ったデータ分析・統計解析らへんの話題をしていくだけ

⭐️Rを使ったモデル構築の最善策を求めて: {dplyr} + {tidyr} + {broom} + {purrr}を使ったアプローチ

RStudioのチーフサイエンティスト、Hadley Wickham(ハドリー)が2月に行った講演のビデオがYouTubeに上がっていたので観た。

"Making Data Analysis Easier"というタイトルでの発表(スライドでは"Managing many models"になっているけど)で、ハドリー自身が考えている、データサイエンスに必要な可視化やモデリングを効率的に行うための手法について、彼の開発してきたパッケージを中心に説明している。

www.youtube.com

分かりやすく、具体例を交えた内容なので、是非YouTubeの動画を観てもらうのが良いと思うが、自分の頭を整理するためにもここでまとめておく。なお、発表スライドはクリエイティブ・コモンズライセンス3.0のもと、表示・非営利のラインセンスで再利用可能となっている。

Hadley Wickham (Chief Scientist, RStudio) 2016. Managing many models. http://wombat2016.org/abstracts/hadley.html

✨ データサイエンスにおいて重要な7つの事柄

ハドリーは自身の発表でしばしば次のような図を提示する。

f:id:u_ribo:20160326224249p:plain

{DiagrammeR}パッケージを使って作図した。ソースコードgistに

これはハドリーがデータサイエンスを行っていく上で重要だと考えているプロセスとその流れを表したもので、ハドリーはこのプロセスに沿ったパッケージを開発してきた*1(そういう風に捉えるとますますハドリーが尊い存在に思えてくる)。上の図で四角い枠に囲われたものがここの過程を示しており(探索的データ解析の中には整形、可視化、モデリングが含まれる)、{package}のように表されているのがそのプロセスを進めていくのに便利なRのパッケージ名となっている。

すなわちデータを取得した後、解析のために利用しやすいデータ形式へと整頓し、データの変形や可視化、モデル化を繰り返して、その情報を伝達したり自動的な処理のプログラムを組むことをハドリーは重視している。そして今回の発表では、この中の「探索的データ解析」に焦点を当てている。これはハドリーのデータサイエンスプロセスでも中心となるものであり、その後の議論を左右する重要な部分である。

ちなみに、以前にもハドリーのこのプロセスに沿ったデータ整形に関する記事を書いたのでそちらも参考に。

uribo.hatenablog.com

🔰 データの階層構造を理解する

現実的な問題: gapminderデータを例に

我々が扱うデータはときに2つ以上の、ときには大量の変数をもっている。例えば{gapminder}パッケージにはgapminderというデータセットがあり、6つの変数をもっている。うち2つの変数は対象の国countryとその地域(大陸)名continentで、これらはカテゴリー変数である。また、年year列はデータが得られた西暦を表している。

library(dplyr)
data("gapminder", package = "gapminder")
# データセットの変数の型を把握する
gapminder %>% glimpse()
## Observations: 1,704
## Variables: 6
## $ country   (fctr) Afghanistan, Afghanistan, Afghanistan, Afghanistan,...
## $ continent (fctr) Asia, Asia, Asia, Asia, Asia, Asia, Asia, Asia, Asi...
## $ year      (int) 1952, 1957, 1962, 1967, 1972, 1977, 1982, 1987, 1992...
## $ lifeExp   (dbl) 28.801, 30.332, 31.997, 34.020, 36.088, 38.438, 39.8...
## $ pop       (int) 8425333, 9240934, 10267083, 11537966, 13079460, 1488...
## $ gdpPercap (dbl) 779.4453, 820.8530, 853.1007, 836.1971, 739.9811, 78...

ここで雑に出生時平均余命lifeExpの変化を年代yearごとにプロットしてみると次の図が得られる。

library(ggplot2)
gapminder %>% ggplot(aes(year, lifeExp, group = country)) + 
    geom_line()

f:id:u_ribo:20160326235655p:plain

20世紀の中盤から21世紀向かって、出生時平均余命が上がっている感じがある。感じがある、を雑に確かめるためには線形回帰を使えば良い。一方で、この図の中には全体的な傾向(右肩上がり)と異なり、時々ガクッと下がったり、それほど上昇していない線もあるようだ。つまり、調査が行われた国ごとに傾向が違うのではないだろうか?それを確認する方法の一つとして、国ごとに線形回帰を行ってみる。

...がちょっと待て。gapminderデータセットにはいくつの国が含まれている?

gapminder$country %>% nlevels()
## [1] 142

これらを一つ一つしらみつぶしに調べていくのは大変面倒だ。

ある基準でデータを内包する

やりたいことを整理してみよう。gapminderデータセットがもつ変数、出生時平均余命lifeExpは年yearに従って増加している、という仮説を検証するために国countryごとにその傾向を線形回帰モデルを使って調べてみる、ということだ。

ここで線形回帰を行うために必要な変数、lifeExpyearは変動するが、国と大陸名は固定されていることに注目する。つまりこのカテゴリー変数を用いてデータを入れ子構造に内包することが可能であるということだ。階層構造をもつデータフレームを生成するには{tidyr}パッケージのnest()関数を使う。また{dplyr}パッケージのgroup_by()で基準となる変数を指定しておくことが大事である。

library(tidyr)
nest_by_country <- gapminder %>% 
  group_by(continent, country) %>% 
  nest()

# 日本のデータは次の位置に格納されている
gapminder$country %>% factor() %>% levels() %>% grep("Japan", .)
## [1] 67
nest_by_country$data[[67]]
## Source: local data frame [12 x 4]
## 
##     year lifeExp       pop gdpPercap
##    (int)   (dbl)     (int)     (dbl)
## 1   1952  63.030  86459025  3216.956
## 2   1957  65.500  91563009  4317.694
## 3   1962  68.730  95831757  6576.649
## 4   1967  71.430 100825279  9847.789
## 5   1972  73.420 107188273 14778.786
## 6   1977  75.380 113872473 16610.377
## 7   1982  77.110 118454974 19384.106
## 8   1987  78.670 122091325 22375.942
## 9   1992  79.360 124329269 26824.895
## 10  1997  80.690 125956499 28816.585
## 11  2002  82.000 127065841 28604.592
## 12  2007  82.603 127467972 31656.068
# 元のgapminderデータセットから日本を抽出し、不要な変数を取り除いたデータと一致する
all.equal(gapminder %>% dplyr::filter(country == "Japan") %>% dplyr::select(-continent, -country), 
          nest_by_country$data[[67]])
## [1] TRUE

この日本のデータにのみ絞って、線形回帰を行ってみよう。

nest_by_country$data[[67]] %>% dim()
## [1] 12  4
lm.res <- nest_by_country$data[[67]] %>% lm(lifeExp ~ year, data = .)
lm.res %>% summary() %>% {
  .$r.squared %>% print() # 決定係数(モデルの当てはまりの良さ)
  .$coefficients %>% print()  # 係数(グラフの傾きと切片)
  .$residuals # 残差(モデルによる予測値と観測値のズレを示す)
}
## [1] 0.9595956
##                 Estimate  Std. Error   t value         Pr(>|t|)
## (Intercept) -623.7469389 45.33137200 -13.75972 0.00000007989525
## year           0.3529042  0.02289954  15.41097 0.00000002695784

##           1           2           3           4           5           6 
## -2.09205128 -1.38657226  0.07890676  1.01438578  1.23986480  1.43534382 
##           7           8           9          10          11          12 
##  1.40082284  1.19630186  0.12178089 -0.31274009 -0.76726107 -1.92878205

Rではこのようにすれば良い。これだけで十分な情報が得られて、Rは素晴らしいな、と思う反面、やはりすべての国に対してこれをやるのは辛い。

{broom}パッケージでRによるモデルの結果をよしなに

先の例では、特に線形回帰を実施する関数lm()の結果を要約summary()したものについて、その要素にアクセスすることでモデルの結果を確認した。Rのモデリング用の関数の多くはこのように、解析の結果を要素として(全体ではリストとして)もつことがほとんどなので、要素名を指定すると必要な情報を得ることができるがいささか面倒だ。そこで{broom}パッケージを利用する。

{broom}パッケージはRの統計解析用の関数の出力を整形し、利用者が再利用しやすい形(データフレームクラスオブジェクト)で出力する。

library(broom)
nest_by_country$data[[67]] %>% lm(lifeExp ~ year, data = .) %>% {
  glance(.) %>% print()
  tidy(.) %>% print()
  augment(.) %>% .$.resid
}
##   r.squared adj.r.squared    sigma statistic          p.value df  logLik
## 1 0.9595956     0.9555552 1.369194   237.498 0.00000002695784  2 -19.704
##        AIC      BIC deviance df.residual
## 1 45.40799 46.86271 18.74691          10
##          term     estimate   std.error statistic          p.value
## 1 (Intercept) -623.7469389 45.33137200 -13.75972 0.00000007989525
## 2        year    0.3529042  0.02289954  15.41097 0.00000002695784

##  [1] -2.09205128 -1.38657226  0.07890676  1.01438578  1.23986480
##  [6]  1.43534382  1.40082284  1.19630186  0.12178089 -0.31274009
## [11] -0.76726107 -1.92878205

出力する値自体は変わりがないが、{broom}を使った統計解析の結果の出力はデータフレームクラスのオブジェクトとなっている点が大きく異なっている。このことが、次の効率的なモデリングを行うのに大切になる。

{dplyr}と{purrr}および{broom}による効率的なモデリング

国別の傾向を知りたいが、問題はデータを分割する基準となる国の数が多いことだ。こうした問題はデータ分析を行う際には頻繁に生じる。性別や都道府県ごとに処理を行う、というのは誰もがやったことがあるのではないだろうか。

ハドリーは、こうした作業をよりプログラミング的に行うべく、do()という関数を用意した。do()関数を利用すると同じく{dplyr}パッケージの関数group_dy()で基準を指定した後、目的の処理を記述するだけで、基準を適用した処理の結果が特殊な階層に保存される。

[http://qiita.com/uri/items/27f1497778b4385acb25:embed:cite]

(do_by_country <- gapminder %>% 
  group_by(country) %>% 
  do(data = lm(lifeExp ~ year, data = .) %>% tidy()))
## Source: local data frame [142 x 2]
## Groups: <by row>
## 
##        country               data
##         (fctr)             (list)
## 1  Afghanistan <data.frame [2,5]>
## 2      Albania <data.frame [2,5]>
## 3      Algeria <data.frame [2,5]>
## 4       Angola <data.frame [2,5]>
## 5    Argentina <data.frame [2,5]>
## 6    Australia <data.frame [2,5]>
## 7      Austria <data.frame [2,5]>
## 8      Bahrain <data.frame [2,5]>
## 9   Bangladesh <data.frame [2,5]>
## 10     Belgium <data.frame [2,5]>
## ..         ...                ...

do()関数内で記述した処理した内容は次のようにして取り出せる。group_by()do()を処理の間に挟むだけで対象の国ごとの線形回帰の結果が得ることができた。素晴らしい!

do_by_country$data[[1]]
##          term     estimate   std.error statistic          p.value
## 1 (Intercept) -507.5342716 40.48416195 -12.53661 0.00000019340553
## 2        year    0.2753287  0.02045093  13.46289 0.00000009835213
do_by_country$data %>% length()
## [1] 142

すべての国について見ていくのはまた面倒なのでtidyr::unnest()を使う。unnest()関数は階層構造のであるデータ形式を展開するものであり、nest()と対をなす関数である。do()の中で実行した処理を再びデータフレームとしてまとめるのに便利である。

do_by_country %>% unnest()
## Source: local data frame [284 x 6]
## 
##        country        term      estimate    std.error  statistic
##         (fctr)       (chr)         (dbl)        (dbl)      (dbl)
## 1  Afghanistan (Intercept)  -507.5342716 40.484161954 -12.536613
## 2  Afghanistan        year     0.2753287  0.020450934  13.462890
## 3      Albania (Intercept)  -594.0725110 65.655359062  -9.048348
## 4      Albania        year     0.3346832  0.033166387  10.091036
## 5      Algeria (Intercept) -1067.8590396 43.802200843 -24.379118
## 6      Algeria        year     0.5692797  0.022127070  25.727749
## 7       Angola (Intercept)  -376.5047531 46.583370599  -8.082385
## 8       Angola        year     0.2093399  0.023532003   8.895964
## 9    Argentina (Intercept)  -389.6063445  9.677729641 -40.258031
## 10   Argentina        year     0.2317084  0.004888791  47.395847
## ..         ...         ...           ...          ...        ...
## Variables not shown: p.value (dbl).

これまでの一連の処理を{purrr}パッケージを使ってよりプログラミング的に記述することができる。

library(purrr)
map_by_country <- gapminder %>% 
  split(.$country) %>% 
  map(., ~ lm(lifeExp ~ year, data = .) %>% tidy())

all.equal(do_by_country$data[[1]], map_by_country[[1]])
## [1] TRUE

どちらを使っても良いが、ハドリーは{purrr}を推している感じがある(講演ではdo()は出てこずmap()の例が出ている)。{purrr}を使った応用例がYouTubeの動画で見れるのできになる人はチェックを。

⭐ まとめ

この記事の内容であまり触れられていない部分もあるが、ハドリーは以下の項目をまとめとしてあげている。

  1. 関連するオブジェクトをリスト - 列の階層にまとめることで、モデリングの結果を取得、整理しやすくなる
  2. {purrr}を使った)関数型プログラミングを体系的に学ぶことで、オブジェクトではなく関数の働きに焦点を当てた処理ができるようになる
  3. {broom}パッケージを使うと、統計解析モデルの結果が整形されたデータフレームとして得られ、可視化や集計に便利である

f:id:u_ribo:20160327011908p:plain

この講演の中でも、彼の著書の中でも、ハドリーはR利用者に自作関数を作ることを勧めている気がする(確かに効率を上げるには関数にしたほうが良いが、関数への理解が必要だ)。それはハドリーがRの標準グラフィックスに対抗するようなggplot2パッケージを開発したり、文字列や日付操作の改善に努め、Rにおける処理のコアとなるデータの読み込みや処理・操作に重点を置いたパッケージを作っていく過程で、Rの標準関数のイマイチな挙動にうんざりしたから、という気がしなくもない。そしていよいよ関数型言語としての利用を推進している、ように感じる。

ハドリーはまた別の発表をするらしい... BARUG at Strata: "Official"­ March Meeting - Bay Area useR Group (R Programming Language) (San Francisco, CA) - Meetup

こちらも資料が公開されたら必見(内容が近い感じがするけど)。

💻 実行環境

package version source
broom 0.4.0 CRAN (R 3.2.2)
dplyr 0.4.3.9000 Github (<hadley/dplyr@9bae2aa>)
magrittr 1.5 Github (<smbache/magrittr@00a1fe3>)
purrr 0.2.1 CRAN (R 3.2.3)
remoji 0.1.0 Github (<richfitz/remoji@dc00779>)
tidyr 0.4.1 CRAN (R 3.2.3)

*1:Data manipulation with dplyr, Getting your data into R https://vimeo.com/125174198