2017-09-24

第1回 Kobe.R annex（犬4匹本輪読会）開催しました！

2017年9月24日大阪府立大学・なんばサテライトにて Kobe.R annex 第1回を開催しました。第1部『ベイズ統計モデリング―R,JAGS,Stanによるチュートリアル― 原著第2版』の輪読，第2部分析ネタなんでもという構成でした。参加者は5名（＋会場をご提供いただいた大阪府立大学の先生）。

f:id:hikaru1122:20170924191307j:plain

　
予定は第6章〜9章の予定だったのですが，8章までとしました。もちろん全員，犬4匹本を購入済みです。以下，輪読会の簡単メモです。

第7章

JAGS の情報が載っていてわかりやすいウェブサイトがある。

http://ecological-stats.com/tag/jags/

MCMCスライド集

hikaru1122.hatenadiary.jp

stan のモデルをコンパイルするとき次のエラーが出たら，Xcode を1度起動して使用許諾に同意する（Xcode のアップデートをした直後に出る恐れ）

追加情報:  警告メッセージ: 
 命令 ''/Library/Frameworks/R.framework/Resources/bin/R' CMD config CXX 2>/dev/null' の実行は状態 69 を持ちました

追加情報:  警告メッセージ: 
 命令 '/Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB file8dae4ce1404e.cpp 2> file8dae4ce1404e.cpp.err.txt' の実行は状態 1 を持ちました

疑問点

Gelman-Rubin統計量（184頁）はRハットのことだろうか？

→どうやらそのようです。

間引き（thinning）はすべき？　すべきでない？

→犬4匹本は基本thinningしない的な立場（191頁）のようだが，thinning したほうがよいという文献もあるらしい。

thinning をした場合，サンプリング時間はどう変わるだろうか？

→あるモデルで「iteration 1000，thin 1」「iteration 1000，thin 2」「iteration 2000，thin 2」を stan で実行してみると，90秒・90秒・180秒となった。保存されたオブジェクトのデータは8MB, 4MB, 8MB。thinning の設定でサンプリング時間は変わらないが，保存されるオブジェクトの容量は変わる。

第8章

runjagsのほうが早くて便利！？
統計学エンドユーザーのための JAGS 入門

hikaru1122.hatenadiary.jp

参加者の方が8章の便利関数をパッケージ化されました！

参加させていただき、ありがとうございました。

8章の便利コマンドですが、パッケージとしてgithubに上げました。
devtools::install_github("kyusque/DBDA2E-utilities", args = "--no-multiarch")
— ヨゼ (@kyusque) 2017年9月24日

LT枠（分析ネタなんでも）

状態空間モデルについて
共変量バランシング傾向スコアを使った統計的因果推論について

　
以上です。次回もいずれ開催したいと思います。

2017-09-20

大きい iPad Pro で論文を快適に読む方法

研究に役立つツール

自分の中でだいたい固めることができたので、大きい iPad Pro で快適に論文を読む方法を書いておきます。最近、質的研究ばかりやっていて、R とかあんまり触っていないため、こういうネタしかないんです。

結論から

　
まず先に結論を言うと、2017年9月20日現在、次の組み合わせが私にはピッタリです。

　
1. iPad Pro 12.9インチ
2. アップルペンシル
3. アップルペンシルのキャップ紛失防止ツール
4. 適度に滑りにくくするための保護シート
5. PDF Expert
6. スリーブケース

私の電子リーダー端末遍歴

　
タブレットやスマホはいろいろ触ってきましたが、論文を読むに関しては、Retina ディスプレイの iPad（第3世代）→ソニーのデジタルペーパー（旧モデル・DPT-S1）→ 小さい iPad Pro（9.7インチ・2015年モデル）→大きい iPad Pro（12.9インチ・2017年モデル）とたどってきました。他にもKindleとかソニーのリーダーとかでもチャレンジしましたが、目には優しいものの動作が重くダメでした。

なぜ紙で論文を読まないのか？

　
私が紙で論文を読みたくないと思ったのは単純に「かさばるから」です。仕事柄出張が多く（昨年は224泊）、論文をたくさん持ち歩きたくありません。社会科学の論文は（少なくとも私の領域では）だいたい1本20ページくらいで、5本でA４・100ページ、10本で200ページ。出張で持ち歩くには思ったよりかさばるんです…。そのため電子化しかないと数年前からいろいろ工夫していました。

　
中には論文は「紙で読む派」という方もいらっしゃると思いますが、差し迫った問題がなければ、私はどちらでもいいと思います。お好きなほうを選べばいいです。

