ここでは因子分析・主成分分析について簡単にメモる。
参考文献: RとPythonで学ぶ[実践的]データサイエンス&機械学習[増補改訂版]
因子分析
簡単な例として英国数理社の5科目のテストについて考えてみる。
まず仮定として、数学や理科では計算を多くこなすとする。一方で社会や国語は読解力が大いに関係するとする。また日本語が不自由であれば数学や理科の題意を理解できない。さらに社会でもグラフや表の読み取るために数理能力が必要になるかもしれない。これらのことを踏まえて、理科の得点を計算式で表そうとすると、
理科の得点 = a_1*計算力 + a_2 * 読解力 + 理科に求められる力
とできる。$a_1$及び$a_2$は各々の変数が理科の得点に与える影響力を示す。上記の仮定を踏まえると、計算力に関する$a_1$の値は大きく$a_2$の値はそれほど大きくないと予想できる。
X_1 = a_{11}f_1 + a_{12}f_2 + e_1 \\
X_2 = a_{21}f_1 + a_{22}f_2 + e_2 \\
X_3 = a_{31}f_1 + a_{32}f_2 + e_3 \\
X_4 = a_{41}f_1 + a_{42}f_2 + e_4 \\
X_5 = a_{51}f_1 + a_{52}f_2 + e_5
$X_i$はi(科目名)のテストの得点を示す。aは5科目と共通能力(計算力、読解力)の2つであるため合計10個のパラメータが存在する。
$a_{ij}$ 因子負荷量(共通因子の重み)
$f_i$ 共通因子(抽出される因子)
$e_i$ 独自因子
主成分分析
主成分分析の考え方は因子分析とほぼ同じであり、共通因子→主成分となるイメージ。因子分析は共通因子で説明できない箇所を独自因子として残すが、主成分分析は全て主成分で事象を説明する。つまり上記の例(5科目のテストの点数)であれば、主成分は5つ。ここから上位の主成分を選択することで次元圧縮を行う。
X_1 = a_{11}f_1 + a_{12}f_2 + a_{13}f_3 + a_{14}f_4 + a_{15}f_5\\
X_2 = a_{21}f_1 + a_{22}f_2 + a_{23}f_3 + a_{24}f_4 + a_{25}f_5 \\
X_3 = a_{31}f_1 + a_{32}f_2 + a_{33}f_3 + a_{34}f_4 + a_{35}f_5 \\
X_4 = a_{41}f_1 + a_{42}f_2 + a_{43}f_3 + a_{44}f_4 + a_{45}f_5 \\
X_5 = a_{51}f_1 + a_{52}f_2 + a_{53}f_3 + a_{54}f_4 + a_{55}f_5
$a_{ij}$ 固有ベクトル (主成分に対する重み)
$f_i$ 主成分 測定項目から抽出される主成分
次元削除
まずは今回使用するデータセットの説明。
1.市町村 東京都の自治体
2.行政CD
3.人気度 (住みたい場所ほど点数が高い)
・
・
と続く。また説明変数となる行政指数(今回は・・で示しているところ)は25存在する。これらを人気度を目的変数とする回帰分析を行うと、項目間の相関が高いくなる、多重共線性が発生する可能性がある。また対象となる自治体が50以上あるためか学習を起こす可能性もある。これらを避けるために因子分析や主成分分析を使用し、次元圧縮する必要がある。ここではまずいくつの共通因子を取り出すかを決めるため、並行分析と言われる手法を使って目安となる因子数や主成分数を計算する。
> result.pr1 <- fa.parallel(DF[, -c(1:3)], fm="ml")
Parallel analysis suggests that the number of factors = 3 and the number of components = 2
出力結果を見ると"the number of factors = 3"で因子数は3つ、 "the number of components = 2 "主成分数は2つが適切という結果が見られる。またグラフからもわかるように第1因子や第1主成分は多くのデータの分散が含まれているが、因子数や主成分数が大きくなるにつれて分散の説明率は低下していっている。目安として赤い点線を見ればいくつ因子数及び主成分数を選べば良いかわかる。
因子分析の実装
#fa:因子分析
#fm:因子抽出法
#最尤法
#nfactors:因子数
#rotate:回転砲(直行回転:varimax, 斜交回転:promax)
#scores:因子得点算出法
#1~2列目は削除
#3列目は対象外
resultFA <- fa(DF[, -c(1:3)],
nfactors=3,
fm="ml",
rotate="varimax",
scores="regression")
> print(resultFA, digits=2, sort=TRUE)
Factor Analysis using method = ml
Call: fa(r = DF[, -c(1:3)], nfactors = 3, rotate = "varimax", scores = "regression",
fm = "ml")
Standardized loadings (pattern matrix) based upon correlation matrix
item ML1 ML2 ML3 h2 u2 com
千人あたり事業所数 18 0.97 0.16 0.14 1.00 0.0046 1.1
昼間人口比_per 6 0.97 0.11 0.12 0.97 0.0310 1.1
千人あたり大型小売店数 21 0.95 0.12 0.20 0.97 0.0327 1.1
千人あたり飲食店数 20 0.95 0.22 0.19 0.98 0.0171 1.2
千人あたり交通事故発生件数 24 0.95 0.00 0.14 0.92 0.0828 1.0
千人あたり刑法犯認知件数 25 0.91 0.23 0.05 0.88 0.1202 1.1
千人あたり幼稚園数 19 0.80 0.24 0.14 0.71 0.2866 1.2
千人あたり病院数 22 0.79 -0.01 -0.26 0.70 0.3041 1.2
転入者_対人口比 4 0.65 0.62 0.41 0.98 0.0240 2.7
課税所得_就業者1人あたり_千円 13 0.50 0.43 0.43 0.62 0.3762 2.9
小売業販売額_事業所あたり_百万円 14 0.49 0.17 0.45 0.47 0.5291 2.2
世帯あたり人数 1 -0.16 -0.92 -0.28 0.95 0.0504 1.3
年齢15未満比率 2 -0.06 -0.91 -0.09 0.83 0.1684 1.0
高齢単身世帯比率 7 0.03 0.74 -0.65 0.96 0.0376 2.0
転出者_対人口比 5 0.48 0.71 0.49 0.97 0.0282 2.6
耕地面積_対可住面積比 12 -0.13 -0.64 -0.11 0.45 0.5546 1.1
千人あたり老人ホーム数 23 0.01 -0.63 -0.16 0.43 0.5712 1.1
可住地面積比率 11 0.01 0.60 0.06 0.36 0.6363 1.0
小売業販売額_売場面積あたり_万円m2 15 0.44 0.58 0.34 0.65 0.3502 2.5
第3次産業従業者数比 10 0.09 0.54 0.23 0.35 0.6497 1.4
第1次産業従業者数比 8 -0.20 -0.54 0.08 0.34 0.6635 1.3
第2次産業従業者数比 9 -0.09 -0.54 -0.23 0.35 0.6536 1.4
ごみリサイクル率_pct 17 -0.24 -0.53 0.11 0.36 0.6444 1.5
年齢65以上比率 3 -0.10 -0.19 -0.85 0.77 0.2340 1.1
国民健保一人あたり診療費_円 16 -0.15 -0.34 -0.48 0.36 0.6382 2.0
ML1 ML2 ML3
SS loadings 8.22 6.34 2.75
Proportion Var 0.33 0.25 0.11
Cumulative Var 0.33 0.58 0.69
Proportion Explained 0.47 0.37 0.16
Cumulative Proportion 0.47 0.84 1.00
Mean item complexity = 1.5
Test of the hypothesis that 3 factors are sufficient.
The degrees of freedom for the null model are 300 and the objective function was 62.46 with Chi Square of 2488.15
The degrees of freedom for the model are 228 and the objective function was 29.15
The root mean square of the residuals (RMSR) is 0.07
The df corrected root mean square of the residuals is 0.08
The harmonic number of observations is 50 with the empirical chi square 152.35 with prob < 1
The total number of observations was 50 with Likelihood Chi Square = 1102.82 with prob < 1.2e-114
Tucker Lewis Index of factoring reliability = 0.442
RMSEA index = 0.276 and the 90 % confidence intervals are 0.263 0.296
BIC = 210.88
Fit based upon off diagonal values = 0.98
Measures of factortscore adequacy
ML1 ML2 ML3
Correlation of (regression) scores with factors 1.00 0.99 0.98
Multiple R square of scores with factors 1.00 0.98 0.96
Minimum correlation of possible factor scores 0.99 0.97 0.93
item 変数のもとの順序(順番がかわているため)
ML1 ~ ML3 共通因子ごとの因子負荷量(共通因数に対する各変数の寄与の度合い)
h2 共通性(各変数が共通因子で説明できる度合い 1-u2)
u2 独自性(各変数が共通因子で説明できない度合い(大きいと望ましくない))
com 複雑性(各変数の因子負荷量が複数の共通因子にまたがる度合い)
まずはML1とML2の絶対値を見比べ、ML1が大きい値を確認する。"小売業販売額事業所あたり百万円"まではML1が大きい。これより上の項目はML1に寄与する項目であり、これより下はML2に寄与する項目である。またML2とML3を比べると、"ごみリサイクル率_pct"まではML2の方が大きい。
ML1へのへの因子負荷量が高い項目{千人あたり事業所数,昼間人口比per,千人あたり大型小売店数,千人あたり飲食店数など}は商業の盛ん度を示していると考えられるため、ML1をビジネスとする。ML2の項目を見ると{千人あたり老人ホーム数(-),高齢単身世帯比率,年齢15未満比率(-), 世帯あたり人数(-), 第3次産業従業者数比}が特徴的であり、単身者が多そうであると予想できる。また三次産業従業者数比が多いことからホワイトカラーの職種が多そうである、以上よりML2を都会生活とする。最後にML3であるが{年齢65以上比率(-), 国民健保一人あたり診療費円(-)}の2項目である。これは高齢者が少ない社会を示しており、非高齢化とする。
次に確認するところとしては"Cumulative Var"であり、今回の例では第3因子までを抽出すると、全体の69%を説明できていることを示す。
#因子負荷量をプロット
#因子負荷量はresultFAのloadingsに格納されている
par(family = "HiraKakuProN-W3")
plot(resultFA$loadings[,1], #ビジネス
resultFA$loadings[,2], #都会生活
type="n") #tyoe="n"でプロットをしない
#枠内にテキストを表示
text(resultFA$loadings[,1],
resultFA$loadings[,2],
rownames(resultFA$loadings), col="steelblue")
#y=0の直線を引く
lines(c(-1,1), c(0,0), col="gray")
#x=0の直線を引く
lines(c(0,0), c(-1,1), col="gray")
上の散布図は各変数の因子負荷量をプロットすることで、どの変数がどの共通因子に寄与しているかを可視化している。軸を変えることで他の変数についても調べることができる。
#因子得点を基にクラスタリング
kmFA <- kmeans(DFfa, 4, iter.max=40)
#色ラベルの配列のために番号をコピー
color_kmFA <- kmFA$cluster
head(color_kmFA)
#factorを文字列に
color_kmFA <- as.factor(color_kmFA)
head(color_kmFA)
#クラスタ番号に色名
cokor_kmFA <- as.factor(color_kmFA)
levels(color_kmFA) <- c("red","orange", "blue", "green")
head(color_kmFA)
#factorを文字列に
color_kmFA <- as.character(color_kmFA)
#因子得点の値を色分け
#ラベルをrownames()
plot(DFfa$ビジネス,
DFfa$都会生活,
type="n") #tyoe="n"でプロットをしない
#枠内にテキストを表示
text(DFfa$ビジネス,
DFfa$都会生活,
rownames(DFfa), col=color_kmFA)
lines(c(-10,10), c(0,0), col="gray")
lines(c(0,0), c(-10,10), col="gray")
これにより市区町ごとのクラスタリングを可視化できる。図を見ると千代田区はビジネスに特化しているように思われる。また都会生活度の低い地域(オレンジ)は東京の西側に多く存在している。
主成分分析の実装
Rで主成分分析を行う際はprcomp()関数を使う。またscale=TRUEでで標準化を行うことができる。
resultPCA <- prcomp(DF[, -(1:3)], scale=TRUE)
summary(resultPCA)
>summary(resultPCA)
Importance of components:
PC1 PC2
Standard deviation 3.3971 2.1271
Proportion of Variance 0.4616 0.1810
Cumulative Proportion 0.4616 0.6426
#各変数ごとの主成分(固有ベクトル)
#第3主成分まで表示
resultPCA$rotation[, 1:3]
#各ケースごとの主成分得点
resultPCA$x[1:5, 1:3]
#主成分得点をプロット
biplot(resultPCA)
summary()のCumulative Proportionを見ると第2因子までで、64%の分散が説明できる。bioplotで主成分の値を出力することができる。また第1主成分と第2主成分を軸に自治体の番号を出力できる。
#主成分得点をベースにクラスタリング
#第3主成分までを使う
kmPCA <- kmeans(DFpca[, 1:2], 4, iter.max = 50)
color_kmPCA <- kmPCA$cluster
#色ラベル
color_kmPCA <- as.factor(color_kmPCA)
levels(color_kmPCA) <- c("red","orange", "blue", "green")
head(color_kmPCA)
color_kmPCA <- as.character(color_kmPCA)
head(color_kmPCA)
#主成分分析をプロット
plot(DFpca$PC1, DFpca$PC2, type="n")
text(DFpca$PC1, DFpca$PC2, rownames(DFpca), col=color_kmPCA)
lines(c(-20, 20), c(0,0), col="grey")
lines(c(0,0), c(-20,20), col="grey")
因子分析の時と同様に今回も主成分分析をもとにクラスタリングを行う。
