まだ厨二病

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

⭐️手を動かしながら学ぶモダンなデータ操作とtidyなデータ(2015年版)

R Advent Calendar 2015の第一日目です。

今日はタイトルの通り、{dplyr}{tidyr}パッケージを使ったデータの整形と集計処理について、実際のデータを交えながら紹介したいと思います(タイトルは流行りの本からとってきました。オマージュです)。

{dplyr}{tidyr}パッケージ、関数の使い方を紹介した記事はあっても、実際のデータを扱った記事を検索しても、日本語の記事がほとんど見つからなかったので、tidyなデータ形式について普及させるために記事を書こうというところです。

一応、自分が集められたtidyデータについての記事へのリンクを貼っておきます。

ちょっと前に海外の方がtidyデータについての記事を書かれていましたので、こちらも参考になります。

🔰 tidyデータとは

Rだけに限った話ではありませんが、データ解析のために用いられる時間の大半は、手持ちのデータを解析処理を行う関数やツールが扱えるようにするためのデータの整形作業です。いわゆる前処理と呼ばれるこの作業は、データ解析の最初の段階としての作業であり、のちの解析にも大きな影響を及ぼす重要な過程と成っています。

R界の「神」ことHadley Wickhamの提唱するtidyデータはデータを規則正しい形に整形することで、データの操作、解析、可視化を柔軟に扱えるようにすることを目指しています。より詳しい内容について知りたい方はHadleyの論文を読んでください。また、日本語での記事として、市川(2014)の解説がとても参考になります。

tidyデータの要点は次のようになります。

  • 一つの変数に対して一つの列が与えられるべし(同一の変数を複数の列にまたがせない)
  • それぞれの観測値は一つの行に収めるべし(同一の個体に対して複数の行を与えない)
  • 観測データの集合は表形式で表現することができる

📦 {tidyr}パッケージ

Hadleyが中心になって開発されるRパッケージの一つに{tidyr}というものがあります。これはRでのデータ処理に優れた{reshape2}をベースにしており、かつ上記のtidyデータの概念を反映しており、手持ちのデータをtidy形式に(柔軟な操作ができるように)する関数や、その他、データの処理のための関数を備えています。

主要な関数

  • gather(): 複数の変数に分かれている共通の変数を1つの変数内の観測値としてまとめます(spread()と対応)
  • spread(): 項目とその値を元にして、変数を分割します(gather()と対応)
  • separate(): 1つの列を複数列に分割します

補助的な関数

  • extract_numeric(): 値から数値を取り出す(¥1000や173cm、-2%のような単位を含むデータから数値だけを抽出します)
  • unite(): 複数列を1列に結合

これらの関数と、すでにRでのデータ操作に関して確固たる地位をしめているであろう{dplyr}パッケージを使っていきます。

🔬 扱うデータについて: モニタリングサイト1000

せっかくなので公開されているデータを使用します。

環境省が行っている業務の一つに「モニタリングサイト1000(通称: モニセン)」というものがあります。これはさまざまな生態系を対象として日本全国各地にモニタリングサイトを設置し、長期的な観測によりその変化や構造について明らかにしようという壮大なプロジェクトです。みなさんどこかで耳にしたことがあるのではないでしょうか。

このモニセン業務では、収集されたデータを公開しており、利用用途などのアンケートに答えることでそのデータをダウンロードできます。

今回は「モニタリングサイト1000森林・草原調査」で行われている調査の1つのである毎木調査(特定の面積内に生育する樹木を対象としたサイズの調査)から、北海道の「雨龍」プロットのデータを例にしてtidyデータの説明を行いたいと思います。

名前がいいですよね。親しみを感じます。行ってみたいプロットの一つです。

このデータでは、2005年から2012年の8年間分の毎木のデータが収められています。また、幹の状態を記録するためのメモが追加されています。

💾 データの読み込みからtidy形式への変形

今回利用するRパッケージを読み込みます。

library(ggplot2)
library(readxl)
library(broom)
library(tidyr)
library(dplyr)
quartzFonts(YuGo = quartzFont(rep("YuGo-Medium", 4)))
theme_set(theme_classic(base_size = 12, base_family = "YuGo"))
df_uryu <- read_excel(path = "UR-BC1-TreeGbh-2005-2012-ver1.xls",
                    sheet  = 1,
                    skip   = 86,
                    na     = "na")

