主成分分析と因子分析

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

主成分分析

ここではデータとして2022年度全国学力・学習状況調査の結果を使う：

df = read.csv("http://okumuralab.org/~okumura/python/data/atest2022.csv")

頭の部分だけ表示してみる：

head(df)

      小国     小算     小理     中国     中数     中理
1 64.44456 61.07105 62.87208 68.59639 48.93763 48.96912
2 67.81161 63.19436 65.83762 69.13618 51.55864 48.98470
3 66.98455 61.59387 63.19816 69.80850 48.52725 47.56724
4 63.68711 60.15438 61.49521 69.14642 48.57422 49.72042
5 70.78273 66.45425 70.61405 72.50219 54.37741 52.11749
6 64.57035 61.05593 63.34265 69.51811 51.91599 50.19502

わかりやすいように、都道府県名を付けておく：

kenmei = c("北海道", "青森", "岩手", "宮城", "秋田", "山形", "福島",
  "茨城", "栃木", "群馬", "埼玉", "千葉", "東京", "神奈川", "新潟",
  "富山", "石川", "福井", "山梨", "長野", "岐阜", "静岡", "愛知", "三重",
  "滋賀", "京都", "大阪", "兵庫", "奈良", "和歌山", "鳥取", "島根", "岡山",
  "広島", "山口", "徳島", "香川", "愛媛", "高知", "福岡", "佐賀", "長崎",
  "熊本", "大分", "宮崎", "鹿児島", "沖縄")
row.names(df) = kenmei
head(df)

           小国     小算     小理     中国     中数     中理
北海道 64.44456 61.07105 62.87208 68.59639 48.93763 48.96912
青森   67.81161 63.19436 65.83762 69.13618 51.55864 48.98470
岩手   66.98455 61.59387 63.19816 69.80850 48.52725 47.56724
宮城   63.68711 60.15438 61.49521 69.14642 48.57422 49.72042
秋田   70.78273 66.45425 70.61405 72.50219 54.37741 52.11749
山形   64.57035 61.05593 63.34265 69.51811 51.91599 50.19502

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

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

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

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

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

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

biplot(prcomp(df))

字が一部欠ける場合は par(xpd=TRUE) というオマジナイをしてから描くとよい。

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

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

biplot(prcomp(df, scale=TRUE))

Rの主成分分析の関数は prcomp() と princomp() の二つがあるが、より新しい prcomp() のほうが推奨。

因子分析

ここでは都道府県は考えずに、6科目の相関係数だけを考えることにする。 

cor(df)

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

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

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

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

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

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