単純集計に関する覚え書き/まとめ(関数化)
単純集計に関する覚え書き/まとめ(関数化)
〜 rubyとRを用いた処理 〜
最終更新日: 2013/08/11
以下に掲げるドキュメントは、
tjs_d.zip
に同梱されている tjs_d.txt と同じ内容です。
tjs_d.zipには、euc-jp | utf-8 で書かれたサンプルスクリプトも同梱されています。
サンプルスクリプトは、MS-Windowsだけでなく、linuxなどでも動作します(CentOS6で確認)。
<はじめに>
「単純集計に関する覚え書き」を3回にわたって記しました。予想していたよりも長い解説になってしまい、読み返すのが面倒です。
具体的な集計の必要に迫られたとき、どのサンプルスクリプトを利用したらいいのか、どんなふうに書き換えたらいいのか、私自身が迷うくらいにだらだらとした解説になってしまいました。
そこで、「最低限これくらいは」というものをRの関数として定義しました。例えば、単一回答の集計の一つ(無回答を除く集計)の sa1() は次のように呼び出します(saは single answer のつもり)。
xx <- sa1(dtf, 2, 3) print(xx)
こうすると、データフレームdtfの2列目と3列目に着目して、クロス集計表を作成します。
戻り値のxxには4つの表が入っています。2つの度数分布表と2つのパーセントの表です。
こうした sa1() のような関数をいくつか作って tjsf.r に入れました。以下で、それらの解説を記します。
なお、この tjsf.r は、rubyに依存する部分はありません。rubyを使わない場合でも利用可能です。
--------
- 1. 単一回答の集計
- 2. 複数回答の集計
- (1) 複数回答と単一回答の組合せ集計 ma1
- (2) 複数回答の列をfactorとして再設定 ma.refactor
- (3) 複数回答と単一回答の分解的集計
- (4) 複数回答の選択肢相互の関係をみる ma2
- (5) 複数回答の選択肢相互の関係を分解的にみる
- (6) 複数回答の選択個数をみる ma.count
- (7) 複数回答にまつわる関数の一覧
- 3. 数値データの集計
--------
1. 単一回答の集計
ここでは、「性別」(男・女)、「意見」(賛成・中立・反対)のように、選択肢の中から1つだけ選んでもらう単一回答の集計を取り上げます。
適切な選択肢がないため回答できない|回答したくないというケースもありますが、それは無回答という形になって表れると思います。そうした無回答を含めて集計する場合についても取り上げます。
単一回答にまつわる主な関数は sa1(), sa2(), sa3() の3つです。
(1) 単一回答の集計に関する関数 sa1
ここで扱い素材データは、ID(整理番号)、性別(男・女)、意見(賛成・中立・反対)の3列からなる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() が返す値はリストで、次の4つのテーブルを含みます。
・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() の第2引数("性別")と第3引数("意見")は、列名でなく列の番号で指定してもかまいません。次のように書いても同じ結果になります。
xx <- sa1(dtf, 2, 3)
tbl1だけを取り出したい時は、xx$tbl1 などのようにします。これには「計」がないので、カイ2乗検定等にかけることができます。例えば次のように書きます。
xx <- sa1(dtf, 2, 3) zz <- chisq.test(xx$tbl1)
こうすると、zzにカイ2乗検定の結果が代入されます。
tbl1, tbl2 などの各々のテーブルには、縦軸のタイトル「性別」と横軸のタイトル「意見」が付いています。これは、data01.csvの列名をそのまま割り当てたものです。
これらタイトル(正式には「次元名」というんでしょうか)を得るには、次のようにします。
dnn <- names(dimnames(xx$tbl1))
こうすると、dnnに「c("性別", "意見")」が代入されます。
sa1() の引数指定で、第3引数がなければ、1つの項目のみに着目した表が作成されます。
xx <- sa1(dtf, 2)
とした場合は、「性別」だけに着目した次のような表になります。tbl1, tbl2 の2つを示します。
−−−− $tbl1 男 女 41 50 $tbl2 男 女 総数 41 50 91 −−−−
上のような1つの項目だけに着目した表を作った場合、横方向には複数の列ができますが、縦方向は1行のみです。
したがって、縦軸に沿ったパーセント(pct2)は意味がないのでNULLになります。
(2) 単一回答の集計に関する関数 sa2, sa3
先にsa1() について述べましたが、sa2(), sa3() も同じように用います。
呼び出す時の引数は2つか3つで、その戻り値は、tbl1など4つのテーブルからなるリストです。
これら 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によって得られる表には「有効回答」という一種の小計欄があるので、カイ2乗検定等の統計検定の材料にはなりません。
(3) 単一回答集計のサンプルスクリプト
data01.csvを素材データとして、「性別」と「意見」のクロス集計を行うrubyスクリプトを掲げます。
Rプログラムの出力結果をそのまま標準出力に出力する他に、集計結果をrubyで受け取って、csvファイルとして書き出します。
−−−− tjs31.rb ここから #! ruby -Ks # 単一回答の集計(2項目のクロス集計) (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 ここまで
上は、2つの項目に着目したクロス集計(いわば2次元)の例ですが、1つの項目のみに着目した集計(1次元)にもそのまま使えます。
「xx$dnn <- names(dimnames(xx$tbl1))」というのは、縦・横の軸のタイトル名(次元名)を取得して、ruby側でそれを参照できるようにするためのものです。
次元名が2つあれば2次元(クロス集計)であり、1つだけなら1次元(1行のみの表)ということになります。
作られた表が2次元だとマトリクス(matrix)であり、1次元の場合は配列(array)です。それをruby側で受け取ると、どちらも配列(rubyのArray)の形になっています。しかし、こまったことに、単純に2次元配列|1次元配列として識別できるわけではありません。
Rの配列(1次元)は、ruby側では2次元配列にになります。「配列の要素番号・項目名・度数」が複数入っている2次元配列です。Rのマトリクスもrubyでは2次元配列なので、その区別がつきにくい。
そこで、dnnに次元名をセットしてrubyに引きわたしています。ruby側において、dnnが配列なら2次元であり、単なる文字列であれば1次元と判断できます。
少々くどくなりましたが、ともあれ tjs31.rb は、2次元の表にも1次元の表にも対応します。ただし、3次元以上の表は作れません。
(4) 補完的関数
tjsf.rには補完的な役割の関数も入っているので、いくつか紹介します。
以下では、素材データとして data02.csv, data02.mem を用います。
data02.csvの中身は data01.csv と同じですが、列名が「性別→sex」、「意見→view」と英語名になっています。また、性別欄は「男→m」、「女→f」。意見欄は「賛成:1, 中立:2, 反対:3」の対応関係で数字が書き込まれています。
data02.mem は、その対応関係をメモ書きしたもので、中身は次の2行です。
sex | m,f | 男,女 view | 1,2,3 | 賛成,中立,反対
半角1文字の '|' は区切り文字です。タブコードも区切り文字になります。区切り文字の前後にある半角スペースは無視されます。
区切り文字によって3つの部分に分けられていますが、各々の部分が何を意味するかは見当がつくと思います。最初が列名、あとの2つは、カンマ区切り形式の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として再設定します。
「remake.factor(dtf, "data02.mem")」は、data02.memの内容に即して、データフレームdtfの各列をfactorとして再設定するものです。
この関数を用いずに、あえて別の形で記述するなら次の2行になります。
dtf$sex <- factor(dtf$sex, levels=c("m","f"), labels=c("男","女")) dtf$view <- factor(dtf$view, levels=1:3, labels=c("賛成","中立","反対"))
・戻り値:再設定された結果(データフレーム)を返します。指定されたメモファイルが存在しないなどの理由で再設定できない場合は、再設定せずに(つまり素のままで)データフレームを返します。
・引数:第1はデータフレーム、第2は再設定用メモ書きファイルの名前です。
第3引数として、読み込むファイルの文字コードを指定することができます。data02.memがutf-8で書かれている場合は次のようにします。
dtf <- remake.factor(dtf, "data02.mem", "utf-8")
第3引数を省略すると、Rのデフォルトの文字コードになります。日本語版のWindowsだと "cp932" です。
◇ rename.dnn
・利用例:xx <- rename.dnn(xx, "性別", "意見")
・機能:テーブルの次元名を変構します。
サンプルスクリプトにおいて、xx(リスト)には tbl1, tbl2, pct1, pct2 の4つのテーブルが入っているわけですが、各々のテーブルの次元名は、「sex, view」になっています。これを「性別, 意見」に変更するのが rename.dnn() です。4つのテーブルすべてについて変更を行います。
・戻り値:第1引数xxの次元名を変更した結果を返します。変更できない時は、xxをそのまま返します。
・引数:第1引数は、テーブルからなるリストを指定します。あるいは、テーブル1つを指定することもできます。
第2引数と第3引数は次元名の指定です。第2が行名、第3が列名。2次元なら「"性別", "意見"」のように2つを指定します。1次元であれば "性別" のように1つだけを指定します(第3引数を省略します)。
◇ lround
・利用例:xx <- lround(xx, 1)
・機能:リストの各構成要素に対して、まるめ処理を行います。
リストの構成要素として、テーブル・マトリクス・ベクトルを想定しますが、それ以外が含まれていても支障ありません。
リストの構成要素が実数値(double)か複素数値(complex)でなければ、まるめ処理を施しません。
・戻り値:まるめ処理を施した結果を返します。
・引数:第1引数はリストを想定していますが、テーブル・マトリクス・ベクトルなどでもかまいません。
第2引数は、まるめ処理に係る小数点の桁数を指定するものです。デフォルトは0です。つまり整数値になるようにまるめます。これを1にすると、小数点1桁までの数値になります。
◇ get.dnn
・利用例:dnn <- get.dnn(xx)
・機能:引数xxの次元名を取得します。
xxがリストであれば、その中の最初にみつかったテーブルの次元名を返します。xxがテーブルなら、その次元名を返します。
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. 複数回答の集計
ここでは、「主な情報源」として、新聞・雑誌・テレビ・ラジオの4つの選択肢から、当てはまるものを複数個選んでもらうような複数回答の集計を取り上げます。
素材データとして data03.csv を用います。これは、ID(整数値), 年齢層(20代・30代・40代のどれか1つ), 主な情報源(新聞・雑誌・テレビ・ラジオの複数回答)から構成されています。98人分のデータです。具体的には次のようになっています。
ID,年齢層,新聞,雑誌,テレビ,ラジオ 1,30代,1,,, 2,40代,1,1,, 3,40代,,1,,
年齢層の列は単一回答です。
主な情報源の列は、4項目のうち選択された項目のところに数字の1を入力し、選択されなかったところは空欄にしてあります。
data03.csv は、乱数発生によって作成したデータです。実態とは関係ありません。
複数回答にまつわる主な関数は、ma1(), ma2() の2つです。他に 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 の2つの表からなるリストです。前者は度数(選択者数)、後者は比率(パーセント)の表です。
上のプログラムの実行結果は次のとおり。
−−−− $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("新聞","雑誌","テレビ","ラジオ"))
第1引数がデータフレーム、第2引数と第3引数が列の指定です。第2・第3は、順序を入れ換えても差し支えありません。
この ma1() は、次のサイトにある ma() を基にして作りました(というより、変更を加えました)。作者に感謝します。
(2) 複数回答の列をfactorとして再設定 ma.refactor
複数回答の列の値について少し触れます。
data03.csvの場合、選択されたら数値の1、そうでなければ空欄になっています。空欄は、read.csvで読み込んだ後に欠損値のNAとなります。
もし空欄でなく数字の0が入力されていたとすると、前述のRプログラムは、意図したような結果を出してくれません。
数字の1と0を機械的に列べると、当然ながら0の方が先にきます。そうすると、ma1() は、0の方を「選択あり」であるかのように処理してしまいます。
前もって、0よりも1の方が先にくるものであることをRに知らせてやる必要があります。そのためには、複数回答の列をfactorとして再設定します。
例えば、0と1からなる3列目を再設定するには次のようにします。
dtf[,3] <- factor(dtf[,3], levels=1:0, labels=1:0)
上の記述を3列目だけでなく3〜6列すべてに対して実行すればいいわけです。
これを少し簡単に記述できるようにするため、ma.refactor() という関数を設けました。次のように呼び出します。
dtf <- ma.refactor(dtf, 3:6, 1:0, 1:0)
こうすると、3〜6列のいずれも factor として再設定されます。
第1引数はデータフレーム、第2が列の指定(番号でなく名前でも可)、第3がlevels、第4はlabelsです。
第4引数を省略すると、「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() は、新聞・雑誌・テレビ・ラジオの4項目からなる複数回答の状況を一括して集計するので、全体像をみるのに向いています。ただ、その集計結果をそのまま統計検定にかけることはできません。
そこで、「新聞」にまるを付けたか否か、「雑誌」にまるを付けたか否かというように、各々の選択肢について二者択一でとらえ、集計表を4つ作ります。そうすると、各々の表をカイ2乗検定等にかけることができます。
前述したように、ma.refactor() の第4引数(labelsの指定)を省略すると、「新聞○」と「新聞×」のように設定されるので、あとは単一回答の集計を行う sa1() を適用してやれば、意図した表を作れます。
下に、それら4つの表を作成するスクリプトを掲げます。
単に表を作るだけでは何なので、カイ2乗検定にかけて、その結果を出力してみます。有意性が認められた時は、調整済み残差の表も出力します。
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 ここまで
参考まで、上のスクリプトを実行した結果の一部を掲げておきます。
「年齢層」と「新聞」のクロス集計表で、1つのセルが「13(41.9)」のように、度数と括弧書きのパーセントという形になっています。
クロス集計表の次は、カイ2乗検定の結果(カイ2乗値・自由度・有意確率・クラメール係数)を列記したものです。
−−−− 新聞 年齢層 新聞○ 新聞× 総数 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は乱数発生によって作ったものなので、カイ2乗検定で有意性が検出されることはないだろうと予想しましたが、上の「新聞」については有意性が認められました。有意性が検出されたのはこれだけです。
調整済み残差をみるまでもなく、「30代」の「新聞○」が多く、「新聞×」が少ない様子が分かります。
tjs34.rbで新たに出てきた関数について説明します。2つあります。
append.cramer(zz) は、引数zz(リスト)にクラメール係数を追加して返すものです。zzには、予めカイ2乗検定の結果を代入しておく必要があります。つまり chisq.test() の戻り値を代入しておきます。
「zz <- append.cramer(zz)」としたとき、zz$cramer にクラメール係数の値がセットされています。
もう1つの join.tbl() は、2つの表を結合するものです。典型的な例として、度数分布の表と構成比の表の結合があります。上のサンプルスクリプトもその一例です。
tt <- join.tbl(tbl, pct, "%d(%.1f)")
とすると、tbl, pct の2つの表が結合されて、その結果が tt に代入されます。ここで結合というのは、対応するセル1つづつを組み合わせるという意味です。
第3引数は、結合セルの書式を指定するものです。これが内部処理においてsprintfのフォーマットとして使われます。"%d(%.1f)" の場合、tblのセルが整数値として書き出され、pctのセルが括弧書きの小数点数(有効桁数1桁)で書き出されます。
もし pct が表として tbl よりも小さくて、行数または列数が少なければ、その少ない部分については tbl のセルだけが書き出されます。その場合、第3引数の書式は用いません。
rubyスクリプトの部分については説明を省略します。
拙作rrxライブラリのメソッドをいくつか使っているので分かりにくいと思いますが、必要に応じ関連のマニュアルを参照して下さい。
◇ append.cramer
・利用例:zz <- chisq.test(xx); zz <- append.cramer(zz)
・機能:引数zzにクラメール係数を追加して返します。zzには予めカイ2乗検定の結果を代入しておきます。
・戻り値:クラメール係数を含むカイ2乗検定結果を返します。戻り値zzのzz$cramerがクラメール係数の値です。
・引数:カイ2乗検定結果を引数として与えます。
◇ join.tbl
・利用例:tt <- join.tbl(tbl, pct, "%d(%.1f)")
・機能:2つの表を結合します。対応するセルを組み合わせます。
・戻り値:結合した表を返します。
・引数:第1引数と第2引数は表データ、第3引数は結合セルの書式で、内部処理においてsprintfのフォーマットとして使われます。省略時は "%s(%s)" とみなされます。
(4) 複数回答の選択肢相互の関係をみる ma2
「主な情報源」の4つの選択肢(新聞・雑誌・テレビ・ラジオ)の相互の関係をみるには、「新聞」と「雑誌」の両方にまるを付けた人がどれくらいいるか、「新聞」と「テレビ」の両方にまるを付けた人はどうか、といった確認をすることになると思います。そうすると、例えば、「新聞」を選ぶ人は「ラジオ」も選ぶ傾向が高いといったことが分かるかもしれません。
この「両方にまるを付けた人がどれくらいか」を一覧表示するのが 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() は、このような表を返しますが、4種類の表を生成して返します。
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 を呼び出す時の引数ですが、第1引数がデータフレーム、第2引数は複数回答の列の指定です。列の指定は、列番号でも列名でも、どちらでもかまいません。
ma2() の戻り値は、次の4つの表を含むリストです。
- 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%とした時に、各セルの値が何%かを確認するのに使います。全体の中で「両方にまるを付けた人」が何%いるかが分かります。
4つの表のうち 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() が生成する表は、統計検定にかけることができません。そこで、各選択肢をばらばらにして集計してみます。
「新聞」と「雑誌」の2つに注目した時の表は次のとおり。
−−−−
雑誌○ | 雑誌× | |
新聞○ | 29 | 23 |
新聞× | 27 | 19 |
−−−−
選択肢が4つなので、そこから2つを抜き出して表を作るやり方は6通りあります。
上のような2行・2列の表を6つ生成して、各々を統計検定にかけるためのサンプルを下に掲げます。カイ2乗検定を用います。
−−−− 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 ここまで
複数回答の答えは「○」か「×」かの2通りなので、各々の表は2行・2列になります。なので、カイ2乗検定の他に、フィッシャーの直接確率検定を適用することができます。
4つのセルの中に5未満のものがあるなど度数が少数のケースでは、カイ2乗検定を行った時に、警告メッセージが出ることがあります。そうした時はフィッシャーの検定を使うのがいいかもしれません。
zz <- fisher.test(xx$tbl1) print(zz$p.value)
このようにすると、p値が表示されます。これが 0.05 未満であれば、有意性が認められることになります。
有意性が認められれば、例えば、「新聞を主な情報源にする人は、ラジオも情報源にする傾向がある。」といったことが言えるかもしれません。しかし、有意性が認められないのであれば、少なくとも統計検定の結果からは、そうした傾向を指摘する根拠が得られない、ということになります。ちなみに、6つの表の中で有意性が認められるものはありませんでした。
このような分解的集計・検定の他に、数量化IIIや対応分析を用いることもできます。そのやり方について「単純集計に関する覚え書き/複数回答の集計」でごく簡単に触れましたが、関連の書籍やWebがいろいろあるので参考にして下さい。
(6) 複数回答の選択個数をみる ma.count
複数回答の選択肢について、回答者各人がまるを何個つけたかを検出しておくと便利です。
IDが1番の人は2個、2番の人は0個、3番の人は3個などのように、各人の選択個数を把握します。そうすれば、選択個数0の人(いわば無回答の人)を除いて集計したりすることが容易になります。
この選択個数の検出を行う関数 ma.count() を設けました。
これを呼び出す時の引数は4つありますが、第1引数がデータフレーム、第2が複数回答の列番号(または列名)です。多くの場合、第3と第4は省略可能です。
とりあえず 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 −−−−
選択個数0の人が4人、選択個数1個の人が24人いることが分かります。
ma.count() の戻り値はベクトルです。dtfには98人分のデータが入っているので、戻り値である count にも98人分の選択個数が入ります。
前述の例では「ma.count(dtf, 3:6)」のように、引数を2つしか指定しませんでしたが、第3引数と第4引数を省略せずに指定すると、次のようになります。
count <- ma.count(dtf, 3:6, 1, FALSE)
第3引数の1は、素材データであるdtfにおいて、複数回答の「選択あり」を示す入力値です。
dtfでは、選択された時に数字の1、選択されなかった時は空欄になっています。そこで、選択個数をカウントする際に数字の1日着目するというのが第3引数の意味です。
dtf <- ma.refactor(dtf, 3:6, c(1,NA), c("yes","no"))
上のように複数回答の列をfactorとして再設定していると、「選択あり」を示すのは "yes" になるので、ma.count() の呼び出し方は次のようになります。
count <- ma.count(dtf, 3:6, "yes", FALSE)
第3引数が "yes" のように文字である場合、実は、それを正規表現とみなします。内部処理では、Rのgrep()によってマッチングを行います。なので、厳密に「yes」の3文字を示すなら "^yes$" とすべきかもしれません。しかし、もう1つの値が "no" なので、そこまでしなくても意図どおりの結果が得られます。
第4引数は、選択個数を素材データのdtfに挿入・追加するか否かの指定です。デフォルトは FALSE です。つまり、挿入・追加をせず、選択個数を記録したベクトルを戻り値として返します。
TRUEにすると、dtfに nnn という名前の列を設けて、それに選択個数を記録します。戻り値として、nnnを追加・挿入したデータフレームを返します。
戻り値がdtf2に代入されている場合、dtf2$nnn または dtf2[,"nnn"] によって選択個数を参照できます。既に素材データのdtfに nnn という名前の列がある時は、そのnnnを上書きします。
選択個数0の無回答の人を除いたデータフレームを得たい時は、次のようにします。
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 の2つの表からなるリストを返します。前者は度数、後者はパーセント。
・引数:第1引数はデータフレーム、第2と第3は列の指定です。どちらかが複数回答の列、もう一方は単一回答の列です。その順番は問いません。
◇ ma2
・利用例:xx <- ma2(dtf, 3:6, 1, FALSE)
・機能:複数回答の選択肢相互の関係をみるための表を生成します。
・戻り値:tbl1, tbl2, pct1, pct2 の4つの表からなるリストを返します。4つの表が何かは、本文を参照して下さい。
・引数:第1引数はデータフレーム、第2引数は複数回答の列の指定です。
◇ ma.refactor
・利用例:dtf <- ma.refactor(dtf, 3:6, c(1,NA), 1:0)
・機能:複数回答の列をfactorとして再設定します。
・戻り値:再設定された結果(データフレーム)を返します。
・引数:第1引数はデータフレーム、第2引数は列の指定、第3はlevels、第4はlabelsです。
第4引数を省略すると、「c("新聞○", "新聞×")」が指定されたものとみなされます。「新聞」のところは各々の列名に置き換えられます。
◇ ma.count
・利用例:count <- ma.count(dtf, 3:6, 1, FALSE)
dtf2 <- ma.count(dtf, 3:6, "○", TRUE)
・機能:複数回答の選択個数を検出します。
・戻り値:第4引数がFALSEの時は、選択個数を記録したベクトルを返し、TRUEの時は、選択個数を記録した列(nnn)が追加・挿入されたデータフレームを返します。
・引数:第1はデータフレーム、第2は複数回答の列の指定、第3が「選択あり」を示す値(文字だと正規表現パターン)、第4は列の追加・挿入の有無です。
------------------------------------------------------------------------
3. 数値データの集計
ここでは、数値データとして「年齢」と「支出」(月当たりの支出額)を含むcsvデータを材料にします。
data04.csvは、ID(整数値)、性別(m,f)、年齢(数値)、支出(数値)から構成されています。98人分のデータです。具体的には次のようになっています。
ID,性別,年齢,支出 1,m,29,82249 2,m,49,93181 3,f,42,93647
数値データの集計にまつわる主な関数は、va1(), va2() の2つです。どちらも、単一回答と数値データの組合せの集計に用います。単一回答の「性別」(男女別)で分類しながら、数値データの平均値等を算出します。
他に date.split(), var2(), sd2() をtjsf.rに盛り込んであります。
date.split() は、「単純集計に関する覚え書き/数値データの集計」で触れたので今回は言及しません。仕様は同じです。
var2(), sd2() は、それぞれ標本分散と標本標準偏差を返す関数です。欠損値NAを取り除いた上で値を算出します。わたす引数は1つのみです。
ちなみに、Rに最初から組み込まれている var(), sd() は、それぞれ不偏分散と不偏標準偏差を返す関数です。こちらは引数が1つとはかぎりません。
(1) 1つの値のみを返す関数を用いた集計 va1
mean() は平均値を返す関数ですが、このように1つの値のみを返す関数によって集計を行う際、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() の第1引数は、素材となるデータフレームです。
第2引数は、分類の手がかりとなる単一回答の列を指定します。名前でも番号でもどちらでもかまいません。
第3引数では数値データの列を指定します。「c("年齢","支出")」とか「3:4」などのように複数指定してもかまいません。ちなみに、第2引数の方は1つだけです。
最後の第4引数は、適用したい関数の名前を指定します。省略すると mean が指定されたものとみなされます。
関数は、既存のものでなく自作のものでもかまいません。例えば、1つのセルを「平均値(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 # 1つの値を返す関数による集計 (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() と同じです。
第4引数で指定する関数が、1つの値だけでなく複数の値を返す場合に va2() を用います。
例えば、数値データに関して summary() は6つの値を返しますが、これを関数として指定したい時は次のようにします。
xx <- va2(dtf, 2, 3:4, summary)
実は、第4引数を省略すると、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.
Keyword(s):
References: