主成分分析は、データの変数が多い時に、出来るだけ情報を減らさずに、次元を圧縮するテクニックとして用いられます。マーケティングにおいても、何かの商品に対する評価項目が複数ある場合に、それを少数の評価項目に圧縮し、総合力のような観点で評価することもできるでしょう。
今回はそんな主成分分析を使って、様々な食品の栄養素に対して主成分分析を行い、総合的に色々な栄養素を含んでいる、つまり、忙しいあなたに、これ一食でも大丈夫という食品を見つけたいと思います。
後半ではRとPythonによる主成分分析実践法も紹介します。
主成分分析とは
ところでデータを主成分に圧縮するとはどのようなことでしょうか？例えば、効率よく一食でタンパク質と炭水化物の両方を摂取できる食品を選びたいとしましょう。筋力トレーニングが大好きな人はまさに理想的な食品ですね。
まずは仮想的に作成した全２００品目の食品に含まれるタンパク質を横軸、炭水化物を縦軸に、散布図を書いてみます。
するとこの緑矢印の先の方にある赤丸で囲った食品がタンパク質と炭水化物の両方を多めに含んでいる食品のようです。つまり、この緑矢印はタンパク質と炭水化物の両方の情報を圧縮した新たな変数、第一主成分となり、この緑矢印を見つけるのが主成分分析なのです。
さて、タンパク質と炭水化物情報を一つの変数で圧縮してしまったために、失われてしまった情報もあります。緑矢印と直角になる矢印、この青色矢印の情報ですね。
この方向のばらつきは緑矢印では説明できませんが、タンパク質を多めに含んでいるか、炭水化物を多めに含んでいるのかという情報が含まれていそうです。これを第二主成分とします。
この第一主成分をx軸に、第二主成分をy軸にしてグラフを書き直してみます。
主成分分析のメカニズムの関係で、値が平均0、分散1に標準化されていますが、関係性としては変わりません。ただし、二つの変数から二つの成分を抽出するということは、全く情報を圧縮していません。ただ単に、データの座標を回転させただけ、つまり見かたを変えただけということです。もう一度強調しておくと、主成分分析はデータ情報を主成分に圧縮することが目的であるということを忘れないようにしましょう。
主成分の求め方
ここで、主成分とはどのような意味を持つのでしょうか？データを圧縮できるということは、主成分がデータの情報を説明できる成分であるということですが、
例えば、タンパク質と炭水化物の関係性がこのようであったら、主成分はどうなるでしょうか？
このデータのバラつきは緑矢印の主成分だけでほとんど説明できそうですね。先ほどとは違い、青矢印方向にはほとんどバラツキはなさそうです。
このように変数ごとのバラツキを一番よく説明できるように、情報を圧縮できる成分が主成分となります。バラツキを使うということは、変数同士でばらつきが大きく違っては、うまく主成分を求めることができません。（ばらつきの大きな変数に引っ張られてしまう。）そこで、平均値を0、分散を1に標準化する必要があったのです。具体的にどう求めるのか気になる方は、数学的な説明な説明が詳しく載っている参考サイトのリンクがページ下部にありますので、そちらをご覧ください。
主成分得点と因子負荷量について
最後にもう一つ主成分得点と因子負荷量についてご紹介します。まず、因子負荷量とは、変数と主成分との関係を表し、-1~1の値を取ります。変数と主成分との相関係数とも言ってもよいでしょう。先ほどのグラフでは、炭水化物やタンパク質の値が大きくなるほど第一主成分も大きくなりました。実際に因子負荷量を求めてみると
- 炭水化物：0.7
- タンパク質：0.7
となります。また第二主成分の因子負荷量は
- 炭水化物：0.7
- タンパク質：-0.7
となります。第一主成分はどちらも正の値、第二主成分はタンパク質が負の値になっているので、グラフの見た目や相関係数の概念とも一致しています。
次に主成分得点とは、個々のサンプルが主成分の要因をどれぐらい持っているのかを表します。つまりは、主成分を横軸にした先ほどのグラフでは、右や左に行くほど第一主成分の要因を多くもつ食品ということになり、この主成分の大きさのことを主成分得点と言います。主成分得点はどのように求めるかというと、因子負荷量と各サンプルの値とを掛け算することで求めることができます。つまり
です。例えば、プラス方向で第一主成分得点の一番高い食品を例にすると、この食品の標準化したタンパク質と炭水化物の量は
- タンパク質：2.4
- 炭水化物： 1.4
です。これにそれぞれの因子負荷量をかけ合わせて足すと
となり、これがこの食品の主成分得点となります。先ほどの主成分を縦軸、横軸にしたグラフでも確認できます。
完全食を探して
さて、主成分分析の心を理解したところで、実際に完全食を探してみましょう。今回扱うデータは、全8618種類の食品栄養素がまとめられたデータです。
データはこちらから入手できます。
中身を確認してみましょう。
| ID | Energy_kcal | Protein_g | Fat_g | Carb_g |
|---|---|---|---|---|
| 01001 | 717 | 0.85 | 81.11 | 0.06 |
| 01002 | 717 | 0.85 | 81.11 | 0.06 |
| 01003 | 876 | 0.28 | 99.48 | 0.00 |
| 01004 | 353 | 21.40 | 28.74 | 2.34 |
| 01005 | 371 | 23.24 | 29.68 | 2.79 |
| 01006 | 334 | 20.75 | 27.68 | 0.45 |
それぞれの行がある一食品に対応し、様々な栄養素（ビタミン、ミネラルなど）全 22種類の情報が含まれています。ちなみに、
| ID | FoodGroup | Descrip |
|---|---|---|
| 01001 | Dairy and Egg Products | Butter, salted |
| 01002 | Dairy and Egg Products | Butter, whipped, with salt |
| 01003 | Dairy and Egg Products | Butter oil, anhydrous |
| 01004 | Dairy and Egg Products | Cheese, blue |
| 01005 | Dairy and Egg Products | Cheese, brick |
| 01006 | Dairy and Egg Products | Cheese, brie |
その食品の種類や、簡単な説明などの情報を見ることもできます。
さっそく、全ての栄養素を網羅的に含んでいる食品、"完全食"、をこの中から探していきます。つまり、第一主成分との関連が一番大きい食品ですね。ここからはRとPythonで記述が異なるので分けて書いていきます。R使いの方は以下から、Python使いの方はこちらです！
R使いの方はこちら
まずはデータの読み込みです。USRDAと名のつく変数はアメリカの推奨摂取量であり、本解析の趣旨とは関係ないので除きます。
library(tidyverse) library(ggfortify)#可視化を楽にするためのパッケージ food <- read_csv("nndb_flat.csv") food <-food %>% select(-matches("_USRDA"))
次に栄養素の分布を確認してみます。カロリーは栄養素ではなく、ある意味で変数の集約情報であるため分析からは外しました。ちなみにコードで出力されるヒストグラムと以下のヒストグラム は画像の関係で異なっています。
#栄養素の情報のみ抜き出す。 food_nut <- food[,-c(1,2,3,4,5,6,7,8)] food_nut %>% gather(key="Name",value="Nutrition") %>% ggplot()+ facet_wrap(.~Name)+ geom_histogram(aes(x=Nutrition),bins=50,fill="green")+ theme(axis.text = element_blank(), axis.title = element_blank(), axis.ticks = element_blank(), strip.text = element_text(size=10))
かなり右裾の長い分布になっていますね。分布を調整するため対数変換を行います。本来であればBox-Cox変換のような高度な変換を行う方が良いのですが、今回は主成分分析がメインということで簡単な対数変換にしました。（0の値には微小な値を足して対数変換しています。）
#ログ変換 food_nut_log <- food_nut %>% map_df(~log(.+0.0000001))
0が多かったためか左端が気になりますが、先ほどよりは分布が整いました。このデータに対して主成分分析を行います。Rで主成分分析を行うにはprcomp{stats}を使います。prcompの引数であるcenterとscaleは標準化するかどうかの指定です。
pc <- prcomp(food_nut_log,center = TRUE, scale.= TRUE) Standard deviations (1, .., p=22): [1] 3.0063931 1.7042684 1.3105218 1.0978138 1.0465825 0.9707619 0.8440906 0.7702426 [9] 0.7095602 0.6951428 0.6614608 0.6120544 0.5830755 0.5541138 0.5446296 0.4749278 [17] 0.4503909 0.4120321 0.3809755 0.3597714 0.3104601 0.2816171 Rotation (n x k) = (22 x 22): PC1 PC2 PC3 PC4 PC5 Protein_g 0.19630383 -0.041307363 -0.383389327 0.013522247 -0.346330236 Fat_g 0.12242931 -0.112798763 -0.397037212 0.322553636 -0.205278110 Carb_g -0.04777346 0.502810434 -0.073820020 -0.037369785 0.000623154 Sugar_g -0.05286341 0.442373501 -0.092043487 0.345584348 0.102885617 Fiber_g -0.01363876 0.493778068 -0.172950526 0.014269967 0.053453524
結果をみやすくするため始めの方だけ出力しましたが、22この栄養素全てで結果が出力されており、22この主成分が抽出されています(PC22まである)。主成分はデータのバラツキを100%説明できるまで抽出されるため、変数の数だけ存在することを覚えておきましょう。
結果の見方ですが、まず上のStandard deviationsは、その主成分がどれぐらいデータの標準偏差を説明しているかを表しています。今回は22この変数を全て標準化し分散を1にしたので、これを二乗して全て足すと22になります。
> sum(pc$sdev^2) [1] 22
つまり、第1主成分は22のうち9（3の二乗）の分散を説明できることになるので、9/22=0.41で、41%の情報を持っていることになります。それぞれの成分がどれぐらいの情報を持っているか、縦軸を成分が持っている分散、横軸を成分としてグラフに表してみます。
大体第4主成分ぐらいまでで、全データの7割ぐらいを説明することができるようです。もし、このデータで他の解析をする場合は、最小で4つの変数まで次元を圧縮できそうです。
次に結果のRotationとあるのは、それぞれの栄養素が、主成分からどれほどの影響を受けているのか、つまり因子負荷量です。
先ほど同様にデータを横軸を第一主成分で、縦軸を第二主成分として、グラフにしてみます。可視化にはautplot{ggfortify}を使います。
autoplot(pc,loadings = TRUE, loadings.label = TRUE)
赤の矢印はそれぞれの栄養素の因子負荷量を表し、多くの栄養素は第一主成分で説明できそうですが、炭水化物系は方向が違いますね（ちなみに食物繊維と砂糖を合わせたものが炭水化物です）。代わりに第二主成分では炭水化物系と強く関連していそうなので、完全食を探すには第二主成分も考慮した方が良いかもしれません。
Python使いの方はこちら
まずはデータの読み込みです。USRDAと名のつく変数はアメリカの推奨摂取量であり、本解析の趣旨とは関係ないので除きます。
import pandas as pd import numpy as np from sklearn.decomposition import PCA from sklearn.preprocessing import StandardScaler import matplotlib.pyplot as plt food <- read_csv("nndb_flat.csv") nut = pd.read_csv("../nndb_flat.csv") nut.drop(nut.columns[nut.columns.str.contains("_USRDA")],axis=1,inplace=True)
次には栄養素の分布を確認してみます。カロリーは栄養素ではなく、ある意味で変数の集約情報であるため分析からは外しました。
#栄養素の情報のみ抜き出す。 food_nut = food.iloc[:,8:] food_nut.hist(bins=50, xlabelsize=-1, ylabelsize=-1, figsize=(11,11)) plt.show()
かなり右裾の長い分布になっていますね。分布を調整するため対数変換を行います。本来であればBox-Cox変換のような高度な変換を行う方が良いのですが、今回は主成分分析がメインということで簡単な対数変換にしました。（0の値には微小な値を足して対数変換しています。）
#ログ変換 f = lambda x: np.log(x+0.0000001) food_nut_log = food_nut.applymap(f) food_nut_log.hist(bins=50, xlabelsize=-1, ylabelsize=-1, figsize=(11,11)) plt.show()
0が多かったためか左端が気になりますが、先ほどよりは分布が整いました。このデータに対して主成分分析を行います。Pythonでは最初にデータを標準化した後に、主成分分析を行います。sklearnを使えばどちらも簡単に実行可能です。
fit_sc = StandardScaler()
food_nut_log_sc = fit_sc.fit_transform(food_nut_log)
fit_pca = PCA()
pca = fit_pca.fit_transform(food_nut_log_sc)
pd.DataFrame({"Variace" : fit_pca.explained_variance_,
"Variace_ratio " : fit_pca.explained_variance_ratio_},
index=["PC{}".format(x) for x in range(1,23)]).iloc[:5,]
pd.DataFrame(fit_pca.components_,columns=food_nut.columns,index=["PC{}".format(x) for x in range(1,23) ]).iloc[:5,:4]
| Variace | Variace_ratio | |
|---|---|---|
| PC1 | 9.039449 | 0.410836 |
| PC2 | 2.904868 | 0.132024 |
| PC3 | 1.717667 | 0.078067 |
| PC4 | 1.205335 | 0.054782 |
| PC5 | 1.095462 | 0.049788 |
| Protein_g | Fat_g | Carb_g | Fiber_g | |
|---|---|---|---|---|
| PC1 | -0.196304 | -0.122429 | 0.047773 | 0.052863 |
| PC2 | 0.041307 | 0.112799 | -0.502810 | -0.442374 |
| PC3 | -0.383389 | -0.397037 | -0.073820 | -0.092043 |
| PC4 | -0.013522 | -0.322554 | 0.037370 | -0.345584 |
| PC5 | -0.346330 | -0.205278 | 0.000623 | 0.102886 |
結果をみやすくするため始めの方だけ出力しましたが、22この栄養素全てで結果が出力されており、22この主成分が抽出されています(PC22まである)。主成分はデータのバラツキを100%説明できるまで抽出されるため、変数の数だけ存在することを覚えておきましょう。
結果の見方ですが、まずVarianceとVariance_ratioは、その主成分がどれぐらいデータの分散を説明しているかを表しています。今回は22この変数を全て標準化し分散を1にしたので、これを全て足すと22になります。
>fit_pca.explained_variance_.sum()
22.00255309272368
つまり、第1主成分は22のうち9の分散を説明できることになるので、9/22=0.41で、41%の情報を持っていることになります。それぞれの成分がどれぐらいの情報を持っているか、縦軸を成分が持っている分散、横軸を成分としてグラフに表してみます。
plt.figure(figsize=(8, 6)) plt.plot(range(1,23),fit_pca.explained_variance_) plt.tick_params(labelsize=15) plt.grid() for x ,y , prop in zip(range(1,23),fit_pca.explained_variance_,np.cumsum(fit_pca.explained_variance_ratio_*100).round(0)): plt.text(x,y,prop,fontsize=13)
少しみにくいですが、大体第4主成分ぐらいまでで、全データの7割ぐらいを説明することができるようです。もし、このデータで他の解析をする場合は、最小で4つの変数まで次元を圧縮できそうです。
その次の結果は、それぞれの栄養素が主成分からどれほどの影響を受けているのか、つまり因子負荷量です。この因子負荷量と第一、第二主成分との関連をグラフにしてみます。
plt.figure(figsize=(8, 6)) for x, y, name in zip(fit_pca.components_[0], fit_pca.components_[1], food_nut.columns): plt.text(x, y, name) plt.scatter(fit_pca.components_[0], fit_pca.components_[1], alpha=0.8) plt.grid() plt.xlabel("PC1",fontsize=19) plt.ylabel("PC2",fontsize=19) plt.tick_params(labelsize=15) plt.show()
多くの栄養素が第一主成分のマイナス方向に偏っています。このように直感的な理解と因子負荷量の方向が異なることがありますが、データのバラツキを説明するという観点から言えば方向の解釈はプラスマイナスを逆転させても問題ありません。 ただし、三大栄養素の一つである炭水化物や砂糖、食物繊維は少し方向性が違うみたいです（ちなみに食物繊維と砂糖を合わせたものが炭水化物です）。代わりに第二主成分では炭水化物系は強く関連していそうなので、完全食を探すには第二主成分も考慮した方が良いかもしれません。
完全食の王者を決める
さて、これで準備が整いました。完全食の王者を決めましょう。ここからはコードが少ないのでRとPythonを同時に書いていきます。とりあえずは第一主成分でデータをソートし、上位top10の食品を抜き出してみます。(Pythonの方は主成分得点の符号が逆転しています。)
Rの方
food %>% mutate(PC1 = pc$x[,1], PC2 = pc$x[,2]) %>% arrange(desc(PC1),desc(PC2)) %>% select(2,3) %>% head(10)
Pythonの方
food["PC1"] = pca[:,0] food["PC2"] = pca[:,1] food[["ID","FoodGroup","ShortDescrip","PC1","PC2"]].sort_values(by="PC1").head(10)
| ID | FoodGroup | ShortDescrip | PC1 | PC2 |
|---|---|---|---|---|
| 14654 | Beverages | BEVERAGES,NUTRITIONAL SHAKE MIX, HI PROT,PDR |
4.0 | 0.9 |
| 14055 | Beverages | BEVERAGES,UNILEVER,SLIMFAST SHAKE MIX, HI PROT,PDR,3-2-1 PLAN |
3.8 | 2.2 |
| 08077 | Breakfast Cereals | CEREALS RTE,GENERAL MILLS, WHL GRAIN TOTAL |
3.8 | 2.2 |
| 08028 | Breakfast Cereals | CEREALS RTE,KELLOGG, KELLOGG'S ALL-BRAN COMPLETE WHEAT FLAKES |
3.8 | 2.1 |
| 08001 | Breakfast Cereals | CEREALS RTE,KELLOGG, KELLOGG'S ALL-BRAN ORIGINAL |
3.6 | 2.1 |
| 25016 | Snacks | FORMUL BAR,MA SNACKFO US, SNICKE MARATH ENER BAR,ALL FLAVO |
3.6 | 2.1 |
| 25021 | Snacks | FORMULATED BAR,LUNA BAR, NUTZ OVER CHOC |
3.5 | 1.9 |
| 08058 | Breakfast Cereals | CEREALS RTE,KELLOGG, KELLOGG'S PRODUCT 19 |
3.4 | 2.1 |
| 25015 | Snacks | FORMBAR,MARSSNACK,SNICKERS MARATHONPROTPERFBAR,CARMELNUTRUSH |
3.4 | 2.0 |
| 08247 | Breakfast Cereals | CEREALS RTE,GENERAL MILLS, TOTAL RAISIN BRAN |
3.3 | 2.0 |
第一主成分トップ10には、飲み物グループ (Beverages)や朝食シリアルグループ(Breakfast Cereals)、スナックグループ(Snacks)がランクインしました。ShortDescriptionをみてみると、飲み物グループは野菜ジュース系、朝食シリアルグループはブランシリアル、スナックグループはエナジーバーなどみたいです。これは確かに完全食に近いですね！
続いて、同じように第二主成分でソートしてみます。
| ID | FoodGroup | ShortDescrip | PC1 | PC2 |
|---|---|---|---|---|
| 11634 | Vegetables and Vegetable Products |
PEPPERS,SWT,GRN, FREEZE-DRIED |
2.6 | 3.0 |
| 11931 | Vegetables and Vegetable Products |
PEPPERS,SWT,RED, FREEZE-DRIED |
2.7 | 3.0 |
| 02012 | Spices and Herbs | CORIANDER LEAF,DRIED | 3.0 | 2.9 |
| 11548 | Vegetables and Vegetable Products |
TOMATO POWDER | 2.4 | 2.9 |
| 14409 | Beverages | BEVERAGES,ORANGE-FLAVOR DRK, BRKFST TYPE,LO CAL,PDR |
-2.0 | 2.9 |
| 14426 | Beverages | ORANGE DRK,BRKFST TYPE, W/ JUC & PULP,FRZ CONC |
-0.3 | 2.9 |
| 14436 | Beverages | ORANGE BRKFST DRK,RTD ,W/ ADDED NUTR |
-1.9 | 2.9 |
| 02023 | Spices and Herbs | MARJORAM,DRIED | 2.7 | 2.8 |
| 02029 | Spices and Herbs | PARSLEY,DRIED | 2.9 | 2.8 |
なんだかスパイス系の食品が多く入ってきています。これは趣旨と離れている気がします...ということで、第一主成分トップ10の中から第二主成分も考慮して完全食を決めたいと思います！
もう一度第一主成分の結果に戻ってみると、第一主成分ランキング2位、3位が、第二主成分も考慮すると僅差です。この二つの最も大きな違いはズバリ液体か固体かです！笑甲乙付け難いので、朝とっても忙しい人は2位の食品を、液体じゃ物足りんという人は3位の食品を食べると良いでしょう！
ということで、この二つを"完全食"と任命します。（ちなみに会社名もデータに含まれているので実際に購入できるかもしれません。）
補足：因子分析と主成分分析
主成分分析と因子分析はよく混同されて説明されている記事をたまに見かけますが、全く違います。主成分分析とは、今回みてきたように、データを主成分に圧縮する手法であるのに対して、因子分析は、現在得られているデータが何らかの因子から発生していると仮定して、そのデータの因子を分析する手法です。
例えば、食品の栄養素というのは何かの因子から発生しているわけではありません。ここに因子分析を用いるということは、食品に含まれる栄養素に何かの力が働いて、その含有量が決めらることを仮定することになります。
式で見てみると主成分分析が
とすることで、二つの変数を一つに圧縮するのに対して、因子解析は
となる因子を求めます。とは言うもののこの二つの違いは実際に分かりにくく、同時にやると混乱するので、まずは理解しやすい主成分分析からマスターすることをオススメします。
まとめ
今回は多変量解析をする上でも重宝する主成分分析を紹介しました。今回のように主成分分析の結果自体からデータをみることもできますし、圧縮した後の主成分を使って、新たに他の解析をすることもできます。
色々と応用範囲が広いので、ざっくりとした概要と、具体的な解析手法は知っておいても損はありませんね。
※本記事は筆者が個人的に学んだことをまとめた記事なります。所属する組織の意見・見解とは無関係です。また、数学の記法や詳細な理論、用語等で誤りがあった際はご指摘頂けると幸いです。
参考
主成分分析をわかりやすく解説したサイト等はりますが、因子分析の解説はなかなか見当たりません。”が”、こちらの書籍は、回帰分析から因子分析などの高度な解析手法まで、とてもわかりやすくまとまっています。多変量解析を最初に学ぶにはとっておきの一冊だと思います。
英語のサイトですが、主成分分析と因子分析の違いについはおそらく一番わかりやすいです。
Principal Components (PCA) and Exploratory Factor Analysis (EFA) with SPSS
Pythonはこちらを参考にさせて頂きました。
数学的な導出に関しては、以下二つがわかりやすいです。
