毎回やる儀式

補遺

前回時間が足りなくてできなかったところがあります。探索的データプロットというところです。今回はそこから始めたいと思います。

今回改めてデータファイルを用意しましたので，baseball2015a.csv というデータファイルをプロジェクトフォルダに入れ，次の一行で読み込んでください。

b.data <- read.csv("baseball2015a.csv",head=T)

探索的データプロット

データは図にする。プロットする。 これが分析の基本です。まずデータがどういう特徴を持っているかを，しっかりと目で把握しましょう。 それができていないと，データにモデルを当てはめるという話が始まりませんので。

ggplot2の導入

探索的データプロットのパッケージ，ggplot2をインストールします。

install.packages("ggplot2")

インストールしただけでは使えません。library関数（またはrequire関数）で装備します。

library(ggplot2)

まずはデータと全体設定を持ったggplotオブジェクトを作ります。

gp <- ggplot(b.data,aes(x=year,y=pay,color=team))

点を打ってみましょう。

gp_point <- gp + geom_point()
plot(gp_point)

色だけじゃなく，大きさや形状を変えてプロットすることもできます。

gp <- ggplot(b.data,aes(x=year,y=pay,color=team,shape=position,size=blood))
gp_point <- gp + geom_point()
plot(gp_point)

タイトルや軸ラベルも書いておきましょう。

gp <- ggplot(b.data,aes(x=year,y=pay,color=team,shape=position,size=blood))
gp_point <- gp + geom_point()+xlab("プレイ期間")+ylab("年俸")+ggtitle("野球選手の年棒とプレイ期間の関係")
plot(gp_point)

画を複数に分けることもできます。

gp <- ggplot(b.data,aes(x=year,y=pay,color=team,shape=position))
gp_point <- gp + geom_point() + facet_wrap(~blood)
plot(gp_point)

横軸がカテゴリカルであれば，箱ひげ図や棒グラフになります。

gp <- ggplot(b.data,aes(x=team,y=pay,fill=position))
gp_box <- gp + geom_boxplot()
plot(gp_box)

gp_bar <- gp + geom_bar(stat="identity",position="dodge")
plot(gp_bar)

# 平均が見たい
gp_bar <- gp + stat_summary(fun.y=mean,geom="bar",position="dodge")
plot(gp_bar)

# 中央値が見たい
gp_bar <- gp + stat_summary(fun.y=median,geom="bar",position="dodge")
plot(gp_bar)

エラーバーを追加する例。

gp <- ggplot(b.data,aes(x=position,y=pay,fill=position))
gp_bar <- gp + stat_summary(fun.y=mean,geom="bar",position="dodge") + 
                stat_summary(fun.data=mean_se,geom="errorbar",width=0.5)
plot(gp_bar)

様々な角度からデータをプロットし，分布の形状や変数間関係をしっかり見極める癖をつけましょう。

本題

今週のテーマは クラスタリング です。

多変量解析の基本

多変量解析とは，要約すること，まとめあげること，でもあるわけです。 まとめたものを同じ「クラスター」に入れる。これを機械的にするのがクラスタリングです。

どうやってまとめるか

基本は似ているもの同士をまとめていきます。

似ているって何

数字と数字が似ている＝距離が近い，です。距離は引き算で出します。データで見ていきましょう。先ほどと同じデータセットを例にします。

例えば「身長」変数を例にします。

b.data[1,] #case1
##         name team   pay position year age height weight blood home BirthY
## 1 杉内　俊哉 巨人 50000     投手   14  34    175     82   A型 福岡   1980
##   BirthM BirthD throw hit
## 1     10     30    左  左
b.data[2,] #case2
##         name       team   pay position year age height weight blood home
## 2 金子　千尋 オリックス 50000     投手   11  31    180     77   O型 新潟
##   BirthY BirthM BirthD throw hit
## 2   1983     11      8    右  左
b.data[3,] #case3
##         name team   pay position year age height weight blood home BirthY
## 3 内海　哲也 巨人 40000     投手   12  33    186     93   A型 京都   1982
##   BirthM BirthD throw hit
## 3      4     29    左  左
b.data$height[1] #case1の身長
## [1] 175
b.data$height[2] #case2の身長
## [1] 180
b.data$height[3] #case3の身長
## [1] 186
b.data$height[1]-b.data$height[2]　#case1とcase2の差
## [1] -5
b.data$height[1]-b.data$height[3]　#case1とcase2の差
## [1] -11
b.data$height[2]-b.data$height[3]　#case1とcase2の差
## [1] -6
sqrt((b.data$height[1]-b.data$height[2])^2)　#case1とcase2の距離
## [1] 5
sqrt((b.data$height[1]-b.data$height[3])^2)　#case1とcase2の距離
## [1] 11
sqrt((b.data$height[2]-b.data$height[3])^2)　#case1とcase2の距離
## [1] 6

これをみると，case1-2の距離はcase1-3より近い。 こういうデータを全部の組み合わせについて見て行きたいわけです。 とても大変な計算量ですが，統計ソフトだと簡単。

こうした組み合わせをひとつずつ計算するのは面倒だけど，dist関数を使うと簡単。 ** head()はデータの上の部分だけを取り出す関数です。全部やるとたくさん出力されるので，便宜的にそうしています。実際はデータ全部を使うのが普通です。**

