「大人の教養・知識・気付き」を伸ばすブログ

一流の大人(ビジネスマン、政治家、リーダー…)として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。目下、データ分析・語学に力点を置いています。

MENU

Rによるデータサイエンス(04/XX)

 Rについて

をベースに学んでいく。
 今回は主成分分析(PP.60-73)を扱う。

5. 主成分分析

 多くの変数により記述された量的データの変数間にある相関関係を利用し、情報の損失を最小限に抑えつつ、もとの変数の線形結合で表される少数個の無相関な合成変数に縮約して分析する手法を主成分分析という。高次元データをより少数個の変数へ要約して次元圧縮を行い、1次元直線や2次元平面、3次元空間に射影してデータ構造を視覚的に把握するための手法として用いることができる*1

5.1 主成分導出の過程

 一般次元の変数および一般の個数の標本に対する主成分導出の過程を説明しよう。各個体の特徴を捉えるp個の変数を\boldsymbol{X}=(X_1,\cdots,X_p)^{T}とする*2。このp個の変数に関して観測されたn個のp次元データ\boldsymbol{X}_1,\cdots,\boldsymbol{X}_n,\boldsymbol{X}_{i}=(X_{i1},\cdots,X_{ip})^{T},i=1,2,\cdots,nに基づいて標本分散共分散行列


\begin{aligned}
S&=(s_{jk})=\displaystyle{\frac{1}{n}\sum_{i=1}^{n}(\boldsymbol{X}_i-\bar{\boldsymbol{X}})(\boldsymbol{X}_i-\bar{\boldsymbol{X}})^{T}}\\
s_{jk}&=\displaystyle{\frac{1}{n}\sum_{i=1}^{n}(X_{ij}-\bar{X}_j)(X_{ik}-\bar{X}_{k})^{T}}
\end{aligned}

を計算する。ここで\bar{\boldsymbol{X}}=(\bar{X}_1,\cdots,\bar{X}_p)^{T},\bar{X}_i=\displaystyle{\frac{1}{n}\sum_{j=1}^{n}X_{ji}}p次元標本平均ベクトルである。
 まずp個の変数の線形結合で表される射影軸


\begin{aligned}
Y=\omega_1X_1+\cdots+\omega_pX_p=\boldsymbol{\omega}^{T}\boldsymbol{X}
\end{aligned}

上へn個のp次元データを射影して、1次元データY_i=\boldsymbol{\omega}^{T}\boldsymbol{X}_i,i=1,2,\cdots,nを得る。これらの平均\bar{y}\bar{Y}=\displaystyle{\frac{1}{n}\sum_{i=1}^{n}\boldsymbol{\omega}^{T}\boldsymbol{X}_i}=\boldsymbol{\omega}^{T}\bar{\boldsymbol{X}}であるから、その分散s_Y^2


\begin{aligned}
s_Y^2=V[Y]&=\displaystyle{\frac{1}{n}\sum_{i=1}^{n}(Y_i-\bar{Y})^2}\\
&=\displaystyle{\frac{1}{n}\sum_{i=1}^{n}(\boldsymbol{\omega}^{T}\boldsymbol{X}_i-\boldsymbol{\omega}^{T}\bar{\boldsymbol{X}})^2}\\
&=\boldsymbol{\omega}^{T}\left\{\displaystyle{\frac{1}{n}\sum_{i=1}^{n}(\boldsymbol{X}_i-\bar{\boldsymbol{X}})^2}\right\}\boldsymbol{\omega}\\
&=\boldsymbol{\omega}^{T}S\boldsymbol{\omega}
\end{aligned}

で与えられる。したがって射影したデータの分散を最大にするような係数ベクトルは、標本分散共分散行列Sの最も大きい固有値\lambda_1に対応する固有ベクトル\boldsymbol{\omega}_1として与えられる。この\boldsymbol{\omega}_1を係数とする射影軸Y_1=\boldsymbol{\omega}_1^{T}\boldsymbol{X}が第1主成分である。またその射影軸上での分散は最大固有値\lambda_1に等しい。
 第2主成分は第1主成分と直交するという条件の下で射影したp次元データの分散を最大にする軸として与えられ、これは分散共分散行列Sの2番目の大きさの固有値\lambda_2に対応する固有ベクトル\boldsymbol{\omega}_2を係数とする射影軸である。更に第3主成分は、第1、第2主成分と直交するという条件の下で射影したp次元データの分散を最大とする軸として定義される。
 この過程を順次繰り返すことによって理論上、もとの変数の線形結合で表されるp個の主成分が導出される。以上の過程から明らかなように、主成分分析は標本分散共分散行列の固有値問題に帰着できる

