実行環境

Mac OSX Mojave
jupyter notebook
Python 3.6.8
statsmodels 0.9.0
R 3.6.0

参考図書

オーム社入門統計学
ファーストブック　多変量解析がわかる

分からないこと・やりたいこと

Python の statsmodels.multivariate.pca.PCA を使って主成分分析をしたときに、主成分負荷量を算出したいのですが、固有ベクトルと主成分負荷量の意味がこんがらがってしまっており、正しい結果を出力するにはどうすればいいのか分からなくなってきました。

参考図書『多変量解析がわかる』P.78には、

p=ax+by+cu+dv+ew
ただし、a^2+b^2+c^2+d^2+e^2=1

合成変量pが最大の分散値を持つとき、そのpを主成分と呼びます。このときの定数a,b,・・・,eの値を主成分負荷量と呼びます。この主成分pを用いてデータを分析するのが主成分分析です。

と書かれています。
そのため、主成分負荷量は、主成分の式の各係数を二乗した合計値は1になるものと思っていました。

知りたいこと

主成分負荷量とは、上記のように各係数を二乗した合計値は1になるものなのでしょうか。それとも、それは固有ベクトルなのでしょうか。
statsmodels.multivariate.pca.PCA を使って主成分負荷量を出力する方法を知りたいです。

試したこと・参考にしたこと・そしてよく分からなくなったこと

ファーストブック　多変量解析がわかる　3-2 にある主成分分析の例題をもとに実施

出席番号	数学 x	理科 y	社会 u	英語 v	国語 w
1	71	64	83	100	71
2	34	48	67	57	68
3	58	59	78	87	66
4	41	51	70	60	72
5	69	56	74	81	66
6	64	65	82	100	71
7	16	45	63	7	59
8	59	59	78	59	62
9	57	54	84	73	72
10	46	54	71	43	62
11	23	49	64	33	70
12	39	48	71	29	66
13	46	55	68	42	61
14	52	56	82	67	60
15	39	53	78	52	72
16	23	43	63	35	59
17	37	45	67	39	70
18	52	51	74	65	69
19	63	56	79	91	70
20	39	49	73	64	60

この表を元に主成分分析をする。

python
1import numpy as np
2from statsmodels.multivariate.pca import PCA
3
4data = np.array([
5    [71,34,58,41,69,64,16,59,57,46,23,39,46,52,39,23,37,52,63,39],
6    [64,48,59,51,56,65,45,59,54,54,49,48,55,56,53,43,45,51,56,49],
7    [83,67,78,70,74,82,63,78,84,71,64,71,68,82,78,63,67,74,79,73],
8    [100,57,87,60,81,100,7,59,73,43,33,29,42,67,52,35,39,65,91,64],
9    [71,68,66,72,66,71,59,62,72,62,70,66,61,60,72,59,70,69,70,60]
10])
11
12ret = PCA(data.T, ncomp=1, standardize=False, normalize=False)
13
14print("主成分", 100 - ret.factors)
15print("主成分負荷量", ret.coeff*-1)

主成分 [[149.96661766]
 [ 90.07331542]
 [130.6235937 ]
 [ 97.38181772]
 [129.76064147]
 [146.50319377]
 [ 37.90790262]
 [107.65862887]
 [119.26918818]
 [ 85.78744696]
 [ 64.52332319]
 [ 69.9967705 ]
 [ 84.47561479]
 [110.96612123]
 [ 91.68946286]
 [ 64.18121726]
 [ 76.26684286]
 [107.50284536]
 [136.34679483]
 [ 99.11866074]]
主成分負荷量 [[0.491513   0.1730531  0.19588597 0.82781535 0.06941201]]

それからこの主成分負荷量は、それぞれ二乗して合計すると1となりました。

python
1sum(ret.coeff[0]**2)

1.0000000000000002

オーム社入門統計学の例題 14章を元に主成分分析を実施

ところが、別の参考書『オーム社入門統計学』の例題 14章をやったときに分からなくなりました。

