前回までは、主成分分析の仕組みについて整理してきた。
主成分分析というのは、分散共分散行列と相関係数行列と深い関係があって、
これらの見通しをより良くするために、ある種の変換を掛けるであった。
重要なことは、変数が多いと変数間の関係が複雑になりすぎるので、
変数間の関係を本質を変えずに取り除き、対象間の関係のみに焦点を絞ること。
ところで、主成分分析を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主成分で観察し直す必要がある。この話はいずれ。
今回は、主成分負荷量と主成分得点を可視化する方法を実践し、
さらに、その結果を横に並べて表示した。
実は、この二つの図を重ね合わせる可視化方法があって、
その方法のことをバイプロットと呼ぶ。
次回は、そのバイプロットについて整理する。
主成分分析というのは、分散共分散行列と相関係数行列と深い関係があって、
これらの見通しをより良くするために、ある種の変換を掛けるであった。
重要なことは、変数が多いと変数間の関係が複雑になりすぎるので、
変数間の関係を本質を変えずに取り除き、対象間の関係のみに焦点を絞ること。
ところで、主成分分析を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主成分」における「主成分負荷量」を可視化してみる。
今は、量の大きさを可視化したいので、棒グラフを用いる。詳しくは、棒グラフの話。
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
なるほど、確かに、大分県では国際協力に関わるボランティア活動の割合が、
他の都道府県と比較して、群を抜いて高い。では、福井県はどうであるか？
うん？今度は、あまり、何かの値が突出しているようには見えない。
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主成分の主成分負荷量を同時に表現してみる。
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")
主成分負荷量の見方は、棒グラフで見た方法と同じ。要するに、一次元で観察する。
主成分負荷量というのは、あくまで、各次元における主成分の方向の強さなので、
幾何的なベクトルとして考えて矢印で表現する。
今度は、主成分得点の図も出してみる。
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主成分で観察し直す必要がある。この話はいずれ。
今回は、主成分負荷量と主成分得点を可視化する方法を実践し、
さらに、その結果を横に並べて表示した。
実は、この二つの図を重ね合わせる可視化方法があって、
その方法のことをバイプロットと呼ぶ。
次回は、そのバイプロットについて整理する。
以下のコマンドを実行すると、下のような実行結果が現れました。このページで示されている実行結果と違うのですが何故でしょうか？
返信削除> A.svd.score[(order(A.svd.score[,1])),][,1]
[1] -84.3659879 -27.2237009 -17.5878902 -17.4725893 -15.6203251 -15.3652856 -14.7809511 -12.8015618
[9] -11.9976794 -11.9089262 -9.1302220 -8.3305146 -7.8976958 -6.4295937 -5.7279074 -4.9173176
[17] -4.1723528 -3.8796902 -2.8219903 -1.7026102 -1.2568486 -0.8619497 0.3182618 2.6431475
[25] 2.6834044 2.8809400 3.6513143 4.1912833 4.8209478 5.4525580 6.0702498 9.7733572
[33] 9.8101319 13.1450657 13.4666745 13.5595689 13.8047619 14.2913455 14.4812483 15.6424106
[41] 15.8493474 15.8798751 18.3468538 19.7863922 21.2399576 21.3140915 23.1504013
>
文系Aさん、いつも、コメントありがとうございます。最近、バタバタで、返事が遅れております。
削除文系Aさんは、ひょっとして、Mac ユーザでしようか？私のMac上で同様の結果が再現されました。
さて、この値をよく見ると、値そのものは同じなのですが、
どうやら、ラベル（県名）が抜け落ちているようですね。
計算上は、問題は無いと思いますが、解釈するには不便です。
OS依存の問題なのか、あるいは、私が実験している間に、
知らず知らずに、データを上書きしていたのか...。
ちゃんと、確認する必要があります。
現在、諸事情で実験機が落ちているので、
復活し次第、確認してみたいと思います。
もうしばらく、もうしばらくお待ちください。
はい、私はMacユーザーです。
削除では、間違いではないのですね。安心しました。ありがとうございます。
以下のコマンドを実行すると、
返信削除下のようなエラーメッセージが出てきます。
> # 中心に円を描く
> x.circle <- seq(-1, 1, by = 0.01)
> y.circle <- sqrt(1 - x.circle^2)
> lines(x.cercle, y = y.circle)
以下にエラー lines(x.cercle, y = y.circle) :
オブジェクト 'x.cercle' がありません
> lines(x.cercle, y = -y.circle)
以下にエラー lines(x.cercle, y = -y.circle) :
オブジェクト 'x.cercle' がありません
これは、スペルミスですね...。「cercle」 になってます。
削除この部分、正確には「circle」です。
本文も修正しておきました。