統計解析ソフトR(以下、統計R)をpythonから利用する場合、
rpy2 とか pypeR というライブラリを利用できます。
ただ、マルチバイト文字の日本語を扱おうとすると苦労します。
そこで、rpy2 を利用しつつ、日本語を扱いやすくするため
pyrcmd.py を作りました。
その pyrcmd の使い方について記します。
今回は基本編です。
統計処理した結果を html や docx(Wordファイル)にするのは
次回にしたいとおもいます。
なお、当サイトで紹介したスクリプトは
pyrcmd02.zip に入れてあります。
スクリプトを動かすのに必要な環境については
pythonから統計解析ソフトRを利用する方法・第１回を参照してください。
pythonライブラリの rpy2 を使えば、統計Rのスクリプトを実行できます。
ただ、スクリプトに間違いがあって正常に実行できない場合、
rpy2が出力するエラー報告を見て、
統計Rのスクリプトの間違いをみつけるのは苦労です。
ちなみに、pythonライブラリの pypeR だと、エラー報告がなく、
pythonスクリプトが終了しないままになります。
やはり、統計Rのログを残して後からチェックできれば便利です。
pythonスクリプトに import pyrcmd as prc という１行を書くと、
prc.rexec(……) により統計Rのスクリプトを実行できます。
引数としてスクリプトファイルの名前を渡してもいいし、
スクリプトを文字列として渡すこともできます。
たとえば下のようにします。
-------- ここから
# encoding: cp932
import pyrcmd as prc
rscript = u'''\
x <- c("大阪", "名古屋", "東京")
x
'''
prc.rexec(rscript)
-------- ここまで
上の pythonスクリプトを実行すると、カレントディレクトリに
pyrcmd_log01.txt というファイルが書き出されます。
その中身は下のとおり。
> options(encoding='cp932')
> sink(file='temp01.dat', split=T)
> x <- c("大阪", "名古屋", "東京")
> x
[1] "大阪" "名古屋" "東京"
>
> save(list = ls(), file="temp02.dat")
>
> proc.time()
ユーザ システム 経過
0.32 0.07 0.39
もとの統計Rのスクリプトには書かれていなかった
options(……)、sink(……)、save(……) が付け加えられています。
それについては後述するとして、
統計Rの実行過程と結果がちゃんと記録されています。
統計Rの警告メッセージ、エラーメッセージも記録されるので
間違いの修正に便利です。
なお、同じ pythonスクリプト内で rexec() を複数回呼び出すと
２回目のログファイルは pyrcmd_log02.txt という名前になります。
ログファイルを残したくない場合は、
log=False というオプションを指定します。
prc.rexec('test.R', log=False)
上のようにすると、ログファイルが残りません。
pyrcmd がどうやって統計Rのスクリプトを実行するかを記します。
Rコマンド(Windowsなら R.exe)にはバッチモードで起動する方法があります。
test.R というスクリプトを実行するなら次のとおり。
R CMD BATCH --no-save -q --encoding=cp932 test.R log.txt [enter]
上のコマンドを実行すれば log.txt にログが記録されます。
(オプションは、あくまで単なる例です。)
pyrcmd では、上記のようなコマンドを子プロセスとして実行しています。
統計Rには reticulate というパッケージがあります。
統計Rから python を利用するためのパッケージです。
それをインストールするための pythonスクリプトを掲げておきます。
-------- ここから
import pyrcmd as prc
url = "https://cran.r-project.org/"
packageName = "reticulate"
rscript = '''\
options(CRAN="%s")
options(repos="%s")
install.packages("%s", dependencies = TRUE)
library(%s)
search()
''' % (url, url, packageName, packageName)
prc.rexec(rscript)
-------- ここまで
Windowsでは proxy が設定されている環境下でも実行可能だとおもいますが、
セキュリティが厳しく設定されている場合は、もしかするとダメかもしれません。
ログを見れば、正常にインストールされたかどうかチェックできます。
統計Rのログには、どんなステートメントが実行されたかを含め
基本的にすべての記録が残されます。
それとは別に、変数の値や検定結果を出力したとき、
その出力だけを得たいことがあります。
x <- 123
print(x)
上の統計Rスクリプトを実行すると
[1] 123
というのが出力されますが、この出力をデータとして得たい場合です。
prc.rexec(……) は、戻り値としてそうした出力を返します。
sink_data = prc.rexec('test.R')
上のようにすると、test.R での出力が変数 sink_data に代入されます。
その仕組みは次のとおりです。
統計Rのスクリプト内に下の１行を挿入します。
sink(file='temp01.dat', split=T)
すると、出力が temp01.dat というファイルに書き出されます。
それを読み取ってデータとして返す訳です。
temp01.dat というファイルは、読み取った後に消去します。
temp01.dat という名前は例えとして出したもので、
実際は python の tempfileライブラリを使います。
このような出力を得る必要がないときは、
prc.rexec() のオプション指定として sink=False を与えます。
sink_data = prc.rexec('test.R', sink=False)
こうすると、変数 sink_data には None が代入されます。
統計Rのスクリプト内に sink(……) という１行が挿入されることもありません。
pyrcmd は、文字エンコーディングとして次ぎのことを仮定します。
別のエンコーディング(euc-jpなど)にしたいなら、pythonスクリプト内で
prc.set_charset('euc-jp')
上の１行を記述します。
そうすると、pythonスクリプトが euc-jp で書かれていると想定します。
統計Rのスクリプトが pythonスクリプト内に書かれていれば、
それも euc-jp だと仮定します。
統計Rのスクリプトが test.R のように別ファイルになっているときは、
そのファイルの文字エンコーディングを自動判別してそれを採用します。
少々ややこしいですが、
要は、Windows環境下であれば cp932 を使い、
それ以外の環境では utf-8 を使うのが最も無難です。
文字エンコーディングの話をもう少し。
統計Rは、低レベルのところで OSのlocaleを参照しているとおもいます。
Windows(日本語版)だと cp932 です。
(マイクロソフトのWebを見ると、残念なことに utf-8 に変更できないようです。)
そして、統計Rのスクリプトがどのエンコーディングで書かれていても
低レベルのエンコーディングには影響しない、と推測します。
pyrcmd は、この低レベルのエンコーディングを utf-8 と仮定しますが、
Windows環境下であれば cp932 を仮定します。
もし OSのlocaleが別の設定になっている場合は、
prc.set_charset('euc-jp', 'euc-jp')
上のようにして第２引数を与えます。
そうすると、euc-jp が低レベルのエンコーディングであると仮定されます。
pyrcmd は、統計Rのスクリプトを加工した上で実行しますが、
文字エンコーディングに関していうと次のようになります。
統計Rスクリプトの文字エンコーディングが cp932 と判別された場合でいうと、
R CMD BATCH …… のオプションとして --encoding=cp932 を指定。options(encoding="cp932") を挿入。
pythonスクリプト内で pyrcmd.py を import すると、
python ver 2.x の場合ですが、デフォルトの文字エンコーディングが設定されます。
Windows の場合は cp932、それ以外なら utf-8 です。
また、set_charset() が呼び出されると、
指定された文字エンコーディングがデフォルトのエンコーディングになります。
prc.set_charset('euc-jp')
上のようにすると、pythonのデフォルトのエンコーディングが euc-jp になります。
これは python ver 2.x の場合の話で、python ver 3.x では関係ありません。
おせっかいな機能とはおもいましたが、
python2 でデフォルトのエンコーディングを設定するのが面倒に感じられ、
こんな仕様にしてみました。
prc._charset = 'euc-jp'
prc._default_charset = 'euc-jp'
上のように記述すれば、
デフォルトのエンコーディングが操作されることなく
pyrcmd が仮定するエンコーディングが設定されます。
pyrcmdで統計Rのスクリプトを実行した場合も
当然ながら python側で各種変数の値を参照できます。
rpy2を利用することでそれを実現しています。
prc.rexec() を呼び出して統計Rのスクリプトを実行した後は、
x = r['x'] などとして統計R側の変数をpython側で取得できます。
たとえば下のとおり。
-------- ここから
# encoding: cp932
from rpy2.robjects import r
import pyrcmd as prc
rscript = u'''\
x <- c("大阪", "名古屋", "東京")
x
'''
prc.rexec(rscript)
x = r['x']
print(x)
-------- ここまで
python側の変数 x は rpy2のオブジェクトなので、下のように表示されます。
[1] "大阪" "名古屋" "東京"
次の項で、背景にある仕組みについて述べます。
統計Rには、各種変数をファイルに書き出すための save() と
そのファイルを読み込んで変数を復元するための load() があります。
変数 x をファイルに書き出して、それを復元する例を掲げます。
-------- ここから
# encoding: cp932
from rpy2.robjects import r
import pyrcmd as prc
rscript = u'''\
x <- c("大阪", "名古屋", "東京")
save(x, file="data01.dat")
x <- NULL
print(x)
load("data01.dat")
print(x)
'''
sink_data = prc.rexec(rscript, rdata=False)
print(sink_data)
-------- ここまで
変数 x に文字列からなるベクトルを代入し、それをファイルに書き出します。
その後、変数 x に NULL を代入して print文で表示し、
次に load() で変数 x を復元・表示しています。
pythonスクリプトの出力は下のとおり。
NULL
[1] "大阪" "名古屋" "東京"
上記は一つの変数だけをファイルに書き出す例ですが、
save(list = ls(), file = "data01.dat")
上のようにすると、すべての変数や関数をファイルに書き出すことになります。
pyrcmd は、この save, load を利用しています。
rexec() で統計Rのスクリプトを実行する際、
暗黙のうちに save() を挿入・実行し、
Rのスクリプトの実行が終了したら、
rpy2を使って load() を実行します。
そうすれば、統計Rのスクリプトに出てきた全変数が
rpy2 のオブジェクトとして使える状態になります。
save, load の利用を抑制したいときは、
rexec() に rdata=False というオプションを与えます。
prc.rexec('test.R', rdata=False)
前述の pythonスクリプトでは、実は、このオプションを指定しました。
巨大なデータフレームが一つの変数に入っていて、
それを含め save, load するのは大変だというようなケースで、
自前で save, load の処理を行う場合は rdata=False を指定して下さい。
自前で処理する例を掲げておきます。
三つの変数 x, y, z のうち、x, y の二つを save する例です。
-------- ここから
# encoding: cp932
import os
from rpy2.robjects import r
import pyrcmd as prc
rscript = u'''\
x <- c("大阪", "名古屋", "東京")
y <- 456
z <- c("cat", "dog", "fox")
save(list=c("x", "y"), file="data01.dat")
'''
prc.rexec(rscript, rdata=False)
r('load("data01.dat")')
os.remove("data01.dat")
x = r['x']
y = r['y']
try:
z = r['z']
except:
z = None
print(x)
print(y)
print(z)
-------- ここまで
変数 z は復元の対象から外したので
pythonスクリプトの側で z = r['z'] とするとエラーが発生します。
上記では try, except を使って、
エラー発生によるスクリプトの実行中断を回避しています。
下の１行を pythonスクリプト内に書くと
r というオブジェクトを使えるようになります。
from rpy2.robjects import r
一方、pyrcmd.py の中にも同じ１行が書かれています。
あるpythonスクリプトが pyrcmd.py を import した場合、
r というオブジェクトが pyrcmd.py との間で共有されるのかどうか。
試してみると共有されるようです。
たとえば a.py の中で r('x <- 123') と書くと
b.py の中で x = r['x'] のように xを参照することができます。
もちろん、その場合、b.py の中で a.py をimportする必要があります。
このように共有される rオブジェクトには、
統計Rの変数とか関数がどんなふうに蓄積されるのか気になります。
一つの pythonスクリプト内で prc.rexec() を複数回 呼び出したとき、
１回ごとに記録がリセットされるのか、それとも継承されるのか……
結論からいうと、継承されます。
意図的に記録を消去しないかぎり、
rオブジェクトにどんどんデータが蓄積されるようです。
そこで、prc.r_clear() とすれば記録が消去されるようにしました。
r('rm(list = ls())')
r('gc()')
上の２行で、統計Rに関連する変数や関数の記録を消去できます。
ただ、これだけだとメモリーが整理された形で解放されないようで、
結局、次の４行で記録を消去します。
(import gc が必要)
gc.collect()
r('rm(list = ls())')
r('gc()')
gc.collect()
prc.r_clear() の中身は上の４行です。
これまで rexec() のオプションをばらばらに説明してきました。
まだ触れていないオプションもあるので、それを含めて取りまとめてみます。
また、pyrcmd.py の中で定義している関数について記します。
rexec(rsf, options=None, log=True, sink=True, rdata=True)
統計Rのスクリプトを実行し、sink_dataを返す。
R CMD BAT …… に与えるオプション。--no-save -q --internet2 が使われる。
引数はない。
統計Rに関する変数・関数をすべて消去。
set_charset(charset = None, default_charset = None)
仮定すべき文字エンコーディングを設定。
pyrcmdモジュール内では、
それぞれ _charset, _default_charset というグローバル変数に値を代入する。
python2の場合、_charset にセットされた文字エンコーディングを
デフォルトのエンコーディングとして設定する。
get_charset(fname)
引数として与えられたファイル名のファイルの文字エンコーディングを返す。
戻り値は文字列。
pythonライブラリの chardet を利用。
read_text(filename, charset='')
テキストファイルの読み込み、unicode文字列を返す。
write_text(filename, text, charset='')
テキスト文字列をファイルに書き出す。
引数はない。
pyrcmd_log01.txt, pyrcmd_log02.txt などログファイルをすべて削除。
print_dict(dct)
引数として渡された辞書(dict, OrderedDict)の中身をprint文で表示。
key, data type, value を空白行を区切りとして表示する。
以上です。
他に下請け的な処理を行う関数がありますが省力します。
rp2pd() などについては、この先で触れます。
rpy2のオブジェクトは、統計Rと類似の処理方法に向いていますが、
python側で扱いやすいかというと、かならずしもそうとはいえません。
そこで、rpy2のオブジェクトをpandasの Series や DataFrame に変換するための
rp2pd() を設けました。
単純な数値や文字列に変換可能なオブジェクトは、
pythonの数値や文字列に変換し、
そうでないオブジェクトを Series, DataFrame に変換します。
なお、rpy2のオブジェクトでないものを rp2pd() に渡すと
None を返すのでご注意ください。
圧縮ファイル pyrcmd02.zip に入っているのは、この先で紹介するスクリプトです。
いずれも data.csv を読み込んで処理します。
「ＩＤ、性別、身長、体重」の４列×400行のデータで
ＩＤ以外の列には、若干の欠損値が含まれています。
ＩＤ,性別,身長,体重
C3,女性,159.1,57.8
W5,男性,163.8,78.2
W11,女性,162.7,59.5
…………
統計Rの summary() は、
数値データを与えると最小値、平均値、最大値などを返します。
数値データ以外だと、該当の度数(人数・個数)を返します。
性別と身長の summary を pandas オブジェクトに変換して表示してみます。
まず、結果として得られる表示は下のとおり。
「性別」のsummary(該当人数)は下のとおり.
NA 2
女性 194
男性 204
dtype: int64
身長のsummaryは下のとおり.
Min. 144.80000
1st Qu. 158.67500
Median 163.60000
Mean 163.46199
3rd Qu. 168.02500
Max. 185.20000
NA's 8.00000
dtype: float64
性別のところに出てくる NA は、性別の記載のない人を示します。
身長の NA's は、身長が欠損している人の数を示します。
上の表示を得るための pythonスクリプトは次のとおり。
1# prc01s.py (encoding: cp932) 2import pandas as pd 3from rpy2.robjects import r 4from textwrap import dedent 5import pyrcmd as prc 6prc.set_charset('cp932') 7 8src_file = 'data.csv' 9rscript = dedent(u'''\ 10 dtf <- read.csv("%s", header=T, fileEncoding="%s") 11 gender_summary <- summary(dtf$性別) 12 height_summary <- summary(dtf$身長) 13 ''') % (src_file, prc.get_charset(src_file)) 14prc.rexec(rscript, sink=False) # Rコマンドの実行 15 16gender_summary = r['gender_summary'] # rp2オブジェクトとして取得 17ser = prc.rp2pd(gender_summary) # pandasオブジェクトに変換 18print(u"「性別」のsummary(該当人数)は下のとおり.") 19print(ser) 20print('') 21 22height_summary = r['height_summary'] 23ser = prc.rp2pd(height_summary) 24print(u"身長のsummaryは下のとおり.") 25print(ser)
python側の変数 gender_summary は rpy2オブジェクトの IntVector です。
height_summary は FloatVector という型です。
どちらも１次元のベクトルなので
prc.rp2pd() で変換すると pandas.Series になります。
もし２次元の Matrix とか DataFrame であれば
pandas.DataFrame に変換されます。
gender_summary, height_summary の両方とも「名前付きベクトル」ですが、
Series にその名前がちゃんと引き継がれます。
性別のところに出てくる NA は、統計R, rpy2 では "" の空文字列です。
でも、"" のままだと処理しにくい点がいろいろ出てくるので
rp2pd() で変換したときは "NA" という２文字に変換します。
統計Rのスクリプトに関連して少し補足します。
スクリプト全体をpythonの関数 dedent() に渡していますが、
これは indent の反対で、実行時に行頭のインデント空白を除去します。
統計Rのスクリプトを字下げして書けるので、見やすくなるとおもいました。
dtf <- read.csv("%s", header=T, fileEncoding="%s")
上の %s は、今回、次のように置き換えられます。
dtf <- read.csv("data.csv", header=T, fileEncoding="utf-8")
utf-8 は、prc.get_charset() で判別した値です。
data.csv の文字エンコーディングを別のものに変更したとしても、
上記スクリプトの実行には支障ないとおもいます。
統計Rには table() という集計用の関数があります。
引数を一つだけ与えたときは１次元の Array を返し、
二つ与えると２次元の Array を返します(クロス集計)。
２次元の Array は、rpy2オブジェクトとして受け取ると Matrix型になるようです。
１次元、２次元の Array を受け取って pandasオブジェクトに変換してみます。
まず、数値データの身長をカテゴリカルデータとして扱えるようにします。
155未満、155〜165、165〜175、175以上のカテゴリに区分けします。
その上で、その「身長区分」だけを table() に渡して集計します。
また、「身長区分」と「性別」を渡したクロス集計も行います。
まず、pythonスクリプトを実行した結果の表示を掲げます。
＊身長区分の集計結果
155未満 46
155〜165 182
165〜175 142
175以上 22
不明 8
dtype: int64
＊身長区分×性別の集計結果
男性 女性 記載なし
155未満 6 40 0
155〜165 72 109 1
165〜175 100 41 1
175以上 20 2 0
不明 6 2 0
次に pythonスクリプトは下のとおり。
クロス集計表は Excelファイルとしても書き出します。
1# prc02s.py (encoding: cp932) 2import pandas as pd 3from rpy2.robjects import r 4from textwrap import dedent 5import xlwt 6import pyrcmd as prc 7prc.set_charset('cp932') 8 9src_file = 'data.csv' 10rscript = dedent(u'''\ 11 dtf <- read.csv("%s", header=T, na.string="", fileEncoding="%s") 12 dtf$性別 <- factor(dtf$性別, levels=c("男性", "女性", NA), 13 labels=c("男性", "女性", "記載なし"), exclude=NULL) 14 clevels <- c("155未満", "155〜165", "165〜175", "175以上") 15 cc <- cut(dtf$身長, breaks=c(-Inf, 155, 165, 175, Inf), 16 labels=clevels, right=FALSE) 17 cc <- factor(cc, levels=c(clevels, NA), 18 labels=c(clevels, "不明"), exclude=NULL) 19 dtf <- cbind(dtf, 身長区分=cc) 20 tbl1 <- table(dtf$身長区分) 21 tbl2 <- table(dtf$身長区分, dtf$性別) 22 ''') % (src_file, prc.get_charset(src_file)) 23prc.rexec(rscript, sink=False) # Rコマンドの実行 24tbl1 = r['tbl1'] # rp2オブジェクトとして取得 25tbl2 = r['tbl2'] 26ptbl1 = prc.rp2pd(tbl1) # pandasオブジェクトに変換 27ptbl2 = prc.rp2pd(tbl2) 28 29print(u'＊身長区分の集計結果') 30print(ptbl1) 31print('') 32print(u'＊身長区分×性別の集計結果') 33print(ptbl2) 34 # ptbl2をExcelに出力 35writer = pd.ExcelWriter("prc02s.xls") 36ptbl2.to_excel(writer, startrow=1, startcol=0) 37worksheet = writer.sheets['Sheet1'] 38worksheet.write(0, 0, # 引数は row, column の順番 39 u"身長区分×性別の集計表") 40writer.save()
統計Rスクリプトの詳細は省略しますが、
pyrcmd.py における Array の扱いについて少々。
統計Rの Array は、３次元データ、４次元データなども表現できます。
ですが、pyrcmd.py では２次元までしか扱えません。
rpy2が３次元以上にどのように対応しているのか確認してませんが、
統計Rスクリプトの側で、２次元データに落とし込むなどして
対応するのが無難だとおもいます。
３次元データは２次元データの集合体、
４次元データは３次元データの集合体です。
1〜9の数字が電話の数字ボタンと同じ３行×３列に並ぶ matrix を考えます。
m <- t(matrix(1:9, nrow = 3))
上の統計Rスクリプトで目的の matrix を生成できます。
これをpython側の変数 m に代入し、
さらに m2 = prc.rp2pd(m) と変換した場合、
m2を print文で表示すると下のようになります。
V1 V2 V3
0 1 2 3
1 4 5 6
2 7 8 9
m2は pandas.DataFrame型で、行ラベルは 0, 1, 2 の通し番号になります。
列ラベルの方は、V1, V2, V3 という名前になります。
行ラベルや列ラベルのないオブジェクトを rp2pd() で変換すると、
行ラベルには 0 から始まる通し番号、
列ラベルには V1, V2, V3, …… が割り当てられます(Vは大文字)。
身長の平均値が男女間で異なるかどうか検証するため、t検定をおこないます。
検定の結果を python側で受け取って表示します。
検定結果は、rpy2のオブジェクトでいうと ListVector型です。
それを rp2pd() で変換すると、OrderedDict型のデータになります。
それを表示すると下のとおり。
前半は統計Rの表示そのまま、後半は各項目の詳細表示です。
＊t検定の結果に関する統計Rの表示
Welch Two Sample t-test
data: male$身長 and female$身長
t = 10.374, df = 386.78, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
5.362199 7.870155
sample estimates:
mean of x mean of y
166.7157 160.0995
＊t検定の結果に関する詳細表示
statistic
(Series)
t 10.373527
dtype: float64
parameter
(Series)
df 386.779601
dtype: float64
p.value
(float64)
2.08272110497e-22
conf.int
(Series)
0 5.362199
1 7.870155
dtype: float64
estimate
(Series)
mean of x 166.715657
mean of y 160.099479
dtype: float64
null.value
(Series)
difference in means 0.0
dtype: float64
alternative
(unicode)
two.sided
method
(unicode)
Welch Two Sample t-test
data.name
(unicode)
male$身長 and female$身長
データ型として Series と unicode文字列が多いですが、
p.value は数値(float型)です。
p.value が 0.05 とか 0.01 より小さいので男女間に有意差があるといえます。
今回は両側検定の結果なので「差がある」という判断ですが、
大きいか否か、小さいか否かを見るには片側検定の結果を参照します。
下に、pythonスクリプトを掲げます。
1# prc03s.py (encoding: cp932) 2import pandas as pd 3from rpy2.robjects import r 4from textwrap import dedent 5import pyrcmd as prc 6prc.set_charset('cp932') 7 8src_file = 'data.csv' 9rscript = dedent(u'''\ 10 dtf <- read.csv("%s", header=T, fileEncoding="%s") 11 male = subset(dtf, 性別 == "男性") 12 female = subset(dtf, 性別 == "女性") 13 flag <- FALSE # 不等分散を仮定 14 t_test <- t.test(male$身長, female$身長, var.equal=flag) 15 t_test 16 ''') % (src_file, prc.get_charset(src_file)) 17sink_data = prc.rexec(rscript) # Rコマンドの実行 18t_test = r['t_test'] # rp2オブジェクト(ListVector)として取得 19pod = prc.rp2pd(t_test) # ListVectorを変換(OrderedDictになる) 20 21print(u'＊t検定の結果に関する統計Rの表示') 22print(sink_data) 23print('') 24print(u'＊t検定の結果に関する詳細表示') 25prc.print_dict(pod)
prc.print_dict() は、辞書型のデータを
key, data type, value を単位として表示する関数です。
だらだらと長い表示になってしまいますが、
中身を取りこぼしなく一通りチェックできます。
今度はカイ２乗検定の結果を受け取って表示します。
検定結果には DataFrame が含まれるので、print文で出力するほか
Excelファイルにも書き出してみます。
まず、身長を155未満、155〜165、165〜175、175以上のカテゴリに区分けします。
その上で、「身長区分」×「性別」のクロス集計を行い、
その集計表をカイ２乗検定にかけます。
今回は欠損値を対象外とします。
スクリプトを実行して得られる表示は下のとおり。
＊カイ２乗検定の結果
chi2: 72.033928, df: 3, p-value: 1.56548e-15
＊観測度数
男性 女性
155未満 6 40
155〜165 72 109
165〜175 100 41
175以上 20 2
＊調整済み残差
男性 女性
155未満 -5.449433 5.449433
155〜165 -4.040047 4.040047
165〜175 5.990430 -5.990430
175以上 3.876832 -3.876832
上と同じ内容が Excelファイルにも書き出されます。
今回、有意確率の値が 0.05 より小さい(というより 0に近い)ので
身長区分と性別の間に有意な関連性があると判断して差し支えないことになります。
となると、気になるのが調整済み残差です。
調整済み残差の表をみれば、表の中のどのセルが
有意に大きいか、あるいは小さいかを見ることができます。
有意水準を５％とした場合、
各セルの値が 1.96 より大きければ有意に大きいと判断でき、
-1.96 より小さければ有意に小さいと判断できます。
有意水準が１％なら 2.58, -2.58 が目安の値です。
今回は、どのセルも有意水準１％で有意性が検出されます。
次に、pythonスクリプトです。
1# prc04s.py (encoding: cp932) 2import pandas as pd 3from rpy2.robjects import r 4from textwrap import dedent 5import xlwt 6import pyrcmd as prc 7prc.set_charset('cp932') 8 9src_file = 'data.csv' 10rscript = dedent(u'''\ 11 dtf <- read.csv("%s", header=T, fileEncoding="%s") 12 dtf$性別 <- factor(dtf$性別, levels=c("男性", "女性"), 13 labels=c("男性", "女性")) 14 cc <- cut(dtf$身長, breaks=c(-Inf, 155, 165, 175, Inf), 15 labels=c("155未満", "155〜165", "165〜175", "175以上"), right=FALSE) 16 dtf <- cbind(dtf, 身長区分=cc) 17 tbl <- table(dtf$身長区分, dtf$性別) 18 chi_test <- chisq.test(tbl) 19 chi_test 20 ''') % (src_file, prc.get_charset(src_file)) 21prc.rexec(rscript, sink=False) # Rコマンドの実行 22chi_test = r['chi_test'] # rp2オブジェクト(ListVector)として取得 23pod = prc.rp2pd(chi_test) # ListVectorを変換(OrderedDictになる) 24chi2 = pod['statistic'][0] # カイ２乗値 25df = pod['parameter'][0] # 自由度 26p_value = pod['p.value'] # 有意確率の値 27 28print(u'＊カイ２乗検定の結果') 29print("chi2: %f, df: %d, p-value: %g\n" % (chi2, df, p_value)) 30print(u"＊観測度数") 31print(pod['observed']) 32print('') 33print(u"＊調整済み残差") 34print(pod['stdres']) 35 36writer = pd.ExcelWriter("prc04s.xls") 37pod['observed'].to_excel(writer, startrow=4, startcol=0) # 観測度数 38worksheet = writer.sheets['Sheet1'] 39worksheet.write(0, 0, u"＊カイ２乗検定の結果") 40worksheet.write(1, 0, "chi2: %f, df: %d, p-value: %g" % (chi2, df, p_value)) 41worksheet.write(3, 0, u"＊観測度数") 42nr = worksheet.last_used_row + 2 # 最下行の次の次 43worksheet.write(nr, 0, u"＊調整済み残差") 44nr = nr + 1 45pod['stdres'].to_excel(writer, startrow=nr, startcol=0) 46writer.save()
t.test() のところを下のようにすれば片側検定となり、
「男性＞女性」と判断していいかどうかの判断材料を得ることができます。
t.test(male$身長, female$身長, var.equal=F, alternative="greater")
Excelファイルの書き出し部分は、順序がちょっと不自然ですが、
writer = pd.ExcelWriter(……) としている関係で、
とりあえず DataFrame.to_excel() を実行しないと
ワークシートオブジェクトを取得できないようです。
そのため、１行目を書き出す前に、まず観測度数の表を書き出しています。
今回は、この辺でおわりにします。
ここでは cp932 のスクリプトを掲げましたが、
zip圧縮ファイルには utf-8 のスクリプトも入っています。
〜 以上 〜
Copyright (C) T. Yoshiizumi, 2018 All rights reserved.