このデータシートでは、欠損値が「na」という文字列で与えられているため、read_excel()na引数に指定しておきます。こうすることでRでも欠損値として扱ってくれるようになります。

names(df_uryu)
##  [1] "mesh_xcord"          "mesh_ycord"          "tag_no"             
##  [4] "indv_no"             "stem_xcord"          "stem_ycord"         
##  [7] "spc_japan"           "gbh05"               "gbh06"              
## [10] "gbh07"               "gbh08"               "gbh09"              
## [13] "gbh10"               "gbh11"               "gbh12"              
## [16] "note05"              "note06"              "note07"             
## [19] "note08"              "note09"              "note10"             
## [22] "note11"              "note12"              "s_date05"           
## [25] "s_date06"            "s_date07"            "s_date08"           
## [28] "s_date09"            "s_date10"            "s_date11"           
## [31] "s_date12"            "stem_id"             "indv_id"            
## [34] "stem_xcord_original" "stem_ycord_original" "tag_no05"           
## [37] "indv_no05"           "tag_no09"            "indv_no09"

読み込んだデータの変数名を確認しました。ちょっと多いので、分析に利用しない列をデータから除いてデータをちょっと軽くします。変数の選択にはdplyr::select()が便利です。正規表現や条件分岐を利用して列の選択ができます。この時に使用できるオプションについては\@hoxo_m さんの記事が詳しいので参考にしてください。

df_uryu %<>% dplyr::select(-ends_with("cord"), -contains("s_date"), 
                           -starts_with("stem_"), -starts_with("indv_no"), 
                           -tag_no05, -tag_no09)

変数が多い(19列)のでglimpse()を使って横方向にすべての変数のデータのいくつかを表示します。

glimpse(df_uryu)
## Observations: 849
## Variables: 19
## $ tag_no    (dbl) 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 1...
## $ spc_japan (chr) "ナナカマド", "トドマツ", "トドマツ", "シナノキ", "イタヤカエデ", "ナナカマド", ...
## $ gbh05     (dbl) 17.0, 161.9, 105.0, 28.6, 89.7, 23.2, 84.9, 223.7, 1...
## $ gbh06     (dbl) 17.5, 163.2, 109.3, 30.3, 90.7, 24.5, 84.9, 224.2, 1...
## $ gbh07     (dbl) 18.4, 163.8, 109.3, 31.8, 93.1, 25.5, 85.0, 224.8, 2...
## $ gbh08     (dbl) 19.2, 165.2, 114.3, 33.3, 92.9, 25.7, 85.1, 224.8, 2...
## $ gbh09     (dbl) 20.2, 165.2, 114.3, 35.4, 95.9, 27.5, 85.1, 224.8, 2...
## $ gbh10     (chr) "21.000000", "165.200000", "117.000000", "37.000000"...
## $ gbh11     (dbl) 22.0, 165.4, 117.3, 39.0, 96.2, 29.0, NA, 224.9, 200...
## $ gbh12     (chr) "22.000000", "167.400000", "121.500000", "40.000000"...
## $ note05    (dbl) NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ note06    (dbl) NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ note07    (dbl) NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ note08    (dbl) NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ note09    (dbl) NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ note10    (chr) NA, NA, NA, NA, NA, NA, "枯死", NA, NA, NA, NA, NA, NA...
## $ note11    (dbl) NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ note12    (chr) NA, NA, NA, NA, NA, NA, "NA", NA, NA, NA, NA, NA, NA...
## $ indv_id   (chr) "2005_374", "2005_163", "2005_164", "2005_373", "200...

ここで変数について説明します。

  • tag_no... 調査する幹を識別するためのタグ番号
  • spc_japan... 種名(和名)
  • gbhで始まる変数... 胸高周囲長。地上からおよそ1.3mの高さで計測される幹の周囲です。森林樹木の調査では樹木のサイズを評価するためにこの指標が標準的に利用されます。末尾の数字は調査した年を示しています。
  • noteで始まる変数... 幹の状態や調査状況について記録したメモ。これも末尾の数字は調査した年を示しています。
  • indv_id... 個体を識別するための番号

gbh10gbh12の列が本来は数値であって欲しいのですが、文字となっていますね。これを修正するために次の処理を行います。dplyr::mutate()はデータの変数に対して、計算を行ったり、関数による処理を適用する関数で、以下のdplyr::mutate_each()複数の変数に対して同様の処理を施す関数です。次の処理で、gbhで始まる変数をすべて変数の値を数値に直します。

