文系のための「多次元データの要約」（２）

分散共分散行列と相関係数行列を用いることで、様々な定量的手法が可能となる。
実際の手法については「中級レベル」で解説するとして、要するに、
これら行列の意味を理解することは極めて重要である。

そもそも、この二つの行列は、どのようなものであったか？

たしか、分散共分散行列は、各変数の分散が対角成分に並び、
それ以外の成分には「共分散」が入っているような「対称行列」であった。

一方、相関係数行列は、対角成分に「1」が並び、
それ以外の部分には「相関係数」が入っているような行列であった。
相関係数行列は、変数間の関係を一目で確認できるので、色々と重宝するのだった。

ふむふむ。そのような話を確かにした。何か「問題」でもあるのか？

実は、2変数の相関係数の話では気づかなかったのであるが、
多変数の相関係数行列にしてみると、ある種の疑惑が生まれる。

本当に、二つの変数の関係「のみ」を表しているのか？
実は、他の変数の影響を受けているのではないか？
という疑惑である。

例えば、3つの変数、「」、「」、「」の3つの変数が存在した場合、
この3つの変数の相関係数行列は以下のようになるはず。



ふむ。何も不審な点は見当たらないように見えるが、
ひょっとすると、「」と「」の関係の中には、
「」の影響が「潜在的」に存在するかもしれない。

2変数の間の相関係数の話では、意識する必要は無いが、
たしかに、多変数となる相関係数行列では、そういった懸念が生じる。
では、そのような別の変数の影響を取り除くにはどうすれば良いか？

この問題が、本日のメインテーマ。

とりあえず、いつものように、いつもの如く、チュートリアルデータを準備する。
今回のデータは、総務省統計局のデータ（e-stat）のデータを用いて、
チュートリアル用の「擬似」データを作成したもの。

• 居住世帯の有無別住宅数及び建物の種類別住宅以外で人が居住する建物数
• １か月当たり家賃・間代別借家数(平成10年～20年)

今回のデータは、居住世帯の有無と一ヶ月辺りの平均家賃の関係を見るためのデータ。
家賃のデータは、19階級に区分し家賃を集計したものであったので、
各階級の平均を一ヶ月の家賃とし、世帯数を掛けたもので都道府県別の平均家賃を算出。

かなり「大雑把な加工」をしているので、正式な分析には用いてはならない。
あくまで、手法を学ぶための「擬似データ」であるので、
ここから実際の解釈を導いても、「深い意味」はない。

一応、注意事項を整理したところで、いつもどおりにデータを読み込む。

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

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

今回のチュートリアルのために用意したデータは以下の通り。

OCU,UNO,PRI
Hokkaido,2340300,390100,40687
Aomori-ken,493500,87400,111566
Iwate-ken,470700,78800,17002
Miyagi-ken,869700,144200,16084
Akita-ken,380300,57100,41254
Yamagata-ken,383000,49700,10523
Fukushima-ken,699700,108600,11879
Ibaraki-ken,1036200,187600,25207
Tochigi-ken,708700,131300,36003
Gumma-ken,725300,130400,27494
Saitama-ken,2688000,341100,22045
Chiba-ken,2344500,373100,94198
Tokyo-to,5939900,840500,60863
Kanagawa-ken,3612200,455600,246996
Niigata-ken,810700,119000,113336
Toyama-ken,368800,55500,24592
Ishikawa-ken,421600,76400,10485
Fukui-ken,259700,49000,15203
Yamanashi-ken,314600,83600,6005
Nagano-ken,758300,188000,10698
Gifu-ken,712600,123100,25608
Shizuoka-ken,1359400,238500,22521
Aichi-ken,2764400,368400,54336
Mie-ken,680900,110100,122595
Shiga-ken,491300,76300,17786
Kyoto-fu,1086800,183400,14849
Osaka-fu,3685100,660900,39764
Hyogo-ken,2169400,351300,150179
Nara-ken,502500,90100,68843
Wakayama-ken,382100,85700,12267
Tottori-ken,208600,38600,10868
Shimane-ken,249900,45900,6686
Okayama-ken,734700,131900,7459
Hiroshima-ken,1147600,208700,27600
Yamaguchi-ken,584100,107600,48448
Tokushima-ken,297000,58600,21579
Kagawa-ken,372700,73700,10305
Ehime-ken,574000,107100,12877
Kochi-ken,312800,64900,22953
Fukuoka-ken,2034000,340800,11205
Saga-ken,286100,36800,103313
Nagasaki-ken,539200,91900,9725
Kumamoto-ken,663800,105700,20633
Oita-ken,467200,79300,29546
Miyazaki-ken,443800,65900,22566
Kagoshima-ken,718200,133100,17694
Okinawa-ken,504400,62100,29919

このデータにおいて、

• OCUは居住世帯、
• UNOは非居住世帯、
• PRIは一ヶ月の平均家賃。

まずは、分散共分散行列と相関係数行列を出してみる。

cov(X)    # 分散共分散行列
cor(X)    # 相関係数行列

以下が、その結果。

> cov(X)    # 分散共分散行列
             OCU          UNO         PRI
OCU 1.297147e+12 185407240994 21146986138
UNO 1.854072e+11  27363809917  2614407571
PRI 2.114699e+10   2614407571  2151559501

> cor(X)    # 相関係数行列
          OCU       UNO       PRI
OCU 1.0000000 0.9841103 0.4002922
UNO 0.9841103 1.0000000 0.3407284
PRI 0.4002922 0.3407284 1.0000000