デジタルペーパー VS 小さいiPad Pro VS 大きい iPad Pro

　
電子化された論文を読む勝ちパターンがこの1年で固まりました。まず昨年、ソニーデジタルペーパー（DPT-S1）を大阪のソニーストアで実機を確認。さすがに定価で買うには抵抗があり、中古をヤフオクでゲット。それでもそれなりのお値段がしました。

デジタルペーパーのいいところ

○　軽い、薄い。
○　なんか目に優しい感じがする。
○　基本、読み書きに特化しているので、論文読みに集中できる。

デジタルペーパーのイマイチだったところ

△　動作がもったり（新モデルは改善しているそうです）
△　文字がギザギザで読みにくい（これも新モデルでは改善しているそうです）
△　思ったよりペンが引っかかる。もっと滑らかなのが個人的には好き。
△　英語論文を読んでいて、知らない単語が出てきたとき、別の端末で調べる必要がある。
△　論文を読んでいて調べ物をしたい場合、ネットで快適に調べることができない。

　
数ヶ月ほど使って「ナンカチガウ」感が出てきて、だんだんと使わなくなっていきました。「これではあんまり論文が読めないなぁ」と困り、iPad Pro はどうだろうと思い立ちました。iPad Pro には大きいサイズ（12.9インチ）と小さいサイズ（当時9.7インチだった）があります。

　
大きい iPad Pro を仕事に使っている知人がいたので見せてもらって、ちょっと大きいかなと思ったので、小さい iPad Pro を使うことに決め、アップルの整備済製品で購入しました。お得です。

iPad整備済製品 - Apple（日本）

小さい iPad Pro のいいところ

　
○　字がきれい。
○　知らない英単語が出てきたとき、付属の英語辞書ですぐわかる。
○　wi-fi につなげて、快適に調べ物ができる。
○　アップルペンシルはなかなかいい感じ。

　
iPad の辞書機能は特に便利です。私は英語に不自由しているので、よく知らない単語が出てきます。これを他の電子辞書を使わずに、iPad 内蔵の辞書で調べることができるのでストレスがとても減り、QOLアップ。　

www.youtube.com

　
PDFを読むとき、PDF Expert を使うのですが、上のように知らない単語を選択して「定義」をタップすれば、内蔵辞書ですぐに意味を確認できます。超便利！

　
アップルペンシルの書き心地は動画では伝わりにくいですが、こんな感じ。

youtu.be

小さい iPad Pro のイマイチなところ

　
△　使っていると、だんだん小さく感じてくる。
△　ケースに入れて使っていると、思ったより重く感じる。
△　アップルペンシルはいいけど、ちょっと滑りすぎる。

　
また半年ほど使って、また「ナンカチガウ」感が出てきたので、大きい iPad Pro を本格検討。新型の大きい iPad が出て、かなり迷いました。清水の舞台から飛び降りる気持ちで移行を決意。在庫があるなら全国どのアップルストアにも行く（出張ついでに）覚悟でした。

　
結果、銀座にアップルストアにあったので、すぐにゲット。そして現在に至ります。まだそんなに経っていないですけど、とても満足して使っています。

大きい iPad Pro のいいところ

　
○　小さい iPad Pro のいいところは全部ある。
○　さらに画面が大きく、見やすい。

大きい iPad Pro のイマイチなところ

　
△　デカい。

大きい iPad Pro をより快適に使う方法

　
以下、より楽しく、快適に iPad Pro 12.9 インチを使って論文を読むために試した方法を書きます。

アップルペンシルのキャップカバー（必須）

　
まずアップルペンシルは必須です。しかし、アップルペンシルの最大の弱点はキャップ。小さい iPad Pro を使っているときにアップルペンシルはすでに活用していましたが、すぐにキャップをなくしました。

　
キャップがなくても普通に使えますが、なんかカッコ悪い。仕方がないので、大きい iPad Pro 導入時にキャップだけ再度購入しました。アップルストアで修理扱いです。修理の予約が必要だし、修理代は1000円くらい。

　
また紛失したら悲しいので、キャップカバーを導入。修理代より安いので、導入したほうがいいです。どのキャップカバーでも構いません。

LENSEN Apple Pencil キャップカバーアダプタ連絡用カバー Pencil Cap Cover Lighting Adapter Cover シリコン１セット (白 LSAPT-008)

出版社/メーカー: LENSEN
メディア: エレクトロニクス
この商品を含むブログを見る

www.youtube.com

　
ちょっと不格好なのですが、紛失のリスクを考えて、ここはガマンですかね。

