JuliaからRCallを使ってANOVA君を呼び出す

ANOVA君は、井関龍太氏が開発している、Rで動作する分散分析プログラムです。様々な実験計画に対応できるようになっており、分散分析に加えて多重比較や単純主効果の検定もやってくれます。自分ではlm()とかaov()なども使ってはいますが、ANOVA君は一回の関数呼び出しですべてを面倒見てくれるのでとても重宝します。

riseki.php.xdomain.jp

数年前に僕がPythonに移行できなかった理由の一つに「まともな分散分析ツールがない」ということがありました。1要因や簡単な2要因のものくらいなら見つかるのですが、3要因以上だったり2要因でも混合計画のものは見つからず、かといって自分で作るのも車輪の再発明をしているかんじでイヤだったのです。Juliaに移行するにあたっても同じ理由でRを捨てられずにいます。

しかし、JuliaはPythonともRとも親和性が高く、それらで書かれたプログラムも簡単に呼び出すことができます。今日は、JuliaでRを呼び出してANOVA君を使う例を紹介します。

まずは必要なパッケージを読み込みます。

julia> using DataFrames
julia> using CSV
julia> using RCall

このRCallが秀逸で、julia>プロンプト時に$を押すとR>プロンプトに変わり、Rの命令がそのまま実行できるようになります(]を押すとPkgモードに入れるのと似たような感じです)。Pkg同様にバックスペースを押すとjulia >プロンプトに戻れます。

julia> 【$を入力】

R> rep(c(1,2,3), 4)
 [1] 1 2 3 1 2 3 1 2 3 1 2 3

R> 【バックスペースを入力】

julia> 

直接Rの命令を実行するにはR命令を使います。Rのあとにダブルクォートで挟んだRコマンドを入れると実行されます。ただし、Rコマンドの中にダブルクォートを入れたい場合はバックスラッシュでエスケープする必要があります。

julia> R"rep(c(1,2,3), 4)"
RObject{RealSxp}
 [1] 1 2 3 1 2 3 1 2 3 1 2 3

julia> R"rep(c(\"A\",\"B\"), 4)"
RObject{StrSxp}
[1] "A" "B" "A" "B" "A" "B" "A" "B"


では、ここからANOVA君を使った分散分析。現時点でANOVA君のページからダウンロードできる最新版がanovakun_484.txtなので、それをローカルに保存して読み込みます。(バックスラッシュでのエスケープを使っています)

julia> R"source(\"anovakun_484.txt\", encoding=\"utf8\")";

より高度な入力方式」ページにある「ロング形式によるデータ入力方式」以下にあるデータを使ってみます。同ページではきれいにフォーマットされていますが、ここでは以下のようなCSVファイル(study_data.csv)として保存してあるとします。

使用するデータは、実験参加者ID、要因WMとStudy、従属変数としてtestがあります。WMは被験者間要因、Studyは被験者内要因の、混合計画になっています。

ID,WM,Study,test
s1,high,lecture,4
s2,high,lecture,6
s3,high,lecture,8
s4,high,lecture,7
s5,low,lecture,4
s6,low,lecture,2
(以下略)

以下のようにCSVファイルを読み込みます。expanduser()はホームディレクトリのパス名を補完してくれるものですので、必要ないかもしれません。

julia> filename = expanduser("~/hoge/study_data.csv");
julia> dat = CSV.read(filename)

datの中身を確認すると、次のようになっています。

julia> dat
30×4 DataFrame
│ Row │ ID     │ WM     │ Study   │ test  │
│     │ String │ String │ String  │ Int64 │
├─────┼────────┼────────┼─────────┼───────┤
│ 1   │ s1     │ high   │ lecture │ 4     │
│ 2   │ s2     │ high   │ lecture │ 6     │
│ 3   │ s3     │ high   │ lecture │ 8     │
⋮
│ 27  │ s7     │ low    │ peer    │ 8     │
│ 28  │ s8     │ low    │ peer    │ 7     │
│ 29  │ s9     │ low    │ peer    │ 9     │
│ 30  │ s10    │ low    │ peer    │ 7     │

ANOVA君の実行は以下のようにします。Julia側で読み込んだdat変数をR側に渡すのに$(dat)としています(文字列への変数内容の埋め込みと同様の文法でできるのがおもしろいですね)。

julia> R"anovakun($(dat), \"AsB\", long=T)"

結果は省略!