5.2 標本分散共分散行列の固有値問題と主成分

 観測されたn個のp次元データに基づく標本分散共分散行列をSとする。標本分散共分散行列はp次元対称行列である。そこでSの固有方程式\mathrm{det}(S-\lambda I_p)=0の解として与えられるp個の固有値


\begin{aligned}
\lambda_1\geq \lambda_2\geq\cdots\lambda_p\geq0
\end{aligned}

とする。また各固有値に対応する、長さが1になるように規格化したp次元固有ベクトルをそれぞれ


\begin{aligned}
\boldsymbol{\omega}_1=\begin{bmatrix}\omega_{11}\\\omega_{12}\\\vdots\\\omega_{1p}\end{bmatrix},\cdots,\boldsymbol{\omega}_p=\begin{bmatrix}\omega_{p1}\\\omega_{p2}\\\vdots\\\omega_{pp}\end{bmatrix}
\end{aligned}

とおく。これらの固有ベクトルは規格化と直交化、すなわち


\begin{aligned}
\|\boldsymbol{\omega}_i\|=1,\boldsymbol{\omega}_i^{T}\boldsymbol{\omega}_j=0,i\neq j
\end{aligned}

が成り立っていることを注意しておく。
 このときもとの変数の線形結合で表されるp個の主成分とその分散は


\begin{aligned}
第1主成分\ \ y_1&=\boldsymbol{\omega}_1^{T}\boldsymbol{X},&V[Y_1]=\lambda_1,\\
第2主成分\ \ y_2&=\boldsymbol{\omega}_2^{T}\boldsymbol{X},&V[Y_2]=\lambda_2,\\
&\vdots&\\
第p主成分\ \ y_p&=\boldsymbol{\omega}_p^{T}\boldsymbol{X},&V[Y_p]=\lambda_p\\
\end{aligned}

で与えられる。
 主成分分析を適用することで、もとのp個の変数について観測されたn個のp次元データ\{\boldsymbol{X}_i=(X_{i1},\cdots,X_{ip})^{T}; i=1,2,\cdots,n\}k(\leq n)次元へと縮約できる。高次元データを1次元直線、2次元平面または3次元空間へと射影してデータ構造を視覚的に把握するとともに、元の変数の線形結合として合成された新たな変数のもつ意味を見出すことによって有益な情報抽出が可能になる。
 合成された主成分の持つ意味は、各変数の係数\omega_{ij}の大きさや政府を通して把握できる。他方で定量的な指標として主成分と各変数の相関を求めておけば、主成分に影響する変数を特定するのに有用である。
 一般に第i主成分Y_iと第j番目の変数X_jとの間の相関は


\begin{aligned}
r_{Y_i,X_j}&=\displaystyle{\frac{\mathrm{Cov}[Y_i,X_j]}{\sqrt{V[Y_i]}\sqrt{V[X_j]}}}\\
&=\displaystyle{\frac{\lambda_i\omega_{ij}}{\sqrt{\lambda_i}\sqrt{s_{jj}}}}\\
&=\displaystyle{\frac{\sqrt{\lambda_i}\omega_{ij}}{\sqrt{s_{jj}}}},i,j=1,2,\cdots,p
\end{aligned}

で与えられる。ここで\lambda_iは第i主成分の分散、\omega_{ij}は第i主成分を構成する変数X_jの係数、s_{jj}は変数X_jの分散である。

5.3 対称行列の固有値固有ベクトル

 標本分散共分散行列


\begin{aligned}
S&=(s_{jk})=\displaystyle{\frac{1}{n}\sum_{i=1}^{n}(\boldsymbol{X}_i-\bar{\boldsymbol{X}})(\boldsymbol{X}_i-\bar{\boldsymbol{X}})^{T}}\\
s_{jk}&=\displaystyle{\frac{1}{n}\sum_{i=1}^{n}(X_{ij}-\bar{X}_j)(X_{ik}-\bar{X}_{k})^{T}}
\end{aligned}

