ラベル 箱ヒゲ図 の投稿を表示しています。 すべての投稿を表示
ラベル 箱ヒゲ図 の投稿を表示しています。 すべての投稿を表示

2012/09/02

文系のための「数の可視化」(8)

分析の基本は、分析の前にデータの性質を理解することである。
分析対象の性質を理解した上で、適切な分析手法を選択する必要がある。

残念ながら有名な分析手法の名前に踊らされる人は多く、
「○○の研究での実績があるから...」とか、「○○先生が使っていたし...」とか、
不可解かつ非論理的な理由で分析する人まで居る。理解できない。

信じ難いのであるが、そんな話が常識となっている分野も存在する。
まぁ、このブログを読んでいる人には、そのような人は居まい。

さて、そのような愚痴はさておき、
本日は多次元データの要約を可視化する方法について考える。

ケトレーの話にも通じるが、データを観察する上で最も基本となるのは、
データの中心からのバラツキであって、
分散標準偏差共分散相関係数といった指標があった。
また、これらの指標に加えて、四分位ヒンジというのもあった。

さらに、それを可視化する方法として、
ヒストグラム散布図箱ヒゲ図、などがあった。

ところで、分散共分散あるいは相関係数を同時に表現する方法として、
分散共分散行列相関係数行列があった。

そこで、今回は分散共分散行列と相関係数行列を可視化する方法について述べる。

分散共分散行列相関係数行列はどのような構造であったか?

確か、分散共分散行列の場合は、対角成分に分散が入っていて、
その他の成分には共分散が入っているような対象行列であった。

一方、相関係数行列は、対角成分が「1」であり、
その他の成分には、相関係数が入っているような対象行列であった。

つまり、結論から述べると散布図行列というのは、
対角成分にバラツキを表すグラフが入り、
それ以外の成分には、散布図が入るような視覚化方法である。

散布図行列には、様々なバリエーションが存在し得るので、
今回は、基本的な見方基本的な形式のみを紹介する。

では、さっそく、データを読み込んでみる。
今回は、いつもとは異なった方法でデータの読み込みを行なってみる

今回のデータは少し大きいので、Google Spreadsheet にデータを置いてある。
このデータを利用するためには、RCurlというパッケージを必要とする。

インストールした記憶の無い人は、まずは、このパッケージをインストールする。
インストールの説明は、新しい箱ヒゲ図の話でしているので、そちらを参照。

install.packages("RCurl", dep=TRUE)

インストールできたら、とりあえず、ライブラリを読み込む

library(RCurl)

ライブラリを無事に読み込めたら、次は、以下のコマンドを実行する。

# Google Spreadsheet からデータをダウンロードする
data <- getURL("https://docs.google.com/spreadsheet/pub?key=0AtOtIs5BjRVhdHo0c0pjb29OTU9aV3BmQUFJUWJQa0E&single=true&gid=0&range=A1%3AK48&output=csv")

# ダウンロードしたデータを読み込む
X <- read.table(textConnection(data), header=TRUE, row.names=1, sep=",")

エラーが出なければ問題は無いはず。念のために、データの確認を行う。

head(X)

head()関数は、データの先頭6行を読み込んで表示するための関数。
データの中身を確認するために用いる。実行結果は以下の通り。

> head(X)
               HM  ELD   HC  CHR   SC   LI   SP  ENV  DIS   IC
Hokkaido     11.3 35.1 34.1 21.2 44.3 14.4 15.5 24.1  5.4 17.9
Aomori-ken   15.2 21.0 22.4 13.3 55.2  7.9  8.4 19.9  8.3  8.0
Iwate-ken     6.0 26.1 12.5 15.9 33.1 10.6 15.6 16.9 11.8 26.9
Miyagi-ken   22.8 40.6 32.8 20.7 43.4 13.4 14.3 22.6 15.2 40.4
Akita-ken     4.9 19.5 25.1 16.1 41.2  7.4  9.8 12.8 10.4 33.2
Yamagata-ken  7.4 26.5 13.4 15.8 33.3  9.5 11.4 13.0  5.9  8.7

