主成分分析

1 はじめに

＜主成分分析の考え＞
　物事は見る方向によって見え方が変わります。鳥になって空から眺めれば、地上の全体像が分かります。
地上にいて見回しても、さまざまなものが重なり合って、手前のものしか見えません。
２つの対象を区別するには、２つが離れて見えることが重要です。
離れて見えるということは、その方向から見たデータの分散が大きいということです。
　「分散が最大となる方向から見る」
これが主成分分析の考え方です。

＜本書の特徴＞
　主成分分析には線形代数の知識が必要とされ、高校での扱いは紹介にとどまっています。
「数理を理解してから応用を学ぶ」というのが従来のスタイルであるからです。
本書では、数理的事実を簡潔に提示し、シミュレーションで確かめるという方法をとります。
これは、統計数理をシミュレーションによる実験・実証で代替・補完しようという試みです。

＜自分で確かめる学習＞
　高校生のみなさんは、グラフの観察や関数による計算を活用して主成分分析に挑んでください。実際の計算は R の関数がやってくれます。様々なモデルにおいて、Rの関数がどのように働くか確かめましょう。
尚、将来の学習のため、行列を用いた数理的説明も簡潔に記しておきます。

　乱数シミュレーションを用いた実証的学習方法を、筆者は「ひょっこり方式」と呼んでいます。この学習方法が、皆さんの探究的学習の一助となれば幸いです。

2 二変量正規分布

簡単のため2変量データから始めましょう。
2変量正規乱数を用いると、指定した平均と分散共分散行列を持つ2変量データを生成することができます。

2.1 二変量正規乱数

• \(I=\begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}\)、\(\mu=\begin{pmatrix} 0 \\ 0 \end{pmatrix}\)　である正規乱数 data1.
  
• \(\Sigma=\begin{pmatrix} 1 & 0.7 \\ 0.7 & 1 \end{pmatrix}\)、\(\mu=\begin{pmatrix} 0 \\ 0 \end{pmatrix}\)　である正規乱数 data.

を各々1000個出力します。
（以降、平均は \(\mu=0\) としておきます）

library(mvtnorm)
n <- 1000
I <- matrix(c(1,0,0,1), ncol=2)
z <- rmvnorm(n, mean=c(0,0), sigma=I)  
data1 <- data.frame(X=z[,1],Y=z[,2])
Sigma <- matrix(c(1,0.7,0.7,1), ncol=2)
x <- rmvnorm(n, mean=c(0,0), sigma=Sigma)  
data <- data.frame(X=x[,1],Y=x[,2])

• data1 を左図にプロットします。
  （data1 は2変量標準正規乱数）
  
• data を右図にプロットします。
  （回帰直線 y=0.7x を赤線で描きます）

par(mfrow=c(1,2))
plot(data1,asp=1,col="lightblue",main="data1",xlim=c(-4,4))
grid()
plot(data,asp=1,col="lightblue",main="data",xlim=c(-4,4))
abline(a=0,b=0.7,col="red")
grid()

グラフ平面上でこれらのデータを観察するとしましょう。

• 左図は、どの方向から観察しても変わりはなく、ほぼ円状に分布しています。 （2変量標準正規分布）
  
• 右図は、円状分布をある方向に拡大したように見えます。
  しかし、その方向は y=0.7x 方向ではないようです。
  どの方向に拡大しているのでしょうか？

拡大している方向を第１主成分方向といい、分散共分散行列の固有ベクトル方向がこれに相当します。
（固有ベクトルは複数あるので、最大固有値に対する固有ベクトル方向です）.

2次正方行列で表される一次変換の固有ベクトルは、次のように定義されます。

2.2 固有値・固有ベクトル

　　　行列 \(A=\begin{pmatrix} a & b \\ c & d \end{pmatrix}\) に対し、 \(A\begin{pmatrix} x \\ y \end{pmatrix}=\lambda \begin{pmatrix} x \\ y \end{pmatrix}\)
　　　を満たす \(\lambda\) を固有値、\(\begin{pmatrix} x \\ y \end{pmatrix}\) を固有ベクトルといいます。

• 固有ベクトルは、その一次変換によって方向が変わらないベクトルをいいます。
  
• 固有値は、その方向に沿った拡大・縮小の倍率を表します。

固有ベクトルは、eigen( )関数で求めることができます。

• 分散共分散行列 \(\Sigma=\begin{pmatrix} 1 & 0.7 \\ 0.7 & 1 \end{pmatrix}\)　に eigen( )関数を適用します。

eigen(Sigma)

## eigen() decomposition
## $values
## [1] 1.7 0.3
## 
## $vectors
##           [,1]       [,2]
## [1,] 0.7071068 -0.7071068
## [2,] 0.7071068  0.7071068

• $values の値 (1.7,0.3) は固有値です。
• $vectors行列の第１列が第１固有ベクトルで、方向は y= x （第１主成分方向）
  
• $vectors行列の第２列が第２固有ベクトルで、方向は y=-x （第２主成分方向）

固有値の平方根が、拡大縮小の倍率になっています。
固有ベクトルは正規化され、大きさは 1 になっています。

［注意］

\(\overrightarrow{p}\) が \(\Sigma\)の固有ベクトルのとき、その定数倍 \(a\overrightarrow{p}\) も固有ベクトルです。
そのため、eigen( )関数で出力される固有ベクトルは、大きさが１に正規化されています。

\(\overrightarrow{p}\) と \(-\overrightarrow{p}\) どちらが出力されるかは関数に委ねられています。
prcomp( ) 関数では、$rotatin として固有ベクトルが出力されますが、その向きは一意ではありません。

出力される固有ベクトルの向きによって、主成分軸の正方向が変わることに注意しましょう。

2.3 拡大縮小行列

左図のdata1 を、x軸方向に \(\sqrt{1.7}\) 倍、y軸方向に \(\sqrt{0.3}\) 倍、
固有値の平方根倍に拡大縮小したデータを d1_new として左に図示します。

拡大縮小の行列は、

• \(A=\begin{pmatrix} \sqrt{1.7} & 0 \\ 0 & \sqrt{0.3} \end{pmatrix}\)　で表されます。

A <- matrix(c(sqrt(1.7),0,0,sqrt(0.3)), ncol=2)
d1_new <- A %*% t(z)
d1_new <- data.frame(X=d1_new[1,],Y=d1_new[2,])
par(mfrow=c(1,2))
plot(d1_new,asp=1,col="lightblue",main="d1_new")
abline(h=0,col="blue")
abline(v=0,col="pink")
plot(data,asp=1,col="lightblue",main="data",xlim=c(-4,4))
abline(a=0,b=1,col="blue")
abline(a=0,b=-1,col="pink")

左右の図は同じ形状の楕円分布ですが、長軸・単軸の方向が異なっています。
右図の長軸は第１主成分方向 y=x、単軸は第２主成分方向 y=-x です。

従って右図は、2変量標準正規分布 data1の円状の分布を、

• 第１主成分（青）方向に \(\sqrt{1.7}\) 倍に拡大、
  
• 第２主成分（桃）方向に \(\sqrt{0.3}\) 倍に縮小した楕円状の分布

と考えることができます。

＜主成分方向の分散＞

　１変量では　\(V(kX)=k^2V(X)\)　　（\(k\) 倍すると分散は \(k^2\) 倍）

この類比で考えると、2変量でも同様のことが成立すると予想されます。

2変量標準正規分布 data1 の分散は、原点を通るどの方向に対しても 1ですが、
右図 data の主成分方向の分散については、

• 第１主成分方向の分散は 1 の 1.7 倍
• 第２主成分方向の分散は 1 の 0.3 倍

固有ベクトルの固有値が、その方向の分散になっていると考えられます。 　

［まとめ］