p\times p次対称行列で、その固有値および固有ベクトルの間には以下が成り立つ:


\begin{aligned}
S\boldsymbol{\omega}_i&=\lambda_i\boldsymbol{\omega}_i,\\
\boldsymbol{\omega}_i^{T}\boldsymbol{\omega}_i&=1,\\
\boldsymbol{\omega}_i^{T}\boldsymbol{\omega}_j&=0,i\neq j,i,j=1,2,\cdots,p
\end{aligned}

 いまp個の固有ベクトルを列ベクトルとするp\times p行列をW固有値を対角成分に持つp\times p対角行列を\Lambdaとする。すなわち


\begin{aligned}
W&=(\boldsymbol{\omega}_1,\cdots,\boldsymbol{\omega}_p),\\
\Lambda&=\begin{bmatrix}\lambda_1&0&\cdots&0\\0&\lambda_2&\cdots&0\\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\lambda_p\end{bmatrix}
\end{aligned}

とする。このとき標本分散共分散行列S固有値固有ベクトルの関係は


\begin{aligned}
(1)&\ \ SW=W\Lambda,W^{T}W=I_p,\\
(2)&\ \ W^{T}SW=\Lambda,\\
(3)&\ \ S=W\Lambda W=\displaystyle{\sum_{i=1}^{p}\lambda_i\boldsymbol{\omega}_i\boldsymbol{\omega}_i^{T}},\\
(4)&\ \ \mathrm{tr}S=\mathrm{tr}\left(W\Lambda W^{T}\right)=\mathrm{tr}\Lambda=\displaystyle{\sum_{i=1}^{p}\lambda_i}
\end{aligned}

である。(2)式は対称行列Sは直交行列Wによって対角化可能であることを示しており、(3)式は対称行列Sのスペクトル分解と呼ばれている。また(4)式はもとの変数X_1,\cdots,X_pの分散の総和\mathrm{tr}S=\displaystyle{\sum_{i=1}^{p}s_{ii}}が新たに構成したp個の主成分の分散の総和\mathrm{tr}\Lambda=\displaystyle{\sum_{i=1}^{p}\lambda_i}に等しいことを表す。

5.4 多次元データの基準化と標本相関係数

 観測されたn個のp次元データを


\begin{aligned}
\boldsymbol{X}_i=(X_{i1},\cdots,X_{ip})^T,i=1,2,\cdots,n
\end{aligned}

とするとき、n個のp次元データに基づく標本平均ベクトル\bar{\boldsymbol{X}}=(\bar{X}_1,\cdots,\bar{X}_p)^{T}と標本分散共分散行列S=(s_{jk})を求める。次に上式の第i番目のp次元データを


\begin{aligned}
\boldsymbol{Z}_i=(z_{i1},\cdots,z_{ip})^{T},z_{ij}=\displaystyle{\frac{x_{ij}-\bar{x}_j}{\sqrt{s_{jj}}}},j=1,2,\cdots,p
\end{aligned}

により規格化する。このように規格化されたn個のp次元データに基づいて標本分散(s_{jj})と標本共分散(s_{jk})を計算すると、


\begin{aligned}
s_{jk}^{*}&=\displaystyle{\frac{1}{n}\sum_{i=1}^{n}z_{ij}z_{jk}}\\
&=\displaystyle{\frac{1}{n}\sum_{i=1}^{n}\frac{(x_{ij}-\bar{x}_j)(x_{ik}-\bar{x}_k)}{\sqrt{s_{jj}}\sqrt{s_{kk}}}}\\
&=\displaystyle{\frac{s_{jk}}{\sqrt{s_{jj}}\sqrt{s_{kk}}}}\equiv r_{jk}
\end{aligned}

である。したがって規格化したデータに基づく標本分散共分散行列は、対角成分r_{jj}がすべて1で、非対角成分r_{jk}j番目とk番目の変数間の標本相関係数からなるp次正方行列


\begin{aligned}
R=(r_{jk})=\begin{bmatrix}
1&r_{12}&\cdots&r_{1p}\\
r_{21}&1&\cdots&r_{2p}\\
\vdots&\vdots&\ddots&\vdots\\
r_{p1}&r_{p2}&\cdots&1
\end{bmatrix}
\end{aligned}