このデータは、総務省統計局からダウンロードしたデータを加工したもので、
都道府県別のボランティア活動に従事した人の年平均従事時間(日/年)のデータ。

書く変数の意味は次の通り。
  • HM:健康や医療サービスに関係した活動
  • ELD:高齢者を対象とした活動
  • HC:障害者を対象とした活動
  • CHR:子供を対象とした活動
  • SC:スポーツ・文化・芸術・学術に関係した活動
  • LI:まちづくりのための活動
  • SP:安全な生活のための活動
  • ENV:自然や環境を守るための活動
  • DIS:災害に関係した活動
  • IC:国際強力に関係した活動
当然、全ての人に聞くことはできないので、サンプル調査の結果であり、
男女を合わせた総数での結果となっている。

さて、まずは、分散共分散行列相関係数行列を確認。
小数点が大きすぎるので、小数点第三位までで丸め込む。

# 分散共分散行列
round(cov(X), 3)

# 相関係数行列
round(cor(X), 3)

実行結果は以下の通り。

> # 分散共分散行列
> round(cov(X), 3)
         HM    ELD      HC    CHR     SC     LI     SP    ENV    DIS      IC
HM   27.370  0.750   5.470 -2.063  7.600  2.305 -5.029  5.586 -1.925 -12.463
ELD   0.750 38.669  22.176 10.859 -2.048  9.253  4.594 16.180 -0.199  44.576
HC    5.470 22.176 103.638 11.745  6.029  9.985  3.673 12.545 -0.205  64.775
CHR  -2.063 10.859  11.745 13.820 -0.021  7.465  8.619 12.185  0.559  11.328
SC    7.600 -2.048   6.029 -0.021 52.300 -2.108  0.040  1.468  1.434   7.519
LI    2.305  9.253   9.985  7.465 -2.108 12.677  5.382 20.473  0.950   0.394
SP   -5.029  4.594   3.673  8.619  0.040  5.382 21.932 12.175  1.550  13.332
ENV   5.586 16.180  12.545 12.185  1.468 20.473 12.175 58.090  0.927  -2.488
DIS  -1.925 -0.199  -0.205  0.559  1.434  0.950  1.550  0.927  7.540  12.452
IC  -12.463 44.576  64.775 11.328  7.519  0.394 13.332 -2.488 12.452 279.265
> # 相関係数行列
> round(cor(X), 3)
        HM    ELD     HC    CHR     SC     LI     SP    ENV    DIS     IC
HM   1.000  0.023  0.103 -0.106  0.201  0.124 -0.205  0.140 -0.134 -0.143
ELD  0.023  1.000  0.350  0.470 -0.046  0.418  0.158  0.341 -0.012  0.429
HC   0.103  0.350  1.000  0.310  0.082  0.275  0.077  0.162 -0.007  0.381
CHR -0.106  0.470  0.310  1.000 -0.001  0.564  0.495  0.430  0.055  0.182
SC   0.201 -0.046  0.082 -0.001  1.000 -0.082  0.001  0.027  0.072  0.062
LI   0.124  0.418  0.275  0.564 -0.082  1.000  0.323  0.754  0.097  0.007
SP  -0.205  0.158  0.077  0.495  0.001  0.323  1.000  0.341  0.121  0.170
ENV  0.140  0.341  0.162  0.430  0.027  0.754  0.341  1.000  0.044 -0.020
DIS -0.134 -0.012 -0.007  0.055  0.072  0.097  0.121  0.044  1.000  0.271
IC  -0.143  0.429  0.381  0.182  0.062  0.007  0.170 -0.020  0.271  1.000

さて、この値のどの部分を注意して観察するべきであったか?
確か、相関係数行列の「1」に近いものと「-1」に近いものに注意するべきであった。
この中では、「0.754」が最も高い。マイナス方向では、最も低いのが「-0.205」である。

では、正負の相関の有無の断定となると色々と手続きが必要であるが、
自然や環境を守るための活動(ENV)」と「まちづくりのための活動(LI)」は、
関係しているように思える。逆に、強い負の関係というのは無さそう。

ふむ。この状況、一つ一つを丁寧に見れば、解らなくは無いのだが、どうも見難い。
そこで、この状況を散布図行列を使って可視化してみる。
pairs(X)

