cucumber flesh

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

質問お待ちしております！できればreprex使ってね 😸

tidyverse R team

最近、チームで分析プロジェクトを進行していることもあり、人のコードをみてコメントしたり、自分もみてもらう、というやりとりが増えてきました。メンバーの中で私は、どちらかといえば「R言語チョットデキル」人間で、時々発生するトラブルや質問について答えることもあります。

よくあるQ&A時の話ですが、回答に答える側になってみると質問してくる人に対してあれこれ言いたいことが出てきました。そんな気持ちを詠んだポエムです。Rの話になりますが、色々な言語で、強いてはプログラミング以外のことにでも言えるかもしれません。

Help me help you

一番に言いたいのは「 私はあなたの問題を解決したいと思っている 」ということ。個人的には、質問は自分の思いつかなかった手法や問題を提示してくれるので歓迎です(一プログラマとしての成長が望まれるので私にも見返りがある。感謝の気持ちもあります)。質問をためらわないでください（もちろん自力で解決できるに越したことはないし、落ち着いて見直したら・再度実行したら解決するもの、調べて済む問題は除く）。

問題に直面したあなたは、やや混乱しているかもしれません。ですが落ち着いて、問題がどこにあるのかを確認してみてください。文字の間違い、括弧や引用符の閉じ忘れ、などはよくある話です。そして質問してください。

あなたからの質問

注: ここでは、メールではなくてSlackのチャットやGitHub issuesでの質問を想定します。

質問はコードと一緒に

あなたからの質問が来ました。

「こんなエラーが出た」「これができない」

... お、おう。

これではちょっと困ってしまいます。あなたが何をしたいのか、何をしたら問題が発生してしまったのか、どこでつまづいているのか、詳細な情報とコードと共に示してくれると助かります。

個人的に、世の中、説明不足な質問で溢れているように感じます。説明が足りていないと、状況を把握しにくい、問題を突き詰めにくい、など悪いことが多いです。過多であることよりも不足していることが困ります。あれそれについての情報・状況を知らせてくださいと、あーだこーだ何度もやりとりを行うのはお互い良くないでしょう。最小限で問題が解決させるためには、あなたからの最初の情報提供が肝心です。

ありがちなパターン

  • 問題の途中から始めてしまう... 問題の箇所だけ、というのも困ります。ここでエラーになった、その情報は大事ですが、実は問題はもっと以前から、根本的な箇所で生じている可能性もあるのです。

コードはテキストで（ハイライトさせてくれるとさらに嬉しい）

次にあなたから実行したい処理の内容と問題の箇所を示したコードが示されました。みてみましょう。... ちょっと待って、もしかしてコンピュータの画面をスマートフォンで撮影した画像を貼り付けました? スクリーンキャプチャにしました? どうしてそんなことを!!!

コードを画像で撮影するのが好ましくないのは、画像をコピペしても問題の処理を再現することができないためです。再現、それは質問に答えるためにとても大事なことです。私が例え長年の経験と勘に優れたプログラマであったとしても、適当なコメントだけで質問に答えようとはしないでしょう。自分の考えが正しいかを確認するために、コードを自分の環境で試したいのです。しかし画像では、画像とエディタを相互に見比べながら画像に書かれている「コード」を書き起こす必要が生じてしまいます。ですので、コードは画像や文章ではなくて「コード」で示してください。ソースコードをファイルに保存してメールで添付してくれても構いません。

ソースコードですが、コピペで貼り付けました?それでも良い... い、いや、問題を説明している他の文章とは区別がつくようにしてください。お願いします。

ソースコードを記述するプログラミング言語は、数値と記号、文字列で表現されます。普通にみると見辛いです。これはあなたと一緒です。私たちはプログラミング言語を扱います。プログラミング言語が自然言語と異なるのはあなたも知っていることでしょうが、SlackやGitHubではMarkdownと呼ばれるマークアップ言語（簡単な文字装飾機能と思ってください）が利用でき、その中にコードを他の文と区別させるための機能があります。

Creating and highlighting code blocks - User Documentation