集落番号	農家数の増減率	農家人口増減率	農業就業人口率	農業従事者率	男子農業就業人口率	農業就業人口のうち生産年齢人口率	60歳未満男子農業専従者率	農業従事者のうちの基幹的農業従事者率	兼業従事者率
1	14.3	-9.7	50.0	73.1	38.5	23.1	0.0	47.4	38.5
2	-14.3	-21.7	55.6	77.8	62.5	30.0	0.0	50.0	38.9
3	-17.6	-22.7	60.0	89.1	50.0	48.5	10.2	51.0	40.0
4	-8.3	3.8	58.7	76.1	44.4	51.9	8.6	57.1	32.6
5	-3.4	-10.6	61.9	85.6	66.7	38.3	6.0	56.6	37.1
6	0.0	-7.5	63.6	87.9	60.0	47.6	13.8	37.9	36.4
7	-20.0	-38.2	45.2	88.1	38.1	52.6	0.0	27.0	52.4
8	-15.0	-8.3	45.2	77.4	39.3	42.9	2.1	45.8	37.1
9	0.0	-21.4	61.9	71.4	63.6	30.8	0.0	73.3	28.6
10	-11.1	-16.7	54.2	79.2	35.7	7.7	0.0	57.9	33.3
11	-20.0	-18.0	42.2	77.8	45.5	36.8	8.6	48.6	55.6
12	-7.7	-24.0	57.9	81.6	55.0	27.3	3.2	64.5	34.2
13	-11.8	-28.4	64.7	92.2	50.0	45.5	2.1	55.3	37.3
14	0.0	15.4	44.4	92.6	38.5	41.7	0.0	24.0	59.3
15	-33.3	-51.3	42.1	89.5	22.2	25.0	0.0	35.3	52.6

データは、14-1.xlsxより

この表を元に主成分分析を実施しました。

python
1import pandas as pd
2from statsmodels.multivariate.pca import PCA
3
4df = pd.read_excel('14-1.xlsx')
5df = df.set_index("集落番号")
6
7ret = PCA(df, ncomp=2, standardize=False, normalize=False)
8
9display(ret.factors)
10display(ret.coeff)

出力結果

集落番号	comp_0	comp_1
1	-9.588784	5.633295
2	-5.765015	-9.884144
3	1.972712	-2.103843
4	-17.528468	14.211070
5	-21.226152	0.077743
6	-12.981545	16.437488
7	35.192034	4.361661
8	3.799128	8.998005
9	-29.143404	-21.182749
10	-2.305750	-16.643283
11	12.697244	2.040215
12	-12.495731	-18.045632
13	-1.230358	-9.314251
14	8.327407	45.738271
15	50.276684	-20.323845

|農家数の増減率|農家人口増減率|農業就業人口率|農業従事者率|男子農業就業人口率|農業就業人口のうち生産年齢人口率|60歳未満男子農業専従者率|農業従事者のうちの基幹的農業従事者率|兼業従事者率|
|------------|-----------|---|----|----|----|----|----|
|comp_0|-0.398821|-0.432523|-0.295694|0.163024|-0.457945|0.021984|-0.059636|-0.454728|0.345431|
|comp_1|0.226974|0.697246|-0.116084|0.110778|-0.082904|0.378380|0.065972|-0.478187|0.231353|

この出力結果のうち comp_0　を二乗して合計してみると

python
1rint(sum(np.array(ret.coeff.T['comp_0'])**2))

1.0000000000000004

正しい結果となっているようですが、
参考図書の解答結果と異なるのです。

参考図書の解答結果は、エクセル統計2010で算出したとのことですが、次のような結果となっていました。

主成分負荷量

変数	第１主成分	第２主成分
農家数の増減率	0.61	0.05
農家人口増減率	0.33	0.35
農業就業人口率	0.83	0.28
農業従事者率	-0.51	0.58
男子農業就業人口率	0.78	0.31
農業就業人口のうち生産年齢人口率	-0.08	0.85
60歳未満男子農業専従者率	0.29	0.73
農業従事者のうちの基幹的農業従事者率	0.82	-0.37
兼業従事者率	-0.89	0.20