pairs()関数は、散布図行列を作成するための関数であり、
分散共分散行列あるいは相関係数行列の状況を可視化したようなもの。

したがって、相関係数行列の各成分に相当する散布図が、
散布図行列の各要素に相当する散布図になっている。

また、散布図行列の対角を挟んで上側下側散布図は、
x軸とy軸が入れ替わっているだけで、結局は同じ状況を示している。
ここでは、とりあえず、下側にだけ注目してみることにする

まずは、先程の比較的関係が強そうなENVとLIの関係に着目してみる。

ENVと他の変数との相関係数がどこに入っていたかと言うと、
相関係数行列では8行目に入っていた。LIとの関係は、6列目

したがって、散布図行列においても8行目を見てみる。
すると、y軸の目盛りが右側に書いてある。これが、ENVの値

LIとの関係は、相関係数行列における6列目であったので、
そのまま、右方向に進み、6番目の散布図を見てみる。これが、ENVとLIの散布図。
LIの値はというと、6列目に並んでいる散布図の一番上の散布図のx軸に書いてある。

y軸の目盛りは左右交互に、x軸の目盛りは上下交互に書いてあるが、
この見方は、相関係数行列と対比しながら見れば解るであろう。

さて、この8行目の6列目の散布図を見てみると、
他の散布図と比較すると、対角線上に綺麗に点が並んでいるように見える。

LIの値が右側に増えるに従って、ENVの値は上に向かって増加しているので、
この二つの変数は連動している

散布図行列だけを見てみると、LI:CHRSP:CHRも、なんとなく右肩上がり。
この散布図の場所は、それぞれ、6行目4列目7行目4列目なので、
そこに該当する相関係数を見てみると、それぞれ、「0.564」と「0.495」となっていて、
今回のデータで見てみると、相対的に相関係数は高い値となっている。

では、逆に対角線に集まって「なさそうな」ものはどれか?
対角を挟んで下側だけで考えてみる。

SC:CHR(5行目4列目)と、SP:SC(7行目5列目)が、対角線上には並んでいなさそう。
相関係数は、それぞれ、「-0.001」と「0.001」であり、やはり、関連は無さそうである。

この二つは、特に関係してい「なさそうな」組み合わせであるが、
散布図が全体的に散らばっているならば、
散布図に対応する相関係数は±1」よりも「0」に近いハズ。

地域に関する活動は、環境や自然に関わる活動や子供に関する活動と関係していて、
さらに、子供に関する活動と安全な生活に関わる活動が関係している。
直感的に、地域住民の意識として、住環境と子供の安全が重要であるのは納得できる。

逆に、スポーツ・文化・芸術・学術に関係した活動は、健康に関わりそうな話であって、
子供に関する活動や安全な生活に関わる活動とは関係しない、というのも頷ける。

なるほど。直感的に納得できそうな結果である。

さて、散布図行列相関係数行列を可視化したようなものであり、
データの全体を俯瞰する上で、非常に強力なツールではあるが、注意すべき点がある

相関係数行列の場合は、変数が多いと値を見落とす危険性があり、
散布図行列の場合は、目の錯覚で、間違った解釈をする場合がある。

この二つの結果を相互に見比べることが重要!!

今度は、散布図行列をもう少し見やすく整えてみる。
散布図行列には、様々な描き方が存在しているが、
ここでは、psych というパッケージを用いた方法を紹介する。

とりあえず、パッケージのインストールし、ライブラリを読み込む。
既に、インストール済みの場合は、ライブラリの読み込みのみ。

install.packages("psych", dep=TRUE)
library(psych)

psych散布図行列のメソッドは、pairs.panels()関数を用いる。
pairs.panels(X)

psych散布図行列を用いると、対角線にはヒストグラムが配置され、
下側には散布図、上側には相関係数が配置される。

それにしても、色々と煩雑なものが入っていて見難い。
それぞれに意味はあるのだが、必要なもの「のみ」を表示することも重要。
pairs.panels(X, smooth=FALSE, density=FALSE, ellipses=FALSE, scale=TRUE)