さて、相関係数行列を見てみると、
居住世帯数（OCU）と非居住世帯数（UNO）の相関が非常に高く、
平均家賃（PRI）と他の変数の関係は強く無さそうに見える。

相関の有無は、ひとまずは置いておいて、
平均家賃に注目してやると、居住世帯数が多いと賃料が高く、
同様に、非居住世帯数が多くて「も」平均賃料が高くなる傾向がある。

うん？非居住世帯が多くなるほど、月の平均賃料が高くなるというのは、
直感的に合わないような気がする。
居住・非居住の世帯数の関係が影響しているかも。

問題は、2変数間の関係に、別の変数の影響が存在している可能性があること。
では、この「潜在的」な影響は、どのようにしたら見えてくるのか？

ここで、相関係数とは別の方法で2変数の関係を考えてみる。
別の方法とは、すなわち、単回帰と呼ばれる方法である。
何のことか解らない人は、もう一度、やり直し。

たしか、単回帰は、2つの変数の関係について、
一方の変数で、もう一方の変数を「説明」するという考え方だった。
この方法を用いることで、他の変数の潜在的な影響を知ることができるかもしれない。

OCUで「説明できない部分」のPRIの値と、
OCUで「説明できない部分」のUNOの値との、
関係を「比較」すれば良いのではないか？とにかくやってみる。

手を動かした方が解りやすいので、まずは、普通に回帰分析を行なってみる。
Rでは、単回帰を行うためのlm()関数があるので、今回は、この関数を使ってみる。

# 平均家賃を固定し、単回帰を行う
X.lm.3_1 <- lm(X[,3]~X[,1])    # PRIをOCUで説明する
X.lm.2_1 <- lm(X[,2]~X[,1])    # UNOをOCUで説明する

本当は、PRIとUNOを比較したいが、OCUの影響が気になる。
したがって、PRIとUNOをそれぞれ「y軸」に固定し、OCUからの説明を考える。

# グラフを並べて描くためのオマジナイ
par(mfrow=c(1,2))

# PRIをy軸とし、OCUをx軸とした散布図
plot(X[,1], X[,3], xlab="OCU", ylab="PRI", ylim=c(0, max(X[,3])))

# 回帰直線を重ねて描く
par(new=TRUE)
plot(X[,1], X.lm.3_1$fitted.value, type="l", ylim=c(0, max(X[,3])), xaxt="n", yaxt="n", xlab="", ylab="")

# UNOをy軸とし、OCUをx軸とした散布図
plot(X[,1], X[,2], xlab="OCU", ylab="UNO", ylim=c(0, max(X[,2])))

# 回帰直線を重ねて描く
par(new=TRUE)
plot(X[,1], X.lm.2_1$fitted.value, type="l", ylim=c(0, max(X[,2])), xaxt="n", yaxt="n", xlab="", ylab="")

次に、各々の「残差の状況」を確認してみる。

# グラフを並べて描くためのオマジナイ
par(mfrow=c(1,2))

# PRIとOCUの残差の検討（標準化残差）
plot(X[,1], scale(X.lm.3_1$residuals), ylim=c(-4,4), xlab="OCU", ylab="PRI")
abline(h=0)
abline(h=c(-2,2), lty=2, col="grey")

# UNOとOCUの残差の検討（標準化残差）
plot(X[,1], scale(X.lm.2_1$residuals), ylim=c(-4,4), xlab="OCU", ylab="UNO")
abline(h=0)
abline(h=c(-2,2), lty=2, col="grey")

この図では、残差の状況を見やすくするために残差を標準化している。
これは「標準化残差」と呼ばれる。「±2」の線は、おおよその「外れ値」の基準線。
詳しい話は、正規分布の話で述べる。

さてさて、たしか、理屈上は、説明できない部分が「残差」と呼ばれるものだった。
つまり、OCUの影響を取り除くには、「残差同士の相関」を見れば良い。

# 残差同士の相関係数
cor(X.lm.3_1$residuals,X.lm.2_1$residuals)

以下は実行結果。

> # 残差同士の相関係数
> cor(X.lm.3_1$residuals,X.lm.2_1$residuals)
[1] -0.3269775

なるほど。OCUの影響を抜いた状況で相関係数を見てみると、
月の平均賃料と非居住世帯数の関係は、「負の関係」となった。
つまり、非世帯数が増加するほど、平均賃料は低くなる傾向がある。

この結果は、影響を抜かない相関係数よりも直感に即しているように思える。
すなわち、人が住んでいない住宅が多いということは、
需要に対する過剰供給が背後に存在し得るので賃料が低くなる。

試しに、残差同士の関係を散布図で可視化すると以下のようになる。

# とりあえず、変数を定義し直す
pri.res <- X.lm.3_1$residuals
uno.res <- X.lm.2_1$residuals

# 残差同士の単回帰（偏回帰）
res.lm <- lm(uno.res ~ pri.res)

# 残差同士の関係を散布図で表す
plot(pri.res, uno.res, ylim=c(min(uno.res),max(uno.res)), xlab="", ylab="")

# 一応、回帰直線も追加
par(new=TRUE)
plot(pri.res, res.lm$fitted.value, type="l", xlab="", ylab="",  xaxt="n", yaxt="n",  ylim=c(min(uno.res),max(uno.res)), col="red")

