主成分分析

---

主成分分析は，高次元（3次元以上）のデータを 2次元のグラフで可視化するのに良く使用される。 場合よっては主成分に分析上有益な意味を見いだすことができる。 多次元データを分析する際は一度この分析手法を試すこと推奨する。

参考：Rのprcomp()関数で主成分分析をするときの注意点

以下は，高校生の成績データの次元縮約を行い，生徒の成績を可視化する例。

データ

出典：「成績データ」『主成分分析』統計科学研究所

d <- read.csv("https://stats.dip.jp/01_ds/data/seiseki_jp.csv")
head(d)

##   学籍番号 国 社 数  理 音 美 体 技 英
## 1       s1 30 43 51  63 60 66 37 44 20
## 2       s2 39 21 49  56 70 72 56 63 16
## 3       s3 29 30 23  57 69 76 33 54  6
## 4       s4 95 87 77 100 77 82 78 96 87
## 5       s5 70 71 78  67 72 82 46 63 44
## 6       s6 67 53 56  61 61 76 70 66 40

library(DT)
datatable(d, caption = "成績データ")

主成分分析

９科目の得点を組み合わせ，できるだけ少ない変数で学生の学力の特徴を調べる。 主成分が何を意味するのか解釈は分析者が考える（実際は意味不明なことが多い）。

summary関数で要約される内容は次のとおり，通常累積割合が8割程度あれば，残りの 主成分は些末であり不要と考え取り除く。このようにして次元縮約ができる。

Importance of components: 成分の重要性

英語名	内容
Standard deviation	標準偏差
Proportion of Variance	変動割合（データ変動を説明する割合=\(\frac{分散}{総分散}\)）
Cumulative Proportion	累積割合（データ変動を説明する割合の累積値）

r <- prcomp(d[, -1], scale = T) # scale = T: 相関行列， F: 分散共分散行列を利用
# 【注意】学籍番号のカラム（1番目）を除いているd[, -1]（科目データだけで分析）

summary(r)

## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.4508 1.0479 0.70060 0.63795 0.54796 0.47059 0.42754
## Proportion of Variance 0.6674 0.1220 0.05454 0.04522 0.03336 0.02461 0.02031
## Cumulative Proportion  0.6674 0.7894 0.84394 0.88916 0.92252 0.94713 0.96744
##                            PC8     PC9
## Standard deviation     0.41376 0.34909
## Proportion of Variance 0.01902 0.01354
## Cumulative Proportion  0.98646 1.00000

# 参考
options(digits = 1) # 表示有効数字2桁
(variance <- r$sdev^2) # 分散（変動），固有値

## [1] 6.0 1.1 0.5 0.4 0.3 0.2 0.2 0.2 0.1

(proportion_variance <- variance / sum(variance)) # 変動割合

## [1] 0.67 0.12 0.05 0.05 0.03 0.02 0.02 0.02 0.01

(cumulative_propotion <- cumsum(proportion_variance)) # 累積変動割合

## [1] 0.7 0.8 0.8 0.9 0.9 0.9 1.0 1.0 1.0

#library(factoextra)
#get_eigenvalue(r)

固有ベクトル（主成分係数）

主成分（PC＊）の固有ベクトルを並べたもの。 主成分得点を計算するときは，変数毎に対応する係数値（表の値）掛けて足し算する。
この表の値は，回帰モデルで例えると，PC1を目的変数，国，社，…を説明変数としたときの 回帰係数に相当する。主成分係数ともいう。

evec <- r$rotation

datatable(round(evec, 2))

主成分得点（スコア）

データ点の主成分座標での値

生徒ごとの主成分の値が分かる。

# レコード名を入力（省くとスコア表とbiplotでレコード連番表示になる）
rownames(r$x) <- d$学籍番号

datatable(round(r$x, 2))

グラフ

library(factoextra)

factoextraパッケージ説明

Scree Plot

データ変動を各主成分が説明する割合％を表すグラフ（scree：がれき，がれ場）
説明割合の低い主成分は無視する →　次元削減
第1主成分と第2主成分だけでデータ変動の8割程度（66.7%＋12.2％）を説明できることが分かる。

fviz_screeplot(r, addlabels = T)

第１主成分への各変数の貢献度（上位５）

第1主成分は，体育を除いた教科がほぼ同程度貢献している。

赤点線は各変数が一様に貢献した場合の％（9科目なので1／9; 11.1%） この線以下の変数は貢献度が低いと言える。

fviz_contrib(r, choice = "var", axes = 1, top = 5)

第２主成分への各変数の貢献度（上位５）