パラメータに関しては以下の通り。

  • smooth: 平滑線の描画の有無(今回は表示しない)
  • density: ヒストグラムにカーネルを重ねるかの有無(今回は重ねない)
  • ellipses: 散布図に相関を円で表したものを表示するか否か(今回は表示しない)
  • scale: 上側の相関係数の表示の大きさを相関の強さで変えるか(今回は変える)
このように表示すると、分散共分散行列相関係数行列を見事に可視化できる。
分散に相当するヒストグラムが対角線上に配され、上側には相関係数行列もある。

最初の方法では、相関係数行列散布図行列を見比べる必要があったが、
この方法を用いることで、必要な情報を一目で確認することができる。

今回は、対角線上にヒストグラムを配したが、対角成分に箱ヒゲ図を配する方法もある。
本ブログでは、使用を避けているが、Rコマンダーと呼ばれる、
Rのグラフィカルインタフェースを導入すると、散布図行列の対角は箱ヒゲ図になる。

また、Rの初期状態で利用可能なpairs()関数を拡張することもでき、
その方法は、pairs()関数のヘルプの例に記載されている。
余力のある人は、試してみると良い。それほど、難しくはない。

ところで、相関係数行列を可視化する方法は、散布図行列だけではない。
参考までに、相関プロットも紹介する。


cor.plot(cor(X))



cor.plot()関数を用いると、相関係数行列色分けで表示することができる。
この図において、白色のセルは相関係数が低く、青色のセルは値が高い。
逆に、赤色のセルは赤色で表示されている。

変数の数が多くなると、どうしても、相関係数行列散布図行列が見難くなる。
そのような場合には、この方法で可視化して全体を把握することもできる。

2012/08/31

文系のための「数の可視化」(7)

ヒストグラム箱ヒゲ図は、両方とも分散あるいは標準偏差の状況を
視覚的に解りやすく表現したものである。

ヒストグラムは、全体的な分布の形を視覚化し、
データの偏り方を視覚化する

一方、箱ヒゲ図は、全体的な分布の範囲を四分位で表し
データの偏り方を視覚化する

分散標準偏差というのは、全体のバラツキの程度を指標化しているだけ
それに対して、箱ヒゲ図ヒストグラムは、全体のデータの偏り方を観察できる。

また、標準偏差からの乖離の程度を見たり、
あるいは、対応する分布を考えるためにも用いる。

これらは、統計においては、非常に重要なことである。

さて、今回は、この辺りの話には深入りせず、
もう少し、データのバラツキの可視化方法について考えてみる。

「箱ヒゲ図」と「ヒストグラム」の両方の特性を併せ持った可視化方法がある。

ここで紹介する方法は、かなり魅力的。利用価値も高い
それにも関わらず、あまり、出回っていないから不思議である。
未だに、箱ヒゲ図を使った図しか見ない。

悪しきソフトウェア依存の習慣が原因だろうが。
まぁ、そういったことはどうでも良いことである。

さて、今回は以下の3つの方法を紹介する。
  • ヴァイオリン・プロット(violin plot)
  • ビーンズ・プロット(bean plot)
  • ビースウォーム(beeswarm)
基本的には、箱ヒゲ図であるが、ヒストグラムの特徴も持っている。
慣れると、箱ヒゲ図よりも得られる情報量が多く、可読性も高い

これらは、Rの初期状態では使えないパッケージのインストールから始める
Rコンソール上で、次のようなコマンドを入力する。

ダウンロード先の指定の画面が出てきたら、最も近い場所を選ぶ。
別に、どこを選んでも良いのだが、遠いとダウンロードに時間がかかる。

そうそう。ダウンロードだから、インターネットの接続環境が必要

install.packages("vioplot", dep=TRUE)        # ヴァイオリンプロット
install.packages("beanplot", dep=TRUE)     # ビーンプロット
install.packages("beeswarm", dep=TRUE)   # ビースウォーム

インストール作業は、最初の一回のみ。install.packages()関数を用いる。
引数にdep=TRUEとあるが、これは依存関係の確認の有無

Rでは、様々なパッケージを組み合わせるので、
目的のパッケージだけでは動かない場合がある。

TRUEにしておくと、足りないパッケージは自動的にインストール