# 最後に図を整える
pcor <-  round(cor(pri.res, uno.res),3)
plm <- round(res.lm$coefficient[2],3)
leg1 <- paste("Partial Correlation", pcor, sep=": ")
leg2 <- paste("Partial Regression Coefficient", plm, sep=": ")
legend("topright", legend=c(leg1, leg2))

title(xlab="Residuals for building rental price", ylab="Residuals for unoccupied buildings")
abline(h=0, v=0)

さて、このように、ある2つの変数の関係を観察するために、
別の変数の影響を取り除いた相関係数のことを、
「偏相関係数（Partial Correlation）」と呼ぶ。

また、ここでは、残差同士の単回帰を行い、その回帰直線も描いているが、
残差同士の単回帰の係数のことを「偏回帰係数（Partial Regression Coefficients）」と呼ぶ。
ちなみに、残差同士の関係なので切片は「0」となる。

さて、今回は、単回帰の残差の関係から偏相関を導いたが、
一般的に、偏相関係数は元の相関係数から、以下のように計算できる。



なお、は、の影響を抜いた と との偏相関を意味し、
、、 はそれぞれ、元の相関係数を意味している。
「」は、ギリシア文字で「ロー」と発音する。「ピー」では無い。

さて、要するに、知りたい2つの変数の相関係数から、
潜在的に影響しているかもしれない別の1つの変数の影響を取り除き、
結果の値が、「±1」の間に収まるように分母の部分で基準化している。

今回のように、3変数の場合は、直感的に理解しやすいのであるが、
4変数以上になると計算が非常に厄介になってくる。行列を使えば、
式自体はコンパクトになるが...余裕の無い人は以下の式の話は飛ばしても良い。

まず、元の相関係数行列における、変数「i」と「j」との相関係数を「」 とし、
その相関係数行列の逆行列における「i」と「j」の関係を 「」する。
ちなみに、逆行列の添字を右下ではなく「右上」に置く。ちょっと、注意。

このとき、4変数以上の偏相関係数行列は、次のように求める。



逆行列の対応する要素を、2つの対角要素の平方根で割って基準化し、
さらに、「符号を反転」させている。逆行列を使った数学の魔法!!
この方法で、偏相関係数を計算することができるのである。

まぁ、実際には、手計算で偏相関係数を求める人は少ないであろう。

Rでは、この偏相関係数を簡単に計算できるため、通常はこれを用いるのであるが、
残念ながら、偏相関行列を計算してくれるパッケージは最初からは入っていないので、
「corpcor」パッケージをインストールする。

# corpcor パッケージのインストール
install.packages("corpcor", dep=TRUE)

上手くインストールできたら、実際に計算してみる。

# ライブラリを読み込む
library(corpcor)

# 偏相関係数行列を計算する。
X.pcor <- cor2pcor(cor(X))

# 変数名が付かないので付け直す。
rownames(X.pcor) <- rownames(cor(X))
colnames(X.pcor) <- colnames(cor(X))

以下が、この実行結果。

> X.pcor
          OCU        UNO        PRI
OCU 1.0000000  0.9839439  0.3892439
UNO 0.9839439  1.0000000 -0.3269775
PRI 0.3892439 -0.3269775  1.0000000

PRIとUNOの関係を見てみると、確かに、残差同士の相関係数と一致している。

相関係数行列は、変数間の関係を観察する上で非常に重要であるが、
さらに、細かい状況を把握するためには、
偏相関係数行列で見た方が、状況が良く解る場合もある。

文系のための「主成分分析の可視化」（１）

前回までは、主成分分析の仕組みについて整理してきた。
主成分分析というのは、分散共分散行列と相関係数行列と深い関係があって、
これらの見通しをより良くするために、ある種の変換を掛けるであった。

重要なことは、変数が多いと変数間の関係が複雑になりすぎるので、
変数間の関係を本質を変えずに取り除き、対象間の関係のみに焦点を絞ること。

ところで、主成分分析を2〜3変数のデータに対して用いている例を見かけるが、
実は、少ない変数の場合は、分散共分散行列と相関係数行列を見ても混乱しない。

したがって、散布図行列を観察するだけで十分に解釈できる場合が多い。
状況に合わせて、適切な方法を選択することが重要である。
重回帰分析や偏相関係数を観察するという方法もある。

さて、少し話が逸れたが、本日は、主成分分析の解釈について整理する。
正確には、主成分得点と主成分負荷量を可視化し、
それらから適切な解釈を導く方法について整理する。

では、さっそく、前回までに用いたデータを用いて、実際の処理を始める。
今回もGoogle Spreadsheet にデータを置いてあるので、それを利用する。
なお、今回は、散布図行列を作成するステップは省略する。本当は必要。

library(RCurl)

データの詳細は、散布図行列の話にあるので、そちらを参照する。
ライブラリを無事に読み込めたら、次は、以下のコマンドを実行する。

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

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

主成分分析の結果を可視化する場合には、
バイプロットと呼ばれる方法を用いるのが一般的であり、
prcomp()関数の結果も、その方法で可視化する。

しかしながら、最初からバイプロットで可視化すると、
初学者の多くは、間違った解釈を行う危険性がある。
そこで、本日は、特異値分解から直接的に分析し、その結果の可視化を行う。

ところで、このチュートリアルで用いているデータは、
分散共分散行列と相関係数行列のうち、どちらが適切であったか？

確か、単位が揃っていて、その単位における値の大小にも意味があったので、
分散共分散行列の方が適切であった。前回の話では、そのように解説した。
っということは、つまりは、偏差行列を準備する必要があった。

