T. Yoshiizumi - tjs_c_guide Diff
- Added parts are displayed like this.
- Deleted parts are displayed
like this.
単純集計に関する覚え書き／数値データの集計
〜 rubyとRを用いた処理 〜
最終更新日: 2013/07/15
以下に掲げるドキュメントは、
[[tjs_c.zip|http://cup.sakura.ne.jp/tjs_c.zip]]
に同梱されている tjs_c.txt と同じ内容です。
＜はじめに＞
「単純集計に関する覚え書き」の第３段として、数値の集計について記します。年齢、身長や体重、価格のように数値として表現できるものであって、平均値とか最大値等を求めることができるものを対象とします。
ちなみに、第１段が[[「単一回答の集計」|http://cup.sakura.ne.jp/hiki/hiki.cgi?tjs_a_guide]]、第２段は[[「複数回答の集計」|http://cup.sakura.ne.jp/hiki/hiki.cgi?tjs_b_guide]]でした。
ソフトウェアとしては R と ruby を用います。拙作 rrxwin(rubyからRを利用するためのライブラリ) も使います。
サンプルのスクリプトを試した環境は次のとおり。
* Windows XP | VISTA
* ruby ver 1.8.7 | 1.9.3
* R ver 2.15.2
--------
{{toc_here}}
--------
!1. 数値データの分布やばらつきを確認
ここでは、「ID・性別・年齢・支出」の４項目からなるcsvデータ data06.csv を取り扱います。100人分のデータです。
IDは、個々人を識別するための整数値です。性別はアルファベットのm, fで記入されており、前者が「男」、後者が「女」です。
年齢と支出が数値データで、どちらも整数値です。支出は、月あたりの支出額のつもりです。乱数発生させた数値に少し手を加えたものです。
以下、このデータを素材にしてサンプルを示します。
!!(1) Rプログラムのsummary関数
○ 要点１：最小値・最大値・平均値等を一括取得するsummary()
○ 要点２：パーセンタイルの考え方による分布の把握
カテゴリー別の集計(男女別の集計)の前に、数値データそのものの分布をみることにします。Rプログラムではsummary関数が用意されており、最小値・最大値・平均値等を一括して取得できます。
csvデータを読み込んだ結果が dtf というデータフレームに収納されている場合、summary(dtf[,"年齢"]) とすれば、年齢に関する分布情報を得られます。
分布情報には次の６つがあります。年齢に関する具体的な値と合わせて記します。
||Min.||最小値||20.00
||1st Qu.||第１の４分位点||28.00
||Median||中央値(第２の４分位点)||34.00
||Mean||平均値||34.29
||3rd Qu.||第３の４分位点||41.00
||Max.||最大値||49.00
上の年齢に関する具体的な値をみると、最も若い人が20歳、最も年長の人が49歳であることが分かります。
平均年齢は 34.29歳です。
第１の４分位点が28なので、ざっくりした言い方ですが、28歳以下が25％いて、28歳以上が75％いることになります。サンプルの総数が100人なので、28歳以下が25人、28歳以上が75人いる(それくらいになる)という目安になります。
実際にデータを数えてみると、25歳「未満」だと24人、25歳「以下」だと27人です。28歳を基準線にして、25％と75％にきちんと厳密に分けられるわけではありません。この辺の詳細について述べるためには「パーセンタイル」の考え方とその計算方法に踏み込まなければなりませんが、ここでは素通りします。
パーセンタイルに関連する値として、中央値が34歳、第３の４分位点が41歳と出ています。34歳以下が50％(50人)くらい、41歳以下が75％くらいいることになります。
20〜49歳の集団としては、34.29歳という平均値は「なるほど」という典型的な値にかなり近いです。中央値(34歳)のところで50％づつに分かれるのも「まさに」という値です。そして、28歳以下が25％、41歳以上が25％となると、特定の年齢層に偏ることなく全体的にばらついていることが推測できます。乱数発生で作ったデータなので当然といえばそれまでですが。
ちなみに、Rプログラムでは、パーセンタイルに関する代表的な値を算出する関数として quantile() というのがあります。使い方は summary() と同じです。
ここで、summary関数を用いるサンプルを掲げておきます。data06.csvの「年齢」と「支出」の２つの項目をsummary() にかけています。
−−−− tjs20.rb ここから
#! ruby -Ks
# summary関数の利用 (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data06.csv", header=T)
items <- c("年齢", "支出")
for (i in items) {
cat(i, "の分布状況\n", sep="")
print(summary(dtf[,i]))
cat("\n")
}
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs20.rb ここまで
!!(2) 平均値と標準偏差
○ 要点１：不偏標準偏差と標本標準偏差の違い
○ 要点２：sapply関数の利用
集団の数値データの代表値として、「平均値」と「標準偏差」を並記する方法があります。平均値だけだと、ばらつきの度合いが分からないので標準偏差(sd:Standard Deviation)を合わせて記します。標準偏差が大きいほど、ばらつきが大きいということになります。
平均値をμ、標準偏差をσとしたとき、おおよそ次のような分布になるようです。
(μ−σ) 〜 (μ＋σ) の範囲に、全構成要素の約68％が入る。
(μ−σ×２) 〜 (μ＋σ×２) の範囲に、全構成要素の約95％が入る。
「標準偏差」を２乗したものを「分散」といいます。というより、「分散」の平方根が「標準偏差」です。
Rプログラムでは、分散を求める関数として var(), 標準偏差を求める関数として sd() があります。どちらも引数として、数値から構成されるベクトルを与えます。
ここで注意しなければならないのは、var(), sd() は、それぞれ不偏分散と不偏標準偏差を算出するものだということです。これらは、標本分散と標本標準偏差と少し違います。
細かな説明は省略しますが、100人のデータから分散を求める場合、計算の最終段階で、100で割り算する方法と、99(100から1を差し引いた値)で割り算する方法の２通りがあります。100で割るのが標本分散、99で割るのが不偏分散です。Rのvar関数は、不偏分散(99で割る方)です。標本分散より少し大きな値になります。
不偏分散(99で割ったもの)の平方根が不偏標準偏差、標本分散(100で割ったもの)の平方根が標本標準偏差です。
Rのvar()が不偏分散を返すのであれば、その値を基にして標本分散を求めるには、99を掛けた上で、改めて100で割ればいいことになります。つまり 99/100(0.99)を掛け算してやります。もちろん 99, 100 という数値は、実際の集団の構成要素数に置き換える必要があります。
不偏分散は、素直に100で割るのではなく、99で割ります。なんでこんなことをするのかという理由について、きちんと説明するのは私には難しいですが、おおよそ次のようなことのようです。
仮に、１万人の集団があったとして、その集団の年齢分布を知りたいとします。お金や時間が潤沢で、１万人全員を調査できたとすれば、集められたデータから分散を求める場合、素直に１万で割り算します。つまり標本分散を使います。標準偏差も標本標準偏差を用います。不偏分散や不偏標準偏差は使いません。
しかし、１万人全員の調査をやる余裕がなく、100人だけ無作為抽出して調査した場合は、１万人分の状況を推測する際、100で割り算するのでなく 99で割ります。100で割るよりも99で割った方が、１万人の母集団の分散に近い、というのが統計学の知見のようです。
要するに、部分から全体を推測する時は「不偏」の分散と標準偏差を用い、調査した集団そのものだけに着目する時は(つまり、より大きな集団について推定したりしない場合)「標本」の分散と標準偏差を用いるということのようです。
平均と標準偏差を並記する場合、適当な方の標準偏差を選んで記すということになるのだと思います。
なお、Wikipediaの「標準偏差」のところを読むと、「不偏」と「標本」の表現の使い分けが必ずしも統一されているわけではないようです。ご注意ください。
話を戻します。
Rプログラムで、標本分散や標本標準偏差を算出するには、例えば次のようにします。xが数値データからなるベクトルだとします。
n <- length(x) # xの要素個数
var100 <- var(x)*(n-1)/n # 標本分散
sd100 <- sqrt(var100) # 標本標準偏差
もしベクトルxの中にNAが含まれている可能性があるのであれば、それを取り除いた上で計算するのが無難です。最初の方に次の１行を挿入します。
x <- x[!is.na(x)] # NAを除去
下にサンプルを掲げてみます。「年齢」と「支出」について、平均と標本標準偏差を出力します。
−−−− tjs21.rb ここから
#! ruby -Ks
# 平均と標本標準偏差 (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data06.csv", header=T)
items <- c("年齢", "支出")
for (i in items) {
x <- dtf[,i]
x <- x[!is.na(x)]
n <- length(x)
sd100 <- sqrt(var(x)*(n-1)/n)
cat(sprintf("%sの平均値: %.1f, 標本標準偏差: %.3f\n", i, mean(x), sd100))
}
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs21.rb ここまで
ここで少し横道にそれますが、上の tjs21.rb の出力結果を表形式で出力することを考えてみます。
my.test() という関数を定義して、与えられたベクトルから平均と標本標準偏差を算出するようにします。
また、Rプログラムのapply関数群の１つである sapply() を用います。
出力結果を先に示すと次のとおり。
|| ||年齢||支出
||平均||34.290000||99673.78
||標本標準偏差||7.750219||28716.33
上の出力を得るためのスクリプトは下のとおりです。
−−−− tjs21b.rb ここから
#! ruby -Ks
# 平均と標本標準偏差(sapply利用) (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
my.test <- function(x) { # mean, var, sd の算出
x <- x[!is.na(x)]
n <- length(x)
m <- mean(x) # 平均値
var99 <- var(x) # 不偏分散
sd99 <- sd(x) # 不偏標準偏差
var100 <- var99 * (n-1) / n # 標本分散
sd100 <- sqrt(var100) # 標本標準偏差
return(c(平均=m, 標本標準偏差=sd100))
}
dtf <- read.csv("data06.csv", header=T)
items <- c("年齢", "支出")
xx <- sapply(items, function(i) my.test(dtf[,i]))
xx
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs21b.rb ここまて
上のスクリプトにおいて、my.test() の中で「不偏」と「標本」の２種類の分散と標準偏差を一通り算出していますが、戻り値として返しているのは「平均」と「標本標準偏差」の２つだけです。
sapply() という関数を使っていますが、これについてちゃんと説明する自信がないので、同じことを for を使って書くとどうなるか示しておきたいと思います。該当部分のみ掲げます。
−−−− sapplyの代替 ここから
items <- c("年齢", "支出")
xx <- NULL
for (i in items) xx <- cbind(xx, my.test(dtf[,i]))
colnames(xx) <- items
xx
−−−− sapplyの代替 ここまで
sapply() は、第１引数 items の各々の要素 i について、function(i) を適用するということだと思います。そして、その function の実態は、my.test(dtf[,i]) であるというわけです。
ちなみに、function(i) の定義は、sapply() の引数指定の中に長々と盛り込むことができます。my.test() の中身をsapply() の引数の一部として書き込むことができます。なので、my.test関数の定義は不要になります。例えば次のように書きます。
−−−− sapply関数の記述例 ここから
xx <- sapply(items, function(i) {
x <- dtf[,i]
x <- x[!is.na(x)]
n <- length(x)
m <- mean(x) # 平均値
sd100 <- sqrt(var(x) * (n-1) / n) # 標本標準偏差
return(c(平均=m, 標本標準偏差=sd100))})
print(xx)
−−−− sapply関数の記述例 ここまで
apply関数群は、他のところで使うことになるので、ちょっと横道にそれましたが取り上げてみました。
!!(3) 数値データのカテゴリー化
○ 要点１：cut() によって数値データにカテゴリー名を割り当てる
○ 要点２：seq() によって等差数列を生成
data06.csvの「年齢」は、20〜49歳の範囲にあります。これを「20代・30代・40代」の３分類でカテゴリー化することを考えます。そして、各々のカテゴリーに属する人の人数を出力してみます。
まず、出力結果は次のようになります。
年齢層別の人数
20代 30代 40代
32 36 32
この出力を得るためのスクリプトは、下のとおりです。
−−−− tjs22.rb ここから
#! ruby -Ks
# 数値データのカテゴリー化(cut関数の利用) (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data06.csv", header=T)
xx <- cut(dtf[,"年齢"], breaks=c(20, 30, 40, 50),
labels=c("20代", "30代", "40代"),
right=FALSE)
table(年齢層別の人数=xx)
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs22.rb ここまで
data06.csvの「年齢」の列をcut() に引きわたして、各々の人の年齢にカテゴリー名を割り当てます。年齢が20以上・30未満なら「20代」、30以上・40未満なら「30代」というカテゴリー名です。
cut() の戻り値 xx には100人分のデータが入っています。28 とか 36 とか 43 といった具体的な年齢ではなく、「20代・30代・40代」の３種類のどれかです。それが100人分収納されています。
といっても、xx は factor なので、カテゴリー名が文字列の形で記録されているわけではなく、整数値の 1, 2, 3 が100人分記録されています。数字の1には「20代」、2には「30代」、3には「40代」が割り振られています。
cutの引数 breaks は、年齢をカテゴリー区分する時の境目となる値を指定するものです。ベクトルで与えます。
サンプルでは c(20, 30, 40, 50) となっているので、20以上・30未満、30以上・40未満、40以上・50未満という区分になります。
このとき、１つの区間について下限の値は含むけれども上限の値は含まない、つまり「20以上・30未満」の20は含むけれども30は含まない、となるのは、常にそうであるわけではなく、cut() の引数 right を FALSE にしているためそうなるものです。
サンプルには出てきていませんが、include.lowest というオプションがあり、それと rightオプションの組合せによって、下限と上限を含めるかどうか指定できます。この辺の話は少々ややこしいので省略します。
breaksオプションとして与えるベクトルにおいて、全体の下限(20)が判然としない時は、マイナス無限大の -Inf を与えます。また、全体の上限(50)が判然としない時は、無限大の Inf を与えれば大丈夫だと思います。
labelsオプションは、カテゴリー名を列挙したベクトルを与えます。factor() の labels と同じような意味合いです。
labelsとして NULL を指定するか、あるいは省略すると、「[20,30)」 「[30,40)」 「[40,50)」というカテゴリー名が割り当てられます。
ところで、「20代・30代・40代」の３分類だと大雑把すぎると思われるかもしれません。そこで、５歳刻みで分類するケースを考えてみます。
breaksオプションとして c(20, 25, 30, …… 50) と記述するのは面倒なので、等差数列を生成する seq() を用います。
brk <- seq(20, 50, by=5)
とすれば、brkに５歳刻みのベクトルがセットされます。seq() の引数には、整数値だけでなく実数値も使えるので、いろいろ応用できると思います。
labelsオプションとして与えるベクトルの生成には、先述した sapply() を用いてみます。次のようにします。
lbl <- sapply(1:(length(brk)-1), function(i)
sprintf("%s〜%s", brk[i], brk[i+1]))
上の２点を盛り込んだスクリプト tjs23.rb を念のため下に掲げておきます。
−−−− tjs23.rb ここから
#! ruby -Ks
# 数値データのカテゴリー化(seq関数の利用) (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data06.csv", header=T)
brk <- seq(20, 50, by=5)
lbl <- sapply(1:(length(brk)-1), function(i)
sprintf("%s〜%s", brk[i], brk[i+1]))
xx <- cut(dtf[,"年齢"], breaks=brk, labels=lbl, right=FALSE)
table(年齢層別の人数=xx)
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs23.rb ここまで
seq() に関して少し補足します。
20〜50歳の範囲を５歳刻みなどではなく、４等分したいといったことがあると思います。４等分だと「20〜27.5、27.5〜35、35〜42.5、42.5〜50」となります。
この場合、seq() の byオプションではなく、length.outオプションを用います。
brk <- seq(20, 50, length.out=5)
とすれば、４等分にできます。length.out の値として、4ではなく5を指定する点に注意して下さい。nを指定すると「n-1」等分されます。
実際の年齢は整数値なので、27.5 とか 42.5 などの小数点数だと違和感があるし、「以上・未満」の境目も分かりにくくなるので、まるめ処理をして整数にしたいという場合は、次の関数を使えます。
floor(x) 少数点以下を切り捨て
ceiling(x) 少数点以下を切り上げ
round(x,0) 少数点以下を四捨五入
これら３つの関数とも、引数にベクトルを与えることができます。ベクトルの要素すべてについて、切り捨て・切り上げ・四捨五入を施します。
最後のroundの四捨五入は、単純な四捨五入ではなく、JISやISOで定められた方式に従った四捨五入のようなので注意して下さい。これに関しては次のサイトが参考になります。
[[JIS,ISO式四捨五入 - RjpWiki|http://www.okada.jp.org/RWiki/?JIS%2CISO%BC%B0%BB%CD%BC%CE%B8%DE%C6%FE]]
数値データの分布を確認する方法は、この辺にしておきます。
--------
!2. 数値データとカテゴリーデータの組合せによる集計
data06.csvには「年齢」と「支出」という数値データの他に、「性別」というカテゴリーデータがあります。
ここでは、男女別に年齢や支出の分布をみたり、カテゴリー化した「年齢」や「支出」と「性別」とのクロス集計を取り上げます。
既に述べた「1. 数値データの分布を確認」と、当連載の第１回目「単一回答の集計」のノウハウを組み合わせたものになります。
!!(1) １つの値(平均値)をカテゴリー別に算出
○ 要点：tapply() によって、分類しつつ関数を適用
data06.csvを読み込んだ結果が dtf に記録されているものとします。dtf はデータフレームです。
男女別に「年齢」と「支出」の平均値を得ようと思ったら、まず dtf の中から「男」に関するデータを male にセットし、「女」に関するデータを female にセットして、その上で各々の平均値を求めることになります。
これを比較的素直に実行しているのが下のプログラムです。実行結果とRプログラムの部分を示します。
|| ||年齢||支出
||男||34.38596||104923.07
||女||34.16279||92715.42
−−−− Rプログラムここから
dtf <- read.csv("data06.csv", header=T)
dtf$性別 <- factor(dtf$性別, levels=c("m","f"), labels=c("男","女"))
nums <- c("年齢", "支出") # 処理する数値欄の列名
male <- subset(dtf, 性別 == "男")
female <- subset(dtf, 性別 == "女")
result <- NULL
for (i in nums) {
m <- mean(male[,i]) # 「男」の平均値
f <- mean(female[,i]) # 「女」の平均値
result <- cbind(result, c(m,f))
}
colnames(result) <- nums
rownames(result) <- c("男", "女")
result
−−−− Rプログラムここまで
「dtf$性別 <- factor(……)」は、「性別」の欄に記入されているのが m, f なので、それを「男・女」に置換・調整するものです。
上のプログラムは、別のケースに応用しやすい形とはいえません。例えば、男女別ではなく地域別のデータがあった場合、それに対応させるためには、かなり書き換える必要が出てきます。地域別は、２分類とは限らず５分類とか10分類かもしれません。「男・女」という２分類前提の書き方をやめなければなりません。
levels(dtf[,"性別"]) が c("男","女") というベクトルを返すことを利用して、少し工夫してみたのが下のスクリプトです。分かりにくく、かつ、短くできたわけでもないので工夫のしがいがないですが、参考まで掲げてみます。
−−−− tjs24.rb ここから
#! ruby -Ks
# 男女別の年齢・支出の平均値 (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data06.csv", header=T)
dtf$性別 <- factor(dtf$性別, levels=c("m","f"), labels=c("男","女"))
key <- "性別" # 分類の手がかりとなる列の名前
nums <- c("年齢", "支出") # 処理対象の数値欄の列名
ctg <- levels(dtf[,key]) # == c("男", "女")
list.dat <- list() # 空のリスト
for (i in ctg) # dtfを男・女それぞれに分けてlist.datに追加
list.dat <- c(list.dat, list(subset(dtf, dtf[,key] == i)))
result <- NULL
for (i in nums) {
v <- NULL
for (j in 1:length(list.dat)) {
dd <- list.dat[[j]] # 男または女それぞれのデータフレームをddにセット
v <- append(v, mean(dd[,i])) # 年齢・支出それぞれの平均をvに追加
}
result <- cbind(result, v)
}
colnames(result) <- nums
rownames(result) <- ctg
result
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs24.rb ここまで
空のリスト list.dat には、１番目に「男」のデータフレーム、２番目に「女」のデータフレームを入れます。なので、list.dat［［1］］ が「男」のデータフレーム、list.dat［［2］］ が「女」のデータフレームです。
あとは、「年齢」と「支出」の平均を男・女それぞれについて求めているだけです。
さて、ここからが本論です。apply関数群を使って、上のプログラムを短くします。
まず、tapply() に着目します。これは、factorである列を含むデータフレームがあるとき、factorの要素の種類(性別欄の男・女)を手がかりにして分類しながら、指定の列(年齢・支出)を指定の関数(mean)に引きわたして、その結果を返します。
例えば次のように記述します。
xx <- tapply(dtf[,"年齢"], dtf[,"性別"], mean)
こうして得られたxxは、34.38596, 34.16279 という２つの要素からなる配列(array)です。１番目の要素が「男」の平均年齢、２番目が「女」の平均年齢です。これらは名前付きで記録されているので「nn <- names(xx)」とすれば、c("男","女") が得られます。
この tapply() を「年齢」と「支出」のそれぞれについて呼び出して、その２つの結果をマトリクスにすれば、前述のプログラムと同じ結果が得られるはずです。
tapply() を２回分記述するところに、for を用いてもいいのですが、既に言及した sapply() を用いると、更に簡潔になります。
ということで、だらだら説明してきましたが、その「だらだら」を踏まえて下のようなスクリプトを作りました。分かりやすいとはいえませんが、短くなったのは確かです。
−−−− tjs25.rb ここから
#! ruby -Ks
# 男女別の年齢・支出の平均値(tapply利用) (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data06.csv", header=T)
dtf$性別 <- factor(dtf$性別, levels=c("m","f"), labels=c("男","女"))
key <- "性別" # 分類の手がかりとなる列の名前
nums <- c("年齢", "支出") # 処理対象の数値欄の列名
result <- sapply(nums, function(i) tapply(dtf[,i], dtf[,key], mean))
result
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs25.rb ここまで
上のスクリプトが、前述した tjs24.rb よりも作業効率がいいのかどうかは分かりません。
「年齢」の平均値を出す時と、「支出」の平均値を出す時の各々で、男女別に分類する処理を実行しているような気がします。つまり、男女別の分類処理を２回実行しているのではないかと思います。
一方、長いスクリプトの tjs24.rb の方は、男女別の分類処理は１回のみです。分類した結果をlist.datにセットして、それを用いて「年齢」と「支出」の平均値を算出しています。
Rプログラムの内部処理がどうなっているか分からないので、上に書いたことが正しいかどうか自信はありません。また、もし正しいとしても、それが作業効率に影響するのかどうか、それも分かりません。データの規模が大きくなるほど作業効率が問題になるので、気になるところではありますが、ちゃんと確認していません。
なお、男女別の他に、男女を合わせた「全体」も付けた形の表を得たい時は、例えば次のようにします。もっと適切なやり方があるような気はしますが、私が思いついたものを掲げておきます。
−−−− tjs25b.rb ここから
#! ruby -Ks
# 男女別の年齢・支出の平均値(全体欄を付加) (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data06.csv", header=T)
dtf$性別 <- factor(dtf$性別, levels=c("m","f"), labels=c("男","女"))
key <- "性別" # 分類の手がかりとなる列の名前
nums <- c("年齢", "支出") # 処理対象の数値欄の列名
xx <- sapply(nums, function(i) tapply(dtf[,i], dtf[,key], mean))
ww <- sapply(nums, function(i) mean(dtf[,i]))
result <- rbind(xx, 全体=ww)
result
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs25b.rb ここまで
上のスクリプトでは、「全体」を別枠で算出して、後から男女別のものに付け加えています。
これを別枠ではなく「男・女・全体」を一括して算出する方法があるような気がしますが、分かりませんでした。
それから、平均値を示す場合、その算出に用いたサンプルの数が気になります。「男」の人数、「女」の人数、「全体」の人数がどれくらいかです。
男女それぞれの人数は、「nn <- table(dtf[,key])」で求めることができます。単一回答や複数回答の集計のところで多用した方法です。nnは、２つの要素からなる配列で、１番目が「男」の人数、２番目が「女」の人数です。
また、男女別の人数を「nn <- tapply(dtf[,key], dtf[,key], length)」として求めることもできます。
「全体」の人数は、length(dtf[,key]) で取得できます。
ここには掲げませんが、最も右の列に「人数」を加えて表示するスクリプト tjs25c.rb を圧縮ファイルに同梱してあるので、よかったら参考にして下さい。
!!(2) 複数の値(最小値・平均値・最大値等)をカテゴリー別に算出
○ 要点：lapply() の利用
平均値という１つの値だけを算出するのであれば、前述した tjs24.rb, tjs25.rb のやり方でいいわけですが、最小値・平均値・最大値といった複数の値の算出を考えると、もう少し工夫してやる必要があります。
平均値を求める mean() ではなく、分布情報を得るための summary() を用いる場合、「年齢」と「支出」の各々について、分布情報が男女間でどのように違うかをみるためには、次のような出力結果が必要になります。
||年齢
|| ||男||女
||Min.||20.00||21.00
||1st Qu.||28.00||27.50
||Median||35.00||33.00
||Mean||34.39||34.16
||3rd Qu.||41.00||41.00
||Max.||49.00||49.00
||支出
|| ||男||女
||Min.||50430||52490
||1st Qu.||82250||70020
||Median||113600||92460
||Mean||104900||92720
||3rd Qu.||129400||116900
||Max.||149400||148400
平均値だけを求める tjs25.rb では、「result <- sapply(……)」として、最終的な結果を sapply() で取得しました。
この sapply は、戻り値としてベクトルやマトリクスを返します。中身の計算結果が複雑化したとしても、それをなるべく単純なベクトルやマトリクスに変換して返します。
しかし、mean ではなく summary だと、複数の値を算出するので、sapply では適切な結果を得られません。このような時は lapply() を用います。lapply() は、Rプログラムでいうところのリストを返します。なので、複雑な階層構造のデータでも意図どおりの形で取得できます。使い方は、基本的に sapply() と同じです。
tjs25.rb の sapply を lapply に変更したプログラムを下に掲げてみます。Rプログラム部分です。
−−−− Rプログラムここから
dtf <- read.csv("data06.csv", header=T)
dtf$性別 <- factor(dtf$性別, levels=c("m","f"), labels=c("男","女"))
key <- "性別" # 分類の手がかりとなる列の名前
nums <- c("年齢", "支出") # 処理対象の数値欄の列名
result <- lapply(nums, function(i) tapply(dtf[,i], dtf[,key], summary))
names(result) <- nums
result
−−−− Rプログラムここまで
上のプログラムを実行した結果は、次のとおりです。「年齢」に関する部分のみ示します。
||$年齢
||$年齢$男
|| ||Min.||1st Qu.||Median||Mean||3rd Qu.||Max.
|| ||20.00||28.00||35.00||34.39||41.00||49.00
||$年齢$女
|| ||Min.||1st Qu.||Median||Mean||3rd Qu.||Max.
|| ||21.00||27.50||33.00||34.16||41.00||49.00
上の出力結果だと男女間の比較がやりにくいので、もう少し工夫してみます。「年齢」と「支出」の各々について、「男・女」と「最小値・平均値・最大値等」とのクロス表を作るようにします。
出力結果は次のとおり。
||$年齢
|| ||Min.||1st Qu.||Median||Mean||3rd Qu.||Max.
||男||20||28.0||35||34.39||41||49
||女||21||27.5||33||34.16||41||49
||$支出
|| ||Min.||1st Qu.||Median||Mean||3rd Qu.||Max.
||男||50430||82250||113600||104900||129400||149400
||女||52490||70020||92460||92720||116900||148400
上のように表示すれば、男女間の違いを比較しやすくなります。
この表示結果を得るためのスクリプトは下のとおりです。
−−−− tjs26.rb ここから
#! ruby -Ks
# 年齢・支出の男女別分布情報(lapply利用) (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data06.csv", header=T)
dtf$性別 <- factor(dtf$性別, levels=c("m","f"), labels=c("男","女"))
key <- "性別" # 分類の手がかりとなる列の名前
nums <- c("年齢", "支出") # 処理対象の数値欄の列名
xx <- lapply(nums, function(i) tapply(dtf[,i], dtf[,key], summary))
names(xx) <- nums
ctg <- levels(dtf[,key])
result <- lapply(nums, function(i) t(sapply(ctg, function(j) xx[[i]][[j]])))
names(result) <- nums
result
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs26.rb ここまで
上のスクリプトで、lapply() が２回出てきますが、１回目は tjs25.rb の sapply を lapply に変更しただけです。
２回目の lapply() の引数部分が分かりにくいかもしれません。この部分では、変数iに「年齢」と「支出」が逐次セットされ、変数jには「男」と「女」が逐次セットされます。
まず「年齢」に関して、summary() の出力結果を取りまとめます。「男」と「女」がばらばらになっていたものを１つの表にします。そして、「支出」に関しても同じように処理します。
sapply(……) を t() で囲んでますが、この t() は、表の縦軸と横軸を入れ換えるものです。今回のケースでは、これがないと縦に長い表になるので入れてみました。
参考まで、男女別の他に「全体」の欄を付加するスクリプトを下に掲げておきます。
−−−− tjs26b.rb ここから
#! ruby -Ks
# 年齢・支出の男女別分布情報(全体欄を付加) (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data06.csv", header=T)
dtf$性別 <- factor(dtf$性別, levels=c("m","f"), labels=c("男","女"))
key <- "性別" # 分類の手がかりとなる列の名前
nums <- c("年齢", "支出") # 処理対象の数値欄の列名
xx <- lapply(nums, function(i) tapply(dtf[,i], dtf[,key], summary))
names(xx) <- nums
ctg <- levels(dtf[,key])
yy <- lapply(nums, function(i) sapply(ctg, function(j) xx[[i]][[j]]))
names(yy) <- nums
ww <- lapply(nums, function(i) summary(dtf[,i]))
names(ww) <- nums
result <- lapply(nums, function(i) t(cbind(yy[[i]], 全体=ww[[i]])))
names(result) <- nums
result
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs26b.rb ここまで
ここには掲げませんが、上のスクリプトと同じ結果をタブ区切りテキストの形で出力する tjs26c.rb を圧縮ファイルに同梱しておきます。
男女別の結果(変数xx)と全体の結果(変数ww)をrubyの側に引きわたして、ruby側でその結合と出力を行います。
すべてRプログラムで記述する場合よるも長くなりますが、引きわたされたデータがどのような形になっているのかを確認する材料になります。よかったら覗いてみて下さい。
なお、tjs26c.rb を動かす際、rrxwin ver 1.02 が必要です。ver 1.01 までだと正常に動かないので注意して下さい。
!!(3) 数値データをカテゴリー化した上でクロス集計
○ 要点：数値データをcut()でカテゴリー化した上でクロス集計
数値データのカテゴリー化は、前述したように cut() で行うことができます。カテゴリー化すれば、他のカテゴリーデータと組み合わせてクロス集計できます。
例えば、「性別」と「年齢」のクロス集計を出力するRプログラムは次のとおりです。「年齢」を「20代・30代・40代」の３分類にします。
−−−− Rプログラムここから
dtf <- read.csv("data06.csv", header=T)
dtf$性別 <- factor(dtf$性別, levels=c("m","f"), labels=c("男","女"))
brk <- c(20, 30, 40, 50)
lbl <- c("20代", "30代", "40代")
age <- cut(dtf[,"年齢"], breaks=brk, labels=lbl, right=FALSE)
tt <- table(dtf[,"性別"], age)
names(dimnames(tt)) <- c("性別", "年齢層")
tt
−−−− Rプログラムここまで
上のプログラムに少し拡張性を持たせてみます。
「性別」 「年齢」 「支出」の３つをそれぞれ組み合わせてクロス集計するスクリプトを下に掲げます。
「年齢」は「20代・30代・40代」の３分類、「支出」は５万〜15万を２万刻みで分類します。
−−−− tjs27.rb ここから
#! ruby -Ks
# カテゴリー化した数値データを含むクロス集計 (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data06.csv", header=T)
dtf$性別 <- factor(dtf$性別, levels=c("m","f"), labels=c("男","女"))
item <- c("性別", "年齢", "支出")
dtf2 <- data.frame(dtf[,item[1]]) # 「性別」の列をdtf2に入れる
brk <- c(20, 30, 40, 50)
lbl <- c("20代", "30代", "40代")
xx <- cut(dtf[,item[2]], breaks=brk, labels=lbl, right=FALSE)
dtf2 <- cbind(dtf2, xx) # カテゴリー化した「年齢」をdtf2に追加
brk <- seq(50000, 150000, by=20000)
lbl <- sapply(1:(length(brk)-1), function(i)
sprintf("%s〜%s万", brk[i]/10000, brk[i+1]/10000))
xx <- cut(dtf[,item[3]], breaks=brk, labels=lbl, right=FALSE)
dtf2 <- cbind(dtf2, xx) # カテゴリー化した「支出」をdtf2に追加
colnames(dtf2) <- item
for (i in 1:(length(item)-1))
for (j in (i+1):length(item)) {
tt <- table(dtf2[,i], dtf2[,j])
names(dimnames(tt)) <- c(item[i], item[j])
print(tt)
cat("\n")
}
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs27.rb ここまで
上のスクリプトを実行すると、３つの表が出力されます。いずれも単一回答を組み合わせたクロス集計表とみなすことができます。合計欄を付加したり、パーセントの表にしたりする場合は、「単一回答の集計」で取り上げたノウハウを利用できます。
--------
!3. 日付データの取扱い
これまで取り上げてきた数値データは、「年齢」と「支出」の２つでした。どちらも整数値ですが、取り上げた処理ノウハウは実数値にも適用できます。
ただ、数値データの１種である日付データ(あるいは日付・時刻データ)を扱う場合は、それなりの工夫が必要になります。その詳細には立ち入りませんが、関連のサンプルをいくつか示してみたいと思います。
扱うcsvデータ data07.csv は、「ID・長男の誕生日・次男の誕生日」の３列から構成されています。100組のデータが入っています。次のような内容です。
ID,長男の誕生日,次男の誕生日
1,1988-10-30,1997-04-13
2,1992-02-04,1997-07-22
3,1995-11-22,1998-05-30
!!(1) 日付の分布状況を確認
○ 要点１：日付データは Dataオブジェクトに変換してから処理
○ 要点２：Dateオブジェクトを文字列に変換するのは format()
Rプログラムにおいて、data07.csv を read.csv() で読み込むと、「1992-02-04」などの日付が文字列として認識され、データフレームになった段階では factor になります。
dtf <- read.csv("data07.csv", header=T)
として読み込むと、dtf$長男の誕生日, dtf$次男の誕生日 のどちらもfactorです。
このままだと日付として扱うのは難しく、例えば、その日付の曜日を簡単に求めることができません。長男と次男の誕生日が何日隔たっているかといった計算もできません。
そこで、このfactorをDate型に変換します。次のようにすると変換できます。
dtf$長男の誕生日 <- as.Date(dtf$長男の誕生日)
dtf$次男の誕生日 <- as.Date(dtf$次男の誕生日)
こうすると、Date型になり、その実態は、1970年１月１日からの経過日数になるようです。つまり数値になります。もっと前の日付だとマイナスの値になります。
1970年１月１日の値が0、翌日の値が1、翌々日の値が2です。
といっても、xxがDate型であるとき、print(xx) で表示されるのは「1970-01-01」などの文字形式です。数値として表示させたい時は、print(as.numeric(xx)) とか print(as.vector(xx)) とします。
ちなみに、数値をDate型に変換する時は、「as.Date(c(0,7,14), origin="1970-01-01")」のように、originオプションによって起点となる日付を指定する必要があります。
日付を文字で表現する場合、「2013-07-01」のように、区切り文字としてマイナス記号を用いるのが基本のようですが、スラッシュを用いてもトラブルなくDate型に変換できるようです。そして、Date型に変換したものをprint() で表示すると、スラッシゅではなくマイナスで表示されます。
また、「2013-07-01」のように、月や日を２桁にするため0を挿入するのが一般的かもしれませんが、「2013-7-1」のように0を省略しても、Date型に変換すれば同じ数値になります。そして、print() で表示すると、0が挿入されて「2013-07-01」のようになります。
なお、「年」を省略して "05-15" などのように「月・日」だけの文字列は、Date型に変換しようとしてもエラーになります。今年の西暦年を自動的に補ってくれる、ということはないようです。
Date型は、基本的に数値なのでいろいろな関数で扱うことができます。例えば、分布状況を表示する summary() もその１つです。最小値(最も早い誕生日)、最大値(最も遅い誕生日)などを確認できます。
下にサンプルスクリプトを掲げておきます。
−−−− tjs28.rb ここから
#! ruby -Ks
# Date型への変換とsummaryの適用 (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data07.csv", header=T)
dtf$長男の誕生日 <- as.Date(dtf$長男の誕生日)
dtf$次男の誕生日 <- as.Date(dtf$次男の誕生日)
summary(dtf$長男の誕生日)
summary(dtf$次男の誕生日)
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs28.rb ここまで
上の実行結果は次のとおりです。長男のもののみ示します。横に長くなるので縦・横逆にして掲げます。
||Min.||"1988-08-25"
||1st Qu.||"1990-06-07"
||Median||"1992-08-09"
||Mean||"1992-07-21"
||3rd Qu.||"1994-09-04"
||Max.||"1996-12-18"
ここで、このDate型のsummary() の結果をrubyの側で受け取ることを考えます。
実は、単純にsummaryの結果を引きわたすと、「1996-12-18」のような文字列ではなく、9848のような数値になります。この数値は、1970年１月１日からの経過日数なので、rubyの側で日付に変換することができます。しかし、Rプログラムで算出されたものは、Rで文字列に変換した方が無難です。
Date型の数値を文字列に変換するのは、format() 関数で行います。例えば次のように呼び出します。
format(xx, "%Y-%m-%d")
こうすると、Date型の数値が「1996-12-18」といった文字列に変換されます。
以下に、rubyにデータを引きわたすサンプルを掲げます。
−−−− tjs28b.rb ここから
#! ruby -Ks
# Date型データをruby側に渡す (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data07.csv", header=T)
item <- c("長男の誕生日", "次男の誕生日")
for (i in item) {
dtf[[i]] <- as.Date(dtf[[i]])
xx <- format(summary(dtf[[i]]), "%Y-%m-%d")
xx$memo <- i
robj(xx)
}
EOS
##
hs = Rrx.rexec(rpro)
n = 0
while (xx = hs["robj#{n+=1}"]) != nil
memo = xx.delete("memo")
ary = []
xx.each do |key, val|
ary << [key, val]
end
puts memo
print Rrx.ary2str(ary, "\t")
printf("\n")
end
−−−− tjs28b.rb ここまで
余談ですが、Rプログラムのsummary() が返す値の構造は、ちょっと特殊です。単純な matrix とか table ではありません。
そして、summary() が返す値に format() を適用すると、構造が変化します。少し単純になります。rubyで受け取ると、素直なハッシュになっています。
summary() が返す値そのままだと、rubyで受け取ったとき、もっと複雑な構造なので注意して下さい。
!!(2) ２つの日付の隔たりを調べる
○ 要点：二つのDateオブジェクトの隔たりはdifftime()で取得
長男の誕生日と次男の誕生日がどれくらい隔たっているかを調べる場合、difftime() を用います。
dd <- difftime(dtf$次男の誕生日, dtf$長男の誕生日, units="days")
とすれば、各々の隔たりの期間が変数ddにセットされます。ddには100組のデータが収納されます。
unitsオプションで "days" を指定しているので、とりあえず隔たりが日数として記録されますが、ddの中身は、単なる日数(整数値)ではなく、difftimeという独自のクラスです。
unitsオプションには、"secs", "mins", "hours", "days", "weeks" が指定できます。サンプルでは日付・時刻のdate-time型ではなくDate型を扱っているので、"days", "weeks" のどちらかということになるでしょうか。
下に difftime() の利用例を掲げます。最も隔たりの大きい兄弟と、最も隔たりの小さい兄弟を表示します。
−−−− tjs29.rb ここから
#! ruby -Ks
# 誕生日の隔たりを算出 (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data07.csv", header=T)
dtf$長男の誕生日 <- as.Date(dtf$長男の誕生日)
dtf$次男の誕生日 <- as.Date(dtf$次男の誕生日)
dd <- difftime(dtf$次男の誕生日, dtf$長男の誕生日, units="days")
dtf <- cbind(dtf, 隔たり=dd)
dmax <- max(dd)
dmin <- min(dd)
xx <- subset(dtf, 隔たり == dmax)
cat("最も隔たりの大きい兄弟\n")
xx
cat("\n")
xx <- subset(dtf, 隔たり == dmin)
cat("最も隔たりの小さい兄弟\n")
xx
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs29.rb ここまで
上を実行すると、隔たりが日数で表示されます。これを「週」の数で示したい時は、difftimeオブジェクト「xx$隔たり」のunitsを次のようにして変更します。
units(xx$隔たり) <- "weeks"
こうすると、「460.7143 weeks」とか「59.57143 weeks」のように表示されます。
残念ながら「年」単位と「月」単位の指定はできないので、少し手間をかけて算出することになります。formatの書式指定(%Y, %m)によって「年」と「月」を文字として得られるので、それを数値に変換します。そうすれば、何年何ヶ月の隔たりかを計算できます。
xxに「最も隔たりの大きい兄弟」のデータが入っている場合、次のようにします。
y1 <- as.numeric(format(xx$長男の誕生日, "%Y"))
y2 <- as.numeric(format(xx$次男の誕生日, "%Y"))
m1 <- as.numeric(format(xx$長男の誕生日, "%m"))
m2 <- as.numeric(format(xx$次男の誕生日, "%m"))
yy <- y2 - y1
mm <- m2 - m1
if (mm < 0) {
mm <- 12 + mm
yy <- yy - 1
}
上のようにして求めた yy, mm は、「yy年mmヶ月の隔たり」ということになります。
ここには掲げませんが、長男と次男の誕生日の隔たりを「月数差」と「年月差」で求めて、最も隔たりが大きい組と、最も隔たりの小さい組を表示するスクリプト tjs29b.rb を圧縮ファイルに同梱したので、よかったら参考にして下さい。
!!(3) 日付データのカテゴリー化(年・月・日・四半期・曜日)
○ 要点１：cut() によって日付データを四半期単位で分類
○ 要点２：unlist(strsplit("2013-01-01", "-"))によって文字列分割
数値データをcut() を用いてカテゴリー化する方法については既に触れました。
日付データの「月」を数値で取得すれば、四半期単位(1〜3月、4〜6月、7〜9月、10〜12月の４分類)にするのは容易です。
xxがベクトルで、100個の日付データが入っている場合でいうと、次のようにして四半期単位の分類を行えます。
mm <- as.numeric(format(xx, "%m")) # 「月」の数を抽出
qq <- as.numeric(cut(mm, breaks=c(1,4,7,10,13), labels=1:4, right=F))
上のようにして求めた qq には100個の数値が入っています。その数値は 1〜4 のどれかです。第１四半期が1、第２四半期が2などです。
cut() が返す値は factor なので、それを as.numeric() で数値に変換します。
xxとqqを対応づければ、該当の日付が第何四半期に当たるかを確認できます。
日付データから「月」を取得する関数 months() もありますが、これが返す値は、"1月" とか "4月" などの文字列です。
xxが複数の日付データからなるベクトルのとき、「mm <- months(xx)」とすれば、mmは、「月」を表す文字列が複数入ったベクトルになります。factorではありません。
この mm を材料にして四半期単位の分類を行うこともできますが、「月」を数値にした方が分類しやすいので、先述の方法を採ることにします。
なお、weekdays() を使うと、「月」ではなく「曜日」に変換することができます。これも months() と同じように、文字列からなるベクトルを返します。
ここで、四半期単位での分類だけでなく、「年・月・日・四半期・曜日」に分けるための関数 date.split() を考えてみます。
ss が "2013-01-01" のような文字列である場合、「strsplit(ss, "-")」によって、マイナス記号を区切り文字としながら「年・月・日」の３つに分割できます。
ただし、そのstrsplit() が返す値はリストなので、unlist() によってベクトルに変換する必要があります。ベクトルに変換すれば、as.numeric() で数値化できます。
こうしたノウハウを利用しながら date.split() を記述します。
＊ date.split() の仕様
date.split() が返す値はデータフレームです。素材として与えられたデータフレームに、year, month, mday, quarter, wday という名前の列を付け加えて返します。５つの列を付加するわけです。
これらの追加された列は、それぞれ数値からなるvectorか、または、factorです。date.split() に与えるオプションによって、vector, factor のどちらにするかを選択できます。factorなら「2013年」とか「8月」などの文字列として扱えるし、一方、vectorなら 2013 とか 8 の数値として扱えます。
dd <- date.split(dtf, "長男の誕生日", "c")
とすると、dtf に year, month などの列を追加した dd を得られます。
第１引数のdtfは、素材となるデータフレームです。
第２引数の "長男の誕生日" は、日付データの列名です。列の番号で指定してもかまいません。指定した列が日付データでない時は、警告メッセージが出て、戻り値の dd は NULL になります。
第３引数が vector, factor の選択で、"n" か "numeric" だと vector、"c" か "character" だと factor になります。省略すると "c" が指定されたものとみなされます。つまり factor になります。
数値のvectorにおける「曜日」の値は 1〜7 で、1が日曜、2が月曜などです。
もし、素材となるデータフレーム dtf に year, month などの５つの列が既にあれば、戻り値である dd には year3, month3, mday3 などのように番号が付いた列名が追加されます。
この番号は、指定された日付データの列番号です。上に掲げた例でいうと、"長男の誕生日" が２番目の列に当たるので、year2, month2, mday2 などになります。
素材である dtf に year, month, mday などの５つの列がなければ、この番号は付きません。
以下、少し長くなりますが date.split() の定義を含むサンプルを掲げます。
「長男の誕生日」を四半期単位で数え上げます。第1四半期が何人、第2四半期が何人というふうに集計した結果を出力します。
−−−− tjs30.rb ここから
#! ruby -Ks
# 日付データのカテゴリー化 (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
date.split <- function(dtf, cc, xtype="c") {
cn <- colnames(dtf)
if (is.character(cc)) cc <- seq_along(cn)[cn %in% cc]
date.dat <- dtf[,cc]
if (class(date.dat) != "Date") {
warning("Date型を指定して下さい.")
return(NULL)
}
if (xtype == "character") xtype <- "c"
else if (xtype == "numeric") xtype <- "n"
if (xtype != "c" && xtype != "n") xtype <- "c"
cname <- c("year", "month", "mday", "quarter", "wday")
nn <- seq_along(cn)[cn %in% cname]
if (length(nn) > 0) # dtf中に既に year, month などの列がある
cname <- sapply(cname, function(i) sprintf("%s%d", i, cc))
fmt <- c("%s年", "%s月", "%s日", "第%s四半期")
wname <- c("日曜日","月曜日","火曜日","水曜日","木曜日","金曜日","土曜日")
ymd <- t(sapply(format(date.dat, "%Y-%m-%d"), function(ss)
as.numeric(unlist(strsplit(ss, "-")))))
qq <- as.numeric(cut(ymd[,2], breaks=c(1,4,7,10,13), labels=1:4, right=F))
ww <- weekdays(date.dat)
for (i in 1:length(wname)) ww[ww==wname[i]] <- i
dd <- data.frame(cbind(ymd, qq, ww))
colnames(dd) <- cname
if (xtype == "c") {
for (i in 1:(ncol(dd)-1)) { # 年・月・日・四半期のfactor化
dd[,i] <- as.factor(dd[,i])
lvl <- levels(dd[,i])
lbl <- sapply(lvl, function(j) sprintf(fmt[i],j))
dd[,i] <- factor(dd[,i], levels=lvl, labels=lbl)
}
dd[,ncol(dd)] <- factor(ww, levels=1:7, labels=wname) # 曜日のfactor化
}
for (i in cname) dtf[,i] <- NULL # 重複する列名を削除
return(cbind(dtf, dd))
}
dtf <- read.csv("data07.csv", header=T)
dtf$長男の誕生日 <- as.Date(dtf$長男の誕生日)
dtf$次男の誕生日 <- as.Date(dtf$次男の誕生日)
dd <- date.split(dtf, "長男の誕生日", "c")
table(dd$quarter)
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs30.rb ここまで
date.split() の定義の最初のところで、日付データであるはずの date.dat が本当にDate型か否かをチェックしています。class(date.dat) の戻り値が "Date" であればDate型です。そうでない時は警告メッセージを出して NULL を返します。
sapply(format(……)) のところは分かりにくいですが、結果として、マトリクス(ymd)を返します。ymdの１列目には「年」、２列目には「月」、３列目には「日」がセットされます。いずれも数値データです。date.datが100個から構成されている場合、ymdの行数が100になります。
date.split() の定義に関するところは、数値のvectorをfactorに変換する処理などがあるため少しごちゃごちゃしていますが、説明は省略します。
念のため、スクリプトの出力結果を掲げておきます。長男の誕生日を四半期単位の分布でみたものです。
||第1四半期||第2四半期||第3四半期||第4四半期
||22||15||32||31
！！［補足］ 日付データへの cut(), seq() の適用について
前述した四半期単位でのカテゴリー化は、日付データの「月」を数値として取り出した上で cut() にかけました。つまり、日付データそのものをcut() にかけたわけではありません。
cut() は、Date型の日付データを直接取り扱うことができます。その場合、breaksオプションにはDate型からなるベクトルを指定します。
例えば、日付データ100個からなるxxに、2010〜2012年の３年間のデータが入っている場合、それを年単位でカテゴリー化するには次のようにします。
brk <- as.Date(c("2010-01-01", "2011-01-01", "2012-01-01", "2013-01-01"))
lbl <- c("2010年", "2011年", "2012年")
yy <- cut(xx, breaks=brk, labels=lbl, right=F)
上のようにすると、「2010年」〜「2012年」の３つの値から構成されるyyがfactorとして生成されます。yyの要素は100個あります。
yyはカテゴリーデータなので、それを基にしてクロス集計等を行うことができます。
それから、2013年の365日分の日付データを生成するのに seq() を使うことができます。
xx <- seq(as.Date("2013-01-01"), as.Date("2013-12-31"), by="days")
こうすると、xxに365日分の日付データがセットされます。
"days" のところには "months", "weeks" を書くこともできます。
"months" にすれば、毎月１日の日付データが12ヶ月分生成されます。"weeks" だと、2013-01-01, 2013-01-08, 2013-01-15, …… のように展開する日付データが生成されます。
xx <- seq(as.Date("2009-01-01"), as.Date("2013-01-01"), by="years")
とすれば、2009-01-01, 2010-01-01, …… 2013-01-01 の５つの日付データが生成されます。
日付データの説明は、この辺でおわりにしようと思います。
「xx <- Sys.Date()」とすれば、xxに今日の日付(Date型)がセットされることなど、他にもいろいろな話題があります。
また、日付と合わせて時刻も記録・保持するdate-time型に言及してませんが、これは、世界標準時の1970年１月１日 0時0分0秒からの経過秒数を実態とするクラスです。日本標準時は世界標準時より９時間早いので、-32400秒のずれがあります。
「xx <- as.POSIXct("1970-01-01 00:00:00")」とすれば、普通は(環境変数TZがちゃんと設定されていれば)日本標準時の秒数がxxにセットされます。秒数といっても独自クラスなので print(xx) とすれば、「1970-01-01 JST」と表示されます。as.numeric(xx) の戻り値は -32400 です。
一方、「xx <- as.POSIXlt("1970-01-01 00:00:00")」とすれば、xxがlist型になり、xx$year, xx$mon などとして、年・月・日・曜日・時・分・秒などを数値として簡単に取り出すことができます。ただ、1でなく0をスタートの値とするものがあるので注意して下さい。
difftime() は、その名前から分かるように、本来、date-time型の差分を算出するものです。
他に時系列を扱うためのクラス(ts)などもあり、日付や時刻にまつわる話題だけでも多岐にわたりますが、私には手に余るので省略します。
!4. その他
数値データの扱いについて、これまで触れなかったことで気になる点を簡単に書いておきます。
!!(1) 欠損値の扱い
○ 要点１：na.rmオプションでNAを無視
○ 要点２：subset() によって、NAを取り除いたデータフレームを生成
数値データを扱う場合、基本的には欠損値を予め取り除いておくのが無難です。
最大値を求める max(), 最小値を求める min(), 平均値を求める mean() などでは、na.rm というオプションがあります。これは、与えられた引数の中にNAがあるとき、それを無視して計算するか否かを指定するものです。
mean(xx, na.rm=TRUE)
とすれば、NAを無視します。
デフォルトは「na.rm = FALSE」となっているので、NAがあると計算できず、mean() はNAを返します。意味のある答えを得られません。
欠損値を扱う場合、na.rmのようなオプションを利用するのも一つの方法ですが、やはり予め欠損値を取り除いておくのがいいと思います。
仮に data07.csv の「長男の誕生日」に欠損値がある場合、次のようにしてNAを取り除くことができます。
dtf <- read.csv("data07.csv", header=T, na.strings="")
dtf2 <- subset(dtf, !is.na(dtf[["長男の誕生日"]]))
こうすると、「長男の誕生日」の中のNAである行が取り除かれて、その結果がdtf2にセットされます。「長男の誕生日」の列だけが行数を減らすのではなく、すべての列の行数が減ります。
「nrow(dtf) - nrow(dtf2)」とすれば、取り除かれた行が何行あるか、その数を確認できます。
「長男の誕生日」と「次男の誕生日」の両方からNAを取り除いた結果を得たい時は、次のようにします。
dtf3 <- subset(dtf, !is.na(dtf[["長男の誕生日"]]) &
!is.na(dtf[["次男の誕生日"]]))
こうして得られたdtf3には、NAが含まれていないことになります。
数値データの列がもっと沢山あって、そのすべてのNAを削除したい場合は、例えば次のようにします。
nums <- c("長男の誕生日", "次男の誕生日", ……)
xx <- dtf
for (i in nums) xx <- subset(xx, !is.na(xx[[i]]))
上のようにして得られたxxは、指定の列(numsで指定した列)のいずれにもNAがない状態になっています。
列を指定するためのnumsは、「nums <- 2:5」のように、数字で指定しても大丈夫です。
!!(2) 文字列と数値が混在している場合の処理
○ 要点１：as.numeric() によって、数値化できるもののみを数値に変換
○ 要点２：as.Date() によって日付データに変換する時は、ダミーの第１要素を挿入
○ 要点３：grep() によって、パターンにマッチする要素のみを抽出
変数xxに、文字列と数値が混在して入っているとします。例えば次のようなものです。
xx <- c("3", "abc", "1.55", "xyz", "-3.14e-001")
このようなxxに対して as.numeric(xx) を適用して数値化を試みると、次のような結果になります。
3.000 NA 1.550 NA -0.314
つまり、数値化できない "abc" とか "xyz" は、NAになります。
このことを利用すれば、文字列と数値が混在している中から数値を取り出すことができます。
ここで、csvを読み込む場合について考えてみます。
ある列が、上のxxのように混在型だとします。具体的には次のようなcsvです。
ID,data
1,3
2,abc
3,1.55
4,xyz
5,-3.14e-001
このcsvをread.csv() で読み込んで、データフレームdtfにセットしたとします。
このとき、dtf$data は、内容的に前述のxxと同じになるはずです。
ただし、dtf$data は、xxのように文字列から構成されるvectorではなく、factorです。なので、as.numeric(dtf$data) としても、意図したような結果が得られません。
dtf$data は５つの要素から構成されていますが、各要素には整理番号が割り振られています。同じ値の要素には同じ整理番号が付けられますが、今回の例では同じものがないので1〜5のどれかが割り振られます。
as.numeric(dtf$data) とすると、その整理番号がvectorになって返されます。具体的には c(3,4,2,5,1) が返されます。
xxを数値化した時と同じ結果を得るためには、dtf$data を一度 as.character() で文字化して、その上で数値化します。具体的には次のようにします。
nn <- as.numeric(as.character(dtf$data))
このようにすると、nn が c(3, NA, 1.55, NA, -0.314) となります。
念のため、先のcsvを読み込んで数値化するRプログラムを掲げておきます。
−−−− Rプログラムここから
dtf <- read.csv("test.csv", header=T)
nn <- as.numeric(as.character(dtf$data))
print(nn)
−−−− Rプログラムここまで
次に、日付データに話を移します。
xxが文字列から構成されるベクトルの場合、その第１要素が日付データに変換可能なものであれば、as.Date() によってベクトル全体を日付型に変換できます。日付データに変換できる要素は変換され、できない要素はNAに変換されます。
ただし、そのような変換が行われるのは、第１要素が "2013-07-13" のように日付データに変換可能な場合に限られます。"abc" とか "123" などのように変換できない文字列が第１要素だと、エラーが発生します。
「それなら、第１要素として変換可能な文字列をあえて挿入してやり、変換後に第１要素を削除すればいい。」という逆手の発想がうかびます。この方法を採ると次のようなプログラムになります。
xx <- c("123", "2012-6-8", "abc", "1998-10-17")
dd <- as.Date(c("1970-1-1", xx))[-1]
print(dd)
「yy <- xx[-1]」とすると、xxの第１要素が削除されたものがyyにセットされますが、上は、その記法を用いたものです。
上のやり方だと、日付データに変換できない要素がNAになるので、xxとddの要素数が同じになります。
そうではなく、日付データに変換可能なものだけを取り出してDate型にしたい場合は、正規表現とgrep() を用いて次のようにします。
xx <- c("123", "2012-6-8", "abc", "1998-10-17")
ptn <- "[0-9]{4}-[0-9]{1,2]-[0-9]{1,2}"
dd <- as.Date(xx[grep(ptn,xx)])
print(dd)
上のやり方だと、日付データに変換可能なものだけがddにセットされるので、要素数が２になります。
grep(ptn,xx) は、正規表現ptnにマッチするものの要素番号を返します。具体的には c(2,4) を返します。つまり、日付データの書き方にマッチするものが２番目と４番目ということです。
「yy <- xx[c(2,4)]」とすると、xxの２番目の要素と４番目の要素の２つがyyにセットされます。このyyをDate型に変換すれば、目的が達成できます。
なお、このやり方だと、xxの中に日付データに変換可能な要素が１つもない場合、ddには character(0) がセットされます。length(dd) が 0 となります。
--------
以上で「単純集計に関する覚え書き／数値データの集計」をおわりにします。
数値データの処理については、まだまだ様々な事柄がありますが、ごく基本的な部分には言及できたのではないかと思います。
Copyright (C) T. Yoshiizumi, 2013 All rights reserved.
〜 rubyとRを用いた処理 〜
最終更新日: 2013/07/15
以下に掲げるドキュメントは、
[[tjs_c.zip|http://cup.sakura.ne.jp/tjs_c.zip]]
に同梱されている tjs_c.txt と同じ内容です。
＜はじめに＞
「単純集計に関する覚え書き」の第３段として、数値の集計について記します。年齢、身長や体重、価格のように数値として表現できるものであって、平均値とか最大値等を求めることができるものを対象とします。
ちなみに、第１段が[[「単一回答の集計」|http://cup.sakura.ne.jp/hiki/hiki.cgi?tjs_a_guide]]、第２段は[[「複数回答の集計」|http://cup.sakura.ne.jp/hiki/hiki.cgi?tjs_b_guide]]でした。
ソフトウェアとしては R と ruby を用います。拙作 rrxwin(rubyからRを利用するためのライブラリ) も使います。
サンプルのスクリプトを試した環境は次のとおり。
* Windows XP | VISTA
* ruby ver 1.8.7 | 1.9.3
* R ver 2.15.2
--------
{{toc_here}}
--------
!1. 数値データの分布やばらつきを確認
ここでは、「ID・性別・年齢・支出」の４項目からなるcsvデータ data06.csv を取り扱います。100人分のデータです。
IDは、個々人を識別するための整数値です。性別はアルファベットのm, fで記入されており、前者が「男」、後者が「女」です。
年齢と支出が数値データで、どちらも整数値です。支出は、月あたりの支出額のつもりです。乱数発生させた数値に少し手を加えたものです。
以下、このデータを素材にしてサンプルを示します。
!!(1) Rプログラムのsummary関数
○ 要点１：最小値・最大値・平均値等を一括取得するsummary()
○ 要点２：パーセンタイルの考え方による分布の把握
カテゴリー別の集計(男女別の集計)の前に、数値データそのものの分布をみることにします。Rプログラムではsummary関数が用意されており、最小値・最大値・平均値等を一括して取得できます。
csvデータを読み込んだ結果が dtf というデータフレームに収納されている場合、summary(dtf[,"年齢"]) とすれば、年齢に関する分布情報を得られます。
分布情報には次の６つがあります。年齢に関する具体的な値と合わせて記します。
||Min.||最小値||20.00
||1st Qu.||第１の４分位点||28.00
||Median||中央値(第２の４分位点)||34.00
||Mean||平均値||34.29
||3rd Qu.||第３の４分位点||41.00
||Max.||最大値||49.00
上の年齢に関する具体的な値をみると、最も若い人が20歳、最も年長の人が49歳であることが分かります。
平均年齢は 34.29歳です。
第１の４分位点が28なので、ざっくりした言い方ですが、28歳以下が25％いて、28歳以上が75％いることになります。サンプルの総数が100人なので、28歳以下が25人、28歳以上が75人いる(それくらいになる)という目安になります。
実際にデータを数えてみると、25歳「未満」だと24人、25歳「以下」だと27人です。28歳を基準線にして、25％と75％にきちんと厳密に分けられるわけではありません。この辺の詳細について述べるためには「パーセンタイル」の考え方とその計算方法に踏み込まなければなりませんが、ここでは素通りします。
パーセンタイルに関連する値として、中央値が34歳、第３の４分位点が41歳と出ています。34歳以下が50％(50人)くらい、41歳以下が75％くらいいることになります。
20〜49歳の集団としては、34.29歳という平均値は「なるほど」という典型的な値にかなり近いです。中央値(34歳)のところで50％づつに分かれるのも「まさに」という値です。そして、28歳以下が25％、41歳以上が25％となると、特定の年齢層に偏ることなく全体的にばらついていることが推測できます。乱数発生で作ったデータなので当然といえばそれまでですが。
ちなみに、Rプログラムでは、パーセンタイルに関する代表的な値を算出する関数として quantile() というのがあります。使い方は summary() と同じです。
ここで、summary関数を用いるサンプルを掲げておきます。data06.csvの「年齢」と「支出」の２つの項目をsummary() にかけています。
−−−− tjs20.rb ここから
#! ruby -Ks
# summary関数の利用 (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data06.csv", header=T)
items <- c("年齢", "支出")
for (i in items) {
cat(i, "の分布状況\n", sep="")
print(summary(dtf[,i]))
cat("\n")
}
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs20.rb ここまで
!!(2) 平均値と標準偏差
○ 要点１：不偏標準偏差と標本標準偏差の違い
○ 要点２：sapply関数の利用
集団の数値データの代表値として、「平均値」と「標準偏差」を並記する方法があります。平均値だけだと、ばらつきの度合いが分からないので標準偏差(sd:Standard Deviation)を合わせて記します。標準偏差が大きいほど、ばらつきが大きいということになります。
平均値をμ、標準偏差をσとしたとき、おおよそ次のような分布になるようです。
(μ−σ) 〜 (μ＋σ) の範囲に、全構成要素の約68％が入る。
(μ−σ×２) 〜 (μ＋σ×２) の範囲に、全構成要素の約95％が入る。
「標準偏差」を２乗したものを「分散」といいます。というより、「分散」の平方根が「標準偏差」です。
Rプログラムでは、分散を求める関数として var(), 標準偏差を求める関数として sd() があります。どちらも引数として、数値から構成されるベクトルを与えます。
ここで注意しなければならないのは、var(), sd() は、それぞれ不偏分散と不偏標準偏差を算出するものだということです。これらは、標本分散と標本標準偏差と少し違います。
細かな説明は省略しますが、100人のデータから分散を求める場合、計算の最終段階で、100で割り算する方法と、99(100から1を差し引いた値)で割り算する方法の２通りがあります。100で割るのが標本分散、99で割るのが不偏分散です。Rのvar関数は、不偏分散(99で割る方)です。標本分散より少し大きな値になります。
不偏分散(99で割ったもの)の平方根が不偏標準偏差、標本分散(100で割ったもの)の平方根が標本標準偏差です。
Rのvar()が不偏分散を返すのであれば、その値を基にして標本分散を求めるには、99を掛けた上で、改めて100で割ればいいことになります。つまり 99/100(0.99)を掛け算してやります。もちろん 99, 100 という数値は、実際の集団の構成要素数に置き換える必要があります。
不偏分散は、素直に100で割るのではなく、99で割ります。なんでこんなことをするのかという理由について、きちんと説明するのは私には難しいですが、おおよそ次のようなことのようです。
仮に、１万人の集団があったとして、その集団の年齢分布を知りたいとします。お金や時間が潤沢で、１万人全員を調査できたとすれば、集められたデータから分散を求める場合、素直に１万で割り算します。つまり標本分散を使います。標準偏差も標本標準偏差を用います。不偏分散や不偏標準偏差は使いません。
しかし、１万人全員の調査をやる余裕がなく、100人だけ無作為抽出して調査した場合は、１万人分の状況を推測する際、100で割り算するのでなく 99で割ります。100で割るよりも99で割った方が、１万人の母集団の分散に近い、というのが統計学の知見のようです。
要するに、部分から全体を推測する時は「不偏」の分散と標準偏差を用い、調査した集団そのものだけに着目する時は(つまり、より大きな集団について推定したりしない場合)「標本」の分散と標準偏差を用いるということのようです。
平均と標準偏差を並記する場合、適当な方の標準偏差を選んで記すということになるのだと思います。
なお、Wikipediaの「標準偏差」のところを読むと、「不偏」と「標本」の表現の使い分けが必ずしも統一されているわけではないようです。ご注意ください。
話を戻します。
Rプログラムで、標本分散や標本標準偏差を算出するには、例えば次のようにします。xが数値データからなるベクトルだとします。
n <- length(x) # xの要素個数
var100 <- var(x)*(n-1)/n # 標本分散
sd100 <- sqrt(var100) # 標本標準偏差
もしベクトルxの中にNAが含まれている可能性があるのであれば、それを取り除いた上で計算するのが無難です。最初の方に次の１行を挿入します。
x <- x[!is.na(x)] # NAを除去
下にサンプルを掲げてみます。「年齢」と「支出」について、平均と標本標準偏差を出力します。
−−−− tjs21.rb ここから
#! ruby -Ks
# 平均と標本標準偏差 (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data06.csv", header=T)
items <- c("年齢", "支出")
for (i in items) {
x <- dtf[,i]
x <- x[!is.na(x)]
n <- length(x)
sd100 <- sqrt(var(x)*(n-1)/n)
cat(sprintf("%sの平均値: %.1f, 標本標準偏差: %.3f\n", i, mean(x), sd100))
}
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs21.rb ここまで
ここで少し横道にそれますが、上の tjs21.rb の出力結果を表形式で出力することを考えてみます。
my.test() という関数を定義して、与えられたベクトルから平均と標本標準偏差を算出するようにします。
また、Rプログラムのapply関数群の１つである sapply() を用います。
出力結果を先に示すと次のとおり。
|| ||年齢||支出
||平均||34.290000||99673.78
||標本標準偏差||7.750219||28716.33
上の出力を得るためのスクリプトは下のとおりです。
−−−− tjs21b.rb ここから
#! ruby -Ks
# 平均と標本標準偏差(sapply利用) (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
my.test <- function(x) { # mean, var, sd の算出
x <- x[!is.na(x)]
n <- length(x)
m <- mean(x) # 平均値
var99 <- var(x) # 不偏分散
sd99 <- sd(x) # 不偏標準偏差
var100 <- var99 * (n-1) / n # 標本分散
sd100 <- sqrt(var100) # 標本標準偏差
return(c(平均=m, 標本標準偏差=sd100))
}
dtf <- read.csv("data06.csv", header=T)
items <- c("年齢", "支出")
xx <- sapply(items, function(i) my.test(dtf[,i]))
xx
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs21b.rb ここまて
上のスクリプトにおいて、my.test() の中で「不偏」と「標本」の２種類の分散と標準偏差を一通り算出していますが、戻り値として返しているのは「平均」と「標本標準偏差」の２つだけです。
sapply() という関数を使っていますが、これについてちゃんと説明する自信がないので、同じことを for を使って書くとどうなるか示しておきたいと思います。該当部分のみ掲げます。
−−−− sapplyの代替 ここから
items <- c("年齢", "支出")
xx <- NULL
for (i in items) xx <- cbind(xx, my.test(dtf[,i]))
colnames(xx) <- items
xx
−−−− sapplyの代替 ここまで
sapply() は、第１引数 items の各々の要素 i について、function(i) を適用するということだと思います。そして、その function の実態は、my.test(dtf[,i]) であるというわけです。
ちなみに、function(i) の定義は、sapply() の引数指定の中に長々と盛り込むことができます。my.test() の中身をsapply() の引数の一部として書き込むことができます。なので、my.test関数の定義は不要になります。例えば次のように書きます。
−−−− sapply関数の記述例 ここから
xx <- sapply(items, function(i) {
x <- dtf[,i]
x <- x[!is.na(x)]
n <- length(x)
m <- mean(x) # 平均値
sd100 <- sqrt(var(x) * (n-1) / n) # 標本標準偏差
return(c(平均=m, 標本標準偏差=sd100))})
print(xx)
−−−− sapply関数の記述例 ここまで
apply関数群は、他のところで使うことになるので、ちょっと横道にそれましたが取り上げてみました。
!!(3) 数値データのカテゴリー化
○ 要点１：cut() によって数値データにカテゴリー名を割り当てる
○ 要点２：seq() によって等差数列を生成
data06.csvの「年齢」は、20〜49歳の範囲にあります。これを「20代・30代・40代」の３分類でカテゴリー化することを考えます。そして、各々のカテゴリーに属する人の人数を出力してみます。
まず、出力結果は次のようになります。
年齢層別の人数
20代 30代 40代
32 36 32
この出力を得るためのスクリプトは、下のとおりです。
−−−− tjs22.rb ここから
#! ruby -Ks
# 数値データのカテゴリー化(cut関数の利用) (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data06.csv", header=T)
xx <- cut(dtf[,"年齢"], breaks=c(20, 30, 40, 50),
labels=c("20代", "30代", "40代"),
right=FALSE)
table(年齢層別の人数=xx)
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs22.rb ここまで
data06.csvの「年齢」の列をcut() に引きわたして、各々の人の年齢にカテゴリー名を割り当てます。年齢が20以上・30未満なら「20代」、30以上・40未満なら「30代」というカテゴリー名です。
cut() の戻り値 xx には100人分のデータが入っています。28 とか 36 とか 43 といった具体的な年齢ではなく、「20代・30代・40代」の３種類のどれかです。それが100人分収納されています。
といっても、xx は factor なので、カテゴリー名が文字列の形で記録されているわけではなく、整数値の 1, 2, 3 が100人分記録されています。数字の1には「20代」、2には「30代」、3には「40代」が割り振られています。
cutの引数 breaks は、年齢をカテゴリー区分する時の境目となる値を指定するものです。ベクトルで与えます。
サンプルでは c(20, 30, 40, 50) となっているので、20以上・30未満、30以上・40未満、40以上・50未満という区分になります。
このとき、１つの区間について下限の値は含むけれども上限の値は含まない、つまり「20以上・30未満」の20は含むけれども30は含まない、となるのは、常にそうであるわけではなく、cut() の引数 right を FALSE にしているためそうなるものです。
サンプルには出てきていませんが、include.lowest というオプションがあり、それと rightオプションの組合せによって、下限と上限を含めるかどうか指定できます。この辺の話は少々ややこしいので省略します。
breaksオプションとして与えるベクトルにおいて、全体の下限(20)が判然としない時は、マイナス無限大の -Inf を与えます。また、全体の上限(50)が判然としない時は、無限大の Inf を与えれば大丈夫だと思います。
labelsオプションは、カテゴリー名を列挙したベクトルを与えます。factor() の labels と同じような意味合いです。
labelsとして NULL を指定するか、あるいは省略すると、「[20,30)」 「[30,40)」 「[40,50)」というカテゴリー名が割り当てられます。
ところで、「20代・30代・40代」の３分類だと大雑把すぎると思われるかもしれません。そこで、５歳刻みで分類するケースを考えてみます。
breaksオプションとして c(20, 25, 30, …… 50) と記述するのは面倒なので、等差数列を生成する seq() を用います。
brk <- seq(20, 50, by=5)
とすれば、brkに５歳刻みのベクトルがセットされます。seq() の引数には、整数値だけでなく実数値も使えるので、いろいろ応用できると思います。
labelsオプションとして与えるベクトルの生成には、先述した sapply() を用いてみます。次のようにします。
lbl <- sapply(1:(length(brk)-1), function(i)
sprintf("%s〜%s", brk[i], brk[i+1]))
上の２点を盛り込んだスクリプト tjs23.rb を念のため下に掲げておきます。
−−−− tjs23.rb ここから
#! ruby -Ks
# 数値データのカテゴリー化(seq関数の利用) (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data06.csv", header=T)
brk <- seq(20, 50, by=5)
lbl <- sapply(1:(length(brk)-1), function(i)
sprintf("%s〜%s", brk[i], brk[i+1]))
xx <- cut(dtf[,"年齢"], breaks=brk, labels=lbl, right=FALSE)
table(年齢層別の人数=xx)
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs23.rb ここまで
seq() に関して少し補足します。
20〜50歳の範囲を５歳刻みなどではなく、４等分したいといったことがあると思います。４等分だと「20〜27.5、27.5〜35、35〜42.5、42.5〜50」となります。
この場合、seq() の byオプションではなく、length.outオプションを用います。
brk <- seq(20, 50, length.out=5)
とすれば、４等分にできます。length.out の値として、4ではなく5を指定する点に注意して下さい。nを指定すると「n-1」等分されます。
実際の年齢は整数値なので、27.5 とか 42.5 などの小数点数だと違和感があるし、「以上・未満」の境目も分かりにくくなるので、まるめ処理をして整数にしたいという場合は、次の関数を使えます。
floor(x) 少数点以下を切り捨て
ceiling(x) 少数点以下を切り上げ
round(x,0) 少数点以下を四捨五入
これら３つの関数とも、引数にベクトルを与えることができます。ベクトルの要素すべてについて、切り捨て・切り上げ・四捨五入を施します。
最後のroundの四捨五入は、単純な四捨五入ではなく、JISやISOで定められた方式に従った四捨五入のようなので注意して下さい。これに関しては次のサイトが参考になります。
[[JIS,ISO式四捨五入 - RjpWiki|http://www.okada.jp.org/RWiki/?JIS%2CISO%BC%B0%BB%CD%BC%CE%B8%DE%C6%FE]]
数値データの分布を確認する方法は、この辺にしておきます。
--------
!2. 数値データとカテゴリーデータの組合せによる集計
data06.csvには「年齢」と「支出」という数値データの他に、「性別」というカテゴリーデータがあります。
ここでは、男女別に年齢や支出の分布をみたり、カテゴリー化した「年齢」や「支出」と「性別」とのクロス集計を取り上げます。
既に述べた「1. 数値データの分布を確認」と、当連載の第１回目「単一回答の集計」のノウハウを組み合わせたものになります。
!!(1) １つの値(平均値)をカテゴリー別に算出
○ 要点：tapply() によって、分類しつつ関数を適用
data06.csvを読み込んだ結果が dtf に記録されているものとします。dtf はデータフレームです。
男女別に「年齢」と「支出」の平均値を得ようと思ったら、まず dtf の中から「男」に関するデータを male にセットし、「女」に関するデータを female にセットして、その上で各々の平均値を求めることになります。
これを比較的素直に実行しているのが下のプログラムです。実行結果とRプログラムの部分を示します。
|| ||年齢||支出
||男||34.38596||104923.07
||女||34.16279||92715.42
−−−− Rプログラムここから
dtf <- read.csv("data06.csv", header=T)
dtf$性別 <- factor(dtf$性別, levels=c("m","f"), labels=c("男","女"))
nums <- c("年齢", "支出") # 処理する数値欄の列名
male <- subset(dtf, 性別 == "男")
female <- subset(dtf, 性別 == "女")
result <- NULL
for (i in nums) {
m <- mean(male[,i]) # 「男」の平均値
f <- mean(female[,i]) # 「女」の平均値
result <- cbind(result, c(m,f))
}
colnames(result) <- nums
rownames(result) <- c("男", "女")
result
−−−− Rプログラムここまで
「dtf$性別 <- factor(……)」は、「性別」の欄に記入されているのが m, f なので、それを「男・女」に置換・調整するものです。
上のプログラムは、別のケースに応用しやすい形とはいえません。例えば、男女別ではなく地域別のデータがあった場合、それに対応させるためには、かなり書き換える必要が出てきます。地域別は、２分類とは限らず５分類とか10分類かもしれません。「男・女」という２分類前提の書き方をやめなければなりません。
levels(dtf[,"性別"]) が c("男","女") というベクトルを返すことを利用して、少し工夫してみたのが下のスクリプトです。分かりにくく、かつ、短くできたわけでもないので工夫のしがいがないですが、参考まで掲げてみます。
−−−− tjs24.rb ここから
#! ruby -Ks
# 男女別の年齢・支出の平均値 (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data06.csv", header=T)
dtf$性別 <- factor(dtf$性別, levels=c("m","f"), labels=c("男","女"))
key <- "性別" # 分類の手がかりとなる列の名前
nums <- c("年齢", "支出") # 処理対象の数値欄の列名
ctg <- levels(dtf[,key]) # == c("男", "女")
list.dat <- list() # 空のリスト
for (i in ctg) # dtfを男・女それぞれに分けてlist.datに追加
list.dat <- c(list.dat, list(subset(dtf, dtf[,key] == i)))
result <- NULL
for (i in nums) {
v <- NULL
for (j in 1:length(list.dat)) {
dd <- list.dat[[j]] # 男または女それぞれのデータフレームをddにセット
v <- append(v, mean(dd[,i])) # 年齢・支出それぞれの平均をvに追加
}
result <- cbind(result, v)
}
colnames(result) <- nums
rownames(result) <- ctg
result
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs24.rb ここまで
空のリスト list.dat には、１番目に「男」のデータフレーム、２番目に「女」のデータフレームを入れます。なので、list.dat［［1］］ が「男」のデータフレーム、list.dat［［2］］ が「女」のデータフレームです。
あとは、「年齢」と「支出」の平均を男・女それぞれについて求めているだけです。
さて、ここからが本論です。apply関数群を使って、上のプログラムを短くします。
まず、tapply() に着目します。これは、factorである列を含むデータフレームがあるとき、factorの要素の種類(性別欄の男・女)を手がかりにして分類しながら、指定の列(年齢・支出)を指定の関数(mean)に引きわたして、その結果を返します。
例えば次のように記述します。
xx <- tapply(dtf[,"年齢"], dtf[,"性別"], mean)
こうして得られたxxは、34.38596, 34.16279 という２つの要素からなる配列(array)です。１番目の要素が「男」の平均年齢、２番目が「女」の平均年齢です。これらは名前付きで記録されているので「nn <- names(xx)」とすれば、c("男","女") が得られます。
この tapply() を「年齢」と「支出」のそれぞれについて呼び出して、その２つの結果をマトリクスにすれば、前述のプログラムと同じ結果が得られるはずです。
tapply() を２回分記述するところに、for を用いてもいいのですが、既に言及した sapply() を用いると、更に簡潔になります。
ということで、だらだら説明してきましたが、その「だらだら」を踏まえて下のようなスクリプトを作りました。分かりやすいとはいえませんが、短くなったのは確かです。
−−−− tjs25.rb ここから
#! ruby -Ks
# 男女別の年齢・支出の平均値(tapply利用) (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data06.csv", header=T)
dtf$性別 <- factor(dtf$性別, levels=c("m","f"), labels=c("男","女"))
key <- "性別" # 分類の手がかりとなる列の名前
nums <- c("年齢", "支出") # 処理対象の数値欄の列名
result <- sapply(nums, function(i) tapply(dtf[,i], dtf[,key], mean))
result
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs25.rb ここまで
上のスクリプトが、前述した tjs24.rb よりも作業効率がいいのかどうかは分かりません。
「年齢」の平均値を出す時と、「支出」の平均値を出す時の各々で、男女別に分類する処理を実行しているような気がします。つまり、男女別の分類処理を２回実行しているのではないかと思います。
一方、長いスクリプトの tjs24.rb の方は、男女別の分類処理は１回のみです。分類した結果をlist.datにセットして、それを用いて「年齢」と「支出」の平均値を算出しています。
Rプログラムの内部処理がどうなっているか分からないので、上に書いたことが正しいかどうか自信はありません。また、もし正しいとしても、それが作業効率に影響するのかどうか、それも分かりません。データの規模が大きくなるほど作業効率が問題になるので、気になるところではありますが、ちゃんと確認していません。
なお、男女別の他に、男女を合わせた「全体」も付けた形の表を得たい時は、例えば次のようにします。もっと適切なやり方があるような気はしますが、私が思いついたものを掲げておきます。
−−−− tjs25b.rb ここから
#! ruby -Ks
# 男女別の年齢・支出の平均値(全体欄を付加) (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data06.csv", header=T)
dtf$性別 <- factor(dtf$性別, levels=c("m","f"), labels=c("男","女"))
key <- "性別" # 分類の手がかりとなる列の名前
nums <- c("年齢", "支出") # 処理対象の数値欄の列名
xx <- sapply(nums, function(i) tapply(dtf[,i], dtf[,key], mean))
ww <- sapply(nums, function(i) mean(dtf[,i]))
result <- rbind(xx, 全体=ww)
result
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs25b.rb ここまで
上のスクリプトでは、「全体」を別枠で算出して、後から男女別のものに付け加えています。
これを別枠ではなく「男・女・全体」を一括して算出する方法があるような気がしますが、分かりませんでした。
それから、平均値を示す場合、その算出に用いたサンプルの数が気になります。「男」の人数、「女」の人数、「全体」の人数がどれくらいかです。
男女それぞれの人数は、「nn <- table(dtf[,key])」で求めることができます。単一回答や複数回答の集計のところで多用した方法です。nnは、２つの要素からなる配列で、１番目が「男」の人数、２番目が「女」の人数です。
また、男女別の人数を「nn <- tapply(dtf[,key], dtf[,key], length)」として求めることもできます。
「全体」の人数は、length(dtf[,key]) で取得できます。
ここには掲げませんが、最も右の列に「人数」を加えて表示するスクリプト tjs25c.rb を圧縮ファイルに同梱してあるので、よかったら参考にして下さい。
!!(2) 複数の値(最小値・平均値・最大値等)をカテゴリー別に算出
○ 要点：lapply() の利用
平均値という１つの値だけを算出するのであれば、前述した tjs24.rb, tjs25.rb のやり方でいいわけですが、最小値・平均値・最大値といった複数の値の算出を考えると、もう少し工夫してやる必要があります。
平均値を求める mean() ではなく、分布情報を得るための summary() を用いる場合、「年齢」と「支出」の各々について、分布情報が男女間でどのように違うかをみるためには、次のような出力結果が必要になります。
||年齢
|| ||男||女
||Min.||20.00||21.00
||1st Qu.||28.00||27.50
||Median||35.00||33.00
||Mean||34.39||34.16
||3rd Qu.||41.00||41.00
||Max.||49.00||49.00
||支出
|| ||男||女
||Min.||50430||52490
||1st Qu.||82250||70020
||Median||113600||92460
||Mean||104900||92720
||3rd Qu.||129400||116900
||Max.||149400||148400
平均値だけを求める tjs25.rb では、「result <- sapply(……)」として、最終的な結果を sapply() で取得しました。
この sapply は、戻り値としてベクトルやマトリクスを返します。中身の計算結果が複雑化したとしても、それをなるべく単純なベクトルやマトリクスに変換して返します。
しかし、mean ではなく summary だと、複数の値を算出するので、sapply では適切な結果を得られません。このような時は lapply() を用います。lapply() は、Rプログラムでいうところのリストを返します。なので、複雑な階層構造のデータでも意図どおりの形で取得できます。使い方は、基本的に sapply() と同じです。
tjs25.rb の sapply を lapply に変更したプログラムを下に掲げてみます。Rプログラム部分です。
−−−− Rプログラムここから
dtf <- read.csv("data06.csv", header=T)
dtf$性別 <- factor(dtf$性別, levels=c("m","f"), labels=c("男","女"))
key <- "性別" # 分類の手がかりとなる列の名前
nums <- c("年齢", "支出") # 処理対象の数値欄の列名
result <- lapply(nums, function(i) tapply(dtf[,i], dtf[,key], summary))
names(result) <- nums
result
−−−− Rプログラムここまで
上のプログラムを実行した結果は、次のとおりです。「年齢」に関する部分のみ示します。
||$年齢
||$年齢$男
|| ||Min.||1st Qu.||Median||Mean||3rd Qu.||Max.
|| ||20.00||28.00||35.00||34.39||41.00||49.00
||$年齢$女
|| ||Min.||1st Qu.||Median||Mean||3rd Qu.||Max.
|| ||21.00||27.50||33.00||34.16||41.00||49.00
上の出力結果だと男女間の比較がやりにくいので、もう少し工夫してみます。「年齢」と「支出」の各々について、「男・女」と「最小値・平均値・最大値等」とのクロス表を作るようにします。
出力結果は次のとおり。
||$年齢
|| ||Min.||1st Qu.||Median||Mean||3rd Qu.||Max.
||男||20||28.0||35||34.39||41||49
||女||21||27.5||33||34.16||41||49
||$支出
|| ||Min.||1st Qu.||Median||Mean||3rd Qu.||Max.
||男||50430||82250||113600||104900||129400||149400
||女||52490||70020||92460||92720||116900||148400
上のように表示すれば、男女間の違いを比較しやすくなります。
この表示結果を得るためのスクリプトは下のとおりです。
−−−− tjs26.rb ここから
#! ruby -Ks
# 年齢・支出の男女別分布情報(lapply利用) (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data06.csv", header=T)
dtf$性別 <- factor(dtf$性別, levels=c("m","f"), labels=c("男","女"))
key <- "性別" # 分類の手がかりとなる列の名前
nums <- c("年齢", "支出") # 処理対象の数値欄の列名
xx <- lapply(nums, function(i) tapply(dtf[,i], dtf[,key], summary))
names(xx) <- nums
ctg <- levels(dtf[,key])
result <- lapply(nums, function(i) t(sapply(ctg, function(j) xx[[i]][[j]])))
names(result) <- nums
result
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs26.rb ここまで
上のスクリプトで、lapply() が２回出てきますが、１回目は tjs25.rb の sapply を lapply に変更しただけです。
２回目の lapply() の引数部分が分かりにくいかもしれません。この部分では、変数iに「年齢」と「支出」が逐次セットされ、変数jには「男」と「女」が逐次セットされます。
まず「年齢」に関して、summary() の出力結果を取りまとめます。「男」と「女」がばらばらになっていたものを１つの表にします。そして、「支出」に関しても同じように処理します。
sapply(……) を t() で囲んでますが、この t() は、表の縦軸と横軸を入れ換えるものです。今回のケースでは、これがないと縦に長い表になるので入れてみました。
参考まで、男女別の他に「全体」の欄を付加するスクリプトを下に掲げておきます。
−−−− tjs26b.rb ここから
#! ruby -Ks
# 年齢・支出の男女別分布情報(全体欄を付加) (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data06.csv", header=T)
dtf$性別 <- factor(dtf$性別, levels=c("m","f"), labels=c("男","女"))
key <- "性別" # 分類の手がかりとなる列の名前
nums <- c("年齢", "支出") # 処理対象の数値欄の列名
xx <- lapply(nums, function(i) tapply(dtf[,i], dtf[,key], summary))
names(xx) <- nums
ctg <- levels(dtf[,key])
yy <- lapply(nums, function(i) sapply(ctg, function(j) xx[[i]][[j]]))
names(yy) <- nums
ww <- lapply(nums, function(i) summary(dtf[,i]))
names(ww) <- nums
result <- lapply(nums, function(i) t(cbind(yy[[i]], 全体=ww[[i]])))
names(result) <- nums
result
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs26b.rb ここまで
ここには掲げませんが、上のスクリプトと同じ結果をタブ区切りテキストの形で出力する tjs26c.rb を圧縮ファイルに同梱しておきます。
男女別の結果(変数xx)と全体の結果(変数ww)をrubyの側に引きわたして、ruby側でその結合と出力を行います。
すべてRプログラムで記述する場合よるも長くなりますが、引きわたされたデータがどのような形になっているのかを確認する材料になります。よかったら覗いてみて下さい。
なお、tjs26c.rb を動かす際、rrxwin ver 1.02 が必要です。ver 1.01 までだと正常に動かないので注意して下さい。
!!(3) 数値データをカテゴリー化した上でクロス集計
○ 要点：数値データをcut()でカテゴリー化した上でクロス集計
数値データのカテゴリー化は、前述したように cut() で行うことができます。カテゴリー化すれば、他のカテゴリーデータと組み合わせてクロス集計できます。
例えば、「性別」と「年齢」のクロス集計を出力するRプログラムは次のとおりです。「年齢」を「20代・30代・40代」の３分類にします。
−−−− Rプログラムここから
dtf <- read.csv("data06.csv", header=T)
dtf$性別 <- factor(dtf$性別, levels=c("m","f"), labels=c("男","女"))
brk <- c(20, 30, 40, 50)
lbl <- c("20代", "30代", "40代")
age <- cut(dtf[,"年齢"], breaks=brk, labels=lbl, right=FALSE)
tt <- table(dtf[,"性別"], age)
names(dimnames(tt)) <- c("性別", "年齢層")
tt
−−−− Rプログラムここまで
上のプログラムに少し拡張性を持たせてみます。
「性別」 「年齢」 「支出」の３つをそれぞれ組み合わせてクロス集計するスクリプトを下に掲げます。
「年齢」は「20代・30代・40代」の３分類、「支出」は５万〜15万を２万刻みで分類します。
−−−− tjs27.rb ここから
#! ruby -Ks
# カテゴリー化した数値データを含むクロス集計 (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data06.csv", header=T)
dtf$性別 <- factor(dtf$性別, levels=c("m","f"), labels=c("男","女"))
item <- c("性別", "年齢", "支出")
dtf2 <- data.frame(dtf[,item[1]]) # 「性別」の列をdtf2に入れる
brk <- c(20, 30, 40, 50)
lbl <- c("20代", "30代", "40代")
xx <- cut(dtf[,item[2]], breaks=brk, labels=lbl, right=FALSE)
dtf2 <- cbind(dtf2, xx) # カテゴリー化した「年齢」をdtf2に追加
brk <- seq(50000, 150000, by=20000)
lbl <- sapply(1:(length(brk)-1), function(i)
sprintf("%s〜%s万", brk[i]/10000, brk[i+1]/10000))
xx <- cut(dtf[,item[3]], breaks=brk, labels=lbl, right=FALSE)
dtf2 <- cbind(dtf2, xx) # カテゴリー化した「支出」をdtf2に追加
colnames(dtf2) <- item
for (i in 1:(length(item)-1))
for (j in (i+1):length(item)) {
tt <- table(dtf2[,i], dtf2[,j])
names(dimnames(tt)) <- c(item[i], item[j])
print(tt)
cat("\n")
}
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs27.rb ここまで
上のスクリプトを実行すると、３つの表が出力されます。いずれも単一回答を組み合わせたクロス集計表とみなすことができます。合計欄を付加したり、パーセントの表にしたりする場合は、「単一回答の集計」で取り上げたノウハウを利用できます。
--------
!3. 日付データの取扱い
これまで取り上げてきた数値データは、「年齢」と「支出」の２つでした。どちらも整数値ですが、取り上げた処理ノウハウは実数値にも適用できます。
ただ、数値データの１種である日付データ(あるいは日付・時刻データ)を扱う場合は、それなりの工夫が必要になります。その詳細には立ち入りませんが、関連のサンプルをいくつか示してみたいと思います。
扱うcsvデータ data07.csv は、「ID・長男の誕生日・次男の誕生日」の３列から構成されています。100組のデータが入っています。次のような内容です。
ID,長男の誕生日,次男の誕生日
1,1988-10-30,1997-04-13
2,1992-02-04,1997-07-22
3,1995-11-22,1998-05-30
!!(1) 日付の分布状況を確認
○ 要点１：日付データは Dataオブジェクトに変換してから処理
○ 要点２：Dateオブジェクトを文字列に変換するのは format()
Rプログラムにおいて、data07.csv を read.csv() で読み込むと、「1992-02-04」などの日付が文字列として認識され、データフレームになった段階では factor になります。
dtf <- read.csv("data07.csv", header=T)
として読み込むと、dtf$長男の誕生日, dtf$次男の誕生日 のどちらもfactorです。
このままだと日付として扱うのは難しく、例えば、その日付の曜日を簡単に求めることができません。長男と次男の誕生日が何日隔たっているかといった計算もできません。
そこで、このfactorをDate型に変換します。次のようにすると変換できます。
dtf$長男の誕生日 <- as.Date(dtf$長男の誕生日)
dtf$次男の誕生日 <- as.Date(dtf$次男の誕生日)
こうすると、Date型になり、その実態は、1970年１月１日からの経過日数になるようです。つまり数値になります。もっと前の日付だとマイナスの値になります。
1970年１月１日の値が0、翌日の値が1、翌々日の値が2です。
といっても、xxがDate型であるとき、print(xx) で表示されるのは「1970-01-01」などの文字形式です。数値として表示させたい時は、print(as.numeric(xx)) とか print(as.vector(xx)) とします。
ちなみに、数値をDate型に変換する時は、「as.Date(c(0,7,14), origin="1970-01-01")」のように、originオプションによって起点となる日付を指定する必要があります。
日付を文字で表現する場合、「2013-07-01」のように、区切り文字としてマイナス記号を用いるのが基本のようですが、スラッシュを用いてもトラブルなくDate型に変換できるようです。そして、Date型に変換したものをprint() で表示すると、スラッシゅではなくマイナスで表示されます。
また、「2013-07-01」のように、月や日を２桁にするため0を挿入するのが一般的かもしれませんが、「2013-7-1」のように0を省略しても、Date型に変換すれば同じ数値になります。そして、print() で表示すると、0が挿入されて「2013-07-01」のようになります。
なお、「年」を省略して "05-15" などのように「月・日」だけの文字列は、Date型に変換しようとしてもエラーになります。今年の西暦年を自動的に補ってくれる、ということはないようです。
Date型は、基本的に数値なのでいろいろな関数で扱うことができます。例えば、分布状況を表示する summary() もその１つです。最小値(最も早い誕生日)、最大値(最も遅い誕生日)などを確認できます。
下にサンプルスクリプトを掲げておきます。
−−−− tjs28.rb ここから
#! ruby -Ks
# Date型への変換とsummaryの適用 (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data07.csv", header=T)
dtf$長男の誕生日 <- as.Date(dtf$長男の誕生日)
dtf$次男の誕生日 <- as.Date(dtf$次男の誕生日)
summary(dtf$長男の誕生日)
summary(dtf$次男の誕生日)
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs28.rb ここまで
上の実行結果は次のとおりです。長男のもののみ示します。横に長くなるので縦・横逆にして掲げます。
||Min.||"1988-08-25"
||1st Qu.||"1990-06-07"
||Median||"1992-08-09"
||Mean||"1992-07-21"
||3rd Qu.||"1994-09-04"
||Max.||"1996-12-18"
ここで、このDate型のsummary() の結果をrubyの側で受け取ることを考えます。
実は、単純にsummaryの結果を引きわたすと、「1996-12-18」のような文字列ではなく、9848のような数値になります。この数値は、1970年１月１日からの経過日数なので、rubyの側で日付に変換することができます。しかし、Rプログラムで算出されたものは、Rで文字列に変換した方が無難です。
Date型の数値を文字列に変換するのは、format() 関数で行います。例えば次のように呼び出します。
format(xx, "%Y-%m-%d")
こうすると、Date型の数値が「1996-12-18」といった文字列に変換されます。
以下に、rubyにデータを引きわたすサンプルを掲げます。
−−−− tjs28b.rb ここから
#! ruby -Ks
# Date型データをruby側に渡す (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data07.csv", header=T)
item <- c("長男の誕生日", "次男の誕生日")
for (i in item) {
dtf[[i]] <- as.Date(dtf[[i]])
xx <- format(summary(dtf[[i]]), "%Y-%m-%d")
xx$memo <- i
robj(xx)
}
EOS
##
hs = Rrx.rexec(rpro)
n = 0
while (xx = hs["robj#{n+=1}"]) != nil
memo = xx.delete("memo")
ary = []
xx.each do |key, val|
ary << [key, val]
end
puts memo
print Rrx.ary2str(ary, "\t")
printf("\n")
end
−−−− tjs28b.rb ここまで
余談ですが、Rプログラムのsummary() が返す値の構造は、ちょっと特殊です。単純な matrix とか table ではありません。
そして、summary() が返す値に format() を適用すると、構造が変化します。少し単純になります。rubyで受け取ると、素直なハッシュになっています。
summary() が返す値そのままだと、rubyで受け取ったとき、もっと複雑な構造なので注意して下さい。
!!(2) ２つの日付の隔たりを調べる
○ 要点：二つのDateオブジェクトの隔たりはdifftime()で取得
長男の誕生日と次男の誕生日がどれくらい隔たっているかを調べる場合、difftime() を用います。
dd <- difftime(dtf$次男の誕生日, dtf$長男の誕生日, units="days")
とすれば、各々の隔たりの期間が変数ddにセットされます。ddには100組のデータが収納されます。
unitsオプションで "days" を指定しているので、とりあえず隔たりが日数として記録されますが、ddの中身は、単なる日数(整数値)ではなく、difftimeという独自のクラスです。
unitsオプションには、"secs", "mins", "hours", "days", "weeks" が指定できます。サンプルでは日付・時刻のdate-time型ではなくDate型を扱っているので、"days", "weeks" のどちらかということになるでしょうか。
下に difftime() の利用例を掲げます。最も隔たりの大きい兄弟と、最も隔たりの小さい兄弟を表示します。
−−−− tjs29.rb ここから
#! ruby -Ks
# 誕生日の隔たりを算出 (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
dtf <- read.csv("data07.csv", header=T)
dtf$長男の誕生日 <- as.Date(dtf$長男の誕生日)
dtf$次男の誕生日 <- as.Date(dtf$次男の誕生日)
dd <- difftime(dtf$次男の誕生日, dtf$長男の誕生日, units="days")
dtf <- cbind(dtf, 隔たり=dd)
dmax <- max(dd)
dmin <- min(dd)
xx <- subset(dtf, 隔たり == dmax)
cat("最も隔たりの大きい兄弟\n")
xx
cat("\n")
xx <- subset(dtf, 隔たり == dmin)
cat("最も隔たりの小さい兄弟\n")
xx
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs29.rb ここまで
上を実行すると、隔たりが日数で表示されます。これを「週」の数で示したい時は、difftimeオブジェクト「xx$隔たり」のunitsを次のようにして変更します。
units(xx$隔たり) <- "weeks"
こうすると、「460.7143 weeks」とか「59.57143 weeks」のように表示されます。
残念ながら「年」単位と「月」単位の指定はできないので、少し手間をかけて算出することになります。formatの書式指定(%Y, %m)によって「年」と「月」を文字として得られるので、それを数値に変換します。そうすれば、何年何ヶ月の隔たりかを計算できます。
xxに「最も隔たりの大きい兄弟」のデータが入っている場合、次のようにします。
y1 <- as.numeric(format(xx$長男の誕生日, "%Y"))
y2 <- as.numeric(format(xx$次男の誕生日, "%Y"))
m1 <- as.numeric(format(xx$長男の誕生日, "%m"))
m2 <- as.numeric(format(xx$次男の誕生日, "%m"))
yy <- y2 - y1
mm <- m2 - m1
if (mm < 0) {
mm <- 12 + mm
yy <- yy - 1
}
上のようにして求めた yy, mm は、「yy年mmヶ月の隔たり」ということになります。
ここには掲げませんが、長男と次男の誕生日の隔たりを「月数差」と「年月差」で求めて、最も隔たりが大きい組と、最も隔たりの小さい組を表示するスクリプト tjs29b.rb を圧縮ファイルに同梱したので、よかったら参考にして下さい。
!!(3) 日付データのカテゴリー化(年・月・日・四半期・曜日)
○ 要点１：cut() によって日付データを四半期単位で分類
○ 要点２：unlist(strsplit("2013-01-01", "-"))によって文字列分割
数値データをcut() を用いてカテゴリー化する方法については既に触れました。
日付データの「月」を数値で取得すれば、四半期単位(1〜3月、4〜6月、7〜9月、10〜12月の４分類)にするのは容易です。
xxがベクトルで、100個の日付データが入っている場合でいうと、次のようにして四半期単位の分類を行えます。
mm <- as.numeric(format(xx, "%m")) # 「月」の数を抽出
qq <- as.numeric(cut(mm, breaks=c(1,4,7,10,13), labels=1:4, right=F))
上のようにして求めた qq には100個の数値が入っています。その数値は 1〜4 のどれかです。第１四半期が1、第２四半期が2などです。
cut() が返す値は factor なので、それを as.numeric() で数値に変換します。
xxとqqを対応づければ、該当の日付が第何四半期に当たるかを確認できます。
日付データから「月」を取得する関数 months() もありますが、これが返す値は、"1月" とか "4月" などの文字列です。
xxが複数の日付データからなるベクトルのとき、「mm <- months(xx)」とすれば、mmは、「月」を表す文字列が複数入ったベクトルになります。factorではありません。
この mm を材料にして四半期単位の分類を行うこともできますが、「月」を数値にした方が分類しやすいので、先述の方法を採ることにします。
なお、weekdays() を使うと、「月」ではなく「曜日」に変換することができます。これも months() と同じように、文字列からなるベクトルを返します。
ここで、四半期単位での分類だけでなく、「年・月・日・四半期・曜日」に分けるための関数 date.split() を考えてみます。
ss が "2013-01-01" のような文字列である場合、「strsplit(ss, "-")」によって、マイナス記号を区切り文字としながら「年・月・日」の３つに分割できます。
ただし、そのstrsplit() が返す値はリストなので、unlist() によってベクトルに変換する必要があります。ベクトルに変換すれば、as.numeric() で数値化できます。
こうしたノウハウを利用しながら date.split() を記述します。
＊ date.split() の仕様
date.split() が返す値はデータフレームです。素材として与えられたデータフレームに、year, month, mday, quarter, wday という名前の列を付け加えて返します。５つの列を付加するわけです。
これらの追加された列は、それぞれ数値からなるvectorか、または、factorです。date.split() に与えるオプションによって、vector, factor のどちらにするかを選択できます。factorなら「2013年」とか「8月」などの文字列として扱えるし、一方、vectorなら 2013 とか 8 の数値として扱えます。
dd <- date.split(dtf, "長男の誕生日", "c")
とすると、dtf に year, month などの列を追加した dd を得られます。
第１引数のdtfは、素材となるデータフレームです。
第２引数の "長男の誕生日" は、日付データの列名です。列の番号で指定してもかまいません。指定した列が日付データでない時は、警告メッセージが出て、戻り値の dd は NULL になります。
第３引数が vector, factor の選択で、"n" か "numeric" だと vector、"c" か "character" だと factor になります。省略すると "c" が指定されたものとみなされます。つまり factor になります。
数値のvectorにおける「曜日」の値は 1〜7 で、1が日曜、2が月曜などです。
もし、素材となるデータフレーム dtf に year, month などの５つの列が既にあれば、戻り値である dd には year3, month3, mday3 などのように番号が付いた列名が追加されます。
この番号は、指定された日付データの列番号です。上に掲げた例でいうと、"長男の誕生日" が２番目の列に当たるので、year2, month2, mday2 などになります。
素材である dtf に year, month, mday などの５つの列がなければ、この番号は付きません。
以下、少し長くなりますが date.split() の定義を含むサンプルを掲げます。
「長男の誕生日」を四半期単位で数え上げます。第1四半期が何人、第2四半期が何人というふうに集計した結果を出力します。
−−−− tjs30.rb ここから
#! ruby -Ks
# 日付データのカテゴリー化 (coding: Windows-31J)
require "rrxwin"
##
rpro = <<EOS
date.split <- function(dtf, cc, xtype="c") {
cn <- colnames(dtf)
if (is.character(cc)) cc <- seq_along(cn)[cn %in% cc]
date.dat <- dtf[,cc]
if (class(date.dat) != "Date") {
warning("Date型を指定して下さい.")
return(NULL)
}
if (xtype == "character") xtype <- "c"
else if (xtype == "numeric") xtype <- "n"
if (xtype != "c" && xtype != "n") xtype <- "c"
cname <- c("year", "month", "mday", "quarter", "wday")
nn <- seq_along(cn)[cn %in% cname]
if (length(nn) > 0) # dtf中に既に year, month などの列がある
cname <- sapply(cname, function(i) sprintf("%s%d", i, cc))
fmt <- c("%s年", "%s月", "%s日", "第%s四半期")
wname <- c("日曜日","月曜日","火曜日","水曜日","木曜日","金曜日","土曜日")
ymd <- t(sapply(format(date.dat, "%Y-%m-%d"), function(ss)
as.numeric(unlist(strsplit(ss, "-")))))
qq <- as.numeric(cut(ymd[,2], breaks=c(1,4,7,10,13), labels=1:4, right=F))
ww <- weekdays(date.dat)
for (i in 1:length(wname)) ww[ww==wname[i]] <- i
dd <- data.frame(cbind(ymd, qq, ww))
colnames(dd) <- cname
if (xtype == "c") {
for (i in 1:(ncol(dd)-1)) { # 年・月・日・四半期のfactor化
dd[,i] <- as.factor(dd[,i])
lvl <- levels(dd[,i])
lbl <- sapply(lvl, function(j) sprintf(fmt[i],j))
dd[,i] <- factor(dd[,i], levels=lvl, labels=lbl)
}
dd[,ncol(dd)] <- factor(ww, levels=1:7, labels=wname) # 曜日のfactor化
}
for (i in cname) dtf[,i] <- NULL # 重複する列名を削除
return(cbind(dtf, dd))
}
dtf <- read.csv("data07.csv", header=T)
dtf$長男の誕生日 <- as.Date(dtf$長男の誕生日)
dtf$次男の誕生日 <- as.Date(dtf$次男の誕生日)
dd <- date.split(dtf, "長男の誕生日", "c")
table(dd$quarter)
EOS
##
hs = Rrx.rexec(rpro)
print hs[:sink]
−−−− tjs30.rb ここまで
date.split() の定義の最初のところで、日付データであるはずの date.dat が本当にDate型か否かをチェックしています。class(date.dat) の戻り値が "Date" であればDate型です。そうでない時は警告メッセージを出して NULL を返します。
sapply(format(……)) のところは分かりにくいですが、結果として、マトリクス(ymd)を返します。ymdの１列目には「年」、２列目には「月」、３列目には「日」がセットされます。いずれも数値データです。date.datが100個から構成されている場合、ymdの行数が100になります。
date.split() の定義に関するところは、数値のvectorをfactorに変換する処理などがあるため少しごちゃごちゃしていますが、説明は省略します。
念のため、スクリプトの出力結果を掲げておきます。長男の誕生日を四半期単位の分布でみたものです。
||第1四半期||第2四半期||第3四半期||第4四半期
||22||15||32||31
！！［補足］ 日付データへの cut(), seq() の適用について
前述した四半期単位でのカテゴリー化は、日付データの「月」を数値として取り出した上で cut() にかけました。つまり、日付データそのものをcut() にかけたわけではありません。
cut() は、Date型の日付データを直接取り扱うことができます。その場合、breaksオプションにはDate型からなるベクトルを指定します。
例えば、日付データ100個からなるxxに、2010〜2012年の３年間のデータが入っている場合、それを年単位でカテゴリー化するには次のようにします。
brk <- as.Date(c("2010-01-01", "2011-01-01", "2012-01-01", "2013-01-01"))
lbl <- c("2010年", "2011年", "2012年")
yy <- cut(xx, breaks=brk, labels=lbl, right=F)
上のようにすると、「2010年」〜「2012年」の３つの値から構成されるyyがfactorとして生成されます。yyの要素は100個あります。
yyはカテゴリーデータなので、それを基にしてクロス集計等を行うことができます。
それから、2013年の365日分の日付データを生成するのに seq() を使うことができます。
xx <- seq(as.Date("2013-01-01"), as.Date("2013-12-31"), by="days")
こうすると、xxに365日分の日付データがセットされます。
"days" のところには "months", "weeks" を書くこともできます。
"months" にすれば、毎月１日の日付データが12ヶ月分生成されます。"weeks" だと、2013-01-01, 2013-01-08, 2013-01-15, …… のように展開する日付データが生成されます。
xx <- seq(as.Date("2009-01-01"), as.Date("2013-01-01"), by="years")
とすれば、2009-01-01, 2010-01-01, …… 2013-01-01 の５つの日付データが生成されます。
日付データの説明は、この辺でおわりにしようと思います。
「xx <- Sys.Date()」とすれば、xxに今日の日付(Date型)がセットされることなど、他にもいろいろな話題があります。
また、日付と合わせて時刻も記録・保持するdate-time型に言及してませんが、これは、世界標準時の1970年１月１日 0時0分0秒からの経過秒数を実態とするクラスです。日本標準時は世界標準時より９時間早いので、-32400秒のずれがあります。
「xx <- as.POSIXct("1970-01-01 00:00:00")」とすれば、普通は(環境変数TZがちゃんと設定されていれば)日本標準時の秒数がxxにセットされます。秒数といっても独自クラスなので print(xx) とすれば、「1970-01-01 JST」と表示されます。as.numeric(xx) の戻り値は -32400 です。
一方、「xx <- as.POSIXlt("1970-01-01 00:00:00")」とすれば、xxがlist型になり、xx$year, xx$mon などとして、年・月・日・曜日・時・分・秒などを数値として簡単に取り出すことができます。ただ、1でなく0をスタートの値とするものがあるので注意して下さい。
difftime() は、その名前から分かるように、本来、date-time型の差分を算出するものです。
他に時系列を扱うためのクラス(ts)などもあり、日付や時刻にまつわる話題だけでも多岐にわたりますが、私には手に余るので省略します。
!4. その他
数値データの扱いについて、これまで触れなかったことで気になる点を簡単に書いておきます。
!!(1) 欠損値の扱い
○ 要点１：na.rmオプションでNAを無視
○ 要点２：subset() によって、NAを取り除いたデータフレームを生成
数値データを扱う場合、基本的には欠損値を予め取り除いておくのが無難です。
最大値を求める max(), 最小値を求める min(), 平均値を求める mean() などでは、na.rm というオプションがあります。これは、与えられた引数の中にNAがあるとき、それを無視して計算するか否かを指定するものです。
mean(xx, na.rm=TRUE)
とすれば、NAを無視します。
デフォルトは「na.rm = FALSE」となっているので、NAがあると計算できず、mean() はNAを返します。意味のある答えを得られません。
欠損値を扱う場合、na.rmのようなオプションを利用するのも一つの方法ですが、やはり予め欠損値を取り除いておくのがいいと思います。
仮に data07.csv の「長男の誕生日」に欠損値がある場合、次のようにしてNAを取り除くことができます。
dtf <- read.csv("data07.csv", header=T, na.strings="")
dtf2 <- subset(dtf, !is.na(dtf[["長男の誕生日"]]))
こうすると、「長男の誕生日」の中のNAである行が取り除かれて、その結果がdtf2にセットされます。「長男の誕生日」の列だけが行数を減らすのではなく、すべての列の行数が減ります。
「nrow(dtf) - nrow(dtf2)」とすれば、取り除かれた行が何行あるか、その数を確認できます。
「長男の誕生日」と「次男の誕生日」の両方からNAを取り除いた結果を得たい時は、次のようにします。
dtf3 <- subset(dtf, !is.na(dtf[["長男の誕生日"]]) &
!is.na(dtf[["次男の誕生日"]]))
こうして得られたdtf3には、NAが含まれていないことになります。
数値データの列がもっと沢山あって、そのすべてのNAを削除したい場合は、例えば次のようにします。
nums <- c("長男の誕生日", "次男の誕生日", ……)
xx <- dtf
for (i in nums) xx <- subset(xx, !is.na(xx[[i]]))
上のようにして得られたxxは、指定の列(numsで指定した列)のいずれにもNAがない状態になっています。
列を指定するためのnumsは、「nums <- 2:5」のように、数字で指定しても大丈夫です。
!!(2) 文字列と数値が混在している場合の処理
○ 要点１：as.numeric() によって、数値化できるもののみを数値に変換
○ 要点２：as.Date() によって日付データに変換する時は、ダミーの第１要素を挿入
○ 要点３：grep() によって、パターンにマッチする要素のみを抽出
変数xxに、文字列と数値が混在して入っているとします。例えば次のようなものです。
xx <- c("3", "abc", "1.55", "xyz", "-3.14e-001")
このようなxxに対して as.numeric(xx) を適用して数値化を試みると、次のような結果になります。
3.000 NA 1.550 NA -0.314
つまり、数値化できない "abc" とか "xyz" は、NAになります。
このことを利用すれば、文字列と数値が混在している中から数値を取り出すことができます。
ここで、csvを読み込む場合について考えてみます。
ある列が、上のxxのように混在型だとします。具体的には次のようなcsvです。
ID,data
1,3
2,abc
3,1.55
4,xyz
5,-3.14e-001
このcsvをread.csv() で読み込んで、データフレームdtfにセットしたとします。
このとき、dtf$data は、内容的に前述のxxと同じになるはずです。
ただし、dtf$data は、xxのように文字列から構成されるvectorではなく、factorです。なので、as.numeric(dtf$data) としても、意図したような結果が得られません。
dtf$data は５つの要素から構成されていますが、各要素には整理番号が割り振られています。同じ値の要素には同じ整理番号が付けられますが、今回の例では同じものがないので1〜5のどれかが割り振られます。
as.numeric(dtf$data) とすると、その整理番号がvectorになって返されます。具体的には c(3,4,2,5,1) が返されます。
xxを数値化した時と同じ結果を得るためには、dtf$data を一度 as.character() で文字化して、その上で数値化します。具体的には次のようにします。
nn <- as.numeric(as.character(dtf$data))
このようにすると、nn が c(3, NA, 1.55, NA, -0.314) となります。
念のため、先のcsvを読み込んで数値化するRプログラムを掲げておきます。
−−−− Rプログラムここから
dtf <- read.csv("test.csv", header=T)
nn <- as.numeric(as.character(dtf$data))
print(nn)
−−−− Rプログラムここまで
次に、日付データに話を移します。
xxが文字列から構成されるベクトルの場合、その第１要素が日付データに変換可能なものであれば、as.Date() によってベクトル全体を日付型に変換できます。日付データに変換できる要素は変換され、できない要素はNAに変換されます。
ただし、そのような変換が行われるのは、第１要素が "2013-07-13" のように日付データに変換可能な場合に限られます。"abc" とか "123" などのように変換できない文字列が第１要素だと、エラーが発生します。
「それなら、第１要素として変換可能な文字列をあえて挿入してやり、変換後に第１要素を削除すればいい。」という逆手の発想がうかびます。この方法を採ると次のようなプログラムになります。
xx <- c("123", "2012-6-8", "abc", "1998-10-17")
dd <- as.Date(c("1970-1-1", xx))[-1]
print(dd)
「yy <- xx[-1]」とすると、xxの第１要素が削除されたものがyyにセットされますが、上は、その記法を用いたものです。
上のやり方だと、日付データに変換できない要素がNAになるので、xxとddの要素数が同じになります。
そうではなく、日付データに変換可能なものだけを取り出してDate型にしたい場合は、正規表現とgrep() を用いて次のようにします。
xx <- c("123", "2012-6-8", "abc", "1998-10-17")
ptn <- "[0-9]{4}-[0-9]{1,2]-[0-9]{1,2}"
dd <- as.Date(xx[grep(ptn,xx)])
print(dd)
上のやり方だと、日付データに変換可能なものだけがddにセットされるので、要素数が２になります。
grep(ptn,xx) は、正規表現ptnにマッチするものの要素番号を返します。具体的には c(2,4) を返します。つまり、日付データの書き方にマッチするものが２番目と４番目ということです。
「yy <- xx[c(2,4)]」とすると、xxの２番目の要素と４番目の要素の２つがyyにセットされます。このyyをDate型に変換すれば、目的が達成できます。
なお、このやり方だと、xxの中に日付データに変換可能な要素が１つもない場合、ddには character(0) がセットされます。length(dd) が 0 となります。
--------
以上で「単純集計に関する覚え書き／数値データの集計」をおわりにします。
数値データの処理については、まだまだ様々な事柄がありますが、ごく基本的な部分には言及できたのではないかと思います。
Copyright (C) T. Yoshiizumi, 2013 All rights reserved.