T. Yoshiizumi - rrx_guide Diff
- Added parts are displayed like this.
- Deleted parts are displayed
like this.
rubyから統計解析ソフトRを使う 〜 rrxwin.rb
最終更新日： 2013/07/03
以下に掲げるドキュメントは、
[[rrx102.zip|http://cup.sakura.ne.jp/rrx102.zip]]
に同梱されている rrxwin.txt と同じ内容です。
--------
rubyから統計解析ソフトRを利用する方法について記します。
rubyとRの連携を取りやすくするためのライブラリ rrxwin.rb を作ったので、それを用います。
MS-Windows用のrrxwin.rbの他に、linuxなどで使えるrrx.rbもあります。それについては同梱の rrx.txt を参照して下さい。
以降で掲げるサンプルスクリプトを私が実行した環境は次のとおり。
* MS-Windows xp | vista
* ruby ver 1.8.7 | 1.9.3
* R ver 2.15.2
--------
{{toc_here}}
----
!＜はじめに＞
統計解析ソフトRは、とても高機能で統計解析のためのいろいろな仕組みを装備しています。また、Rにはプログラミング言語として様々な仕組みがあり、それに熟達している人にとっては、Rだけで処理できることが大半だろうと思います。
ただ、私は ruby を愛用しており、それと連携させてRを使いたいと考えました。rubyを経由することで、Rの解析結果を容易にExcelに貼り付けることができます。また、sqliteなど様々なデータベースに蓄積されたデータをRで解析するのも容易になります。
rubyからRを利用する強力なライブラリとして rsruby があります。これについても機会があれば書いてみたいと思いますが、もう少し手軽に使えるものとして rrxwin.rb を作ったので、その利用方法を記します。
１つのcsvファイルから複数の表(行列)を切り取るためのメソッドなども用意し、R活用の時の周辺的な便利さを追求したつもりです。
使用環境は MS-Windows です。
Rとrubyの両方が予めインストールされている必要があります。
文字コードは Windows-31J(cp932) を前提にしています。
Rの使い方には深入りしません。Rについては様々な解説がWebにあるのでそれらを参照して下さい。
サンプルではカイ２乗検定(クロス集計表の検定)、ｔ検定(平均値の差に関する検定)、相関係数の３つを取り上げます。また、相関係数に関連して散布図の描画も取り上げます。
--------
!1. rrxwin.rb のインストール
!!(1) インストールせずに使う方法
rrxwin.rb がカレントディレクトリにある場合は、rubyスクリプトの先頭の方に「require "./rrxwin"」を書いておけば、以下のインストール作業をしなくても使えます。
どのディレクトリからでも使えるようにしたい時は、下に掲げるインストール作業を行って下さい。
!!(2) インストール方法
rrxwin.rb, setup.rb の二つのファイルがカレントディレクトリにあることを確認して、次のようにrubyを実行します。
ruby.exe setup.rb [enter]
こうすると、rubyのライブラリパスが通ったディレクトリに rrxwin.rb がコピーされます。
このインストールの際に、統計解析ソフトRの実行ファイルである Rcmd.exe がどのディレクトリにあるかを検出すると、そのフルパス名を標準出力に出力します。また、検出できない時はその旨のメッセージを出力します。
検出できなかった場合、Rcmd.exeがパスの通ったディレクトリにないのであれば、そのフルパス名を利用者がセットしてやる必要があります。それは次のようなrubyスクリプトを実行して行います。
−−−− set_path.rb ここから
#! ruby -Ks
require "rrxwin"
pathname = "C:/usr/R/R-2.15.2/bin/i386/Rcmd.exe"
Rrx.path = pathname
−−−− set_path.rb ここまで
"C:/usr/R/……/Rcmd.exe" の部分は、ご自分のものに置き換えた上で実行して下さい。これを実行すると、そのフルパス名が登録されます(rrxwin.rbと同じディレクトリにあるrrxwin.cfgに書き込まれます。)。
Rをインストールする際、デフォルトのままで行うと次のようなフルパスになると思います。
C:/Program Files/R/R-2.15.2/bin/i386/Rcmd.exe
「2.15.2」は Rのバージョン番号です。どのRをインストールしたかによって違ってきます。
ユーザー権限等の関係で、"C:/Program Files" の下にディレクトリを作れない状況であれば、Documentディレクトリの下にあるかもしれません。例えば次のようなものです。
C:/Users/Taro/Documents/R/R-2.15.2/bin/i386/Rcmd.exe
!!(3) アンインストールの方法
上のようにしてインストールしたrrxwin.rbを削除したい時は unset.rb を実行します。
ruby.exe unset.rb [enter]
とすれば、rrxwin.rbとその関係ファイル(rrxwin.cfg)が削除されます。
unset.rbを実行する時は、カレントディレクトリに rrxwin.rb がなくてもかまいません。
!!［参考１］ rrxwin.rbが Rcmd.exe のパスを検索するやり方
rrxwin.rbは、アクセス可能なドライブ(ハードディスク、USBメモリーなど)を順番に検索しますが、各ドライブについて次のディレクトリがあるかどうか調べます。該当のものがみつかれば、そこで検索を中断します。
* /R (ルートディレクトリ直下の R というディレクトリ)
* /Program Files/R
* /Program Files (x86)/R
* /Users/Taro/Documents/R (Taro はユーザー名の一例)
上の４種類のディレクトリがどのドライブでもみつからないと、rrxwin.rbは、Rcmd.exeのパスを検出できません。その時は、前述した set_path.rb でパスを設定して下さい。
Rというディレクトリの下にバージョンの異なるRが複数ある場合は、最も新しいバージョンを自動的に選びます。
!!［参考２］ Rのインストール方法
Windows版のRは、管理者権限がないユーザーでもインストールできます。また、USBメモリーなどに入れて持ち歩き、常用パソコン以外でも利用できます。
「USBメモリへのRのインストール方法」は次のサイトで詳しく紹介されています。
[[http://humansci.let.hokudai.ac.jp/m/terao/stat/r/usb_install/usb_install.html|http://humansci.let.hokudai.ac.jp/m/terao/stat/r/usb_install/usb_install.html]]
拙作「rubyなしでExcel操縦rubyスクリプトを実行するためのexl.exe」に rrxwin.rb も組み込んだので、Rとexl.exeをUSBメモリーに入れておけば、常用パソコン以外でも利用できます。
[[http://cup.sakura.ne.jp/exl.htm|http://cup.sakura.ne.jp/exl.htm]]
--------
!2. RプログラムとRcmd.exeの使い方(rubyなしで行う方法)
以下では、クロス集計表(分割表)に対してカイ二乗検定を行う例を取り上げます。
サンプルの data01.csv には、結婚の意思に関して男女別に集計した表が盛り込まれています。次のような表です。
|| ||いずれ結婚するつもり||一生結婚しないつもり||不詳
||男子||3902||386||383
||女子||3402||284||268
この集計表から、結婚の意思について男女間で違いがあるかを確認するのにカイ二乗検定を用います。
分割表が data01.csv に書き込まれている場合、次の３行からなるRプログラムを実行すれば、カイ二乗検定の結果を得ることができます。
−−−− Rプログラムここから
dtf <- read.csv("data01.csv", header=T, row.names=1)
xx <- chisq.test(dtf)
xx
−−−− Rプログラムここまで
このRプログラム(名前をrprog01.rとします)を実行するには、Rcmd.exeを次のように起動します。
Rcmd.exe BATCH --no-save -q rprog01.r log.txt [enter]
こうすると、log.txtに検定結果などが書き出されます。
log.txtの中には、次のような行があるはずです。
X-squared = 10.5394, df = 2, p-value = 0.005145
これは、カイ二乗値、自由度、有意確立の値(p値)を示したものです。検定結果についてレポートなどに記載する際、これらの値を記します。
p値が0.05より小さければ、一般的に、有意性があるとみなされます。今回は0.005145と小さい値なので、有意性が認められることになります。つまり、結婚の意思に関して男女間で有意な違いがあるということになります。
以上のような Rcmd.exe の使い方を知っていれば、実行過程が記録されているlog.txtをrubyで読み込んで、必要な情報を得ることができます。
!!［参考］ 構成比(パーセント)を求める
いきなりカイ二乗検定を持ち出しましたが、本来、観測度数と構成比をきちんと確認した上で統計的検定を行うべきところです。
観測度数はcsvファイルに書かれているので、構成比を求めるRプログラムを参考まで載せておきます。
横方向の合計に対する比率と、縦方向の合計に対する比率の両方を表示します。
−−−− rprog02.r ここから
dtf <- read.csv("data01.csv", header=T, row.names=1)
mt <- as.matrix(dtf) # データフレームを行列に変換
xx1 <- prop.table(mt,1)*100 # 横方向の合計に対する比率
xx1
xx2 <- prop.table(mt,2)*100 # 縦方向の合計に対する比率
xx2
−−−− rprog02.r ここまで
Rプログラムは、行列を一括して処理できるところが魅力です。行列の各々の要素について反復的に計算する必要がありません。
余談ですが、上のRプログラムで構成比を出してみると、男女間で結婚の意思にそれほど違いがないように思われます。差があるとしても、僅かな差に過ぎないのでは？との印象を受けます。しかし、カイ２乗検定の結果では有意性が認められます。
こうした疑問がある時に、クラメール係数が２つの変数(「結婚の意思」と「性別」)の間の連関の強さに示唆を与えてくれます。クラメール係数については後述します。
ちなみに、男子・女子それぞれについてパーセンテージを算出すると次のとおり。
|| ||いずれ結婚するつもり||一生結婚しないつもり||不詳||計
||男子||83.5||8.3||8.2||100.0
||女子||86.0||7.2||6.8||100.0
--------
!3. rrxwin.rbを使ってRプログラムを実行(カイ２乗検定)
ここでは、Rプログラムによるカイ２乗検定を素材にして、rrxwin.rbの利用例を示します。
検定の結果、有意性が認められた時は残差分析も行います。
!!(1) 単純な実行例
Rプログラム rprog01.r と、データファイル data01.csv がカレントディレクトリにあるという前提で、その rprog01.r をrubyによって実行する例を示します。chi01.rbというrubyスクリプトです。
ruby.exe chi01.rb [enter]
のように実行すると、log.txtが書き出されます。これにRプログラムの実行過程＆結果が記録されています。
−−−− chi01.rb ここから
#! ruby -Ks
# 単純な実行例 (coding: Windows-31J)
require "rrxwin"
Rrx.rexec( "rprog01.r", "log.txt" )
puts "'rprog01.r' を実行し, 'log.txt' を作成しました."
−−−− chi01.rb ここまで
Rrx.rexec() がRプログラムを実行するためのメソッドです。
第１引数にはRプログラムを指定します。今回はファイル名をわたしましたが、Rプログラム本体(文字列)をわたすこともできます。
第２引数はログファイル名です。Rプログラムの実行過程＆結果の記録簿の名前です。これを省略するか、あるいは nil を指定すると、ログファイルは作成されません。
今回は指定していませんが、第３引数として Rcmd.exe の起動オプションをわたすことができます。"--vanilla --no-timing" のような文字列を第３引数としてわたすと、それが Rcmd.exe の起動オプションになります。第３引数を指定しない場合、"--no-save -q" が指定されたものとみなされます。
!!(2) rubyスクリプト中にRプログラムを埋め込む
chi02.rbは、rubyスクリプト中にRプログラムを埋め込む例です。
rprog01.rのようなRプログラムを別ファイルとして予め用意する必要がなくなります。
データファイルの data01.csv は、カレントディレクトリに置いてあるものとします。
−−−− chi02.rb ここから
#! ruby -Ks
# rubyスクリプト中にRプログラムを埋め込む (coding: Windows-31J)
require "rrxwin"
rpro = <<EOS
dtf <- read.csv("data01.csv", header=T, row.names=1)
xx <- chisq.test(dtf)
xx
EOS
Rrx.rexec(rpro, "log.txt")
puts "'log.txt' を作成しました."
−−−− chi02.rb ここまで
上のスクリプトをもう少し簡略化して、rproという変数を用いない形にすると次のようになります。
−−−− chi02b.rb ここから
#! ruby -Ks
# rubyスクリプト中にRプログラムを埋め込む・その２ (coding: Windows-31J)
require "rrxwin"
Rrx.rexec(<<EOS, "log.txt")
dtf <- read.csv("data01.csv", header=T, row.names=1)
xx <- chisq.test(dtf)
xx
EOS
puts "'log.txt' を作成しました."
−−−− chi02b.rb ここまで
!!(3) 検定結果を表示する(ログファイルの参照)
これまでのサンプルを実行した場合、log.txtの中に検定結果が書き出されています。
その検定結果に関する情報を表示する例を掲げます。ログファイルは残しません。
−−−− chi03.rb ここから
#! ruby -Ks
# 検定結果を表示する(ログファイルの参照) (coding: Windows-31J)
require "rrxwin"
rpro = <<EOS
dtf <- read.csv("data01.csv", header=T, row.names=1)
xx <- chisq.test(dtf)
xx
EOS
hs = Rrx.rexec(rpro)
print hs[:sink] if hs[:sink] != nil
−−−− chi03.rb ここまで
上のスクリプトを実行すると、カイ２乗検定の結果が標準出力に出力されます。ログファイルはできません。
「hs = Rrx.rexec( rpro )」は、Rプログラムの実行結果を変数 hs に代入するものです。
hsはハッシュです。何らかの理由でRプログラムが正常に実行されなかった時は、hsは、空のハッシュ {} になります。
hs[:sink] にはRプログラムの実行結果が文字列として記録されています。
また、hs[:log] というのもあります。この中には、Rプログラムとしてどんなコマンドが実行されたか(どんな関数が呼び出されたか)の記録と、その実行結果の記録の両方が含まれています。いわばRプログラムの実行過程が全部記録されているわけです。
それに対し、hs[:sink] にはコマンドの実行結果だけが記録されています。
Rプログラムの中で、例えば、sink("result.txt") という１行を置くと、それ以降のコマンド実行に関して、その実行結果だけが result.txt に書き出されます。どんなコマンドが実行されたかの記録は含まれません。
実は、Rrx.rexec() によってRプログラムを実行するとき、指定されたRプログラムの先頭に sink(……) を置いて、実行結果をテンポラリファイルに書き出すようにしています。そして、そのテンポラリファイルの内容を hs[:sink] に代入します。
!!(4) 調整済み残差を得る(暗黙のうちのテンポラリファイル利用)
カイ２乗検定を行った結果、p値が 0.05 より小さければ、一般的にいって有意性が認められることになります。
次のステップとして、表のどの欄の値が有意に大きい｜小さいかをみることになります。男女それぞれの「いずれ結婚するつもり」 「一生結婚しないつもり」 「不詳」のうち、どれが多いのか、どれが少ないのかをみます。
そのために、調整済み残差を確認します。これまでのRプログラムに即していうと、xx$stdres で調整済み残差を参照できます。
この調整済み残差をrubyスクリプトで参照するサンプル chi04.rb を掲げます。
−−−− chi04.rb ここから
#! ruby -Ks
# 調整済み残差の取得(暗黙のテンポラリファイル利用) (coding: Windows-31J)
require "rrxwin"
rpro = <<EOS
dtf <- read.csv("data01.csv", header=T, row.names=1)
xx <- chisq.test(dtf)
xx
write.csv(xx$stdres, file="@@stdres@@")
EOS
hs = Rrx.rexec(rpro)
if !hs[:sink] or !hs["stdres"]
exit
end
regex = /^\w\S* = (.+?), df = (.+?), p-value [<>=]+ (.+)$/
chi,df,pval = regex =~ hs[:sink] ? [$1.to_f, $2.to_i, $3.to_f] : [nil]*3
File.open("chitest.csv", "w") do |ff|
ff.printf("\"χ2 %g\",\"df %d\",\"p-value %g\"\n", chi, df, pval)
if pval < 0.05
ff.puts "\"調整済み残差\""
ff.print hs["stdres"]
end
end
puts "'chitest.csv' を作成しました."
−−−− chi04.rb ここまで
上のスクリプトでは、Rプログラムの最後のところに
write.csv(xx$stdres, file="@@stdres@@")
と書きました。「write.csv」は、指定された表データをcsv形式で書き出すためのRの関数です。本来の記述スタイルで記すなら次のようになります。
write.csv(xx$stdres, file="stdres.csv")
こうすると、調盛済み残差の表が stdres.csv という名前のファイルとして書き出されます。
今回、「"stdres.csv"」でなく「"@@stdres@@"」としているのは、次のような事情によります。
rrxwin.rbは、Rプログラムの中に「@@……@@」という記述があると、それをランダムな名前(テンポラリファイルの名前)に置き換えます。その上でそのRプログラムを実行します。
Rプログラムの実行が終了すると、テンポラリファイルの中身を読み込みます。
読み込んだ中身は、ハッシュに格納されます。ハッシュのキーは「@@……@@」の「……」の部分です。ハッシュは、rexec()の戻り値として返されます。
サンプルの場合でいうと、「@@stdres@@」の中身が hs["stdres"] として得られることになります。
テンポラリファイルのうち、暗黙のうちに作成されたものは、Rプログラムの実行が終了した時点で削除されます。
ユーザーが意図的に作成したテンポラリファイル(temp_makeメソッドで作成したもの)は、rubyスクリプトが終了する時点で削除されます。
テンポラリファイルの扱いは、rubyの標準添付ライブラリ tempfile に委ねています。
!!［参考］ 調整済み残差、期待値、残差、各欄のp値
調整済み残差の各欄の値が 1.96 より大きければ、その欄の値は、期待値よりも有意に大きいことになります。一方、-1.96より小さければ、有意に小さいことになります。その場合の有意水準は 5% です。有意性があると考えられるが、その結論が誤っている危険率が5%だという意味です。
2.58よりも大きい、あるいは -2.58よりも小さければ、やはり期待値よりも有意に大きい｜小さいということになりますが、その場合の有意水準は 1% です。1.96の時よりも危険率が低いことになります。
今回のサンプルでは、男子の「いずれ結婚するつもり」が -3.22 なので有意に小さく、一方、女子のそれは 3.22 なので有意に大きいと考えられます。
「一生結婚しないつもり」は、男子が 1.87、女子が -1.87 なので、男子の方が多い傾向ではありますが、有意性は認められません。
期待値は、男女間で結婚の意思に差がないと仮定した時に論理的に導き出される値です。Rプログラムでは xx$expected で参照できます。
xx$observed が観測度数、xx$residuals が残差(調整済み残差ではない)です。
また、各欄のp値を得たい時は、Rプログラムの中で
write.csv(pnorm(abs(xx$stdres),lower.tail=F)*2, file="@@p-values@@")
などのように記述すれば参照できます。これで得られるp値は絶対値なので、期待値よりも大きいか小さいかを判断する材料にはなりませんが、0.05より小さければ有意水準5%で有意性が認められる、などと判断する手がかりになります。
!!(5) csvファイルなどから表(行列)を切り取る cut_table, cut_num
Rプログラムで行列を処理する場合、その中に空欄があったり、あるいは、合計欄があると、トラブルの原因になることがあります。
そこで、rrxwin.rbでは、空欄を含まない行列を切り取るための cut_table() を設けました。矩形選択を自動化するものです。
このメソッドは、行列の左と上に項目名があると、それも一緒に切り取ります。
また、行列の中に合計欄があると判断される場合は、その合計欄の行または列を削除します。
項目名を一緒に切り取るのが難しそうな時は cut_num() を使うのが無難です。これは、数値だけから構成される行列を切り取ります。
２つの表(項目名の付いた行列)が含まれている data02.csv を材料にして、サンプルを示します。
data02.csvには、結婚の意思を年齢階層別に集計した表が入っています。男子と女子を別々に集計してあります。したがって、２つの表があります。
この２つの表をカイ２乗検定にかけるrubyスクリプトを下に掲げます。
−−−− chi05.rb ここから
#! ruby -Ks
# １つのcsvファイルから複数の行列を切り取って処理 (coding: Windows-31J)
require "rrxwin"
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", header=T, row.names=1)
dtf
xx <- chisq.test(dtf)
xx
EOS
# 行列を一つづつ取り出して処理
count = 0
Rrx.cut_table("data02.csv") do |tbl|
count += 1
Rrx.temp_make(tbl) # csvデータの一時ファイルを作成
logname = sprintf("log%02d.txt", count)
hs = Rrx.rexec(rpro, logname)
next unless hs[:sink]
regex = /^\w\S* = (.+?), df = (.+?), p-value [<>=]+ (.+)$/
chi,df,pval = regex =~ hs[:sink] ? [$1.to_f, $2.to_i, $3.to_f] : [nil]*3
printf("検定#%d\n", count)
printf("χ2 %g, df %d, p-value %g\n\n", chi, df, pval)
end
−−−− chi05.rb ここまで
Rプログラムの中で、「read.csv("@@last@@", ……」というのがあります。
この「@@last@@」は、rubyスクリプト内で作られたテンポラリファイルの名前に置き換えられます。rrxwin.rbの temp_make() メソッドを利用して作られたテンポラリファイルであって、最新のものの名前が当てはめられます。
サンプルの場合、「@@last@@」は、一つの行列が書き込まれたファイルの名前に置換されます。
「Rrx.cut_table("data02.csv")」は、data02.csvを読み込んで、その中から数値の行列を切り取ります。項目名が付いていればそれも合わせて切り取ります。
data02.csvの中身をみると、左端が空欄になっていますが、それらは含めません。
また、男子に関する表の場合、最初の行と左端の列が「総数」になっていますが、これは除去されます。
女子の場合は、最後の行と右端の列が「総数」ですが、これも除去されます。
cut_table() は、与えられたデータの中から、まず数値によって構成される矩形(２×２以上の正方形または長方形)の行列を探索します。そして、そのすぐ上にある行と、すぐ左側の列を項目名とみなして一緒に切り取ります。項目名が検出できない場合(数値の行列が上端とか左端に位置する場合)は、項目名を付けずに行列を切り取ります。
次に、そうやって切り取った表に対して合計欄の削除を試みます。合計欄が右端の列あるいは下端の行にあるかどうか調べます。実際に各要素を足し算してみて、確かに合計欄ということになれば削除します。
右端も下端も合計欄でない場合、今度は左端と上端が合計欄かどうかをチェックします。そして、合計欄であれば削除します。
１つの行列を検出する度に、行列本体(項目名を含む)を返します。サンプルではブロック付きで呼び出しているので、レシーバーの tbl にそれが代入されます。
ブロックなしで呼び出すた時は、複数の行列からなる配列が返されます。具体的には次のようなものです。
[tbl1, tbl2, tbl3, ……]
cut_num() も同じ挙動をしますが、項目名があることを想定せず、数値の行列(２×２以上)だけを切り取ります。１行だけとか１列だけの行列は切り取りません。
!!［補足１］ cut_table(), cut_num() の引数
これまでのサンプルでは、cut_table() の引数に "data02.csv" というファイル名をわたしていましたが、その他に、データそのもの(文字列)をわたすこともできます。
データは、csv形式またはタブ区切りテキストの２種類をわたすことができます。csvをわたした時は、戻り値の行列もcsvの文字列になり、タブ区切りテキストをわたした場合は、やはりタブ区切りテキスト(文字列)が戻り値になります。ブロック付きで呼び出した時は、そのレシーバーにそれらが代入されます。
csvかタブ区切りかの識別は、データの中にカンマが多く含まれているか、タブコードの方が多く含まれているかによって判別します。
引数にファイル名をわたした時は、そのファイルの中身がcsvかタブ区切りテキストかを同じように識別して、それに応じて戻り値を設定します。
また、引数に配列(Array)をわたすこともできます。その場合、戻り値の行列も配列になります。文字列でなく配列が返されます。
引数は、一つだけでなく複数わたすことができます。
その他、rubyのライブラリである spreadsheet または拙作 exlap がrequireされている場合は、Excelファイルの名前をわたすこともできます。例えば次のように書けます。
Rrx.cut_table("data.xls")
spreadsheetの場合は、パソコンにExcelがインストールされていなくても大丈夫ですが、扱えるのは *.xls ファイルのみです。
exlapの場合は、Excelがインストールされていないと使えませんが、*.xls, *.xml, *.htm, *.xlsx, *.xlsm, *.ods を扱えます。*.xmlは Excel2003用のxmlスプレッドシート、*.odsはオープンドキュメントスプレッドシートです。パソコンに導入されているのがExcel2007よりも前のものだと、扱えない拡張子のファイルがあるので注意して下さい。
!!［補足２］ cut_num() を使った時の項目名の付加
cut_num() は、項目名なしで数値の行列だけを切り取ります。項目名を付加したい場合は、Rプログラムの中で colnames(), rownames() を使うのが便利だと思います。
dtfがデータフレームである場合、次のようにします。
colnames(dtf) <- c("横1", "横2", "横3")
rownames(dtf) <- c("縦1", "縦2", "縦3")
上のようにすれば、例えば次のような表(データフレーム)になります。
|| ||横1||横2||横3
||縦1||10||20||30
||縦2||40||50||60
||縦3||70||80||90
!!(6) クリップボードの利用(str2clp, clp2str)
Rプログラムにおいて、ファイル名を記述するところに clipboard と書けば、ファイル経由でなくクリップボード経由で入出力を行えます。
read.csv("clipboard", ……)
write.csv(xx, "clipboard")
などのように記述します。
一方、rrxwin.rbには、文字列をクリップボードに書き出すメソッド str2clp(別名 s2c)と、その逆に、クリップボードのデータを文字列として取得するメソッド clp2str(別名 c2s)があります。例えば次のように用います。
str1 = "abc"
Rrx.str2clp(str1) # クリップボードへの書き出し
str2 = Rrx.clp2str() # クリップボードの読み込み
p str2 # => "abc"
chi05.rbをクリップボード利用の形に書き換えたのが chi06.rb です。
cut_table() の２番目の引数として数値の 2 をわたしていますが、これについてはスクリプトの後で説明します。
−−−− chi06.rb ここから
#! ruby -Ks
# クリップボードの利用 (coding: Windows-31J)
require "rrxwin"
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("clipboard", header=T, row.names=1)
dtf
xx <- chisq.test(dtf)
xx
EOS
# 行列を一つづつ取り出して処理
count = 0
Rrx.cut_table("data02.csv", 2) do |tbl, title|
count += 1
Rrx.str2clp(tbl) # 抽出した行列をクリップボードに書き出す
logname = sprintf("log%02d.txt", count)
hs = Rrx.rexec(rpro, logname)
next unless hs[:sink]
regex = /^\w\S* = (.+?), df = (.+?), p-value [<>=]+ (.+)$/
chi,df,pval = regex =~ hs[:sink] ? [$1.to_f, $2.to_i, $3.to_f] : [nil]*3
printf("検定#%d(%s)\n", count, title)
printf("χ2 %g, df %d, p-value %g\n\n", chi, df, pval)
end
−−−− chi06.rb ここまで
今回、cut_table() の２番目の引数として、数値の 2 を与えています。
実は、cut_table() の引数として、最後に 1, 2, 3 のどれかを与えると、戻り値の個数、またはブロック指定時のレシーバーの個数が変わります。
res = Rrx.cut_table("data02.csv", 2)
とした場合の戻り値は次のようになります。
[[tbl1, title1], [tbl2, title2], ……]
このとき、title1 には tbl1 の左上端の辺りにある表のタイトルと思われる文字列がセットされます。それらしい文字列がみつからない場合は "" の空文字列がセットされます。
サンプルの data02.csv の場合でいうと、title1 には "結婚の意思|男子" が入り、title2 には "女子" が入ります。
一方、数値の 3 を与えた場合
res = Rrx.cut_table("data02.csv", 3)
のようにすると、次のような戻り値になります。
[[tbl1, title1, ropt1], [tbl2, title2, ropt2], ……]
このとき、ropt1 には Rプログラムの read.csv() 用のオプションの一部がセットされます。"header=T, row.names=1" などの文字列がセットされます。
行列を切り取った時に、上端の項目名が検出された場合は「header=T」、検出されなかった時は「header=F」がセットされます。「row.names=1」は、左端の項目名が検出された時に入ります。
最後の引数として 1, 2, 3 の数値を指定しなかった時は、1 が指定されたものとみなされます。つまり title, ropt なしで tbl だけが返されます。
!!(7) グローバル変数によるRプログラムの置換
cut_table() で行列を切り取った場合、行列本体の他に title, ropt も返してくることを前述しました。その ropt を利用する例を掲げます。
roptは、Rプログラムの read.csv() にわたす引数の一部です。具体的には "header=T, row.names=1" のような文字列です。行列を切り取る際の項目名の検出具合によって変ってきます。
これをRプログラムに埋め込むことができれば、Rプログラムにおける行列の取り込みを自動化することにつながります。
rubyの sub() などの文字列置換機能を使えば、容易にそれを実現できます。例えば次のような形です。
rpro = "dtf <- read.csv("clipboard", __ropt__)\n ……"
rprog = rpro.sub(/__ropt__/, ropt)
Rrx.rexec(rprog, "log.txt")
上のような形でもいいのですが、rrxwin.rbでは、もう少しだけ簡略化できるよう、「グローバル変数によるRプログラムの置換」を行えるようにしてあります。
Rプログラムの中に「@@$ropt@@」という記述があると、その部分は、rubyスクリプト側のグローバル変数 $ropt の値に置き換えられます。
上述の例を次のように書くことができます。
rpro = "dtf <- read.csv("clipboard", @@$ropt@@)\n ……"
$ropt = ropt
Rrx.rexec(rpro, "log.txt")
さほど簡略化になっていないような気がしますが、もっと複雑なRプログラムで、複数個の置換文字列を使う時は、それなりに便利だと思います。
この場合のグローバル変数の名前は、ドル記号の次の文字が必ず半角アルファベット小文字でなければなりません。そのような制約条件を設けています。
グローバル変数の値は、配列とかハッシュでもかまいません。つまり「@@$ary[0]@@」とか「@@$hs{'x'}@@」といった記述でもかまいません。
もちろん、Rrx.rexec()が呼び出される直前の時点で、該当するグローバル変数の値がきちんと定義されている必要があります。
以下に、この「グローバル変数によるRプログラムの置換」のサンプルを掲げておきます。
−−−− chi07.rb ここから
#! ruby -Ks
# グローバル変数によるRプログラムの置換 (coding: Windows-31J)
require "rrxwin"
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", @@$ropt@@)
dtf
xx <- chisq.test(dtf)
xx
EOS
# 行列を一つづつ取り出して処理
count = 0
Rrx.cut_table("data02.csv", 3) do |tbl, title, ropt|
count += 1
Rrx.temp_make(tbl) # 抽出した行列を一時ファイルに書き出す
$ropt = ropt # read.csv() のオプションをグローバル変数に代入
logname = sprintf("log%02d.txt", count)
hs = Rrx.rexec(rpro, logname)
next unless hs[:sink]
regex = /^\w\S* = (.+?), df = (.+?), p-value [<>=]+ (.+)$/
chi,df,pval = regex =~ hs[:sink] ? [$1.to_f, $2.to_i, $3.to_f] : [nil]*3
printf("検定#%d(%s)\n", count, title)
printf("χ2 %g, df %d, p-value %g\n\n", chi, df, pval)
end
−−−− chi07.rb ここまで
!!(8) 文字列と配列の相互変換(str2ary, ary2str)
これまでのサンプルでは、Rプログラムの出力結果を文字列の形で取得してきました。例えば、調整済み残差も文字列として得ました。この文字列を配列に変換できれば、rubyスクリプトの側でいろいろ処理することが可能になります。
そこで、文字列を配列に変換するメソッド str2ary() を設けました。別名 s2a() です。次のように用います。
ary = Rrx.str2ary(str)
ary = Rrx.str2ary(str, "\t")
２番目の例は第２引数を指定していますが、タブ区切りテキストの文字列を配列に変換するという意味です。
１番目のように第２引数を指定しなければ、カンマ区切りなのかタブ区切りなのかを適当に判断します。文字列中にカンマとタブのどちらが多く含まれているかで判断します。
逆に、配列を文字列に変換するメソッド ary2str() もあります。別名 a2s() です。次のように用います。
str = Rrx.ary2str(ary)
str = Rrx.ary2str(ary, "\t")
２番目の例は第２引数を指定していますが、タブ区切りテキストに変換するためのものです。
１番目のように第２引数を指定しない時は、カンマ区切りのcsvの文字列に変換します。
これら str2ary() と ary2str() は、rubyの標準添付ライブラリcsvを用いて変換を行っています。
ただし、str2ary() の場合、ライブラリcsvで変換した結果(どのセルも文字列)に対して、数値化できるセルは数値に変換します。整数と考えられるものは to_i、浮動小数点数は to_f で変換を施した上で配列を返します。
どのセルも文字列で得たい時は、ライブラリcsvを用いて変換して下さい。例えば次のようにします。
sep = "," # フィールドの区切り文字
ary = (RUBY_VERSION < "1.9") ? CSV.parse(str, sep) :
CSV.parse(str, {:col_sep=>sep})
以上は Rrx のメソッドの説明でした。これを利用するサンプルを下に掲げます。
data02.csvを読み込んで、調整済み残差を確認します。そして、有意に多い｜少ない欄を抽出し、ちょっとしたレポート風に出力します。
−−−− chi08.rb ここから
#! ruby -Ks
# 文字列を配列に変換：調整済み残差のレポート (coding: Windows-31J)
require "rrxwin"
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", @@$ropt@@)
xx <- chisq.test(dtf)
xx
write.csv(xx$stdres, file="@@stdres@@")
EOS
# 行列を一つづつ取り出して処理
count = 0
Rrx.cut_table("data02.csv", 3) do |tbl, title, ropt|
count += 1
Rrx.temp_make(tbl) # 抽出した行列を一時ファイルに書き出す
$ropt = ropt
hs = Rrx.rexec(rpro)
next if !hs[:sink] or !hs["stdres"]
regex = /^\w\S* = (.+?), df = (.+?), p-value [<>=]+ (.+)$/
chi,df,pval = regex =~ hs[:sink] ? [$1.to_f, $2.to_i, $3.to_f] : [nil]*3
printf("検定#%d(%s)\n", count, title)
printf("χ2 %g\tdf %d\tp-value %g\n", chi, df, pval)
if pval < 0.05 # 有意性が認められる場合
ary = Rrx.str2ary(hs["stdres"]) # 調整済み残差を配列に変換
col = ary[0] # 上端(横方向の項目名)
row = Rrx.table_turn(ary)[0] # 左端(縦方向の項目名)
more = [] # 有意に多いセルをこれに記録
less = [] # 有意に少ないセルをこれに記録
for i in 1...ary.size
for j in 1...ary[i].size
val = ary[i][j]
more << [row[i], col[j], val] if val > 1.96
less << [row[i], col[j], val] if val < -1.96
end
end
puts "* 有意に多い欄"
more.each {|a| puts "\t" + a.join("\t")}
puts "* 有意に少ない欄"
less.each {|a| puts "\t" + a.join("\t")}
end
printf("\n")
end
−−−− chi08.rb ここまで
Rrx.table_turn() というのが出てきますが、これは、引数として与えられた配列を縦・横逆転させた結果を返すメソッドです。左端の縦１列が上端の横１行になる形で返ってきます。引数として与えた元の配列は変更されません。元のままです。
一般に、rubyスクリプトで配列を扱う場合、縦方向を扱うのは少々面倒です。それに対し、横方向を扱うのは容易です。なので、縦・横逆転させて、縦方向のものを横方向に転換すると、処理が簡単になることがあります。そこで、rrxwin.rbに table_turn() というメソッドを設けました。
後半は、調整済み残差の配列の要素(セル)を一つづつチェックして、有意に多いものと、有意に少ないものを別々の配列に記録し、最後にそれらを出力しています。
!!(9) Rプログラムのオブジェクトをrubyスクリプト側で一括して取得するための robj
Rプログラムの中で「xx <- chisq.test(dtf)」とすれば、カイ２乗検定の結果がxxにセットされます。xxの中にはカイ２乗値、p値、調整済み残差等が記録されています。
このxxのようなRプログラムのオブジェクトをrubyスクリプトの方で簡単に参照できるようにするための仕組みを設けてあります。Rプログラムの中で robj() という関数を用います。
data02.csvを読み込んでカイ２乗検定を行う例を示してみます。検定結果が収められているxxの情報をrubyスクリプト側で一覧表示します。
−−−− chi09.rb ここから
#! ruby -Ks
# robj() の利用例 (coding: Windows-31J)
require "rrxwin"
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", @@$ropt@@)
xx <- chisq.test(dtf)
robj(xx)
EOS
# 各々の表についてRプログラムを実行
count = 0
Rrx.cut_table("data02.csv", 3) do |tbl, title, ropt|
count += 1
Rrx.temp_make(tbl) # 抽出した表を一時ファイルに書き出す
$ropt = ropt
hs = Rrx.rexec(rpro)
next unless hs["robj1"]
printf("検定#%d(%s)\n", count, title)
xx = hs["robj1"]
xx.each do |key, val|
printf("%s\t", key)
p val
end
printf("\n")
end
−−−− chi09.rb ここまで
Rプログラムの中の robj(xx) は、xxをテンポラリファイルに書き出すものです。robj() の関数定義は、rrxwin.rbが自動的に付加します。
テンポラリファイルの中身は、rubyスクリプト側の rexec() の中で読み込まれて、その戻り値 hs にセットされます。hs["robj1"] でそれを参照できます。
「xx = hs["robj1"]」によって、Rプログラムにおけるxxオブジェクトが、rubyスクリプトの変数xxにセットされます。
Rプログラムの中では、例えば xx$p.value でp値を参照できますが、rubyスクリプトでは xx["p.value"] で参照できます。
rubyスクリプト側で調整済み残差をみたい時は、xx["stdres"] を確認します。これには配列(２次元配列)がセットされているはずです。
xx["statistic"] がカイ２乗値、xx["parameter"] が自由度、xx["method"] が検定の方法です。検定の方法には "Pearson's Chi-squared test" という文字列が代入されているはずです。
Rプログラムの中で robj() が何度か呼び出されている場合、その中身をrubyの方で参照する時は、hs["robj1"], hs["robj2"], hs["robj3"] …… のように参照します。
Rプログラムの中で「robj(xx, "chitest")」のように書くと、つまり第２引数を指定すると、rubyの側では hs["chitest"] でそれを参照できます。
「robj(xx, "@@chitest@@")」などのように「@@」を付けると正しく機能しません。read.csv() とか write.csv() の場合は「@@」を付けましたが、それとは違うので注意して下さい。
robj(xx) の引数xxには、リスト、データフレーム、マトリクス、ベクトル、数値、文字列のどのデータ型をわたしても大丈夫です。
引数xxに何もname属性がない場合や、name属性があるとしてもデータフレームとかマトリクスのように一つのまとまりとして扱う方がいい場合(ruby側で２次元配列として扱える場合)、xxを一つのまとまりとして出力します。
この場合、rubyの方で「xx = hs["robj1"]」とすると、ハッシュであるxxには xx["self"] という要素しかありません。この xx["self"] には、文字列、数値、配列のどれかがセットされているはずです。
ハッシュxxのキーとして "self" がある場合は、それ以外の要素がないことを意味しています。
例えば、Rプログラムの側で robj(xx$stdres) のように調整済み残差を引きわたした場合、rubyの方でその調整済み残差をみるには次のようにします。
xx = hs["robj1"]
ary = xx["self"]
こうすると、調整済み残差(２次元配列)が ary にセットされます。
!!(10) クラメール係数
カイ２乗検定を材料にしたrrxwin.rb利用の締めくくりとして、クラメール係数に触れておきます。
クラメール係数は、クロス集計表の横方向の項目と縦方向の項目の間の連関の強さを示すものです。0〜1の間の値をとります。
その数値をどう解釈すべきか一概にはいえませんが、一例を示すと次のようになります。
0.8≦v： 非常に強い相関、関連がある
0.5≦v＜0.8： やや強い相関
0.25≦v＜0.5： やや弱い相関
v＜0.25： 非常に弱い相関、関連がない
クラメール係数の値は、論文によっては「効果量」として記載されることがあるようです。
Rプログラムにおいて、クラメール係数は次のように求めることができます。
cramer_v <- sqrt(chisq.test(x)$statistic/(min(nrow(x), ncol(x)-1)/sum(x))
上の計算式で、xは、Rプログラムでいうところのmatrixでなければなりません。
以上のことを踏まえて、data02.csvを素材にして、カイ２乗検定の結果とクラメール係数を出力するスクリプトを示します。
−−−− chi10.rb ここから
#! ruby -Ks
# カイ２乗検定の結果とクラメール係数を表示 (coding: Windows-31J)
require "rrxwin"
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", @@$ropt@@)
mt <- as.matrix(dtf) # データフレームを行列に変換
xx <- chisq.test(mt)
xx$cramer <- sqrt(xx$statistic / (min(nrow(mt), ncol(mt)) - 1) / sum(mt))
robj(xx)
EOS
# 各々の表についてRプログラムを実行
count = 0
Rrx.cut_table("data02.csv", 3) do |tbl, title, ropt|
count += 1
Rrx.temp_make(tbl) # 抽出した表を一時ファイルに書き出す
$ropt = ropt
hs = Rrx.rexec(rpro)
next unless hs["robj1"]
xx = hs["robj1"]
printf("検定#%d(%s)\n", count, title)
names = %w(chi df p-value cramer)
%w(statistic parameter p.value cramer).each do |key|
printf("\t%s\t%g\n", names.shift, xx[key])
end
printf("\n")
end
−−−− chi10.rb ここまで
Rプログラムにおけるxxオブジェクトは、いわゆるリスト型です。データフレーム、マトリクス、ベクトル、数値、文字列等のいろいろな型のデータをこれに入れることができます。
クラメール係数の値をxxに追加したい時は、例えば次のようにします。cramer_vにクラメール係数の値が既にセットされているものとします。２通りの記述形式を示します。
xx$cramer <- cramer_v
xx[["cramer"]] <- cramer_v
上のようにしてから robj(xx) を書いておけば、rubyの側で xx["cramer"] でクラメール係数の値を参照できます。
余談ですが、検定の結果をみると、結婚の意思と年齢との連関の強さは、男子のクラメール係数が 0.158、女子は 0.242 です。どちらも弱い連関ですが、比較すると女子の方が高い値になっています。
ちなみに、最初の方で取り上げた data01.csv (男女別と結婚の意思)のクラメール係数は 0.035 という非常に低い値です。p値が 0.0051 なので有意性が検出されるわけですが、性別と結婚の意思の間に関連性があるといっていいのかどうか。パーセンテージをみても僅かな違いに過ぎないとの印象です。これをどんなニュアンスでレポートするかは、書き手の判断ということになるでしょうか。
!!［補足］ Rプログラムにおいてオブジェクトの属性等を知る方法
カイ２乗検定の結果が xx に入っている場合、xx$stdres で調整済み残差、xx$expected で期待値を参照できるわけですが、ここに出てくる「stdres」とか「expected」という名前をどうやって知ることができるかというと、Rプログラムの中で attributes(xx) を実行すると、それら名前が表示されます。これは、オブジェクトの属性一覧を表示するものです。
また、属性の名前だけでなくその値も知りたいという場合は、「dput(xx, "output.txt")」を実行して、後から output.txt の中身をみれば、そこに属性や値がずらずらと書かれています。
dput() は、本来、オブジェクトをテキストファイルで保存するためのものです。dput() で保存したオブジェクトは、別の機会に dget() で取り込んで、オブジェクトを復元することができます。
テキストファイルといっても、人が読んで分かりやすい形ではありませんが、「こんな構造になっているのか」と理解する手がかりになります。
参考まで attributes(), dput() の利用例を下に掲げておきます。data01.csv を読み込んでカイ２乗検定を行うケースです。
−−−−
#! ruby -Ks
require "rrxwin"
rpro = <<EOS
dtf <- read.csv("data01.csv", header=T, row.names=1)
xx <- chisq.test(dtf)
attributes(xx)
dput(xx, "@@dput@@")
EOS
hs = Rrx.rexec( rpro )
puts "* 属性一覧"
print hs[:sink] + "\n"
puts "* dputの出力結果"
puts hs["dput"]
−−−−
以上、かなり回りくどい説明になってしまいましたが、カイ２乗検定を材料にした説明は終わりにします。
rrxwin.rbの使い方の大半は、紹介できたと思います。
この後は、ｔ検定と相関係数を材料にして、説明しきれていないrrxwin.rbの使い方に触れたいと思います。
--------
!4. データの分類と平均値の差の検定(ｔ検定)
ここでは、２つのグループの平均値(例えば男女それぞれの身長など)について差が認められるかどうかの検定を取り上げます。
なお、以下のrubyスクリプトでは、最初の方に「include Rrx」という１行を置いて、「Rrx.」という記述を省略するようにしました。これまで Rrx.rexec() と書いていたのを単に rexec() と書きます。
!!(1) 数値以外のデータを含む行列の切り取り cut_data
data03.csv は、氏名、性別、身長、体重の４項目からなる表です。このうち氏名と性別は数値ではありません。具体的なデータの中身は次のようなものです。
氏名,性別,身長,体重
安部,男,158.9,55.6
伊藤,女,151.0,45.7
これまで表(行列)の切り取りに使ってきた cut_table() は、数値からなる行列を対象にしていました。その上と左に項目名があれば一緒に切り取りますが、あくまで数値の行列が対象です。しかし、今回は数値でない欄を含む行列データなので、このメソッドは使えません。
このような場合は cut_data() を用います。これは、数値でない欄も含めて、矩形(２×２以上の正方形または長方形)の行列を切り取ります。空欄を含まない行列を切り取ります。
行列の中に、数値を含まない行や列があるかもしれないので、その行または列が項目名なのかデータなのかを機械的に判別することはできません。その辺は利用者が判断して、適当に処理する必要があります。
data03.csv の中にある行列は、上端に横１行の項目名がありますが、左端の縦１列は項目名でなくデータです。Rプログラムでこれを読み込む時は下のようにします。
dtf <- read.csv("@@last@@", header=T)
下に data03.csv を読み込んで、そのデータフレームの構造を表示するサンプルを掲げます。
−−−− ft01.rb ここから
#! ruby -Ks
# データフレームの構造を表示 (coding: Windows-31J)
require "rrxwin"
include Rrx
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", header=T)
str(dtf) # dtfの構造を出力
EOS
# 行列の切り取りとRプログラムの実行
cut_data("data03.csv") do |data|
temp_make(data) # 一時ファイルへの行列の書き出し
hs = rexec(rpro)
next unless hs[:sink]
print hs[:sink]
end
−−−− ft01.rb ここまで
上のRプログラムの中に出てくる「str(dtf)」は、データフレームの構造を表示するものです。strは、structureの最初の３文字だと思います。
dtfの構造が次のように表示されます。
'data.frame': 10 obs. of 4 variables:
$ 氏名: Factor w/ 10 levels "安部","伊藤",..: 1 2 7 3 6 4 5 8 9 10
$ 性別: Factor w/ 2 levels "女","男": 2 1 2 1 2 1 2 1 2 1
$ 身長: num 160 151 178 164 170 ...
$ 体重: num 56.6 50.7 73.9 70.6 67.7 62.1 71.5 51.5 68.2 51.3
氏名と性別は、素材データでは数値でなく文字列なわけですが、Rプログラムでデータフレームとして読み込むと、levelsとして数が割り当てられていることが分かります。
!!［補足］ cut_data() の処理速度
cut_data() あるいは、カイ２乗検定のところで取り上げた cut_table(), cut_num() もそうですが、行列の切り取りには時間がかかるケースがあります。行列の規模が大きくなるほど、要する時間が相乗的に増えます。切り取りのアルゴリズムを単純なままにしていて、効率性を追求していないためです。
カイ２乗検定で扱うのはクロス集計表なので、行列の規模がそれほど大きくなることはないと思いますが、ｆ検定やｔ検定で扱うデータは、大きなものになることがあります。
そのような場合は、cut_data() を用いず、素材データをRプログラムで扱えるよう予め調整し、それを read.csv() で直接読み込むのがいいと思います。
!!(2) データの分類
ここで目的としているのは、男性の身長と女性の身長に差があるといえるかどうかを検定することです。
その検定を行う前提として、素材データの中から、男性だけのデータを取り出して１つのグループにし、また、それとは別に、女性のデータを取り出して別のグループに分けなければなりません。
この分類作業をrubyで行ってもいいのですが、Rプログラムだと簡単に行えます。次のようにします。
dtf <- read.csv("@@last@@", header=T)
male <- subset(dtf, 性別 == "男")
female <- subset(dtf, 性別 == "女")
上をRプログラムとして実行すると、男性だけからなるデータフレームの maleと、女性だけからなるデータフレームの female が作られます。
男性のデータと女性のデータをrubyの側で受け取って、csvデータとして出力するサンプルを掲げてみます。
−−−− ft02.rb ここから
#! ruby -Ks
# 性別による分類 (coding: Windows-31J)
require "rrxwin"
include Rrx
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", header=T)
male <- subset(dtf, 性別 == "男")
female <- subset(dtf, 性別 == "女")
write.csv(male, file="@@male@@")
write.csv(female, file="@@female@@")
EOS
# 行列を切り出してRプログラムを実行
cut_data("data03.csv") do |data|
temp_make(data) # 抽出した表を一時ファイルに書き出す
hs = rexec(rpro)
next if !hs["male"] or !hs["female"]
File.open("output.csv", "w") do |ff|
ff.print hs["male"] + "\n"
ff.print hs["female"] + "\n"
end
end
puts "'output.csv' を作成しました."
−−−− ft02.rb ここまで
上のスクリプトを実行すると、左端の縦１列に通し番号が入ってしまいます。それを削除したい時は、rrxwin.rbの中で定義されている del_col() を用います。これは、指定された列を削除するもので、次のように用います。
del_col(ary, 0) # 第１列目(左端の縦１列)を削除
del_col(ary, 1, 2) # 第２列目と第３列目を削除
第１引数は、２次元配列以外に、csvデータ(文字列)、タブ区切りテキスト(文字列)を指定することができます。配列を与えれば配列が、文字列を与えれば文字列が返されます。
２番目以降の引数には数値を与えます。複数の列を削除したい時に、必要な個数だけ列番号を指定します。列番号は 0 がスタートの値です。
del_col には del_tate, del_vertical という別名が割り当てられています。
del_row() というのもあり、こちらは指定行を削除します。del_yoko, del_horizontal という別名が割り当てられています。
ここには掲げませんが、del_col() の利用サンプルをft02b.rb として同梱しておきます。
!!(3) ｆ検定とｔ検定
本題の統計検定に話を進めます。
これまでのRプログラムで、男性のデータは male、女性のデータは female に入れることができました。それぞれの身長のデータは「male$身長」 「female$身長」で参照できます。
ここで、男女の身長に関して、まずは等分散かどうかを確認します。等分散か否かによって、平均値の検定の仕方が違ってくるためです。
等分散かどうかは「var.test(male$身長, female$身長)」で行います。これをｆ検定というようです。その実行結果を下に記しておきます。
−−−−
F test to compare two variances
data: male$身長 and female$身長
F = 0.9098, num df = 4, denom df = 4, p-value = 0.9292
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.09472154 8.73777711
sample estimates:
ratio of variances
0.9097558
−−−−
詳細は省きますが、上の表示の中で、とりあえず注目されるのは p-value です。この値が 0.05 より小さいと、「二つのグループの分散が等しいとみなすのは難しい。」ということになります。等分散とはいいきれない、つまり、ばらつきの度合いが異なるだろうということになります。
今回の表示をみると、p値は 0.9292 とかなり高い値なので、等分散と判断できます。
等分散である場合、平均値に差があるか否かの検定(ｔ検定)は、次のようにして行います。
t.test(male$身長, female$身長, var.equal=T) # Studentのt検定
ちなみに、等分散でない時は次のようにします。
t.test(male$身長, female$身長, var.equal=F) # Aspin-welchのt検定
以下に、ｆ検定とｔ検定を行うサンプルを掲げます。
−−−− ft03.rb ここから
#! ruby -Ks
# ｆ検定とｔ検定 (coding: Windows-31J)
require "rrxwin"
include Rrx
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", header=T)
male <- subset(dtf, 性別 == "男")
female <- subset(dtf, 性別 == "女")
ff <- var.test(male$身長, female$身長)
ff
flag <- (ff$p.value >= 0.05)
t.test(male$身長, female$身長, var.equal=flag)
EOS
# 行列を切り出してRプログラムを実行
cut_data("data03.csv") do |data|
temp_make(data) # 抽出した表を一時ファイルに書き出す
hs = rexec(rpro)
next unless hs[:sink]
print hs[:sink]
end
−−−− ft03.rb ここまで
上の実行結果をみると、男性の身長の平均値が 170.66、女性は 161.68 です。
ただ、ｔ検定のp値が 0.1887 です。これは 0.05 より大きいので、「男女の身長に有意な差があるとは言いきれない」ということになります。
!!(4) 身長と体重の検定を行う
data03.csvには、身長だけでなく体重のデータもあります。この両方についてｆ検定とｔ検定を行うことを考えます。
Rプログラムの中に「身長」と書き込んでしまうと、それを体重に応用することができません。そこで、「グローバル変数によるRプログラムの置換」のところで紹介したやり方を応用します。
まず、Rプログラムの中で「male$身長」と書いたところを「male$@@$item@@」と記述しておきます。
そして、rubyスクリプトの方で rexec() によってRプログラムを実行する直前に、グローバル変数 $item に "身長" とか "体重" を代入します。そうすると、Rプログラムの中の「@@$item@@」が「身長」とか「体重」に置き換えられた上で実行されます。
以下、そのサンプルです。
−−−− ft04.rb ここから
#! ruby -Ks
# 身長と体重の両方を検定 (coding: Windows-31J)
require "rrxwin"
include Rrx
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", header=T)
male <- subset(dtf, 性別 == "男")
female <- subset(dtf, 性別 == "女")
ff <- var.test(male$@@$item@@, female$@@$item@@)
ff
flag <- (ff$p.value >= 0.05)
t.test(male$@@$item@@, female$@@$item@@, var.equal=flag)
EOS
# 行列を切り出してRプログラムを実行
cut_data("data03.csv") do |data|
temp_make(data) # 抽出した表を一時ファイルに書き出す
%w(身長 体重).each do |item|
$item = item
hs = rexec(rpro)
next unless hs[:sink]
print hs[:sink]
end
end
−−−− ft04.rb ここまで
上のサンプルを実行してみると、身長に関するｔ検定の結果は、先述したように有意差が認められませんが、体重については、p値が 0.02153 と 0.05 より小さいので、男女の体重に有意な違いがあるということになります。
ちなみに、平均体重は、男性が 66.18、女性が 53.04 です。
ここで少し効率性を問題にしてみます。
上のスクリプトでは、Rプログラムを２度実行しています。したがって、当然ながら read.csv() も２度実行されます。
本来、データの読み込みは１回だけ行い、検定は２度行う(身長と体重について行う)というのが妥当なやり方です。その意味で、上のスクリプトは効率的でないといえます。
ということで、妥当と思われる方法の一例を下に掲げておきます。
−−−− ft05.rb ここから
#! ruby -Ks
# 身長と体重の両方を検定・その２ (coding: Windows-31J)
require "rrxwin"
include Rrx
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", header=T)
male <- subset(dtf, 性別 == "男")
female <- subset(dtf, 性別 == "女")
for (i in 1:ncol(dtf)) {
item <- colnames(dtf)[i]
if (item != "身長" && item != "体重") next
cat("**「", item, "」についての検定結果\n", sep="")
ff <- var.test(male[[item]], female[[item]])
print(ff)
flag <- (ff$p.value >= 0.05)
tt <- t.test(male[[item]], female[[item]], var.equal=flag)
print(tt)
}
EOS
# 行列を切り出してRプログラムを実行
cut_data("data03.csv") do |data|
temp_make(data) # 抽出した表を一時ファイルに書き出す
hs = rexec(rpro)
next unless hs[:sink]
print hs[:sink]
end
−−−− ft05.rb ここまで
ここには掲げませんが、上の ft05.rb を robj() 利用の形に書き換えた ft05b.rb を同梱しておきます。必要に応じて参考にして下さい。
以上でｆ検定とｔ検定を題材にした説明をおわりにします。
--------
!5. 相関係数と散布図の作成
ここでは、先でも利用した10人分の身長と体重のデータを使って、相関係数を算出します。また、身長と体重の散布図を描く方法について記します。
!!(1) 相関係数の算出
相関係数は、例えば、身長の高い人ほど体重も重い傾向がみられるか、といったことを検証する時に利用します。
0〜1 または -1〜0 の間の値を取り、相関係数の絶対値が0に近いほど二つの間の関連性が弱く、1に近いほど強いことになります。1に近ければ正比例のような関係、-1に近ければ逆比例のような関係といえます。
Rプログラムでは cor.test() で相関係数を求めることができます。
「xx <- cor.test(dtf$身長, dtf$体重)」とすれば、身長と体重の関係が得られます。
xx$estimate で相関係数を、xx$p.value でp値を参照できます。
cor.test() の結果を一通り出力するrubyスクリプトを掲げておきます。
−−−− cor01.rb ここから
#! ruby -Ks
# 相関係数の算出 (coding: Windows-31J)
require "rrxwin"
include Rrx
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", header=T)
xx <- cor.test(dtf$身長, dtf$体重)
robj(xx)
EOS
# 行列の切り取りとRプログラムの実行
cut_data("data03.csv") do |data|
temp_make(data) # 一時ファイルへの行列の書き出し
hs = rexec(rpro)
xx = hs["robj1"]
next unless xx
xx.each do |key, val|
printf("%s\t", key)
p val
end
end
−−−− cor01.rb ここまで
上のスクリプトを実行してみると、相関係数 xx["estimate"] が 0.8605302 と高い値です。
p値 xx["p.value"] は 0.001394111 となっており、0.05より小さいので、相関があると考えて差し支えなさそうです。
cor.test() には第３引数を与えることができます。相関係数の算出手法を指定するものです。次の３種類があります。
* ピアソンの積率相関係数 （pearson）
* ケンドールの順位相関係数 （kendall）
* スピアマンの順位相関係数 （spearman）
例えば次のように記述します。
cor.test(dtf$身長, dtf$体重, method="pearson")
cor.test(dtf$身長, dtf$体重, method="p") # 頭文字１文字のみ指定
第３引数を省略すると、「method="pearson"」が指定されたものとみなされます。
!!(2) 身長と体重の散布図の作成(pdfファイル)
data03.csvの身長と体重の関係を視覚的に把握するため、散布図を描いてみます。pdfファイルを作成します。
横軸を体重、縦軸を身長にして描きます。
Rプログラムでは plot() によってグラフや図を描くようです。以下にそのサンプルを掲げます。
−−−− cor02.rb ここから
#! ruby -Ks
# グラフをpdfに出力 (coding: Windows-31J)
require "rrxwin"
include Rrx
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", header=T)
pdf(file="cor02.pdf")
par(family = "Japan1GothicBBB") # 日本語フォントにゴチックを指定
par(mar = c(5, 5, 3, 3)) # 余白
plot(dtf$体重, dtf$身長, type="p",
main="身長と体重の散布図", cex.main=2, # メインタイトルとその大きさ
xlab="体重", ylab="身長", cex.lab=2)
dev.off()
EOS
# 行列の切り取りとRプログラムの実行
cut_data("data03.csv") do |data|
temp_make(data) # 一時ファイルへの行列の書き出し
hs = rexec(rpro)
print hs[:log]
end
−−−− cor02.rb ここまで
Rプログラムの pdf() は、pdfファイルの名前を指定して、それを出力先デバイスとして開くものです。
par() は、pdfファイルのパラメータを設定するもので、「family = ……」による日本語フォントの指定は欠かせません。これがないと、pdfに日本語が正しく書き込まれません。
plot() では様々なオプション指定ができるようですが、サンプルではメインタイトル、横軸と縦軸の項目名を指定しています。
dev.off() は、出力先デバイスとして開いたpdfファイルを閉じるものです。上のサンプルでは必要ないかもしれませんが、念のため置きました。
!!(3) MS-Windowsのメタファイル作成とExcelでの取り込み
MS-Windows専用の話題になりますが(linuxなどでは使えない)、Windowsの画像ファイルであるメタファイルを作成します。先のpdfファイルと同じ内容のものを作成します。
メタファイルは、拡張子が .emf のファイルで、Microsoft Office Picture Manager でそれを閲覧することができます。
このemfファイルは、Excel, Word, PowerPoint などに貼り付けることができます。下に掲げる cor03.rb では、Excelファイルに貼り付けています。
RプログラムでWindowsのメタファイルを作成するのは簡単で、pdfの時と同じように行います。パラメータの設定方法が少し違うだけです。
Excelを動かすのに、拙作 exlap.rb を使っています。
以下、サンプルです。
−−−− cor03.rb ここから
#! ruby -Ks
# グラフをemfに出力 (coding: Windows-31J)
require "exlap"
require "rrxwin"
include Rrx
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", header=T)
win.metafile(file="cor03.emf")
windowsFonts(JP1=windowsFont("MS Gothic"), JP2=windowsFont("MS Mincho"))
par(family="JP1")
plot(dtf$体重, dtf$身長, type="p",
main="身長と体重の散布図", cex.main=2, # メインタイトルとその大きさ
xlab="体重", ylab="身長", cex.lab=2)
dev.off()
EOS
# 行列の切り取りとRプログラムの実行
cut_data("data03.csv") do |data|
temp_make(data) # 一時ファイルへの行列の書き出し
hs = rexec(rpro)
print hs[:log]
end
# Excelファイルへの画像の貼り付け
emf_file = Exl::fullpath("cor03.emf")
Exlap.new("cor03.xls") do |wb|
ss = wb.fes
ss.Range("A1:B1").Value = "統計解析ソフトRの画像取込み",
"(男女10人の身長と体重)"
rng = ss.Range("A3:J13") # 描画領域。大きさはどうでもいい
shape = ss.Shapes.AddPicture({
'Filename'=>emf_file,
'LinkToFile'=>false, # 素材画像にリンクをはらない
'SaveWithDocument'=>true, # 素材画像を取り込んで保存
'Left'=>rng.Left,
'Top'=>rng.Top,
'Width'=>0,
'Height'=>0})
shape.ScaleWidth(1.0, true) # 素材画像と同じ幅にする(1.0倍)
shape.ScaleHeight(1.0, true) # 素材画像と同じ高さにする(1.0倍)
wb.save
end
puts "'cor03.xls' を作成しました."
−−−− cor03.rb ここまで
「rng = ss.Range("A3:J13")」というのは、画像を貼り付ける領域を示すものですが、実は、その大きさは適当でかまいません。ここでは、一応、10行×10列の合計100個のセルを画像領域にしていますが、大事なのは A3 の始点です。
AddPicture() で画像の貼り付けを行っているわけですが、画像領域の左端を rng.Left、上端を rng.Top にしています。これは A3 の位置を指すものです。Width, Height については、とりあえず0にしておいて、後から素材画像と同じ幅・高さになるよう調整しています。
--------
!6. 覚え書き
覚え書きを補足的に記しておきます。
!!(1) erbの利用
rubyには、テキストファイルにrubyスクリプトを埋め込むためのerbというライブラリが標準添付されています。そのため、雛形となるテキストを場面に応じて柔軟に変更できます。Rプログラムを雛形テキストに見立てれば、このerbを活用できます。
例えば、同じRプログラムを何度も実行したいけれども、その都度 read.csv() に与える引数を切り替えたい(読み込むcsvファイルを次々と変更したい)といったケースがあります。
この場合、erbを用いると次のようなスクリプトになります。サンプルの chi02.rb を変更したものを掲げてみます。
−−−− chi02c.rb ここから
#! ruby -Ks
# erbの利用例 (coding: Windows-31J)
require "erb"
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv(<%= read_param %>)
xx <- chisq.test(dtf)
xx
EOS
##
read_param = "\"data01.csv\", header=T, row.names=1"
erb = ERB.new(rpro)
Rrx.rexec(erb.result(binding), "log.txt")
puts "'log.txt' を作成しました."
−−−− chi02c.rb ここまで
上記は、ごく簡単なRプログラムなので erb を用いるまでもありません。rrxwin.rb の「グローバル変数によるRプログラムの置換」でも十分に対応できます。
ただ、もっと複雑なRプログラムになると、erbの利用が効果を発揮するケースがあると思います。
!!(2) read.table(), write.table() の利用
これまでのサンプルでは、read.csv(), write.csv() を使ってきました。ですが、read.table(), write.table() の方がより汎用的です。
参考まで chi04.rb の変更版を掲げてみます。カイ２乗検定の調整済み残差を出力するものです。
「sep=","」によって、フィールドセパレータ(区切り文字)を指定している点が異なります。
タブ区切りテキストを扱いたい時は、「sep="\t"」とします。
また、write.table() のオプションに、「col.names=NA」を付け加えています。これがないと、第１行目の横軸の項目名の列びが、左端の空欄がない状態になり、第２行目以降に比べて、フィールドの数が１個不足することになります。
各々のデータを引用符で囲まなくてもいい場合は、write.table() のオプションとして「quote=F」を加えます。
−−−− chi04b.rb ここから
#! ruby -Ks
# read.table, write.table の利用 (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.table("data01.csv", sep=",", header=T, row.names=1)
dtf
xx <- chisq.test(dtf)
xx
write.table(xx$stdres, file="stdres.csv", sep=",", col.names=NA)
EOS
Rrx.rexec(rpro, "log.txt")
puts "'log.txt' と 'stdres.csv' を作成しました."
−−−− chi04b.rb ここまで
!!(3) 文字コードの調整
仮にRプログラムがeuc-jpで書かれているとします。その場合、Rにそのことを伝えるためには、起動オプションとして --encoding=euc-jp を指定します。例えば次のようにします。
Rcmd.exe BATCH --encoding=euc-jp test.r log.txt [enter]
linuxなどでは次のようになります。
R CMD BATCH --encoding=euc-jp test.r log.txt [enter]
また、念のため Rプログラムの先頭のところに次の１行を置いておきます。どのOSにも共通です。
options(encoding="euc-jp")
このように文字コードを指定しておくと、Rプログラムの中に日本語が交じっていても正しく処理されます。読み込むファイルなどの文字コードもeuc-jpが想定されることになります。
rrxwin.rbは、実行中のrubyスクリプトの文字コードに従って、上の文字コード調整を自動的に行います。
実行中のrubyスクリプトの文字コードは、ruby 1.8系であれば $KCODE をみて判断し、ruby 1.9系であれば、rubyスクリプトのencode指定によって判断します。
Rrx.rencode() というメソッドがあり、これは前述のようにして判断した文字コードに変換するものです。例えば次のように用います。
str2 = Rrx.rencode(str)
こうすると、rubyスクリプトの文字コードがeuc-jpであれば、str2がeuc-jpの文字列として生成されます。
どの文字コードで書かれているか不明なファイルを読み込むような場合、この rencode() を適当に適用してやれば、トラブル回避になるかもしれません。
rencode()での文字コードの変換は nkf を用いています。
ちなみに、cut_table(), cut_num(), cut_data() では内部的に rencode() を呼び出しているので、わざわざ呼び出す必要はありません。
--------
以上で rrxwin.rb についての説明をおわりにします。
rrxwin.rb で定義されている各種メソッドなどについては、同梱の[[rrxwin02.txt|http://cup.sakura.ne.jp/rrxwin02.htm]]を参照して下さい。
〜 以上 〜
Copyright (C) T. Yoshiizumi, 2013 All rights reserved.
最終更新日： 2013/07/03
以下に掲げるドキュメントは、
[[rrx102.zip|http://cup.sakura.ne.jp/rrx102.zip]]
に同梱されている rrxwin.txt と同じ内容です。
--------
rubyから統計解析ソフトRを利用する方法について記します。
rubyとRの連携を取りやすくするためのライブラリ rrxwin.rb を作ったので、それを用います。
MS-Windows用のrrxwin.rbの他に、linuxなどで使えるrrx.rbもあります。それについては同梱の rrx.txt を参照して下さい。
以降で掲げるサンプルスクリプトを私が実行した環境は次のとおり。
* MS-Windows xp | vista
* ruby ver 1.8.7 | 1.9.3
* R ver 2.15.2
--------
{{toc_here}}
----
!＜はじめに＞
統計解析ソフトRは、とても高機能で統計解析のためのいろいろな仕組みを装備しています。また、Rにはプログラミング言語として様々な仕組みがあり、それに熟達している人にとっては、Rだけで処理できることが大半だろうと思います。
ただ、私は ruby を愛用しており、それと連携させてRを使いたいと考えました。rubyを経由することで、Rの解析結果を容易にExcelに貼り付けることができます。また、sqliteなど様々なデータベースに蓄積されたデータをRで解析するのも容易になります。
rubyからRを利用する強力なライブラリとして rsruby があります。これについても機会があれば書いてみたいと思いますが、もう少し手軽に使えるものとして rrxwin.rb を作ったので、その利用方法を記します。
１つのcsvファイルから複数の表(行列)を切り取るためのメソッドなども用意し、R活用の時の周辺的な便利さを追求したつもりです。
使用環境は MS-Windows です。
Rとrubyの両方が予めインストールされている必要があります。
文字コードは Windows-31J(cp932) を前提にしています。
Rの使い方には深入りしません。Rについては様々な解説がWebにあるのでそれらを参照して下さい。
サンプルではカイ２乗検定(クロス集計表の検定)、ｔ検定(平均値の差に関する検定)、相関係数の３つを取り上げます。また、相関係数に関連して散布図の描画も取り上げます。
--------
!1. rrxwin.rb のインストール
!!(1) インストールせずに使う方法
rrxwin.rb がカレントディレクトリにある場合は、rubyスクリプトの先頭の方に「require "./rrxwin"」を書いておけば、以下のインストール作業をしなくても使えます。
どのディレクトリからでも使えるようにしたい時は、下に掲げるインストール作業を行って下さい。
!!(2) インストール方法
rrxwin.rb, setup.rb の二つのファイルがカレントディレクトリにあることを確認して、次のようにrubyを実行します。
ruby.exe setup.rb [enter]
こうすると、rubyのライブラリパスが通ったディレクトリに rrxwin.rb がコピーされます。
このインストールの際に、統計解析ソフトRの実行ファイルである Rcmd.exe がどのディレクトリにあるかを検出すると、そのフルパス名を標準出力に出力します。また、検出できない時はその旨のメッセージを出力します。
検出できなかった場合、Rcmd.exeがパスの通ったディレクトリにないのであれば、そのフルパス名を利用者がセットしてやる必要があります。それは次のようなrubyスクリプトを実行して行います。
−−−− set_path.rb ここから
#! ruby -Ks
require "rrxwin"
pathname = "C:/usr/R/R-2.15.2/bin/i386/Rcmd.exe"
Rrx.path = pathname
−−−− set_path.rb ここまで
"C:/usr/R/……/Rcmd.exe" の部分は、ご自分のものに置き換えた上で実行して下さい。これを実行すると、そのフルパス名が登録されます(rrxwin.rbと同じディレクトリにあるrrxwin.cfgに書き込まれます。)。
Rをインストールする際、デフォルトのままで行うと次のようなフルパスになると思います。
C:/Program Files/R/R-2.15.2/bin/i386/Rcmd.exe
「2.15.2」は Rのバージョン番号です。どのRをインストールしたかによって違ってきます。
ユーザー権限等の関係で、"C:/Program Files" の下にディレクトリを作れない状況であれば、Documentディレクトリの下にあるかもしれません。例えば次のようなものです。
C:/Users/Taro/Documents/R/R-2.15.2/bin/i386/Rcmd.exe
!!(3) アンインストールの方法
上のようにしてインストールしたrrxwin.rbを削除したい時は unset.rb を実行します。
ruby.exe unset.rb [enter]
とすれば、rrxwin.rbとその関係ファイル(rrxwin.cfg)が削除されます。
unset.rbを実行する時は、カレントディレクトリに rrxwin.rb がなくてもかまいません。
!!［参考１］ rrxwin.rbが Rcmd.exe のパスを検索するやり方
rrxwin.rbは、アクセス可能なドライブ(ハードディスク、USBメモリーなど)を順番に検索しますが、各ドライブについて次のディレクトリがあるかどうか調べます。該当のものがみつかれば、そこで検索を中断します。
* /R (ルートディレクトリ直下の R というディレクトリ)
* /Program Files/R
* /Program Files (x86)/R
* /Users/Taro/Documents/R (Taro はユーザー名の一例)
上の４種類のディレクトリがどのドライブでもみつからないと、rrxwin.rbは、Rcmd.exeのパスを検出できません。その時は、前述した set_path.rb でパスを設定して下さい。
Rというディレクトリの下にバージョンの異なるRが複数ある場合は、最も新しいバージョンを自動的に選びます。
!!［参考２］ Rのインストール方法
Windows版のRは、管理者権限がないユーザーでもインストールできます。また、USBメモリーなどに入れて持ち歩き、常用パソコン以外でも利用できます。
「USBメモリへのRのインストール方法」は次のサイトで詳しく紹介されています。
[[http://humansci.let.hokudai.ac.jp/m/terao/stat/r/usb_install/usb_install.html|http://humansci.let.hokudai.ac.jp/m/terao/stat/r/usb_install/usb_install.html]]
拙作「rubyなしでExcel操縦rubyスクリプトを実行するためのexl.exe」に rrxwin.rb も組み込んだので、Rとexl.exeをUSBメモリーに入れておけば、常用パソコン以外でも利用できます。
[[http://cup.sakura.ne.jp/exl.htm|http://cup.sakura.ne.jp/exl.htm]]
--------
!2. RプログラムとRcmd.exeの使い方(rubyなしで行う方法)
以下では、クロス集計表(分割表)に対してカイ二乗検定を行う例を取り上げます。
サンプルの data01.csv には、結婚の意思に関して男女別に集計した表が盛り込まれています。次のような表です。
|| ||いずれ結婚するつもり||一生結婚しないつもり||不詳
||男子||3902||386||383
||女子||3402||284||268
この集計表から、結婚の意思について男女間で違いがあるかを確認するのにカイ二乗検定を用います。
分割表が data01.csv に書き込まれている場合、次の３行からなるRプログラムを実行すれば、カイ二乗検定の結果を得ることができます。
−−−− Rプログラムここから
dtf <- read.csv("data01.csv", header=T, row.names=1)
xx <- chisq.test(dtf)
xx
−−−− Rプログラムここまで
このRプログラム(名前をrprog01.rとします)を実行するには、Rcmd.exeを次のように起動します。
Rcmd.exe BATCH --no-save -q rprog01.r log.txt [enter]
こうすると、log.txtに検定結果などが書き出されます。
log.txtの中には、次のような行があるはずです。
X-squared = 10.5394, df = 2, p-value = 0.005145
これは、カイ二乗値、自由度、有意確立の値(p値)を示したものです。検定結果についてレポートなどに記載する際、これらの値を記します。
p値が0.05より小さければ、一般的に、有意性があるとみなされます。今回は0.005145と小さい値なので、有意性が認められることになります。つまり、結婚の意思に関して男女間で有意な違いがあるということになります。
以上のような Rcmd.exe の使い方を知っていれば、実行過程が記録されているlog.txtをrubyで読み込んで、必要な情報を得ることができます。
!!［参考］ 構成比(パーセント)を求める
いきなりカイ二乗検定を持ち出しましたが、本来、観測度数と構成比をきちんと確認した上で統計的検定を行うべきところです。
観測度数はcsvファイルに書かれているので、構成比を求めるRプログラムを参考まで載せておきます。
横方向の合計に対する比率と、縦方向の合計に対する比率の両方を表示します。
−−−− rprog02.r ここから
dtf <- read.csv("data01.csv", header=T, row.names=1)
mt <- as.matrix(dtf) # データフレームを行列に変換
xx1 <- prop.table(mt,1)*100 # 横方向の合計に対する比率
xx1
xx2 <- prop.table(mt,2)*100 # 縦方向の合計に対する比率
xx2
−−−− rprog02.r ここまで
Rプログラムは、行列を一括して処理できるところが魅力です。行列の各々の要素について反復的に計算する必要がありません。
余談ですが、上のRプログラムで構成比を出してみると、男女間で結婚の意思にそれほど違いがないように思われます。差があるとしても、僅かな差に過ぎないのでは？との印象を受けます。しかし、カイ２乗検定の結果では有意性が認められます。
こうした疑問がある時に、クラメール係数が２つの変数(「結婚の意思」と「性別」)の間の連関の強さに示唆を与えてくれます。クラメール係数については後述します。
ちなみに、男子・女子それぞれについてパーセンテージを算出すると次のとおり。
|| ||いずれ結婚するつもり||一生結婚しないつもり||不詳||計
||男子||83.5||8.3||8.2||100.0
||女子||86.0||7.2||6.8||100.0
--------
!3. rrxwin.rbを使ってRプログラムを実行(カイ２乗検定)
ここでは、Rプログラムによるカイ２乗検定を素材にして、rrxwin.rbの利用例を示します。
検定の結果、有意性が認められた時は残差分析も行います。
!!(1) 単純な実行例
Rプログラム rprog01.r と、データファイル data01.csv がカレントディレクトリにあるという前提で、その rprog01.r をrubyによって実行する例を示します。chi01.rbというrubyスクリプトです。
ruby.exe chi01.rb [enter]
のように実行すると、log.txtが書き出されます。これにRプログラムの実行過程＆結果が記録されています。
−−−− chi01.rb ここから
#! ruby -Ks
# 単純な実行例 (coding: Windows-31J)
require "rrxwin"
Rrx.rexec( "rprog01.r", "log.txt" )
puts "'rprog01.r' を実行し, 'log.txt' を作成しました."
−−−− chi01.rb ここまで
Rrx.rexec() がRプログラムを実行するためのメソッドです。
第１引数にはRプログラムを指定します。今回はファイル名をわたしましたが、Rプログラム本体(文字列)をわたすこともできます。
第２引数はログファイル名です。Rプログラムの実行過程＆結果の記録簿の名前です。これを省略するか、あるいは nil を指定すると、ログファイルは作成されません。
今回は指定していませんが、第３引数として Rcmd.exe の起動オプションをわたすことができます。"--vanilla --no-timing" のような文字列を第３引数としてわたすと、それが Rcmd.exe の起動オプションになります。第３引数を指定しない場合、"--no-save -q" が指定されたものとみなされます。
!!(2) rubyスクリプト中にRプログラムを埋め込む
chi02.rbは、rubyスクリプト中にRプログラムを埋め込む例です。
rprog01.rのようなRプログラムを別ファイルとして予め用意する必要がなくなります。
データファイルの data01.csv は、カレントディレクトリに置いてあるものとします。
−−−− chi02.rb ここから
#! ruby -Ks
# rubyスクリプト中にRプログラムを埋め込む (coding: Windows-31J)
require "rrxwin"
rpro = <<EOS
dtf <- read.csv("data01.csv", header=T, row.names=1)
xx <- chisq.test(dtf)
xx
EOS
Rrx.rexec(rpro, "log.txt")
puts "'log.txt' を作成しました."
−−−− chi02.rb ここまで
上のスクリプトをもう少し簡略化して、rproという変数を用いない形にすると次のようになります。
−−−− chi02b.rb ここから
#! ruby -Ks
# rubyスクリプト中にRプログラムを埋め込む・その２ (coding: Windows-31J)
require "rrxwin"
Rrx.rexec(<<EOS, "log.txt")
dtf <- read.csv("data01.csv", header=T, row.names=1)
xx <- chisq.test(dtf)
xx
EOS
puts "'log.txt' を作成しました."
−−−− chi02b.rb ここまで
!!(3) 検定結果を表示する(ログファイルの参照)
これまでのサンプルを実行した場合、log.txtの中に検定結果が書き出されています。
その検定結果に関する情報を表示する例を掲げます。ログファイルは残しません。
−−−− chi03.rb ここから
#! ruby -Ks
# 検定結果を表示する(ログファイルの参照) (coding: Windows-31J)
require "rrxwin"
rpro = <<EOS
dtf <- read.csv("data01.csv", header=T, row.names=1)
xx <- chisq.test(dtf)
xx
EOS
hs = Rrx.rexec(rpro)
print hs[:sink] if hs[:sink] != nil
−−−− chi03.rb ここまで
上のスクリプトを実行すると、カイ２乗検定の結果が標準出力に出力されます。ログファイルはできません。
「hs = Rrx.rexec( rpro )」は、Rプログラムの実行結果を変数 hs に代入するものです。
hsはハッシュです。何らかの理由でRプログラムが正常に実行されなかった時は、hsは、空のハッシュ {} になります。
hs[:sink] にはRプログラムの実行結果が文字列として記録されています。
また、hs[:log] というのもあります。この中には、Rプログラムとしてどんなコマンドが実行されたか(どんな関数が呼び出されたか)の記録と、その実行結果の記録の両方が含まれています。いわばRプログラムの実行過程が全部記録されているわけです。
それに対し、hs[:sink] にはコマンドの実行結果だけが記録されています。
Rプログラムの中で、例えば、sink("result.txt") という１行を置くと、それ以降のコマンド実行に関して、その実行結果だけが result.txt に書き出されます。どんなコマンドが実行されたかの記録は含まれません。
実は、Rrx.rexec() によってRプログラムを実行するとき、指定されたRプログラムの先頭に sink(……) を置いて、実行結果をテンポラリファイルに書き出すようにしています。そして、そのテンポラリファイルの内容を hs[:sink] に代入します。
!!(4) 調整済み残差を得る(暗黙のうちのテンポラリファイル利用)
カイ２乗検定を行った結果、p値が 0.05 より小さければ、一般的にいって有意性が認められることになります。
次のステップとして、表のどの欄の値が有意に大きい｜小さいかをみることになります。男女それぞれの「いずれ結婚するつもり」 「一生結婚しないつもり」 「不詳」のうち、どれが多いのか、どれが少ないのかをみます。
そのために、調整済み残差を確認します。これまでのRプログラムに即していうと、xx$stdres で調整済み残差を参照できます。
この調整済み残差をrubyスクリプトで参照するサンプル chi04.rb を掲げます。
−−−− chi04.rb ここから
#! ruby -Ks
# 調整済み残差の取得(暗黙のテンポラリファイル利用) (coding: Windows-31J)
require "rrxwin"
rpro = <<EOS
dtf <- read.csv("data01.csv", header=T, row.names=1)
xx <- chisq.test(dtf)
xx
write.csv(xx$stdres, file="@@stdres@@")
EOS
hs = Rrx.rexec(rpro)
if !hs[:sink] or !hs["stdres"]
exit
end
regex = /^\w\S* = (.+?), df = (.+?), p-value [<>=]+ (.+)$/
chi,df,pval = regex =~ hs[:sink] ? [$1.to_f, $2.to_i, $3.to_f] : [nil]*3
File.open("chitest.csv", "w") do |ff|
ff.printf("\"χ2 %g\",\"df %d\",\"p-value %g\"\n", chi, df, pval)
if pval < 0.05
ff.puts "\"調整済み残差\""
ff.print hs["stdres"]
end
end
puts "'chitest.csv' を作成しました."
−−−− chi04.rb ここまで
上のスクリプトでは、Rプログラムの最後のところに
write.csv(xx$stdres, file="@@stdres@@")
と書きました。「write.csv」は、指定された表データをcsv形式で書き出すためのRの関数です。本来の記述スタイルで記すなら次のようになります。
write.csv(xx$stdres, file="stdres.csv")
こうすると、調盛済み残差の表が stdres.csv という名前のファイルとして書き出されます。
今回、「"stdres.csv"」でなく「"@@stdres@@"」としているのは、次のような事情によります。
rrxwin.rbは、Rプログラムの中に「@@……@@」という記述があると、それをランダムな名前(テンポラリファイルの名前)に置き換えます。その上でそのRプログラムを実行します。
Rプログラムの実行が終了すると、テンポラリファイルの中身を読み込みます。
読み込んだ中身は、ハッシュに格納されます。ハッシュのキーは「@@……@@」の「……」の部分です。ハッシュは、rexec()の戻り値として返されます。
サンプルの場合でいうと、「@@stdres@@」の中身が hs["stdres"] として得られることになります。
テンポラリファイルのうち、暗黙のうちに作成されたものは、Rプログラムの実行が終了した時点で削除されます。
ユーザーが意図的に作成したテンポラリファイル(temp_makeメソッドで作成したもの)は、rubyスクリプトが終了する時点で削除されます。
テンポラリファイルの扱いは、rubyの標準添付ライブラリ tempfile に委ねています。
!!［参考］ 調整済み残差、期待値、残差、各欄のp値
調整済み残差の各欄の値が 1.96 より大きければ、その欄の値は、期待値よりも有意に大きいことになります。一方、-1.96より小さければ、有意に小さいことになります。その場合の有意水準は 5% です。有意性があると考えられるが、その結論が誤っている危険率が5%だという意味です。
2.58よりも大きい、あるいは -2.58よりも小さければ、やはり期待値よりも有意に大きい｜小さいということになりますが、その場合の有意水準は 1% です。1.96の時よりも危険率が低いことになります。
今回のサンプルでは、男子の「いずれ結婚するつもり」が -3.22 なので有意に小さく、一方、女子のそれは 3.22 なので有意に大きいと考えられます。
「一生結婚しないつもり」は、男子が 1.87、女子が -1.87 なので、男子の方が多い傾向ではありますが、有意性は認められません。
期待値は、男女間で結婚の意思に差がないと仮定した時に論理的に導き出される値です。Rプログラムでは xx$expected で参照できます。
xx$observed が観測度数、xx$residuals が残差(調整済み残差ではない)です。
また、各欄のp値を得たい時は、Rプログラムの中で
write.csv(pnorm(abs(xx$stdres),lower.tail=F)*2, file="@@p-values@@")
などのように記述すれば参照できます。これで得られるp値は絶対値なので、期待値よりも大きいか小さいかを判断する材料にはなりませんが、0.05より小さければ有意水準5%で有意性が認められる、などと判断する手がかりになります。
!!(5) csvファイルなどから表(行列)を切り取る cut_table, cut_num
Rプログラムで行列を処理する場合、その中に空欄があったり、あるいは、合計欄があると、トラブルの原因になることがあります。
そこで、rrxwin.rbでは、空欄を含まない行列を切り取るための cut_table() を設けました。矩形選択を自動化するものです。
このメソッドは、行列の左と上に項目名があると、それも一緒に切り取ります。
また、行列の中に合計欄があると判断される場合は、その合計欄の行または列を削除します。
項目名を一緒に切り取るのが難しそうな時は cut_num() を使うのが無難です。これは、数値だけから構成される行列を切り取ります。
２つの表(項目名の付いた行列)が含まれている data02.csv を材料にして、サンプルを示します。
data02.csvには、結婚の意思を年齢階層別に集計した表が入っています。男子と女子を別々に集計してあります。したがって、２つの表があります。
この２つの表をカイ２乗検定にかけるrubyスクリプトを下に掲げます。
−−−− chi05.rb ここから
#! ruby -Ks
# １つのcsvファイルから複数の行列を切り取って処理 (coding: Windows-31J)
require "rrxwin"
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", header=T, row.names=1)
dtf
xx <- chisq.test(dtf)
xx
EOS
# 行列を一つづつ取り出して処理
count = 0
Rrx.cut_table("data02.csv") do |tbl|
count += 1
Rrx.temp_make(tbl) # csvデータの一時ファイルを作成
logname = sprintf("log%02d.txt", count)
hs = Rrx.rexec(rpro, logname)
next unless hs[:sink]
regex = /^\w\S* = (.+?), df = (.+?), p-value [<>=]+ (.+)$/
chi,df,pval = regex =~ hs[:sink] ? [$1.to_f, $2.to_i, $3.to_f] : [nil]*3
printf("検定#%d\n", count)
printf("χ2 %g, df %d, p-value %g\n\n", chi, df, pval)
end
−−−− chi05.rb ここまで
Rプログラムの中で、「read.csv("@@last@@", ……」というのがあります。
この「@@last@@」は、rubyスクリプト内で作られたテンポラリファイルの名前に置き換えられます。rrxwin.rbの temp_make() メソッドを利用して作られたテンポラリファイルであって、最新のものの名前が当てはめられます。
サンプルの場合、「@@last@@」は、一つの行列が書き込まれたファイルの名前に置換されます。
「Rrx.cut_table("data02.csv")」は、data02.csvを読み込んで、その中から数値の行列を切り取ります。項目名が付いていればそれも合わせて切り取ります。
data02.csvの中身をみると、左端が空欄になっていますが、それらは含めません。
また、男子に関する表の場合、最初の行と左端の列が「総数」になっていますが、これは除去されます。
女子の場合は、最後の行と右端の列が「総数」ですが、これも除去されます。
cut_table() は、与えられたデータの中から、まず数値によって構成される矩形(２×２以上の正方形または長方形)の行列を探索します。そして、そのすぐ上にある行と、すぐ左側の列を項目名とみなして一緒に切り取ります。項目名が検出できない場合(数値の行列が上端とか左端に位置する場合)は、項目名を付けずに行列を切り取ります。
次に、そうやって切り取った表に対して合計欄の削除を試みます。合計欄が右端の列あるいは下端の行にあるかどうか調べます。実際に各要素を足し算してみて、確かに合計欄ということになれば削除します。
右端も下端も合計欄でない場合、今度は左端と上端が合計欄かどうかをチェックします。そして、合計欄であれば削除します。
１つの行列を検出する度に、行列本体(項目名を含む)を返します。サンプルではブロック付きで呼び出しているので、レシーバーの tbl にそれが代入されます。
ブロックなしで呼び出すた時は、複数の行列からなる配列が返されます。具体的には次のようなものです。
[tbl1, tbl2, tbl3, ……]
cut_num() も同じ挙動をしますが、項目名があることを想定せず、数値の行列(２×２以上)だけを切り取ります。１行だけとか１列だけの行列は切り取りません。
!!［補足１］ cut_table(), cut_num() の引数
これまでのサンプルでは、cut_table() の引数に "data02.csv" というファイル名をわたしていましたが、その他に、データそのもの(文字列)をわたすこともできます。
データは、csv形式またはタブ区切りテキストの２種類をわたすことができます。csvをわたした時は、戻り値の行列もcsvの文字列になり、タブ区切りテキストをわたした場合は、やはりタブ区切りテキスト(文字列)が戻り値になります。ブロック付きで呼び出した時は、そのレシーバーにそれらが代入されます。
csvかタブ区切りかの識別は、データの中にカンマが多く含まれているか、タブコードの方が多く含まれているかによって判別します。
引数にファイル名をわたした時は、そのファイルの中身がcsvかタブ区切りテキストかを同じように識別して、それに応じて戻り値を設定します。
また、引数に配列(Array)をわたすこともできます。その場合、戻り値の行列も配列になります。文字列でなく配列が返されます。
引数は、一つだけでなく複数わたすことができます。
その他、rubyのライブラリである spreadsheet または拙作 exlap がrequireされている場合は、Excelファイルの名前をわたすこともできます。例えば次のように書けます。
Rrx.cut_table("data.xls")
spreadsheetの場合は、パソコンにExcelがインストールされていなくても大丈夫ですが、扱えるのは *.xls ファイルのみです。
exlapの場合は、Excelがインストールされていないと使えませんが、*.xls, *.xml, *.htm, *.xlsx, *.xlsm, *.ods を扱えます。*.xmlは Excel2003用のxmlスプレッドシート、*.odsはオープンドキュメントスプレッドシートです。パソコンに導入されているのがExcel2007よりも前のものだと、扱えない拡張子のファイルがあるので注意して下さい。
!!［補足２］ cut_num() を使った時の項目名の付加
cut_num() は、項目名なしで数値の行列だけを切り取ります。項目名を付加したい場合は、Rプログラムの中で colnames(), rownames() を使うのが便利だと思います。
dtfがデータフレームである場合、次のようにします。
colnames(dtf) <- c("横1", "横2", "横3")
rownames(dtf) <- c("縦1", "縦2", "縦3")
上のようにすれば、例えば次のような表(データフレーム)になります。
|| ||横1||横2||横3
||縦1||10||20||30
||縦2||40||50||60
||縦3||70||80||90
!!(6) クリップボードの利用(str2clp, clp2str)
Rプログラムにおいて、ファイル名を記述するところに clipboard と書けば、ファイル経由でなくクリップボード経由で入出力を行えます。
read.csv("clipboard", ……)
write.csv(xx, "clipboard")
などのように記述します。
一方、rrxwin.rbには、文字列をクリップボードに書き出すメソッド str2clp(別名 s2c)と、その逆に、クリップボードのデータを文字列として取得するメソッド clp2str(別名 c2s)があります。例えば次のように用います。
str1 = "abc"
Rrx.str2clp(str1) # クリップボードへの書き出し
str2 = Rrx.clp2str() # クリップボードの読み込み
p str2 # => "abc"
chi05.rbをクリップボード利用の形に書き換えたのが chi06.rb です。
cut_table() の２番目の引数として数値の 2 をわたしていますが、これについてはスクリプトの後で説明します。
−−−− chi06.rb ここから
#! ruby -Ks
# クリップボードの利用 (coding: Windows-31J)
require "rrxwin"
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("clipboard", header=T, row.names=1)
dtf
xx <- chisq.test(dtf)
xx
EOS
# 行列を一つづつ取り出して処理
count = 0
Rrx.cut_table("data02.csv", 2) do |tbl, title|
count += 1
Rrx.str2clp(tbl) # 抽出した行列をクリップボードに書き出す
logname = sprintf("log%02d.txt", count)
hs = Rrx.rexec(rpro, logname)
next unless hs[:sink]
regex = /^\w\S* = (.+?), df = (.+?), p-value [<>=]+ (.+)$/
chi,df,pval = regex =~ hs[:sink] ? [$1.to_f, $2.to_i, $3.to_f] : [nil]*3
printf("検定#%d(%s)\n", count, title)
printf("χ2 %g, df %d, p-value %g\n\n", chi, df, pval)
end
−−−− chi06.rb ここまで
今回、cut_table() の２番目の引数として、数値の 2 を与えています。
実は、cut_table() の引数として、最後に 1, 2, 3 のどれかを与えると、戻り値の個数、またはブロック指定時のレシーバーの個数が変わります。
res = Rrx.cut_table("data02.csv", 2)
とした場合の戻り値は次のようになります。
[[tbl1, title1], [tbl2, title2], ……]
このとき、title1 には tbl1 の左上端の辺りにある表のタイトルと思われる文字列がセットされます。それらしい文字列がみつからない場合は "" の空文字列がセットされます。
サンプルの data02.csv の場合でいうと、title1 には "結婚の意思|男子" が入り、title2 には "女子" が入ります。
一方、数値の 3 を与えた場合
res = Rrx.cut_table("data02.csv", 3)
のようにすると、次のような戻り値になります。
[[tbl1, title1, ropt1], [tbl2, title2, ropt2], ……]
このとき、ropt1 には Rプログラムの read.csv() 用のオプションの一部がセットされます。"header=T, row.names=1" などの文字列がセットされます。
行列を切り取った時に、上端の項目名が検出された場合は「header=T」、検出されなかった時は「header=F」がセットされます。「row.names=1」は、左端の項目名が検出された時に入ります。
最後の引数として 1, 2, 3 の数値を指定しなかった時は、1 が指定されたものとみなされます。つまり title, ropt なしで tbl だけが返されます。
!!(7) グローバル変数によるRプログラムの置換
cut_table() で行列を切り取った場合、行列本体の他に title, ropt も返してくることを前述しました。その ropt を利用する例を掲げます。
roptは、Rプログラムの read.csv() にわたす引数の一部です。具体的には "header=T, row.names=1" のような文字列です。行列を切り取る際の項目名の検出具合によって変ってきます。
これをRプログラムに埋め込むことができれば、Rプログラムにおける行列の取り込みを自動化することにつながります。
rubyの sub() などの文字列置換機能を使えば、容易にそれを実現できます。例えば次のような形です。
rpro = "dtf <- read.csv("clipboard", __ropt__)\n ……"
rprog = rpro.sub(/__ropt__/, ropt)
Rrx.rexec(rprog, "log.txt")
上のような形でもいいのですが、rrxwin.rbでは、もう少しだけ簡略化できるよう、「グローバル変数によるRプログラムの置換」を行えるようにしてあります。
Rプログラムの中に「@@$ropt@@」という記述があると、その部分は、rubyスクリプト側のグローバル変数 $ropt の値に置き換えられます。
上述の例を次のように書くことができます。
rpro = "dtf <- read.csv("clipboard", @@$ropt@@)\n ……"
$ropt = ropt
Rrx.rexec(rpro, "log.txt")
さほど簡略化になっていないような気がしますが、もっと複雑なRプログラムで、複数個の置換文字列を使う時は、それなりに便利だと思います。
この場合のグローバル変数の名前は、ドル記号の次の文字が必ず半角アルファベット小文字でなければなりません。そのような制約条件を設けています。
グローバル変数の値は、配列とかハッシュでもかまいません。つまり「@@$ary[0]@@」とか「@@$hs{'x'}@@」といった記述でもかまいません。
もちろん、Rrx.rexec()が呼び出される直前の時点で、該当するグローバル変数の値がきちんと定義されている必要があります。
以下に、この「グローバル変数によるRプログラムの置換」のサンプルを掲げておきます。
−−−− chi07.rb ここから
#! ruby -Ks
# グローバル変数によるRプログラムの置換 (coding: Windows-31J)
require "rrxwin"
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", @@$ropt@@)
dtf
xx <- chisq.test(dtf)
xx
EOS
# 行列を一つづつ取り出して処理
count = 0
Rrx.cut_table("data02.csv", 3) do |tbl, title, ropt|
count += 1
Rrx.temp_make(tbl) # 抽出した行列を一時ファイルに書き出す
$ropt = ropt # read.csv() のオプションをグローバル変数に代入
logname = sprintf("log%02d.txt", count)
hs = Rrx.rexec(rpro, logname)
next unless hs[:sink]
regex = /^\w\S* = (.+?), df = (.+?), p-value [<>=]+ (.+)$/
chi,df,pval = regex =~ hs[:sink] ? [$1.to_f, $2.to_i, $3.to_f] : [nil]*3
printf("検定#%d(%s)\n", count, title)
printf("χ2 %g, df %d, p-value %g\n\n", chi, df, pval)
end
−−−− chi07.rb ここまで
!!(8) 文字列と配列の相互変換(str2ary, ary2str)
これまでのサンプルでは、Rプログラムの出力結果を文字列の形で取得してきました。例えば、調整済み残差も文字列として得ました。この文字列を配列に変換できれば、rubyスクリプトの側でいろいろ処理することが可能になります。
そこで、文字列を配列に変換するメソッド str2ary() を設けました。別名 s2a() です。次のように用います。
ary = Rrx.str2ary(str)
ary = Rrx.str2ary(str, "\t")
２番目の例は第２引数を指定していますが、タブ区切りテキストの文字列を配列に変換するという意味です。
１番目のように第２引数を指定しなければ、カンマ区切りなのかタブ区切りなのかを適当に判断します。文字列中にカンマとタブのどちらが多く含まれているかで判断します。
逆に、配列を文字列に変換するメソッド ary2str() もあります。別名 a2s() です。次のように用います。
str = Rrx.ary2str(ary)
str = Rrx.ary2str(ary, "\t")
２番目の例は第２引数を指定していますが、タブ区切りテキストに変換するためのものです。
１番目のように第２引数を指定しない時は、カンマ区切りのcsvの文字列に変換します。
これら str2ary() と ary2str() は、rubyの標準添付ライブラリcsvを用いて変換を行っています。
ただし、str2ary() の場合、ライブラリcsvで変換した結果(どのセルも文字列)に対して、数値化できるセルは数値に変換します。整数と考えられるものは to_i、浮動小数点数は to_f で変換を施した上で配列を返します。
どのセルも文字列で得たい時は、ライブラリcsvを用いて変換して下さい。例えば次のようにします。
sep = "," # フィールドの区切り文字
ary = (RUBY_VERSION < "1.9") ? CSV.parse(str, sep) :
CSV.parse(str, {:col_sep=>sep})
以上は Rrx のメソッドの説明でした。これを利用するサンプルを下に掲げます。
data02.csvを読み込んで、調整済み残差を確認します。そして、有意に多い｜少ない欄を抽出し、ちょっとしたレポート風に出力します。
−−−− chi08.rb ここから
#! ruby -Ks
# 文字列を配列に変換：調整済み残差のレポート (coding: Windows-31J)
require "rrxwin"
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", @@$ropt@@)
xx <- chisq.test(dtf)
xx
write.csv(xx$stdres, file="@@stdres@@")
EOS
# 行列を一つづつ取り出して処理
count = 0
Rrx.cut_table("data02.csv", 3) do |tbl, title, ropt|
count += 1
Rrx.temp_make(tbl) # 抽出した行列を一時ファイルに書き出す
$ropt = ropt
hs = Rrx.rexec(rpro)
next if !hs[:sink] or !hs["stdres"]
regex = /^\w\S* = (.+?), df = (.+?), p-value [<>=]+ (.+)$/
chi,df,pval = regex =~ hs[:sink] ? [$1.to_f, $2.to_i, $3.to_f] : [nil]*3
printf("検定#%d(%s)\n", count, title)
printf("χ2 %g\tdf %d\tp-value %g\n", chi, df, pval)
if pval < 0.05 # 有意性が認められる場合
ary = Rrx.str2ary(hs["stdres"]) # 調整済み残差を配列に変換
col = ary[0] # 上端(横方向の項目名)
row = Rrx.table_turn(ary)[0] # 左端(縦方向の項目名)
more = [] # 有意に多いセルをこれに記録
less = [] # 有意に少ないセルをこれに記録
for i in 1...ary.size
for j in 1...ary[i].size
val = ary[i][j]
more << [row[i], col[j], val] if val > 1.96
less << [row[i], col[j], val] if val < -1.96
end
end
puts "* 有意に多い欄"
more.each {|a| puts "\t" + a.join("\t")}
puts "* 有意に少ない欄"
less.each {|a| puts "\t" + a.join("\t")}
end
printf("\n")
end
−−−− chi08.rb ここまで
Rrx.table_turn() というのが出てきますが、これは、引数として与えられた配列を縦・横逆転させた結果を返すメソッドです。左端の縦１列が上端の横１行になる形で返ってきます。引数として与えた元の配列は変更されません。元のままです。
一般に、rubyスクリプトで配列を扱う場合、縦方向を扱うのは少々面倒です。それに対し、横方向を扱うのは容易です。なので、縦・横逆転させて、縦方向のものを横方向に転換すると、処理が簡単になることがあります。そこで、rrxwin.rbに table_turn() というメソッドを設けました。
後半は、調整済み残差の配列の要素(セル)を一つづつチェックして、有意に多いものと、有意に少ないものを別々の配列に記録し、最後にそれらを出力しています。
!!(9) Rプログラムのオブジェクトをrubyスクリプト側で一括して取得するための robj
Rプログラムの中で「xx <- chisq.test(dtf)」とすれば、カイ２乗検定の結果がxxにセットされます。xxの中にはカイ２乗値、p値、調整済み残差等が記録されています。
このxxのようなRプログラムのオブジェクトをrubyスクリプトの方で簡単に参照できるようにするための仕組みを設けてあります。Rプログラムの中で robj() という関数を用います。
data02.csvを読み込んでカイ２乗検定を行う例を示してみます。検定結果が収められているxxの情報をrubyスクリプト側で一覧表示します。
−−−− chi09.rb ここから
#! ruby -Ks
# robj() の利用例 (coding: Windows-31J)
require "rrxwin"
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", @@$ropt@@)
xx <- chisq.test(dtf)
robj(xx)
EOS
# 各々の表についてRプログラムを実行
count = 0
Rrx.cut_table("data02.csv", 3) do |tbl, title, ropt|
count += 1
Rrx.temp_make(tbl) # 抽出した表を一時ファイルに書き出す
$ropt = ropt
hs = Rrx.rexec(rpro)
next unless hs["robj1"]
printf("検定#%d(%s)\n", count, title)
xx = hs["robj1"]
xx.each do |key, val|
printf("%s\t", key)
p val
end
printf("\n")
end
−−−− chi09.rb ここまで
Rプログラムの中の robj(xx) は、xxをテンポラリファイルに書き出すものです。robj() の関数定義は、rrxwin.rbが自動的に付加します。
テンポラリファイルの中身は、rubyスクリプト側の rexec() の中で読み込まれて、その戻り値 hs にセットされます。hs["robj1"] でそれを参照できます。
「xx = hs["robj1"]」によって、Rプログラムにおけるxxオブジェクトが、rubyスクリプトの変数xxにセットされます。
Rプログラムの中では、例えば xx$p.value でp値を参照できますが、rubyスクリプトでは xx["p.value"] で参照できます。
rubyスクリプト側で調整済み残差をみたい時は、xx["stdres"] を確認します。これには配列(２次元配列)がセットされているはずです。
xx["statistic"] がカイ２乗値、xx["parameter"] が自由度、xx["method"] が検定の方法です。検定の方法には "Pearson's Chi-squared test" という文字列が代入されているはずです。
Rプログラムの中で robj() が何度か呼び出されている場合、その中身をrubyの方で参照する時は、hs["robj1"], hs["robj2"], hs["robj3"] …… のように参照します。
Rプログラムの中で「robj(xx, "chitest")」のように書くと、つまり第２引数を指定すると、rubyの側では hs["chitest"] でそれを参照できます。
「robj(xx, "@@chitest@@")」などのように「@@」を付けると正しく機能しません。read.csv() とか write.csv() の場合は「@@」を付けましたが、それとは違うので注意して下さい。
robj(xx) の引数xxには、リスト、データフレーム、マトリクス、ベクトル、数値、文字列のどのデータ型をわたしても大丈夫です。
引数xxに何もname属性がない場合や、name属性があるとしてもデータフレームとかマトリクスのように一つのまとまりとして扱う方がいい場合(ruby側で２次元配列として扱える場合)、xxを一つのまとまりとして出力します。
この場合、rubyの方で「xx = hs["robj1"]」とすると、ハッシュであるxxには xx["self"] という要素しかありません。この xx["self"] には、文字列、数値、配列のどれかがセットされているはずです。
ハッシュxxのキーとして "self" がある場合は、それ以外の要素がないことを意味しています。
例えば、Rプログラムの側で robj(xx$stdres) のように調整済み残差を引きわたした場合、rubyの方でその調整済み残差をみるには次のようにします。
xx = hs["robj1"]
ary = xx["self"]
こうすると、調整済み残差(２次元配列)が ary にセットされます。
!!(10) クラメール係数
カイ２乗検定を材料にしたrrxwin.rb利用の締めくくりとして、クラメール係数に触れておきます。
クラメール係数は、クロス集計表の横方向の項目と縦方向の項目の間の連関の強さを示すものです。0〜1の間の値をとります。
その数値をどう解釈すべきか一概にはいえませんが、一例を示すと次のようになります。
0.8≦v： 非常に強い相関、関連がある
0.5≦v＜0.8： やや強い相関
0.25≦v＜0.5： やや弱い相関
v＜0.25： 非常に弱い相関、関連がない
クラメール係数の値は、論文によっては「効果量」として記載されることがあるようです。
Rプログラムにおいて、クラメール係数は次のように求めることができます。
cramer_v <- sqrt(chisq.test(x)$statistic/(min(nrow(x), ncol(x)-1)/sum(x))
上の計算式で、xは、Rプログラムでいうところのmatrixでなければなりません。
以上のことを踏まえて、data02.csvを素材にして、カイ２乗検定の結果とクラメール係数を出力するスクリプトを示します。
−−−− chi10.rb ここから
#! ruby -Ks
# カイ２乗検定の結果とクラメール係数を表示 (coding: Windows-31J)
require "rrxwin"
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", @@$ropt@@)
mt <- as.matrix(dtf) # データフレームを行列に変換
xx <- chisq.test(mt)
xx$cramer <- sqrt(xx$statistic / (min(nrow(mt), ncol(mt)) - 1) / sum(mt))
robj(xx)
EOS
# 各々の表についてRプログラムを実行
count = 0
Rrx.cut_table("data02.csv", 3) do |tbl, title, ropt|
count += 1
Rrx.temp_make(tbl) # 抽出した表を一時ファイルに書き出す
$ropt = ropt
hs = Rrx.rexec(rpro)
next unless hs["robj1"]
xx = hs["robj1"]
printf("検定#%d(%s)\n", count, title)
names = %w(chi df p-value cramer)
%w(statistic parameter p.value cramer).each do |key|
printf("\t%s\t%g\n", names.shift, xx[key])
end
printf("\n")
end
−−−− chi10.rb ここまで
Rプログラムにおけるxxオブジェクトは、いわゆるリスト型です。データフレーム、マトリクス、ベクトル、数値、文字列等のいろいろな型のデータをこれに入れることができます。
クラメール係数の値をxxに追加したい時は、例えば次のようにします。cramer_vにクラメール係数の値が既にセットされているものとします。２通りの記述形式を示します。
xx$cramer <- cramer_v
xx[["cramer"]] <- cramer_v
上のようにしてから robj(xx) を書いておけば、rubyの側で xx["cramer"] でクラメール係数の値を参照できます。
余談ですが、検定の結果をみると、結婚の意思と年齢との連関の強さは、男子のクラメール係数が 0.158、女子は 0.242 です。どちらも弱い連関ですが、比較すると女子の方が高い値になっています。
ちなみに、最初の方で取り上げた data01.csv (男女別と結婚の意思)のクラメール係数は 0.035 という非常に低い値です。p値が 0.0051 なので有意性が検出されるわけですが、性別と結婚の意思の間に関連性があるといっていいのかどうか。パーセンテージをみても僅かな違いに過ぎないとの印象です。これをどんなニュアンスでレポートするかは、書き手の判断ということになるでしょうか。
!!［補足］ Rプログラムにおいてオブジェクトの属性等を知る方法
カイ２乗検定の結果が xx に入っている場合、xx$stdres で調整済み残差、xx$expected で期待値を参照できるわけですが、ここに出てくる「stdres」とか「expected」という名前をどうやって知ることができるかというと、Rプログラムの中で attributes(xx) を実行すると、それら名前が表示されます。これは、オブジェクトの属性一覧を表示するものです。
また、属性の名前だけでなくその値も知りたいという場合は、「dput(xx, "output.txt")」を実行して、後から output.txt の中身をみれば、そこに属性や値がずらずらと書かれています。
dput() は、本来、オブジェクトをテキストファイルで保存するためのものです。dput() で保存したオブジェクトは、別の機会に dget() で取り込んで、オブジェクトを復元することができます。
テキストファイルといっても、人が読んで分かりやすい形ではありませんが、「こんな構造になっているのか」と理解する手がかりになります。
参考まで attributes(), dput() の利用例を下に掲げておきます。data01.csv を読み込んでカイ２乗検定を行うケースです。
−−−−
#! ruby -Ks
require "rrxwin"
rpro = <<EOS
dtf <- read.csv("data01.csv", header=T, row.names=1)
xx <- chisq.test(dtf)
attributes(xx)
dput(xx, "@@dput@@")
EOS
hs = Rrx.rexec( rpro )
puts "* 属性一覧"
print hs[:sink] + "\n"
puts "* dputの出力結果"
puts hs["dput"]
−−−−
以上、かなり回りくどい説明になってしまいましたが、カイ２乗検定を材料にした説明は終わりにします。
rrxwin.rbの使い方の大半は、紹介できたと思います。
この後は、ｔ検定と相関係数を材料にして、説明しきれていないrrxwin.rbの使い方に触れたいと思います。
--------
!4. データの分類と平均値の差の検定(ｔ検定)
ここでは、２つのグループの平均値(例えば男女それぞれの身長など)について差が認められるかどうかの検定を取り上げます。
なお、以下のrubyスクリプトでは、最初の方に「include Rrx」という１行を置いて、「Rrx.」という記述を省略するようにしました。これまで Rrx.rexec() と書いていたのを単に rexec() と書きます。
!!(1) 数値以外のデータを含む行列の切り取り cut_data
data03.csv は、氏名、性別、身長、体重の４項目からなる表です。このうち氏名と性別は数値ではありません。具体的なデータの中身は次のようなものです。
氏名,性別,身長,体重
安部,男,158.9,55.6
伊藤,女,151.0,45.7
これまで表(行列)の切り取りに使ってきた cut_table() は、数値からなる行列を対象にしていました。その上と左に項目名があれば一緒に切り取りますが、あくまで数値の行列が対象です。しかし、今回は数値でない欄を含む行列データなので、このメソッドは使えません。
このような場合は cut_data() を用います。これは、数値でない欄も含めて、矩形(２×２以上の正方形または長方形)の行列を切り取ります。空欄を含まない行列を切り取ります。
行列の中に、数値を含まない行や列があるかもしれないので、その行または列が項目名なのかデータなのかを機械的に判別することはできません。その辺は利用者が判断して、適当に処理する必要があります。
data03.csv の中にある行列は、上端に横１行の項目名がありますが、左端の縦１列は項目名でなくデータです。Rプログラムでこれを読み込む時は下のようにします。
dtf <- read.csv("@@last@@", header=T)
下に data03.csv を読み込んで、そのデータフレームの構造を表示するサンプルを掲げます。
−−−− ft01.rb ここから
#! ruby -Ks
# データフレームの構造を表示 (coding: Windows-31J)
require "rrxwin"
include Rrx
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", header=T)
str(dtf) # dtfの構造を出力
EOS
# 行列の切り取りとRプログラムの実行
cut_data("data03.csv") do |data|
temp_make(data) # 一時ファイルへの行列の書き出し
hs = rexec(rpro)
next unless hs[:sink]
print hs[:sink]
end
−−−− ft01.rb ここまで
上のRプログラムの中に出てくる「str(dtf)」は、データフレームの構造を表示するものです。strは、structureの最初の３文字だと思います。
dtfの構造が次のように表示されます。
'data.frame': 10 obs. of 4 variables:
$ 氏名: Factor w/ 10 levels "安部","伊藤",..: 1 2 7 3 6 4 5 8 9 10
$ 性別: Factor w/ 2 levels "女","男": 2 1 2 1 2 1 2 1 2 1
$ 身長: num 160 151 178 164 170 ...
$ 体重: num 56.6 50.7 73.9 70.6 67.7 62.1 71.5 51.5 68.2 51.3
氏名と性別は、素材データでは数値でなく文字列なわけですが、Rプログラムでデータフレームとして読み込むと、levelsとして数が割り当てられていることが分かります。
!!［補足］ cut_data() の処理速度
cut_data() あるいは、カイ２乗検定のところで取り上げた cut_table(), cut_num() もそうですが、行列の切り取りには時間がかかるケースがあります。行列の規模が大きくなるほど、要する時間が相乗的に増えます。切り取りのアルゴリズムを単純なままにしていて、効率性を追求していないためです。
カイ２乗検定で扱うのはクロス集計表なので、行列の規模がそれほど大きくなることはないと思いますが、ｆ検定やｔ検定で扱うデータは、大きなものになることがあります。
そのような場合は、cut_data() を用いず、素材データをRプログラムで扱えるよう予め調整し、それを read.csv() で直接読み込むのがいいと思います。
!!(2) データの分類
ここで目的としているのは、男性の身長と女性の身長に差があるといえるかどうかを検定することです。
その検定を行う前提として、素材データの中から、男性だけのデータを取り出して１つのグループにし、また、それとは別に、女性のデータを取り出して別のグループに分けなければなりません。
この分類作業をrubyで行ってもいいのですが、Rプログラムだと簡単に行えます。次のようにします。
dtf <- read.csv("@@last@@", header=T)
male <- subset(dtf, 性別 == "男")
female <- subset(dtf, 性別 == "女")
上をRプログラムとして実行すると、男性だけからなるデータフレームの maleと、女性だけからなるデータフレームの female が作られます。
男性のデータと女性のデータをrubyの側で受け取って、csvデータとして出力するサンプルを掲げてみます。
−−−− ft02.rb ここから
#! ruby -Ks
# 性別による分類 (coding: Windows-31J)
require "rrxwin"
include Rrx
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", header=T)
male <- subset(dtf, 性別 == "男")
female <- subset(dtf, 性別 == "女")
write.csv(male, file="@@male@@")
write.csv(female, file="@@female@@")
EOS
# 行列を切り出してRプログラムを実行
cut_data("data03.csv") do |data|
temp_make(data) # 抽出した表を一時ファイルに書き出す
hs = rexec(rpro)
next if !hs["male"] or !hs["female"]
File.open("output.csv", "w") do |ff|
ff.print hs["male"] + "\n"
ff.print hs["female"] + "\n"
end
end
puts "'output.csv' を作成しました."
−−−− ft02.rb ここまで
上のスクリプトを実行すると、左端の縦１列に通し番号が入ってしまいます。それを削除したい時は、rrxwin.rbの中で定義されている del_col() を用います。これは、指定された列を削除するもので、次のように用います。
del_col(ary, 0) # 第１列目(左端の縦１列)を削除
del_col(ary, 1, 2) # 第２列目と第３列目を削除
第１引数は、２次元配列以外に、csvデータ(文字列)、タブ区切りテキスト(文字列)を指定することができます。配列を与えれば配列が、文字列を与えれば文字列が返されます。
２番目以降の引数には数値を与えます。複数の列を削除したい時に、必要な個数だけ列番号を指定します。列番号は 0 がスタートの値です。
del_col には del_tate, del_vertical という別名が割り当てられています。
del_row() というのもあり、こちらは指定行を削除します。del_yoko, del_horizontal という別名が割り当てられています。
ここには掲げませんが、del_col() の利用サンプルをft02b.rb として同梱しておきます。
!!(3) ｆ検定とｔ検定
本題の統計検定に話を進めます。
これまでのRプログラムで、男性のデータは male、女性のデータは female に入れることができました。それぞれの身長のデータは「male$身長」 「female$身長」で参照できます。
ここで、男女の身長に関して、まずは等分散かどうかを確認します。等分散か否かによって、平均値の検定の仕方が違ってくるためです。
等分散かどうかは「var.test(male$身長, female$身長)」で行います。これをｆ検定というようです。その実行結果を下に記しておきます。
−−−−
F test to compare two variances
data: male$身長 and female$身長
F = 0.9098, num df = 4, denom df = 4, p-value = 0.9292
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.09472154 8.73777711
sample estimates:
ratio of variances
0.9097558
−−−−
詳細は省きますが、上の表示の中で、とりあえず注目されるのは p-value です。この値が 0.05 より小さいと、「二つのグループの分散が等しいとみなすのは難しい。」ということになります。等分散とはいいきれない、つまり、ばらつきの度合いが異なるだろうということになります。
今回の表示をみると、p値は 0.9292 とかなり高い値なので、等分散と判断できます。
等分散である場合、平均値に差があるか否かの検定(ｔ検定)は、次のようにして行います。
t.test(male$身長, female$身長, var.equal=T) # Studentのt検定
ちなみに、等分散でない時は次のようにします。
t.test(male$身長, female$身長, var.equal=F) # Aspin-welchのt検定
以下に、ｆ検定とｔ検定を行うサンプルを掲げます。
−−−− ft03.rb ここから
#! ruby -Ks
# ｆ検定とｔ検定 (coding: Windows-31J)
require "rrxwin"
include Rrx
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", header=T)
male <- subset(dtf, 性別 == "男")
female <- subset(dtf, 性別 == "女")
ff <- var.test(male$身長, female$身長)
ff
flag <- (ff$p.value >= 0.05)
t.test(male$身長, female$身長, var.equal=flag)
EOS
# 行列を切り出してRプログラムを実行
cut_data("data03.csv") do |data|
temp_make(data) # 抽出した表を一時ファイルに書き出す
hs = rexec(rpro)
next unless hs[:sink]
print hs[:sink]
end
−−−− ft03.rb ここまで
上の実行結果をみると、男性の身長の平均値が 170.66、女性は 161.68 です。
ただ、ｔ検定のp値が 0.1887 です。これは 0.05 より大きいので、「男女の身長に有意な差があるとは言いきれない」ということになります。
!!(4) 身長と体重の検定を行う
data03.csvには、身長だけでなく体重のデータもあります。この両方についてｆ検定とｔ検定を行うことを考えます。
Rプログラムの中に「身長」と書き込んでしまうと、それを体重に応用することができません。そこで、「グローバル変数によるRプログラムの置換」のところで紹介したやり方を応用します。
まず、Rプログラムの中で「male$身長」と書いたところを「male$@@$item@@」と記述しておきます。
そして、rubyスクリプトの方で rexec() によってRプログラムを実行する直前に、グローバル変数 $item に "身長" とか "体重" を代入します。そうすると、Rプログラムの中の「@@$item@@」が「身長」とか「体重」に置き換えられた上で実行されます。
以下、そのサンプルです。
−−−− ft04.rb ここから
#! ruby -Ks
# 身長と体重の両方を検定 (coding: Windows-31J)
require "rrxwin"
include Rrx
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", header=T)
male <- subset(dtf, 性別 == "男")
female <- subset(dtf, 性別 == "女")
ff <- var.test(male$@@$item@@, female$@@$item@@)
ff
flag <- (ff$p.value >= 0.05)
t.test(male$@@$item@@, female$@@$item@@, var.equal=flag)
EOS
# 行列を切り出してRプログラムを実行
cut_data("data03.csv") do |data|
temp_make(data) # 抽出した表を一時ファイルに書き出す
%w(身長 体重).each do |item|
$item = item
hs = rexec(rpro)
next unless hs[:sink]
print hs[:sink]
end
end
−−−− ft04.rb ここまで
上のサンプルを実行してみると、身長に関するｔ検定の結果は、先述したように有意差が認められませんが、体重については、p値が 0.02153 と 0.05 より小さいので、男女の体重に有意な違いがあるということになります。
ちなみに、平均体重は、男性が 66.18、女性が 53.04 です。
ここで少し効率性を問題にしてみます。
上のスクリプトでは、Rプログラムを２度実行しています。したがって、当然ながら read.csv() も２度実行されます。
本来、データの読み込みは１回だけ行い、検定は２度行う(身長と体重について行う)というのが妥当なやり方です。その意味で、上のスクリプトは効率的でないといえます。
ということで、妥当と思われる方法の一例を下に掲げておきます。
−−−− ft05.rb ここから
#! ruby -Ks
# 身長と体重の両方を検定・その２ (coding: Windows-31J)
require "rrxwin"
include Rrx
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", header=T)
male <- subset(dtf, 性別 == "男")
female <- subset(dtf, 性別 == "女")
for (i in 1:ncol(dtf)) {
item <- colnames(dtf)[i]
if (item != "身長" && item != "体重") next
cat("**「", item, "」についての検定結果\n", sep="")
ff <- var.test(male[[item]], female[[item]])
print(ff)
flag <- (ff$p.value >= 0.05)
tt <- t.test(male[[item]], female[[item]], var.equal=flag)
print(tt)
}
EOS
# 行列を切り出してRプログラムを実行
cut_data("data03.csv") do |data|
temp_make(data) # 抽出した表を一時ファイルに書き出す
hs = rexec(rpro)
next unless hs[:sink]
print hs[:sink]
end
−−−− ft05.rb ここまで
ここには掲げませんが、上の ft05.rb を robj() 利用の形に書き換えた ft05b.rb を同梱しておきます。必要に応じて参考にして下さい。
以上でｆ検定とｔ検定を題材にした説明をおわりにします。
--------
!5. 相関係数と散布図の作成
ここでは、先でも利用した10人分の身長と体重のデータを使って、相関係数を算出します。また、身長と体重の散布図を描く方法について記します。
!!(1) 相関係数の算出
相関係数は、例えば、身長の高い人ほど体重も重い傾向がみられるか、といったことを検証する時に利用します。
0〜1 または -1〜0 の間の値を取り、相関係数の絶対値が0に近いほど二つの間の関連性が弱く、1に近いほど強いことになります。1に近ければ正比例のような関係、-1に近ければ逆比例のような関係といえます。
Rプログラムでは cor.test() で相関係数を求めることができます。
「xx <- cor.test(dtf$身長, dtf$体重)」とすれば、身長と体重の関係が得られます。
xx$estimate で相関係数を、xx$p.value でp値を参照できます。
cor.test() の結果を一通り出力するrubyスクリプトを掲げておきます。
−−−− cor01.rb ここから
#! ruby -Ks
# 相関係数の算出 (coding: Windows-31J)
require "rrxwin"
include Rrx
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", header=T)
xx <- cor.test(dtf$身長, dtf$体重)
robj(xx)
EOS
# 行列の切り取りとRプログラムの実行
cut_data("data03.csv") do |data|
temp_make(data) # 一時ファイルへの行列の書き出し
hs = rexec(rpro)
xx = hs["robj1"]
next unless xx
xx.each do |key, val|
printf("%s\t", key)
p val
end
end
−−−− cor01.rb ここまで
上のスクリプトを実行してみると、相関係数 xx["estimate"] が 0.8605302 と高い値です。
p値 xx["p.value"] は 0.001394111 となっており、0.05より小さいので、相関があると考えて差し支えなさそうです。
cor.test() には第３引数を与えることができます。相関係数の算出手法を指定するものです。次の３種類があります。
* ピアソンの積率相関係数 （pearson）
* ケンドールの順位相関係数 （kendall）
* スピアマンの順位相関係数 （spearman）
例えば次のように記述します。
cor.test(dtf$身長, dtf$体重, method="pearson")
cor.test(dtf$身長, dtf$体重, method="p") # 頭文字１文字のみ指定
第３引数を省略すると、「method="pearson"」が指定されたものとみなされます。
!!(2) 身長と体重の散布図の作成(pdfファイル)
data03.csvの身長と体重の関係を視覚的に把握するため、散布図を描いてみます。pdfファイルを作成します。
横軸を体重、縦軸を身長にして描きます。
Rプログラムでは plot() によってグラフや図を描くようです。以下にそのサンプルを掲げます。
−−−− cor02.rb ここから
#! ruby -Ks
# グラフをpdfに出力 (coding: Windows-31J)
require "rrxwin"
include Rrx
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", header=T)
pdf(file="cor02.pdf")
par(family = "Japan1GothicBBB") # 日本語フォントにゴチックを指定
par(mar = c(5, 5, 3, 3)) # 余白
plot(dtf$体重, dtf$身長, type="p",
main="身長と体重の散布図", cex.main=2, # メインタイトルとその大きさ
xlab="体重", ylab="身長", cex.lab=2)
dev.off()
EOS
# 行列の切り取りとRプログラムの実行
cut_data("data03.csv") do |data|
temp_make(data) # 一時ファイルへの行列の書き出し
hs = rexec(rpro)
print hs[:log]
end
−−−− cor02.rb ここまで
Rプログラムの pdf() は、pdfファイルの名前を指定して、それを出力先デバイスとして開くものです。
par() は、pdfファイルのパラメータを設定するもので、「family = ……」による日本語フォントの指定は欠かせません。これがないと、pdfに日本語が正しく書き込まれません。
plot() では様々なオプション指定ができるようですが、サンプルではメインタイトル、横軸と縦軸の項目名を指定しています。
dev.off() は、出力先デバイスとして開いたpdfファイルを閉じるものです。上のサンプルでは必要ないかもしれませんが、念のため置きました。
!!(3) MS-Windowsのメタファイル作成とExcelでの取り込み
MS-Windows専用の話題になりますが(linuxなどでは使えない)、Windowsの画像ファイルであるメタファイルを作成します。先のpdfファイルと同じ内容のものを作成します。
メタファイルは、拡張子が .emf のファイルで、Microsoft Office Picture Manager でそれを閲覧することができます。
このemfファイルは、Excel, Word, PowerPoint などに貼り付けることができます。下に掲げる cor03.rb では、Excelファイルに貼り付けています。
RプログラムでWindowsのメタファイルを作成するのは簡単で、pdfの時と同じように行います。パラメータの設定方法が少し違うだけです。
Excelを動かすのに、拙作 exlap.rb を使っています。
以下、サンプルです。
−−−− cor03.rb ここから
#! ruby -Ks
# グラフをemfに出力 (coding: Windows-31J)
require "exlap"
require "rrxwin"
include Rrx
# Rプログラムの定義
rpro = <<EOS
dtf <- read.csv("@@last@@", header=T)
win.metafile(file="cor03.emf")
windowsFonts(JP1=windowsFont("MS Gothic"), JP2=windowsFont("MS Mincho"))
par(family="JP1")
plot(dtf$体重, dtf$身長, type="p",
main="身長と体重の散布図", cex.main=2, # メインタイトルとその大きさ
xlab="体重", ylab="身長", cex.lab=2)
dev.off()
EOS
# 行列の切り取りとRプログラムの実行
cut_data("data03.csv") do |data|
temp_make(data) # 一時ファイルへの行列の書き出し
hs = rexec(rpro)
print hs[:log]
end
# Excelファイルへの画像の貼り付け
emf_file = Exl::fullpath("cor03.emf")
Exlap.new("cor03.xls") do |wb|
ss = wb.fes
ss.Range("A1:B1").Value = "統計解析ソフトRの画像取込み",
"(男女10人の身長と体重)"
rng = ss.Range("A3:J13") # 描画領域。大きさはどうでもいい
shape = ss.Shapes.AddPicture({
'Filename'=>emf_file,
'LinkToFile'=>false, # 素材画像にリンクをはらない
'SaveWithDocument'=>true, # 素材画像を取り込んで保存
'Left'=>rng.Left,
'Top'=>rng.Top,
'Width'=>0,
'Height'=>0})
shape.ScaleWidth(1.0, true) # 素材画像と同じ幅にする(1.0倍)
shape.ScaleHeight(1.0, true) # 素材画像と同じ高さにする(1.0倍)
wb.save
end
puts "'cor03.xls' を作成しました."
−−−− cor03.rb ここまで
「rng = ss.Range("A3:J13")」というのは、画像を貼り付ける領域を示すものですが、実は、その大きさは適当でかまいません。ここでは、一応、10行×10列の合計100個のセルを画像領域にしていますが、大事なのは A3 の始点です。
AddPicture() で画像の貼り付けを行っているわけですが、画像領域の左端を rng.Left、上端を rng.Top にしています。これは A3 の位置を指すものです。Width, Height については、とりあえず0にしておいて、後から素材画像と同じ幅・高さになるよう調整しています。
--------
!6. 覚え書き
覚え書きを補足的に記しておきます。
!!(1) erbの利用
rubyには、テキストファイルにrubyスクリプトを埋め込むためのerbというライブラリが標準添付されています。そのため、雛形となるテキストを場面に応じて柔軟に変更できます。Rプログラムを雛形テキストに見立てれば、このerbを活用できます。
例えば、同じRプログラムを何度も実行したいけれども、その都度 read.csv() に与える引数を切り替えたい(読み込むcsvファイルを次々と変更したい)といったケースがあります。
この場合、erbを用いると次のようなスクリプトになります。サンプルの chi02.rb を変更したものを掲げてみます。
−−−− chi02c.rb ここから
#! ruby -Ks
# erbの利用例 (coding: Windows-31J)
require "erb"
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv(<%= read_param %>)
xx <- chisq.test(dtf)
xx
EOS
##
read_param = "\"data01.csv\", header=T, row.names=1"
erb = ERB.new(rpro)
Rrx.rexec(erb.result(binding), "log.txt")
puts "'log.txt' を作成しました."
−−−− chi02c.rb ここまで
上記は、ごく簡単なRプログラムなので erb を用いるまでもありません。rrxwin.rb の「グローバル変数によるRプログラムの置換」でも十分に対応できます。
ただ、もっと複雑なRプログラムになると、erbの利用が効果を発揮するケースがあると思います。
!!(2) read.table(), write.table() の利用
これまでのサンプルでは、read.csv(), write.csv() を使ってきました。ですが、read.table(), write.table() の方がより汎用的です。
参考まで chi04.rb の変更版を掲げてみます。カイ２乗検定の調整済み残差を出力するものです。
「sep=","」によって、フィールドセパレータ(区切り文字)を指定している点が異なります。
タブ区切りテキストを扱いたい時は、「sep="\t"」とします。
また、write.table() のオプションに、「col.names=NA」を付け加えています。これがないと、第１行目の横軸の項目名の列びが、左端の空欄がない状態になり、第２行目以降に比べて、フィールドの数が１個不足することになります。
各々のデータを引用符で囲まなくてもいい場合は、write.table() のオプションとして「quote=F」を加えます。
−−−− chi04b.rb ここから
#! ruby -Ks
# read.table, write.table の利用 (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.table("data01.csv", sep=",", header=T, row.names=1)
dtf
xx <- chisq.test(dtf)
xx
write.table(xx$stdres, file="stdres.csv", sep=",", col.names=NA)
EOS
Rrx.rexec(rpro, "log.txt")
puts "'log.txt' と 'stdres.csv' を作成しました."
−−−− chi04b.rb ここまで
!!(3) 文字コードの調整
仮にRプログラムがeuc-jpで書かれているとします。その場合、Rにそのことを伝えるためには、起動オプションとして --encoding=euc-jp を指定します。例えば次のようにします。
Rcmd.exe BATCH --encoding=euc-jp test.r log.txt [enter]
linuxなどでは次のようになります。
R CMD BATCH --encoding=euc-jp test.r log.txt [enter]
また、念のため Rプログラムの先頭のところに次の１行を置いておきます。どのOSにも共通です。
options(encoding="euc-jp")
このように文字コードを指定しておくと、Rプログラムの中に日本語が交じっていても正しく処理されます。読み込むファイルなどの文字コードもeuc-jpが想定されることになります。
rrxwin.rbは、実行中のrubyスクリプトの文字コードに従って、上の文字コード調整を自動的に行います。
実行中のrubyスクリプトの文字コードは、ruby 1.8系であれば $KCODE をみて判断し、ruby 1.9系であれば、rubyスクリプトのencode指定によって判断します。
Rrx.rencode() というメソッドがあり、これは前述のようにして判断した文字コードに変換するものです。例えば次のように用います。
str2 = Rrx.rencode(str)
こうすると、rubyスクリプトの文字コードがeuc-jpであれば、str2がeuc-jpの文字列として生成されます。
どの文字コードで書かれているか不明なファイルを読み込むような場合、この rencode() を適当に適用してやれば、トラブル回避になるかもしれません。
rencode()での文字コードの変換は nkf を用いています。
ちなみに、cut_table(), cut_num(), cut_data() では内部的に rencode() を呼び出しているので、わざわざ呼び出す必要はありません。
--------
以上で rrxwin.rb についての説明をおわりにします。
rrxwin.rb で定義されている各種メソッドなどについては、同梱の[[rrxwin02.txt|http://cup.sakura.ne.jp/rrxwin02.htm]]を参照して下さい。
〜 以上 〜
Copyright (C) T. Yoshiizumi, 2013 All rights reserved.