# 行列の引き算ができるように、平均値を総数個複製する。
M <- matrix(rep(m <- t(colMeans(X)), each=nrow(X)), nrow=nrow(X))

# 元の行列から、偏差行列を作成する。
A <- as.matrix(X-M)

次に、作成された偏差行列に対して特異値分解を行い、
主成分得点と主成分負荷量を計算する。

# 特異値分解を行う
A.svd <- svd(A)

# 主成分スコアと主成分負荷量を計算する
A.svd.score <- A.svd$u %*% diag(A.svd$d)
A.svd.loadings <- A.svd$v %*% diag(A.svd$d)

各主成分の説明力についての情報も欲しい。説明力は何に対応していたか？
確か、主成分得点の行列の分散の割合に相当していた。

# 割合を100分率で表し、さらに、小数点第二位で切り詰める。
A.svd.per <- round(diag((var(A.svd.score)/sum(diag(var(A.svd.score)))))*100,2)

割合なので、小数点以下をダラダラと書いても仕方がない。
小数点第二位〜第三位当たりで切る。この場合は、第二位までで十分。

では、さっそく、「第1主成分」における「主成分負荷量」を可視化してみる。
今は、量の大きさを可視化したいので、棒グラフを用いる。詳しくは、棒グラフの話。

# 第一主成分というのは、主成分の第1次元ということ
dim <- 1

# 各変数の分散の割合の大きさを表すだけなので、とりあえず棒グラフで表現
barplot(as.vector(A.svd.loadings[,dim]), names=colnames(X))

# タイトルは必ず必要。タイトル文字を作成する。少々面倒。
pca.loadings.title <- paste("Factors of principal component in dimension", dim, sep="")
pca.loadings.title <- paste(pca.loadings.title , A.svd.per[1],sep="\n")
pca.loadings.title <- paste(pca.loadings.title , "%", sep="")

# タイトルをグラフに与える。
title(main=pca.loadings.title, ylab="variance")

まずは、主成分負荷量の見方から整理する。

この図では、第1主成分における主成分負荷量は、正と負の値を取っている。
実は、主成分分析の正と負は計算方法に依存しソフトや関数によって逆転する。

したがって、分散の大きさの差と、「0」を中心に二つの相反する性質を表す。
正と負は関係無いので、正負に捕らわれずに判断することが重要。

この図の場合、第1主成分では、国際協力に関する活動と医療や健康に関する活動、
の二つの変数に相当する主成分負荷量が両極にあって、国際協力に近い方から見ていくと、
障害者に関する活動、高齢者に関する活動・・・という順番で並んでいる。

医療や健康に関する活動「のみ」が反対側に向いていて、
しかも、出ている度合いがかなり小さいのが気になるが...とりあえず、次に進む。

こうした状況を観察してみると、第1主成分における二つの対極というのは、
一方がグローバルな活動を表し、もう一方がローカルな活動を表している「らしく」、
この状況は、元のデータの50.79％を説明していることが解る。

「らしく」というのは、この状況だけでは、まだ、なんとも言えない。焦ってはいけない。
主成分得点に関しても同様の図を作成し、
意味については、両方を比較しながら再度検討する。

# 第一主成分というのは、主成分の第1次元ということ
dim <- 1

# 各変数の分散の割合の大きさを表すだけなので、とりあえず棒グラフで表現
barplot(as.vector(A.svd.score[,dim]), names=rownames(X), las=2)

# タイトルは必ず必要。タイトル文字を作成する。少々面倒。
pca.loadings.title <- paste("Scores of principal component in dimension", dim, sep="")
pca.loadings.title <- paste(pca.loadings.title , A.svd.per[dim],sep="\n")
pca.loadings.title <- paste(pca.loadings.title , "%", sep="")

# タイトルをグラフに与える。
title(main=pca.loadings.title, ylab="variance")

第1主成分における主成分得点は、第1主成分における主成分負荷量に対応している。
つまり、この図においては、「0」を中心に下側はグローバルな活動の強さを表し、
一方、上側はローカルな活動の強さを表している。

なお、「0」付近というのは、どちらにも寄らないことを意味するので、
比較的に平均「的」な対象が集まっている。

そうそう、具体的な数値もあった方が良い。

> A.svd.score[(order(A.svd.score[,1])),][,1]
      Oita       Aichi       Miyagi     Ibaraki    Ishikawa    Kanagawa
-84.3659879 -27.2237009 -17.5878902 -17.4725893 -15.6203251 -15.3652856
       Gifu        Nara   Fukushima     Saitama       Ehime       Hyogo
-14.7809511 -12.8015618 -11.9976794 -11.9089262  -9.1302220  -8.3305146
   Kumamoto       Kochi      Toyama    Tokyo-to       Akita     Niigata
 -7.8976958  -6.4295937  -5.7279074  -4.9173176  -4.1723528  -3.8796902
   Shizuoka     Okinawa       Shiga    Osaka-fu     Fukuoka   Tokushima
 -2.8219903  -1.7026102  -1.2568486  -0.8619497   0.3182618   2.6431475
     Nagano       Gumma    Hokkaido      Iwate        Chiba     Okayama
  2.6834044   2.8809400   3.6513143   4.1912833   4.8209478   5.4525580
     Kagawa     Tochigi   Kagoshima   Hiroshima        Saga        Mie
  6.0702498   9.7733572   9.8101319  13.1450657  13.4666745  13.5595689
  Yamanashi     Tottori    Wakayama    Miyazaki    Kyoto-fu   Yamaguchi
 13.8047619  14.2913455  14.4812483  15.6424106  15.8493474  15.8798751
   Nagasaki     Aomori      Shimane    Yamagata       Fukui
 18.3468538  19.7863922  21.2399576  21.3140915  23.1504013

