Technically Impossible

Lets look at the weak link in your statement. Anything "Technically Impossible" basically means we haven't figured out how yet.

数式なしの主成分分析ーRでの分析から予測までの手順

主成分分析と言うデータ分析手法がある。複数の評価項目のあるデータの全体像をつかむ際に、その評価項目の特性を分析して、より小数の新たな評価項目を見出し、データの全体像をより簡単に把握、可視化しようという試みだ。

例えば、中学生の成績を評価するとしよう。いわゆる5大教科(国、数、英、理、社)に加えて体育、音楽、技術など、それぞれの評価項目から全体から把握するよりも、これらに起因するより簡単な評価項目(評価軸)を見出すことができれば、おぼろげながらも、より簡単な把握が可能になる。例えば、高校入試科目とそれ以外とか、筆記試験を伴う科目と実技試験を伴う科目とか。このような評価軸を統計的に見出だし、リサーチャーが解釈し、データの傾向を把握するのが主成分分析だ。

ここでは主成分分析そのものを紹介するのではなく、Rを用いてその分析手順を紹介していく。分析対象データはRに含まれているものを用いるので、RStudioだけで手順を再現することができる。ここでは次の環境を用いている。

  • R 3.5.2
  • RStudio 1.2.5001

ちなみに現時点でのR、最新版は3.6.2だ。RStudioと共に参照先は、この投稿末尾にまとめている。

ライブラリとデータの準備

この投稿では、次のライブラリを用いる。

FactoMineR 主成分分析に、関数PCAを用いる。
Factoextra 主成分分析の視覚化に用いる。
corrplot 相関行列の視覚化に用いる。
LifeCycleSavings 今回の分析対象データ。

いずれについても、この投稿末尾に参照先をまとめている。

主成分分析に用いる関数はPCA以外にもprcompやprincompなどを用いる場合もあるが、ここではPCAを用いている。端的には、私が気に入っているから、と言うのがその理由だ。特に分析結果を構造的に出力してくれる点を気に入っている。例えば「mypca = PCA(...)」の出力結果の一部は、次のような具合だ。

mypca$eig 固有値
mypca$var 変数(評価項目) $coord 座標
$cor 主成分との相関
$cos2 スコア
$contrib 寄与率
mypca$ind 個々のデータ(分析対象) $coord 座標
$cos2 スコア
$contrib 寄与率
mypca$$ind.sup 個々のデータ(予測対象) $coord 座標
cos2 スコア

LifeCycleSavingsは可処分所得に対する貯蓄率のデータだ。国別に集計されている。変数は次のように定義されている。数値は、1960~1970年の平均だ。

sr 個人の貯蓄率
pop15 15歳未満の人口割合
pop75 75歳以上の人口割合
dpi 可処分所得に対する金額
ddpi dpiの成長率

データの読み込みを含め、これらを利用するためのコードは次のようになる。

library(FactoMineR)
library(factoextra)
library(corrplot)

mydata = LifeCycleSavings

主成分の解釈

まず予測無しに、主成分分析の流れを追ってみる。
関数PCAを用いて、mydataに対する主成分分析を実施する。分析結果をmypcaに格納する。
mypca$eigで、分析結果から各主成分の固有値を表示させると、

mypca = PCA(mydata, scale.unit = TRUE, graph = FALSE)
mypca$eig
       eigenvalue percentage of variance cumulative percentage of variance
comp 1 2.82207781              56.441556                          56.44156
comp 2 1.25606655              25.121331                          81.56289
comp 3 0.60452546              12.090509                          93.65340
comp 4 0.23964496               4.792899                          98.44630
comp 5 0.07768522               1.553704                         100.00000

「cumulative percentage of variance」を参照すると、第1、2主成分で、データ全体の81.6%を説明できることが分かる。

fviz_pca_varで各変数と主成分との相関を確認する。スコア(cos2)、寄与率(contrib)と共に確認してみる。矢印の長さが主成分との相関の表している。矢印が長ければ(円に近ければ)、相関は高い。pop15~dpiは第1主成分との相関が高く、ddpiが第2主成分との相関が高いことが分かる。
ちなみに第1主成分が横軸、第2主成分が縦軸だ。

fviz_pca_var(mypca, col.var = "cos2", gradient.cols = c("blue", "red"), repel = TRUE)
fviz_pca_var(mypca, col.var = "contrib", gradient.cols = c("blue", "red"), repel = TRUE)

f:id:espio999:20191006032315p:plain
f:id:espio999:20191006032331p:plain

具体的な数値で把握するため、corrplotを用いると、第1、2主成分が形成する4象限を解釈しやすい。

corrplot(mypca$var$cos2, addCoef.col = "gray")
corrplot(mypca$var$contrib, is.corr = FALSE, addCoef.col = "gray")

f:id:espio999:20191006032413p:plain
f:id:espio999:20191006032426p:plain

第1主成分左側はpop15、右側はpop75、dpiの領域だ。第1主成分の左から右にかけて可処分所得が大きくなる。左へ行くほど若い人の領域、右に行くほど年寄りの領域、と解釈できる。
第2主成分上側はddpiの領域だ、上に行くほど可処分所得の成長率が高く、下に行くほど低いと解釈できる。

主成分の解釈はリサーチャーに委ねられる。今回のデータは、たまたま素直な出力を返してくれたので、主成分の解釈はかなり率直だった。データによっては解釈に苦しむ出力を返すこともある。考えることが好きな人には面白いところでもある。

個々のデータのプロット

まず主成分上のすべてのデータをプロットしてみる。fviz_pca_indを用いる。