ペーパーライク保護フィルム（必須）

　
ソニーデジタルペーパーと比べ、iPad Pro はツルッツルです。そのためアップルペンシルが滑りすぎます。多少の引っかかりが必要なため、ディスプレイが紙のような感じになる「ペーパーライク保護フィルム」を導入しました。

エレコム 2017年新型 iPad Pro 12.9 / 2015年発売 iPad pro 12.9 液晶保護フィルムペーパーライク反射防止 TB-A17LFLAPL

出版社/メーカー: エレコム
発売日: 2017/06/07
メディア: Personal Computers
この商品を含むブログを見る

　
最近は安いのが出ているので、好きなのを選べばいいかと思います。ただ、気をつけてほしいことがあります。それは、せっかくの美麗なディスプレイが台無しになることです。

　
画面はちょっとざらついた見た目になりますし、指紋はむっちゃ付着します。きれいな画像・映像を見たい人にはガマンできないと思います。

　
しかし、この iPad Pro は論文を読むためにあるのです！　割り切ってエンタメは諦めましょう…。次第に目は慣れていきます。なんか違和感のある手触りも慣れます。ざらつきはありますが、画面はヌルヌル動くので、操作は問題ありません。

www.youtube.com

PDFリーダー（必須。お好きなのを）

　
先ほど書いたように、PDF Expert を使っています。他のPDFリーダーも試して、いちばんしっくりきたためです。かなり感覚的でちゃんとした理由はありません。ドロップボックスに読むべき論文のフォルダを作り、そのフォルダと PDF Expert を同期させます。

　
f:id:hikaru1122:20170920015459p:plain

格安SIM（お好きにどうぞ）

　
wi-fi でつながっていないときに（電車の中とか）で論文を読んでいて、ネットで調べ物をしたいときに便利です。ただ、つい余計なものを見てしまうときもあるので、導入するかどうかは好き好きで。私は楽天ポイント（ホテル手配は楽天トラベルだからけっこう貯まる）で支払える楽天SIMにしています。

www.youtube.com

　
アンテナが2本しか立っていなかったのでスピードが少し遅いですね。大きな問題はないです。

スリーブケース（お好きにどうぞ）

　
ケースは重いので、スリーブケースにしました。デスクにどかんと置いて使うか、手に持って使うだけなので、別に立てて使うことはないからです。アマゾンとかで好きなモノを選んでください。

　
私は純正品を使っていますが、キャップ付きのアップルペンシルだと収まりが悪く、ちょっとカッコわるいです。

f:id:hikaru1122:20170920012251j:plain

　
というわけで、大きい iPad Pro であなたも快適な論文読みライフを送りませんか？　コスパには目をつむってになりますけど…。

2017-08-24

第1回 Kobe.R 別館を開催します！

Kobe.R

f:id:hikaru1122:20170824004200j:plain

来る2017年9月24日に Kobe.R の弟分的な存在の勉強会「Kobe.R annex（別館）」第1回を開催します。小さく、緩く（いい意味で）、楽しくやっていきたいと思います。２部制にして、第1部は『ベイズ統計モデリング: R,JAGS, Stanによるチュートリアル原著第2版』（犬4匹本）を輪読します。第2部では、初心者セッションおよびＬＴ大会です。第1部から参加、もしくは第2部のみ参加をお選びいただけます。

　
参加者数は最大11名と小さいため、ＬＴ枠などが埋まらない場合は、参加者同士で時間の許す限りＲ関連の話題やＲを使った研究ネタについて（アイデアレベルでも）と話し合おうと思っています。初心者セッションは “R for Datascience” または “Quantitative Social Science: An Introduction” を題材にしてもいいかな～とも思っています。

　
日曜の大阪・難波での開催です。すでに滋賀県から参加される方もいらっしゃいます。関西在住の方（もちろんその他の地域の方も）は奮ってご参加ください！

　
kobe-r-annex.connpass.com

2017-07-30

Proximity Matching による因果効果の推定？

ノンプログラマーのためのR入門参加報告統計分析統計的因果推論

2017年7月30日のKazutan.Rで「Proximity Matching による因果効果の推定」というタイトルで LT をしてきました。スライドは↓ですが，R のスクリプトは少なめなので，その補足説明をここに残しておこうと思います。

https://hikarugoto.github.io/kazutanR20170730/

LT の動機

　
今回の発表のモチベーションは，たまたま次の論考（何かの予稿集かな？）を見かけ，傾向スコアに代わる何かがあるのだったら，調べてみようかなと思った次第です。ランダムフォレストと傾向スコアについては，＠TJO さんや＠nakamichi さんなどがすでにブログで書いていらっしゃいますが，この論考では Proximity というのを使ってマッチングを行うものなので，ちょっと方法が違うようです。