加えて、シンタックスハイライトという文字に色をつける機能が用意されていて、これは多くの言語で利用できます。ソースコードの中の役割のようなものごとに異なる色が与えられ、その色を識別することで'なんとなく'読めるようになっているのです。普段シンタックスハイライトに慣れていると、装飾されていないコードを読むのは正直辛いのです。

Markdownでは` を三回繰り返しで示し、同じく ` の三回繰り返しで囲んだ部分がコードブロックとして扱われます。コードブロックをハイライトさせるために例えばRであれば次のようにします。

```r
library(reprex)

x <- mean(c(1, NA, 3:4))
mean(x)
# [1] NA
```

pythonであればrの部分をpythonとしてください。このブログにもシンタックスハイライトは実装されていて、上記のコードを鮮やかにすることができます。

library(reprex)

x <- mean(c(1, NA, 3:4))
mean(x)
# [1] NA

色をつける。単純なことですがこれなら少し読む気になりません?

reprexパッケージを使う

さて、いよいよRならではの話に入ってきます。これまでに説明した、質問へのお願いを実現するために reprex パッケージが便利です。パッケージの機能を説明した id:yutannihilation さんの記事がすでにありますが、ちょっと内容の追記と更新をしておきます。なおreprexは rep roducible ex ample からなる造語で、Romain Francoisによって名付けられました。

notchained.hatenablog.com

Prepare Reproducible Example Code for Sharing • reprex

reprex はGitHub issuesやStackOverflow、Slackでの再現可能なコードの例を用意してくれるパッケージおよびその関数名です。現在CRANに登録されているバージョンは0.1.2で、ここでもそのバージョンを利用しています。

reprex()ないで対象のコードを次のいずれかの形式で与えるか、ファイルを指定するかで、一時フォルダの中でコードが実行されます。RStudioを利用しているのであれば、専用のアドインも利用可能です（個人的にはこちらがおすすめ）。

reprex(mean(rnorm(10)))

reprex(input = "mean(rnorm(10))\n")

reprex(input = "my_reprex.R")

ここで肝心なのは、reprex()を実行している環境とは別の場所でコードが再実行されることと、コピーペーストでGitHub等のサービスにシンタックスハイライト可能なMarkdown形式で貼り付け可能なテキストが用意されることかと思います。

reprex()は、質問への回答者がコピーペーストで問題を再現できるように調整されています。ですのでreprex()は、対象のコードを実行して、読み込まれていないパッケージの関数やオブジェクトがある場合、出力もエラーメッセージを含んだものになります。これは実行したコードでは示されていないパッケージの読み込みやオブジェクトの作成がコードとは別の場所（例えばコンソールで実行しただけでコードには書かれていない）で行われていたことに起因します。人間のコピペだと、うっかり問題の箇所を見落としたり、情報を不足させる危険があるので、reprex()での実行環境を変えた再実行が効果的なのです。

ここで注意なのが、csvデータ等の読み込みです。read.csv()などを伴う処理を行う場合、パスを明確に指定しないとファイルがない、とエラーになることがあります。例えば現在の作業ディレクトリ（Rprojファイルがあるディレクトリ）で次のコードをreprex()に渡すと再現が失敗します。

d <- read.csv("iris.csv")
table(d$Species)

d <- read.csv("iris.csv")
#> Warning in file(file, "rt"): cannot open file 'iris.csv': No such file or
#> directory
#> Error in file(file, "rt"): cannot open the connection
table(d$Species)
#> Error in table(d$Species): object 'd' not found

ファイルの入出力を伴う場合はパスはフルパスで示すようにしましょう。

d <- read.csv("/home/rstudio/iris.csv")
table(d$Species)
#>
#>     setosa versicolor  virginica
#>         50         50         50

問題がデータではなくて処理にある場合、データはRに最初から用意されているデータやダミーデータを作成してくれると助かります（ローカルファイルの結果を示すことはできますが、再現にはファイルそのものが必要になるためです）。reprex公式では、irisやmtcarsが良いとされています。またデータはすべてでなくて構わない時があります。この時はhead()やsample()を使って一部のデータを示してくれるだけで良いです。sample()などランダムに結果が異なる処理では、set.seed()により乱数の固定をすることも再現性の確保に必要になります。