参考までに第１主成分の各値を二乗して合計してみました。

python
1t1 = np.array([0.61,0.33,0.83,-0.51,0.78,-0.08,0.29,0.82,-0.89])
2print(sum(t1**2))

3.5934000000000004

値が１になりません。。。

この辺りからこんがらがってきました。
参考図書やソフトウェアによって主成分負荷量は異なってくるのか。
そしたら一体どれが正しいのか。

そこで、Rでも試してみました。
14-1.csv は先ほどの表から集落番号のみ削ったcsvです。

R
1> data = read.csv('14-1.csv')
2> pca = prcomp(data, scale=T)
3> pca$rotation[,1:2]

                                             PC1         PC2
農家数の増減率                       -0.32022098  0.03432674
農家人口増減率                       -0.17526318  0.24398192
農業就業人口率                       -0.44142820  0.19710407
農業従事者率                          0.26724962  0.40110345
男子農業就業人口率                   -0.41161307  0.21355265
農業就業人口のうち生産年齢人口率      0.04087309  0.59333011
X60歳未満男子農業専従者率            -0.15337729  0.50718952
農業従事者のうちの基幹的農業従事者率 -0.43172740 -0.25672953
兼業従事者率                          0.46845603  0.13698137

この値は、固有ベクトルになるとRで主成分分析を説明してあるサイトにかかれていました。
また、この固有ベクトルから主成分負荷量を求めるには、固有値の正の平方根を掛け合わせれば良いとも。

R
1> t(t(pca$rotation)*pca$sdev)[,1:2]

                                             PC1         PC2
農家数の増減率                       -0.60524876  0.04945746
農家人口増減率                       -0.33126443  0.35152552
農業就業人口率                       -0.83434219  0.28398460
農業従事者率                          0.50512775  0.57790386
男子農業就業人口率                   -0.77798869  0.30768347
農業就業人口のうち生産年齢人口率      0.07725411  0.85486117
X60歳未満男子農業専従者率            -0.28989798  0.73075109
農業従事者のうちの基幹的農業従事者率 -0.81600673 -0.36989208
兼業従事者率                          0.88542741  0.19736072

これによって得た値は、参考図書に掲載されているエクセル統計2010で算出した主計分負荷量とほぼ同じ値になりました。
（符号は違ったりしますが）

というわけで、pythonのstatsmodelsで同じことをやってみょうとしたのですが、

python
1import pandas as pd
2from statsmodels.multivariate.pca import PCA
3
4df = pd.read_excel('14-1.xlsx')
5df = df.set_index("集落番号")
6
7ret = PCA(df, ncomp=2, standardize=True, normalize=False)
8
9display(ret.coeff)

出力結果

|農家数の増減率|農家人口増減率|農業就業人口率|農業従事者率|男子農業就業人口率|農業就業人口のうち生産年齢人口率|60歳未満男子農業専従者率|農業従事者のうちの基幹的農業従事者率|兼業従事者率|
|------------|--|--|-|-|-|-|
|comp_0|-0.320221|-0.175263|-0.441428|0.267250|-0.411613|0.040873|-0.153377|-0.431727|0.468456|
|comp_1|0.034327|0.243982|0.197104|0.401103|0.213553|0.593330|0.507190|-0.256730|0.136981|

PCAのstandarizeをTrueに変えて、固有ベクトルの値を同じにすることができたのですが、
そのあとの固有値の正の平方根を掛け合わせることができません。（わかりません）

このように、参考図書によって違ってきたりするので、何が正しいことなのかよく分からなくなってしまいました。

どうかよろしくお願いします。

行動規範の内容に同意します

回答1件

ベストアンサー

まず、固有ベクトルは理屈の上では長さがいくつであってもいいのです。固有ベクトルですから。