　\(N_2(0, \Sigma\)）の分布の特徴を整理しておきます。

• \(\Sigma\) の固有ベクトルが主成分方向
  （長軸が第１主成分方向、短軸が第２主成分方向）
  
• \(\Sigma\) の固有値が主成分方向の分散
  （大きい固有値が第１主成分方向、小さい固有値が第２主成分方向）

3 主成分座標

二変量正規分布 \(N_2(0, \Sigma\)）は、楕円状の分布をしています。
楕円の（長軸, 短軸）を座標軸とした新しい座標を、主成分座標といいます。
長軸座標 PC1 が第１主成分座標、短軸座標 PC2 が第２主成分座標です。

分散最大の方向の座標が PC1、分散最小の方向の座標が PC2 です。
主成分座標は、主成分得点（スコア）とも呼ばれます。
（数理については、高校生の線形代数「5.内積と正射影」参照）

新座標 (PC1, PC2) を求める数理は後述しますが、Rではprcomp( )関数を用いて算出することができます。

3.1 内積の定義

ベクトル \(\overrightarrow{a}=\begin{pmatrix} a_1 \\ a_2 \end{pmatrix}\)　 \(\overrightarrow{b}=\begin{pmatrix} b_1 \\ b_2 \end{pmatrix}\)　の内積は、

• \(\overrightarrow{a}\cdot \overrightarrow{b}=|\overrightarrow{a}||\overrightarrow{b}|\cos\theta\)　（なす角\(\theta\)による定義）

• \(\overrightarrow{a}\cdot \overrightarrow{b}=a_1b_1+a_2b_2\)　（成分による定義）

この２つの定義は余弦定理によって相互に繋がっています。

特に、\(\overrightarrow{p}\) が単位ベクトルのとき、

• \(\overrightarrow{p}\cdot \overrightarrow{x}=|\overrightarrow{x}|\cos\theta\)は、\(\overrightarrow{x}\)の\(\overrightarrow{p}\)への正射影の座標を表します。
  正射影の座標は、\(\overrightarrow{p}\) を座標軸としたときの\(\overrightarrow{x}\)の座標に相当します。

3.2 主成分座標の計算

主成分座標は、固有ベクトルを基底とする座標であり、固有ベクトルへの正射影です。
固有ベクトルへの正射影は固有ベクトルとの内積で求まります。

計算方法は、次のようになります。
\(\Sigma\) の固有ベクトルについて、固有値の大きい順に
　　　\(p=\begin{pmatrix} p_1\\ p_2 \end{pmatrix}\)　\(q=\begin{pmatrix} q_1\\ q_2 \end{pmatrix}\)　　（ \(|p|=1,|q|=1\)）
とすると、\(p\) は第１主成分方向、\(q\) は第２主成分方向の単位ベクトルです。

データ \(x=\begin{pmatrix} x_1\\ x_2 \end{pmatrix}\) に対し、

• p方向への正射影の座標は　\(\overrightarrow{p} \cdot \overrightarrow{x}=p_1x_1+p_2x_2\) （ \(p\)と\(x\)の内積）
• q方向への正射影の座標は　\(\overrightarrow{q} \cdot \overrightarrow{x}=q_1x_1+q_2x_2\) （ \(q\)と\(x\)の内積）

主成分座標は、各固有ベクトルへの正射影の座標なので、

• x の第１主成分座標は、\(PC1=p_1x_1+p_2x_2\)

• x の第２主成分座標は、\(PC2=q_1x_1+q_2x_2\)

主成分を構成する係数は、分散共分散行列の固有ベクトルの成分です。
これを主成分係数（Principal component coefficient）といいます。
つまり、主成分係数とは変量から主成分を求める係数のことです。

この式は、固有ベクトルを基底とする座標への座標変換を意味します。

3.3 主成分プロット

prcomp( )関数を用いて主成分座標をプロットしましょう。
基本的な流れは次のとおりです。

• \(N_2(0, \Sigma)\)に従う data を生成.
  標本分散共分散行列から、eigen関数で固有値と固有ベクトルを算出.
  

• dataを標準化しないでprcomp( )関数を適用（scale = FALSE）.
  固有ベクトル ve と主成分座標 data_PC を抽出.

• 左図に data をプロット（固有ベクトルは青、赤で表示）.
  右図に主成分座標 data_PC をプロット.

主成分プロットは、固有ベクトルの向きに応じ上下反転・左右反転することがあります。

3.3.1 大標本データ

• \(\Sigma=\begin{pmatrix} 1 & 0.7 \\ 0.7 & 1 \end{pmatrix}\)、\(\mu=\begin{pmatrix} 0 \\ 0 \end{pmatrix}\)　である正規乱数 data.

Sigma <- matrix(c(1,0.7,0.7,1), ncol=2)
x <- rmvnorm(500, mean=c(0,0), sigma=Sigma)  
data <- data.frame(X=x[,1],Y=x[,2])
var(data)

##           X         Y
## X 1.0079280 0.6818276
## Y 0.6818276 1.0223493

標本分散共分散行列は \(\Sigma\) に近いですが、完全に一致するわけではありません。

pca_result <- prcomp(data, scale = FALSE)
ve <-pca_result$rotation#固有ベクトル
ve

##          PC1        PC2
## X -0.7033581  0.7108357
## Y -0.7108357 -0.7033581

固有ベクトルの方向が、座標軸の正方向になります。

主成分座標を抽出します。

data_PC <- pca_result$x
head(data_PC)

##             PC1        PC2
## [1,] -1.4418490 -1.1397085
## [2,] -0.7087794 -0.3441445
## [3,] -0.2992616  0.5705260
## [4,] -0.4160481  1.0365751
## [5,]  0.9282066 -0.5823057
## [6,] -0.2099791  0.3428498

data（左図）、data_PC（右図）をプロットします。

par(mfrow=c(1,2))
plot(data,asp=1,col="lightblue",main="data",xlim=c(-4,4))
abline(a=0,b=ve[2,1]/ve[1,1],col="lightblue")
abline(a=0,b=ve[2,2]/ve[1,2],col="pink")
arrows(x0=0, y0=0, x1=ve[1,1],  y1=ve[2,1],col="blue",length = 0.1)
arrows(x0=0, y0=0, x1=ve[1,2],  y1=ve[2,2],col="red",length = 0.1)
plot(pca_result$x,,asp=1,col="lightblue",main="data_PC",xlim=c(-4,4))
abline(h=0,col="lightblue")
abline(v=0,col="pink")
arrows(x0=0, y0=0, x1=1,  y1=0, code=2, col="blue",length = 0.1)
arrows(x0=0, y0=0, x1=0,  y1=1, code=2, col="red", length = 0.1)

右図の data_PC は、固有ベクトルの向きに応じ左右反転・上下反転します。

主成分座標は、dataの長軸をx軸、短軸をy軸とした座標です。

3.3.2 小標本データ

小標本のほうが、主成分座標の分布や反転の状況がよく分かります。

データ数を20として data（左図）、data_PC（右図）をプロットします。

Sigma <- matrix(c(1,0.7,0.7,1), ncol=2)
x <- rmvnorm(20, mean=c(0,0), sigma=Sigma)  
data <- data.frame(X=x[,1],Y=x[,2])
pca_result <- prcomp(data, scale = FALSE)
ve <-pca_result$rotation
data_PC <- pca_result$x
par(mfrow=c(1,2))
plot(data,asp=1,col="lightblue",main="data",xlim=c(-4,4))
arrows(x0=0, y0=0, x1=ve[1,1],  y1=ve[2,1], code=2, col="blue", length = 0.1,lwd = 2)
arrows(x0=0, y0=0, x1=ve[1,2],  y1=ve[2,2], code=2, col="red", length = 0.1,lwd = 2)
abline(a=0,b=ve[2,1]/ve[1,1],col="lightblue")
abline(a=0,b=ve[2,2]/ve[1,2],col="pink")
plot(data_PC,asp=1,col="lightblue",xlim=c(-4,4),main="data_PC")
abline(h=0,col="lightblue")
abline(v=0,col="pink")
arrows(x0=0, y0=0, x1=1,  y1=0, col="blue",length = 0.1)
arrows(x0=0, y0=0, x1=0,  y1=1, col="red", length = 0.1)

小標本の方が分布の特徴が見えやすくなります。

data_PC は、固有ベクトルの向きに応じ左右反転・上下反転します。

3.3.3 データ番号表示

データ番号を表示することも、dataプロットとdata_PCプロットの対応の理解に役立ちます。

x <- rmvnorm(20, mean=c(0,0), sigma=Sigma)  
data <- data.frame(X=x[,1],Y=x[,2])
pca_result <- prcomp(data, scale = FALSE)
ve <-pca_result$rotation
par(mfrow=c(1,2))
plot(data,asp=1,type="n",xlim=c(-3,3),main="data")
text(1:20,x=data[,1],y=data[,2],col="lightblue")
arrows(x0=0, y0=0, x1=ve[1,1],  y1=ve[2,1], code=2, col="blue", length = 0.1,lwd = 2)
arrows(x0=0, y0=0, x1=ve[1,2],  y1=ve[2,2], code=2, col="red", length = 0.1,lwd = 2)
abline(a=0,b=ve[2,1]/ve[1,1],col="lightblue")
abline(a=0,b=ve[2,2]/ve[1,2],col="pink")
plot(pca_result$x,,asp=1,type="n",xlim=c(-3,3),main="data_PC")
text(1:20,x=pca_result$x[,1],y=pca_result$x[,2],col="lightblue")
abline(h=0,col="blue")
abline(v=0,col="pink")
arrows(x0=0, y0=0, x1=1,  y1=0, col="blue",length = 0.1)
arrows(x0=0, y0=0, x1=0,  y1=1, col="red", length = 0.1)

data_PC は、固有ベクトルの向きに応じ左右反転・上下反転します。

4 負荷量

負荷量は、主成分分析や因子分析において用いられる用語です。 私見ですが、その意味について考察します。

負荷量は、荷を負わせる量です。
複数の原因が結果に関係するとき、原因が結果に及ぼす影響量の荷重が負荷量です。

ひとつの例えですが、
（父, 母）という複数の原因によって、（子供）という結果が出力されると考えましょう。

• （子供）= 0.5（父）+ 0.5（母）
  （子供）に対する負荷量は、（父, 母）＝（0.5, 0.5）

数式としては、その原因に乗ぜられる係数が負荷量です。
この負荷量からは、男女平等であることが感じられます。

（子供）＝（男, 女）という多変量で分析しましょう。

• （男）= 0.8（父）+ 0.6（母）
  （女）= 0.6（父）+ 0.8（母）

であったとしましょう。この場合、

• （男）に対する負荷量は、（父, 母）＝（0.8, 0.6）
  （女）に対する負荷量は、（父, 母）＝（0.6, 0.8）

父親は男の子供により大きな影響を与え、母親は女の子供により大きな影響を与えるように感じられます。

［注意］
　結果と原因が変われば負荷量も変わります。
「3.2 主成分座標の計算」の \((x_1,x_2)\) から \((PC1,PC2)\) を求める式に注目します。

• \(PC1=p_1x_1+p_2x_2\)　では、
  \(PC1\) に対する \((x_1,x_2)\) の負荷量は \((p_1,p_2)\)
  
• \(PC2=q_1x_1+q_2x_2\)　では、
  \(PC2\) に対する \((x_1,x_2)\) の負荷量は \((q_1,q_2)\)

このように、何を結果とするかによって負荷量は変わります。

分析の方向を変えてみましょう。

4.1 主成分負荷量

\(N_2(0, \Sigma\)）について、「3.2 主成分座標の計算」の設定を継続します。

\(\Sigma\) の固有値を大きい順に \(\lambda_1,~\lambda_2\)、対応する固有ベクトルを

　　　\(p=\begin{pmatrix} p_1\\ p_2 \end{pmatrix}\)　\(q=\begin{pmatrix} q_1\\ q_2 \end{pmatrix}\)　　（ \(|p|=1,|q|=1\)）

とすると、\(p\) は第１主成分方向、\(q\) は第２主成分方向の単位ベクトルです。

データ \(x=\begin{pmatrix} x_1\\ x_2 \end{pmatrix}\) に対し、

　　　\(\begin{pmatrix} PC_1\\ PC_2 \end{pmatrix} =\begin{pmatrix} p_1 & p_2 \\ q_1 & q_2 \end{pmatrix}\begin{pmatrix} x_1 \\ x_2 \end{pmatrix}\)

係数は、\((PC1,PC2)\) を求めるための \((x_1,x_2)\) の負荷量です。

分析の方向を変えます。

分散共分散行列の固有ベクトルは直交しており、係数行列は直交行列です。
従って、逆行列は転置行列として直ちに求まり、

　　　\(\begin{pmatrix} x_1\\ x_2 \end{pmatrix} =\begin{pmatrix} p_1 & q_1 \\ p_2 & q_2 \end{pmatrix}\begin{pmatrix} PC_1 \\ PC_2 \end{pmatrix}\)　

係数は、\((x_1,x_2)\) を求めるための \((PC1,PC2)\) の負荷量で、主成分負荷量といいます。
係数行列は固有ベクトルを列とする直交行列であり、rotation に一致します。

• \(x_1\) に対する \((PC1,PC2)\) の主成分負荷量は \((p_1,q_1)\)
  
• \(x_2\) に対する \((PC1,PC2)\) の主成分負荷量は \((p_2,q_2)\)

これらは、変量を構成する主成分の割合と考えることができます。
後述の biplot では、\(x_1,~x_2\) の主成分負荷量としてプロットされます。

4.2 負荷量ベクトル

　 \(\begin{pmatrix} x_1\\ x_2 \end{pmatrix} = \begin{pmatrix} p_1\\ p_2 \end{pmatrix}PC1 + \begin{pmatrix} q_1\\ q_2 \end{pmatrix}PC2\)　と考えると、

• \(\begin{pmatrix} x_1\\ x_2 \end{pmatrix}\) に対する \(PC1\) の負荷量ベクトルは \(\begin{pmatrix} p_1\\ p_2 \end{pmatrix}\)

• \(\begin{pmatrix} x_1\\ x_2 \end{pmatrix}\) に対する \(PC2\) の負荷量ベクトルは \(\begin{pmatrix} q_1\\ q_2 \end{pmatrix}\)

負荷量ベクトルの大きさはどちらも１であり、係数ベクトルに主成分の影響力の強さは現れません。
\(PC1\) の分散は \(\lambda_1\)、\(PC2\) の分散は \(\lambda_2\) であり、影響力は主成分の値のなかに含まれています。

そこで、
主成分の分散を１に直すという工夫をすると、係数ベクトルに主成分の影響力が現れます。

　 　\(\begin{pmatrix} x_1\\ x_2 \end{pmatrix} = \begin{pmatrix} \sqrt{\lambda_1}p_1\\ \sqrt{\lambda_1}p_2 \end{pmatrix}\dfrac{PC1}{\sqrt{\lambda_1}} + \begin{pmatrix} \sqrt{\lambda_2}q_1\\ \sqrt{\lambda_2}q_2 \end{pmatrix}\dfrac{PC2}{\sqrt{\lambda_2}}\)　

\(z_1=\dfrac{PC1}{\sqrt{\lambda_1}}\)、\(z_2=\dfrac{PC2}{\sqrt{\lambda_2}}\) とおくと、

　 　\(\begin{pmatrix} x_1\\ x_2 \end{pmatrix} = \begin{pmatrix} \sqrt{\lambda_1}p_1\\ \sqrt{\lambda_1}p_2 \end{pmatrix}z_1 + \begin{pmatrix} \sqrt{\lambda_2}q_1\\ \sqrt{\lambda_2}q_2 \end{pmatrix}z_2\)　

　 　\(\begin{pmatrix} x_1\\ x_2 \end{pmatrix} = \begin{pmatrix} \sqrt{\lambda_1}p_1 & \sqrt{\lambda_2}q_1\\ \sqrt{\lambda_1}p_2 & \sqrt{\lambda_2}q_2 \end{pmatrix} \begin{pmatrix} z_1\\ z_2 \end{pmatrix}\)　

\(A=\begin{pmatrix} \sqrt{\lambda_1}p_1 & \sqrt{\lambda_2}q_1\\ \sqrt{\lambda_1}p_2 & \sqrt{\lambda_2}q_2 \end{pmatrix}\) を 主成分負荷行列 といいます。

これは、標準化された主成分による負荷量であり、

• \(x_1\) に対する \((z_1,z_2)\) の主成分負荷量は \((\sqrt{\lambda_1}p_1,\sqrt{\lambda_2}q_1)\)
  これを \(x_1\) の主成分負荷ベクトルと呼びます。
  
• \(x_2\) に対する \((z_1,z_2)\) の主成分負荷量は \((\sqrt{\lambda_1}p_2,\sqrt{\lambda_2}q_2)\)
  これを \(x_2\) の主成分負荷ベクトルと呼びます。

これをプロットすると、影響力の大きさがベクトルの長さとして表現されます。
このため、主成分を標準化する操作がデフォルトになっています。

4.3 主成分負荷行列

＜主成分を正規化しないとき＞.

　　\(\begin{pmatrix} x_1\\ x_2 \end{pmatrix} =\begin{pmatrix} p_1 & q_1 \\ p_2 & q_2 \end{pmatrix}\begin{pmatrix} PC_1 \\ PC_2 \end{pmatrix}\)　

　　主成分負荷行列は　\(A=\begin{pmatrix} p_1 & q_1 \\ p_2 & q_2 \end{pmatrix}\)

• \(x_1\) の主成分負荷ベクトル \((p_1, q_1)\) と 、
  
• \(x_2\) の主成分負荷ベクトル \((p_2, q_2)\) は直交する単位ベクトルです。
  （ A は直交行列なので、行ベクトルも直交するため）

＜主成分を正規化するとき＞.

　 　\(\begin{pmatrix} x_1\\ x_2 \end{pmatrix} = \begin{pmatrix} \sqrt{\lambda_1}p_1 & \sqrt{\lambda_2}q_1\\ \sqrt{\lambda_1}p_2 & \sqrt{\lambda_2}q_2 \end{pmatrix} \begin{pmatrix} z_1\\ z_2 \end{pmatrix}\)　

　　主成分負荷行列は　\(A=\begin{pmatrix} \sqrt{\lambda_1}p_1 & \sqrt{\lambda_2}q_1\\ \sqrt{\lambda_1}p_2 & \sqrt{\lambda_2}q_2 \end{pmatrix}\)

　\(V(x_1)=\lambda_1p_1^2+\lambda_2q_1^2\) 、\(V(x_2)=\lambda_1p_2^2+\lambda_2q_2^2\) であるので、

• \(x_1\) の主成分負荷ベクトル \((\sqrt{\lambda_1}p_1, \sqrt{\lambda_2}q_1)\) の大きさは、\(x_1\)の標準偏差に等しく、
  
• \(x_2\) の主成分負荷ベクトル \((\sqrt{\lambda_1}p_2, \sqrt{\lambda_2}q_2)\) の大きさは、\(x_2\)の標準偏差に等しくなります。

　　（A は直交行列ではなく、\(x_1,~x_2\) の主成分負荷ベクトルは直交しません）
　

＜発展：二変量正規分布＞

固有ベクトルに固有値の平方根を乗ずる理由は、二変量正規分布の定義とつながっています。
将来の学習のため、行列を用いて簡潔に記します。

　　\(\begin{pmatrix} \sqrt{\lambda_1}p_1 & \sqrt{\lambda_2}q_1\\ \sqrt{\lambda_1}p_2 & \sqrt{\lambda_2}q_2 \end{pmatrix} = \begin{pmatrix} p_1 & q_1 \\ p_2 & q_2 \end{pmatrix} \begin{pmatrix} \sqrt{\lambda_{1}} & 0 \\ 0 & \sqrt{\lambda_{2}} \end{pmatrix}\)

であるので、

　　\(\begin{pmatrix} x_1\\ x_2 \end{pmatrix} = \begin{pmatrix} p_1 & q_1 \\ p_2 & q_2 \end{pmatrix} \begin{pmatrix} \sqrt{\lambda_{1}} & 0 \\ 0 & \sqrt{\lambda_{2}} \end{pmatrix} \begin{pmatrix} z_1\\ z_2 \end{pmatrix}\)

• 2変量標準正規分布 \((z_1,z_2)\) の円状分布を x軸方向に\(\sqrt{\lambda_1}\)倍、y方向に\(\sqrt{\lambda_2}\)倍に拡大縮小することで楕円状分布となり、
  
• 直交行列で回転することで楕円を回転した正規分布 \((x_1,x_2)\) となります。
  （固有ベクトル方向を適切にすると直交行列は回転になります）

簡潔に表現すると、

　　\(X\sim N_2(m,\Sigma)\) から構成した主成分負荷行列 \(A\) に対し、

　　　　\(X=AZ\) 、\(Z\sim N_2(0,I_2)\)　　が成立します。

\(X\) から \(Z\) を求めることは \(X\) を標準化することであり、\(Z\) から \(X\) を求めることは 標準正規分布から一般の正規分布を構成することです。
主成分座標 \((z_1,~z_2)\) は、２変量正規変量 \((x_1,~x_2)\)の標準化と言えます。

これらのことは、線形代数で二変量正規分布を学んでから、再度学習してください。

4.4 biplot

biplot は、データの主成分座標と変数の主成分負荷ベクトルを一つのグラフに重ねて表示する関数です。

biplot 関数による描画は容易ですが、図の意味を説明できるでしょうか。
ここでは、主成分座標と主成分負荷ベクトルをplot( ) 関数で描画することで理解を深めます。

設定には、prcomp のオプションと biplot のオプションがあります。

＜prcomp のオプション＞

• prcomp(data,scale=FALSE) ：dataを標準化しないで処理。
  
• prcomp(data,scale=TRUE) ：dataを標準化して処理。

＜biplot のオプション＞

• biplot(pca_result, scale=0) ：主成分座標を標準化しないで表示。
  
• biplot(pca_result, scale=1) ：主成分座標を標準化して表示。

これらの組み合わせで、４通りの出力が可能です。
それぞれの場合の特徴を把握しましょう。
一つ一つ確かめます。

5 データを標準化せず分析

• \(\Sigma=\begin{pmatrix} 1 & 1.4 \\ 1.4 & 4 \end{pmatrix}\)　である正規乱数 data
  \(s_x=1,~s_y=2,~r=0.7\) として生成します。

library(mvtnorm)
r <- 0.7
sx <- 1
sy <- 2
Sigma <- matrix(c(sx^2,r*sx*sy,r*sx*sy,sy^2), ncol=2)
Sigma

##      [,1] [,2]
## [1,]  1.0  1.4
## [2,]  1.4  4.0

データを20個生成します。
prcomp( )のオプションを “scale=FALSE” とすると、データを標準化しません。
固有ベクトルは、rotation で表示されます。

x <- rmvnorm(20, mean=c(0,0), sigma=Sigma)  
data <- data.frame(X=x[,1],Y=x[,2])
pca_result <- prcomp(data, scale = FALSE)
ve <-pca_result$rotation#固有ベクトル
ve

##         PC1        PC2
## X 0.3715088  0.9284294
## Y 0.9284294 -0.3715088

固有ベクトルは、直交する単位ベクトルです。
データは標準正規分布を固有ベクトル方向へ固有値の平方根倍に拡大縮小した楕円状の分布になります。
第１固有ベクトルが楕円の長軸方向、第２固有ベクトルが単軸方向になります。

5.1 主成分を正規化しないとき

• 左図に data をプロット（固有ベクトルは第１を青、第2を赤で表示）.
  
• 右図に主成分座標 data_PC をプロット（正規化しない）.

par(mfrow=c(1,2))
plot(data,asp=1,type="n",xlim=c(-3,3),main="data")
text(1:20,x=data[,1],y=data[,2],col="lightblue")
abline(a=0,b=ve[2,1]/ve[1,1],col="lightblue")
abline(a=0,b=ve[2,2]/ve[1,2],col="pink")
arrows(x0=0, y0=0, x1=ve[1,1],  y1=ve[2,1],col="blue",length = 0.1)
arrows(x0=0, y0=0, x1=ve[1,2],  y1=ve[2,2],col="red",length = 0.1)
plot(pca_result$x,,asp=1,type="n",xlim=c(-5,5),main="data_PC")
text(1:20,x=pca_result$x[,1],y=pca_result$x[,2],col="lightblue")
abline(h=0,col="lightblue")
abline(v=0,col="pink")
arrows(x0=0, y0=0, x1=1,  y1=0, col="blue",length = 0.1)
arrows(x0=0, y0=0, x1=0,  y1=1, col="red", length = 0.1)

• 主成分プロットは、楕円の標準形の形状に分布します。
  （但し、固有ベクトルの向きに応じ、左右反転・上下反転します）

＜主成分負荷ベクトル＞

\(X= (p_1, q_1)\)、\(Y= (p_2, q_2)\) を青矢印で記入します。
主成分負荷ベクトルは直交する単位ベクトルです。

• 同じ尺度では負荷量ベクトルは小さくなるので（左図）、
  
• 尺度を変えて拡大して表示します（右図）。

par(mfrow=c(1,2))
plot(pca_result$x,,asp=1,type="n",xlim=c(-5,5))
text(1:20,x=pca_result$x[,1],y=pca_result$x[,2],col="lightblue")
abline(h=0,col="lightblue")
abline(v=0,col="pink")
arrows(x0=0, y0=0, x1=ve[1,1],  y1=ve[1,2],col="blue",length = 0.2)
text(x = 1.5*ve[1,1], y = 1.5*ve[1,2], labels = "X", col = "blue")
arrows(x0=0, y0=0, x1=ve[2,1],  y1=ve[2,2],col="blue",length = 0.2)
text(x = 1.5*ve[2,1], y = 1.5*ve[2,2], labels = "Y", col = "blue")
plot(pca_result$x,,asp=1,type="n",xlim=c(-5,5))
text(1:20,x=pca_result$x[,1],y=pca_result$x[,2],col="lightblue")
abline(h=0,col="lightblue")
abline(v=0,col="pink")
arrows(x0=0, y0=0, x1=4*ve[1,1],  y1=4*ve[1,2],col="blue",length = 0.2)
text(x = 4.5*ve[1,1], y = 4.5*ve[1,2], labels = "X", col = "blue")
arrows(x0=0, y0=0, x1=4*ve[2,1],  y1=4*ve[2,2],col="blue",length = 0.2)
text(x = 4.5*ve[2,1], y = 4.5*ve[2,2], labels = "Y", col = "blue")

• 主成分負荷ベクトルは直交する単位ベクトルです。
  直交しているのは、\(A=\begin{pmatrix} p_1 & q_1 \\ p_2 & q_2 \end{pmatrix}\) が直交行列であり、行ベクトルも直交するからです。
  
• 主成分負荷ベクトルは \(A\) の行ベクトルですが、 \(A\)が対称行列のときは列ベクトル（dataの固有ベクトル）にも一致します。これは、固有ベクトルの向きのとりかたに依存します。

主成分負荷ベクトルの成分は、変数を構成する主成分の割合を表します。

ve

##         PC1        PC2
## X 0.3715088  0.9284294
## Y 0.9284294 -0.3715088

主成分プロットは、元の dataプロットを.

　　　\(\begin{pmatrix} PC_1\\ PC_2 \end{pmatrix} =\begin{pmatrix} p_1 & p_2 \\ q_1 & q_2 \end{pmatrix}\begin{pmatrix} x_1 \\ x_2 \end{pmatrix}\)

によって直交変換（回転と対称移動）したものです。

pca_result をもとにしたプロット（左図）と、biplot 関数の出力（右図）を比べてみましょう。

par(mfrow=c(1,2))
plot(pca_result$x,,asp=1,type="n",xlim=c(-5,5))
text(1:20,x=pca_result$x[,1],y=pca_result$x[,2],col="lightblue")
abline(h=0,col="lightblue")
abline(v=0,col="pink")
arrows(x0=0, y0=0, x1=4*ve[1,1],  y1=4*ve[1,2],col="blue",length = 0.2)
text(x = 4.5*ve[1,1], y = 4.5*ve[1,2], labels = "X", col = "blue")
arrows(x0=0, y0=0, x1=4*ve[2,1],  y1=4*ve[2,2],col="blue",length = 0.2)
text(x = 4.5*ve[2,1], y = 4.5*ve[2,2], labels = "Y", col = "blue")
biplot(pca_result,col = c("lightblue", "blue"),scale = 0) 
abline(h=0,col="lightblue")
abline(v=0,col="pink")

biplot では、主成分負荷ベクトルの目盛りは上・右にあり、見やすいように自動調整されています。

＜研究＞

主成分負荷ベクトルと主成分プロットの関係を考察します。
\(\begin{pmatrix} PC_1\\ PC_2 \end{pmatrix} =\begin{pmatrix} p_1 & p_2 \\ q_1 & q_2 \end{pmatrix}\begin{pmatrix} x_1 \\ x_2 \end{pmatrix}\)　より、

　　　\(\begin{pmatrix} PC_1\\ PC_2 \end{pmatrix} = x_1\begin{pmatrix} p_1\\ q_1 \end{pmatrix} + x_2\begin{pmatrix} p_2\\ q_2 \end{pmatrix}\)

\((PC_1,~PC_2) =k(p_1,~q_1)\)　のとき、

　　　\(k\begin{pmatrix} p_1\\ q_1 \end{pmatrix} = x_1\begin{pmatrix} p_1\\ q_1 \end{pmatrix} + x_2\begin{pmatrix} p_2\\ q_2 \end{pmatrix}\)

　　\(\therefore~(x_1,x_2)=(k,0)\).

即ち、\(x_1\) の主成分負荷ベクトルの延長上に主成分点があるとき、データは \(x_1\) にのみ依存し \(x_2\) は0に近くなります。

これは、「主成分負荷ベクトルを基底としたとき、主成分点は座標 \((x_1,x_2)\) の位置にある」ということが背景にあります。

5.2 主成分を正規化するとき

＜主成分座標を正規化すると円状の分布になります＞

• \(\begin{pmatrix} z_1 \\ z_2 \end{pmatrix} = \begin{pmatrix} 1/\sqrt{\lambda_{1}} & 0 \\ 0 & 1/\sqrt{\lambda_{2}} \end{pmatrix} \begin{pmatrix} PC_1\\ PC_2 \end{pmatrix}\)

この主成分座標は、正規化しない場合の主成分座標の固有値平方根逆数倍です。

固有値の平方根 sdev を用いて、標準化された主成分座標を計算すると、

sdev <- pca_result$sdev
norm_result <- pca_result$x %*% diag(1/sdev)
head(norm_result)

##             [,1]        [,2]
## [1,]  1.11295726 -1.43513218
## [2,] -0.03036659  1.69621260
## [3,]  0.44371026 -0.07494734
## [4,] -0.27747317 -0.05252919
## [5,] -0.94384836  1.13719407
## [6,] -0.01544462 -1.04225550

＜主成分負荷行列は、rotationに対角行列を乗じて得られます＞

　　\(\begin{pmatrix} \sqrt{\lambda_1}p_1 & \sqrt{\lambda_2}q_1\\ \sqrt{\lambda_1}p_2 & \sqrt{\lambda_2}q_2 \end{pmatrix} = \begin{pmatrix} p_1 & q_1 \\ p_2 & q_2 \end{pmatrix} \begin{pmatrix} \sqrt{\lambda_{1}} & 0 \\ 0 & \sqrt{\lambda_{2}} \end{pmatrix}\)

転置して列を取り出すと、

• \(\begin{pmatrix} \sqrt{\lambda_1}p_1\\ \sqrt{\lambda_2}q_1 \end{pmatrix} = \begin{pmatrix} \sqrt{\lambda_{1}} & 0 \\ 0 & \sqrt{\lambda_{2}} \end{pmatrix} \begin{pmatrix} p_1\\ q_1 \end{pmatrix}\)　

• \(\begin{pmatrix} \sqrt{\lambda_1}p_1\\ \sqrt{\lambda_2}q_2 \end{pmatrix} = \begin{pmatrix} \sqrt{\lambda_{1}} & 0 \\ 0 & \sqrt{\lambda_{2}} \end{pmatrix} \begin{pmatrix} p_2 \\ q_2 \end{pmatrix}\)　

この主成分負荷ベクトルは、正規化しない場合の負荷ベクトルの固有値平方根倍です。

＜主成分負荷行列＞

sdev <- pca_result$sdev
vn <- ve %*% diag(sdev);vn

##        [,1]       [,2]
## X 0.8809357  0.5810462
## Y 2.2015266 -0.2325043

主成分負荷ベクトルは、主成分負荷行列の行ベクトルです。

＜正規化しない場合のbiplotと正規化する場合のbiplotの比較＞

　主成分座標を正規化すると、

• 主成分座標は固有値平方根の逆数倍、
  
• 主成分負荷ベクトルは固有値の平方根倍.

主成分座標の伸縮倍率と負荷ベクトルの伸縮倍率は、逆数の関係にあります。

sdev

## [1] 2.3712374 0.6258378

これが固有値の平方根 \(sdev=(\sqrt{\lambda_{1}},\sqrt{\lambda_{2}})\) です。
左図から右図への変化を観察し、理解を深めましょう。

par(family = "HiraKakuProN-W3") 
par(mfrow=c(1,2))
plot(pca_result$x,,asp=1,type="n",xlim=c(-5,5),main="正規化しないbiplot")
text(1:20,x=pca_result$x[,1],y=pca_result$x[,2],col="lightblue")
abline(h=0,col="lightblue")
abline(v=0,col="pink")
arrows(x0=0, y0=0, x1=1*ve[1,1],  y1=1*ve[1,2],col="blue",length = 0.2)
text(x = 1.5*ve[1,1], y = 1.5*ve[1,2], labels = "X", col = "blue")
arrows(x0=0, y0=0, x1=1*ve[2,1],  y1=1*ve[2,2],col="blue",length = 0.2)
text(x = 1.5*ve[2,1], y = 1.5*ve[2,2], labels = "Y", col = "blue")
plot(norm_result,asp=1,type="n",xlim=c(-4,4),main="正規化するbiplot",xlab="PC1",ylab="PC2")
text(1:20,x=norm_result[,1],y=norm_result[,2],col="lightblue")
abline(h=0,col="lightblue")
abline(v=0,col="pink")
arrows(x0=0, y0=0, x1=vn[1,1],  y1=vn[1,2],col="blue",length = 0.2)
text(x = 1.5*vn[1,1], y = 1.5*vn[1,2], labels = "X", col = "blue")
arrows(x0=0, y0=0, x1=vn[2,1],  y1=vn[2,2],col="blue",length = 0.2)
text(x = 1.5*vn[2,1], y = 1.5*vn[2,2], labels = "Y", col = "blue")

• 主成分は横方向へ\(1/\sqrt{\lambda_{1}}\)倍に縮小されるので、主成分負荷ベクトルのx座標は\(\sqrt{\lambda_{1}}\)倍に拡大.
  

• 縦方向へ\(1/\sqrt{\lambda_{2}}\)倍に拡大されるので、主成分負荷ベクトルのy座標は\(\sqrt{\lambda_{2}}\)倍に縮小.

• 左図で直交していた主成分負荷ベクトルは、右図では直交しません。

• 左図では主成分負荷ベクトルの大きさは１ですが、右図では各変量の標準偏差になります。
  右図での主成分負荷ベクトルの大きさは、

　　　　　\(SD(x_1)=\sqrt{\lambda_1p_1^2+\lambda_2q_1^2},~~SD(x_2)=\sqrt{\lambda_1p_2^2+\lambda_2q_2^2}\)

その値を計算すると、

c(sqrt(sum(vn[1,]^2)),sqrt(sum(vn[2,]^2)))

## [1] 1.055302 2.213770

関係式は、

　　　\(\begin{pmatrix} x_1\\ x_2 \end{pmatrix} = \begin{pmatrix} \sqrt{\lambda_1}p_1 & \sqrt{\lambda_2}q_1\\ \sqrt{\lambda_1}p_2 & \sqrt{\lambda_2}q_2 \end{pmatrix} \begin{pmatrix} z_1\\ z_2 \end{pmatrix}\)
　
係数行列は直交行列でないので逆行列は転置だけでは求まらず、

　　　\(\begin{pmatrix} z_1 \\ z_2 \end{pmatrix} = \begin{pmatrix} 1/\sqrt{\lambda_{1}} & 0 \\ 0 & 1/\sqrt{\lambda_{2}} \end{pmatrix} \begin{pmatrix} p_1 & p_2 \\ q_1 & q_2 \end{pmatrix}\begin{pmatrix} x_1 \\ x_2 \end{pmatrix}\)

正規化しない場合に成立していた主成分プロットと主成分負荷ベクトルの関係は成り立たなくなります。

＜biplotの示す情報＞

　dataプロット（左図）とbiplot（右図）を観察しましょう。
　biplotからどんな情報が得られるでしょうか？

par(family = "HiraKakuProN-W3") 
par(mfrow=c(1,2))
plot(data,asp=1,type="n",xlim=c(-3,3),main="data")
text(1:20,x=data[,1],y=data[,2],col="lightblue")
abline(a=0,b=ve[2,1]/ve[1,1],col="lightblue")
abline(a=0,b=ve[2,2]/ve[1,2],col="pink")
arrows(x0=0, y0=0, x1=ve[1,1],  y1=ve[2,1],col="blue",length = 0.1)
arrows(x0=0, y0=0, x1=ve[1,2],  y1=ve[2,2],col="red",length = 0.1)
biplot(pca_result,col = c("lightblue", "blue"),scale = 1,main="biplot") 

• 主成分プロットは円状に分布しますが、点の相対的位置関係は変わりません。
  （標準化されているので、ユークリッド距離とは一致しません）

• ベクトルの長さは 変量\(x_1,x_2\) の標準偏差を表し、
  

• ベクトルのなす角の余弦は 変量\(x_1,x_2\) の相関係数に一致します。

数理的説明の概略を記します。

　\(\begin{pmatrix} x_1\\ x_2 \end{pmatrix} = \begin{pmatrix} \sqrt{\lambda_1}p_1 & \sqrt{\lambda_2}q_1\\ \sqrt{\lambda_1}p_2 & \sqrt{\lambda_2}q_2 \end{pmatrix} \begin{pmatrix} z_1\\ z_2 \end{pmatrix}\) より、

• \(x_1,x_2\)の分散は、
  \(V(x_1)=\lambda_1p_1^2+\lambda_2q_1^2\)
  \(V(x_2)=\lambda_1p_2^2+\lambda_2q_2^2\)

• \(x_1,x_2\)の共分散は、
  \(Cov(x_1,x_2)=\lambda_1p_1p_2+\lambda_2q_1q_2\)

• \(x_1,x_2\)の相関係数は、
  \(cor(x_1,x_2)=\dfrac{\lambda_1p_1p_2+\lambda_2q_1q_2}{{\sqrt{\lambda_1p_1^2+\lambda_2q_1^2}}{\sqrt{\lambda_1p_2^2+\lambda_2q_2^2}}}\)

これは、主成分負荷ベクトル \((\sqrt{\lambda_1}p_1, \sqrt{\lambda_2}q_1)\)　\((\sqrt{\lambda_1}p_2, \sqrt{\lambda_2}q_2)\) のなす角の余弦に一致します。

これらは、数理の学習が進んでから、再度学習してください。

6 データを標準化して分析

標準化されたデータでは、分散共分散行列は相関行列に一致します。
そこで、相関行列を分散共分散行列としてデータを生成します。

• \(s_x=1,~s_y=1,~r=0.7\) である二変量正規乱数 data

「3.3.2 小標本データ」のように、20個のデータを生成します。
得られたデータを scale( ) 関数で完全に標準化します。

分散共分散行列は、

Sigma <- matrix(c(1,0.7,0.7,1), ncol=2)
x <- rmvnorm(20, mean=c(0,0), sigma=Sigma)  
data <- data.frame(X=x[,1],Y=x[,2])
data <- scale(data)
var(data)

##           X         Y
## X 1.0000000 0.9068941
## Y 0.9068941 1.0000000

標準化されても、相関関係は残ります。

標本相関行列の固有ベクトルが主成分方向であり、主成分方向に固有値の平方根倍に伸縮された分布となります。

prcomp( ) 関数において、標準化して分析するオプションは “scale=TRUE” です。
すでに標準化していますが、 “scale=TRUE” とします。

pca_result <- prcomp(data, scale = TRUE)
pca_result

## Standard deviations (1, .., p=2):
## [1] 1.3809033 0.3051327
## 
## Rotation (n x k) = (2 x 2):
##         PC1        PC2
## X 0.7071068 -0.7071068
## Y 0.7071068  0.7071068

• Rotation の列ベクトルが主成分方向であり、
  
• Standard deviations は主成分方向の標準偏差です。

この情報から、標準正規分布を「どの方向にどれだけ伸縮した分布であるか」イメージできます。

実際に、標準化された data（左図）と、主成分を横軸・縦軸にした data_PC（右図）を描画します。

ve <-pca_result$rotation
data_PC <- pca_result$x
par(mfrow=c(1,2))
plot(data,asp=1,col="lightblue",main="data",xlim=c(-4,4))
arrows(x0=0, y0=0, x1=ve[1,1],  y1=ve[2,1], code=2, col="blue", length = 0.1)
arrows(x0=0, y0=0, x1=ve[1,2],  y1=ve[2,2], code=2, col="red", length = 0.1)
abline(a=0,b=ve[2,1]/ve[1,1],col="lightblue")
abline(a=0,b=ve[2,2]/ve[1,2],col="pink")
plot(data_PC,asp=1,col="lightblue",xlim=c(-4,4),main="data_PC")
abline(h=0,col="lightblue")
abline(v=0,col="pink")
arrows(x0=0, y0=0, x1=1,  y1=0, col="blue",length = 0.1)
arrows(x0=0, y0=0, x1=0,  y1=1, col="red", length = 0.1)

また、標準化した分析では、主成分方向は必ず y=x, y=-x となることが分かります。

このことについて補足しておきます。

6.1 相関行列の固有値・固有ベクトル

対称行列タイプ１　\(A= \begin{pmatrix} a & b \\ b & a \end{pmatrix}\)　を考えます。

このタイプの対称行列は、

• 固有値　 \(\lambda_1=a+b\)　　\(\lambda_2=a-b\).
  
• 固有ベクトル　\(\begin{pmatrix} 1 \\ 1 \end{pmatrix}\)　　\(\begin{pmatrix} -1 \\ 1 \end{pmatrix}\)

固有値と固有ベクトルは簡明です。
（拙著「高校生の線形代数 3.5対称行列」参照）

相関行列　\(A= \begin{pmatrix} 1 & r \\ r & 1 \end{pmatrix}\)　の場合、

• 固有値　 \(\lambda_1=1+r\)　　\(\lambda_2=1-r\).
  
• 固有ベクトル　\(\begin{pmatrix} 1 \\ 1 \end{pmatrix}\)　　\(\begin{pmatrix} -1 \\ 1 \end{pmatrix}\)

不動直線は、y=x と y=-x であり、これが主成分軸です。

　　rotationは　\(P=\begin{pmatrix} 1/\sqrt{2} & -1/\sqrt{2} \\ 1/\sqrt{2} & 1/\sqrt{2} \end{pmatrix}\)
列ベクトルが固有ベクトルであり主成分方向です。

＜主成分を正規化しない場合＞

　　主成分負荷行列は　\(P^{-1}=\begin{pmatrix} 1/\sqrt{2} & 1/\sqrt{2} \\ -1/\sqrt{2} & 1/\sqrt{2} \end{pmatrix}\)

　　主成分負荷ベクトル \(\begin{pmatrix} 1/\sqrt{2}\\ -1/\sqrt{2} \end{pmatrix}\)　\(\begin{pmatrix} 1/\sqrt{2}\\ 1/\sqrt{2} \end{pmatrix}\) は直交します。

　主成分負荷ベクトルはいつも同じということになります。
　これでは分析の意味をなしません。

＜主成分を正規化する場合＞

　\(\begin{pmatrix} \sqrt{\lambda_1}p_1 & \sqrt{\lambda_2}q_1\\ \sqrt{\lambda_1}p_2 & \sqrt{\lambda_2}q_2 \end{pmatrix} = \begin{pmatrix} p_1 & q_1 \\ p_2 & q_2 \end{pmatrix} \begin{pmatrix} \sqrt{\lambda_{1}} & 0 \\ 0 & \sqrt{\lambda_{2}} \end{pmatrix}\) から、

主成分負荷量ベクトルは、

• \(\begin{pmatrix} \sqrt{\lambda_1}p_1\\ \sqrt{\lambda_2}q_1 \end{pmatrix} = \begin{pmatrix} \sqrt{\lambda_{1}} & 0 \\ 0 & \sqrt{\lambda_{2}} \end{pmatrix} \begin{pmatrix} p_1\\ q_1 \end{pmatrix} =\dfrac{1}{\sqrt{2}} \begin{pmatrix} \sqrt{1+r}\\ -\sqrt{1-r} \end{pmatrix}\)　

• \(\begin{pmatrix} \sqrt{\lambda_1}p_1\\ \sqrt{\lambda_2}q_2 \end{pmatrix} = \begin{pmatrix} \sqrt{\lambda_{1}} & 0 \\ 0 & \sqrt{\lambda_{2}} \end{pmatrix} \begin{pmatrix} p_2 \\ q_2 \end{pmatrix} = \dfrac{1}{\sqrt{2}} \begin{pmatrix} \sqrt{1+r}\\ \sqrt{1-r} \end{pmatrix}\)　

その内積は \(r\) 、大きさは１なので、\(\cos\theta=\dfrac{r}{1\cdot 1}=r\)

主成分負荷量ベクトルのなす角を \(\theta\) とすると、

　　　\(\cos\theta=r=0.7\)

なす角の余弦は、相関係数と等しくなります。 biplot で変量の負荷ベクトルのなす角の余弦を見れば、変数同士の相関係数の程度がわかることになります。

angle_deg <- acos(0.7) * (180 / pi)
angle_deg

## [1] 45.573

主成分負荷ベクトルのなす角度は、約 45度となります。

6.2 主成分を正規化しないとき

biplot で主成分を正規化しないオプションは “scale=0” です。

par(mfrow=c(1,2))
plot(data,asp=1,type="n",xlim=c(-3,3),main="data")
text(1:20,x=data[,1],y=data[,2],col="lightblue")
abline(a=0,b=ve[2,1]/ve[1,1],col="lightblue")
abline(a=0,b=ve[2,2]/ve[1,2],col="pink")
arrows(x0=0, y0=0, x1=ve[1,1],  y1=ve[2,1],col="blue",length = 0.1)
arrows(x0=0, y0=0, x1=ve[1,2],  y1=ve[2,2],col="red",length = 0.1)
biplot(pca_result,col = c("lightblue", "blue"),scale = 0,main="biplot") 
abline(h=0,col="lightblue")
abline(v=0,col="pink")

biplot の主成分負荷ベクトルは、相関係数に依存せず、必ずこの状態になります。

そこで次に、主成分を正規化する場合を考えます。

6.3 主成分を正規化するとき

biplot で主成分を正規化するオプションは “scale=1” です。

scale=0 の場合を左図に、scale=1 の場合を右図に描画します。

par(mfrow=c(1,2))
plot(data,asp=1,type="n",xlim=c(-3,3),main="data")
text(1:20,x=data[,1],y=data[,2],col="lightblue")
abline(a=0,b=ve[2,1]/ve[1,1],col="lightblue")
abline(a=0,b=ve[2,2]/ve[1,2],col="pink")
arrows(x0=0, y0=0, x1=ve[1,1],  y1=ve[2,1],col="blue",length = 0.1)
arrows(x0=0, y0=0, x1=ve[1,2],  y1=ve[2,2],col="red",length = 0.1)
biplot(pca_result,col = c("lightblue", "blue"),scale = 1,main="biplot") 
abline(h=0,col="lightblue")
abline(v=0,col="pink")

• scale=1 のbiplotでは、データの分布は円状になります。
  
• 主成分負荷ベクトルの大きさは、左右どちらも１です。
  
• scale=1 のbiplotでは、主成分負荷ベクトルのなす角は約45度で、その余弦は変量の相関係数を反映しています。

7 標準化のまとめ.

7.1 データの標準化

　主成分分析ではデータを標準化することが推奨され、特に、単位の異なる変量を扱う場合は標準化が必須とされています。

標準偏差 \(s_x,~s_y\)、相関係数 \(r\) が定まれば、分散共分散行列は.

　　　\(\begin{pmatrix} s_x^2 & s_{xy}\\ s_{xy} & s_y^2 \end{pmatrix} = \begin{pmatrix} s_x^2 & r s_x s_y \\ r s_x s_y & s_y^2 \end{pmatrix}\)　によって求まります。

相関係数を \(r=0.7\) と定め、データを標準化しない場合と標準化する場合を比較しましょう。

• \(s_x=1,~s_y=2,~r=0.7\) である二変量正規乱数 data について.

r <- 0.7;sx <- 1;sy <- 2
Sigma <- matrix(c(sx^2,r*sx*sy,r*sx*sy,sy^2), ncol=2)
x <- rmvnorm(1000, mean=c(0,0), sigma=Sigma)  
data <- data.frame(x1=x[,1],x2=x[,2])
pca_result <- prcomp(data, scale = FALSE)
ve <-pca_result$rotation#固有ベクトル
cor(data$x1,data$x2)#相関係数

## [1] 0.722464

• data を標準化した data1 について.

x <- scale(x)#標準化
data1 <- data.frame(x1=x[,1],x2=x[,2])
pca_result1 <- prcomp(data1, scale = FALSE)
ve1 <-pca_result1$rotation#固有ベクトル
cor(data1$x1,data1$x2)#相関係数

## [1] 0.722464

標準化しても、相関関係が変わるわけではなく、どちらも相関係数は同じです。

左図に標準化しないdata、右図に標準化した data1 をプロットします。

par(mfrow=c(1,2))
plot(data,asp=1,col="lightblue",main="data",xlim=c(-4,4))
abline(a=0,b=ve[2,1]/ve[1,1],col="lightblue")
abline(a=0,b=ve[2,2]/ve[1,2],col="pink")
arrows(x0=0, y0=0, x1=ve[1,1],  y1=ve[2,1],col="blue",length = 0.1)
arrows(x0=0, y0=0, x1=ve[1,2],  y1=ve[2,2],col="red",length = 0.1)
plot(data1,asp=1,col="lightblue",main="data1",xlim=c(-4,4))
abline(a=0,b=ve1[2,1]/ve1[1,1],col="lightblue")
abline(a=0,b=ve1[2,2]/ve1[1,2],col="pink")
arrows(x0=0, y0=0, x1=ve1[1,1],  y1=ve1[2,1],col="blue",length = 0.1)
arrows(x0=0, y0=0, x1=ve1[1,2],  y1=ve1[2,2],col="red",length = 0.1)

• 分散が変わると、固有ベクトルの方向が変わり、主成分方向も変わります。
  
• 標準化した右図では、主成分方向は常に y=x, y=-x です。

7.2 主成分の標準化

　biplot のデフォルトでは主成分は標準化され、主成分プロットは円状に分布します。第１主成分の分散が最大であることが分かるよう表示するには、scale=o としなければなりません。

• 主成分の標準化は何のために行うのでしょうか。

　　　\(\begin{pmatrix} x_1\\ x_2 \end{pmatrix} = \begin{pmatrix} \sqrt{\lambda_1}p_1 & \sqrt{\lambda_2}q_1\\ \sqrt{\lambda_1}p_2 & \sqrt{\lambda_2}q_2 \end{pmatrix} \begin{pmatrix} z_1\\ z_2 \end{pmatrix}\)　‥(1)

主成分を正規化し、主成分の影響力を負荷行列（loading）で表現しています。
これにより、負荷行列（loading）を見れば、主成分が変数にどう影響しているかが分かります。

係数行列を主成分負荷行列といいます。

この主成分負荷行列について、

• 列ベクトルは、その主成分の各変数への影響力を表し、
  
• 行ベクトルは、その変数における各主成分の割合を表します。

主成分負荷行列の成分の平方和（sum of squares）については、

• 行ベクトルの平方和は各変量の分散に等しく、

　　　\(V(x_1)=\lambda_1p_1^2+\lambda_2q_1^2\)、　\(V(x_2)=\lambda_1p_2^2+\lambda_2q_2^2\).

• 列ベクトルの平方和（ss loadings）は、各主成分の分散に等しくなります。

　　　\((\sqrt{\lambda_1}p_1)^2+(\sqrt{\lambda_1}p_2)^2=\lambda_1\)、　\((\sqrt{\lambda_2}q_1)^2+(\sqrt{\lambda_2}q_2)^2=\lambda_2\).

7.3 データと主成分の標準化

データと主成分を標準化した分析がよく用いられるのは、主成分負荷行列がデータと主成分の相関係数になるからです。

行列を用いて説明します。高校生の方は行列を学んでから、再度読み返してみると良いでしょう。
説明には、標準化した変量の共分散が相関係数に一致することを用います。

変量 \(x,~y\) の標準偏差 \(s_x=1,~s_y=1\) とすると、

　　相関係数　\(Cor(x,y)=\dfrac{Cov(x,y)}{s_x s_y}=E(xy)\).
　（変量の積の期待値が相関係数になります）

データ \((x_1,x_2)\) と主成分 \((z_1,z_2)\) が標準化されているとします。

＜データと主成分の相関＞.
行列を用いると、一括して相関係数が求まります。

7.2 の式 (1) の両辺に、右から \(\begin{pmatrix} z_1 & z_2 \end{pmatrix}\) をかけます。

\(\begin{pmatrix} x_1\\ x_2 \end{pmatrix} \begin{pmatrix} z_1 & z_2 \end{pmatrix} = \begin{pmatrix} \sqrt{\lambda_1}p_1 & \sqrt{\lambda_2}q_1\\ \sqrt{\lambda_1}p_2 & \sqrt{\lambda_2}q_2 \end{pmatrix} \begin{pmatrix} z_1\\ z_2 \end{pmatrix} \begin{pmatrix} z_1 & z_2 \end{pmatrix}\) から、

　　\(\begin{pmatrix} x_1z_1 & x_1z_2 \\ x_2z_1 & x_2z_2 \end{pmatrix} = \begin{pmatrix} \sqrt{\lambda_1}p_1 & \sqrt{\lambda_2}q_1\\ \sqrt{\lambda_1}p_2 & \sqrt{\lambda_2}q_2 \end{pmatrix} \begin{pmatrix} z_1z_1 & z_1z_2 \\ z_2z_1 & z_2z_2 \end{pmatrix}\)

両辺の期待値をとります。

　\(E(z_1^2)=E(z_2^2)=1\) 、 \(E(z_1z_2)=E(z_2z_1)=0\) であるので、

　　\(\begin{pmatrix} Cor(x_1,z_1) & Cor(x_1,z_2) \\ Cor(x_2,z_1) & Cor(x_2,z_2) \end{pmatrix} = \begin{pmatrix} \sqrt{\lambda_1}p_1 & \sqrt{\lambda_2}q_1\\ \sqrt{\lambda_1}p_2 & \sqrt{\lambda_2}q_2 \end{pmatrix}\)

主成分負荷行列の成分は相関係数となります。

	\(z_1\)	\(z_2\)	平方和
\(x_1\)	\(Cor(x_1,z_1)\)	\(Cor(x_1,z_2)\)	1
\(x_2\)	\(Cor(x_2,z_1)\)	\(Cor(x_2,z_2)\)	1

biplotのデフォルトで主成分が標準化されているのは、変量と主成分との相関関係を見るためです。
biplot は変数間の相関の視覚的な解釈に適したプロットなのです。

biplot の特徴としては、

• \(x_1,x_2\) を表すベクトルの大きさは両方とも１です。
  
• ベクトルの成分が、各主成分との相関係数を表します。
  通常は第１主成分の影響が大きいので、傾きの絶対値は小さくなります。
  
• ベクトルのなす角の余弦が \(x_1,x_2\) の相関係数を表します。

このような理由で、データと主成分の両方を標準化することが、biplot のデフォルトになっています。

8 分析と評価

主成分分析の目的は、データを要約し特徴を見出すことです。
主成分はデータの個性をまとめた新しい指標であり、データの理解に有用な指標である必要があります。

そのためには、

• 主成分にまとめる際に、失われる情報を少なくすること.
  
• 主成分がデータをよく説明していること.

が重要になります。

8.1 寄与率

主成分の分散がデータの全分散に占める割合を寄与率といいます。
寄与率が大きい主成分は多くの情報を持っています。

＜行列のトレース＞.
　線形代数では、行列 A の対角成分の合計をトレースと呼び \(tr(A)\) と表します。

• データの分散共分散行列 \(\Sigma\) の対角成分は各変量の分散なので、\(tr( \Sigma)\) は全分散の合計になります。
  
• 主成分の寄与率は、\(tr(\Sigma)\)に対するその主成分の分散の割合です。

特に、データと主成分を標準化した分析では、

• 各変量の分散は１であるので、全分散の合計は変量の個数に一致します。
  （下の表で各行の平方和が分散で、その合計は２です）.
  
• 各列の平方和は ss loading と呼ばれ、各主成分の分散になります。
  （下の表では \(\lambda_1,\lambda_2\) であり、\(\lambda_1+\lambda_2=2\)）.

	\(z_1\)	\(z_2\)	平方和
\(x_1\)	\(\sqrt{\lambda_1}p_1\)	\(\sqrt{\lambda_2}q_1\)	1
\(x_2\)	\(\sqrt{\lambda_1}p_2\)	\(\sqrt{\lambda_2}q_2\)	1
平方和	\(\lambda_1\)	\(\lambda_2\)

• \(z_1\) の寄与率は \(\dfrac{\lambda_1}{\lambda_1+\lambda_2}=\dfrac{\lambda_1}{2}\)　　\(z_2\) の寄与率は \(\dfrac{\lambda_2}{\lambda_1+\lambda_2}=\dfrac{\lambda_2}{2}\).
  
• \(\lambda_1>\lambda_2\)であるので、\(z_1\)の寄与率＞\(z_2\)の寄与率.

主成分の分散は降順となっており、データの次元圧縮を合理的に行うことができます。

8.2 単純構造

「主成分がデータをよく説明している」とはどのような状態でしょうか。

データ \((x_1,x_2)\) から主成分 \((z_1,z_2)\) が得られたとします。

• \(x_1=a\cdot z_1\)
  \(x_2=d\cdot z_2\)

このように、変量が一つの主成分に支配される構造を単純構造と言います。

実際、\(z_1\)が決まれば \(x_1\)が決まり、\(z_2\)が決まれば \(x_2\)が決まるので、 主成分はデータをよく説明しています。

しかし、実際のデータには相関があるので、完全な単純構造にはなりません。

• \(x_1=a\cdot z_1+b\cdot z_2~~~~(b\fallingdotseq 0)\)
  \(x_2=c\cdot z_1+d\cdot z_2~~~~(c\fallingdotseq 0)\)

変量が近似的に一つの主成分に支配されていれば単純構造とみなします。

＜単純構造＞

• 単純構造とは、各変数が特定の主成分に対してのみ高い負荷量を持ち、他の主成分への負荷量がゼロに近い状態をいいます。

• biplotでのイメージでは、変量の主成分負荷ベクトルが横軸・縦軸に近い状態のことです。

＜回転＞

• 主成分を単純構造に近づけるには回転を行います。
  
• 主成分主成分が標準化され円状に分布していれば、回転してもデータ分布に影響を与える心配はありません。

psychパッケージの principal( ) 関数には、“rotate” オプションがあり、

• rotate = “none” とすると回転を行わなわず、
  rotate = “varimax” と指定すると単純構造に近づける回転を適用します。

“varimax” は回転の基準のひとつで、各主成分ごとの係数二乗値の分散の合計を最大化します。

8.3 分析例

次の例は、principal( ) 関数による分析例です。

• \(\Sigma=\begin{pmatrix} 1 & 1 \\ 1 & 4 \end{pmatrix}\)　である正規乱数 data を
  \(s_x=1,~s_y=2,~r=0.5\) として生成し、

回転を行わない biplot（左）、回転を行う biplot（右）として表示します。

library(psych)
r <- 0.5;sx <- 1;sy <- 2
Sigma <- matrix(c(sx^2,r*sx*sy,r*sx*sy,sy^2), ncol=2)
set.seed(11145378)
x <- rmvnorm(20, mean=c(0,0), sigma=Sigma)  
data <- data.frame(X=x[,1],Y=x[,2])
pca_result <- principal(data, nfactors = 2, rotate = "none")
vpca_result <- principal(data, nfactors = 2, rotate = "varimax")
par(mfrow=c(1,2))
biplot(pca_result,col = c("lightblue", "blue"),main="Rotate : none") 
biplot(vpca_result,col = c("lightblue", "blue"),main="Rotate : varimax") 

回転なので、主成分プロットと各変量の主成分負荷ベクトルの相対的位置関係は変わりません。

• 右図では、RC2はYを、RC1はXをよく説明していますが、
  左図では、PC1,PC2 の変量 X,Y への影響が中途半端です。

バリマックス回転各主成分ごとの係数二乗値の分散の合計を最大化する回転です。
つまり、

• 単純構造のときの係数二乗値の分散の合計は、そうでないときより大きくなります。

このことを確かめましょう。

二つの数\(A,B\) の分散は、
　　　\(\left(A-\dfrac{A+B}{2}\right)^2+\left(B-\dfrac{A+B}{2}\right)^2=\dfrac{(A-B)^2}{2}\)　で計算されます。

［右図］

• \(X=0.96\cdot RC2+0.27\cdot RC1\)
  \(Y=0.31 \cdot RC2+0.95\cdot RC1\)

RC は、回転された主成分（Rotated Components）を意味します。

片方の係数がゼロに近ければ、各主成分の係数の分散は大きくなります。

• \((0.96,0.31)\)の二乗値の分散は、\((0.96^2-0.31^2)^2/2=0.3407\).
  \((0.27,0.95)\)の二乗値の分散は、\((0.27^2-0.95^2)^2/2=0.3411\).

分散合計は、\(0.6848\)

［左図］

• \(X=0.88 \cdot PC1-0.47 \cdot PC2\)
  \(Y=0.88\cdot PC1+0.47 \cdot PC2\)

係数が似ていれば、各主成分の係数の分散は小さくなります。

• \((0.88,0.88)\)の二乗値の分散は、\((0.88^2-0.88^2)^2/2=0\).
  \((-0.47,0.47)\)の二乗値の分散は、\((0.47^2-0.47^2)^2/2=0\).

分散合計は、\(0\).

右図の分散合計のほうが、左図よりも大きくなっています。

principal( ) 関数の詳細については、因子分析を学ぶときに説明します。

9 分析の出力

9.1 principal( )の出力

principal( )は因子分析のための関数であり、その出力には因子分析の指標が含まれています。
因子分析は、変量を潜在因子と独自因子で表す手法で、社会科学の分野でよく用いられます。
principal( ) は、主成分を因子とみなして分析しますが、独自因子（誤差）はゼロとして扱います。

上記「8.3 分析例」を題材に、Rの出力結果を説明します。

＜回転しない場合＞.
「回転を行わない biplot（左）」の分析出力です。

• “Standardized loadings” は主成分負荷量が標準化されていることを意味します。
  “based upon correlation matrix” はデータが標準化され相関行列を用いた分析であることを意味します。

主成分負荷行列のところの h2,u2,com は因子分析の指標です。

• 因子分析では、変量＝主成分＋誤差　であり、
  主成分の分散：h2（共通性）、誤差分散：u2（独自性）
  主成分分析では、誤差を想定しないので、h2=1,u2=0 です。
  （詳細には立ち入りません）

• com はComplexity（複雑性）であり、“Mean item complexity” の数値と一致します。 これは影響を受ける因子の数に対応します。単純構造であれば1に近く、単純構造でなければ大きくなります。

ss loadinds（sum of squared loadongs）は、主成分負荷行列の列ベクトル成分の平方和であり、各主成分の分散を表します。

• Proportion Var（寄与率）は、主成分の分散の全分散に対する割合、
  
• Cumulative Var（累積寄与率）は、第1主成分から順に寄与率を合計したものです。

pca_result

## Principal Components Analysis
## Call: principal(r = data, nfactors = 2, rotate = "none")
## Standardized loadings (pattern matrix) based upon correlation matrix
##    PC1   PC2 h2       u2 com
## X 0.88 -0.47  1 -2.2e-16 1.5
## Y 0.88  0.47  1 -1.1e-16 1.5
## 
##                        PC1  PC2
## SS loadings           1.55 0.45
## Proportion Var        0.78 0.22
## Cumulative Var        0.78 1.00
## Proportion Explained  0.78 0.22
## Cumulative Proportion 0.78 1.00
## 
## Mean item complexity =  1.5
## Test of the hypothesis that 2 components are sufficient.
## 
## The root mean square of the residuals (RMSR) is  0 
##  with the empirical chi square  0  with prob <  NA 
## 
## Fit based upon off diagonal values = 1

“2 components are sufficient” の部分は、データが２つの主成分で十分説明されているかを検定評価しています。この部分は因子分析の学習の際に学んでください。
　

＜回転する場合＞.

「回転を行う biplot（右）」の分析出力です。

バリマックス回転により単純構造に近づけているので、

• com（複雑性）は1に近くなります。
  これは影響を受ける因子の数と考えることができます。
  
• SS loadings は平均化され、固有値とも一致しなくなります。
  第１主成分に集中していた説明力は平均化され、各々の主成分が各々の変量を説明するようになります。

vpca_result

## Principal Components Analysis
## Call: principal(r = data, nfactors = 2, rotate = "varimax")
## Standardized loadings (pattern matrix) based upon correlation matrix
##    RC2  RC1 h2       u2 com
## X 0.96 0.27  1 -2.2e-16 1.2
## Y 0.31 0.95  1 -1.1e-16 1.2
## 
##                        RC2  RC1
## SS loadings           1.02 0.98
## Proportion Var        0.51 0.49
## Cumulative Var        0.51 1.00
## Proportion Explained  0.51 0.49
## Cumulative Proportion 0.51 1.00
## 
## Mean item complexity =  1.2
## Test of the hypothesis that 2 components are sufficient.
## 
## The root mean square of the residuals (RMSR) is  0 
##  with the empirical chi square  0  with prob <  NA 
## 
## Fit based upon off diagonal values = 1

回転すると、主成分の説明性は向上しますが、

• PCAの特徴である「第1主成分で全分散を最大化する」という性質は、回転によって失われます。
  

• 第一主成分、第二主成分という順序の重要性が低くなります。

• 分散合計（SS loadings の合計）はどちらも２に一致しています。

回転のメリットは、変数と主成分の関連が明確になり、分析結果の解釈が容易になることです。

9.2 prcomp( ) の出力

前節ではバリマックス回転を説明するため、principal( ) 関数を用いましたが、主成分分析で最もよく用いられるのは prcomp( )関数です。
違いを整理しておきます。

＜軸の回転（Rotation）＞

• prcomp(): 固有値・固有ベクトルを用いて分散が最大になる座標軸を順番に見つけるものなので、「回転」という概念はありません。
  
• principal(): 各変数がどの主成分に支配されているか、解釈しやすくするため回転を行うことができます。

＜負荷量（Loadings）の定義＞

• prcomp(): 出力される $rotation は固有ベクトルを列として構成されます。主成分負荷行列は主成分座標を元の座標に変換するための直交行列です。
  
• psych::principal(): 出力される loadings は、主成分と元の変数の相関係数です。 そのため、1つの主成分が各変数をどの程度説明しているかを直感的に把握できます。

「8.3 分析例」のデータを、prconp( ) で分析し、結果を比較します。

prcomp( )の出力は主成分の標準偏差とRotation（主成分方向）.

prcomp_result <- prcomp(data, scale = TRUE)
prcomp_result

## Standard deviations (1, .., p=2):
## [1] 1.2466625 0.6677069
## 
## Rotation (n x k) = (2 x 2):
##         PC1        PC2
## X 0.7071068  0.7071068
## Y 0.7071068 -0.7071068

Rotation の列ベクトルが、主成分の方向ベクトルを表します。

主成分を標準化した主成分負荷量にするには、固有値の平方根を乗じる必要があります.

sdev <- prcomp_result$sdev
ve <-prcomp_result$rotation#固有ベクトル
vn <- ve %*% diag(sdev);vn

##        [,1]       [,2]
## X 0.8815235  0.4721401
## Y 0.8815235 -0.4721401

principalは、分散を計算する際に自由度nを使用するため、prcompと計算結果がわずかに異なる場合があります（特にサンプルサイズが小さい場合）。

標準化したprcompのbiplot（左図） と、回転しない場合のprincipal（右図） を比べてみましょう。

par(mfrow=c(1,2))
biplot(prcomp_result,scale=1,main="Prcomp biplot")
biplot(pca_result,col = c("lightblue", "blue"),main="Principal biplot") 

上下反転していますが、同じものであることが分かります。

軸の目盛りは逆になっていることに注意しましょう。

• prcomp()：主成分座標の目盛りは下・左に、主成分負荷ベクトルの目盛りは上・右にあります。
  
• principal()：主成分座標の目盛りは上・右に、主成分負荷ベクトルの目盛りは下・左にあります。

主成分負荷ベクトルのscaleは見やすいように調整されており、実際の値の何倍かになっています。

10 三変量正規分布

10.1 座標変換と次元圧縮

主成分分析で用いる手法は「座標変換」と「次元圧縮」です。

＜座標変換＞
主成分分析は「分散が最大となる方向から見る」手法です。
そのため、データがどのような方向に広がっているか、その分布状況を把握する必要があります。

解決策を与えてくるのは、分散共分散行列の固有値・固有ベクトル分解です。
データは固有ベクトル方向に固有値の平方根倍に広がっているからです。

従って、固有ベクトルを基底とした座標を知ることが分析の第一歩です。

• 分散共分散行列は対称行列なので固有ベクトルは互いに直交しています。
  
• 固有ベクトルは正規直交基底です。
  
• 固有ベクトルを基底とした座標を求めることは、データの座標変換に他なりません。

＜次元圧縮＞.
データを次元圧縮するのは、データを可視化するためです。
データを見やすいように、地図上に MAP するのです。
しかし、多次元のデータを二次元の地図に MAP すれば、三次元の情報が失われます。
失われる情報を少なくするには、情報量の大きい座標を残す必要があります。

• 分散共分散行列の固有値は固有ベクトル方向の分散を表します。
  
• そこで、固有値の降順に固有ベクトルを選び、その方向の座標を主成分座標 PC1,PC2 ‥ とします。
  
• \((PC1,PC2)\) を地図上に MAP すれば、失われる情報量は最小となります。

二変量の主成分分析
これまでに扱った二変量の主成分分析では、座標変換だけが問題でした。
データは二次元なので次元圧縮の必要がなかったからです。
二変量正規分布の座標変換を振り返り、同じ手法で三変量正規分布の座標変換を考えます。

二変量正規分布は楕円状の分布であり、
　・楕円の長軸が第１主成分方向、短軸が第２主成分方向、
　・主成分方向は、分散共分散行列の固有ベクトルの方向です。
主成分座標を求めることは、固有ベクトルを基底とする座標への座標変換です。

• \(\Sigma=\begin{pmatrix} 1 & 0.7 \\ 0.7 & 1 \end{pmatrix}\)、\(\mu=\begin{pmatrix} 0 \\ 0 \end{pmatrix}\)　である正規乱数 dataを出力

固有ベクトルを基底とした主成分座標 \((PC1,PC2)\) に変換します。

library(mvtnorm)
Sigma <- matrix(c(1,0.7,0.7,1), ncol=2)
x <- rmvnorm(50, mean=c(0,0), sigma=Sigma)  
data <- data.frame(X=x[,1],Y=x[,2])
pca_result <- prcomp(data, scale = FALSE)
ve <-pca_result$rotation
par(mfrow=c(1,2))
plot(data,asp=1,col="lightblue",main="data",xlim=c(-4,4))
abline(a=0,b=ve[2,1]/ve[1,1],col="lightblue")
abline(a=0,b=ve[2,2]/ve[1,2],col="pink")
arrows(x0=0, y0=0, x1=ve[1,1],  y1=ve[2,1],col="blue",length = 0.1)
arrows(x0=0, y0=0, x1=ve[1,2],  y1=ve[2,2],col="red",length = 0.1)
plot(pca_result$x,asp=1,col="lightblue",main="data_PC",xlim=c(-4,4))
abline(h=0,col="lightblue")
abline(v=0,col="pink")
arrows(x0=0, y0=0, x1=1,  y1=0, code=2, col="blue",length = 0.1)
arrows(x0=0, y0=0, x1=0,  y1=1, code=2, col="red", length = 0.1)

矢線は固有ベクトルによる正規直交基底です。
楕円の長軸・短軸が主成分方向になっています。
主成分プロットは、固有ベクトルの向きに応じ上下反転・左右反転します。

10.2 三次元楕円体

三変量正規分布は三次元の楕円体状の分布になります。
三変量正規分布でも、固有ベクトルが主成分方向であることに変わりはありません。
イメージをつかむために、三次元楕円体について説明します。

楕円体の標準形

　　　\(\dfrac{x^2}{a^2} + \dfrac{y^2}{b^2} +\dfrac{z^2}{c^2}= 1 \quad (a > 0, b > 0, c > 0)\)

R の rgl パッケージの ellipse3d( ) 関数を用いて楕円体を描画します。
rgl による描画は、マウス操作によるインタラクティブな観察が可能です。
電子書籍ではできませんが、html上では可能です。
ひょっこり統計「高校生の主成分分析」で体験してください。

［例］
　　　\(\dfrac{x^2}{9} + \dfrac{y^2}{4} +\dfrac{z^2}{1}= 1\)　の３D表示.

options(rgl.printRglwidget = TRUE) 
library(rgl)
library(heplots)

shape_mat <- diag(c(9, 4, 1)) 
ellips <- ellipse3d(shape_mat, centre = c(0, 0, 0), t = 1)
shade3d(ellips, col = "lightblue", alpha = 0.2)
wire3d(ellips, col = "lightblue")
axes3d() 
ellipse3d.axes(shape_mat, centre = c(0, 0, 0), t=1, col = c("blue","red","green"), lwd = 2, labels = TRUE)
title3d(xlab = "X", ylab = "Y", zlab = "Z")
view3d(theta=20, phi=20)

　・軸の長さの順に、主成分軸 \(PC1,~PC2,~PC3\) となり、
　・次元圧縮では、第３主成分方向から見ることになります。

［データの例］
　現実のデータは広がりを持ちます。

• \(\Sigma=\begin{pmatrix} 9 & 0 & 0\\ 0 & 3 & 0\\ 0 & 0 & 1\ \end{pmatrix}\)、\(\mu=\begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix}\)　である正規乱数 data.

options(rgl.printRglwidget = TRUE) 
library(rgl)
library(MASS)
mu <- c(0,0,0)
Sigma <- matrix(c(9, 0,0,0,3,0,0,0, 1), 3, 3)
data <- mvrnorm(1000, mu, Sigma)
data <- data.frame(X=data[,1],Y=data[,2],Z=data[,3])
plot3d(data$X, data$Y, data$Z, type = "s", col = "lightblue", size = 0.3,xlim=c(-8,8),ylim=c(-8,8),zlim=c(-8,8),xlab="X",ylab="Y",zlab="Z")
abclines3d(0, 0, 0, a = diag(3), col = c("blue", "red", "green"))
view3d(theta=20, phi=10)

10.3 三変量の座標変換

一般の三変量正規分布について、固有ベクトルを基底とした主成分座標に変換します。
直交変換により、楕円体は標準形に座標変換されます。

• 分散共分散行列の固有ベクトルが、主成分方向であり、
  
• 固有値の大きさの順に、主成分軸 \(PC1,~PC2,~PC3\) となります。

主成分座標を求めることは、固有ベクトルを基底とする座標への座標変換です。

• \(\Sigma=\begin{pmatrix} 1 & 0.8 & 0.5\\ 0.8 & 1 & 0.3\\ 0.5 & 0.3 & 1\ \end{pmatrix}\)、\(\mu=\begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix}\)　に従う正規乱数 data を1000 個.

データを標準化しないで分析を行います。
（prcomp( )関数で、scale=FALSE と指定）

options(rgl.printRglwidget = TRUE) 
library(rgl)
library(MASS)
library(plot3D)

mu <- c(0,0,0)
Sigma <- matrix(c(1, 0.8,0.5,0.8,1,0.3,0.5,0.3, 1), 3, 3)
data <- mvrnorm(1000, mu, Sigma)
data <- data.frame(X=data[,1],Y=data[,2],Z=data[,3])
result <- prcomp(data,scale=FALSE)
ve <- result$rotation

close3d()
data <- data.frame(X=data[,1],Y=data[,2],Z=data[,3])
plot3d(data$X, data$Y, data$Z, type = "s", col = "lightblue", size = 0.4,xlim=c(-5,5),ylim=c(-5,5),zlim=c(-5,5),xlab="X",ylab="Y",zlab="Z")
abclines3d(0, 0, 0, ve[,1], col="blue")
abclines3d(0, 0, 0, ve[,2], col="red")
abclines3d(0, 0, 0, ve[,3], col="green")

上図をインタアクティブに回転し軸の方向を同じにすると、下図と一致することが確かめられます。
電子書籍ではできませんが、ひょっこり統計「高校生の主成分分析」でやってみましょう。

10.4 次元削減

次は、第３固有ベクトル方向から見た二次元 MAPを作成します。
この MAP は \((PC1,PC2)\)をプロットしたものです。

• 上図：元データ.
  
• 下図：二次元 MAP

固有ベクトルの大きさは１ですが、見やすくするため大きさを2倍にして表示します。
上図をアクティブに動かしてPC3方向から見た図と 二次元MAP が一致することを確かめましょう。

mu <- c(0,0,0)
Sigma <- matrix(c(1, 0.8,0.5,0.8,1,0.3,0.5,0.3, 1), 3, 3)
data <- mvrnorm(100, mu, Sigma)
data <- scale(data)
data <- data.frame(X=data[,1],Y=data[,2],Z=data[,3])
pca_result <- prcomp(data, scale = TRUE)
ve <-pca_result$rotation
data_PC <- pca_result$x[,1:2]

close3d()
plot3d(data$X, data$Y, data$Z, type = "s", col = "lightblue", size =0.5,xlim=c(-4,4),ylim=c(-4,4),zlim=c(-4,4),xlab="X",ylab="Y",zlab="Z")
arrow3d(p0 = c(0, 0, 0), p1 =2*ve[,1] , col = "blue",barblen = 0.3)
arrow3d(p0 = c(0, 0, 0), p1 =2*ve[,2] , col = "red",barblen = 0.3)
arrow3d(p0 = c(0, 0, 0), p1 =2*ve[,3] , col = "green",barblen = 0.3)
view3d(theta=0, phi=30)

plot(data_PC,asp=1,col="skyblue",main="data_PC",xlim=c(-5,5))
abline(v=0, h=0, col = c("lightblue", "red"))
arrows(x0=0, y0=0, x1=1,  y1=0, code=2, col="blue",lwd = 2)
arrows(x0=0, y0=0, x1=0,  y1=1, code=2, col="red", lwd = 2)

アクティブな操作は電子書籍ではできませんが、ひょっこり統計「高校生の主成分分析」でやってみましょう。

10.4.1 小標本

この操作を小標本で行い、データ番号を表示すると、対応がはっきりします。

close3d()
options(rgl.printRglwidget = TRUE) 
library(rgl)
library(MASS)
mu <- c(0,0,0)
Sigma <- matrix(c(1, 0.8,0.5,0.8,1,0.3,0.5,0.3, 1), 3, 3)
data <- mvrnorm(20, mu, Sigma)
data <- scale(data)
data <- data.frame(X=data[,1],Y=data[,2],Z=data[,3])
pca_result <- prcomp(data, scale = TRUE)
ve <-pca_result$rotation#固有ベクトル
data_PC <- pca_result$x[,1:2]

plot3d(data$X, data$Y, data$Z, type = "n",xlab="X",ylab="Y",zlab="Z")
arrow3d(p0 = c(0, 0, 0), p1 =ve[,1] , col = "blue",barblen = 0.3)
arrow3d(p0 = c(0, 0, 0), p1 = ve[,2] , col = "red",barblen = 0.3)
arrow3d(p0 = c(0, 0, 0), p1 =ve[,3] , col = "green",barblen = 0.3)
text3d(data$X, data$Y, data$Z,1:20,col="blue")

plot(data_PC,type="n",asp=1,col="lightblue",main="data_PC",xlim=c(-4,4))
abline(v=0, h=0, col = c("lightblue", "red"))
arrows(x0=0, y0=0, x1=1,  y1=0, code=2, col="blue",lwd = 2)
arrows(x0=0, y0=0, x1=0,  y1=1, code=2, col="red", lwd = 2)
text(1:20,x=pca_result$x[,1],y=pca_result$x[,2],col="blue")

アクティブな操作をひょっこり統計「高校生の主成分分析」でやってみましょう。

11 ３次元biplot

biplotは主成分座標と主成分負荷量ベクトルを一つのMAP に表示したものです。
二次元にすることで、表示される主成分負荷量ベクトルは情報を失っています。
そこで、
三次元のまま主成分負荷量ベクトルを表示し、biplot を二次元にした場合、どのように情報が失われているか観察しましょう。

「4.3 主成分負荷行列」を参照し、二次元の類比で三次元の主成分負荷量ベクトルを考えます。

11.1 主成分負荷行列

＜主成分を正規化しないとき＞.

　　\(\begin{pmatrix} x_1\\ x_2\\ x_3 \end{pmatrix} =\begin{pmatrix} p_1 & q_1 & r_1\\ p_2 & q_2 & r_2\\ p_3 & q_3 & r_3\\ \end{pmatrix}\begin{pmatrix} PC_1 \\ PC_2 \\ PC_3 \end{pmatrix}\)　

　　主成分負荷行列は　\(A=\begin{pmatrix} p_1 & q_1 & r_1\\ p_2 & q_2 & r_2\\ p_3 & q_3 & r_3\\ \end{pmatrix}\)

• \(x_1\) の主成分負荷ベクトル \((p_1, q_1, r_1)\) と 、
  
• \(x_2\) の主成分負荷ベクトル \((p_2, q_2, r_2)\) と 、
  
• \(x_3\) の主成分負荷ベクトル \((p_3, q_3, r_3)\) は互いに直交します。
  （ \(A\) は直交行列なので、行ベクトルも直交するため）

＜主成分を正規化するとき＞.

　 　\(\begin{pmatrix} x_1\\ x_2\\ x_3 \end{pmatrix} = \begin{pmatrix} \sqrt{\lambda_1}p_1 & \sqrt{\lambda_2}q_1 & \sqrt{\lambda_3}r_1\\ \sqrt{\lambda_1}p_2 & \sqrt{\lambda_2}q_2 & \sqrt{\lambda_3}r_2\\ \sqrt{\lambda_1}p_3 & \sqrt{\lambda_2}q_3 & \sqrt{\lambda_3}r_3\\ \end{pmatrix} \begin{pmatrix} z_1\\ z_2\\ z_3 \end{pmatrix}\)　

　　主成分負荷行列は　\(A=\begin{pmatrix} \sqrt{\lambda_1}p_1 & \sqrt{\lambda_2}q_1 & \sqrt{\lambda_3}r_1\\ \sqrt{\lambda_1}p_2 & \sqrt{\lambda_2}q_2 & \sqrt{\lambda_3}r_2\\ \sqrt{\lambda_1}p_3 & \sqrt{\lambda_2}q_3 & \sqrt{\lambda_3}r_3\\ \end{pmatrix}\)

\(V(x_1)=\lambda_1p_1^2+\lambda_2q_1^2+\lambda_3r_1^2\) 、\(V(x_2)=\lambda_1p_2^2+\lambda_2q_2^2+\lambda_3r_2^2\)、
\(V(x_3)=\lambda_1p_3^2+\lambda_2q_3^2+\lambda_3r_3^2\)　 であるので、

• \(x_1\) の主成分負荷ベクトル \((\sqrt{\lambda_1}p_1, \sqrt{\lambda_2}q_1,\sqrt{\lambda_3}r_1)\) の大きさは、\(x_1\)の標準偏差に等しく、
  
• \(x_2\) の主成分負荷ベクトル \((\sqrt{\lambda_1}p_2, \sqrt{\lambda_2}q_2,\sqrt{\lambda_3}r_2)\) の大きさは、\(x_2\)の標準偏差に等しく、
  
• \(x_3\) の主成分負荷ベクトル \((\sqrt{\lambda_1}p_3, \sqrt{\lambda_2}q_3,\sqrt{\lambda_3}r_3)\) の大きさは、\(x_3\)の標準偏差に等しくなります。
  A は直交行列ではなく、\(x_1,~x_2,~x_3\) の主成分負荷ベクトルは直交しません。

11.2 主成分座標を標準化しないとき

• \(\Sigma=\begin{pmatrix} 1 & 0.8 & 0.5\\ 0.8 & 1 & 0.3\\ 0.5 & 0.3 & 1\ \end{pmatrix}\)、\(\mu=\begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix}\)　に従う正規乱数 data を20 個.