Session Information

ソースコードの中身を追ったところ、コードには問題がないようです。ではどうして問題が発生しているのでしょう?

問題を深追いするために、あなたの利用環境について教えてもらう必要があります。RではsessionInfo()が用意されていて、これを実行するとコンピュータのOSや読み込んでいるパッケージの情報が出力されます。

sessionInfo()
#> R version 3.4.3 (2017-11-30)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Debian GNU/Linux 9 (stretch)
#>
#> Matrix products: default
#> BLAS: /usr/lib/openblas-base/libblas.so.3
#> LAPACK: /usr/lib/libopenblasp-r0.2.19.so
#>
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=C             
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#>
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#>
#> loaded via a namespace (and not attached):
#>  [1] compiler_3.4.3  backports_1.1.2 magrittr_1.5    rprojroot_1.3-2
#>  [5] tools_3.4.3     htmltools_0.3.6 yaml_2.1.16     Rcpp_0.12.15   
#>  [9] stringi_1.1.6   rmarkdown_1.8   knitr_1.19      stringr_1.2.0  
#> [13] digest_0.6.15   evaluate_0.10.1

パッケージの違いを追うのにこの情報が必要なことがあります。さらに専用のパッケージ、 sessioninfo が用意されていて、こっちを使ってもらっても良いです（こっちの方が良いかな）。sessioninfo の出力の例です。

sessioninfo::session_info()
#> ─ Session info ──────────────────────────────────────────────────────────
#>  setting  value                       
#>  version  R version 3.4.3 (2017-11-30)
#>  os       Debian GNU/Linux 9 (stretch)
#>  system   x86_64, linux-gnu           
#>  ui       X11                         
#>  language (EN)                        
#>  collate  en_US.UTF-8                 
#>  tz       UTC                         
#>  date     2018-02-28                  
#>
#> ─ Packages ──────────────────────────────────────────────────────────────
#>  package     * version    date       source                      
#>  backports     1.1.2      2017-12-13 CRAN (R 3.4.3)              
#>  clisymbols    1.2.0      2017-05-21 CRAN (R 3.4.3)              
#>  digest        0.6.15     2018-01-28 CRAN (R 3.4.3)              
#>  evaluate      0.10.1     2017-06-24 CRAN (R 3.4.3)              
#>  htmltools     0.3.6      2017-04-28 CRAN (R 3.4.3)              
#>  knitr         1.19       2018-01-29 CRAN (R 3.4.3)              
#>  magrittr      1.5        2014-11-22 CRAN (R 3.4.3)              
#>  Rcpp          0.12.15    2018-01-20 CRAN (R 3.4.3)              
#>  rmarkdown     1.8        2017-11-17 CRAN (R 3.4.3)              
#>  rprojroot     1.3-2      2018-01-03 CRAN (R 3.4.3)              
#>  sessioninfo   1.0.0      2017-06-21 CRAN (R 3.4.3)              
#>  stringi       1.1.6      2017-11-17 CRAN (R 3.4.3)              
#>  stringr       1.2.0      2017-02-18 CRAN (R 3.4.3)              
#>  withr         2.1.1.9000 2018-02-28 Github (r-lib/withr@5d05571)
#>  yaml          2.1.16     2017-12-12 CRAN (R 3.4.3)

sessionInfo()を含めないのはなぜ?

id:yutannihilation さんの記事の中で

session infoは含めない！

と書かれています。うーむ。これについて、私も最初ユタニさん同様、どうしてだろうと思いました。

ですが、sessionInfo()の情報が必要になるのは、問題がコードにない時だろうと思います。問題がコードにあるのであればsessionInfo()は余分な情報となります。

sessionInfo()はコードだけで問題解決ができない時、追加の情報として求めれば良いのだろうと考えています。

reprex 0.2.0 が良い感じ