まず、グローバルな活動の強さを見ると、大分県が圧倒的に強く現れており、
次いで、愛知県、宮城県、茨城県、石川県、神奈川県・・・と続く。

一方、ローカルな活動の強さを見ると、福井県が最も値が大きく、
山形県、島根県、青森県、長崎県、山口県・・・と続く。

全体的な印象として、グローバルな活動は、比較的人口が集中している地域が多く、
ローカルな活動は、過疎化などが問題になっている地域が多いように思える。

ここで、気になるのは、グローバルな活動が大分県で高いこと。
このように「なんで？」という疑問は良く生じる。

そういった時に、元データを参照すると新しいことが見えてくる。
常に、元データと分析結果を交互に行き来することが重要。

> 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    15.2 21.0 22.4 13.3 55.2  7.9  8.4 19.9  8.3   8.0
Iwate      6.0 26.1 12.5 15.9 33.1 10.6 15.6 16.9 11.8  26.9
Miyagi    22.8 40.6 32.8 20.7 43.4 13.4 14.3 22.6 15.2  40.4
Akita      4.9 19.5 25.1 16.1 41.2  7.4  9.8 12.8 10.4  33.2
Yamagata   7.4 26.5 13.4 15.8 33.3  9.5 11.4 13.0  5.9   8.7
Fukushima 14.8 30.5 16.2 15.9 40.9 10.5 16.4 29.6 15.2  41.8
Ibaraki    7.9 29.6 38.6 21.4 63.1  9.5 29.6 26.9 10.2  38.5
Tochigi   16.2 36.2 24.3 20.4 51.2 12.1 14.5 24.9  4.7  14.5
Gumma     14.6 32.4 23.1 18.8 41.9  8.7 16.8 26.5  5.9  23.3
Saitama   10.3 36.3 21.8 18.3 34.5 15.0 24.7 28.9  8.5  38.3
Chiba     10.9 29.4 33.1 20.6 42.2 15.5 24.7 22.6  7.1  17.6
Tokyo-to  11.4 41.9 36.8 24.5 37.9 26.1 18.6 48.1  6.6  23.8
Kanagawa  11.9 31.4 26.0 24.6 44.7 16.0 23.7 42.2  5.4  40.7
Niigata   10.3 35.6 36.8 19.0 49.6  8.4 20.6 23.6  5.2  24.8
Toyama    16.8 25.1 28.3 15.6 37.8 11.5 16.5 19.2  9.0  32.7
Ishikawa   6.9 41.7 13.4 25.3 35.0 10.4 23.3 20.3  9.1  43.9
Fukui      7.5 29.3 12.1 17.4 44.2 10.2 15.1 19.6  4.4   5.8
Yamanashi 12.8 27.4 18.2 17.1 39.2  8.7 12.0 25.4  5.4  14.6
Nagano    12.4 24.9 31.8 16.6 29.9  9.3 14.1 18.5  4.9  22.9
Gifu      18.5 29.9 20.7 16.4 44.6  8.2 18.7 26.3  5.2  43.8
Shizuoka  10.6 32.7 18.5 17.8 50.8 11.1 12.8 21.1  5.4  30.9
Aichi      9.0 42.6 21.1 20.2 36.0 11.5 22.1 20.7  5.2  54.2
Mie       23.6 26.8 14.2 17.6 45.4 12.1 15.8 22.6  4.6  16.3
Shiga      9.5 28.0 15.8 22.3 40.9 10.1 13.3 20.7  4.6  31.1
Kyoto-fu  17.6 32.7 22.3 30.3 38.2 18.0 22.4 34.5  5.8   8.3
Osaka-fu  11.9 28.8 41.4 22.2 36.3 13.6 17.5 38.9  7.4  21.2
Hyogo      7.0 40.0 39.8 26.4 42.2 15.6 16.1 23.3  7.8  27.4
Nara      16.3 34.0 36.7 24.9 39.9 11.1 17.2 23.6  7.6  35.0
Wakayama   8.2 38.2 25.1 16.7 39.6 13.2 10.8 33.7  8.2   8.9
Tottori   12.6 31.6 18.3 15.9 41.9  9.2 12.7 23.7  5.1  13.2
Shimane   17.4 28.1 12.3 18.6 31.9  9.9 20.9 20.5  6.6   8.3
Okayama   17.9 39.6 26.9 23.0 32.7 11.9 12.7 27.5  4.6  18.2
Hiroshima 10.0 33.8 25.2 23.1 36.2 12.5 20.7 24.9  4.2  10.8
Yamaguchi 14.2 26.2 31.4 19.7 32.9 11.2 17.5 24.0  5.5   7.9
Tokushima 19.3 26.0 29.8 15.7 46.5  9.7 12.0 17.5  6.9  23.2
Kagawa    18.2 30.8 21.7 20.1 54.5 11.1 12.6 20.7  8.4  20.5
Ehime      9.0 44.3 26.2 18.7 37.5  8.0 13.4 22.5  7.3  33.1
Kochi     19.1 38.3 27.6 19.7 43.0 14.6 11.9 34.7  5.4  30.7
Fukuoka    9.6 29.7 30.5 27.9 44.0 14.4 18.6 24.0 13.3  22.8
Saga      18.9 31.5 27.8 20.7 52.2 10.2 20.0 17.4  5.2  10.2
Nagasaki  26.7 33.4 25.0 18.8 55.1 11.7 10.7 26.9  4.3   6.1
Kumamoto  21.5 38.1 56.9 17.3 30.4 12.5 13.2 20.0  4.2  23.5
Oita      13.4 44.5 58.9 23.8 47.7 12.8 14.5 21.2  6.7 102.6
Miyazaki  13.6 23.8 26.0 20.1 38.7 13.7 21.9 33.0  6.9   9.6
Kagoshima 21.2 26.6 22.2 14.1 40.3 11.1 10.7 18.7  4.7  18.4
Okinawa   21.7 39.9 26.4 22.4 46.6 22.4 20.7 46.0 10.6  24.3