scale( )関数で標準化したうえで分析を行います。
（標本分散共分散行列は、完全には \(\Sigma\) に一致しません）.

close3d()
options(rgl.printRglwidget = TRUE) 
library(rgl)
library(MASS)
mu <- c(0,0,0)
Sigma <- matrix(c(1, 0.8,0.5,0.8,1,0.3,0.5,0.3, 1), 3, 3)
data <- mvrnorm(20, mu, Sigma)
data <- scale(data)
data <- data.frame(X=data[,1],Y=data[,2],Z=data[,3])
pca_result <- prcomp(data, scale = TRUE)
ve <-pca_result$rotation#固有ベクトル

これは、「1.4.1 小標本」と同じ data を用いています。
三次元 biplot は \(( PC1,~PC2,~PC3)\) の描画であり、楕円体の標準形の形状になっています。

主成分負荷量ベクトルは、rotation（固有ベクトルを列とする行列）の行ベクトルです。
これは大きさが１で、互いに直交しています。 見やすくするため、大きさを4倍にして表示します。

biplot <- pca_result$x
data <- data.frame(PC1=biplot[,1],PC2=biplot[,2],PC3=biplot[,3])
plot3d(data$PC1, data$PC2, data$PC3, type = "n",xlim=c(-4,4),ylim=c(-4,4),zlim=c(-4,4),xlab="X",ylab="Y",zlab="Z")
aspect3d(1, 1, 1)
text3d(data$PC1, data$PC2, data$PC3,1:20,col="blue")
arrow3d(p0 = c(0, 0, 0), p1 =4*ve[1,] , col = "blue")
text3d(x=4.3*ve[1,1], y=4.3*ve[1,2], z=4.3*ve[1,3], texts = "X", col = "blue")
arrow3d(p0 = c(0, 0, 0), p1 =4*ve[2,] , col = "red")
text3d(x=4.3*ve[2,1], y=4.3*ve[2,2], z=4.3*ve[2,3], texts = "Y", col = "red")
arrow3d(p0 = c(0, 0, 0), p1 =4*ve[3,] , col = "green")
text3d(x=4.3*ve[3,1], y=4.3*ve[3,2], z=4.3*ve[3,3], texts = "Z", col = "green")