reprex パッケージの最新版は 0.1.2であると述べましたが、おそらく次期バージョンの 0.2.0 ではさらに使いやすく、より便利になっています。問題がなければこちらのバージョンを先取りし、GitHubからインストールして使うのを勧めます。

インストールは次のコードの実行で行われます。

install.packages("githubinstall")
githubinstall::githubinstall("reprex")

0.2.0では次の画像のような出力になります。reprexを使ったこととそのバージョンも示されていて良いですね。

f:id:u_ribo:20180301073429p:plain

ちなみにパッケージは現在は tidyverseリポジトリに属するようになっています。

それでは、良い質問と良い回答を! Enjoy!

(あとでちょっと追記するかも)

最近お気に入りのRStudio(ver.1.1.423)の機能: ファイル内検索とTODOハイライト

rstudio tips

先日、RStudioのバージョン1.1.423がリリースされました。

v1.1からの更新内容については先日開催されたrstudio::conf 2018にて、開発者の一人であるKevinが発表しているので詳細はこちらをご覧ください。国内ではぞうさんこと@kazutanが記事にしています。

RStudio v1.1 – What’s New

kazutan.github.io

そんな開発が活発なRStudioですが、最新のv1.1.423で利用可能な個人的におすすめしたいプロジェクトの運用がしやすくなる2つの機能について紹介します。

プロジェクト中のファイル内検索

どうしてこの機能に気がつかなかったのだろう。この機能を使う前と後では私のRStudioの使い方は一変してしまいました、そのくらいのお気に入りの機能。

プロジェクトやパッケージが肥大してくると、全てのRファイルや関数について把握するのは難しくなります。「このデータファイルはここで使っている」、「この図はあそこのRファイルで出力している」などなど。対象のファイルの名前やパス、オブジェクト名を変更した際、影響を及ぼす範囲を把握しきれていないと変更が追いつかず、ある日エラーになっている、という経験が皆さんにもあるでしょう。

そんな時にこの機能を使うと、Rプロジェクト中のファイルから、任意の文字列を検索し、どのファイル、どの部分で記述されているかがわかって大変便利です。特に誰かのパッケージやプロジェクトに首を突っ込む時、自分が変更した箇所が他のファイルに反映されているかを確かめるときに重宝します。

ショートカットキーが当てられていて、macOSでは「コマンド + シフト + F」で検索用の画面が起動します。正規表現の利用や対象のパス、ファイル形式について設定できて高機能です。

f:id:u_ribo:20180221071044g:plain

検索結果は、RStudioのコンソールパネルのタブとして表示されます。検索結果を見ると、検索した文字列がどこで使われているかがハイライトされます。また、特定行をクリックすることで、その部分をソースパネルで編集することができます。

全てのファイルの中身を一度に変更することはできませんが、こうして影響のある部分を書き換えていくことで、変更によるトラブルを防げるようになります。

TODOハイライト

こちらはid:niszetさんのつぶやきで知りましたが、Rファイル中で特定の文字列がハイライトされるようになっています。特定の文字列とは"TODO", "FIXME"の2つです。

github.com

コード中にこうしたコメントを残しておくことで、視認性・コードの保守性が向上しますね。

f:id:u_ribo:20180221073136p:plain

todorパッケージ

github.com

この機能と組み合わせて、todorパッケージが便利です。このパッケージは、プロジェクト中の"TODO"や"FIXME"などの特定の文字列を検索します。こちらが対応しているのはRStudioのものよりも多いです。

CRANには登録されていないため、利用の際にはGitHub経由でインストールを行います。

remotes::install_github("dokato/todor")

f:id:u_ribo:20180221073343g:plain

todor::todor()およびRStudio Addinsから実行します。

RStudio v1.2に向けて

GitHubのNEWSを見ると、どうやら次のバージョン(v1.2)ではreticulateパッケージによるPythonとの連携機能の強化が行われている様子。

また、daily build版だけなのかもしれませんが、アイコンが...

Enjoy!

RStanのおさらいをしながら読む 岩波DS1 (ややアップデートした2018年版)

skimr stan modeling ggplot2

