主成分分析と因子分析

[追記] より詳しく説明した特異値分解・主成分分析・バイプロットという記事も書きました。

主成分分析

ここではデータとして平成26年度全国学力・学習状況調査の結果を使う。中学校・小学校の結果をあらかじめ読み込んでおく：

chu = read.csv("http://okumuralab.org/~okumura/stat/data/atest2014chu.csv",
               fileEncoding="CP932")
sho = read.csv("http://okumuralab.org/~okumura/stat/data/atest2014sho.csv",
               fileEncoding="CP932")

中学校の結果の頭の部分だけ表示してみる：

> head(chu)
  都道府県 国語A 国語B 数学A 数学B
1   北海道  79.4  49.9  66.0  59.4
2   青森県  81.0  52.0  69.3  60.7
3   岩手県  80.5  51.8  64.1  57.5
4   宮城県  80.3  52.0  65.6  59.4
5   秋田県  84.4  55.8  73.0  65.5
6   山形県  80.6  52.1  67.6  60.2

1列目 chu[ ,1] は都道府県名だから，いらない。2〜5列目 chu[ ,2:5] だけでいい。でも消すのはもったいないので，都道府県名は行の名前にして，それから都道府県の列を外すことにする：

> row.names(chu) = chu[,1]
> chu = chu[,2:5]
> head(chu)
       国語A 国語B 数学A 数学B
北海道  79.4  49.9  66.0  59.4
青森県  81.0  52.0  69.3  60.7
岩手県  80.5  51.8  64.1  57.5
宮城県  80.3  52.0  65.6  59.4
秋田県  84.4  55.8  73.0  65.5
山形県  80.6  52.1  67.6  60.2

各科目の平均点は気にしないので，以下では点数からその科目の平均点を引いた値（偏差）を使うことにする。ここで，すべての値は「ほげ」と「ふが」という二つの独立な隠れた値で決まるという大胆な仮説を立てる：

北海道の国語Aの偏差 = 北海道のほげ × 国語Aのほげ＋北海道のふが × 国語Aのふが
青森県の国語Aの偏差 = 青森県のほげ × 国語Aのほげ＋青森県のふが × 国語Aのふが
……

「ほげ」と「ふが」が独立なので，北海道から沖縄県まで47成分ある「ほげ」ベクトルと「ふが」ベクトルは直交する（内積0）。また，国語Aから理科まで5成分ある「ほげ」ベクトルと「ふが」ベクトルも直交する。これだけではまだ決まらないので，以下では5成分のベクトルはどれも長さが1で，47成分のベクトルは「ほげ」の長さが最大になるように選ぶ。ここで，ベクトルの長さとは，自分自身との内積の平方根である（つまり成分の2乗和の平方根である）。

もちろん正確に「ほげ」と「ふが」だけに分解することは不可能だが，近似的に上のように分解できればうれしい。これが主成分分析の考え方である。これは数学的には「特異値分解」（singular value decomposition，SVD）というものに相当する。

主成分分析は，「ほげ」「ふが」の二つに限らず，三つでも四つでも好きな数に分解すればいいのだが，後で結果を「バイプロット」（biplot）という2次元の図に描きたいので，ここでは二つに限定した。

Rで実際に計算してみよう。主成分分析し，結果をバイプロットで描く。上のいいかげんな説明の「ほげ」が横軸，「ふが」が縦軸である。「ほげ」「ふが」の正式な名前は「第1主成分」「第2主成分」である。図ではPC1，PC2と書かれているが，PCは主成分（principal component）のことである。

biplot(prcomp(chu))

よく見ると字が一部欠けている。これを防ぐには

par(xpd=TRUE)

というオマジナイをしてから描くとよい。

この図から次のことがわかる：

二つの数学，二つの国語は近い関係にある。
数学のほうが国語より長いのは，標準偏差が大きいことを表す。
福井県や秋田県が優秀だが，福井県は数学，秋田県は国語が特に良い。
沖縄県は残念な結果だった。

上の例では，各科目の標準偏差をわざと揃えなかった。標準偏差の違いを考えないのであれば，各科目の点数から平均を引いた後で標準偏差で割るのがよい。このようにしてから主成分分析を行う場合は，次のようにする：

biplot(prcomp(chu, scale=TRUE))

実際に描いて結果を比べてみるとよい。また，小学校の場合もやってみるとよい。

row.names(sho) = sho[,1]
sho = sho[,2:5]
biplot(prcomp(sho))
biplot(prcomp(sho, scale=TRUE))

Rの主成分分析の関数は prcomp() と princomp() の二つがあるが，より新しい prcomp() のほうを使うべきである。

因子分析

4科目では簡単すぎるので，小学校と合わせて8科目にしてみよう。

> all = data.frame(sho, chu)
> names(all) = c("小国A","小国B","小算A","小算B","中国A","中国B","中数A","中数B")
> head(all)
       小国A 小国B 小算A 小算B 中国A 中国B 中数A 中数B
北海道  71.8  52.9  75.8  55.2  79.4  49.9  66.0  59.4
青森県  76.6  60.5  81.3  60.8  81.0  52.0  69.3  60.7
岩手県  73.7  58.3  78.9  58.7  80.5  51.8  64.1  57.5
宮城県  74.2  54.3  77.3  56.8  80.3  52.0  65.6  59.4
秋田県  77.4  67.3  85.1  66.2  84.4  55.8  73.0  65.5
山形県  74.3  56.7  77.8  57.7  80.6  52.1  67.6  60.2

これを主成分分析・バイプロットしてもいいが，ここでは都道府県は考えずに，8科目の相関係数だけを考えることにする。

cor(all)

と打ち込むと，すべての相関係数の行列が出力される。もちろん対角成分（同じ科目どうしの相関係数）は1であるが，対角成分に意味はない。そこで，各科目に相当するベクトルを作って，異なるベクトル間の内積が相関係数にできるだけ等しくなるようにするという問題を考える。これが因子分析（factor analysis）の問題である。

具体的には次のようにすればよい：

r = factanal(all, factors=2)
plot(NULL, xlim=c(0,1), ylim=c(0,1), xlab="因子1", ylab="因子2")
text(r$loadings, names(all))

因子分析では，各科目を表すベクトルどうしの内積だけに意味があるので，原点を中心に回転しても変わらない。この自由度を使って，通常はできるだけ解釈のしやすい座標になるように計算される。具体的には，座標軸のそばにベクトルが多く集まるように回転した結果を出力する（バリマックス回転）。

これから次のことがわかる：

小学校と中学校がはっきり分かれる。
中学校では国語のほうが小学校の成績に近い。
小学校では算数Bが中学校の成績に近い。
小学校の国語Aは他教科との相関がやや小さい。

なお，因子分析でもバイプロットを描くことがある：

r = factanal(all, factors=2, scores="regression")
biplot(r$scores, r$loadings)