３つの主成分負荷ベクトルの大きさは、

c(sum(ve[1,]^2),sum(ve[2,]^2),sum(ve[3,]^2))

## [1] 1 1 1

主成分負荷ベクトルは互いに直交する単位ベクトルです。

次に、二次元のbiplot を描きます。

biplot(pca_result,scale=0)

二次元のbiplot では、X と Y の主成分負荷ベクトルが直交することは読み取れません。
射影によって情報が失われていることに注意しておく必要があります。
また、長さも異なるように見えますが、これは射影のためで、全て大きさは１です。

11.3 主成分座標を標準化するとき

• 主成分負荷行列を求めるには、rotation（固有ベクトルを列とする行列）に固有値の平方根を乗じる必要があります。
  
• 主成分座標を標準化するには、固有値の平方根 (sdev)で除する必要があります。

行列を用いた計算では、

• 主成分負荷行列（vn）の計算には、固有値の平方根を対角にもつ対角行列 diag(sdev) を右から乗じ、
  
• 標準化した主成分座標（biplotn）の計算には、固有値の平方根の逆数を対角にもつ対角行列 diag(1/sdev) を右から乗じます。

sdev <- pca_result$sdev
vn <- ve %*% diag(sdev);vn

##         [,1]       [,2]         [,3]
## X -0.9316835  0.2700844  0.242940873
## Y -0.9331439  0.2640566 -0.243960359
## Z -0.7377363 -0.6750867  0.001770421