dist(head(b.data$height)) #身長変数について，最初の6行分だけ，距離を計算する
##    1  2  3  4  5
## 2  5            
## 3 11  6         
## 4  8  3  3      
## 5 10  5  1  2   
## 6  6  1  5  2  4

こうして,近い方から「これは同じクラスタだね」とまとめていく。例えばcase5とcase3,case6とcase3は距離が1しかないから同じクラスタにする。簡単ですね。

変数が二つ以上あると

身長だけでなく，体重も考えたい，というときは，ユークリッド距離を算出します。

dist(head(b.data[,7:8]))
##           1         2         3         4         5
## 2  7.071068                                        
## 3 15.556349 17.088007                              
## 4 13.601471 16.278821  3.000000                    
## 5 14.866069 16.763055  1.000000  2.000000          
## 6 16.155494 20.024984  6.403124  4.472136  5.656854

変数が三つ以上になっても基本的に同じです。

もっともこの場合は，身長(cm)と体重(kg)という異なる単位のものを等価に扱っているという問題があります。単位の違いを考慮するならば，データを標準化*しておいたほうがいいですね。標準化する関数はscaleというのがあります。

scale(head(b.data[,7:8]))
##       height     weight
## 1 -1.6736548 -0.9178485
## 2 -0.4184137 -1.5582078
## 3  1.0878756  0.4909422
## 4  0.3347310  0.4909422
## 5  0.8368274  0.4909422
## 6 -0.1673655  1.0032297
## attr(,"scaled:center")
##    height    weight 
## 181.66667  89.16667 
## attr(,"scaled:scale")
##   height   weight 
## 3.983298 7.808115

これを使って距離データを準備しておくことにしましょう。

b.distance <- dist(scale(b.data[,7:8]))

距離を順にまとめていく階層的クラスター分析法

まとまったclusterをAとします。次は，cluster Aとcase2,3,4,6はどこが近いかな？あるいは他のクラスターB(case4と3とか)との距離は？ということを考える必要があります。 まとめ方によっていろいろ変わりますが，代表的なものを三つ挙げておきます。

  • 最も近いもの同士をまとめる 最近隣法
  • 最も遠いもの同士をまとめる 最遠隣法
  • 二つのクラスターを統合したときに差がハッキリ別れるようにする ward法

順次まとめて束ね上げていく

顎の長さデータをつかって，階層的クラスター分析をやってみよう。

result <- hclust(b.distance,method="ward.D2")　#ward法
plot(result)

このような図をデンドログラムという。これを見ながら，いくつのグループに分けるかを考えます。いくつのグループに分けるかはご随意に！

カットしてから考える

2つのグループかしら？というあたりをつけたとしましょう。この木の枝をカットし，元変数に足しておきます。

b.data$class <- cutree(result,2)
b.data$class
##   [1] 1 1 2 2 2 2 2 2 1 2 2 2 1 2 2 2 2 1 1 1 2 2 2 2 2 2 2 2 1 2 1 2 2 1 1
##  [36] 1 2 1 1 2 2 1 2 2 2 1 2 1 1 2 2 2 2 2 2 1 2 1 2 1 1 1 2 2 1 2 1 2 1 2
##  [71] 1 1 1 2 2 2 1 1 2 2 1 1 2 1 2 2 1 2 1 2 2 1 1 2 1 2 1 2 1 1 2 2 2 1 2
## [106] 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 1 1 2 1 2 2 2 2
## [141] 1 2 1 1 1 2 2 2 1 2 2 1 2 1 2 1 2 1 1 2 2 2 1 2 2 1 2 2 1 1 1 1 1 1 1
## [176] 1 2 2 1 1 1 2 1 1 2 2 2 2 2 1 1 1 1 1 1 1 2 2 1 1

さてこの違いはどこから来てるのでしょうか？ 度数を数えるときはクロス集計表というのがいいですね。あるいは検定して調べるとか。

xtabs(~position+class,data=b.data)
##         class
## position  1  2
##   外野手 12 22
##   投手   42 59
##   内野手 24 33
##   捕手    4  4
chisq.test(xtabs(~position+class,data=b.data))
## Warning in chisq.test(xtabs(~position + class, data = b.data)): Chi-squared
## approximation may be incorrect
## 
##  Pearson's Chi-squared test
## 
## data:  xtabs(~position + class, data = b.data)
## X-squared = 0.76852, df = 3, p-value = 0.857
t.test(pay~class,data=b.data)
## 
##  Welch Two Sample t-test
## 
## data:  pay by class
## t = -3.7347, df = 193.13, p-value = 0.0002472
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -7781.285 -2402.961
## sample estimates:
## mean in group 1 mean in group 2 
##        10456.10        15548.22

とりあえず多くのデータがあって困っているときに，いくつかのグループに分けて，グループごとの平均や相関などをみて，このクラスターはこういう特徴だったのだと事後的に解釈することは少なくありません。そのために，検定があったり，群ごとに検証するモデリングがでてきたりします。

本日の課題

  • 年俸と守備位置とBMIの関係をグラフにして表現してください。何グラフを用いても良い。
  • 投手データかバッターのデータ，いずれか一方を使ってクラスター分析をし，選手をいくつかのクラスターに分類してください。
  • 分類したクラスターはどういう基準で分かれているのでしょうか。図表や検定を使って考察してください。