最近、以前書いた岩波DSの久保さんパートのおさらい記事を久保さん本人がRTしていたのを見て（ああ、あの頃はまだ大学院生をやっていたのだと懐かしくなったと同時に）、移り変わりのあるRコードは今も動くのか心配になった。結論を先に言うと、RPubsに書いたコードは現在の環境でも問題なく動作した。が、その頃にはtidyverseも十分に整備されていなかったし、私のR技術力も低く、ヘンテコなコードを書いていたのが気になる。そこで、新しいパッケージも出ていることだし、2018年版ということで再度実行してみた、というのがこの記事の顛末である。

元の記事でもstanコードは id:StatModeling さんの記事に書かれているものをコピーしたが、代入に「<-」を使うとstanが警告文を出すので「=」に書き換えておいた（非推奨なだけで実行自体はしてくれる。<-が非推奨〜云々は「StanとRでベイズ統計モデリング」でも書かれている）。

StanとRでベイズ統計モデリング (Wonderful R)

StanとRでベイズ統計モデリング (Wonderful R)

以下、変更した部分だけを書く。結果は変わらないので内容ではない。

まず例題として用意されているデータ(久保さんのウェブサイトからダウンロード)をdplyr::glimepse()で確認しているが、これと同じデータの確認は、私は最近skimrパッケージのskim()を使ってやることが多い。

skimr::skim(d)
# Skim summary statistics
#  n obs: 20 
#  n variables: 6 
# 
# Variable type: factor 
#  variable missing complete  n n_unique             top_counts ordered
#      pref       0       20 20       10 A: 2, B: 2, C: 2, D: 2   FALSE
# 
# Variable type: integer 
#  variable missing complete  n mean   sd p0   p25 median p75 p100     hist
#         N       0       20 20 54.1 2.77 49 52.75   54    56   59 ▃▃▂▇▆▃▂▅
#         X       0       20 20  0.5 0.51  0  0       0.5   1    1 ▇▁▁▁▁▁▁▇
# 
# Variable type: numeric 
#  variable missing complete  n   mean   sd     p0    p25 median    p75   p100     hist
#       Age       0       20 20   0.5  0.51   0      0      0.5    1      1    ▇▁▁▁▁▁▁▇
#    mean.Y       0       20 20 155.3  2.8  151.36 153.06 155.69 157.22 161.71 ▃▇▁▂▇▁▁▁
#      sd.Y       0       20 20   2.97 0.23   2.45   2.9    3.06   3.11   3.21 ▂▁▂▁▁▂▇▆

これだけでも数値からなる変数のデータのばらつきが把握できて大変良い。なおskimrパッケージの良さは id:niszet 氏が語ってくれている。

niszet.hatenablog.com