biplotn <- pca_result$x %*% diag(1/sdev)
data <- data.frame(PC1=biplotn[,1],PC2=biplotn[,2],PC3=biplotn[,3])
plot3d(data$PC1, data$PC2, data$PC3, type = "n",xlim=c(-4,4),ylim=c(-4,4),zlim=c(-4,4),xlab="X",ylab="Y",zlab="Z")
aspect3d(1, 1, 1)
text3d(data$PC1, data$PC2, data$PC3,1:20,col="blue")
arrow3d(p0 = c(0, 0, 0), p1 =4*vn[1,] , col = "red")
text3d(x=4.2*vn[1,1], y=4.2*vn[1,2], z=4.2*vn[1,3], texts = "X", col = "blue")
arrow3d(p0 = c(0, 0, 0), p1 =4*vn[2,] , col = "blue")
text3d(x=4.2*vn[2,1], y=4.2*vn[2,2], z=4.2*vn[2,3], texts = "Y", col = "blue")
arrow3d(p0 = c(0, 0, 0), p1 =4*vn[3,] , col = "green")
text3d(x=4.2*vn[3,1], y=4.2*vn[3,2], z=4.2*vn[3,3], texts = "Z", col = "blue")

• \(x_1\) の主成分負荷ベクトル \((\sqrt{\lambda_1}p_1, \sqrt{\lambda_2}q_1,\sqrt{\lambda_3}r_1)\) の大きさは、\(x_1\)の標準偏差に等しく、
  
