2017年度版 RStudioを使ったReproducible Research、補足ポエム
この記事はRStudioアドベントカレンダーの10日目の投稿の補足です。私ももう、ゴールしても良いよね、という気になってきました(注)。本体は以前書いた記事で申し訳ないのですが、
になります。古くなったので刷新し、追記をしました。
RによるReproducible Researchの話は@psycle44さんによる5日目の記事でも書かれているのですが、補足として私の所感を述べたいと思います。以下はざっくりとした私の見解です。ポエムです。
ここで扱う、Reproducible Researchの話はRやPythonで記述されたコードにおける再現性の話です。海外では数年前から流れができ始め、現在では市民権を得て大きな流れになっているように感じます。一方日本国内でも少しずつですが認知されつつあり、一部のユーザの間ではかなり浸透しているという印象があります。
R Markdownは研究者ほど使うべき。ぞうさんがいうんだから間違いないですよ。
— kazutan v3.4.3 (@kazutan) June 19, 2016
そもそも再現性とは
学術論文では、調査地や調査対象、時期に始まり、データを得るための詳細や解析に利用したツール、ライブラリとそのバージョンまで書いたりしますよね。良い論文というのはしっかりとした再現性が担保できることだという話も聞いたことがあります。
で、現在ではRやPythonといったスクリプト言語で統計解析やモデリングを実行しますよね。そうなったら、そのコードで書いた内容も「方法」として扱っても良いよね、となるわけです。そうした意味ではバージョン管理をすることは一つの選択肢となり得ると思います。
Git?
Gitはバージョン管理システムの一つで、プログラマやエンジニアの中で使われてきたものですが最近ではデザイナや書籍の編集者も使っています。
Gitの恩恵はいろいろあり、その感じ方も多様ですが、Reproducible Researchでの再現性を保証するのに十分な機能を備えている気がします。Gitを用いる理由は自分自身のためでもあり、共同編集者のためでもあり、第三者のためでもあります。
以前、研究者に対してバージョン管理システムの話をしたことがありましたが、「Dropboxのように意識せずとも使えるようなものでないと使う気にならない」と言われました。Gitは確かにそこまで気がきくものではありません。ですが得られるものは大きいです。ぜひ一度使ってみてほしいものです。またGit以外で他に良いツールがあれば教えてほしいです。
何も全てのGitコマンドに精通する必要はありません。コミット、プッシュ、そしてブランチという概念について少しでもわかっていれば十分にGitを導入することができると思います。私も必要な操作については都度調べたりします。そしてRStudioではこうした最低限必要であろう機能が備わっています。逆に複雑なGitの機能は実装されていません（Shell機能を通して実行することはできます）。
まだRStudio or Gitを使ったことがないという人は今日できることから始めてみてはどうでしょう。まずはコミットで記録していくだけで良いと思います。
Happy Coding with R!
過去の資料とか
オススメ
中級者向けggplot2でこんな図が描きたい - 地図編
どーも。ggplot2は空手の一種として知られているので（要出典）普段の稽古が欠かせまん。今年を振り返り、ggplot2での作図について、いくつかの知見を共有します（書いている余裕がなかったんや...）
library(magrittr) library(jpndistrict)
## Loading required package: jpmesh
## This package provide map data is based on the Digital Map
## 25000(Map Image) published by Geospatial Information Authorityof
## Japan (Approval No.603FY2017 information usage
## <http://www.gsi.go.jp>)
library(sf)
## Linking to GEOS 3.6.1, GDAL 2.1.3, proj.4 4.9.3
library(tidyverse)
## ─ Attaching packages ───────────────── tidyverse 1.2.1 ─
## ✔ ggplot2 2.2.1.9000 ✔ purrr 0.2.4
## ✔ tibble 1.3.4 ✔ dplyr 0.7.4
## ✔ tidyr 0.7.2 ✔ stringr 1.2.0
## ✔ readr 1.1.1 ✔ forcats 0.2.0
## ─ Conflicts ─────────────────── tidyverse_conflicts() ─
## ✖ tidyr::extract() masks magrittr::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ purrr::set_names() masks magrittr::set_names()
library(ggimage)
ggplot2で複数の図を一枚の図にまとめる方法あれこれ
図の中に小さな図（サブプロット）を配置する
こんな図を作りたい。
これは最近あった案件なのですが、やってみると結構手こずりました。また、作業のあとにggimageパッケージの更新があり、こちらの関数を利用するともっと楽になることが判明しました。残念。
自分がやった方法とggimageパッケージの利用例をそれぞれ示します。まずはオレオレから。
地図データは国土交通省 国土数値情報の提供する行政区域データ（国土地理院発行の数値地図（国土基本情報））となります。また、地図と合わせて掲載するデータを環境省 モニタリングサイト1000のページからデータをダウンロードしてきます。
# 北海道のShapefile 国土数値情報 # 行政区域データを利用 sf_pref01 <- jpndistrict::jpn_pref(1, district = FALSE) %>% st_simplify(dTolerance = 0.001)
まずはベースとなる地図を描画します。あとでグラフを追加したいので領域を広げておきます。
p_base <- sf_pref01 %>% ggplot() + geom_sf(fill = "white") + xlim(137, 149) + ylim(41, 46) + theme(panel.border = element_blank(), axis.title = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), plot.background = element_rect(colour = "black"), plot.caption = element_text(size = 6)) + labs(caption = "この地図は、国土地理院長の承認を得て、\n同院発行の数値地図（国土基本情報）電子国土基本図（地図情報）\nを使用したものである (承認番号 平28情使、第603号)")
今回は作図するのが主題なので、データの読み込み等の処理は説明を省略します。ダウンロードしたファイルの読み込みと整形、作図のための関数を書いておきます。
x <- list.files("~/Documents/projects2017/monit1000/data-raw/KOZ02/", pattern = ".xls$", full.names = TRUE) read_temp_xls <- function(path) { readxl::read_excel(path = path, sheet = 2, range = "A2:E10488", col_names = c("datetime", "above_temp", "below_5cm_temp", "below_10cm_temp", "notes"), col_types = c("guess", "text", "text", "text", "text"), skip = 1) %>% select(-below_5cm_temp, -notes) %>% gather(position, temperature, -datetime) %>% filter(!is.na(temperature)) %>% mutate(temperature = as.numeric(temperature)) } plot_temperature <- function(path, ...) { read_temp_xls(path) %>% ggplot(aes(datetime, temperature, color = position)) + geom_line() + ggtitle(...) + theme_set(theme_gray(base_size = 12, base_family = "IPAexGothic")) + theme_transparent() + theme(legend.position = "none", text = element_text(size = 6)) }
これらの関数を使って、一気に作図まで行います。このデータは、 モニタリングサイト1000の高山帯調査で記録された、大雪山の黒岳風衝地における地温・地表面温度です。凡例をつけていませんが、赤い線が地表、青い線が地下10cmでの気温（摂氏）を示しています。
plot_temperature(path = x[1], label = "黒岳風衝地")
ではこの2つの図を組み合わせて行きます。これにはggplotGrob()とannotation_custom()を利用します。
sub_p1 <- ggplotGrob(plot_temperature(path = x[1], label = "黒岳風衝地"))
p_base + annotation_custom(sub_p1, xmin = 136.2, xmax = 140.6, ymin = 43.4, ymax = 45.8)
きちんと理解していないので適当ですが、このように、図中に含ませたい図をggplotGrob()でgtableオブジェクトに変換し、annotation_custom()で上乗せする、という方法です。annotation_custom()の引数で図の大きさ、配置を調整します。
で、これだと図が多くなると手間ですよね。というわけでggimageを使います。私が作業する時には更新が間に合いませんでしたが、今度からはこっちを使いたいと思います。簡単です。
まずは図を用意しましょう。今度は3地点の図を一度に載せます。グラフの描画にもpurrrは便利ですね。
plots <- map2(.x = x[c(1, 5, 9)], .y = c("黒岳風衝地", "黒岳石室", "赤岳第4雪渓"), .f = ~plot_temperature(path = .x, label = .y)) # 個別の図はリストに格納されています # plots[[1]]
ggimageのgeom_subview()はデータフレームに配置する座標や図を格納し他データフレームを展開します。そのためまずは図を配置する座標、サイズ、含めたい図のオブジェクトからなるデータフレームを作りましょう。
# plot変数にlistで図を格納する df_plots <- data_frame(x = c(138.2, 138.2, 146), y = c(45.4, 43.4, 41.8), width = 3.5, height = 1.7, plot = list(plots[[1]], plots[[2]], plots[[3]]))
これをgeom_subview()に渡します。
p_add <- p_base + geom_subview(aes(x = x, y = y, subview = plot, width = width, height = height), data = df_plots) p_add
できたー！格好いい図が描けました。よかったですね。(軸が潰れてしまっていますが、一旦は無視してください)
出典: 「モニタリングサイト1000高山帯調査 地温・地表面温度調査」(環境省生物多様性センター) (http://www.biodic.go.jp/moni1000/findings/data/index_file_earth-temp.html)
この図におけるグラフデータは「モニタリングサイト1000高山帯調査 地温・地表面温度調査」(環境省生物多様性センター)を使用し、瓜生真也が作成・加工したものである。(http://www.biodic.go.jp/moni1000/findings/data/index_file_earth-temp.html)
沖縄を移動、南方諸島を省略した日本地図
さてもう一段。今日はRアドベントカレンダーなので2つ目のtipsです。
日本地図を描画する際、沖縄県や小笠原群島を中心とした南方諸島の島々の位置が気になることがあります。そのため地図の配置を変えるという試みがこれまでされてきているのですが、これのggplot2版というのは見たことがありません。
これもggimage::geom_subview()を使うとできてしまうのでやってみましょう。完成図は上に示したものです。
sf_ja <- 1:47 %>% magrittr::extract(-13) %>% map(~jpndistrict::jpn_pref(pref_code = ., district = FALSE)) %>% reduce(rbind) %>% st_simplify(dTolerance = 0.01) sf_pref13 <- jpn_pref(pref_code = 13, district = TRUE) %>% st_simplify(dTolerance = 0.01) %>% mutate(city_code = as.numeric(city_code)) %>% filter(city_code != 13421) %>% st_union() %>% as.data.frame() %>% mutate(jis_code = "13", prefecture = "東京都") %>% magrittr::set_names(c("geometry", "jis_code", "prefecture")) %>% st_as_sf() sf_ja_omit47 <- sf_ja %>% filter(jis_code != "47") sf_ja_pref47 <- sf_ja %>% filter(jis_code == "47") sf_ja_pref47$geometry %<>% magrittr::add(c(5.6, 17.5)) sf_ja_pref47 %<>% st_set_crs(value = 4326) p <- ggplot(sf_ja_omit47) + geom_sf() + geom_sf(data = sf_pref13, inherit.aes = TRUE) + geom_sf(data = sf_ja_pref47, inherit.aes = TRUE) + geom_segment(aes(x = round(st_bbox(sf_ja_omit47)[1], 0), xend = 132.5, y = 40, yend = 40)) + geom_segment(aes(x = 132.5, xend = 138, y = 40, yend = 42)) + geom_segment(aes(x = 138, xend = 138, y = 42, yend = round(st_bbox(sf_ja_omit47)[4], 0))) + xlab(NULL) + ylab(NULL) + theme(plot.caption = element_text(size = 6)) + labs(caption = "この地図は、国土地理院長の承認を得て、\n同院発行の数値地図（国土基本情報）電子国土基本図（地図情報）\nを使用したものである (承認番号 平28情使、第603号)") p
いいですね。Enjoy! 時間がなくてコードの説明が足りていないので、後ほど追記します。
本日発表!ほくぽえむ大賞2017 俳句の部
ホクソエムといえばポエムです（要出典）。今日はTwitterでのホクソエム氏の投稿から、俳句を探してこようと思います。そして、今年のベスト俳句を独断と偏見により決めます。戦略としては、Twitterから投稿を取得、日本語形態素解析システム JUMAN++により形態素解析、単語の読みの文字数から俳句の定型を判定、というものです。
それでは早速データを取得するところから始めましょう。データの利用に際して、ホクソエム氏から許可をいただいております。
どういうことか良く分からないけど好きにしたらいい
— h(o x o\_)m＜ Fool on the Hill (@hoxo\_m) December 6, 2017
言葉にするのって難しいですね。
データ取得と前処理
ツイッターのAPIを利用するため、rtweetパッケージを用います。認証を済ませ、タイムラインから投稿を取得します。
library(magrittr) library(tidyverse) library(rtweet)
get_timeline()でユーザを指定し、取得件数を決定します。他の引数からAPIのパラメータを変更可能なようですが、下記の例ではうまくいかなかったです。
df_hoxom <- get_timeline("hoxo_m", n = 3000, exclude_replies = TRUE, include_rts = FALSE)
取得したデータから、今年の投稿を抽出します。また、RTやリプライ、URLを含んだ投稿等を除外します。加えて今回の俳句判定には、アルファベットを含むとカウントを正常に行えないという課題がありましたのでこれも除外します。その他の空白や記号の除去も次のコードで実行します。
df_hoxom$created_at %>% range()
## [1] "2016-09-24 14:21:10 UTC" "2017-12-06 13:15:42 UTC"
df_hoxom %<>% # 今年の投稿、RTやメンションを取り除く filter(between(created_at, lubridate::ymd_hms("20170101 00:00:00"), lubridate::ymd_hms("20171231 23:59:59")), is_retweet == FALSE, is.na(reply_to_status_id), is.na(urls_url), !grepl("^@", text)) %>% # 記号、空白の除去 mutate(text = str_replace_all(text, "[:punct:]", "") %>% str_replace_all(., "[:space:]", "")) %>% filter(!grepl("[A-Za-z]", text)) %>% # 余分な列を削除 # 投稿のID、投稿日時、投稿文に制限 select(status_id, created_at, text)
俳句の判定
肝となる俳句の判定は、次の条件で行います。
- 投稿文を平仮名読みに直し、文字数を数える
- 単語の文字数で5,7,5のリズム(17文字)になるものを「俳句」とする
"古池や蛙飛び込む水の音" %>% str_count()
## [1] 11
"ふるいけやかわずとびこむみずのおと" %>% str_count()
## [1] 17
df_hoxom$text[1]
## [1] "謙虚さと学び"
df_hoxom$text[1] %>% str_count()
## [1] 6
「柿食えば鐘が鳴るなり法隆寺」のように17文字ではないものも俳句になりますが、今回は厳密に17文字のものを俳句とみなしています。
肝心の自然言語処理の部分はJUMAN++にやってもらいます。JUMAN++のRラッパパッケージrjumanppをid:songcunyouzai (y__mattu) が開発してくれているのでそれを使います。サンキューマッツ！
[https://github.com/ymattu/rjumanpp:embed:cite]
[http://y-mattu.hatenablog.com/entry/2017/08/19/230432:embed:cite]
text_wakati <- rjumanpp::jum_text("かずたんがフリー素材になっている") text_wakati[[1]]
## [1] "かず" "かず"
## [3] "かず" "名詞"
## [5] "6" "普通名詞"
## [7] "1" "*"
## [9] "0" "*"
## [11] "0" "\"代表表記:下図/かず"
## [13] "自動獲得:EN_Wiktionary\""
text_wakati[[5]]
## [1] "素材" "そざい"
## [3] "素材" "名詞"
## [5] "6" "普通名詞"
## [7] "1" "*"
## [9] "0" "*"
## [11] "0" "\"代表表記:素材/そざい"
## [13] "〜を〜に構成語" "カテゴリ:人工物-その他\""
「かずたん」（もしかして: @kazutan）がうまく分かち書きできていないのが気になりますが、読みや品詞の区分が行われた結果が得られました。
この結果をデータフレームにしておくと色々捗るので関数を用意しましょう。
jum_text_separate <- function(text) { x <- rjumanpp::jum_text(text) res <- x %>% purrr::map_df(~tibble::data_frame(.[2], .[4], .[6])) %>% purrr::set_names(c("yomi", "hinsi_bunrui_dai", "hinsi_bunrui_sai")) return(res) }
次に今作成した関数を利用し、読みの文字数が17文字となるものを抽出します。これが条件1をクリアした「俳句」の候補となります。そして、条件をクリアした投稿文に対して、単語の読みの文字数を数えます。これを足し合わせ、5,7(12),5(17)となっているものを最終的な「俳句」とみなします。
"かずたんがフリー素材になっている" %>% jum_text_separate(text = .) %>% # 17文字となるものに制限 filter_at(vars("yomi"), any_vars(str_count(paste(., collapse = "")) == 17)) %>% # 単語の文字数を累積 mutate(haiku = cumsum(str_count(yomi)))
## # A tibble: 8 x 4
## yomi hinsi_bunrui_dai hinsi_bunrui_sai haiku
## <chr> <chr> <chr> <int>
## 1 かず 名詞 普通名詞 2
## 2 たん 名詞 普通名詞 4
## 3 が 助詞 格助詞 5
## 4 ふりー 形容詞 * 8
## 5 そざい 名詞 普通名詞 11
## 6 に 助詞 格助詞 12
## 7 なって 動詞 * 15
## 8 いる 接尾辞 動詞性接尾辞 17
この「かずたん〜」は5,7,5になっているので条件をクリアしていることになります。ということでこれも関数化しておきましょう。先ほど書いた関数を追加して、一度に全ての処理を完結するようにします。
is_haiku <- function(text, ...) { x <- jum_text_separate(text) %>% filter_at(vars("yomi"), any_vars(str_count(paste(., collapse = "")) == 17)) %>% mutate(haiku = cumsum(str_count(yomi))) %>% filter(haiku %in% c(5, 12, 17)) res <- if_else(nrow(x) == 3, TRUE, FALSE) return(res) }
次のような判定結果を得ます。
c("かずたんがフリー素材になっている", "謙虚さと学び", "なつくさやつわものどもがゆめのあと") %>% map_lgl(is_haiku)
## [1] TRUE FALSE TRUE
is_haiku()の結果を投稿データに適用し、いよいよ「ほくぽえむ大賞2017 俳句の部」の作品発表です！
df_hoxom %<>% filter(pmap_lgl(., ~is_haiku(text, ...)) == TRUE)
ほくぽえむ大賞2017 俳句の部
入選作品の発表です。厳しい審査の結果、今年は4件の発表が選ばれました。
df_hoxom %>% knitr::kable(format = "markdown")
| status_id | created_at | text |
|---|---|---|
| 917606569378017280 | 2017-10-10 04:23:49 | ぞうさんのおすすめテーマにしてみた |
| 874064391338999808 | 2017-06-12 00:42:44 | エレベタやおっさんどもが臭の跡 |
| 852532562635313152 | 2017-04-13 14:42:56 | 多様体いろんなとこに潜んでる |
| 832896107428536320 | 2017-02-18 10:14:41 | かずたんがフリー素材になっている |
佳作
かずたんがフリー素材になっている。
— h(o x o_)m＜ Fool on the Hill (@hoxo_m) February 18, 2017
どういう状況かわかりませんが、声に出して読みたい感が高評価につながりました。おめでとうございます！また、「かずたん」は入選作品に2回も登場することから俳句に適したポエットな単語であることがわかりました。
多様体いろんなとこに潜んでる
— h(o x o_)m＜ Fool on the Hill (@hoxo_m) April 13, 2017
データサイエンティストとしての意識からなのか、ふとしたつぶやきに才能を感じるこの投稿が佳作です。わたしはこの作品から、星新一のショートショートを連想しました。皆さんはどうでしょう。
大賞
エレベタや おっさんどもが 臭の跡
— h(o x o_)m＜ Fool on the Hill (@hoxo_m) June 12, 2017
こちらの俳句は初夏に詠まれたものです。「エレベタ」と「おっさん」の組み合わせ、そして「息の跡」とポエム度が高いです！スペースも入れて、詠み人自身の俳句としての意識が強いこの投稿が大賞となります！おめでとうございます！！
審査委員の言葉
いかがだったでしょう。思ったりもポエム（俳句）が少ないですね。17文字という自由度の低さが原因だったのか、この辺は改良が必要そうです。また、来年は「俳句」で肝心なキレも採点基準に加えたいと思います。
来年はより多くのポエムが投稿されることをねがいます。ホクソエムアドベントカレンダー、明日は未定です！誰か〜
Enjoy!