さて、インストールができたら、今度は、ライブラリを読み込む
ライブラリの読み込みは、毎回必要
読み込まないと使えない。ヘルプも見れない

library(vioplot)
library(beanplot)
library(beeswarm)

では、さっそく。と言いたいところであるが...。

新しく、パッケージをインストールしたら、最初にヘルプを確認すること。

R使いは、ヘルプを見るのが当たり前。いや、見ないと何もできない。
教科書やホームページに頼ってばかりではいけない。

最初は、パッケージついての情報を確認する。

help(package="vioplot")

ここで重要なのは、「Index:」 と書かれている部分から下の部分
左端に「vioplot」とあり、空白の右側に「violin plot」と書いてある。

実は、左端が関数名で、右側はその説明
ふむ。vioplotの説明は、非常に不親切なようである。説明が少ない。

vioplotパッケージには、一つの関数しか含まれていないようであるが、
パッケージによっては、数十もの関数を含む場合もある。

これで、関数名が解ったので、今度はvioplot()関数ヘルプを確認する。

?vioplot        # あるいは、help(vioplot)

最も重要な部分は、「Usage:」の部分。パラメータの確認はここで。
その次に、「Arguments:」の部分で、各パラメータの説明を確認する。
最後に、「Example:」を確認して、使い方を確認する。

まぁ、理解しにくいものや、無理に理解する必要の無いものもあるのだけれど、
その辺りの判断には、ある程度の慣れが必要。

同様に、他のパッケージについても確認。これは各自で。
ひと通り、確認したところで、いよいよ、実際に使ってみる。

今回のデータは、何するか?迷うが、発電所のデータにしてみるか。

このデータは、平成22年度の「地域別発電電力量」のデータ
総務省統計局からダウンロードしたものを編集したもの。詳細は、棒グラフの話。

これをいつものように、

# Windows と Linux の人は次のコマンド
X <- read.table("clipboard", sep=",") 

# Mac の人は次のコマンド
X <- read.table(pipe("pbpaste") , sep=",")

入力して、その後に以下をコピーして、コマンドを実行する。

pnum,max,hydro_pnum,hydro_max,thermal_pnum,thermal_max,atomic_pnum,atomic_max
Hokkaidou,66,7419,53,1234,11,4065,1,2070
Touhoku,228,17206,209,2423,13,11286,2,3274
Tokyo,192,64988,162,8981,25,38696,3,17308
Chubu,197,32828,183,5219,11,23969,1,3617
Hokuriku,138,8057,127,1904,6,4400,1,1746
Kansai,165,34877,149,8196,12,16907,3,9768
Chugoku,110,11986,97,2906,12,7801,1,1280
Shikoku,65,6963,58,1141,4,3797,1,2022
Kyushu,194,20330,139,3279,45,11577,2,5258
Okinawa,22,1919,0,0,21,1919,0,0

前処理として、色を付けるための色の定義を。
これは、棒グラフの話で出てきた色の定義方法。

col.h <- rgb(50, 100, 255, maxColorValue = 255)  # 水力発電所のイメージ色
col.t <- rgb(255, 100, 50, maxColorValue = 255)  # 火力発電所のイメージ色
col.a <- rgb(200, 255, 100, maxColorValue = 255)  # 原子力発電所のイメージ色??

まずは、通常の箱ヒゲ図から。今回は発電所種別ごとの総出力で見てみる。
boxplot(X[,c(4,6,8)], names=c("Hydro","Thermal","Atomic"), col=c(col.h, col.t, col.a))
title(main="Maximum output of electronic plants in Japan", ylab="Maximum output (Kw)")

これは大丈夫のはず。これで、何の変哲もない、普通の箱ヒゲ図が描けた。
次に、ヴァイオリン・プロット。描き方はかなり異なる。
vec.h <- as.vector(X[,4])    # 元のデータから水力発電の総出力をベクトルで抽出
vec.t <- as.vector(X[,6])    # 元のデータから火力発電の総出力をベクトルで抽出
vec.a <- as.vector(X[,8])    # 元のデータから原子力発電の総出力をベクトルで抽出