• \(x_2\) の主成分負荷ベクトル \((\sqrt{\lambda_1}p_2, \sqrt{\lambda_2}q_2,\sqrt{\lambda_3}r_2)\) の大きさは、\(x_2\)の標準偏差に等しく、
  
• \(x_3\) の主成分負荷ベクトル \((\sqrt{\lambda_1}p_3, \sqrt{\lambda_2}q_3,\sqrt{\lambda_3}r_3)\) の大きさは、\(x_2\)の標準偏差に等しくなります。
  今、\(x_1,x_2,x_3\) は標準化されているので、主成分負荷ベクトルは大きさはすべて１です。

c(sum(vn[1,]^2),sum(vn[2,]^2),sum(vn[3,]^2))

## [1] 1 1 1

また、主成分負荷行列は直交行列ではなく、\(x_1,~x_2,~x_3\) の主成分負荷ベクトルは直交しません。
　
次に、二次元のbiplot を描きます。

biplot(pca_result,scale=1)

XとYのなす角は小さく見えますが、三次元で見ると小さくありません。
二次元biplotから読み取れる情報は、二次元に削減された情報です。
XとYの正確な相関係数を求めるには、三次元で計算する必要があります。

「5.2」で確認したように、
変量 X,Y の主成分負荷ベクトルのなす角の余弦は、変量 X,Y の相関係数に一致します。
変量 X,Y の主成分負荷ベクトルは単位ベクトルなので、変量 X,Y の相関係数は主成分負荷ベクトルの内積として計算されます。