に等しく、これを標本相関行列という。
 規格化した多次元データから出発する主成分分析は標本相関行列の固有値問題に帰着され、主成分を導出する過程は標本分散共分散行列から出発する場合と同様である。

5.5 次元圧縮と情報損失

 主成分を用いて高次元空間に散らばるデータの次元を我々の直感で捉えることができる1次元直線、2次元平面または3次元空間へと圧縮することにより、視覚的に高次元データ構造を捉えることができる。ただし高次元データを低次元の空間へと圧縮すると情報の損失が生じることから、次元圧縮による情報損失を定量的に把握しておく必要がある。
 情報損失の例を述べる。2科目の試験において(60点,80点),(80点,60点)というデータをy=x_1+x_2軸へ射影することで圧縮するとどちらも140となるため、区別が出来なくなる。この損失をできるだけ少なくするために射影したデータが最も散らばる軸を探していると言える。
 主成分分析は情報の量を分散で捉えているため、求めた主成分の分散の相対的な大きさを見れば情報の損失を定量的に把握できる。すでに示したように第i主成分の分散V[s_{Y_i}^2]=V[Y_i]=\lambda_iと対応する固有値で表されたから、もとのp個の変数X_1,\cdots,X_pのもつ情報の中で、第i主成分Y_iにはどの程度の情報が含まれている情報の量を測るための基準は、


\begin{aligned}
\displaystyle{\frac{\lambda_i}{\lambda_1+\lambda_2+\cdots+\lambda_p}}
\end{aligned}

で与えられ、寄与率という。また第1主成分から第k主成分までを用いたときの情報の量は


\begin{aligned}
\displaystyle{\frac{\lambda_1+\cdots+\lambda_k}{\lambda_1+\lambda_2+\cdots+\lambda_p}}
\end{aligned}

で与えられ、累積寄与率という。
 固有ベクトルを列ベクトルとするW\in\mathbb{R}^{p\times p}および固有値を対角成分にもつ\Lambda\in\mathbb{R}^{p\times p}、すなわち


\begin{aligned}
W&=(\boldsymbol{\omega}_1,\boldsymbol{\omega}_2,\cdots,\boldsymbol{\omega}_p),\\
\Lambda&=\begin{bmatrix}
\lambda_1&0&\cdots&0\\
0&\lambda_2&\cdots&0\\
\vdots&\vdots\ddots&\vdots\\
0&0&\cdots&\lambda_p
\end{bmatrix}
\end{aligned}

に対して


\begin{aligned}
\mathrm{tr}S=\mathrm{tr}\Lambda
\end{aligned}

が得られるから、もとの変数X_1,X_2,\cdots,X_pのもつ全情報(=各変数の分散の総和\mathrm{tr}S)はp個の主成分Y_1,\cdots,Y_pの全情報\mathrm{tr}\Lambdaへと受け継がれるが、次元圧縮により一部の主成分を用いたときには情報の損失が生じる。そのために用いる主成分にどの程度の情報が含まれているかを定量に測るために、全主成分の分散の総和に対する各主成分の分散の相対的な割合を情報損失の基準としていると言える。

5.6 Rによる分析結果

 データとして\mathrm{iris}を用いて主成分分析を行ってみる。


サンプルデータの分布

# library(stats) # statsは読込不要

# サンプルデータはiris
data(iris)

# まずはサンプルデータの概形を分析
library(ggplot2)
library(ggrepel)

# Speciesごとに統計量を算出
for(species in unique(iris[,"Species"])){
  print(summary(iris[iris[,"Species"]==species,]))
}

# 列ごとに図示
ggplot(iris,aes(x = Species, y = Sepal.Length)) + geom_boxplot()
ggplot(iris,aes(x = Species, y = Sepal.Width)) + geom_boxplot()
ggplot(iris,aes(x = Species, y = Petal.Length)) + geom_boxplot()
ggplot(iris,aes(x = Species, y = Petal.Width)) + geom_boxplot()

# 主成分分析を実施
ls_comp <- princomp(iris[,-5], cor = F, scores = T)

ls_comp$sdev # 標準偏差
ls_comp$loadings # 固有ベクトル:空欄は近似的に0
ls_comp$center # 平均
ls_comp$scale # 各変数に適用されたスケール
ls_comp$n.obs # 観測値のサイズ
ls_comp$scores # 主成分データの得点
summary(ls_comp) # 寄与率・累積寄与率