df_uryu %<>% dplyr::mutate_each(funs(as.numeric), starts_with("gbh"))

ここからデータをtidy形式にしていきます。上述の通り、gbh_*, note_*はそれぞれ年ごとに記録された値ですので、これはtidyデータの概念から外れています(一つの変数に対して一つの列が与えられるべし)。そこで、この年ごとに値の異なる変数を「調査年」と「値」に分けてtidyにしてみましょう。

gbh_*, note_*を一旦分離して結合することを試みます。

df_uryu %>% {
  tmp_df_gbh <<- dplyr::select(., -contains("note")) %>% 
    gather(key = year, value = gbh, select = -c(tag_no, spc_japan, indv_id)) %>% 
  dplyr::mutate(year = extract_numeric(year) + 2000)
  tmp_df_note <<- dplyr::select(., -contains("gbh")) %>% 
    gather(key = year, value = note, select = -c(tag_no, spc_japan, indv_id)) %>% 
  dplyr::mutate(year = extract_numeric(year) + 2000)
}

今、tmp_df_gbhtmp_df_noteというデータフレームができました。以下のように年ごとに記録されていたgbhの値が1列に収まっており、調査された年についての情報はyear列にあります。noteについても同様です。

tidyデータ形式になったtmp_df_gbhtmp_df_noteを結合して、完全なtidyデータを作りましょう。ここではdplyr::inner_join()を使います。この関数は、2つのデータフレームを比較して、相互に値を補いながら結合してくれます。

df_uryu_tidy <- tmp_df_gbh %>% inner_join(tmp_df_note)
head(df_uryu_tidy) %>% kable(format = "markdown")
tag_no spc_japan indv_id year gbh note
1 ナナカマド 2005_374 2005 17.0 NA
2 トドマツ 2005_163 2005 161.9 NA
3 トドマツ 2005_164 2005 105.0 NA
4 シナノキ 2005_373 2005 28.6 NA
5 イタヤカエデ 2005_162 2005 89.7 NA
6 ナナカマド 2005_372 2005 23.2 NA

🔨 tidyデータを操作する

tidyデータは、データ操作、モデリング、可視化するために利用しやすい形式であると最初に述べましたが、それを確かめるために次にあげるいくつか処理を加えてみましょう。

  • データ操作
    • 列の分割... tidyr::separate()
    • 新たな変数の作成... dplyr::mutate()
    • 条件に応じたデータの抽出... dplyr::filter()
    • グループごとの集計... dplyr::group_by(), dplyr::summarise()x
  • 可視化
  • モデリング

さて、これまでに使ってない変数として、indv_idがあります。これは、幹が最初に記録された年と番号からなる変数ですので、2つの変数に分割します。変数の分割にはtidyr::separate()を利用します。col引数に対象とする変数名を、into引数に新たに生成する変数名を与えます。sep引数に与える値は分割の基準(文字)です。

df_uryu_tidy %<>% separate(col    = indv_id, 
                           into   = c("obs_first", "id"), 
                           sep    = "_",
                           remove = TRUE)

次は既存の変数を元にして、新たな変数を追加する、ということを実行します。今回のデータでは、GBH(幹の周囲長)が記録されているので、この幹の肥大成長量を計算して新たな変数にしましょう。

生物の成長量を評価する指標として、相対成長速度 Relative growth rate (RGR)というものがあります。これは以下の式で定義されます。

f:id:u_ribo:20151201063241p:plain

biomass1, biomass2はそれぞれ、ある時間におけるサイズを、t1, t2は時間を意味します。

これをtag_no(幹を識別するタグ番号)ごとに計算します。そのためにdplyr::group_by()tag_noを指定します。またその前に実行しているarrange()は、結果の見やすさのための並び替えです。

df_uryu_tidy %<>% arrange(tag_no, year) %>% 
  group_by(tag_no) %>% 
  mutate(rgr = (log(gbh) - log(gbh[1]) ) / (year - min(year, na.rm = TRUE)) ) %>% 
  ungroup() %>% 
  mutate(rgr = gsub("NaN|Inf", NA, rgr)) %>% 
  mutate(rgr = as.numeric(rgr) %>% round(digits = 4))

今度はデータに基づいた可視化を行います。まず、調査が行われた年・種ごとに何本の記録があるのかを棒グラフで確認します。また、集計した年・種ごとの幹本数を集計した値を元に、本数が50以上ある種をtgt_spというベクトルで保存しておきます(これは後で利用します)。