sum(vn[1,]*vn[2,])#XとYの相関係数

## [1] 0.8814444

相関係数は余弦の値なので、三次元における X,Y のなす角は cosの逆関数 acos( ) によって計算されます。

acos(sum(vn[1,]*vn[2,]))* 180 / pi#XとYのなす角

## [1] 28.18291

二次元biplotにおける X,Y のなす角は、三次元biplotでの実際の角度 より小さくなります。

12 回転と寄与率

Rには主成分分析の関数として、基本関数である prcomp( ) と、psychパッケージの principal( ) があります。
（「9 分析の出力」参照）

• principal( ) は心理学でよく用いられ、変量及び主成分得点を標準化して扱います。
• principal( ) では “rotate” オプションで回転を指定することができ、主成分の意味を解釈するのに役立ちます。
  
• 寄与率は、Proportion Var で出力されます。

3変量データを乱数で出力し、principal( )関数を適用してみましょう。

12.1 principal( )の出力

• \(\Sigma=\begin{pmatrix} 1 & 0.8 & 0.5\\ 0.8 & 1 & 0.2\\ 0.5 & 0.2 & 1\ \end{pmatrix}\)、\(\mu=\begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix}\)　に従う正規乱数 data を
  20 個生成します。

library(MASS)
mu <- c(0, 0, 0) # 平均ベクトル
Sigma <- matrix(c(1.0, 0.8, 0.5,
                  0.8, 1.0, 0.2,
                  0.5, 0.2, 1.0), nrow = 3, byrow = TRUE)