# 第1成分だけで9割近くの寄与度があるため、第1・第2成分のみで主成分を表示
summary(ls_comp) # 寄与率・累積寄与率

df_prin_comp <- data.frame(names = row.names(ls_comp$loadings),ls_comp$loadings[,c(1,2)])
row.names(df_prin_comp) <- NULL

g <- ggplot(df_prin_comp,aes(x = Comp.1,y = Comp.2, label = names)) + 
  geom_point()
g <- g  + geom_label_repel()+ theme(plot.title = element_text(hjust = 0.5),legend.position = "bottom",
                                    legend.title=element_text(size = 7),
                                    legend.text=element_text(size = 7))
plot(g)

# 主成分得点
# 情報を縮約した結果であり、元の変数・個体間の相対的関係に過ぎない
df_prin_comp_2 <- data.frame(labels = iris[,5],ls_comp$scores)

g_pcs <- ggplot(data = df_prin_comp_2, aes(x = Comp.1, y = Comp.2, colour = labels)) + 
  geom_point()
g_pcs <- g_pcs + theme(plot.title = element_text(hjust = 0.5),legend.position = "bottom",
                       legend.title=element_text(size = 7),
                       legend.text=element_text(size = 7))
plot(g_pcs)

# 主成分の数を判定する:ここでは平行分析法
library(psych)
lst_pc_num <- fa.parallel(x = iris[,-5], fa = "pc", cor = "cov", n.iter = 1000, error.bars = F)

# 結果を図示
df_parallell_analysis <- data.frame("NumComp" = 1:length(lst_pc_num$pc.values),
                                    "DataKind" = rep("RealData",length(lst_pc_num$pc.values)),
                                    "Data" = lst_pc_num$pc.values)
df_parallell_analysis <- rbind(df_parallell_analysis,
                               data.frame("NumComp" = 1:length(lst_pc_num$pc.sim),
                                          "DataKind" = rep("SimData",length(lst_pc_num$pc.sim)),
                                          "Data" = lst_pc_num$pc.sim))
df_parallell_analysis <- rbind(df_parallell_analysis,
                               data.frame("NumComp" = 1:length(lst_pc_num$pc.simr),
                                          "DataKind" = rep("Resampled",length(lst_pc_num$pc.simr)),
                                          "Data" = lst_pc_num$pc.simr))

g_pa <- ggplot(data = df_parallell_analysis,aes(x = NumComp, y = Data, colour = DataKind)) + geom_line()
g_pa <- g_pa + ggtitle("Parallel Analysis Scree Plots") + xlab("Component Number") + 
  ylab("eigen values of principal components")
g_pa <- g_pa + geom_point() + theme(plot.title = element_text(hjust = 0.5),legend.position = "bottom",
                                    legend.title=element_text(size = 7),
                                    legend.text=element_text(size = 7))
plot(g_pa)

# 予測:全データの一部を学習データ、残りをテストデータとして予測する
# 学習データの行番号
vc_resampling <- sample(x = 1:nrow(iris),size = as.integer(round(0.5 * nrow(iris))),replace = F)
ls_comp <- princomp(iris[vc_resampling,-5], cor = F, scores = T) # 学習データでの主成分分析

# テストデータで予測
pred_comp <- predict(object = ls_comp, data = iris[-vc_resampling,-5])
pred_comp


主成分(第1成分と第2成分)の散布図


主成分得点の散布図(第1成分と第2成分)


平行分析法*3

補足 スペック情報

エディション Windows 10 Home
バージョン 20H2
プロセッサ Intel(R) Core(TM) i5-1035G4 CPU @ 1.10GHz 1.50 GHz
実装 RAM 8.00 GB
システムの種類 64 ビット オペレーティング システム、x64 ベース プロセッサ
R バージョン 3.6.3 (2020-02-29)
RStudio バージョン 1.2.5033

*1:以降、5.1から5.5までは小西貞則(2010)「多変量解析入門 線形から非線形へ」(岩波書店) PP.225-236参照。

*2:{}^{T}は転置を表す記号とする。

*3:これで見ると、第1成分のみで十分と言える。

プライバシーポリシー お問い合わせ