df_uryu_tidy %>% dplyr::select(-note) %>% 
  dplyr::filter(!is.na(gbh)) %>% 
  group_by(year) %>% 
  summarise(stem_num = n()) %>% 
  ggplot(aes(year, stem_num, fill = factor(year))) + 
  geom_bar(stat = "identity", alpha = 0.5) +
  theme(axis.text.x = element_text(angle = 40, hjust = 1)) +
  guides(fill = FALSE)

f:id:u_ribo:20151201052459p:plain

df_uryu_tidy %>% dplyr::filter(!is.na(gbh)) %>% 
  group_by(year, spc_japan) %>% 
  summarise(stem_num = n()) %>% {
    tgt_sp <<- dplyr::filter(., stem_num > 50) %$% spc_japan %>% unique()
    p <<- ggplot(., aes(year, stem_num, color = spc_japan)) +
      geom_point(aes(shape = ifelse(stem_num > 50, 16, 21))) +
      geom_line(aes(linetype = ifelse(stem_num > 50, 1, 3))) +
      scale_shape_identity() +
      scale_linetype_identity()
  }
p

f:id:u_ribo:20151201052528p:plain

50本以上の記録がある種は次の5種になりました。この5種について一般化線形モデルによる統計モデリングを実行しましょう。

tgt_sp
## [1] "イタヤカエデ" "トドマツ"     "ナナカマド"   "ハリギリ"    
## [5] "ミズナラ"

dplyr::filter()によってデータに制限を与えます。先ほど作成したtgt_spに含まれる種のうち、2012年のデータから、rgrが欠損しておらず、0よりも大きいデータについて対象にしています。

相対成長速度は一般に、サイズが大きくなると低下するので、これを一般化線形モデル GLMで確認します。全データをごちゃ混ぜにした状態で実行した結果が以下になります。

df_uryu_tidy %>% dplyr::filter(spc_japan %in% tgt_sp, 
                               year == 2012, 
                               !is.na(rgr), rgr > 0) %>% 
  glm(formula = rgr ~ gbh, 
      family = "Gamma"(link = "log"),
      data = .) %>% 
  tidy() %>% kable(format = "markdown")
term estimate std.error statistic p.value
(Intercept) -3.4690330 0.0459167 -75.55051 0
gbh -0.0081179 0.0005593 -14.51315 0

種ごとの傾向を確認したい時にはdplyr::do()を使うと良いです。dplyr::do()については、過去に書いていますのでそちらをご参考ください。

uryu_glm_broom <- df_uryu_tidy %>% dplyr::filter(spc_japan %in% tgt_sp,
                                                 year == 2012, 
                                                 !is.na(rgr), rgr > 0) %>% 
  group_by(spc_japan) %>% 
  do(glm(formula = rgr ~ gbh, 
         family = "Gamma"(link = "log"),
         data = .) %>% 
       tidy())
uryu_glm_broom %>% kable(format = "markdown")
spc_japan term estimate std.error statistic p.value
イタヤカエデ (Intercept) -3.2367434 0.1079893 -29.972816 0.0000000
イタヤカエデ gbh -0.0150706 0.0019277 -7.817903 0.0000000
トドマツ (Intercept) -3.6114314 0.1699936 -21.244509 0.0000000
トドマツ gbh -0.0067302 0.0015603 -4.313316 0.0000348
ナナカマド (Intercept) -3.6525532 0.1324492 -27.577018 0.0000000
ナナカマド gbh -0.0072686 0.0034717 -2.093685 0.0385247
ハリギリ (Intercept) -3.4240583 0.1027251 -33.332256 0.0000000
ハリギリ gbh -0.0081614 0.0016080 -5.075394 0.0000060
ミズナラ (Intercept) -3.1651068 0.0745401 -42.461813 0.0000000
ミズナラ gbh -0.0091480 0.0006453 -14.176492 0.0000000

この結果も図示しておきましょう。

df_uryu_tidy %>% dplyr::filter(spc_japan %in% tgt_sp,
                               year == 2012, 
                               !is.na(rgr), rgr > 0, rgr < 0.10) %>% 
  ggplot(aes(gbh, rgr)) +
  geom_point(aes(colour = spc_japan), size = 2, shape = 21, stroke = 1) +
  geom_smooth(method = "glm",
              method.args = list(family = "Gamma"(link = "log")),
              se = FALSE) +
  facet_wrap(~ spc_japan, scales = "free_x", ncol = 5) +
  guides(colour = FALSE)