実際は長さ不定、だと不便だし、数値的に計算するときも解が求まらなかったりします。なので、ノルムが1になるようにする、ということはよくやります。

ここでいうノルムは二乗和の平方根ですが、結果的に1になるように正規化されたのであればただの二乗和を計算しても1になります。

また、主成分負荷量は各変数と主成分の相関係数で、すべての変数に対して計算してノルムが1になる保証はどこにもありません。そう書いてあったなら本の間違いです（とはいえ、この説明もよく見るので、もしかしたらなにか根拠があってそう言っているのかもしれません。ただいずれにせよ、そう理解することはおすすめしません）。

計算方法については、こちらの10ページなどがわかりやすいかと思います。固有ベクトルに対応する固有値をかけるというのは、入力の標準化をした場合に使える便法です。

http://manabukano.brilliant-future.net/document/text-PCA.pdf

なお、主成分負荷量は知りたければ出しておいても良いんじゃないかな程度であって、固有ベクトルだけで解釈してもそんなに困ることはないと思います。状況によりけりですが……

sklearnでやる場合

素直に上の方法で計算可能です。

python
1import numpy as np
2from sklearn.decomposition import PCA as PCA
3from sklearn.preprocessing import StandardScaler
4
5data = np.array([
6    [71,34,58,41,69,64,16,59,57,46,23,39,46,52,39,23,37,52,63,39],
7    [64,48,59,51,56,65,45,59,54,54,49,48,55,56,53,43,45,51,56,49],
8    [83,67,78,70,74,82,63,78,84,71,64,71,68,82,78,63,67,74,79,73],
9    [100,57,87,60,81,100,7,59,73,43,33,29,42,67,52,35,39,65,91,64],
10    [71,68,66,72,66,71,59,62,72,62,70,66,61,60,72,59,70,69,70,60]
11])
12
13pca = PCA(n_components=1, svd_solver="full")
14X = pca.fit_transform(StandardScaler().fit_transform(data.T))
15R = []
16for row in data:
17    R.append(np.corrcoef(X.ravel(), row)[0,1])
18print(np.array(R))
19print(pca.components_[0]*np.sqrt(pca.explained_variance_[0]))
20""" =>
21[-0.94621924 -0.91293882 -0.91759568 -0.94095263 -0.48408647]
22[-0.97080045 -0.93665547 -0.94143331 -0.96539702 -0.49666224]
23"""

少しずれますが、原因は究明していません（自由度の調整かなにかが怪しいような気がする。ただの浮動小数点数の誤差かもしれないけど）。

statsmodelsでやる場合

eigenvalsで取れるかなとも思ったけど、なんだか信じがたい数字が出たので保留。

投稿2019/06/22 05:53

hayataka2049

総合スコア30933

waku_nagoya

2019/06/22 09:16

ご解答ありがとうございます！固有ベクトルさえあれば、主成分負荷量はどちらでも良い程度のものなのですね。まずは教えていただいた内容を元にもう少し調べさせてください。

hayataka2049

2019/06/22 09:17

statsmodelsはあまり詳しくないんですが、整合性のある結果が出ない気がします。こっちでももう少し検討してみます。

waku_nagoya

2019/06/23 23:43

ありがとうございます。まずは、この手の話題を話せる人が身近にいないため、ご相談できて嬉しいです。いただいたsklearnのソースを元に自分なりに落とし込んでみました。 sklearnだとcomponents_によって固有ベクトルとexplained_variance_で固有値が取得できるので、それを掛け合わせると求められるというわけですね。私も調べていた時にこの方法はちらっと試してみたことがありました。近い値はでるのですが、若干ずれるのでなんでだろうと思って、そのままにしておりました。 hayataka2049さんの相関係数の方法は、参考書の値とぴったり一致しました。あと気になるのは、skleanだと機械学習からのアプローチになりますが、統計学からだと、statsmodelsやscipi.statsあたりから算出するといいのかなと思っていました。とはいっても、両方とも、多変量解析になってくるとどうしても機能的に足りない部分がでてくるので、この辺りは別の方法に頼るしかないのか、それとも計算式を理解して、関数を作るしかないのかと悩んでいたりもしました。 statsmodelsのeigenvalsなど、確かに本当に合っているのか分からないような値、固有ベクトルがそのままでているのではないかという点もあって、どこまで使えるのかが怪しいです。ソースは公開されているのですが、何分、計算式の方の知識が足りなくて解読できずにいました。。。