skimrはさらにここからdplyrパッケージの関数と組み合わせて利用することでグループ化したデータの要約や並びかえも行える。例えば、都道府県 pref 、測定回数 Age (1,2回目の測定をそれぞれ0と1で示す）ごとの平均身長 mean.Yの中央値を求め、平均身長の順に並びかえて表示するには次のようにする。dplyr(あるいはtidyverse)は、モデリングを行う際でも使う機会があるので最初に読み込んでおくと良いだろう。

library(dplyr)

# skimr::skim()以外はすべてdplyrパッケージのデータ操作関数を利用
d %>%
  # 必要な変数にのみ制限
  select(pref, Age, mean.Y) %>%
  # 都道府県、測定回数でグループ化し
  group_by(pref, Age) %>%
  # 統計量を求める
  skimr::skim() %>% 
  ungroup() %>% 
  # 統計量から中央値を抽出
  filter(stat == "median") %>% 
  # 平均身長の中央値が低い順に並びをかえる
  arrange(value)
# # A tibble: 20 x 8
#    pref    Age variable type    stat   level value formatted
#    <fct> <dbl> <chr>    <chr>   <chr>  <chr> <dbl> <chr>    
#  1 A      0    mean.Y   numeric median .all    151 151.36   
#  2 B      0    mean.Y   numeric median .all    152 151.56   
#  3 C      0    mean.Y   numeric median .all    152 152.22   
#  4 I      0    mean.Y   numeric median .all    153 152.67   
#  5 G      0    mean.Y   numeric median .all    153 152.98   
#  6 D      0    mean.Y   numeric median .all    153 153.09   
#  7 E      0    mean.Y   numeric median .all    153 153.22   
#  8 H      0    mean.Y   numeric median .all    153 153.27   
#  9 F      0    mean.Y   numeric median .all    153 153.31   
# 10 J      0    mean.Y   numeric median .all    155 155.37   
# 11 D      1.00 mean.Y   numeric median .all    156 156      
# 12 I      1.00 mean.Y   numeric median .all    157 156.82   
# 13 B      1.00 mean.Y   numeric median .all    157 156.83   
# 14 C      1.00 mean.Y   numeric median .all    157 157.08   
# 15 F      1.00 mean.Y   numeric median .all    157 157.22   
# 16 E      1.00 mean.Y   numeric median .all    157 157.24   
# 17 A      1.00 mean.Y   numeric median .all    157 157.27   
# 18 G      1.00 mean.Y   numeric median .all    158 157.81   
# 19 H      1.00 mean.Y   numeric median .all    159 158.95   
# 20 J      1.00 mean.Y   numeric median .all    162 161.71

glm()の部分は依然としてbroomパッケージを使っているので変更なし。

さて以前はstanコードのdataブロックに渡すデータをリストで作成したが、ここも変わりがない。強いて言うなら、 list()の代わりにtibble::lst()にすることや、ユニークな都道府県の数をカウントした値（length(unique(d$pref))）はよく使うので、別なオブジェクトにしておくと良いだろう。また、length(unique(x))の部分はn_distinct(x)に書き換えるのが新しいやり方。もしくはこうした定数は別なファイルに記述して、configパッケージで必要な時に参照する方法もある。

n_pref <- n_distinct(d$pref)

# list
tibble::lst(
  N_r     = 2,
  N_pref  = n_pref,
  Mean_Y  = t(matrix(d$mean.Y, n_pref, 2)),
  Sigma_Y = t(matrix(d$sd.Y / sqrt(d$N), n_pref, 2)),
  X       = t(matrix(d$X, n_pref, 2)),
  Age     = t(matrix(d$Age, n_pref, 2)))

suryu.me

blog.hoxo-m.com

stanコードのビルド、MCMCサンプリングの実行も変わりがない。記事を書いた時からStanは色々なアップデートがあっただろうが、動いてよかった..。

最後にサンプリング結果から推定したパラメータの可視化を行なっているが、ここもそのままのコードが動作する。ただしgridExtraを使った複数のggplot2作図の結合はちょっと古いかつ、若干わかりにくいのでeggあるいはpatchwork (2018年2月現在 GitHub版のみ利用可能 https://github.com/thomasp85/patchwork) パッケージで代替する。eggの方はgridExtraの開発者が中心に開発していて正式な後継っぽいのだが、個人的にはpatchworkの方が好み。

事後分布のヒストグラムやら、ggmcmcで作図するところまでは同じ。まずはeggを使った方法から。

library(egg)
(p <- ggarrange(p1, p2, p3, p4, widths = 1:2))

library(patchwork)
(p <- p1 + p2 + p3 + p4 + plot_layout(ncol = 2))

# egg, patchwork いずれも保存は同じ方法
ggsave(filename = "samping_parameters_result.png",
       p, 
       width = 6, height = 4)

f:id:u_ribo:20180213002849p:plain

ああ、まだ「ベイズモデリングの世界」が買ったまま読めていない...。

他にも、私が追えていないだけで更新できる内容があるかもしれないが、2018年版ということでひとまずは Enjoy!

ベイズモデリングの世界

ベイズモデリングの世界

Rで楽しむベイズ統計入門[しくみから理解するベイズ推定の基礎] (Data Science Library)

Rで楽しむベイズ統計入門[しくみから理解するベイズ推定の基礎] (Data Science Library)

（こっちもよろしくね😼）