Zhao, P., Su, X., Ge, T., & Fan, J. (2016). Propensity score and proximity matching using random forest. Contemporary Clinical Trials, 47, 85-92.

　
ランダムフォレストから Proximity を出して，マッチングを行う例は，すでに “Causal inference with observational data in R: A walk-through from some recent research” という YouTube 動画と資料があったので今回の LT ではそれをなぞってみました。

使用したデータ

　
統計的因果推論といえば，岩波データサイエンス vol.3 ！そちらのサポートページで公開されているデータを使いました。

github.com

岩波データサイエンス Vol.3

作者: 岩波データサイエンス刊行委員会
出版社/メーカー: 岩波書店
発売日: 2016/06/10
メディア: 単行本（ソフトカバー）
この商品を含むブログ (4件) を見る

　
10000行で35変数ありますが，「CM視聴の有無でゲーム利用回数はどの程度異なるか」をアウトカムにして，23変数を使うことにしました。イメージ的にはこんな感じ（あんまりいい図ではないけど）。

f:id:hikaru1122:20170730181643p:plain 　
上の図は DiagrammeR で書きました。

grViz("
      digraph boxes_and_circles {
      # a 'graph' statement
      graph [overlap = true, fontsize = 10]
      
      # several 'node' statements
      node [shape = box,
      fontname = Helvetica]
      原因変数（CM視聴の有無）; ゲーム利用回数

      node [shape = circle,
      width = 1] // sets as circles
      共変量
      
      # several 'edge' statements
      共変量->{原因変数（CM視聴の有無） ゲーム利用回数}
      原因変数（CM視聴の有無）->{ゲーム利用回数}
      }
      ")

使用したパッケージ

　
使ったパッケージは次のとおりです。

library(ggplot2)
library(randomForest)
library(dplyr)
library(magrittr)
library(broom)
library(knitr)
library(DiagrammeR)
library(ggpubr)

分析の手順

　
まず，使用する変数を元のデータから取り出します。select で -（マイナス）を使い，不要な列を削除しました。

CI_data %<>% select(-gamedummy, -area_keihan, -job_dummy8, -fam_str_dummy5, -T:-M3, -gamesecond)
CI_data %>% glimpse()

Observations: 10,000
Variables: 23
$ cm_dummy        <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
$ area_kanto      <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, ...
$ area_tokai      <int> 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
$ area_keihanshin <int> 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, ...
$ age             <dbl> 44.5, 34.5, 24.5, 44.5, 34.5, 24.5, 19.0, 54.5, 44.5, 34.5, 34.5, 34.5, ...

＜以下，省略＞

　
%<>% はあまり見かけませんが，左からデータフレームを右に流していって，その後，同じ名前のオプジェクトに置き換える破壊的な操作をするパイプです。今回は軽い分析なので使ってみました。 2行目はチラ見関数です。str() よりもモダンな感じ？

ランダムフォレストの実行

　
proximity を得るために，ランダムフォレストを実行します。

cm_dummy = CI_data$cm_dummy==1 # 論理ベクトル
rf = randomForest(cm_dummy ~ ., data = CI_data, proximity = TRUE) # 1分くらいかかる

　
1行目は CM視聴した有無を論理値ベクトルで入手し，それを randamForest で使います。なぜ，論理値ベクトルにする必要があるのか，私には不明です。proximity = T のオプションを忘れないのが重要です。

Proximity でマッチングさせたデータを作る

　
参考資料では，以下のようにして，マッチングを行っています。

prox = rf$proximity - diag(nrow(rf$proximity))  # 何やってるの？
prox.true = prox[cm_dummy, !cm_dummy] # 何やってるの？
untre.samples = apply(prox.true, 1, which.max) # 行ごとに proximity が最大の列を抜き出し？

tre = CI_data[cm_dummy,] # 処置群のみ取り出し
untre = CI_data[!cm_dummy,] # 対照群のデータのみ取り出し

untre.samples2 = unique(untre.samples) # 重複している行を削除。
untre2 = untre[untre.samples2,] # 処置群と似ている対照群を抜き出している？
whole.sample.u = rbind(tre, untre2) # それぞれのデータを結合。

　
もっとも重要なところなのでしょうが，勉強不足でよくわからん･･･。目をつぶって先に進みます。

共変量のバランスは取れたのかチェックしてみる。

　
せっかくなので，CM見た群と見なかった群の共変量の分布が近づいたかどうか，TVwatch_day と pmoney について見てみます。マッチング前後で見比べるため，ggplot で作図したものを並べます。ggplot を並べる方法はいくつかありますが，今回は ggpubr の ggarrange() を使います。

a1 = CI_data %>% 
ggplot(aes(x = TVwatch_day)) + geom_density(aes(colour = factor(cm_dummy) ,fill = factor(cm_dummy)), alpha = 0.5)
a2 = whole.sample.u %>% 
ggplot(aes(x = TVwatch_day)) + geom_density(aes(colour = factor(cm_dummy) ,fill = factor(cm_dummy)), alpha = 0.5)

b1 = CI_data %>% 
ggplot(aes(x = pmoney)) + geom_density(aes(colour = factor(cm_dummy) ,fill = factor(cm_dummy)), alpha = 0.5)
b2 = whole.sample.u %>% 
ggplot(aes(x = pmoney)) + geom_density(aes(colour = factor(cm_dummy) ,fill = factor(cm_dummy)), alpha = 0.5)

　
TVwatch_day のマッチング前の密度分布（a1），マッチング後（a2）として，pmoney も同様にマッチング前後（b1 と b2）を ggplot() で作ります。並べるには次のようにggpub パッケージの ggarrange()を使います。

ggarrange(a1, a2, b1, b2, ncol = 2, nrow = 2, 
labels = c("Before Matching","After Matching","Before Matching","After Matching"), 
common.legend = TRUE)

　
図を2行2列で並べて，それぞれマッチング前とマッチング後のラベルをつけ，共通の凡例を表示させます。

　
f:id:hikaru1122:20170730162121p:plain

　
確認しても，これで共変量が調整されているかなんとも言えないです･･･。なお，この並べ方は次のページを参考にしました。

http://www.sthda.com/english/wiki/ggplot2-easy-way-to-mix-multiple-graphs-on-the-same-page#.V5F9h53_d2E.twitter

因果効果の推定結果を確認

　
参考資料には「あとは回帰分析すればいいぜ」と書いてあったので，やってみます。なんか怪しい。まぁ，RCT に近づいたからやっていいという考え方だとは思うのですが･･･。

lm(gamecount ~ ., data = whole.sample.u) %>% 
tidy() %>% 
head(3) %>%
kable()

　
lm() で得た結果を2行目でデータフレームにして，3行目ではじめの数行だけ表示させ，4行目で下の表を出すためのマークダウンを得ました。

term	estimate	std.error	statistic	p.value
(Intercept)	46.467720	7.462751	6.226621	0.0000000
cm_dummy	5.080259	1.499469	3.388039	0.0007088
area_kanto	-3.428637	2.041777	-1.679242	0.0931606

　
上の表から，CMを見た場合の平均処置効果（ATE）は 5.08 であることがわかりました。岩波DS vol.3 でIPWによってATEを出したときの結果が 5.32 ですから（96ページ），近いと言えば近い*1。

	ATE
岩波DS vol.3	5.32
RF で IPW @nakamichi*2	5.32
Proximity でマッチング	5.08

分析してみた感想

　
やっていることをよく理解できず，今回の結果を信じるには至っていません。ちなみに傾向スコアマッチングは使わないほうがいいという論考が出ています*3。今回の方法は傾向スコアの弱点をクリアしているものの*4，IPW法など他のやり方があるので，わざわざマッチングをしなくてもいいように思いました。

最後に

　
やはり，Kazutan.R に参加してよかったです。主催者 kazutan 先生やハンズオンの MrUnadadon氏，LTされた方々に感謝。私にはデータ分析をする人が身近にいないので，情報の入手先は twitter やこういう勉強会がすべてなので（書籍の情報も twitter で知ることが多い），助かります。

*1:ちなみに，スライドで2.32と書いたのは ATT でまちがいでした。スミマセン

*2:https://rpubs.com/nakamichi_/study-iwanami-ds3

*3:King and Nielsen （2016）"Why Propensity Scores Should Not Be Used for Matching"

*4:Zhao らによれば，傾向スコアの限界は「model misspecifications, categorical variables with more than two levels, difficulties in handling missing data, and nonlinear relationships」とのことです。

2017-07-07

テキストマイニングの環境整備つまずきメモ

テキストマイニングノンプログラマーのためのR入門書籍紹介

テキストマイニングをしようと決意して，その準備をしたところ，いくつかつまづいたところがあったので，そのメモを残しておきます。解決策はネット上に散らばっていて，同じようなトラブルにまたあったときのためです。なお，Mac です。

1．環境構築の基礎知識をゲット

　
まず『Rによるやさしいテキストマイニング』（やさテキ）を読み，次に『Rによるテキストマイニング入門（第2版）』を読みました。2冊でざっとテキストマイニングの基礎知識と分析の流れをつかみます。どちらも読みやすいですが，まず『やさテキ』から『Rによる〜』の順がよいと思います。

Rによるやさしいテキストマイニング

作者: 小林雄一郎
出版社/メーカー: オーム社
発売日: 2017/02/17
メディア: 単行本
この商品を含むブログを見る

Rによるテキストマイニング入門

作者: 石田基広
出版社/メーカー: 森北出版
発売日: 2017/06/28
メディア: 単行本
この商品を含むブログを見る

2．MeCab のインストール時のつまずき

　
『Rによる〜』に書かれている方法で MeCab のインストールを行ったところ，うまく行きませんでした。ターミナル上で「mecab すもももももももものうち」と打っても結果が

??,????,,,,,?,???????,???????

　
と文字化けするのです。これは MeCab をもう1度次の方法でインストールして解決できました。

$ ./configure --enable-utf8-only --with-charset=utf8 --with-mecab-config=/usr/local/bin/mecab-config
$ make
$ make install

　
詳細はこちらにあります。

kumagonjp2.blog.fc2.com

3．RMeCab インストール時のつまずき

　
私の環境では『Rによる〜』に記載されている方法でインストールして RMecab を使うと RStudio ごと落ちるという現象が起きました（Rコンソールでもダメ）。tweet したところ，石田先生に助けていただきました。

辞書もソースからインストールされているということですよね。では、Rで　install.packages("RMeCab", repos = "https://t.co/ItGfYL8TAs", type = "source")　を実行してみてください
— RMeCab (@rmecab) 2017年6月30日

　
これで無事解決。ありがとうございます。

4．mecab-ipadic-NEologd を使用するときのつまずき

　
テキストマイニングの対象はちょっと特別な言葉が入っているので，標準の辞書ではない NEologd を使うことにしました。この辞書のインストールは次のページに従って行いました。

qiita.com

　
しかし，RMeCab で NEologd を使うときは，ユーザ辞書が必要とのこと。次のページの「R」セクションの記述に従って，ユーザ辞書を作り直しました。

github.com

　
作った辞書（〜〜.csv.dic という名前のファイル）は，適当なところにコピーして使います。こんな感じ。

> RMeCabC("オールセラミック", 
　　　　　　dic = "/Users/jibun/mecab-user-dict-seed.20170630.csv.dic")
[[1]]
              名詞 
"オールセラミック"

　
標準の辞書だと，

> RMeCabC("オールセラミック")
[[1]]
  接頭詞 
"オール" 

[[2]]
        名詞 
"セラミック"

　
ちゃんと，オールセラミックを1つの言葉として扱ってくれてる。よかった。

　
まだまだテキストマイニングは準備すら難しい印象を持ちました（私の環境では）。でも『やさテキ』と『Rによるテキストマイニング第2版』とネット上の先達のおかげでなんとかなりました。KH Coder に流れたくなったけど，せっかくなんとか準備が整ったので R でがんばりたいと思います。

Rによるやさしいテキストマイニング

作者: 小林雄一郎
出版社/メーカー: オーム社
発売日: 2017/02/17
メディア: 単行本
この商品を含むブログを見る

Rによるテキストマイニング入門

作者: 石田基広
出版社/メーカー: 森北出版
発売日: 2017/06/28
メディア: 単行本
この商品を含むブログを見る

2017-05-20

R による整列ランク変換を使った分散分析

ノンプログラマーのためのR入門参加報告読んだ論文統計分析

整列ランク変換（Aligned Rank Transform：ART）という言葉を聞いたことありますか？　私はなかったです。ART（えーあーるーてぃー。「あーと」ではない）について、Tokyo.R 第61回で LT をしてきました。補足も兼ねて ART についてまとめます。

きっかけ

　
先日、自分の専門分野の先行研究を調べていて論文（Heidenreich, et al. 2015）を読んでいたら次のような記述に出会いました。

“As our data were not normally distributed (D(243)=2.193, p<0.01), we employed the aligned rank transform (ART) procedure for hypothesis testing. For non-normal distributed data, the ART procedure is more robust and powerful than the traditional analysis of variance (ANOVA)”

　
超訳すると「データが正規分布に従わなかったので、aligned rank transform (ART) を使って検定やった。非正規分布のデータなら、ARTを使ったほうが普通の分散分析よりロバストでパワフルだ」。

　
ART について、ネット上で検索してもほとんど日本語の情報がありません。かろうじて「整列ランク変換」という訳語が見つかるくらいです。そこで少し調べてみました。

ART法はタイプIエラーをコントロールできるエクセレントな方法だ。研究者がやりがちな古典的分散分析と比べてパワーがある（Mansouri, et al. 2004）

ART はどんな分布のデータにも使えるわけじゃないけど、ロバストだし一般的なノンパラメトリック検定より望ましい性質を持っている（Higgins, et al. 1990）

ART 検定はパワフルでロバストで使いやすい（Mansouri 1998）

交互作用を調べるとき、いくつかの仮定が成り立たない場合、パラメトリックな分析は必ずしもよいわけではない。代替の策として、調整ランク変換検定（adjusted rank transform test）がある。それはパラメトリックな手法とノンパラメトリックな手法の中間みたいなものだ（Leys & Schumann 2010）

　
つまり、分散分析したい、特に交互作用に注目したい時にデータが正規分布に従わない場合、従来の分散分析は（厳密には）使えないので、整列ランク変換を施してから分散分析すればいいよ、ということだと思います。

　
R パッケージがないか探したところ、ARTool というものがちゃんとありました。CRAN に登録されています。使い方もカンタン。以下、ART を使った分散分析の方法です。例題として、性別と世代の不倫の許容度に対する違い・交互作用を見たいとします。応答変数は「結婚していても場合によっては不倫を許せる？（１：許せない～４：許せる）」です。4件法でスミマセン。

０．データの分布の確認

ggplot(rensyu, aes(furin)) + 
geom_histogram(fill = "red", alpha = 0.7, binwidth = 0.5)

f:id:hikaru1122:20170520195014p:plain

　
あきらかに正規分布には従わないっぽい。もちろん、正規性の検定はすべきですが、ここでは飛ばします。

１．ARTool を使えるようにする。

install.packages("ARTool") # 1回のみでOK。
library(ARTool)

２．データを整列ランク変換する。

henkan = art(furin ~ nendai*sex, data = rensyu)
summary(henkan)

　
art() で変換しますが、式はよくあるタイプです。式は「art(furin ~ nendai + sex + nendai:sex, data = rensyu)」を短縮した書き方にしています。どちらでもOK。lm() のときと同じです。

　
変換後、summary(henkan) の出力結果に注意します。すべての数値が０になっていることが望ましいです。今回は残念ながらそうならなかったので、警告が出ています。しかし、練習なので先に進めます。

## Warning in summary.art(henkan): F values of ANOVAs on aligned responses not
## of interest are not all ~0. ART may not be appropriate.

## Aligned Rank Transform of Factorial Model
## 
## Call:
## art(formula = furin ~ nendai * sex, data = rensyu)
## 
## Column sums of aligned responses (should all be ~0):
##     nendai        sex nendai:sex 
##          0          0          0 
## 
## F values of ANOVAs on aligned responses not of interest (should all be ~0):
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0000000 0.0000000 0.0000000 0.0144000 0.0007531 0.0854000

３．分散分析を実行する。

anova(henkan)

## Analysis of Variance of Aligned Rank Transformed Data
## 
## Table Type: Anova Table (Type III tests) 
## Model: No Repeated Measures (lm)
## Response: art(furin)
## 
##              Df Df.res  F value     Pr(>F)       Sum Sq Sum Sq.res
## 1 nendai      4   1261   3.7856  0.0045704  **  1935456  161178877
## 2 sex         1   1261 100.7097 < 2.22e-16 *** 10944212  137033949
## 3 nendai:sex  4   1261   4.6317  0.0010275  **  2369093  161247964
## ---
## Signif. codes:   0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

　
「anova(変換したデータ)」でOK。実際には、anova.art() が動いて、デフォルトでタイプ３の分散分析をやってくれます。世代と性別の主効果、交互作用は有意になりました（なお、効果量 $\eta^{2}$ は世代がほぼ0、性別が0.02、交互作用が0.01でほとんどインパクトはありません）。いちおう、交互作用のグラフをプロットします。

library(ggplot2)
ggplot(rensyu) +
aes(x = nendai, color = sex, group = sex, y = furin) +
stat_summary(fun.y = mean, geom = "point", size = 2) +
stat_summary(fun.y = mean, geom = "line", size = 1) +
xlab("世代") + ylab("不倫の許容度") +
scale_colour_discrete(name ="性別")

f:id:hikaru1122:20170520224157p:plain

伝統的な分散分析と結果を比べてみる。

aaa = lm(furin ~ nendai*sex, data = rensyu)
library(car)
Anova(aaa, type="III")

## Anova Table (Type III tests)
## 
## Response: furin
##             Sum Sq   Df  F value  Pr(>F)    
## (Intercept) 224.03    1 535.4585 < 2e-16 ***
## nendai        2.39    4   1.4269 0.22278    
## sex           2.22    1   5.3061 0.02141 *  
## nendai:sex    3.02    4   1.8051 0.12541    
## Residuals   527.59 1261                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

　
同じくタイプ３でやってみましたが、こちらでは世代の主効果と交互作用が有意になりませんでした。いい感じに対比できましたが、データが怪しいので、結果はあくまで参考用として受け取ってください。

まとめ

整列ランク変換を利用した分散分析では交互作用が有意になった。
ただし、今回のデータはちゃんと整列ランク変換できてはいないのでなんとも言えない。
今後、データが正規分布に従わない場合、使ったほうがよい手法かもしれない。

さいごに

　
多重比較もできます。繰り返しのデータも分析可能のようです。詳しくはパッケージ作者の Github ページを参考にしてください。今後、必要あれば、使っていきたい分析法です。

参考文献

ARTool パッケージ
https://github.com/mjskay/ARTool
Feys, J. (2016). Nonparametric Tests for the Interaction in Two-way Factorial Designs Using R. The R Journal, 8(1), 367-378.
Heidenreich, S., Wittkowski, K., Handrich, M., & Falk, T. (2015). The dark side of customer co-creation: exploring the consequences of failed co-created services. Journal of the Academy of Marketing Science, 43(3), 279-296.
Higgins, J. J., Blair, R. C., & Tashtoush, S. (1990). The aligned rank transform procedure.
Leys, C., & Schumann, S. (2010). A nonparametric method to analyze interactions: The adjusted rank transform test. Journal of Experimental Social Psychology, 46(4), 684-688.
Mansouri, H. (1998). Multifactor analysis of variance based on the aligned rank transform technique. Computational statistics & data analysis, 29(2), 177-189.
Mansouri, H., Paige, R. L., & Surles, J. G. (2004). Aligned rank transform techniques for analysis of variance and multiple comparisons. Communications in Statistics-Theory and Methods, 33(9), 2217-2232.
Wobbrock, J.O., Findlater, L., Gergle, D. and Higgins, J.J. (2011). The Aligned Rank Transform for nonparametric factorial analyses using only ANOVA procedures. Proceedings of the ACM Conference on Human Factors in Computing Systems (CHI ‘11). Vancouver, British Columbia (May 7-12, 2011). New York: ACM Press, 143-146.

目を通してないけど、重要そうな文献

Salter, K. C., & Fawcett, R. F. (1993). The ART test of interaction: a robust and powerful rank test of interaction in factorial models. Communications in Statistics-Simulation and Computation, 22(1), 137-153.
Sawilowsky, S. S. (1990). Nonparametric tests of interaction in experimental design. Review of Educational Research, 60(1), 91-126.

2017-04-09

RMarkdown で参考文献を作るときにちょっとつまずいた話

ノンプログラマーのためのR入門研究に役立つツール

f:id:hikaru1122:20170409035108j:plain

Rを使うドキュメントはできるだけ RStudio内で完結するように（コピペが面倒なのと，スライドも配布するPDF資料もRStudioでできるので）環境を整えているところですが，参考文献も自動で挿入できるようにしたいと思いました。やってみて，つまずいたところをメモします。

　
RMarkdown チートシートを見ると，2ページ目に「Citations and Bibliographies」というセクションがあります。YAML ヘッダーに簡単な追記をすればできるようなので，やって見たところ，次のように怒られました。

pandoc: Filter pandoc-citeproc not found

　
エラーメッセージをそのまま検索すると，こちらの答えがありました。

　
github.com

　
どうやら，

brew install pandoc pandoc-citeproc

　
でいいそうです（Mac で homebrew が使える場合）。確かに，pandoc-citeproc をインストールした後は，YAMLヘッダーに参考文献のリスト（.bib）を用意すれば，うまく参考文献が表示されるようになります。

　
今回，日本語文献だったので，順番がうまくいきませんでしたが，論文ではないので，まぁOKかなと思います。記述したYAMLヘッダーはこちら。

---
title: "Rの操作の復習と簡単な分散分析"
author: Hikaru GOTO
date: 2017年4月15日（土）
output:
  rmdshower::shower_presentation:
    self_contained: false
    theme: ribbon
    ratio: 4x3
bibliography: mylib.bib
csl: journal-of-business-research.csl
---