なるほど、確かに、大分県では国際協力に関わるボランティア活動の割合が、
他の都道府県と比較して、群を抜いて高い。では、福井県はどうであるか？

うん？今度は、あまり、何かの値が突出しているようには見えない。

ついでに、データのサマリと箱ヒゲ図も確認。

> boxplot(X, main="Average Days for Participation in Volunteer Activities", ylab="Average days in a year")

> summary(X)
       HM             ELD              HC             CHR    
 Min.   : 4.90   Min.   :19.50   Min.   :12.10   Min.   :13.30
 1st Qu.: 9.80   1st Qu.:27.70   1st Qu.:20.90   1st Qu.:16.90
 Median :12.80   Median :31.50   Median :25.20   Median :19.70
 Mean   :13.80   Mean   :32.35   Mean   :26.59   Mean   :19.85
 3rd Qu.:17.75   3rd Qu.:37.20   3rd Qu.:31.60   3rd Qu.:22.25
 Max.   :26.70   Max.   :44.50   Max.   :58.90   Max.   :30.30
       SC              LI              SP             ENV    
 Min.   :29.90   Min.   : 7.40   Min.   : 8.40   Min.   :12.80
 1st Qu.:36.90   1st Qu.: 9.80   1st Qu.:12.75   1st Qu.:20.40
 Median :41.20   Median :11.20   Median :15.80   Median :23.60
 Mean   :41.89   Mean   :12.05   Mean   :16.53   Mean   :24.98
 3rd Qu.:45.05   3rd Qu.:13.50   3rd Qu.:20.30   3rd Qu.:26.90
 Max.   :63.10   Max.   :26.10   Max.   :29.60   Max.   :48.10
      DIS               IC      
 Min.   : 4.200   Min.   :  5.80
 1st Qu.: 5.200   1st Qu.: 13.85
 Median : 5.900   Median : 23.20
 Mean   : 7.028   Mean   : 25.08
 3rd Qu.: 8.250   3rd Qu.: 32.90
 Max.   :15.200   Max.   :102.60

さて、当初はローカルな活動の強さを表していたのかと思ったが、
福井県のデータは、全体的に各値が低く、特に、国際協力に関する活動は著しく低い。

この傾向は、福井県の次の山形県にも現れている。

したがって、ローカルというよりは「反」グローバル的とするか、
ボランティアに活動に対する総合的な取り組み度合いと見ても良かもしれない。

今度は、第2主成分に関しても同じことをしてみる。
今回は、二つの図を並べて同時に出力する。

#並べて描くためのオマジナイ。
par(mfrow=c(1,2))

# 今度は次元が2になる。
dim <- 2

# 一つ目の図には、主成分負荷量の図
# 各変数の分散の割合の大きさを表すだけなので、とりあえず棒グラフで表現
barplot(as.vector(A.svd.loadings[,dim]), names=colnames(X))

# タイトルは必ず必要。タイトル文字を作成する。少々面倒。
pca.loadings.title <- paste("Factors of principal component in dimension", dim, sep=" ")
pca.loadings.title <- paste(pca.loadings.title , A.svd.per[dim],sep="\n")
pca.loadings.title <- paste(pca.loadings.title , "%", sep="")

# タイトルをグラフに与える。
title(main=pca.loadings.title, ylab="variance")

# 次は、主成分得点の可視化
# 各変数の分散の割合の大きさを表すだけなので、とりあえず棒グラフで表現
barplot(as.vector(A.svd.score[,dim]), names=rownames(X), las=2)

# タイトルは必ず必要。タイトル文字を作成する。少々面倒。
pca.loadings.title <- paste("Scores of principal component in dimension", dim, sep=" ")
pca.loadings.title <- paste(pca.loadings.title , A.svd.per[dim],sep="\n")
pca.loadings.title <- paste(pca.loadings.title , "%", sep="")

# タイトルをグラフに与える。
title(main=pca.loadings.title, ylab="variance")

実際の値も表示する。これは、結果のみ。

> A.svd.score[(order(A.svd.score[,2])),][,2]
    Tokyo-to     Kumamoto     Osaka-fu      Okinawa     Kyoto-fu        Hyogo 
