T. Yoshiizumi - tjs_d_guide Diff
- Added parts are displayed like this.
- Deleted parts are displayed
like this.
単純集計に関する覚え書き／まとめ(関数化)
〜 rubyとRを用いた処理 〜
最終更新日: 2013/08/11
以下に掲げるドキュメントは、
[[tjs_d.zip|http://cup.sakura.ne.jp/tjs_d.zip]]
に同梱されている tjs_d.txt と同じ内容です。{{br}}
tjs_d.zipには、euc-jp | utf-8 で書かれたサンプルスクリプトも同梱されています。{{br}}
サンプルスクリプトは、MS-Windowsだけでなく、linuxなどでも動作します(CENT6(CentOS6で確認)。
＜はじめに＞
「単純集計に関する覚え書き」を３回にわたって記しました。予想していたよりも長い解説になってしまい、読み返すのが面倒です。
具体的な集計の必要に迫られたとき、どのサンプルスクリプトを利用したらいいのか、どんなふうに書き換えたらいいのか、私自身が迷うくらいにだらだらとした解説になってしまいました。
そこで、「最低限これくらいは」というものをRの関数として定義しました。例えば、単一回答の集計の一つ(無回答を除く集計)の sa1() は次のように呼び出します(saは single answer のつもり)。
xx <- sa1(dtf, 2, 3)
print(xx)
こうすると、データフレームdtfの２列目と３列目に着目して、クロス集計表を作成します。
戻り値のxxには４つの表が入っています。２つの度数分布表と２つのパーセントの表です。
こうした sa1() のような関数をいくつか作って tjsf.r に入れました。以下で、それらの解説を記します。
なお、この tjsf.r は、rubyに依存する部分はありません。rubyを使わない場合でも利用可能です。
--------
{{toc_here}}
--------
!1. 単一回答の集計
ここでは、「性別」(男・女)、「意見」(賛成・中立・反対)のように、選択肢の中から１つだけ選んでもらう単一回答の集計を取り上げます。
適切な選択肢がないため回答できない｜回答したくないというケースもありますが、それは無回答という形になって表れると思います。そうした無回答を含めて集計する場合についても取り上げます。
単一回答にまつわる主な関数は sa1(), sa2(), sa3() の３つです。
!!(1) 単一回答の集計に関する関数 sa1
ここで扱い素材データは、ID(整理番号)、性別(男・女)、意見(賛成・中立・反対)の３列からなるcsvデータ(data01.csv)です。98人分が記録されています。具体的には次のような内容です。
ID,性別,意見
1,女,賛成
2,男,賛成
3,男,賛成
IDには空欄がありませんが、性別と意見の列には少数ながら空欄があります。つまり、欠損値(無回答)があります。
この素材データを集計します。
まず、欠損値を数に入れずに集計する sa1() のサンプルを示します。Rプログラム部分のみです。
−−−− sa1()の利用例ここから
source("tjsf.r")
dtf <- read.csv("data01.csv", header=T, na.strings="")
dtf$性別 <- factor(dtf$性別, levels=c("男", "女"), labels=c("男", "女"))
xx <- sa1(dtf, "性別", "意見")
xx
−−−− sa1()の利用例ここまで
「dtf$性別 <- factor(……)」は、なくても支障はありませんが、これがないと、作成される表の中の「男・女」の順番が「女・男」になります。
sa1() が返す値はリストで、次の４つのテーブルを含みます。
・tbl1: 総数(計)なしの表
・tbl2: 総数(計)を付加した表
・pct1: 横軸に足し算を行った「計」に対するパーセントの表
・pct2: 縦軸に足し算を行った「計」に対するパーセントの表
Rプログラムの出力結果は次のとおり。
−−−−
$tbl1
意見
性別 賛成 中立 反対
男 14 12 11
女 9 17 20
$tbl2
意見
性別 賛成 中立 反対 総数
男 14 12 11 37
女 9 17 20 46
総数 23 29 31 83
$pct1
意見
性別 賛成 中立 反対 総数
男 37.83784 32.43243 29.72973 100.00000
女 19.56522 36.95652 43.47826 100.00000
総数 27.71084 34.93976 37.34940 100.00000
$pct2
意見
性別 賛成 中立 反対 総数
男 60.86957 41.37931 35.48387 44.57831
女 39.13043 58.62069 64.51613 55.42169
総数 100.00000 100.00000 100.00000 100.00000
−−−−
sa1() の第２引数("性別")と第３引数("意見")は、列名でなく列の番号で指定してもかまいません。次のように書いても同じ結果になります。
xx <- sa1(dtf, 2, 3)
tbl1だけを取り出したい時は、xx$tbl1 などのようにします。これには「計」がないので、カイ２乗検定等にかけることができます。例えば次のように書きます。
xx <- sa1(dtf, 2, 3)
zz <- chisq.test(xx$tbl1)
こうすると、zzにカイ２乗検定の結果が代入されます。
tbl1, tbl2 などの各々のテーブルには、縦軸のタイトル「性別」と横軸のタイトル「意見」が付いています。これは、data01.csvの列名をそのまま割り当てたものです。
これらタイトル(正式には「次元名」というんでしょうか)を得るには、次のようにします。
dnn <- names(dimnames(xx$tbl1))
こうすると、dnnに「c("性別", "意見")」が代入されます。
sa1() の引数指定で、第３引数がなければ、１つの項目のみに着目した表が作成されます。
xx <- sa1(dtf, 2)
とした場合は、「性別」だけに着目した次のような表になります。tbl1, tbl2 の２つを示します。
−−−−
$tbl1
男 女
41 50
$tbl2
男 女 総数
41 50 91
−−−−
上のような１つの項目だけに着目した表を作った場合、横方向には複数の列ができますが、縦方向は１行のみです。
したがって、縦軸に沿ったパーセント(pct2)は意味がないのでNULLになります。
!!(2) 単一回答の集計に関する関数 sa2, sa3
先にsa1() について述べましたが、sa2(), sa3() も同じように用います。
呼び出す時の引数は２つか３つで、その戻り値は、tbl1など４つのテーブルからなるリストです。
これら sa1, sa2, sa3 の違いを下に列記します。
* sa1: 無回答を含まない表の作成
* sa2: 無回答を含む表の作成
* sa3: 有効回答(一種の小計になっている)を含む表の作成
以下に、sa1〜sa3の各々の出力例を示します。合計欄の「総数」が付いた表を掲げます。
欠損値(無回答)は、「無効回答」と表記されます。
−−−−
＊ sa1の出力例
|| ||賛成||中立||反対||総数
||男||14||12||11||37
||女||9||17||20||46
||総数||23||29||31||83
＊ sa2の出力例
|| ||賛成||中立||反対||無効回答||総数
||男||14||12||11||4||41
||女||9||17||20||4||50
||無効回答||2||3||2||0||7
||総数||25||32||33||8||98
＊ sa3の出力例
|| ||賛成||中立||反対||有効回答||無効回答||総数
||男||14||12||11||37||4||41
||女||9||17||20||46||4||50
||有効回答||23||29||31||83||8||91
||無効回答||2||3||2||7||0||7
||総数||25||32||33||90||8||98
−−−−
sa3() の戻り値は、tbl1に「総数」があるので、tbl2の方は NULL になります。この点、sa1, sa2 と少し違うので留意して下さい。
それから、sa3によって得られる表には「有効回答」という一種の小計欄があるので、カイ２乗検定等の統計検定の材料にはなりません。
!!(3) 単一回答集計のサンプルスクリプト
data01.csvを素材データとして、「性別」と「意見」のクロス集計を行うrubyスクリプトを掲げます。
Rプログラムの出力結果をそのまま標準出力に出力する他に、集計結果をrubyで受け取って、csvファイルとして書き出します。
−−−− tjs31.rb ここから
#! ruby -Ks
# 単一回答の集計(２項目のクロス集計) (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
source("tjsf.r")
dtf <- read.csv("data01.csv", header=T, na.strings="")
dtf$性別 <- factor(dtf$性別, levels=c("男", "女"), labels=c("男", "女"))
xx <- sa1(dtf, "性別", "意見")
for (i in 1:length(xx)) if (!is.null(xx[[i]])) xx[[i]] <- round(xx[[i]],1)
xx
xx$dnn <- names(dimnames(xx$tbl1))
robj(xx)
EOS
##
hs = Rrx.rexec(rpro)
puts "** Rプログラムの出力"
print hs[:sink]
printf("\n")
xx = hs["robj1"]
dnn = xx.delete("dnn") # 項目タイトル
dnnstr = (dnn.class == Array) ? dnn.join("/") : dnn
keys = %w(tbl1 tbl2 pct1 pct2)
titles = %w(基本表 「計」付きの表 横軸パーセント 縦軸パーセント)
str = ""
while key = keys.shift
title = titles.shift
val = xx[key]
next unless val
if dnn.class == Array # クロス集計表の場合
val[0][0] = dnnstr
else
val.delete_at(0)
val = Rrx.table_turn(val)
val.delete_at(0)
end
str = str + "** #{title}\n"
str = str + Rrx.ary2str(val, ",") + "\n"
end
File.open("result.csv", "w") {|ff| ff.write str}
puts "'result.csv' を書き出しました."
−−−− tjs31.rb ここまで
上は、２つの項目に着目したクロス集計(いわば２次元)の例ですが、１つの項目のみに着目した集計(１次元)にもそのまま使えます。
「xx$dnn <- names(dimnames(xx$tbl1))」というのは、縦・横の軸のタイトル名(次元名)を取得して、ruby側でそれを参照できるようにするためのものです。
次元名が２つあれば２次元(クロス集計)であり、１つだけなら１次元(１行のみの表)ということになります。
作られた表が２次元だとマトリクス(matrix)であり、１次元の場合は配列(array)です。それをruby側で受け取ると、どちらも配列(rubyのArray)の形になっています。しかし、こまったことに、単純に２次元配列｜１次元配列として識別できるわけではありません。
Rの配列(１次元)は、ruby側では２次元配列にになります。「配列の要素番号・項目名・度数」が複数入っている２次元配列です。Rのマトリクスもrubyでは２次元配列なので、その区別がつきにくい。
そこで、dnnに次元名をセットしてrubyに引きわたしています。ruby側において、dnnが配列なら２次元であり、単なる文字列であれば１次元と判断できます。
少々くどくなりましたが、ともあれ tjs31.rb は、２次元の表にも１次元の表にも対応します。ただし、３次元以上の表は作れません。
!!(4) 補完的関数
tjsf.rには補完的な役割の関数も入っているので、いくつか紹介します。
以下では、素材データとして data02.csv, data02.mem を用います。
data02.csvの中身は data01.csv と同じですが、列名が「性別→sex」、「意見→view」と英語名になっています。また、性別欄は「男→m」、「女→f」。意見欄は「賛成:1, 中立:2, 反対:3」の対応関係で数字が書き込まれています。
data02.mem は、その対応関係をメモ書きしたもので、中身は次の２行です。
sex | m,f | 男,女
view | 1,2,3 | 賛成,中立,反対
半角１文字の '|' は区切り文字です。タブコードも区切り文字になります。区切り文字の前後にある半角スペースは無視されます。
区切り文字によって３つの部分に分けられていますが、各々の部分が何を意味するかは見当がつくと思います。最初が列名、あとの２つは、カンマ区切り形式のlevelsとlabelsです。
このような素材データ(data02.csv, data02.mem)を基にして前掲のスクリプトと同じ結果を得るための tjs32.rb を下に掲げます。
−−−− tjs32.rb ここから
#! ruby -Ks
# 単一回答の集計(補完的関数の利用) (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
source("tjsf.r")
dtf <- read.csv("data02.csv", header=T, na.strings="")
dtf <- remake.factor(dtf, "data02.mem")
xx <- sa2(dtf, 2, 3)
xx <- rename.dnn(xx, "性別", "意見")
xx <- lround(xx, 1)
xx
xx$dnn <- get.dnn(xx)
robj(xx)
EOS
##
hs = Rrx.rexec(rpro)
puts "** Rプログラムの出力"
print hs[:sink]
printf("\n")
xx = hs["robj1"]
dnn = xx.delete("dnn") # 項目タイトル
dnnstr = (dnn.class == Array) ? dnn.join("/") : dnn
keys = %w(tbl1 tbl2 pct1 pct2)
titles = %w(基本表 「計」付きの表 横軸パーセント 縦軸パーセント)
str = ""
while key = keys.shift
title = titles.shift
val = xx[key]
next unless val
if dnn.class == Array # クロス集計表の場合
val[0][0] = dnnstr
else
val.delete_at(0)
val = Rrx.table_turn(val)
val.delete_at(0)
end
str = str + "** #{title}\n"
str = str + Rrx.ary2str(val, ",") + "\n"
end
File.open("result.csv", "w") {|ff| ff.write str}
puts "'result.csv' を書き出しました."
−−−− tjs32.rb ここまで
上のスクリプトに出てくる補完的な関数について、その仕様を記します。
◇ remake.factor
・利用例：dtf <- remake.factor(dtf, "data02.mem")
・機能：メモファイルに基づいてデータフレームの列をfactorとして再設定します。{{br}}
「remake.factor(dtf, "data02.mem")」は、data02.memの内容に即して、データフレームdtfの各列をfactorとして再設定するものです。{{br}}
この関数を用いずに、あえて別の形で記述するなら次の２行になります。{{br}}
dtf$sex <- factor(dtf$sex, levels=c("m","f"), labels=c("男","女"))
dtf$view <- factor(dtf$view, levels=1:3, labels=c("賛成","中立","反対"))
・戻り値：再設定された結果(データフレーム)を返します。指定されたメモファイルが存在しないなどの理由で再設定できない場合は、再設定せずに(つまり素のままで)データフレームを返します。
・引数：第１はデータフレーム、第２は再設定用メモ書きファイルの名前です。{{br}}
第３引数として、読み込むファイルの文字コードを指定することができます。data02.memがutf-8で書かれている場合は次のようにします。{{br}}
dtf <- remake.factor(dtf, "data02.mem", "utf-8"){{br}}
第３引数を省略すると、Rのデフォルトの文字コードになります。日本語版のWindowsだと "cp932" です。
◇ rename.dnn
・利用例：xx <- rename.dnn(xx, "性別", "意見")
・機能：テーブルの次元名を変構します。{{br}}
サンプルスクリプトにおいて、xx(リスト)には tbl1, tbl2, pct1, pct2 の４つのテーブルが入っているわけですが、各々のテーブルの次元名は、「sex, view」になっています。これを「性別, 意見」に変更するのが rename.dnn() です。４つのテーブルすべてについて変更を行います。
・戻り値：第１引数xxの次元名を変更した結果を返します。変更できない時は、xxをそのまま返します。
・引数：第１引数は、テーブルからなるリストを指定します。あるいは、テーブル１つを指定することもできます。{{br}}
第２引数と第３引数は次元名の指定です。第２が行名、第３が列名。２次元なら「"性別", "意見"」のように２つを指定します。１次元であれば "性別" のように１つだけを指定します(第３引数を省略します)。
◇ lround
・利用例：xx <- lround(xx, 1)
・機能：リストの各構成要素に対して、まるめ処理を行います。{{br}}
リストの構成要素として、テーブル・マトリクス・ベクトルを想定しますが、それ以外が含まれていても支障ありません。{{br}}
リストの構成要素が実数値(double)か複素数値(complex)でなければ、まるめ処理を施しません。
・戻り値：まるめ処理を施した結果を返します。
・引数：第１引数はリストを想定していますが、テーブル・マトリクス・ベクトルなどでもかまいません。{{br}}
第２引数は、まるめ処理に係る小数点の桁数を指定するものです。デフォルトは0です。つまり整数値になるようにまるめます。これを1にすると、小数点１桁までの数値になります。
◇ get.dnn
・利用例：dnn <- get.dnn(xx)
・機能：引数xxの次元名を取得します。{{br}}
xxがリストであれば、その中の最初にみつかったテーブルの次元名を返します。xxがテーブルなら、その次元名を返します。{{br}}
sa1()〜sa3()の戻り値(リスト)から次元名を得ることを目的として、この get.dnn() を設けました。
・戻り値：引数xxの次元名。例えば c("sex", "view") など。
・引数：リスト、またはテーブル
以上が補完的関数の説明です。
余談ですが、tjs31.rb, tjs32.rb のどちらも、小数点以下の桁数についてのまるめ処理をRプログラムの側で行っています。各々のテーブルに対して round() を適用しています。
この処理をR側ではなくruby側で行うサンプルスクリプトを tjs32b.rb として同梱しておきます。
参考まで、rubyスクリプトの該当部分を下に掲げます。
−−−− tjs32b.rb 抜粋ここから
str = ""
while key = keys.shift
title = titles.shift
val = xx[key]
next unless val
if dnn.class == Array # クロス集計表の場合
val[0][0] = dnnstr
else
val.delete_at(0)
val = Rrx.table_turn(val)
val.delete_at(0)
end
if key =~ /^pct[12]$/
val.map! {|row|
row.map! {|cell|
[Fixnum,Float].include?(cell.class) ? sprintf("%.1f",cell) : cell
}
}
end
str = str + "** #{title}\n"
str = str + Rrx.ary2str(val, ",") + "\n"
end
−−−− tjs32b.rb 抜粋ここまで
------------------------------------------------------------------------
!2. 複数回答の集計
ここでは、「主な情報源」として、新聞・雑誌・テレビ・ラジオの４つの選択肢から、当てはまるものを複数個選んでもらうような複数回答の集計を取り上げます。
素材データとして data03.csv を用います。これは、ID(整数値), 年齢層(20代・30代・40代のどれか１つ), 主な情報源(新聞・雑誌・テレビ・ラジオの複数回答)から構成されています。98人分のデータです。具体的には次のようになっています。
ID,年齢層,新聞,雑誌,テレビ,ラジオ
1,30代,1,,,
2,40代,1,1,,
3,40代,,1,,
年齢層の列は単一回答です。
主な情報源の列は、４項目のうち選択された項目のところに数字の1を入力し、選択されなかったところは空欄にしてあります。
data03.csv は、乱数発生によって作成したデータです。実態とは関係ありません。
複数回答にまつわる主な関数は、ma1(), ma2() の２つです。他に ma.refactor(), ma.count() があります。
!!(1) 複数回答と単一回答の組合せ集計 ma1
まず、年齢層(単一回答)と主な情報源(複数回答)の組合せ集計を行います。
ma1()という関数をtjsf.rに盛り込んだので、それを用います。そのサンプル(Rプログラム部分のみ)を掲げます。
−−−− Rプログラムここから
source("tjsf.r")
dtf <- read.csv("data03.csv", header=T, na.strings="")
xx <- ma1(dtf, 2, 3:6)
xx
−−−− Rプログラムここまで
ma1() の戻り値は、tbl, pct の２つの表からなるリストです。前者は度数(選択者数)、後者は比率(パーセント)の表です。
上のプログラムの実行結果は次のとおり。
−−−−
$tbl
新聞 雑誌 テレビ ラジオ 該当数
20代 13 17 14 14 31
30代 24 17 20 16 34
40代 15 22 17 16 33
合計 52 56 51 46 98
$pct
新聞 雑誌 テレビ ラジオ 該当数
20代 41.93548 54.83871 45.16129 45.16129 100
30代 70.58824 50.00000 58.82353 47.05882 100
40代 45.45455 66.66667 51.51515 48.48485 100
合計 53.06122 57.14286 52.04082 46.93878 100
−−−−
ma1() を呼び出す場合、データフレームの列を番号でなく名前で指定することができます。次のように書くことができます。
ma1(dtf, "年齢層", c("新聞","雑誌","テレビ","ラジオ"))
第１引数がデータフレーム、第２引数と第３引数が列の指定です。第２・第３は、順序を入れ換えても差し支えありません。
この ma1() は、次のサイトにある ma() を基にして作りました(というより、変更を加えました)。作者に感謝します。
[[R -- マルチアンサーのクロス集計|http://aoki2.si.gunma-u.ac.jp/R/multianswer.html]]
!!(2) 複数回答の列をfactorとして再設定 ma.refactor
複数回答の列の値について少し触れます。
data03.csvの場合、選択されたら数値の1、そうでなければ空欄になっています。空欄は、read.csvで読み込んだ後に欠損値のNAとなります。
もし空欄でなく数字の0が入力されていたとすると、前述のRプログラムは、意図したような結果を出してくれません。
数字の1と0を機械的に列べると、当然ながら0の方が先にきます。そうすると、ma1() は、0の方を「選択あり」であるかのように処理してしまいます。
前もって、0よりも1の方が先にくるものであることをRに知らせてやる必要があります。そのためには、複数回答の列をfactorとして再設定します。
例えば、0と1からなる３列目を再設定するには次のようにします。
dtf[,3] <- factor(dtf[,3], levels=1:0, labels=1:0)
上の記述を３列目だけでなく３〜６列すべてに対して実行すればいいわけです。
これを少し簡単に記述できるようにするため、ma.refactor() という関数を設けました。次のように呼び出します。
dtf <- ma.refactor(dtf, 3:6, 1:0, 1:0)
こうすると、３〜６列のいずれも factor として再設定されます。
第１引数はデータフレーム、第２が列の指定(番号でなく名前でも可)、第３がlevels、第４はlabelsです。
第４引数を省略すると、「c("新聞○", "新聞×")」が指定されたものとみなされます。「新聞」のところは各々の列名に置き換えられます。
data03.csvの場合でいうと、数字の0ではなく空欄なので、再設定は次のようにします。
dtf <- ma.refactor(dtf, 3:6, c(1,NA), 1:0)
ここで、ma1(), ma.refactor() の利用例を掲げておきます。
Rプログラムの出力を標準出力に出力するほか、ma1() の戻り値をrubyの側で受け取って、それをcsvファイルとして書き出します。
−−−− tjs33.rb ここから
#! ruby -Ks
# 複数回答と単一回答の組合せ集計 (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
source("tjsf.r")
dtf <- read.csv("data03.csv", header=T, na.strings="")
dtf <- ma.refactor(dtf, 3:6, c(1,NA))
xx <- ma1(dtf, 2, 3:6)
xx <- lround(xx, 1)
xx
robj(xx)
EOS
##
hs = Rrx.rexec(rpro)
puts "** Rプログラムの出力"
print hs[:sink]
printf("\n")
xx = hs["robj1"]
keys = %w(tbl pct)
titles = %w(回答の度数分布 回答比率)
str = ""
while key = keys.shift
title = titles.shift
val = xx[key]
next unless val
str = str + "** #{title}\n"
str = str + Rrx.ary2str(val, ",") + "\n"
end
File.open("result.csv", "w") {|ff| ff.write str}
puts "'result.csv' を書き出しました."
−−−− tjs33.rb ここまで
!!(3) 複数回答と単一回答の分解的集計
ma1() は、新聞・雑誌・テレビ・ラジオの４項目からなる複数回答の状況を一括して集計するので、全体像をみるのに向いています。ただ、その集計結果をそのまま統計検定にかけることはできません。
そこで、「新聞」にまるを付けたか否か、「雑誌」にまるを付けたか否かというように、各々の選択肢について二者択一でとらえ、集計表を４つ作ります。そうすると、各々の表をカイ２乗検定等にかけることができます。
前述したように、ma.refactor() の第４引数(labelsの指定)を省略すると、「新聞○」と「新聞×」のように設定されるので、あとは単一回答の集計を行う sa1() を適用してやれば、意図した表を作れます。
下に、それら４つの表を作成するスクリプトを掲げます。
単に表を作るだけでは何なので、カイ２乗検定にかけて、その結果を出力してみます。有意性が認められた時は、調整済み残差の表も出力します。
Rプログラムの実行結果は標準出力に出力し、ruby側では同じ内容をcsvファイルに書き出しています。
まだ紹介していないRの関数が出てきますが、後で説明します。
−−−− tjs34.rb ここから
#! ruby -Ks
# 複数回答と単一回答の分解的集計 (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
source("tjsf.r")
dtf <- read.csv("data03.csv", header=T, na.strings="")
dtf <- ma.refactor(dtf, 3:6, c(1,NA))
items <- c("χ2", "df", "p-value", "cramer")
keys <- c("statistic", "parameter", "p.value", "cramer")
for (i in 3:6) {
xx <- sa1(dtf, 2, i)
zz <- chisq.test(xx$tbl1)
zz <- append.cramer(zz)
zz$join <- join.tbl(xx$tbl2, xx$pct1, "%d(%.1f)")
ww <- sapply(keys, function(k) zz[[k]])
names(ww) <- items
print(zz$join); cat("\n")
print(ww); cat("\n")
if (zz$p.value < 0.05) {
cat("☆ 調整済み残差\n")
print(zz$stdres); cat("\n")
}
cat("\n\n")
robj(zz)
}
EOS
##
hs = Rrx.rexec(rpro)
puts "** Rプログラムの出力"
print hs[:sink]
printf("\n")
items = %w(χ2 df p-value cramer)
keys = %w(statistic parameter p.value cramer)
str = ""
Rrx.robj_each() do |kk, zz|
str = str + Rrx.ary2str(zz["join"], ",") + "\n"
ary = []
keys.each {|key| ary << sprintf("%g", zz[key])}
str = str + Rrx.ary2str([items, ary], ",") + "\n"
if zz["p.value"] < 0.05
str = str + "☆ 調整済み残差\n"
str = str + Rrx.ary2str(zz["stdres"], ",") + "\n"
end
str = str + "\n\n"
end
File.open("result.csv", "w") {|ff| ff.write str}
puts "'result.csv' を書き出しました."
−−−− tjs34.rb ここまで
参考まで、上のスクリプトを実行した結果の一部を掲げておきます。
「年齢層」と「新聞」のクロス集計表で、１つのセルが「13(41.9)」のように、度数と括弧書きのパーセントという形になっています。
クロス集計表の次は、カイ２乗検定の結果(カイ２乗値・自由度・有意確率・クラメール係数)を列記したものです。
−−−−
新聞
年齢層 新聞○ 新聞× 総数
20代 13(41.9) 18(58.1) 31(100.0)
30代 24(70.6) 10(29.4) 34(100.0)
40代 15(45.5) 18(54.5) 33(100.0)
総数 52(53.1) 46(46.9) 98(100.0)
χ2 df p-value cramer
6.50090607 2.00000000 0.03875665 0.25755733
−−−−
余談ですが、data03.csvは乱数発生によって作ったものなので、カイ２乗検定で有意性が検出されることはないだろうと予想しましたが、上の「新聞」については有意性が認められました。有意性が検出されたのはこれだけです。
調整済み残差をみるまでもなく、「30代」の「新聞○」が多く、「新聞×」が少ない様子が分かります。
tjs34.rbで新たに出てきた関数について説明します。２つあります。
append.cramer(zz) は、引数zz(リスト)にクラメール係数を追加して返すものです。zzには、予めカイ２乗検定の結果を代入しておく必要があります。つまり chisq.test() の戻り値を代入しておきます。
「zz <- append.cramer(zz)」としたとき、zz$cramer にクラメール係数の値がセットされています。
もう１つの join.tbl() は、２つの表を結合するものです。典型的な例として、度数分布の表と構成比の表の結合があります。上のサンプルスクリプトもその一例です。
tt <- join.tbl(tbl, pct, "%d(%.1f)")
とすると、tbl, pct の２つの表が結合されて、その結果が tt に代入されます。ここで結合というのは、対応するセル１つづつを組み合わせるという意味です。
第３引数は、結合セルの書式を指定するものです。これが内部処理においてsprintfのフォーマットとして使われます。"%d(%.1f)" の場合、tblのセルが整数値として書き出され、pctのセルが括弧書きの小数点数(有効桁数１桁)で書き出されます。
もし pct が表として tbl よりも小さくて、行数または列数が少なければ、その少ない部分については tbl のセルだけが書き出されます。その場合、第３引数の書式は用いません。
rubyスクリプトの部分については説明を省略します。
拙作rrxライブラリのメソッドをいくつか使っているので分かりにくいと思いますが、必要に応じ関連のマニュアルを参照して下さい。
◇ append.cramer
・利用例：zz <- chisq.test(xx); zz <- append.cramer(zz)
・機能：引数zzにクラメール係数を追加して返します。zzには予めカイ２乗検定の結果を代入しておきます。
・戻り値：クラメール係数を含むカイ２乗検定結果を返します。戻り値zzのzz$cramerがクラメール係数の値です。
・引数：カイ２乗検定結果を引数として与えます。
◇ join.tbl
・利用例：tt <- join.tbl(tbl, pct, "%d(%.1f)")
・機能：２つの表を結合します。対応するセルを組み合わせます。
・戻り値：結合した表を返します。
・引数：第１引数と第２引数は表データ、第３引数は結合セルの書式で、内部処理においてsprintfのフォーマットとして使われます。省略時は "%s(%s)" とみなされます。
!!(4) 複数回答の選択肢相互の関係をみる ma2
「主な情報源」の４つの選択肢(新聞・雑誌・テレビ・ラジオ)の相互の関係をみるには、「新聞」と「雑誌」の両方にまるを付けた人がどれくらいいるか、「新聞」と「テレビ」の両方にまるを付けた人はどうか、といった確認をすることになると思います。そうすると、例えば、「新聞」を選ぶ人は「ラジオ」も選ぶ傾向が高いといったことが分かるかもしれません。
この「両方にまるを付けた人がどれくらいか」を一覧表示するのが ma2() です。具体的には次のような表を生成します。
−−−−
新聞 雑誌 テレビ ラジオ 該当数
新聞 52 29 23 28 52
雑誌 29 56 30 25 56
テレビ 23 30 51 22 51
ラジオ 28 25 22 46 46
該当数 52 56 51 46 98
−−−−
上の表で、「該当数」というのは、いわば「総数」に当たるものです。「新聞」にまるを付けた人の総数が何人か、「雑誌」の総数は何人かなどが「該当数」をみれば分かります。「該当数」の最も右下の98は、回答者総数を意味します。
「新聞」と「雑誌」が交差するセルの値は29ですが、これは、この両方にまるを付けた人が29人いることを意味します。
この表は、縦軸と横軸を入れ換えても同じものになります。
ma2() は、このような表を返しますが、４種類の表を生成して返します。
ma2() の仕様について説明する前に、その利用例を下に掲げます。Rプログラム部分のみです。
−−−− Rプログラムここから
source("tjsf.r")
dtf <- read.csv("data03.csv", header=T, na.strings="")
dtf <- ma.refactor(dtf, 3:6, c(1,NA))
xx <- ma2(dtf, 3:6)
xx <- lround(xx, 1)
xx
−−−− Rプログラムここまで
まず ma2 を呼び出す時の引数ですが、第１引数がデータフレーム、第２引数は複数回答の列の指定です。列の指定は、列番号でも列名でも、どちらでもかまいません。
ma2() の戻り値は、次の４つの表を含むリストです。
* tbl1: 「該当数」を含まない表
* tbl2: 「該当数」を含む表
* pct1: 右端の列の「該当数」を100%とする時のパーセント
* pct2: 右下端の「回答者総数」を100%とする時のパーセント
先に掲げた表は、「該当数」付きの tbl2 です。tbl1 は、それから縦と横の「該当数」を取り除いたものです。
pct1, pct2 を下に掲げておきます。
−−−−
$pct1
新聞 雑誌 テレビ ラジオ 該当数
新聞 100.0 55.8 44.2 53.8 100
雑誌 51.8 100.0 53.6 44.6 100
テレビ 45.1 58.8 100.0 43.1 100
ラジオ 60.9 54.3 47.8 100.0 100
該当数 53.1 57.1 52.0 46.9 100
$pct2
新聞 雑誌 テレビ ラジオ 該当数
新聞 53.1 29.6 23.5 28.6 53.1
雑誌 29.6 57.1 30.6 25.5 57.1
テレビ 23.5 30.6 52.0 22.4 52.0
ラジオ 28.6 25.5 22.4 46.9 46.9
該当数 53.1 57.1 52.0 46.9 100.0
−−−−
pct1 は、「新聞」にまるを付けた人の総数を100%とした時に、「雑誌」にもまるを付けた人が何％か、「テレビ」にもまるを付けた人が何％かなどをみるのに使えます。
pct2 の方は、回答者総数98人を100%とした時に、各セルの値が何％かを確認するのに使います。全体の中で「両方にまるを付けた人」が何％いるかが分かります。
４つの表のうち pct1 以外は、縦軸と横軸を入れ換えても同じものになります。pct1 は同じになりません。
こうした表は、きちんと分析するのにはあまり適していませんが、ざっと全体像をつかむ時に、それなりに使えるのではないかと思います。
下に ma2() の利用サンプルを掲げておきます。
−−−− tjs35.rb ここから
#! ruby -Ks
# 複数回答の相互関係をみる (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
source("tjsf.r")
dtf <- read.csv("data03.csv", header=T, na.strings="")
dtf <- ma.refactor(dtf, 3:6, c(1,NA))
xx <- ma2(dtf, 3:6)
xx <- lround(xx, 1)
xx
robj(xx)
EOS
##
hs = Rrx.rexec(rpro)
puts "** Rプログラムの出力"
print hs[:sink]
printf("\n")
xx = hs["robj1"]
keys = %w(tbl1 tbl2 pct1 pct2)
titles = %w(基本表 「該当数」付きの表 横軸パーセント 対総数のパーセント)
str = ""
while key = keys.shift
title = titles.shift
val = xx[key]
next unless val
str = str + "** #{title}\n"
str = str + Rrx.ary2str(val, ",") + "\n"
end
File.open("result.csv", "w") {|ff| ff.write str}
puts "'result.csv' を書き出しました."
−−−− tjs35.rb ここまで
!!(5) 複数回答の選択肢相互の関係を分解的にみる
前述の ma2() が生成する表は、統計検定にかけることができません。そこで、各選択肢をばらばらにして集計してみます。
「新聞」と「雑誌」の２つに注目した時の表は次のとおり。
−−−−
|| ||雑誌○||雑誌×
||新聞○||29||23
||新聞×||27||19
−−−−
選択肢が４つなので、そこから２つを抜き出して表を作るやり方は６通りあります。
上のような２行・２列の表を６つ生成して、各々を統計検定にかけるためのサンプルを下に掲げます。カイ２乗検定を用います。
−−−− tjs36.rb ここから
#! ruby -Ks
# 複数回答の相互関係を分解的にみる (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
source("tjsf.r")
dtf <- read.csv("data03.csv", header=T, na.strings="")
dtf <- ma.refactor(dtf, 3:6, c(1,NA))
items <- c("χ2", "df", "p-value", "cramer")
keys <- c("statistic", "parameter", "p.value", "cramer")
for (i in 3:5) {
for (j in (i+1):6) {
xx <- sa1(dtf, i, j)
zz <- chisq.test(xx$tbl1)
zz <- append.cramer(zz)
zz$join <- join.tbl(xx$tbl2, xx$pct1, "%d(%.1f)")
ww <- sapply(keys, function(k) zz[[k]])
names(ww) <- items
print(zz$join); cat("\n")
print(ww); cat("\n")
if (zz$p.value < 0.05) {
cat("☆ 調整済み残差\n")
print(zz$stdres); cat("\n")
}
cat("\n\n")
robj(zz)
}
}
EOS
##
hs = Rrx.rexec(rpro)
puts "** Rプログラムの出力"
print hs[:sink]
printf("\n")
items = %w(χ2 df p-value cramer)
keys = %w(statistic parameter p.value cramer)
str = ""
Rrx.robj_each() do |kk, zz|
str = str + Rrx.ary2str(zz["join"], ",") + "\n"
ary = []
keys.each {|key| ary << sprintf("%g", zz[key])}
str = str + Rrx.ary2str([items, ary], ",") + "\n"
if zz["p.value"] < 0.05
str = str + "☆ 調整済み残差\n"
str = str + Rrx.ary2str(zz["stdres"], ",") + "\n"
end
str = str + "\n\n"
end
File.open("result.csv", "w") {|ff| ff.write str}
puts "'result.csv' を書き出しました."
−−−− tjs36.rb ここまで
複数回答の答えは「○」か「×」かの２通りなので、各々の表は２行・２列になります。なので、カイ２乗検定の他に、フィッシャーの直接確率検定を適用することができます。
４つのセルの中に５未満のものがあるなど度数が少数のケースでは、カイ２乗検定を行った時に、警告メッセージが出ることがあります。そうした時はフィッシャーの検定を使うのがいいかもしれません。
zz <- fisher.test(xx$tbl1)
print(zz$p.value)
このようにすると、p値が表示されます。これが 0.05 未満であれば、有意性が認められることになります。
有意性が認められれば、例えば、「新聞を主な情報源にする人は、ラジオも情報源にする傾向がある。」といったことが言えるかもしれません。しかし、有意性が認められないのであれば、少なくとも統計検定の結果からは、そうした傾向を指摘する根拠が得られない、ということになります。ちなみに、６つの表の中で有意性が認められるものはありませんでした。
このような分解的集計・検定の他に、数量化IIIや対応分析を用いることもできます。そのやり方について「単純集計に関する覚え書き／複数回答の集計」でごく簡単に触れましたが、関連の書籍やWebがいろいろあるので参考にして下さい。
!!(6) 複数回答の選択個数をみる ma.count
複数回答の選択肢について、回答者各人がまるを何個つけたかを検出しておくと便利です。
IDが１番の人は２個、２番の人は０個、３番の人は３個などのように、各人の選択個数を把握します。そうすれば、選択個数０の人(いわば無回答の人)を除いて集計したりすることが容易になります。
この選択個数の検出を行う関数 ma.count() を設けました。
これを呼び出す時の引数は４つありますが、第１引数がデータフレーム、第２が複数回答の列番号(または列名)です。多くの場合、第３と第４は省略可能です。
とりあえず ma.count の利用例を掲げます。Rプログラム部分のみです。
−−−− Rプログラムここから
source("tjsf.r")
dtf <- read.csv("data03.csv", header=T, na.strings="")
count <- ma.count(dtf, 3:6)
table(count)
−−−− Rプログラムここまで
上のプログラムを実行すると、下のような表示が得られます。
−−−−
count
0 1 2 3 4
4 24 34 31 5
−−−−
選択個数０の人が４人、選択個数１個の人が24人いることが分かります。
ma.count() の戻り値はベクトルです。dtfには98人分のデータが入っているので、戻り値である count にも98人分の選択個数が入ります。
前述の例では「ma.count(dtf, 3:6)」のように、引数を２つしか指定しませんでしたが、第３引数と第４引数を省略せずに指定すると、次のようになります。
count <- ma.count(dtf, 3:6, 1, FALSE)
第３引数の１は、素材データであるdtfにおいて、複数回答の「選択あり」を示す入力値です。
dtfでは、選択された時に数字の１、選択されなかった時は空欄になっています。そこで、選択個数をカウントする際に数字の１日着目するというのが第３引数の意味です。
dtf <- ma.refactor(dtf, 3:6, c(1,NA), c("yes","no"))
上のように複数回答の列をfactorとして再設定していると、「選択あり」を示すのは "yes" になるので、ma.count() の呼び出し方は次のようになります。
count <- ma.count(dtf, 3:6, "yes", FALSE)
第３引数が "yes" のように文字である場合、実は、それを正規表現とみなします。内部処理では、Rのgrep()によってマッチングを行います。なので、厳密に「yes」の３文字を示すなら "^yes$" とすべきかもしれません。しかし、もう１つの値が "no" なので、そこまでしなくても意図どおりの結果が得られます。
第４引数は、選択個数を素材データのdtfに挿入・追加するか否かの指定です。デフォルトは FALSE です。つまり、挿入・追加をせず、選択個数を記録したベクトルを戻り値として返します。
TRUEにすると、dtfに nnn という名前の列を設けて、それに選択個数を記録します。戻り値として、nnnを追加・挿入したデータフレームを返します。
戻り値がdtf2に代入されている場合、dtf2$nnn または dtf2[,"nnn"] によって選択個数を参照できます。既に素材データのdtfに nnn という名前の列がある時は、そのnnnを上書きします。
選択個数０の無回答の人を除いたデータフレームを得たい時は、次のようにします。
dtf2 <- ma.count(dtf, 3:6, 1, TRUE)
dtf2 <- subset(dtf2, nnn > 0)
上のようにすると、無回答の人を除く94人分のデータが dtf2 にセットされます。
以上が ma.count() の仕様です。
nnn という列名には違和感があるかもしれませんが、count だと、よく使われるような感じがして、競合を避けるため nnn にしました。
!!(7) 複数回答にまつわる関数の一覧
複数回答にまつわる関数の一覧を掲げておきます。
◇ ma1
・利用例：xx <- ma1(dtf, 2, 3:6)
・機能：複数回答と単一回答の組合せ集計を行います。
・戻り値：tbl, pct の２つの表からなるリストを返します。前者は度数、後者はパーセント。
・引数：第１引数はデータフレーム、第２と第３は列の指定です。どちらかが複数回答の列、もう一方は単一回答の列です。その順番は問いません。
◇ ma2
・利用例：xx <- ma2(dtf, 3:6, 1, FALSE)
・機能：複数回答の選択肢相互の関係をみるための表を生成します。
・戻り値：tbl1, tbl2, pct1, pct2 の４つの表からなるリストを返します。４つの表が何かは、本文を参照して下さい。
・引数：第１引数はデータフレーム、第２引数は複数回答の列の指定です。
◇ ma.refactor
・利用例：dtf <- ma.refactor(dtf, 3:6, c(1,NA), 1:0)
・機能：複数回答の列をfactorとして再設定します。
・戻り値：再設定された結果(データフレーム)を返します。
・引数：第１引数はデータフレーム、第２引数は列の指定、第３はlevels、第４はlabelsです。{{br}}
第４引数を省略すると、「c("新聞○", "新聞×")」が指定されたものとみなされます。「新聞」のところは各々の列名に置き換えられます。
◇ ma.count
・利用例：count <- ma.count(dtf, 3:6, 1, FALSE){{br}}
dtf2 <- ma.count(dtf, 3:6, "○", TRUE)
・機能：複数回答の選択個数を検出します。
・戻り値：第４引数がFALSEの時は、選択個数を記録したベクトルを返し、TRUEの時は、選択個数を記録した列(nnn)が追加・挿入されたデータフレームを返します。
・引数：第１はデータフレーム、第２は複数回答の列の指定、第３が「選択あり」を示す値(文字だと正規表現パターン)、第４は列の追加・挿入の有無です。
------------------------------------------------------------------------
!3. 数値データの集計
ここでは、数値データとして「年齢」と「支出」(月当たりの支出額)を含むcsvデータを材料にします。
data04.csvは、ID(整数値)、性別(m,f)、年齢(数値)、支出(数値)から構成されています。98人分のデータです。具体的には次のようになっています。
ID,性別,年齢,支出
1,m,29,82249
2,m,49,93181
3,f,42,93647
数値データの集計にまつわる主な関数は、va1(), va2() の２つです。どちらも、単一回答と数値データの組合せの集計に用います。単一回答の「性別」(男女別)で分類しながら、数値データの平均値等を算出します。
他に date.split(), var2(), sd2() をtjsf.rに盛り込んであります。
date.split() は、「単純集計に関する覚え書き／数値データの集計」で触れたので今回は言及しません。仕様は同じです。
var2(), sd2() は、それぞれ標本分散と標本標準偏差を返す関数です。欠損値NAを取り除いた上で値を算出します。わたす引数は１つのみです。
ちなみに、Rに最初から組み込まれている var(), sd() は、それぞれ不偏分散と不偏標準偏差を返す関数です。こちらは引数が１つとはかぎりません。
!!(1) １つの値のみを返す関数を用いた集計 va1
mean() は平均値を返す関数ですが、このように１つの値のみを返す関数によって集計を行う際、va1() を使えます。
data04.csvを素材にして、男女別に平均年齢を算出し、それを表として得たい場合、次のようにします。Rプログラム部分のみ示します。
−−−− Rプログラムここから
source("tjsf.r")
dtf <- read.csv("data04.csv", header=T, na.strings="")
dtf$性別 <- factor(dtf$性別, levels=c("m","f"), labels=c("男","女"))
xx <- va1(dtf, "性別", "年齢", mean)
print(xx)
−−−− Rプログラムここまで
上のプログラムを実行した結果は下のとおりです。va1() は matrix を返します。
−−−−
年齢
男 34.37500
女 34.04762
全体 34.23469
−−−−
va1() の第１引数は、素材となるデータフレームです。
第２引数は、分類の手がかりとなる単一回答の列を指定します。名前でも番号でもどちらでもかまいません。
第３引数では数値データの列を指定します。「c("年齢","支出")」とか「3:4」などのように複数指定してもかまいません。ちなみに、第２引数の方は１つだけです。
最後の第４引数は、適用したい関数の名前を指定します。省略すると mean が指定されたものとみなされます。
関数は、既存のものでなく自作のものでもかまいません。例えば、１つのセルを「平均値(sd:標準偏差)」の書式で出力したい時は次のようにします。
xx <- va1(dtf, 2, 3:4, function(x)
sprintf("%.1f(sd:%.2f)", mean(x), sd(x)))
上のようにすると、下の結果が得られます。
−−−−
年齢 支出
男 "34.4(sd:8.23)" "105605.7(sd:29094.26)"
女 "34.0(sd:7.41)" "92032.0(sd:27015.33)"
全体 "34.2(sd:7.85)" "99788.4(sd:28879.60)"
−−−−
下にサンプルスクリプトを掲げておきます。平均値と中央値を別々に表示するものです。
Rプログラムの出力結果を標準出力に出力し、ruby側では同じ内容をcsvとして書き出します。
−−−− tjs37.rb ここから
#! ruby -Ks
# １つの値を返す関数による集計 (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
source("tjsf.r")
dtf <- read.csv("data04.csv", header=T, na.strings="")
dtf$性別 <- factor(dtf$性別, levels=c("m","f"), labels=c("男","女"))
xx1 <- va1(dtf, 2, 3:4, mean)
cat("平均値\n")
xx1
cat("\n")
xx2 <- va1(dtf, c("年齢","支出"), "性別", median)
cat("中央値\n")
xx2
robj(list(平均値=xx1, 中央値=xx2))
EOS
##
hs = Rrx.rexec(rpro)
puts "** Rプログラムの出力"
print hs[:sink]
printf("\n")
str = ""
Rrx.robj_each(3) do |name, key, val|
str = str + "** #{key}\n"
str = str + Rrx.ary2str(val, ",") + "\n\n"
end
File.open("result.csv", "w") {|ff| ff.write str}
puts "'result.csv' を書き出しました."
−−−− tjs37.rb ここまで
!!(2) 複数の値を返す関数を用いた集計 va2
va2() は、使い方は va1() と同じです。
第４引数で指定する関数が、１つの値だけでなく複数の値を返す場合に va2() を用います。
例えば、数値データに関して summary() は６つの値を返しますが、これを関数として指定したい時は次のようにします。
xx <- va2(dtf, 2, 3:4, summary)
実は、第４引数を省略すると、summary が指定されたものとみなされます。
va2() は、戻り値として list を返します。list の各要素は matrix です。
以下にサンプルスクリプトとその出力結果を示します。
−−−− tjs38.rb ここから
#! ruby -Ks
# 複数の値を返す関数を用いた集計 (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
source("tjsf.r")
dtf <- read.csv("data04.csv", header=T, na.strings="")
dtf$性別 <- factor(dtf$性別, levels=c("m","f"), labels=c("男","女"))
xx <- va2(dtf, 2, 3:4, summary)
xx
robj(xx)
EOS
##
hs = Rrx.rexec(rpro)
puts "** Rプログラムの出力"
print hs[:sink]
printf("\n")
str = ""
Rrx.robj_each(3) do |name, key, val|
str = str + "** #{key}\n"
str = str + Rrx.ary2str(val, ",") + "\n\n"
end
File.open("result.csv", "w") {|ff| ff.write str}
puts "'result.csv' を書き出しました."
−−−− tjs38.rb ここまで
−−−− 出力結果ここから
$年齢
Min. 1st Qu. Median Mean 3rd Qu. Max.
男 20 28.00 34.5 34.38 41.25 49
女 21 27.25 33.0 34.05 41.00 49
全体 20 28.00 34.0 34.23 41.00 49
$支出
Min. 1st Qu. Median Mean 3rd Qu. Max.
男 50430 82960 115200 105600 129600 149400
女 52490 69530 92120 92030 114900 148400
全体 50430 73590 98100 99790 124500 149400
−−−− 出力結果ここまで
数値データについては、例えば「年齢」を「20代・30代・40代」のようにカテゴリー化するのもよくあることだと思いますが、これを関数化するのは難しいのでやってません。
カテゴリー化の方法については「単純集計に関する覚え書き／数値データの集計」のところで触れているので、そちらを参考にして下さい。
------------------------------------------------------------------------
以上で「単純集計に関する覚え書き／まとめ(関数化)」をおわりにします。
tjsf.r の中には、上で紹介していない関数も含まれていますが、それらは内部処理で用いるために設けたサブルーチン的な関数です。興味があったら覗いてみて下さい。
Copyright (C) T. Yoshiizumi, 2013 All rights reserved.
〜 rubyとRを用いた処理 〜
最終更新日: 2013/08/11
以下に掲げるドキュメントは、
[[tjs_d.zip|http://cup.sakura.ne.jp/tjs_d.zip]]
に同梱されている tjs_d.txt と同じ内容です。{{br}}
tjs_d.zipには、euc-jp | utf-8 で書かれたサンプルスクリプトも同梱されています。{{br}}
サンプルスクリプトは、MS-Windowsだけでなく、linuxなどでも動作します
＜はじめに＞
「単純集計に関する覚え書き」を３回にわたって記しました。予想していたよりも長い解説になってしまい、読み返すのが面倒です。
具体的な集計の必要に迫られたとき、どのサンプルスクリプトを利用したらいいのか、どんなふうに書き換えたらいいのか、私自身が迷うくらいにだらだらとした解説になってしまいました。
そこで、「最低限これくらいは」というものをRの関数として定義しました。例えば、単一回答の集計の一つ(無回答を除く集計)の sa1() は次のように呼び出します(saは single answer のつもり)。
xx <- sa1(dtf, 2, 3)
print(xx)
こうすると、データフレームdtfの２列目と３列目に着目して、クロス集計表を作成します。
戻り値のxxには４つの表が入っています。２つの度数分布表と２つのパーセントの表です。
こうした sa1() のような関数をいくつか作って tjsf.r に入れました。以下で、それらの解説を記します。
なお、この tjsf.r は、rubyに依存する部分はありません。rubyを使わない場合でも利用可能です。
--------
{{toc_here}}
--------
!1. 単一回答の集計
ここでは、「性別」(男・女)、「意見」(賛成・中立・反対)のように、選択肢の中から１つだけ選んでもらう単一回答の集計を取り上げます。
適切な選択肢がないため回答できない｜回答したくないというケースもありますが、それは無回答という形になって表れると思います。そうした無回答を含めて集計する場合についても取り上げます。
単一回答にまつわる主な関数は sa1(), sa2(), sa3() の３つです。
!!(1) 単一回答の集計に関する関数 sa1
ここで扱い素材データは、ID(整理番号)、性別(男・女)、意見(賛成・中立・反対)の３列からなるcsvデータ(data01.csv)です。98人分が記録されています。具体的には次のような内容です。
ID,性別,意見
1,女,賛成
2,男,賛成
3,男,賛成
IDには空欄がありませんが、性別と意見の列には少数ながら空欄があります。つまり、欠損値(無回答)があります。
この素材データを集計します。
まず、欠損値を数に入れずに集計する sa1() のサンプルを示します。Rプログラム部分のみです。
−−−− sa1()の利用例ここから
source("tjsf.r")
dtf <- read.csv("data01.csv", header=T, na.strings="")
dtf$性別 <- factor(dtf$性別, levels=c("男", "女"), labels=c("男", "女"))
xx <- sa1(dtf, "性別", "意見")
xx
−−−− sa1()の利用例ここまで
「dtf$性別 <- factor(……)」は、なくても支障はありませんが、これがないと、作成される表の中の「男・女」の順番が「女・男」になります。
sa1() が返す値はリストで、次の４つのテーブルを含みます。
・tbl1: 総数(計)なしの表
・tbl2: 総数(計)を付加した表
・pct1: 横軸に足し算を行った「計」に対するパーセントの表
・pct2: 縦軸に足し算を行った「計」に対するパーセントの表
Rプログラムの出力結果は次のとおり。
−−−−
$tbl1
意見
性別 賛成 中立 反対
男 14 12 11
女 9 17 20
$tbl2
意見
性別 賛成 中立 反対 総数
男 14 12 11 37
女 9 17 20 46
総数 23 29 31 83
$pct1
意見
性別 賛成 中立 反対 総数
男 37.83784 32.43243 29.72973 100.00000
女 19.56522 36.95652 43.47826 100.00000
総数 27.71084 34.93976 37.34940 100.00000
$pct2
意見
性別 賛成 中立 反対 総数
男 60.86957 41.37931 35.48387 44.57831
女 39.13043 58.62069 64.51613 55.42169
総数 100.00000 100.00000 100.00000 100.00000
−−−−
sa1() の第２引数("性別")と第３引数("意見")は、列名でなく列の番号で指定してもかまいません。次のように書いても同じ結果になります。
xx <- sa1(dtf, 2, 3)
tbl1だけを取り出したい時は、xx$tbl1 などのようにします。これには「計」がないので、カイ２乗検定等にかけることができます。例えば次のように書きます。
xx <- sa1(dtf, 2, 3)
zz <- chisq.test(xx$tbl1)
こうすると、zzにカイ２乗検定の結果が代入されます。
tbl1, tbl2 などの各々のテーブルには、縦軸のタイトル「性別」と横軸のタイトル「意見」が付いています。これは、data01.csvの列名をそのまま割り当てたものです。
これらタイトル(正式には「次元名」というんでしょうか)を得るには、次のようにします。
dnn <- names(dimnames(xx$tbl1))
こうすると、dnnに「c("性別", "意見")」が代入されます。
sa1() の引数指定で、第３引数がなければ、１つの項目のみに着目した表が作成されます。
xx <- sa1(dtf, 2)
とした場合は、「性別」だけに着目した次のような表になります。tbl1, tbl2 の２つを示します。
−−−−
$tbl1
男 女
41 50
$tbl2
男 女 総数
41 50 91
−−−−
上のような１つの項目だけに着目した表を作った場合、横方向には複数の列ができますが、縦方向は１行のみです。
したがって、縦軸に沿ったパーセント(pct2)は意味がないのでNULLになります。
!!(2) 単一回答の集計に関する関数 sa2, sa3
先にsa1() について述べましたが、sa2(), sa3() も同じように用います。
呼び出す時の引数は２つか３つで、その戻り値は、tbl1など４つのテーブルからなるリストです。
これら sa1, sa2, sa3 の違いを下に列記します。
* sa1: 無回答を含まない表の作成
* sa2: 無回答を含む表の作成
* sa3: 有効回答(一種の小計になっている)を含む表の作成
以下に、sa1〜sa3の各々の出力例を示します。合計欄の「総数」が付いた表を掲げます。
欠損値(無回答)は、「無効回答」と表記されます。
−−−−
＊ sa1の出力例
|| ||賛成||中立||反対||総数
||男||14||12||11||37
||女||9||17||20||46
||総数||23||29||31||83
＊ sa2の出力例
|| ||賛成||中立||反対||無効回答||総数
||男||14||12||11||4||41
||女||9||17||20||4||50
||無効回答||2||3||2||0||7
||総数||25||32||33||8||98
＊ sa3の出力例
|| ||賛成||中立||反対||有効回答||無効回答||総数
||男||14||12||11||37||4||41
||女||9||17||20||46||4||50
||有効回答||23||29||31||83||8||91
||無効回答||2||3||2||7||0||7
||総数||25||32||33||90||8||98
−−−−
sa3() の戻り値は、tbl1に「総数」があるので、tbl2の方は NULL になります。この点、sa1, sa2 と少し違うので留意して下さい。
それから、sa3によって得られる表には「有効回答」という一種の小計欄があるので、カイ２乗検定等の統計検定の材料にはなりません。
!!(3) 単一回答集計のサンプルスクリプト
data01.csvを素材データとして、「性別」と「意見」のクロス集計を行うrubyスクリプトを掲げます。
Rプログラムの出力結果をそのまま標準出力に出力する他に、集計結果をrubyで受け取って、csvファイルとして書き出します。
−−−− tjs31.rb ここから
#! ruby -Ks
# 単一回答の集計(２項目のクロス集計) (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
source("tjsf.r")
dtf <- read.csv("data01.csv", header=T, na.strings="")
dtf$性別 <- factor(dtf$性別, levels=c("男", "女"), labels=c("男", "女"))
xx <- sa1(dtf, "性別", "意見")
for (i in 1:length(xx)) if (!is.null(xx[[i]])) xx[[i]] <- round(xx[[i]],1)
xx
xx$dnn <- names(dimnames(xx$tbl1))
robj(xx)
EOS
##
hs = Rrx.rexec(rpro)
puts "** Rプログラムの出力"
print hs[:sink]
printf("\n")
xx = hs["robj1"]
dnn = xx.delete("dnn") # 項目タイトル
dnnstr = (dnn.class == Array) ? dnn.join("/") : dnn
keys = %w(tbl1 tbl2 pct1 pct2)
titles = %w(基本表 「計」付きの表 横軸パーセント 縦軸パーセント)
str = ""
while key = keys.shift
title = titles.shift
val = xx[key]
next unless val
if dnn.class == Array # クロス集計表の場合
val[0][0] = dnnstr
else
val.delete_at(0)
val = Rrx.table_turn(val)
val.delete_at(0)
end
str = str + "** #{title}\n"
str = str + Rrx.ary2str(val, ",") + "\n"
end
File.open("result.csv", "w") {|ff| ff.write str}
puts "'result.csv' を書き出しました."
−−−− tjs31.rb ここまで
上は、２つの項目に着目したクロス集計(いわば２次元)の例ですが、１つの項目のみに着目した集計(１次元)にもそのまま使えます。
「xx$dnn <- names(dimnames(xx$tbl1))」というのは、縦・横の軸のタイトル名(次元名)を取得して、ruby側でそれを参照できるようにするためのものです。
次元名が２つあれば２次元(クロス集計)であり、１つだけなら１次元(１行のみの表)ということになります。
作られた表が２次元だとマトリクス(matrix)であり、１次元の場合は配列(array)です。それをruby側で受け取ると、どちらも配列(rubyのArray)の形になっています。しかし、こまったことに、単純に２次元配列｜１次元配列として識別できるわけではありません。
Rの配列(１次元)は、ruby側では２次元配列にになります。「配列の要素番号・項目名・度数」が複数入っている２次元配列です。Rのマトリクスもrubyでは２次元配列なので、その区別がつきにくい。
そこで、dnnに次元名をセットしてrubyに引きわたしています。ruby側において、dnnが配列なら２次元であり、単なる文字列であれば１次元と判断できます。
少々くどくなりましたが、ともあれ tjs31.rb は、２次元の表にも１次元の表にも対応します。ただし、３次元以上の表は作れません。
!!(4) 補完的関数
tjsf.rには補完的な役割の関数も入っているので、いくつか紹介します。
以下では、素材データとして data02.csv, data02.mem を用います。
data02.csvの中身は data01.csv と同じですが、列名が「性別→sex」、「意見→view」と英語名になっています。また、性別欄は「男→m」、「女→f」。意見欄は「賛成:1, 中立:2, 反対:3」の対応関係で数字が書き込まれています。
data02.mem は、その対応関係をメモ書きしたもので、中身は次の２行です。
sex | m,f | 男,女
view | 1,2,3 | 賛成,中立,反対
半角１文字の '|' は区切り文字です。タブコードも区切り文字になります。区切り文字の前後にある半角スペースは無視されます。
区切り文字によって３つの部分に分けられていますが、各々の部分が何を意味するかは見当がつくと思います。最初が列名、あとの２つは、カンマ区切り形式のlevelsとlabelsです。
このような素材データ(data02.csv, data02.mem)を基にして前掲のスクリプトと同じ結果を得るための tjs32.rb を下に掲げます。
−−−− tjs32.rb ここから
#! ruby -Ks
# 単一回答の集計(補完的関数の利用) (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
source("tjsf.r")
dtf <- read.csv("data02.csv", header=T, na.strings="")
dtf <- remake.factor(dtf, "data02.mem")
xx <- sa2(dtf, 2, 3)
xx <- rename.dnn(xx, "性別", "意見")
xx <- lround(xx, 1)
xx
xx$dnn <- get.dnn(xx)
robj(xx)
EOS
##
hs = Rrx.rexec(rpro)
puts "** Rプログラムの出力"
print hs[:sink]
printf("\n")
xx = hs["robj1"]
dnn = xx.delete("dnn") # 項目タイトル
dnnstr = (dnn.class == Array) ? dnn.join("/") : dnn
keys = %w(tbl1 tbl2 pct1 pct2)
titles = %w(基本表 「計」付きの表 横軸パーセント 縦軸パーセント)
str = ""
while key = keys.shift
title = titles.shift
val = xx[key]
next unless val
if dnn.class == Array # クロス集計表の場合
val[0][0] = dnnstr
else
val.delete_at(0)
val = Rrx.table_turn(val)
val.delete_at(0)
end
str = str + "** #{title}\n"
str = str + Rrx.ary2str(val, ",") + "\n"
end
File.open("result.csv", "w") {|ff| ff.write str}
puts "'result.csv' を書き出しました."
−−−− tjs32.rb ここまで
上のスクリプトに出てくる補完的な関数について、その仕様を記します。
◇ remake.factor
・利用例：dtf <- remake.factor(dtf, "data02.mem")
・機能：メモファイルに基づいてデータフレームの列をfactorとして再設定します。{{br}}
「remake.factor(dtf, "data02.mem")」は、data02.memの内容に即して、データフレームdtfの各列をfactorとして再設定するものです。{{br}}
この関数を用いずに、あえて別の形で記述するなら次の２行になります。{{br}}
dtf$sex <- factor(dtf$sex, levels=c("m","f"), labels=c("男","女"))
dtf$view <- factor(dtf$view, levels=1:3, labels=c("賛成","中立","反対"))
・戻り値：再設定された結果(データフレーム)を返します。指定されたメモファイルが存在しないなどの理由で再設定できない場合は、再設定せずに(つまり素のままで)データフレームを返します。
・引数：第１はデータフレーム、第２は再設定用メモ書きファイルの名前です。{{br}}
第３引数として、読み込むファイルの文字コードを指定することができます。data02.memがutf-8で書かれている場合は次のようにします。{{br}}
dtf <- remake.factor(dtf, "data02.mem", "utf-8"){{br}}
第３引数を省略すると、Rのデフォルトの文字コードになります。日本語版のWindowsだと "cp932" です。
◇ rename.dnn
・利用例：xx <- rename.dnn(xx, "性別", "意見")
・機能：テーブルの次元名を変構します。{{br}}
サンプルスクリプトにおいて、xx(リスト)には tbl1, tbl2, pct1, pct2 の４つのテーブルが入っているわけですが、各々のテーブルの次元名は、「sex, view」になっています。これを「性別, 意見」に変更するのが rename.dnn() です。４つのテーブルすべてについて変更を行います。
・戻り値：第１引数xxの次元名を変更した結果を返します。変更できない時は、xxをそのまま返します。
・引数：第１引数は、テーブルからなるリストを指定します。あるいは、テーブル１つを指定することもできます。{{br}}
第２引数と第３引数は次元名の指定です。第２が行名、第３が列名。２次元なら「"性別", "意見"」のように２つを指定します。１次元であれば "性別" のように１つだけを指定します(第３引数を省略します)。
◇ lround
・利用例：xx <- lround(xx, 1)
・機能：リストの各構成要素に対して、まるめ処理を行います。{{br}}
リストの構成要素として、テーブル・マトリクス・ベクトルを想定しますが、それ以外が含まれていても支障ありません。{{br}}
リストの構成要素が実数値(double)か複素数値(complex)でなければ、まるめ処理を施しません。
・戻り値：まるめ処理を施した結果を返します。
・引数：第１引数はリストを想定していますが、テーブル・マトリクス・ベクトルなどでもかまいません。{{br}}
第２引数は、まるめ処理に係る小数点の桁数を指定するものです。デフォルトは0です。つまり整数値になるようにまるめます。これを1にすると、小数点１桁までの数値になります。
◇ get.dnn
・利用例：dnn <- get.dnn(xx)
・機能：引数xxの次元名を取得します。{{br}}
xxがリストであれば、その中の最初にみつかったテーブルの次元名を返します。xxがテーブルなら、その次元名を返します。{{br}}
sa1()〜sa3()の戻り値(リスト)から次元名を得ることを目的として、この get.dnn() を設けました。
・戻り値：引数xxの次元名。例えば c("sex", "view") など。
・引数：リスト、またはテーブル
以上が補完的関数の説明です。
余談ですが、tjs31.rb, tjs32.rb のどちらも、小数点以下の桁数についてのまるめ処理をRプログラムの側で行っています。各々のテーブルに対して round() を適用しています。
この処理をR側ではなくruby側で行うサンプルスクリプトを tjs32b.rb として同梱しておきます。
参考まで、rubyスクリプトの該当部分を下に掲げます。
−−−− tjs32b.rb 抜粋ここから
str = ""
while key = keys.shift
title = titles.shift
val = xx[key]
next unless val
if dnn.class == Array # クロス集計表の場合
val[0][0] = dnnstr
else
val.delete_at(0)
val = Rrx.table_turn(val)
val.delete_at(0)
end
if key =~ /^pct[12]$/
val.map! {|row|
row.map! {|cell|
[Fixnum,Float].include?(cell.class) ? sprintf("%.1f",cell) : cell
}
}
end
str = str + "** #{title}\n"
str = str + Rrx.ary2str(val, ",") + "\n"
end
−−−− tjs32b.rb 抜粋ここまで
------------------------------------------------------------------------
!2. 複数回答の集計
ここでは、「主な情報源」として、新聞・雑誌・テレビ・ラジオの４つの選択肢から、当てはまるものを複数個選んでもらうような複数回答の集計を取り上げます。
素材データとして data03.csv を用います。これは、ID(整数値), 年齢層(20代・30代・40代のどれか１つ), 主な情報源(新聞・雑誌・テレビ・ラジオの複数回答)から構成されています。98人分のデータです。具体的には次のようになっています。
ID,年齢層,新聞,雑誌,テレビ,ラジオ
1,30代,1,,,
2,40代,1,1,,
3,40代,,1,,
年齢層の列は単一回答です。
主な情報源の列は、４項目のうち選択された項目のところに数字の1を入力し、選択されなかったところは空欄にしてあります。
data03.csv は、乱数発生によって作成したデータです。実態とは関係ありません。
複数回答にまつわる主な関数は、ma1(), ma2() の２つです。他に ma.refactor(), ma.count() があります。
!!(1) 複数回答と単一回答の組合せ集計 ma1
まず、年齢層(単一回答)と主な情報源(複数回答)の組合せ集計を行います。
ma1()という関数をtjsf.rに盛り込んだので、それを用います。そのサンプル(Rプログラム部分のみ)を掲げます。
−−−− Rプログラムここから
source("tjsf.r")
dtf <- read.csv("data03.csv", header=T, na.strings="")
xx <- ma1(dtf, 2, 3:6)
xx
−−−− Rプログラムここまで
ma1() の戻り値は、tbl, pct の２つの表からなるリストです。前者は度数(選択者数)、後者は比率(パーセント)の表です。
上のプログラムの実行結果は次のとおり。
−−−−
$tbl
新聞 雑誌 テレビ ラジオ 該当数
20代 13 17 14 14 31
30代 24 17 20 16 34
40代 15 22 17 16 33
合計 52 56 51 46 98
$pct
新聞 雑誌 テレビ ラジオ 該当数
20代 41.93548 54.83871 45.16129 45.16129 100
30代 70.58824 50.00000 58.82353 47.05882 100
40代 45.45455 66.66667 51.51515 48.48485 100
合計 53.06122 57.14286 52.04082 46.93878 100
−−−−
ma1() を呼び出す場合、データフレームの列を番号でなく名前で指定することができます。次のように書くことができます。
ma1(dtf, "年齢層", c("新聞","雑誌","テレビ","ラジオ"))
第１引数がデータフレーム、第２引数と第３引数が列の指定です。第２・第３は、順序を入れ換えても差し支えありません。
この ma1() は、次のサイトにある ma() を基にして作りました(というより、変更を加えました)。作者に感謝します。
[[R -- マルチアンサーのクロス集計|http://aoki2.si.gunma-u.ac.jp/R/multianswer.html]]
!!(2) 複数回答の列をfactorとして再設定 ma.refactor
複数回答の列の値について少し触れます。
data03.csvの場合、選択されたら数値の1、そうでなければ空欄になっています。空欄は、read.csvで読み込んだ後に欠損値のNAとなります。
もし空欄でなく数字の0が入力されていたとすると、前述のRプログラムは、意図したような結果を出してくれません。
数字の1と0を機械的に列べると、当然ながら0の方が先にきます。そうすると、ma1() は、0の方を「選択あり」であるかのように処理してしまいます。
前もって、0よりも1の方が先にくるものであることをRに知らせてやる必要があります。そのためには、複数回答の列をfactorとして再設定します。
例えば、0と1からなる３列目を再設定するには次のようにします。
dtf[,3] <- factor(dtf[,3], levels=1:0, labels=1:0)
上の記述を３列目だけでなく３〜６列すべてに対して実行すればいいわけです。
これを少し簡単に記述できるようにするため、ma.refactor() という関数を設けました。次のように呼び出します。
dtf <- ma.refactor(dtf, 3:6, 1:0, 1:0)
こうすると、３〜６列のいずれも factor として再設定されます。
第１引数はデータフレーム、第２が列の指定(番号でなく名前でも可)、第３がlevels、第４はlabelsです。
第４引数を省略すると、「c("新聞○", "新聞×")」が指定されたものとみなされます。「新聞」のところは各々の列名に置き換えられます。
data03.csvの場合でいうと、数字の0ではなく空欄なので、再設定は次のようにします。
dtf <- ma.refactor(dtf, 3:6, c(1,NA), 1:0)
ここで、ma1(), ma.refactor() の利用例を掲げておきます。
Rプログラムの出力を標準出力に出力するほか、ma1() の戻り値をrubyの側で受け取って、それをcsvファイルとして書き出します。
−−−− tjs33.rb ここから
#! ruby -Ks
# 複数回答と単一回答の組合せ集計 (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
source("tjsf.r")
dtf <- read.csv("data03.csv", header=T, na.strings="")
dtf <- ma.refactor(dtf, 3:6, c(1,NA))
xx <- ma1(dtf, 2, 3:6)
xx <- lround(xx, 1)
xx
robj(xx)
EOS
##
hs = Rrx.rexec(rpro)
puts "** Rプログラムの出力"
print hs[:sink]
printf("\n")
xx = hs["robj1"]
keys = %w(tbl pct)
titles = %w(回答の度数分布 回答比率)
str = ""
while key = keys.shift
title = titles.shift
val = xx[key]
next unless val
str = str + "** #{title}\n"
str = str + Rrx.ary2str(val, ",") + "\n"
end
File.open("result.csv", "w") {|ff| ff.write str}
puts "'result.csv' を書き出しました."
−−−− tjs33.rb ここまで
!!(3) 複数回答と単一回答の分解的集計
ma1() は、新聞・雑誌・テレビ・ラジオの４項目からなる複数回答の状況を一括して集計するので、全体像をみるのに向いています。ただ、その集計結果をそのまま統計検定にかけることはできません。
そこで、「新聞」にまるを付けたか否か、「雑誌」にまるを付けたか否かというように、各々の選択肢について二者択一でとらえ、集計表を４つ作ります。そうすると、各々の表をカイ２乗検定等にかけることができます。
前述したように、ma.refactor() の第４引数(labelsの指定)を省略すると、「新聞○」と「新聞×」のように設定されるので、あとは単一回答の集計を行う sa1() を適用してやれば、意図した表を作れます。
下に、それら４つの表を作成するスクリプトを掲げます。
単に表を作るだけでは何なので、カイ２乗検定にかけて、その結果を出力してみます。有意性が認められた時は、調整済み残差の表も出力します。
Rプログラムの実行結果は標準出力に出力し、ruby側では同じ内容をcsvファイルに書き出しています。
まだ紹介していないRの関数が出てきますが、後で説明します。
−−−− tjs34.rb ここから
#! ruby -Ks
# 複数回答と単一回答の分解的集計 (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
source("tjsf.r")
dtf <- read.csv("data03.csv", header=T, na.strings="")
dtf <- ma.refactor(dtf, 3:6, c(1,NA))
items <- c("χ2", "df", "p-value", "cramer")
keys <- c("statistic", "parameter", "p.value", "cramer")
for (i in 3:6) {
xx <- sa1(dtf, 2, i)
zz <- chisq.test(xx$tbl1)
zz <- append.cramer(zz)
zz$join <- join.tbl(xx$tbl2, xx$pct1, "%d(%.1f)")
ww <- sapply(keys, function(k) zz[[k]])
names(ww) <- items
print(zz$join); cat("\n")
print(ww); cat("\n")
if (zz$p.value < 0.05) {
cat("☆ 調整済み残差\n")
print(zz$stdres); cat("\n")
}
cat("\n\n")
robj(zz)
}
EOS
##
hs = Rrx.rexec(rpro)
puts "** Rプログラムの出力"
print hs[:sink]
printf("\n")
items = %w(χ2 df p-value cramer)
keys = %w(statistic parameter p.value cramer)
str = ""
Rrx.robj_each() do |kk, zz|
str = str + Rrx.ary2str(zz["join"], ",") + "\n"
ary = []
keys.each {|key| ary << sprintf("%g", zz[key])}
str = str + Rrx.ary2str([items, ary], ",") + "\n"
if zz["p.value"] < 0.05
str = str + "☆ 調整済み残差\n"
str = str + Rrx.ary2str(zz["stdres"], ",") + "\n"
end
str = str + "\n\n"
end
File.open("result.csv", "w") {|ff| ff.write str}
puts "'result.csv' を書き出しました."
−−−− tjs34.rb ここまで
参考まで、上のスクリプトを実行した結果の一部を掲げておきます。
「年齢層」と「新聞」のクロス集計表で、１つのセルが「13(41.9)」のように、度数と括弧書きのパーセントという形になっています。
クロス集計表の次は、カイ２乗検定の結果(カイ２乗値・自由度・有意確率・クラメール係数)を列記したものです。
−−−−
新聞
年齢層 新聞○ 新聞× 総数
20代 13(41.9) 18(58.1) 31(100.0)
30代 24(70.6) 10(29.4) 34(100.0)
40代 15(45.5) 18(54.5) 33(100.0)
総数 52(53.1) 46(46.9) 98(100.0)
χ2 df p-value cramer
6.50090607 2.00000000 0.03875665 0.25755733
−−−−
余談ですが、data03.csvは乱数発生によって作ったものなので、カイ２乗検定で有意性が検出されることはないだろうと予想しましたが、上の「新聞」については有意性が認められました。有意性が検出されたのはこれだけです。
調整済み残差をみるまでもなく、「30代」の「新聞○」が多く、「新聞×」が少ない様子が分かります。
tjs34.rbで新たに出てきた関数について説明します。２つあります。
append.cramer(zz) は、引数zz(リスト)にクラメール係数を追加して返すものです。zzには、予めカイ２乗検定の結果を代入しておく必要があります。つまり chisq.test() の戻り値を代入しておきます。
「zz <- append.cramer(zz)」としたとき、zz$cramer にクラメール係数の値がセットされています。
もう１つの join.tbl() は、２つの表を結合するものです。典型的な例として、度数分布の表と構成比の表の結合があります。上のサンプルスクリプトもその一例です。
tt <- join.tbl(tbl, pct, "%d(%.1f)")
とすると、tbl, pct の２つの表が結合されて、その結果が tt に代入されます。ここで結合というのは、対応するセル１つづつを組み合わせるという意味です。
第３引数は、結合セルの書式を指定するものです。これが内部処理においてsprintfのフォーマットとして使われます。"%d(%.1f)" の場合、tblのセルが整数値として書き出され、pctのセルが括弧書きの小数点数(有効桁数１桁)で書き出されます。
もし pct が表として tbl よりも小さくて、行数または列数が少なければ、その少ない部分については tbl のセルだけが書き出されます。その場合、第３引数の書式は用いません。
rubyスクリプトの部分については説明を省略します。
拙作rrxライブラリのメソッドをいくつか使っているので分かりにくいと思いますが、必要に応じ関連のマニュアルを参照して下さい。
◇ append.cramer
・利用例：zz <- chisq.test(xx); zz <- append.cramer(zz)
・機能：引数zzにクラメール係数を追加して返します。zzには予めカイ２乗検定の結果を代入しておきます。
・戻り値：クラメール係数を含むカイ２乗検定結果を返します。戻り値zzのzz$cramerがクラメール係数の値です。
・引数：カイ２乗検定結果を引数として与えます。
◇ join.tbl
・利用例：tt <- join.tbl(tbl, pct, "%d(%.1f)")
・機能：２つの表を結合します。対応するセルを組み合わせます。
・戻り値：結合した表を返します。
・引数：第１引数と第２引数は表データ、第３引数は結合セルの書式で、内部処理においてsprintfのフォーマットとして使われます。省略時は "%s(%s)" とみなされます。
!!(4) 複数回答の選択肢相互の関係をみる ma2
「主な情報源」の４つの選択肢(新聞・雑誌・テレビ・ラジオ)の相互の関係をみるには、「新聞」と「雑誌」の両方にまるを付けた人がどれくらいいるか、「新聞」と「テレビ」の両方にまるを付けた人はどうか、といった確認をすることになると思います。そうすると、例えば、「新聞」を選ぶ人は「ラジオ」も選ぶ傾向が高いといったことが分かるかもしれません。
この「両方にまるを付けた人がどれくらいか」を一覧表示するのが ma2() です。具体的には次のような表を生成します。
−−−−
新聞 雑誌 テレビ ラジオ 該当数
新聞 52 29 23 28 52
雑誌 29 56 30 25 56
テレビ 23 30 51 22 51
ラジオ 28 25 22 46 46
該当数 52 56 51 46 98
−−−−
上の表で、「該当数」というのは、いわば「総数」に当たるものです。「新聞」にまるを付けた人の総数が何人か、「雑誌」の総数は何人かなどが「該当数」をみれば分かります。「該当数」の最も右下の98は、回答者総数を意味します。
「新聞」と「雑誌」が交差するセルの値は29ですが、これは、この両方にまるを付けた人が29人いることを意味します。
この表は、縦軸と横軸を入れ換えても同じものになります。
ma2() は、このような表を返しますが、４種類の表を生成して返します。
ma2() の仕様について説明する前に、その利用例を下に掲げます。Rプログラム部分のみです。
−−−− Rプログラムここから
source("tjsf.r")
dtf <- read.csv("data03.csv", header=T, na.strings="")
dtf <- ma.refactor(dtf, 3:6, c(1,NA))
xx <- ma2(dtf, 3:6)
xx <- lround(xx, 1)
xx
−−−− Rプログラムここまで
まず ma2 を呼び出す時の引数ですが、第１引数がデータフレーム、第２引数は複数回答の列の指定です。列の指定は、列番号でも列名でも、どちらでもかまいません。
ma2() の戻り値は、次の４つの表を含むリストです。
* tbl1: 「該当数」を含まない表
* tbl2: 「該当数」を含む表
* pct1: 右端の列の「該当数」を100%とする時のパーセント
* pct2: 右下端の「回答者総数」を100%とする時のパーセント
先に掲げた表は、「該当数」付きの tbl2 です。tbl1 は、それから縦と横の「該当数」を取り除いたものです。
pct1, pct2 を下に掲げておきます。
−−−−
$pct1
新聞 雑誌 テレビ ラジオ 該当数
新聞 100.0 55.8 44.2 53.8 100
雑誌 51.8 100.0 53.6 44.6 100
テレビ 45.1 58.8 100.0 43.1 100
ラジオ 60.9 54.3 47.8 100.0 100
該当数 53.1 57.1 52.0 46.9 100
$pct2
新聞 雑誌 テレビ ラジオ 該当数
新聞 53.1 29.6 23.5 28.6 53.1
雑誌 29.6 57.1 30.6 25.5 57.1
テレビ 23.5 30.6 52.0 22.4 52.0
ラジオ 28.6 25.5 22.4 46.9 46.9
該当数 53.1 57.1 52.0 46.9 100.0
−−−−
pct1 は、「新聞」にまるを付けた人の総数を100%とした時に、「雑誌」にもまるを付けた人が何％か、「テレビ」にもまるを付けた人が何％かなどをみるのに使えます。
pct2 の方は、回答者総数98人を100%とした時に、各セルの値が何％かを確認するのに使います。全体の中で「両方にまるを付けた人」が何％いるかが分かります。
４つの表のうち pct1 以外は、縦軸と横軸を入れ換えても同じものになります。pct1 は同じになりません。
こうした表は、きちんと分析するのにはあまり適していませんが、ざっと全体像をつかむ時に、それなりに使えるのではないかと思います。
下に ma2() の利用サンプルを掲げておきます。
−−−− tjs35.rb ここから
#! ruby -Ks
# 複数回答の相互関係をみる (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
source("tjsf.r")
dtf <- read.csv("data03.csv", header=T, na.strings="")
dtf <- ma.refactor(dtf, 3:6, c(1,NA))
xx <- ma2(dtf, 3:6)
xx <- lround(xx, 1)
xx
robj(xx)
EOS
##
hs = Rrx.rexec(rpro)
puts "** Rプログラムの出力"
print hs[:sink]
printf("\n")
xx = hs["robj1"]
keys = %w(tbl1 tbl2 pct1 pct2)
titles = %w(基本表 「該当数」付きの表 横軸パーセント 対総数のパーセント)
str = ""
while key = keys.shift
title = titles.shift
val = xx[key]
next unless val
str = str + "** #{title}\n"
str = str + Rrx.ary2str(val, ",") + "\n"
end
File.open("result.csv", "w") {|ff| ff.write str}
puts "'result.csv' を書き出しました."
−−−− tjs35.rb ここまで
!!(5) 複数回答の選択肢相互の関係を分解的にみる
前述の ma2() が生成する表は、統計検定にかけることができません。そこで、各選択肢をばらばらにして集計してみます。
「新聞」と「雑誌」の２つに注目した時の表は次のとおり。
−−−−
|| ||雑誌○||雑誌×
||新聞○||29||23
||新聞×||27||19
−−−−
選択肢が４つなので、そこから２つを抜き出して表を作るやり方は６通りあります。
上のような２行・２列の表を６つ生成して、各々を統計検定にかけるためのサンプルを下に掲げます。カイ２乗検定を用います。
−−−− tjs36.rb ここから
#! ruby -Ks
# 複数回答の相互関係を分解的にみる (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
source("tjsf.r")
dtf <- read.csv("data03.csv", header=T, na.strings="")
dtf <- ma.refactor(dtf, 3:6, c(1,NA))
items <- c("χ2", "df", "p-value", "cramer")
keys <- c("statistic", "parameter", "p.value", "cramer")
for (i in 3:5) {
for (j in (i+1):6) {
xx <- sa1(dtf, i, j)
zz <- chisq.test(xx$tbl1)
zz <- append.cramer(zz)
zz$join <- join.tbl(xx$tbl2, xx$pct1, "%d(%.1f)")
ww <- sapply(keys, function(k) zz[[k]])
names(ww) <- items
print(zz$join); cat("\n")
print(ww); cat("\n")
if (zz$p.value < 0.05) {
cat("☆ 調整済み残差\n")
print(zz$stdres); cat("\n")
}
cat("\n\n")
robj(zz)
}
}
EOS
##
hs = Rrx.rexec(rpro)
puts "** Rプログラムの出力"
print hs[:sink]
printf("\n")
items = %w(χ2 df p-value cramer)
keys = %w(statistic parameter p.value cramer)
str = ""
Rrx.robj_each() do |kk, zz|
str = str + Rrx.ary2str(zz["join"], ",") + "\n"
ary = []
keys.each {|key| ary << sprintf("%g", zz[key])}
str = str + Rrx.ary2str([items, ary], ",") + "\n"
if zz["p.value"] < 0.05
str = str + "☆ 調整済み残差\n"
str = str + Rrx.ary2str(zz["stdres"], ",") + "\n"
end
str = str + "\n\n"
end
File.open("result.csv", "w") {|ff| ff.write str}
puts "'result.csv' を書き出しました."
−−−− tjs36.rb ここまで
複数回答の答えは「○」か「×」かの２通りなので、各々の表は２行・２列になります。なので、カイ２乗検定の他に、フィッシャーの直接確率検定を適用することができます。
４つのセルの中に５未満のものがあるなど度数が少数のケースでは、カイ２乗検定を行った時に、警告メッセージが出ることがあります。そうした時はフィッシャーの検定を使うのがいいかもしれません。
zz <- fisher.test(xx$tbl1)
print(zz$p.value)
このようにすると、p値が表示されます。これが 0.05 未満であれば、有意性が認められることになります。
有意性が認められれば、例えば、「新聞を主な情報源にする人は、ラジオも情報源にする傾向がある。」といったことが言えるかもしれません。しかし、有意性が認められないのであれば、少なくとも統計検定の結果からは、そうした傾向を指摘する根拠が得られない、ということになります。ちなみに、６つの表の中で有意性が認められるものはありませんでした。
このような分解的集計・検定の他に、数量化IIIや対応分析を用いることもできます。そのやり方について「単純集計に関する覚え書き／複数回答の集計」でごく簡単に触れましたが、関連の書籍やWebがいろいろあるので参考にして下さい。
!!(6) 複数回答の選択個数をみる ma.count
複数回答の選択肢について、回答者各人がまるを何個つけたかを検出しておくと便利です。
IDが１番の人は２個、２番の人は０個、３番の人は３個などのように、各人の選択個数を把握します。そうすれば、選択個数０の人(いわば無回答の人)を除いて集計したりすることが容易になります。
この選択個数の検出を行う関数 ma.count() を設けました。
これを呼び出す時の引数は４つありますが、第１引数がデータフレーム、第２が複数回答の列番号(または列名)です。多くの場合、第３と第４は省略可能です。
とりあえず ma.count の利用例を掲げます。Rプログラム部分のみです。
−−−− Rプログラムここから
source("tjsf.r")
dtf <- read.csv("data03.csv", header=T, na.strings="")
count <- ma.count(dtf, 3:6)
table(count)
−−−− Rプログラムここまで
上のプログラムを実行すると、下のような表示が得られます。
−−−−
count
0 1 2 3 4
4 24 34 31 5
−−−−
選択個数０の人が４人、選択個数１個の人が24人いることが分かります。
ma.count() の戻り値はベクトルです。dtfには98人分のデータが入っているので、戻り値である count にも98人分の選択個数が入ります。
前述の例では「ma.count(dtf, 3:6)」のように、引数を２つしか指定しませんでしたが、第３引数と第４引数を省略せずに指定すると、次のようになります。
count <- ma.count(dtf, 3:6, 1, FALSE)
第３引数の１は、素材データであるdtfにおいて、複数回答の「選択あり」を示す入力値です。
dtfでは、選択された時に数字の１、選択されなかった時は空欄になっています。そこで、選択個数をカウントする際に数字の１日着目するというのが第３引数の意味です。
dtf <- ma.refactor(dtf, 3:6, c(1,NA), c("yes","no"))
上のように複数回答の列をfactorとして再設定していると、「選択あり」を示すのは "yes" になるので、ma.count() の呼び出し方は次のようになります。
count <- ma.count(dtf, 3:6, "yes", FALSE)
第３引数が "yes" のように文字である場合、実は、それを正規表現とみなします。内部処理では、Rのgrep()によってマッチングを行います。なので、厳密に「yes」の３文字を示すなら "^yes$" とすべきかもしれません。しかし、もう１つの値が "no" なので、そこまでしなくても意図どおりの結果が得られます。
第４引数は、選択個数を素材データのdtfに挿入・追加するか否かの指定です。デフォルトは FALSE です。つまり、挿入・追加をせず、選択個数を記録したベクトルを戻り値として返します。
TRUEにすると、dtfに nnn という名前の列を設けて、それに選択個数を記録します。戻り値として、nnnを追加・挿入したデータフレームを返します。
戻り値がdtf2に代入されている場合、dtf2$nnn または dtf2[,"nnn"] によって選択個数を参照できます。既に素材データのdtfに nnn という名前の列がある時は、そのnnnを上書きします。
選択個数０の無回答の人を除いたデータフレームを得たい時は、次のようにします。
dtf2 <- ma.count(dtf, 3:6, 1, TRUE)
dtf2 <- subset(dtf2, nnn > 0)
上のようにすると、無回答の人を除く94人分のデータが dtf2 にセットされます。
以上が ma.count() の仕様です。
nnn という列名には違和感があるかもしれませんが、count だと、よく使われるような感じがして、競合を避けるため nnn にしました。
!!(7) 複数回答にまつわる関数の一覧
複数回答にまつわる関数の一覧を掲げておきます。
◇ ma1
・利用例：xx <- ma1(dtf, 2, 3:6)
・機能：複数回答と単一回答の組合せ集計を行います。
・戻り値：tbl, pct の２つの表からなるリストを返します。前者は度数、後者はパーセント。
・引数：第１引数はデータフレーム、第２と第３は列の指定です。どちらかが複数回答の列、もう一方は単一回答の列です。その順番は問いません。
◇ ma2
・利用例：xx <- ma2(dtf, 3:6, 1, FALSE)
・機能：複数回答の選択肢相互の関係をみるための表を生成します。
・戻り値：tbl1, tbl2, pct1, pct2 の４つの表からなるリストを返します。４つの表が何かは、本文を参照して下さい。
・引数：第１引数はデータフレーム、第２引数は複数回答の列の指定です。
◇ ma.refactor
・利用例：dtf <- ma.refactor(dtf, 3:6, c(1,NA), 1:0)
・機能：複数回答の列をfactorとして再設定します。
・戻り値：再設定された結果(データフレーム)を返します。
・引数：第１引数はデータフレーム、第２引数は列の指定、第３はlevels、第４はlabelsです。{{br}}
第４引数を省略すると、「c("新聞○", "新聞×")」が指定されたものとみなされます。「新聞」のところは各々の列名に置き換えられます。
◇ ma.count
・利用例：count <- ma.count(dtf, 3:6, 1, FALSE){{br}}
dtf2 <- ma.count(dtf, 3:6, "○", TRUE)
・機能：複数回答の選択個数を検出します。
・戻り値：第４引数がFALSEの時は、選択個数を記録したベクトルを返し、TRUEの時は、選択個数を記録した列(nnn)が追加・挿入されたデータフレームを返します。
・引数：第１はデータフレーム、第２は複数回答の列の指定、第３が「選択あり」を示す値(文字だと正規表現パターン)、第４は列の追加・挿入の有無です。
------------------------------------------------------------------------
!3. 数値データの集計
ここでは、数値データとして「年齢」と「支出」(月当たりの支出額)を含むcsvデータを材料にします。
data04.csvは、ID(整数値)、性別(m,f)、年齢(数値)、支出(数値)から構成されています。98人分のデータです。具体的には次のようになっています。
ID,性別,年齢,支出
1,m,29,82249
2,m,49,93181
3,f,42,93647
数値データの集計にまつわる主な関数は、va1(), va2() の２つです。どちらも、単一回答と数値データの組合せの集計に用います。単一回答の「性別」(男女別)で分類しながら、数値データの平均値等を算出します。
他に date.split(), var2(), sd2() をtjsf.rに盛り込んであります。
date.split() は、「単純集計に関する覚え書き／数値データの集計」で触れたので今回は言及しません。仕様は同じです。
var2(), sd2() は、それぞれ標本分散と標本標準偏差を返す関数です。欠損値NAを取り除いた上で値を算出します。わたす引数は１つのみです。
ちなみに、Rに最初から組み込まれている var(), sd() は、それぞれ不偏分散と不偏標準偏差を返す関数です。こちらは引数が１つとはかぎりません。
!!(1) １つの値のみを返す関数を用いた集計 va1
mean() は平均値を返す関数ですが、このように１つの値のみを返す関数によって集計を行う際、va1() を使えます。
data04.csvを素材にして、男女別に平均年齢を算出し、それを表として得たい場合、次のようにします。Rプログラム部分のみ示します。
−−−− Rプログラムここから
source("tjsf.r")
dtf <- read.csv("data04.csv", header=T, na.strings="")
dtf$性別 <- factor(dtf$性別, levels=c("m","f"), labels=c("男","女"))
xx <- va1(dtf, "性別", "年齢", mean)
print(xx)
−−−− Rプログラムここまで
上のプログラムを実行した結果は下のとおりです。va1() は matrix を返します。
−−−−
年齢
男 34.37500
女 34.04762
全体 34.23469
−−−−
va1() の第１引数は、素材となるデータフレームです。
第２引数は、分類の手がかりとなる単一回答の列を指定します。名前でも番号でもどちらでもかまいません。
第３引数では数値データの列を指定します。「c("年齢","支出")」とか「3:4」などのように複数指定してもかまいません。ちなみに、第２引数の方は１つだけです。
最後の第４引数は、適用したい関数の名前を指定します。省略すると mean が指定されたものとみなされます。
関数は、既存のものでなく自作のものでもかまいません。例えば、１つのセルを「平均値(sd:標準偏差)」の書式で出力したい時は次のようにします。
xx <- va1(dtf, 2, 3:4, function(x)
sprintf("%.1f(sd:%.2f)", mean(x), sd(x)))
上のようにすると、下の結果が得られます。
−−−−
年齢 支出
男 "34.4(sd:8.23)" "105605.7(sd:29094.26)"
女 "34.0(sd:7.41)" "92032.0(sd:27015.33)"
全体 "34.2(sd:7.85)" "99788.4(sd:28879.60)"
−−−−
下にサンプルスクリプトを掲げておきます。平均値と中央値を別々に表示するものです。
Rプログラムの出力結果を標準出力に出力し、ruby側では同じ内容をcsvとして書き出します。
−−−− tjs37.rb ここから
#! ruby -Ks
# １つの値を返す関数による集計 (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
source("tjsf.r")
dtf <- read.csv("data04.csv", header=T, na.strings="")
dtf$性別 <- factor(dtf$性別, levels=c("m","f"), labels=c("男","女"))
xx1 <- va1(dtf, 2, 3:4, mean)
cat("平均値\n")
xx1
cat("\n")
xx2 <- va1(dtf, c("年齢","支出"), "性別", median)
cat("中央値\n")
xx2
robj(list(平均値=xx1, 中央値=xx2))
EOS
##
hs = Rrx.rexec(rpro)
puts "** Rプログラムの出力"
print hs[:sink]
printf("\n")
str = ""
Rrx.robj_each(3) do |name, key, val|
str = str + "** #{key}\n"
str = str + Rrx.ary2str(val, ",") + "\n\n"
end
File.open("result.csv", "w") {|ff| ff.write str}
puts "'result.csv' を書き出しました."
−−−− tjs37.rb ここまで
!!(2) 複数の値を返す関数を用いた集計 va2
va2() は、使い方は va1() と同じです。
第４引数で指定する関数が、１つの値だけでなく複数の値を返す場合に va2() を用います。
例えば、数値データに関して summary() は６つの値を返しますが、これを関数として指定したい時は次のようにします。
xx <- va2(dtf, 2, 3:4, summary)
実は、第４引数を省略すると、summary が指定されたものとみなされます。
va2() は、戻り値として list を返します。list の各要素は matrix です。
以下にサンプルスクリプトとその出力結果を示します。
−−−− tjs38.rb ここから
#! ruby -Ks
# 複数の値を返す関数を用いた集計 (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
source("tjsf.r")
dtf <- read.csv("data04.csv", header=T, na.strings="")
dtf$性別 <- factor(dtf$性別, levels=c("m","f"), labels=c("男","女"))
xx <- va2(dtf, 2, 3:4, summary)
xx
robj(xx)
EOS
##
hs = Rrx.rexec(rpro)
puts "** Rプログラムの出力"
print hs[:sink]
printf("\n")
str = ""
Rrx.robj_each(3) do |name, key, val|
str = str + "** #{key}\n"
str = str + Rrx.ary2str(val, ",") + "\n\n"
end
File.open("result.csv", "w") {|ff| ff.write str}
puts "'result.csv' を書き出しました."
−−−− tjs38.rb ここまで
−−−− 出力結果ここから
$年齢
Min. 1st Qu. Median Mean 3rd Qu. Max.
男 20 28.00 34.5 34.38 41.25 49
女 21 27.25 33.0 34.05 41.00 49
全体 20 28.00 34.0 34.23 41.00 49
$支出
Min. 1st Qu. Median Mean 3rd Qu. Max.
男 50430 82960 115200 105600 129600 149400
女 52490 69530 92120 92030 114900 148400
全体 50430 73590 98100 99790 124500 149400
−−−− 出力結果ここまで
数値データについては、例えば「年齢」を「20代・30代・40代」のようにカテゴリー化するのもよくあることだと思いますが、これを関数化するのは難しいのでやってません。
カテゴリー化の方法については「単純集計に関する覚え書き／数値データの集計」のところで触れているので、そちらを参考にして下さい。
------------------------------------------------------------------------
以上で「単純集計に関する覚え書き／まとめ(関数化)」をおわりにします。
tjsf.r の中には、上で紹介していない関数も含まれていますが、それらは内部処理で用いるために設けたサブルーチン的な関数です。興味があったら覗いてみて下さい。
Copyright (C) T. Yoshiizumi, 2013 All rights reserved.