行動規範の内容に同意します

あなたの回答

tips

プレビュー

行動規範の内容に同意します

質問の解決につながる回答をしましょう。サンプルコードなど、より具体的な説明があると質問者の理解の助けになります。また、読む側のことを考えた、分かりやすい文章を心がけましょう。

15分調べてもわからないことは
teratailで質問しよう！

ただいまの回答率
85.48%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう！

statsmodelsによる主成分分析での主成分負荷量の算出の方法

実行環境

参考図書

分からないこと・やりたいこと

知りたいこと

試したこと・参考にしたこと・そしてよく分からなくなったこと

ファーストブック　多変量解析がわかる　3-2 にある主成分分析の例題をもとに実施

オーム社入門統計学の例題 14章を元に主成分分析を実施

sklearnでやる場合

statsmodelsでやる場合

関連した質問

出席番号	数学 x	理科 y	社会 u	英語 v	国語 w
1	71	64	83	100	71
2	34	48	67	57	68
3	58	59	78	87	66
4	41	51	70	60	72
5	69	56	74	81	66
6	64	65	82	100	71
7	16	45	63	7	59
8	59	59	78	59	62
9	57	54	84	73	72
10	46	54	71	43	62
11	23	49	64	33	70
12	39	48	71	29	66
13	46	55	68	42	61
14	52	56	82	67	60
15	39	53	78	52	72
16	23	43	63	35	59
17	37	45	67	39	70
18	52	51	74	65	69
19	63	56	79	91	70
20	39	49	73	64	60

出席番号	数学 x	理科 y	社会 u	英語 v	国語 w
1	71	64	83	100	71
2	34	48	67	57	68
3	58	59	78	87	66
4	41	51	70	60	72
5	69	56	74	81	66
6	64	65	82	100	71
7	16	45	63	7	59
8	59	59	78	59	62
9	57	54	84	73	72
10	46	54	71	43	62
11	23	49	64	33	70
12	39	48	71	29	66
13	46	55	68	42	61
14	52	56	82	67	60
15	39	53	78	52	72
16	23	43	63	35	59
17	37	45	67	39	70
18	52	51	74	65	69
19	63	56	79	91	70
20	39	49	73	64	60

質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう！

実行環境

参考図書

分からないこと・やりたいこと

知りたいこと

試したこと・参考にしたこと・そしてよく分からなくなったこと

ファーストブック 多変量解析がわかる 3-2 にある主成分分析の例題をもとに実施

オーム社 入門統計学 の例題 14章 を元に主成分分析を実施

sklearnでやる場合

statsmodelsでやる場合

関連した質問

ファーストブック　多変量解析がわかる　3-2 にある主成分分析の例題をもとに実施

オーム社入門統計学の例題 14章を元に主成分分析を実施

出席番号	数学 x	理科 y	社会 u	英語 v	国語 w
1	71	64	83	100	71
2	34	48	67	57	68
3	58	59	78	87	66
4	41	51	70	60	72
5	69	56	74	81	66
6	64	65	82	100	71
7	16	45	63	7	59
8	59	59	78	59	62
9	57	54	84	73	72
10	46	54	71	43	62
11	23	49	64	33	70
12	39	48	71	29	66
13	46	55	68	42	61
14	52	56	82	67	60
15	39	53	78	52	72
16	23	43	63	35	59
17	37	45	67	39	70
18	52	51	74	65	69
19	63	56	79	91	70
20	39	49	73	64	60