-25.48330834 -19.54068817 -17.89270947 -17.53265633 -11.60217588 -10.09716695 
    Wakayama     Hokkaido     Nagasaki     Miyazaki        Kochi      Niigata 
 -7.99181874  -7.95731799  -7.78153288  -7.29715992  -6.45035170  -6.32440501 
       Chiba     Kanagawa    Yamaguchi      Okayama      Ibaraki         Nara 
 -6.11855778  -6.08651180  -6.01139834  -5.64678869  -5.49745999  -4.43666871 
   Hiroshima      Fukuoka      Tochigi         Saga       Miyagi        Oita  
 -4.06168230  -3.71780644  -3.37521039  -2.60332055  -1.57447800   0.03078564 
       Gumma       Nagano    Tokushima        Ehime       Kagawa      Saitama 
  1.82577203   3.22769906   3.37501834   3.66652229   3.79919709   3.97604392 
     Tottori      Aomori     Yamanashi    Kagoshima       Toyama         Mie  
  4.70810331   5.05798505   5.39623139   6.31970748   6.62959774   6.96469789 
     Shimane        Fukui         Gifu     Shizuoka    Fukushima       Aichi  
  8.14448030   9.14451299   9.99399234  10.06750081  11.39446940  12.86385463 
       Shiga     Yamagata     Ishikawa        Akita       Iwate  
 13.29618683  14.73308634  15.04994751  16.76995612  18.64582591 

第2主成分は、元情報の16.38％の情報を持っている。
主成分負荷量を見ると、障害者に関する活動と国際協力に関する活動が対立している。
全体的に見てみると、恒常的な状況と特殊な状況に分かれているように見える。

一方、主成分得点を見てみると、下側に最も値が高いのが東京都で、
熊本県、大阪府、沖縄県、京都府・・・と続く。
逆に、上側に最も値が高いのは、岩手県で、秋田県、石川県、山形県・・・と続く。

それぞれの上位を見てみると、下側には西日本が多く含まれていて、
上側には東日本が多く含まれているようにも見える。
文化の違いだろうか？それとも、人口密度の問題だろうか？

研究として踏み込むのであれば、より詳細な検討を必要とする。
元データがどのようにして得られたのか、あるいは、外国人の居住人口の割合、
高齢者や障害者の人口比率、犯罪件数、環境汚染といった他の情報も必要とする。
少なくとも、今回のデータでは、これ以上の解釈はできない。

これは、主成分分析に限ったことでは無いが、データの背景を知ることは極めて重要。
また、ある分析から導かれるのが解釈ではなく、新しい疑問であることも多い。

分析結果を安易に鵜呑みにしたり、
強引に分析結果を結びつけることも厳禁！
こういった状況を目にすることは多い。

さて、話が逸れたが、主成分分析の結果を可視化するために、
棒グラフを用いて説明してきたが、これは、各主成分を1次元ごとに観察するため。
次は、二つの次元を同時に可視化する方法を行うが1次元毎に解釈するのは同じ。

では、第1主成分と第2主成分の主成分負荷量を同時に表現してみる。

# 第1主成分と第2主成分の次元を設定する
dim1 <- 1
dim2 <- 2

# 第1主成分と第2主成分の主成分負荷量をベクトルに入れる
x <- A.svd.loadings[,dim1]
y <- A.svd.loadings[,dim2]

# x軸の軸ラベル
label_x <- ""
label_x <- paste("Dim.", dim1)
label_x <- paste(label_x, "(")
label_x <- paste(label_x, A.svd.per[dim1], sep="")
label_x <- paste(label_x, "%)")

# y軸の軸ラベル
label_y <- ""
label_y <- paste("Dim.", dim2)
label_y <- paste(label_y, "(")
label_y <- paste(label_y, A.svd.per[dim2], sep="")
label_y <- paste(label_y, "%)")

# まずは、プロット領域を作成し、ラベルを配置する
plot(x, y, xlab=label_x, ylab=label_y, type="n", main="Loadings of principal component")
text(x, y, colnames(X), cex=0.8)

# 原点0からの方向を矢印で表現する
arrows(x0=0, y0=0, x1=x, y1=y, length=0.1, col="red")

# 最後に、原点軸を入れて整形する
abline(h=0, v=0, lty=2, col="lightgrey")

主成分負荷量の見方は、棒グラフで見た方法と同じ。要するに、一次元で観察する。
主成分負荷量というのは、あくまで、各次元における主成分の方向の強さなので、
幾何的なベクトルとして考えて矢印で表現する。

今度は、主成分得点の図も出してみる。

# 第1主成分と第2主成分の次元を設定する
dim1 <- 1
dim2 <- 2

# 第1主成分と第2主成分の主成分負荷量をベクトルに入れる
x <- A.svd.loadings[,dim1]
y <- A.svd.loadings[,dim2]

# x軸の軸ラベル
label_x <- ""
label_x <- paste("Dim.", dim1)
label_x <- paste(label_x, "(")
label_x <- paste(label_x, A.svd.per[dim1], sep="")
label_x <- paste(label_x, "%)")

# y軸の軸ラベル
label_y <- ""
label_y <- paste("Dim.", dim2)
label_y <- paste(label_y, "(")
label_y <- paste(label_y, A.svd.per[dim2], sep="")
label_y <- paste(label_y, "%)")

# まずは、プロット領域を作成し、ラベルを配置する
plot(x, y, xlab=label_x, ylab=label_y, type="n", main="Loadings of principal component")
text(x, y, colnames(X), cex=0.8)

# 原点0からの方向を矢印で表現する
arrows(x0=0, y0=0, x1=x, y1=y, length=0.1, col="red")

# 最後に、原点軸を入れて整形する
abline(h=0, v=0, lty=2, col="lightgrey")

# 第1主成分と第2主成分の主成分負荷量をベクトルに入れる
x <- A.svd.score[,dim1]
y <- A.svd.score[,dim2]