vioplot(vec.h, vec.t, vec.a, names=c("Hydro","Thermal","Atomic"), col="lightpink")
title(main="Maximum output of electronic plants in Japan", ylab="Maximum output (Kw)")

violot()関数は、ベクトルデータを並べて指定する
そのため、vec.hvec.tvec.a の3つのベクトルを作成し、
それぞれ、水力発電、火力発電、原子力発電の総出力を代入。

実は、色を個別に付けるのは無理なので、ここでは、薄いピンク色。
ソースコードが公開されているので、修正することは難しくないが、
後々になって、説明することが面倒なので、そのままで。

中心の白い「○」が中央値
真ん中の黒い「帯」がIQR(四分位範囲)
真ん中を突っ切る「線」がIQRの±1.5倍のヒゲ

左右に広がる薄いピンクの「山」は、カーネルと呼ばれる密度推定の山
これが、ヒストグラムに相当する。カーネルの話はいずれ。

vioplot()関数は、ヒストグラムの特性と箱ヒゲ図の特性を合わせた特性を持ち、
重要なこととして、箱ヒゲ図で表現されている全て情報も表現している。
問題は、データの指定方法で、この関数だけデータの指定方法が特殊である。

今度は、ビーンプロットを試してみる。

beanplot(X[,c(4,6,8)], names=c("Hydro","Thermal","Atomic"), col="lightblue", 
    main="Maximum output of electronic plants in Japan", ylab="Maximum output (Kw)", cut=0)

beanplot()関数でのデータの選択は、箱ヒゲ図と同じ。

このグラフにおいて、真ん中の太い線は、中央値を表し、
複数の細い「線」は、データが存在している場所を示している。
横に広がる「山」は、vioplot()関数と同様にカーネルの山

vioplot()関数と比較すると beanplot()関数の山は細かく波を打っている。
これは、カーネルのパラメータの違い。カーネルは、様々な描き方ができる。

beaswarm()関数は、他の二つの方法と比較して、
分布に重点を置いているため、より詳細な密度推定が可能となる。
その一方で、箱ヒゲ図と比較して、外れ値やIQRの表現力は低い

注意するべきは、cut パラメータ。この値は「0」にする。
既定のままだと、データが存在しない所まで山裾が広がってしまう

色に関しては、複数のbean(マメ)の色を変えることはできず、
title()関数を使って、後から、タイトルを入れることもできない。

最後はビースウォーム
beeswarm(X[,c(4,6,8)], labels=c("Hydro","Thermal","Atomic"), col=c(col.h, col.t, col.a), pch=16)
title(main="Maximum output of electronic plants in Japan", ylab="Maximum output (Kw)")

beeswarm()関数は、分布をカーネルの山で表現せず、実際の点を布置する。

今回の場合は、対象数が少ないので、見難いが、
対象数がある程度あれば、この方法が最も誤解無く情報を伝えることができる。

boxplot()関数と異なるのは、ラベルの指定パラメータ。
boxplot()関数names であるのに対し、beeswarm()関数では labels となっている。

beeswarm()関数は、他の二つの方法と比較して、推定を用いていないため、
より実態に即した情報を伝えることができるが、外れ値やIQRの表現はできない
ただし、色の指定は可能なので、色分けをしたい場合に使えるというメリットがある。

さてさて。では、今回紹介した方法のうちで、最も適切なものはどれか?
結論から言うと、ケース・バイ・ケースであるが、無難な方法はビースウォーム

うん?色々と欠点があるのに?

実は、ビースウォーム・プロットの欠点は、
箱ヒゲ図と重ね合わせることで克服できる。

boxplot(X[,c(4,6,8)], names=c("Hydro","Thermal","Atomic"))
beeswarm(X[,c(4,6,8)], labels=c("Hydro","Thermal","Atomic"), col="red", pch=16, add=TRUE)
title(main="Maximum output of electronic plants in Japan", ylab="Maximum output (Kw)")

beeswarm()関数のパラメータの最後に add=TRUE としてある。
このようにすることで、前に描画したプロットに、新しいプロットを追加できる

なお、この方法は、他の方法にも適用できるが、
ヴァイオリン・プロットやビーンプロットを箱ヒゲ図に重ねるメリットは無い。