2015年3月17日火曜日

Rのtapply関数を使って、複数の度数分布表・ヒストグラムをまとめて出力する

tapplyを使えば、データフレームのカテゴリの列でグループ分けしつつ、値の列に対して関数を適用するということができます。

・・・と言っても分かりにくいと思いますので、具体的なサンプルをあげてみますと、

> # サンプルデータを作る
> 名前 <- c("A", "B",  "C", "D",  "E", "F")
> 性別 <- c("男", "男", "女", "男", "女", "女" )
> 身長 <- c(175,  165, 165,  170, 160, 155  )
> df <- data.frame(名前, 性別, 身長)
> df

  名前 性別 身長
1    A   男  175
2    B   男  165
3    C   女  165
4    D   男  170
5    E   女  160
6    F   女  155

>
> # ここから本題
> tapply(df$身長, df$性別, mean)
 
 女  男
160 170


性別でグループ分けして、身長の平均を算出する(mean関数を適用)という感じです。

第3引数にhist関数を指定してやれば、グループ分けを適用したあとに、複数の度数分布表を算出したり、複数の度数分布図(ヒストグラム)を書いたりすることもできます。

もう少し大きいサンプルデータを使いましょう。↓この本に載っていたサンプルを使わせていただきます。

データマイニング入門

↓こちらからダウンロードできる、"年収.csv"を使います。

データマイニング入門-Rで学ぶ最新データ解析 - 東京図書


↓こんな感じのサンプルデータです。

> d <- read.csv("年収.csv")
> head(d)
  学位 地位 性別 現職年数 学位年数  年収
1 博士 教授 男性       25       35 36350
2 博士 教授 男性       13       22 35350
3 博士 教授 男性       10       23 28200
4 博士 教授 女性        7       27 26775
5 修士 教授 男性       19       30 33696
6 博士 教授 男性       16       21 28516

tapplyでhistを適用しても、バババっとヒストグラムが次々描画される(そして消えていく)ことになってしまうので、ちょっと工夫が必要です。

工夫と言っても、戻り値を受け取って、所望の形式で出力させるだけですが。↓こんな感じ。

> # 度数の集計のため(描画が目的ではない)
> h <- tapply(d$年収, d$性別, hist, breaks=seq(15000, 40000, 5000))
>
> # 度数分布表の出力
> for ( i in 1:length(h) ) {
+   print ( names(h[i]) )
+   print ( data.frame(階級値=h[[i]]$mids, 度数=h[[i]]$counts) )
+   cat ( "\n") # 改行
+ }
[1] "女性"
  階級値 度数
1  17500    6
2  22500    5
3  27500    2
4  32500    0
5  37500    1

[1] "男性"
  階級値 度数
1  17500    9
2  22500   13
3  27500    9
4  32500    5
5  37500    2

それぞれに適用されたhistの戻り値がリストの形で h に格納されてますので、h[[1]], h[[2]], h[[3]], ... という感じで、それぞれにアクセスします。$midsには$breaksの中間の値が入っているので階級値として使えます。$countsには度数が入っています。

違うグループ分け(地位)でやってみましょう↓

> # 今度は、地位でグループ分け
> h <- tapply(d$年収, d$地位, hist, breaks=seq(15000, 40000, 5000))
>
> # 度数分布表の出力
> for ( i in 1:length(h) ) {
+   print ( names(h[i]) )
+   print ( data.frame(階級値=h[[i]]$mids, 度数=h[[i]]$counts) )
+   cat ( "\n") # 改行
+ }
[1] "教授"
  階級値 度数
1  17500    0
2  22500    2
3  27500   10
4  32500    5
5  37500    3

[1] "講師"
  階級値 度数
1  17500   15
2  22500    3
3  27500    0
4  32500    0
5  37500    0

[1] "助教"
  階級値 度数
1  17500    0
2  22500   13
3  27500    1
4  32500    0
5  37500    0

次は、ヒストグラムを書いてみましょう。前述のスクリプトでも描画自体はされていたのですが、複数なので上書きされてしまうのと、タイトルが入っていないのでどのグループのヒストグラムか分からない、というところをちょっと工夫します↓

# あらかじめグループ名を取得するために
h <- tapply(d$年収, d$性別, hist, breaks=seq(15000, 40000, 5000))
 
# カテゴリの数だけヒストグラムを縦に並べるための設定
par(mfrow=c(length(h),1))
 
# ヒストグラムの描画
for ( i in 1:length(h) ) {
  plot(h[[i]], main=names(h[i]), xlab="年収", col="gray")
}

tapplyでヒストグラムを描く(性別でグループ分け)

別のグループ分けを適用したいときは、tapplyの引数を1つ変えるだけ↓

# 今度は、地位でグループ分け
h <- tapply(d$年収, d$地位, hist, breaks=seq(15000, 40000, 5000))
 
# カテゴリの数だけヒストグラムを縦に並べるための設定
par(mfrow=c(length(h),1))
 
# ヒストグラムの描画
for ( i in 1:length(h) ) {
  plot(h[[i]], main=names(h[i]), xlab="年収", col="gray")
}

tapplyでヒストグラムを描く(地位でグループ分け)

カテゴリの種類が多い場合や、カテゴリが含む水準が多い場合は、いちいち手でやっていられないですからね。ちょっとExcelにはできない芸当だと思います。




0 件のコメント:

コメントを投稿