Sigma

##      [,1] [,2] [,3]
## [1,]  1.0  0.8  0.5
## [2,]  0.8  1.0  0.2
## [3,]  0.5  0.2  1.0

set.seed(1114537)
x <- mvrnorm(n = 20, mu = mu, Sigma = Sigma)
x<- scale(x)
data <- data.frame(x = x[,1], y = x[,2],z = x[,3])

回転について、

• 回転しない場合は、rotate = “none”.
  
• varimax 回転する場合は、rotate = “varimax”.

と指定します。

biplot について、
回転しない場合（左図）と回転する場合（右図）を比較しましょう。

• 左図の軸 PC1, PC2は “Principal Component”（主成分）
  
• 右図の軸 RC1, RC2は “Rotated Component”（回転された主成分）

を意味します。

library(psych)
pca_result <- principal(data, nfactors = 2, rotate = "none")
vpca_result <- principal(data, nfactors = 2, rotate = "varimax")
par(mfrow=c(1,2))
biplot(pca_result,col = c("lightblue", "blue"),main="Rotate : none") 
biplot(vpca_result,col = c("lightblue", "blue"),main="Rotate : varimax") 

変量 X,Y は第１主成分で説明され、変量 Z は第２主成分で説明される単純構造に近いことが観察されます。

• 回転をしない場合の出力 pca_result.
  
• varimax回転をする場合の出力 vpca_result.

を比較しましょう。

＜回転しない場合＞

pca_result

## Principal Components Analysis
## Call: principal(r = data, nfactors = 2, rotate = "none")
## Standardized loadings (pattern matrix) based upon correlation matrix
##    PC1   PC2   h2     u2 com
## x 0.96 -0.15 0.94 0.0599 1.0
## y 0.90 -0.38 0.95 0.0486 1.3
## z 0.61  0.79 1.00 0.0036 1.9
## 
##                        PC1  PC2
## SS loadings           2.10 0.79
## Proportion Var        0.70 0.26
## Cumulative Var        0.70 0.96
## Proportion Explained  0.73 0.27
## Cumulative Proportion 0.73 1.00
## 
## Mean item complexity =  1.4
## Test of the hypothesis that 2 components are sufficient.
## 
## The root mean square of the residuals (RMSR) is  0.03 
##  with the empirical chi square  0.13  with prob <  NA 
## 
## Fit based upon off diagonal values = 1

• ss loading は各主成分の分散、全分散は３です。
  
• 第１主成分の寄与率（Proportion Var）は 70%、第２主成分の寄与率は 26%、
  第２主成分までの累積寄与率（Cumulative Var）は 96% です。
  
• com（影響を受ける因子の数）は、変量 z が 1.9個と多く、3変量の平均個数は “Mean item complexity” で1.4 個です。
  
• 第１主成分の分散が大きく、第１主成分の係数の絶対値が大きいです。

次に、varimax 回転の指標を計算するため、主成分負荷行列の成分を二乗します。

s_loading <- pca_result$loadings[,]^2
s_loading

##         PC1        PC2
## x 0.9184829 0.02157064
## y 0.8078703 0.14353431
## z 0.3694786 0.62691365

列ごとの分散は、単純構造に近づくと大きくなります。
列ごとの分散は、apply( )関数で第２引数を2にすると計算されます。

apply(s_loading, 2, var)

##        PC1        PC2 
## 0.08430468 0.10249515

列ごとの分散を合計します。

sum(apply(s_loading, 2, var))

## [1] 0.1867998

var( )を用いた手計算も可能です。

var(pca_result$loadings[,1]^2)+var(pca_result$loadings[,2]^2)

## [1] 0.1867998

各主成分ごとの係数二乗値の分散の合計は、0.1867998 です。
（varimax回転では、この指標が最大となるよう回転します）

＜回転する場合＞

vpca_result

## Principal Components Analysis
## Call: principal(r = data, nfactors = 2, rotate = "varimax")
## Standardized loadings (pattern matrix) based upon correlation matrix
##    RC1  RC2   h2     u2 com
## x 0.92 0.31 0.94 0.0599 1.2
## y 0.97 0.07 0.95 0.0486 1.0
## z 0.18 0.98 1.00 0.0036 1.1
## 
##                        RC1  RC2
## SS loadings           1.82 1.06
## Proportion Var        0.61 0.35
## Cumulative Var        0.61 0.96
## Proportion Explained  0.63 0.37
## Cumulative Proportion 0.63 1.00
## 
## Mean item complexity =  1.1
## Test of the hypothesis that 2 components are sufficient.
## 
## The root mean square of the residuals (RMSR) is  0.03 
##  with the empirical chi square  0.13  with prob <  NA 
## 
## Fit based upon off diagonal values = 1

• 第１主成分の寄与率（Proportion Var）は 61%、第２主成分の寄与率は 35%、
  第２主成分までの累積寄与率（Cumulative Var）は 96% で変わりません。
  
• com（影響を受ける因子の数）は、どの変量も１に近くなっています。
  3変量の平均個数 “Mean item complexity” は 1.1 個で、単純構造に近づいていることが分かります。
  
• 第１主成分の分散が大きいですが、分散の偏りは小さくなっています。

次に、varimax回転の指標を計算するため、主成分負荷行列の成分を二乗し、列ごとの分散を合計します。

s_vloading <- vpca_result$loadings[,]^2
s_vloading

##          RC1         RC2
## x 0.84613168 0.093921841
## y 0.94609736 0.005307221
## z 0.03231067 0.964081605

apply(s_vloading, 2, var)

##       RC1       RC2 
## 0.2512173 0.2807131

sum(apply(s_vloading, 2, var))

## [1] 0.5319305

var( )を用いた手計算では、

var(vpca_result$loadings[,1]^2)+var(vpca_result$loadings[,2]^2)

## [1] 0.5319305

各主成分ごとの係数二乗値の分散の合計は 0.5319305 で、回転しない場合より大きくなっています。
varimax回転は、この指標が最大となる回転を行っています。

12.2 prcomp()との比較

prcomp( ) は、固有値・固有ベクトルを用いた主成分分析の基本関数です。

prcomp_result <- prcomp(data, scale = TRUE)
prcomp_result

## Standard deviations (1, .., p=3):
## [1] 1.4476988 0.8899543 0.3348875
## 
## Rotation (n x k) = (3 x 3):
##         PC1        PC2        PC3
## x 0.6619990 -0.1650304  0.7311103
## y 0.6208585 -0.4257064 -0.6582620
## z 0.4198715  0.8896848 -0.1793567

prcomp() で出力される $Rotation は固有ベクトルを列として構成される行列です。
Rotation の列ベクトルは、主成分の方向ベクトルを表します。

主成分を標準化した主成分負荷量にするには、固有値の平方根を乗じる必要があります。

sdev <- prcomp_result$sdev
ve <-prcomp_result$rotation#固有ベクトル
vn <- ve %*% diag(sdev);vn

##        [,1]       [,2]        [,3]
## x 0.9583751 -0.1468695  0.24483969
## y 0.8988160 -0.3788592 -0.22044369
## z 0.6078475  0.7917788 -0.06006431

この行列の２列目までを用いれば、principal()で得た主成分負荷行列を得ます。

round(vn[,1:2],2)

##   [,1]  [,2]
## x 0.96 -0.15
## y 0.90 -0.38
## z 0.61  0.79

標準化したprcompのbiplot（左図） と、回転しない場合のprincipal（右図） を比べてみましょう。

par(mfrow=c(1,2))
biplot(prcomp_result,scale=1,main="Prcomp biplot")
biplot(pca_result,col = c("lightblue", "blue"),main="Principal biplot") 

軸の目盛りは、

• prcomp()：主成分座標の目盛りは下・左に、主成分負荷ベクトルの目盛りは上・右にあります。
• principal()：主成分座標の目盛りは上・右に、主成分負荷ベクトルの目盛りは下・左にあります。

プロットの配置は、ほぼ一致していることが分かります。

主成分負荷ベクトルのscaleは見やすいように調整されており、実際の値の何倍かになっています。

13 終わりに

好きな映画に「スパルタカス」があります。
スパルタカスはローマの奴隷で剣闘士でしたが、解放を求めて立ち上がりました。
「俺は自由だ、だが何も知らない、字さえ読めない」「だから知りたい、どうしても」 「なぜ星は落ち、鳥は落ちない」「太陽は夜どこへ行く、なぜ月の形は変わる」 自由を得たスパルタカスは、この世界の不思議を学びたかったのです。

私は退職後に統計学を学び始めました。学生時代、よく分からなかったからです。
gaccoで統計学Ⅲ「多変量データ解析法」を受講しました。
ところが、統計学Ⅲのテキストp151に
「主成分分析は、数学的にいえば分散共分散行列のスペクトル分解に他ならない」
と言う言葉があり考え込んでしまいました。
私は「主成分分析は統計学の特殊な方法」と思い込み、その背景にある基礎を理解していなかったからです。

今日ではICTやAIの進歩により、少しの操作でもっともらしい資料が得られます。
しかし、その意味は理解できているでしょうか。その妥当性を評価できるでしょうか。
本書では、
納得いくまで確かめるため、二次元・三次元でシミュレーションを行いました。
三次元のインタラクティブな操作で分布を回転させ、主成分分析の仕組みを確かめることができました。
　「主成分分析は、データを見えやすい方向から見るための基底変換である」
理解ができず心の中にずっとひっかかっていたことが、教材をつくる中で見えてきて喜んでいます。

高校生は線形代数を学習していませんが、視覚的に理解できるよう工夫しました。
そして、将来の学習の参考となるよう、数理への入り口を紹介しています。線形代数で行列を学んでから読み返せば、きっと腑に落ちてくることでしょう。