f:id:u_ribo:20151201053641p:plain

さらにこうした種の効果を考慮するモデルとして、一般化線形モデルよりも一般化混合線形モデルが適しています、がここは割愛します。

♻ 元に戻す

tidyデータは便利ですが、入力する際にはやはり横方向に記録する形式が良いかと思います。なので、tidyデータを元の形式に戻すコードを書いておきます。

df_uryu_tidy %>% dplyr::select(-rgr) %>% {
  tmp_df_gbh <<- dplyr::select(., -note) %>% 
    unite(col = indv_id, obs_first, id, sep = "_") %>% 
    spread(key = year, value = gbh) %>% 
    dplyr::rename(gbh05 = `2005`,
                  gbh06 = `2006`,
                  gbh07 = `2007`,
                  gbh08 = `2008`,
                  gbh09 = `2009`,
                  gbh10 = `2010`,
                  gbh11 = `2011`,
                  gbh12 = `2012`)
  tmp_df_note <<- dplyr::select(., -gbh) %>% 
    unite(col = indv_id, obs_first, id, sep = "_") %>% 
    spread(key = year, value = note) %>% 
    dplyr::rename(note05 = `2005`,
                  note06 = `2006`,
                  note07 = `2007`,
                  note08 = `2008`,
                  note09 = `2009`,
                  note10 = `2010`,
                  note11 = `2011`,
                  note12 = `2012`)
}
df_uryu_untidy <- tmp_df_gbh %>% inner_join(tmp_df_note)
dim(df_uryu_untidy)
## [1] 849  19

{dplyr}と比べると関数が少ないので覚えやすいと思うので、皆さんもtidyってください。

💻 実行環境

devtools::session_info() %>% {
  print(.$platform)
  .$packages %>% dplyr::filter(`*` == "*") %>% kable(format = "markdown")
}
##  setting  value                       
##  version  R version 3.2.2 (2015-08-14)
##  system   x86_64, darwin13.4.0        
##  ui       X11                         
##  language En                          
##  collate  en_US.UTF-8                 
##  tz       Asia/Tokyo                  
##  date     2015-12-01
package * version date source
broom * 0.4.0 2015-11-30 CRAN (R 3.2.2)
dplyr * 0.4.3.9000 2015-10-28 Github (<hadley/dplyr@52d10f6>)
ggplot2 * 1.0.1.9003 2015-10-17 Github (<hadley/ggplot2@864d64f>)
knitr * 1.11.8 2015-10-19 Github (<yihui/knitr@a1b235d>)
magrittr * 1.5 2015-07-28 Github (<smbache/magrittr@effbadc>)
pipeR * 0.6.0.6 2015-07-08 CRAN (R 3.2.0)
readxl * 0.1.0 2015-04-14 CRAN (R 3.2.0)
remoji * 0.1.0 2015-11-19 Github (<richfitz/remoji@dc00779>)
tidyr * 0.3.1 2015-09-10 CRAN (R 3.2.0)

🍵 おまけ

明日のR Advent Calendar 2015の担当者は@u_ribo です!!!

(連投になりますが、ご了承ください。3日からは豪華なメンバーによる記事が続きます)。

🔖 出典

このページで利用したデータは、モニタリングサイト1000森林・草原調査 「毎木調査」(環境省生物多様性センター)を加工したものである。(http://www.biodic.go.jp/moni1000/findings/data/index_file.html

📚 参考

  • 市川太祐 (2014). 外部パッケージを用いた集計・整形処理 Rによるモダンな集計処理 『データサイエンティスト養成読本 R活用編』. 技術評論社.
  • Wickham, H., Cook, D., & Hofmann, H. (2015). Visualizing statistical models: Removing the blindfold. Statistical Analysis and Data Mining: the ASA Data Science Journal, 8(4), 203–225. http://doi.org/10.1002/sam.11271
  • Wickham, H. (2014). Tidy data. The Journal of Statistical Software.
  • Horton, N. J., Baumer, B. S., & Wickham, H. (2015). Taking a Chance in the Classroom: Setting the Stage for Data Science: Integration of Data Management Skills in Introductory and Second Courses in Statistics. Chance. http://doi.org/10.1080/09332480.2015.1042739
  • Wickham, H. (2011). The split-apply-combine strategy for data analysis. Journal of Statistical Software.