fviz_pca_ind(mypca, repel = TRUE)

f:id:espio999:20191006032445p:plain
大変、見難い。各データを第1、2主成分に反映するにあたり、元のデータから失われている情報がある。各データのスコアを参照し、あまり低いデータは表示しないよう、フィルタをかける。まずスコアを確認する。

fviz_cos2(mypca, choice = "ind", axes = 1:2)

f:id:espio999:20191006032514p:plain

第1、2象限に対するcos2について、0.5を超えるあたりで一区切りするのが妥当に見えるが、それでもまだデータ数は多い。ここではcos2 > 0.8を表示対象にしてみる。かなりすっきりした。

fviz_pca_ind(mypca, repel = TRUE, select.ind = list(cos2 = 0.8))

f:id:espio999:20191006032535p:plain

全体の解釈

fviz_pca_biplotで変数と、個々のデータをプロットしてみる。注意しなければならないのは、個々のデータのポジションと、変数を示す矢印の方向に注目することだ。変数を示す矢印と同じ方向にあるポジションが、変数の傾向を色濃く反映しているデータであり、反対方向にあれば傾向は弱くなる。

fviz_pca_biplot(
  mypca,
  col.var = "blue",
  col.ind = "cos2",
  gradient.cols = c("blue", "red"),
  repel = TRUE, 
  select.ind = list(cos2 = 0.8))

f:id:espio999:20191006032710p:plain

例えばChinaは若い領域におり、可処分所得の金額は控えめなものの、その成長率は高いと解釈できる。60~70年代、まだ中国は今のような超大国ではなかった。
Japanは年寄りの領域におり、可処分所得は高く、その成長率も高い。実際、60~70年代と言えば高度経済成長期だった。

可処分所得が高いものの、成長率の低い領域には北欧、西欧諸国が確認できる。一方、南米諸国は可処分所得が低く、成長率も低いことが分かる。

予測

貯蓄率と日本を除外したデータで主成分分析し、その結果をもとに、それらがどのようにプロットされるのかを確認してみる。つまり、主成分分析の結果に基づいて、それらを予測するということだ。

貯蓄率はデータの1列目、日本はデータの23行目に記録されている。主成分分析の再実施に際し、日本をind.supに、貯蓄率をquanti.supに指定する。

mypca = PCA(mydata, ind.sup = 23, quanti.sup = 1 ,scale.unit = TRUE, graph = FALSE)

PCA実行後、ind.supとquanti.supを確認するとともに、それらをfviz_pca_indにプロットしてみる。

mypca$ind.sup
mypca$quanti.sup
$coord
          Dim.1    Dim.2      Dim.3      Dim.4
Japan 0.3953987 1.620821 0.09006812 -0.8875493

$cos2
           Dim.1     Dim.2       Dim.3     Dim.4
Japan 0.04367951 0.7339682 0.002266467 0.2200859

$dist
   Japan 
1.891892 

$coord
       Dim.1    Dim.2      Dim.3      Dim.4
sr 0.3629719 0.275223 -0.1636127 -0.2067154

$cor
       Dim.1    Dim.2      Dim.3      Dim.4
sr 0.3629719 0.275223 -0.1636127 -0.2067154

$cos2
       Dim.1      Dim.2     Dim.3      Dim.4
sr 0.1317486 0.07574771 0.0267691 0.04273127
myplot = fviz_pca_biplot(
  mypca,
  col.var = "orange",
  col.ind = "cos2",
  gradient.cols = c("gray", "green"),
  repel = TRUE, 
  select.ind = list(cos2 = 0.8))
myplot = fviz_add(myplot, mypca$ind.sup$cos2, color = "blue", repel = TRUE)
myplot

分析結果であるmypcaのindには予測対象が含まれないため、fviz_addでind.supのプロットを指定し、myplotと重ね書きする必要がある。
先ほどとは異なる位置に、貯蓄率srが点線の青矢印で表示され、日本が青でプロットされているのが分かる。
f:id:espio999:20191006051538p:plain

コード

library(FactoMineR)
library(factoextra)
library(corrplot)

mydata = LifeCycleSavings

mypca = PCA(mydata, scale.unit = TRUE, graph = FALSE)
mypca$eig

fviz_pca_var(mypca, col.var = "cos2", gradient.cols = c("blue", "red"), repel = TRUE)
fviz_pca_var(mypca, col.var = "contrib", gradient.cols = c("blue", "red"), repel = TRUE)

corrplot(mypca$var$cos2, addCoef.col = "gray")
corrplot(mypca$var$contrib, is.corr = FALSE, addCoef.col = "gray")

fviz_pca_ind(mypca, repel = TRUE)
fviz_cos2(mypca, choice = "ind", axes = 1:2)
fviz_pca_ind(mypca, repel = TRUE, select.ind = list(cos2 = 0.8))

fviz_pca_biplot(
  mypca,
  col.var = "blue",
  col.ind = "cos2",
  gradient.cols = c("blue", "red"),
  repel = TRUE, 
  select.ind = list(cos2 = 0.8))

mypca = PCA(mydata, ind.sup = 23, quanti.sup = 1 ,scale.unit = TRUE, graph = FALSE)
mypca$ind.sup
mypca$quanti.sup

myplot = fviz_pca_biplot(
  mypca,
  col.var = "orange",
  col.ind = "cos2",
  gradient.cols = c("gray", "green"),
  repel = TRUE, 
  select.ind = list(cos2 = 0.8))
myplot = fviz_add(myplot, mypca$ind.sup$cos2, color = "blue", repel = TRUE)
myplot|