第2主成分は，体育が大きく貢献している。

fviz_contrib(r, choice = "var", axes = 2, top = 5)

主成分負荷量（principal component loading）

変数と主成分との間の相関係数 \(Cor(x_i, y_i)\)で，主成分に寄与する（相関の高い）変数を特定することに使用する。

\[Cor(x_i, y_i) = \sqrt{l_j}h_{ij}\]

ここで，\(x_i\)：\(i\)番目の変数，\(y_j\)：第\(j\)主成分，\(l_j\)：第\(j\)固有値（=分散）， \(h_{ij}\)：主成分係数（\(i\)行，\(j\)列）

library("corrplot")

var <- get_pca_var(r)
corrplot(var$cor, is.corr = T, addCoef.col = "gray") 

#corrplot(var$contrib, is.corr = F, addCoef.col = "gray") 

主成分負荷量ベクトル

主成分軸で主成分負荷量をベクトル表示したもの

理系，文系科目に分かれることを期待したが，実技試験（体育）と筆記試験（その他の教科）になった。

fviz_pca_var(r, 
             col.var = "contrib", # 色分け 
             repel = T) # repel: テキストラベルの重なり防止

Biplot

主成分得点と主成分負荷量ベクトルの2種類を一緒に表示したグラフ

生徒個人の学問的立ち位置（傾向）が分かる。

fviz_pca_biplot(r, col.ind = "contrib", repel = T)

数字画像のグルーピング（参考）

手書きの0〜9までの数字の8×8の画像データを用いて数字のグルーピングを行う例。

【注意】 手書き文字の認識は深層学習が精度的に優れているので実用に際しては， それを利用すること。

データ

出典： Sci-kit learn； 手書きの0〜9までの数字の8×8の画像データ

d0 <- read.csv(file = "https://stats.dip.jp/01_ds/data/hand_writing_numbers0-9.csv")

library(DT)
datatable(d0)

d <- d0[, -1]
number <- d0$number

draw.images <- function(img, i.fr, i.to)
{
  par(mfrow = c(5, 5), 
      mar = c(0, 0, 1, 0)+0.1, 
      cex.main = 0.9)
  
  DX <- 8
  DY <- 8
  BIT <- 16
  
  for (i in i.fr:i.to)
  {
    plot(NA, type = "n",axes = F,
         xlim = c(0, DX),
         ylim = c(0, DY),
         xlab = "",
         ylab = "",
        main = paste("Fig.", i-1))
    
    m <- matrix(unlist(img[i, ])/BIT, nrow = 8, byrow = T)
    
    rasterImage(m, 0, 0, DX, DY) 
  }
}

# 描画
draw.images(img = d, i.fr = 1, i.to = 20)

# 第2主成分までの主成分分析
r <- prcomp(d, rank. = 2)

fviz_screeplot(r, addlabels = T)

# 第１主成分への各変数の貢献度（上位５）
fviz_contrib(r, choice = "var", axes = 1, top = 5)

# 第２主成分への各変数の貢献度（上位５）
fviz_contrib(r, choice = "var", axes = 2, top = 5)

# 変数ベクトル
fviz_pca_var(r, 
             col.var = "contrib", # 色分け 
             repel = T) # repel: テキストラベルの重なり防止

# 個別グラフ（数字を大体グルーピング）
fviz_pca_ind(r, 
             label = d0$number,
             habillage = number,
             addEllipses = T,
             ellipse.level = 0.95)

Python

import pandas as pd
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

d0 = pd.read_csv("https://stats.dip.jp/01_ds/data/hand_writing_numbers0-9.csv")

# 最初の列（正解番号）を除く
d1 = d0.loc[:, d0.columns!="number"]

# 第2主成分まで計算する設定でオブジェクトを作成
pca = PCA(2)

r = pca.fit_transform(d1) # 主成分分析

d = pd.DataFrame(data = r, columns = ['pc1', 'pc2']) # データフレーム化

# 散布図
plt.scatter(d["pc1"], d["pc2"],
            c = d0["number"],
            edgecolor = 'none', alpha = 0.5,
            cmap=plt.cm.get_cmap('Spectral', 10))
            
plt.xlabel('第１主成分')
plt.ylabel('第２主成分')

plt.colorbar()

## <matplotlib.colorbar.Colorbar object at 0x000002BA0EF71C90>

plt.show()

#library(webshot)
#webshot(url = 'https://rpubs.com/tkdhss111/PCA', file = 'PCA.png')
#library(magick) # install.packages('magick')
#image_write(image_read('PCA.png'), 'PCA.pdf', quality = 100)