# x軸の軸ラベル
label_x <- ""
label_x <- paste("Dim.", dim1)
label_x <- paste(label_x, "(")
label_x <- paste(label_x, A.svd.per[dim1], sep="")
label_x <- paste(label_x, "%)")

# y軸の軸ラベル
label_y <- ""
label_y <- paste("Dim.", dim2)
label_y <- paste(label_y, "(")
label_y <- paste(label_y, A.svd.per[dim2], sep="")
label_y <- paste(label_y, "%)")
# まずは、プロット領域を作成し、ラベルを配置する
plot(x, y, xlab=label_x, ylab=label_y, type="n", main="Scores of principal component")
text(x+1, y+1, rownames(X), cex=0.8)

# 原点0からの方向を矢印で表現する
points(x, y, pch=16, col="blue")

# 最後に、原点軸を入れて整形する
abline(h=0, v=0, lty=2, col="lightgrey")

主成分得点は、散布図の「ように」可視化する。あくまで「ように」。
重要なことは、主成分負荷量と同様に、各次元ごとに1次元で解釈すること。

さて、ここでは第1主成分と第2主成分における主成分得点を同時に布置している。
したがって、比較的近い所に布置されている物同士は、
この二つの主成分のみにおいて似た性質を持っていると言える。

また、「0」付近に布置されるもの（今回は存在しない）は、
可視化された二つの軸において、平均「的」な性質を持っていることを意味する。

さて、最後に、相関係数行列に対応した場合でのプロットも作成してみる。
今回の場合、分散共産分散行列に対応した方法が良いのであるが練習だけ。

# 特異値分解を行う
A.svd <- svd(scale(X))

# 主成分スコアと主成分負荷量を計算する
A.svd.score <- A.svd$u %*% diag(A.svd$d)
A.svd.loadings <- (A.svd$v %*% diag(A.svd$d))/sqrt(nrow(X)-1)

# 二つのプロットを並べて描くためのオマジナイ
par(mfrow=c(1,2))

# 第1主成分と第2主成分の次元を設定する
dim1 <- 1
dim2 <- 2

# x軸の軸ラベル
label_x <- ""
label_x <- paste("Dim.", dim1)
label_x <- paste(label_x, "(")
label_x <- paste(label_x, A.svd.per[dim1], sep="")
label_x <- paste(label_x, "%)")

# y軸の軸ラベル
label_y <- ""
label_y <- paste("Dim.", dim2)
label_y <- paste(label_y, "(")
label_y <- paste(label_y, A.svd.per[dim2], sep="")
label_y <- paste(label_y, "%)")

# 主成分負荷量のプロット
# 第1主成分と第2主成分の主成分負荷量をベクトルに入れる
x <- A.svd.loadings[,dim1]
y <- A.svd.loadings[,dim2]

# まずは、プロット領域を作成し、ラベルを配置する
plot(x, y, xlab=label_x, ylab=label_y, type="n", main="Loadings of principal component", xlim=c(-1,1),ylim=c(-1,1),asp=1)
text(x, y, colnames(X), cex=0.8)

# 中心に円を描く
x.circle <- seq(-1, 1, by = 0.01)
y.circle <- sqrt(1 - x.circle^2)
lines(x.circle, y = y.circle)
lines(x.circle, y = -y.circle)

# 原点0からの方向を矢印で表現する
arrows(x0=0, y0=0, x1=x, y1=y, length=0.1, col="red")

# 最後に、原点軸を入れて整形する
abline(h=0, v=0, lty=2, col="lightgrey")

# 主成分得点のプロット
# 第1主成分と第2主成分の主成分負荷量をベクトルに入れる
x <- A.svd.score[,dim1]
y <- A.svd.score[,dim2]

# まずは、プロット領域を作成し、ラベルを配置する
plot(x, y, xlab=label_x, ylab=label_y, type="n", main="Scores of principal component")
text(x+0.1, y+0.1, rownames(X), cex=0.8)

# 原点0からの方向を矢印で表現する
points(x, y, pch=16, col="blue")

# 最後に、原点軸を入れて整形する
abline(h=0, v=0, lty=2, col="lightgrey")

相関係数行列に対応する主成分分析を行う場合、
主成分負荷量は、sqrt(n-1) で基準化する必要がある。

そして、これを可視化する場合には、分散が1であることを示すために半径1の円を描く。
別の言い方をすると、主成分負荷量のプロットの中心に円が描かれている場合には、
その分析結果が相関係数行列に対応していることが解るのである。

さて、今回は第1主成分と第2主成分のみの可視化を行った。
これは、元のデータの67.17%を再現している。
これが十分であるかは、もう少し議論がある。

基本的には、この値が高いほど元のデータに近づくのであるが、
必ずしも、値が高ければ良いとも言えない。

というのは、第1主成分に特定の変数の値が強く反映され、
他の状況が上手く読み取れない場合が存在する。
そのような場合、第1主成分を除いて、第2主成分以下で分析することがある。

また、ある特殊な状況においては、第1主成分と第2主成分の主成分得点のプロットが、
「U」の字に布置されることがあって（馬蹄形問題）、そのような場合には、
第1主成分と第3主成分で観察し直す必要がある。この話はいずれ。

今回は、主成分負荷量と主成分得点を可視化する方法を実践し、
さらに、その結果を横に並べて表示した。
実は、この二つの図を重ね合わせる可視化方法があって、
その方法のことをバイプロットと呼ぶ。

次回は、そのバイプロットについて整理する。