Part 3 主成分分析
Chap.0 全体の流れ
Part3では、主成分分析を行う。主成分分析とは、多くの説明変数1を、より数少ない指標や合成変数2で表すことである。つまり、統計的にデータを要約することで、直感的に理解しづらい多次元のデータ全体の傾向や雰囲気を掴むということである。例えば、体型を表す指標として身長と体重の2変数があるが、これらを1変数のBMIにまとめるようなものである。
ChemTHEATREのデータの場合、一つのサンプルに対して、いくつかの化学物質の濃度を計測している事が多い。つまり、各サンプルについて説明変数が複数あるということである。今回は、各説明変数の持つ情報を可能な限り損なわないよう主成分分析し、データの概要を掴むことを目指す。
Chap.1 ライブラリの読み込み
%matplotlib inline
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
import matplotlib.ticker as ticker
import sklearn
from sklearn.decomposition import PCA
まず、必要なライブラリを読み込む。1行目は例によって、Jupyter Notebook内にmatplotlibで描画したグラフを表示するマジックコマンドである。
続いて、2行目以降が必要なライブラリである。なお、これらのライブラリはAnacondaにインストールされているものである。
今回主成分分析に利用するscikit-learnは、オープンソース3の機械学習ライブラリであり、他に分類4や回帰4、クラスタリング4などの様々なアルゴリズムも実装されている。
| ライブラリ | 概要 | 今回の使用目的 | 公式URL |
|---|---|---|---|
| NumPy | 数値計算ライブラリ | 統計処理上の数値計算に利用 | https://www.numpy.org |
| pandas | データ分析ライブラリ | データ読み込みや整形に利用 | https://pandas.pydata.org |
| Matplotlib | グラフ描画ライブラリ | データの可視化に利用 | https://matplotlib.org |
| scikit-learn | 機械学習ライブラリ | 主成分分析（PCA）の実装に利用 | https://scikit-learn.org |
Chap.2 データの読み込み
Chap.1で必要なライブラリが読み込めたので、次に今回取り扱うデータを準備する。今回取り扱うデータは、スナメリ（Neophocaena phocaenoides）、スジイルカ（Stenella coeruleoalba）、カズハゴンドウ（Peponocephala electra）の三種である。いずれもChemTHEATREのSample Searchからmeasureddataとsamplesをダウンロードして、Part1やPart2と同様に読み込む。
- ChemTHEATREのメニューバーから「Sample Search」を選択する。
- 「Sample Type」で「Biotic - Mammals - Marine mammals」」を選択する。
- 「Scientific Name」から「Neophocaena phocaenoides」を選択する。
- 「Search」ボタンをクリックすると，条件に合致する試料の一覧が出力される。
- 「Export samples」で試料の情報が，「Export measured data」で目的の化学物質の測定値が，タブ区切りのテキストファイルとして出力する。
- 他の2種も同様に繰り返す。
各3つずつのファイルの読み込みが成功したら、measureddata、samplesのそれぞれ3つのファイルを後の処理で使いやすいよう、結合してひとまとめにする。この際、Pandasのconcat関数を使用すると便利である。パラメーターのignore_indexをTrueにすると、結合するPandasのDataFrame5の列番号が振り直される。
data_file1 = "measureddata_20191004033314.tsv" #変数に入力する文字列を、各自のmeasureddataのtsvファイル名に変更する
data_file2 = "measureddata_20191004033325.tsv" #変数に入力する文字列を、各自のmeasureddataのtsvファイル名に変更する
data_file3 = "measureddata_20191004033338.tsv" #変数に入力する文字列を、各自のmeasureddataのtsvファイル名に変更する
data1 = pd.read_csv(data_file1, delimiter="\t")
data2 = pd.read_csv(data_file2, delimiter="\t")
data3 = pd.read_csv(data_file3, delimiter="\t")
data = pd.concat([data1, data2, data3], ignore_index=True) #Pandasのconcat関数で結合する
data
| MeasuredID | ProjectID | SampleID | ScientificName | ChemicalID | ChemicalName | ExperimentID | MeasuredValue | AlternativeData | ... | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 80 | PRA000002 | SAA000087 | Neophocaena phocaenoides | CH0000154 | TBT | EXA000001 | 170.0000 | NaN | ... |
| 1 | 81 | PRA000002 | SAA000087 | Neophocaena phocaenoides | CH0000155 | DBT | EXA000001 | 220.3591 | NaN | ... |
| 2 | 82 | PRA000002 | SAA000087 | Neophocaena phocaenoides | CH0000156 | MBT | EXA000001 | 44.5445 | NaN | ... |
| 3 | 83 | PRA000002 | SAA000087 | Neophocaena phocaenoides | CH0000157 | ΣBTs | EXA000001 | 434.9036 | NaN | ... |
| 4 | 84 | PRA000002 | SAA000087 | Neophocaena phocaenoides | CH0000158 | TPT | EXA000001 | 12.9220 | NaN | ... |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 9774 | 24272 | PRA000036 | SAA002163 | Peponocephala electra | CH0000151 | cis-nonachlor | EXA000001 | 170.0000 | NaN | ... |
| 9775 | 24056 | PRA000036 | SAA002163 | Peponocephala electra | CH0000152 | ΣCHLs | EXA000001 | 1100.0000 | NaN | ... |
| 9776 | 24376 | PRA000036 | SAA002163 | Peponocephala electra | CH0000153 | HCB | EXA000001 | 120.0000 | NaN | ... |
| 9777 | 24646 | PRA000036 | SAA002163 | Peponocephala electra | CH0000533 | TCPMe | EXA000001 | 6.1000 | NaN | ... |
| 9778 | 24700 | PRA000036 | SAA002163 | Peponocephala electra | CH0000534 | TCPMOH | EXA000001 | 16.0000 | NaN | ... |
9779 rows × 13 columns
sample_file1 = "samples_20191004033311.tsv" #変数に入力する文字列を、各自のsamplesのtsvファイル名に変更する
sample_file2 = "samples_20191004033323.tsv" #変数に入力する文字列を、各自のsamplesのtsvファイル名に変更する
sample_file3 = "samples_20191004033334.tsv" #変数に入力する文字列を、各自のsamplesのtsvファイル名に変更する
sample1 = pd.read_csv(sample_file1, delimiter="\t")
sample2 = pd.read_csv(sample_file2, delimiter="\t")
sample3 = pd.read_csv(sample_file3, delimiter="\t")
sample = pd.concat([sample1, sample2, sample3], ignore_index=True) #Pandasのconcat関数で結合する
sample
| ProjectID | SampleID | SampleType | TaxonomyID | UniqCodeType | UniqCode | SampleName | ScientificName | CommonName | CollectionYear | ... | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | PRA000002 | SAA000087 | ST004 | 34892 | es-BANK | EW00884 | NaN | Neophocaena phocaenoides | Finless porpoises | 1996 | ... |
| 1 | PRA000002 | SAA000088 | ST004 | 34892 | es-BANK | EW00812 | NaN | Neophocaena phocaenoides | Finless porpoises | 1999 | ... |
| 2 | PRA000002 | SAA000089 | ST004 | 34892 | es-BANK | EW00873 | NaN | Neophocaena phocaenoides | Finless porpoises | 1995 | ... |
| 3 | PRA000002 | SAA000090 | ST004 | 34892 | es-BANK | EW04787 | NaN | Neophocaena phocaenoides | Finless porpoises | 2000 | ... |
| 4 | PRA000002 | SAA000091 | ST004 | 34892 | es-BANK | EW00867 | NaN | Neophocaena phocaenoides | Finless porpoises | 1998 | ... |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 312 | PRA000036 | SAA002159 | ST004 | 103596 | es-BANK | EW04779 | 060301-1 | Peponocephala electra | Melon-headed whale | 2006 | ... |
| 313 | PRA000036 | SAA002160 | ST004 | 103596 | es-BANK | EW00115 | M32625 | Peponocephala electra | Melon-headed whale | 2001 | ... |
| 314 | PRA000036 | SAA002161 | ST004 | 103596 | es-BANK | EW00122 | M32633 | Peponocephala electra | Melon-headed whale | 2001 | ... |
| 315 | PRA000036 | SAA002162 | ST004 | 103596 | es-BANK | EW00116 | M32626 | Peponocephala electra | Melon-headed whale | 2001 | ... |
| 316 | PRA000036 | SAA002163 | ST004 | 103596 | es-BANK | EW00117 | M32627 | Peponocephala electra | Melon-headed whale | 2001 | ... |
317 rows × 66 columns
Chap.3 データの下処理
データの読み込みが終了したら、次はデータの分析ができる形に下処理をする。まずmeasureddataから今回利用するデータを取り出すところから始める。
今回利用するのは、スナメリ、スジイルカ、カズハゴンドウの3種についての、5種類の化学物質（ΣPCBs, ΣDDTs, ΣPBDEs, ΣCHLs, ΣHCHs）である。
最初に、単位が[ng/g lipid]のデータを抽出する。これは、Part 1でも述べたように単位が異なるデータ同士は、単純に比較できないからである。
その次に、化学物質名が先述の5種のいずれかであるデータを抽出して、reset_indexメソッドを利用して行番号を振り直す。
data_lipid = data[(data["Unit"] == "ng/g lipid")]
data_lipid = data_lipid[(data_lipid["ChemicalName"] == "ΣPCBs") | (data_lipid["ChemicalName"] == "ΣDDTs") |
(data_lipid["ChemicalName"] == "ΣPBDEs") | (data_lipid["ChemicalName"] == "ΣCHLs") |
(data_lipid["ChemicalName"] == "ΣHCHs")]
data_lipid = data_lipid.reset_index(drop=True)
data_lipid[data_lipid["SampleID"] == "SAA001941"] #SampleIDがSAA001941のデータ一覧
| MeasuredID | ProjectID | SampleID | ScientificName | ChemicalID | ChemicalName | ExperimentID | MeasuredValue | AlternativeData | ... | |
|---|---|---|---|---|---|---|---|---|---|---|
| 115 | 25112 | PRA000030 | SAA001941 | Neophocaena phocaenoides | CH0000033 | ΣDDTs | EXA000001 | 9400.0 | NaN | ... |
| 116 | 25098 | PRA000030 | SAA001941 | Neophocaena phocaenoides | CH0000096 | ΣPCBs | EXA000001 | 1100.0 | NaN | ... |
| 117 | 25140 | PRA000030 | SAA001941 | Neophocaena phocaenoides | CH0000146 | ΣHCHs | EXA000001 | 41.0 | NaN | ... |
| 118 | 25126 | PRA000030 | SAA001941 | Neophocaena phocaenoides | CH0000152 | ΣCHLs | EXA000001 | 290.0 | NaN | ... |
measureddataのデータ抽出が終了したら、次は、samplesとmeasureddataのDataFrameの統合を行う。
Part1、Part2では各化学物質のmeasureddataにsamplesの情報（学名や採集地など）を付け加えたが、今回は、各sampleに化学物質の計測値を付与するので処理が今までと異なることに注意する。
# SampleIDと学名だけのDataFrameを作成
df = sample[["SampleID", "ScientificName"]]
for chem in ["ΣCHLs", "ΣDDTs", "ΣHCHs", "ΣPBDEs", "ΣPCBs"]: #化学物質名でループ処理
# リストchem番目の化学物質のみを対象としたDataFrameを作成
data_lipid_chem = data_lipid[data_lipid["ChemicalName"] == chem]
# SampleIDと測定値だけを残して、列名"MeasuredValue"を化学物質名に変更
data_lipid_chem = data_lipid_chem[["SampleID", "MeasuredValue"]].rename(columns={"MeasuredValue": chem})
data_lipid_chem = data_lipid_chem.drop_duplicates(subset='SampleID')
# SampleIDでdf(ID+学名)とdata_lipid_chem(ID+測定値)をマージ
df = pd.merge(df, data_lipid_chem, on="SampleID")
df = df.dropna(how="any") #NaNが含まれる行を除去
df = df.drop("SampleID", axis=1)
df
| ScientificName | ΣCHLs | ΣDDTs | ΣHCHs | ΣPBDEs | ΣPCBs | |
|---|---|---|---|---|---|---|
| 0 | Neophocaena phocaenoides | 770.0 | 68000.0 | 1100.0 | 170.0 | 5700.0 |
| 1 | Neophocaena phocaenoides | 1200.0 | 140000.0 | 3500.0 | 120.0 | 6500.0 |
| 2 | Neophocaena phocaenoides | 1000.0 | 140000.0 | 5500.0 | 86.0 | 5500.0 |
| 3 | Neophocaena phocaenoides | 950.0 | 130000.0 | 6600.0 | 100.0 | 5800.0 |
| 4 | Neophocaena phocaenoides | 1400.0 | 280000.0 | 6100.0 | 140.0 | 11000.0 |
| ... | ... | ... | ... | ... | ... | ... |
| 89 | Peponocephala electra | 3000.0 | 15000.0 | 170.0 | 260.0 | 12000.0 |
| 90 | Peponocephala electra | 5100.0 | 23000.0 | 380.0 | 490.0 | 19000.0 |
| 91 | Peponocephala electra | 5700.0 | 33000.0 | 240.0 | 300.0 | 25000.0 |
| 92 | Peponocephala electra | 2800.0 | 12000.0 | 220.0 | 230.0 | 9300.0 |
| 93 | Peponocephala electra | 5700.0 | 27000.0 | 240.0 | 180.0 | 21000.0 |
94 rows × 6 columns
上で表示されている通り、今の状態では学名がそのままDataFrameに格納されている。種名は今後可視化する際に利用するので、学名を表す変数を便宜的に用意する。最初に、各学名に値（スナメリ：0、スジイルカ：1、カズハゴンドウ：2）を対応させた辞書を作る。
species = df["ScientificName"].unique().tolist() #dfにある学名の一覧をリストで出力する
class_dic = {}
for i in range(len(species)):
class_dic[species[i]] = i
class_dic
{'Neophocaena phocaenoides': 0,
'Stenella coeruleoalba': 1,
'Peponocephala electra': 2}
次に、dfにClass列を追加し、各行の学名に対応した値をClass列に入力する。
```python df["Class"] = 0 for irow in range(len(df)): df.at[irow, "Class"] = class_dic[df.at[irow, "ScientificName"]] df = df.loc[:, ["Class", "ScientificName", "ΣPCBs", "ΣDDTs", "ΣPBDEs", "ΣCHLs", "ΣHCHs"]] df ```| Class | ScientificName | ΣPCBs | ΣDDTs | ΣPBDEs | ΣCHLs | ΣHCHs | |
|---|---|---|---|---|---|---|---|
| 0 | 0 | Neophocaena phocaenoides | 5700.0 | 68000.0 | 170.0 | 770.0 | 1100.0 |
| 1 | 0 | Neophocaena phocaenoides | 6500.0 | 140000.0 | 120.0 | 1200.0 | 3500.0 |
| 2 | 0 | Neophocaena phocaenoides | 5500.0 | 140000.0 | 86.0 | 1000.0 | 5500.0 |
| 3 | 0 | Neophocaena phocaenoides | 5800.0 | 130000.0 | 100.0 | 950.0 | 6600.0 |
| 4 | 0 | Neophocaena phocaenoides | 11000.0 | 280000.0 | 140.0 | 1400.0 | 6100.0 |
| ... | ... | ... | ... | ... | ... | ... | ... |
| 89 | 2 | Peponocephala electra | 12000.0 | 15000.0 | 260.0 | 3000.0 | 170.0 |
| 90 | 2 | Peponocephala electra | 19000.0 | 23000.0 | 490.0 | 5100.0 | 380.0 |
| 91 | 2 | Peponocephala electra | 25000.0 | 33000.0 | 300.0 | 5700.0 | 240.0 |
| 92 | 2 | Peponocephala electra | 9300.0 | 12000.0 | 230.0 | 2800.0 | 220.0 |
| 93 | 2 | Peponocephala electra | 21000.0 | 27000.0 | 180.0 | 5700.0 | 240.0 |
94 rows × 7 columns
ピボットテーブルを使った方が簡単だったので，data_lipid以下を下記のように修正。
df = pd.pivot_table(data_lipid, index=['ScientificName', 'SampleID'], columns='ChemicalName', values=['MeasuredValue'])
df = df.dropna(how="any") #NaNが含まれる行を除去
df
| MeasuredValue | ||||||
|---|---|---|---|---|---|---|
| ChemicalName | ΣCHLs | ΣDDTs | ΣHCHs | ΣPBDEs | ΣPCBs | |
| ScientificName | SampleID | |||||
| Neophocaena phocaenoides | SAA001903 | 770.0 | 68000.0 | 1100.0 | 170.0 | 5700.0 |
| SAA001904 | 1200.0 | 140000.0 | 3500.0 | 120.0 | 6500.0 | |
| SAA001905 | 1000.0 | 140000.0 | 5500.0 | 86.0 | 5500.0 | |
| SAA001906 | 950.0 | 130000.0 | 6600.0 | 100.0 | 5800.0 | |
| SAA001907 | 1400.0 | 280000.0 | 6100.0 | 140.0 | 11000.0 | |
| ... | ... | ... | ... | ... | ... | ... |
| Stenella coeruleoalba | SAA001984 | 730.0 | 3500.0 | 40.0 | 84.0 | 3200.0 |
| SAA001985 | 7300.0 | 44000.0 | 300.0 | 530.0 | 26000.0 | |
| SAA001986 | 12000.0 | 88000.0 | 490.0 | 850.0 | 42000.0 | |
| SAA001987 | 9100.0 | 56000.0 | 380.0 | 650.0 | 32000.0 | |
| SAA001988 | 11000.0 | 67000.0 | 520.0 | 610.0 | 36000.0 | |
94 rows × 5 columns
これ以降，ピボットテーブルのデータフレームを使って解析。
Chap.4 データ概要の可視化
ここまででデータの下準備は完了したので、ひとまずデータ全体を可視化して確認してみる。今回取り扱っているデータは、5次元（説明変数が5つある）である。人間は5次元空間を把握できないので、この5次元データすべての分布をそのまま可視化するのは不可能である。ここでは、${}_5 C_2 = 10$通りの分布図が生成する。
複数変数の散布図の可視化にはpandasのpandas.plotting.scatter_matrix()関数を利用するのが便利である。scatter_matrix()関数では、可視化する範囲をframeパラメーターで指定する。今回は、dfの全行・2列目6以降すべてである。
pd.plotting.scatter_matrix(
frame=df["MeasuredValue"],
c=list(df["MeasuredValue"].iloc[:, 0]),
alpha=0.5)
plt.show()
Chap.5 主成分分析
Chap.4では、データ全体の概要として、${}_5 C_2 = 10$通りの分布図が生成された。しかしこれでは、2種類の化学物質の組み合わせごとのデータの分布しか確認できない。つまり、5次元のデータ全体の概形は掴めないのである。そこで主成分分析で5次元データを要約して可視化する。
まず、このままでは各軸の長さが異なるので、各変数について標準化7し、5次元すべての軸の長さを統一する。
dfs = df.apply(lambda x: (x-x.mean())/x.std(), axis=0)
dfs
| MeasuredValue | ||||||
|---|---|---|---|---|---|---|
| ChemicalName | ΣCHLs | ΣDDTs | ΣHCHs | ΣPBDEs | ΣPCBs | |
| ScientificName | SampleID | |||||
| Neophocaena phocaenoides | SAA001903 | -0.743053 | 0.320011 | 0.179064 | -0.308268 | -0.768856 |
| SAA001904 | -0.588580 | 1.574027 | 1.297534 | -0.528241 | -0.693375 | |
| SAA001905 | -0.660428 | 1.574027 | 2.229592 | -0.677822 | -0.787727 | |
| SAA001906 | -0.678390 | 1.399859 | 2.742224 | -0.616230 | -0.759421 | |
| SAA001907 | -0.516732 | 4.012392 | 2.509209 | -0.440252 | -0.268792 | |
| ... | ... | ... | ... | ... | ... | ... |
| Stenella coeruleoalba | SAA001984 | -0.757422 | -0.803378 | -0.314927 | -0.686621 | -1.004736 |
| SAA001985 | 1.602782 | -0.097994 | -0.193759 | 1.275535 | 1.146484 | |
| SAA001986 | 3.291208 | 0.668349 | -0.105213 | 2.683360 | 2.656112 | |
| SAA001987 | 2.249413 | 0.111009 | -0.156477 | 1.803469 | 1.712595 | |
| SAA001988 | 2.931969 | 0.302594 | -0.091233 | 1.627491 | 2.090001 | |
94 rows × 5 columns
あるいは標準化にStandardScalerを使う方法もある。
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
scaler.fit(df)
scaler.transform(df)
dfs = pd.DataFrame(scaler.transform(df), columns=df.columns)
dfs
| MeasuredValue | |||||
|---|---|---|---|---|---|
| ChemicalName | ΣCHLs | ΣDDTs | ΣHCHs | ΣPBDEs | ΣPCBs |
| 0 | -0.747037 | 0.321727 | 0.180024 | -0.309921 | -0.772979 |
| 1 | -0.591736 | 1.582467 | 1.304491 | -0.531073 | -0.697093 |
| 2 | -0.663969 | 1.582467 | 2.241547 | -0.681457 | -0.791950 |
| 3 | -0.682027 | 1.407365 | 2.756927 | -0.619534 | -0.763493 |
| 4 | -0.519503 | 4.033907 | 2.522663 | -0.442612 | -0.270233 |
| ... | ... | ... | ... | ... | ... |
| 89 | -0.761484 | -0.807686 | -0.316615 | -0.690303 | -1.010123 |
| 90 | 1.611376 | -0.098520 | -0.194798 | 1.282374 | 1.152631 |
| 91 | 3.308856 | 0.671933 | -0.105778 | 2.697748 | 2.670354 |
| 92 | 2.261475 | 0.111604 | -0.157316 | 1.813139 | 1.721777 |
| 93 | 2.947690 | 0.304217 | -0.091722 | 1.636218 | 2.101208 |
94 rows × 5 columns
標準化したら、主成分分析を行う。ここではStandardScalerの結果を使う。
主成分分析には、scikit-learnのsklearn.decomposition.PCA()を利用する。
主成分分析では、まずデータの分布のばらつきが最大となる向きに軸を引くように合成変数（第1主成分）を作り、その合成変数で失われるデータを補完するように第2主成分、第3主成分と合成変数を作っていく。その後、各主成分を軸とした空間にデータを転写すれば、後は解釈するのみである。
pca = PCA()
feature = pca.fit(dfs) #主成分分析
feature = pca.transform(dfs) #主成分空間にデータを転写
pd.DataFrame(feature, columns=["PC{}".format(x + 1) for x in range(len(dfs.columns))])
| PC1 | PC2 | PC3 | PC4 | PC5 | |
|---|---|---|---|---|---|
| 0 | -0.931614 | 0.574985 | 0.348535 | -0.188477 | -0.178934 |
| 1 | -0.525203 | 2.177125 | 0.038871 | -0.242288 | -0.498098 |
| 2 | -0.700197 | 2.892482 | -0.248495 | 0.251040 | -0.282242 |
| 3 | -0.720222 | 3.155379 | -0.366880 | 0.634978 | -0.017083 |
| 4 | 0.554487 | 4.585518 | 0.039662 | -0.952484 | -0.965309 |
| ... | ... | ... | ... | ... | ... |
| 89 | -1.607963 | -0.471519 | 0.077579 | 0.103413 | -0.068820 |
| 90 | 2.190969 | -0.704970 | 0.009171 | 0.563020 | 0.053168 |
| 91 | 4.974287 | -0.699630 | 0.023300 | 0.765578 | 0.013150 |
| 92 | 3.218248 | -0.753396 | -0.005256 | 0.688830 | 0.065009 |
| 93 | 3.804695 | -0.750210 | -0.513720 | 0.727646 | -0.224109 |
94 rows × 5 columns
固有ベクトルについても求めてみる。固有ベクトルとは、元データを主成分軸空間に転写する際に加える「軸を傾ける力」のようなものである。
Pythonでの固有ベクトルの計算は、scikit-learnのPCAに含まれる、components_で実装されている。
pd.DataFrame(pca.components_, columns=df.columns[2:], index=["PC{}".format(x + 1) for x in range(len(dfs.columns))])
| MeasuredValue | |||||
|---|---|---|---|---|---|
| ChemicalName | ΣCHLs | ΣDDTs | ΣHCHs | ΣPBDEs | ΣPCBs |
| PC1 | 0.572482 | 0.307088 | -0.005288 | 0.488078 | 0.582848 |
| PC2 | -0.212874 | 0.641425 | 0.730929 | -0.048824 | -0.081347 |
| PC3 | -0.328297 | 0.180897 | -0.249982 | 0.782267 | -0.430192 |
| PC4 | 0.427788 | -0.529630 | 0.566847 | 0.279604 | -0.370130 |
| PC5 | -0.579799 | -0.425488 | 0.286195 | 0.263203 | 0.575856 |
固有ベクトルがわかったので、ここで各主成分の固有値について求める。固有値とは、固有ベクトルに付随する係数で、各主成分軸の大きさを示す指標である。
Pythonでは、scikit-learnのPCAに含まれる、explained_variance_を利用すれば、寄与率を求めることができる。
pd.DataFrame(pca.explained_variance_, index=["PC{}".format(x + 1) for x in range(len(dfs.columns))])
| 0 | |
|---|---|
| PC1 | 2.301828 |
| PC2 | 1.526633 |
| PC3 | 0.681585 |
| PC4 | 0.450884 |
| PC5 | 0.092834 |
次に、寄与率を求める。寄与率とは、各主成分軸がデータのどれだけの情報を説明しているかを示す指標である。つまり寄与率60%の主成分はデータ全体の60%を説明できており、40%を損失しているということである。
Pythonでは、scikit-learnのPCAに含まれる、explained_variance_ratio_を利用すれば、寄与率を求めることができる。
pd.DataFrame(pca.explained_variance_ratio_, index=["PC{}".format(x + 1) for x in range(len(dfs.columns))])
| 0 | |
|---|---|
| PC1 | 0.455468 |
| PC2 | 0.302078 |
| PC3 | 0.134867 |
| PC4 | 0.089218 |
| PC5 | 0.018369 |
ここで累積寄与率も求める。なぜ累積寄与率が必要かというと、「あと何割のデータが失われているか」を見るためである。つまり、ある一定以上の累積寄与率となる主成分まで見れば、次元数を削減してデータ全体の概形をつかめるというわけである。
ここでは、numpyのcumsum関数を利用して累積した合計値を求めたものを利用している。
plt.plot([0] + list( np.cumsum(pca.explained_variance_ratio_)), "-o")
plt.xlabel("Number of principal components")
plt.ylabel("Cumulative contribution rate")
plt.grid()
plt.show()
生成した各主成分軸についてだが、このままではどのような意味を持つかよくわからない。しかし、主成分における各変数の寄与度を可視化することで各主成分のもつ意味が見えてくる。ここでは、第1主成分と第2主成分における各変数の寄与度を散布図にプロットしてみる。
```python plt.figure() for x, y, name in zip(pca.components_[0], pca.components_[1], df["MeasuredValue"].columns[:]): plt.text(x, y, name) plt.scatter(pca.components_[0], pca.components_[1], alpha=0.8) plt.grid() plt.xlabel("PC1") plt.ylabel("PC2") plt.show() ```
最後に主成分分析した結果を可視化する。一般に固有値が1以上または、累積寄与率が80％以上までの主成分を採用し、採用した主成分の組み合わせで分布をみる。今回の場合は、第2主成分までを採用するとして、グラフに可視化する。一般の散布図同様、matplotlibのscatter関数を利用し、パラメーターのX,Yそれぞれに第1主成分（featureの0列目）、第2主成分（featureの1列目）を指定すれば良い。
plt.figure()
plt.scatter(x=feature[:, 0], y=feature[:, 1], alpha=0.8, c=list(df.iloc[:, 0])) #cにplotの色分けを記述
plt.grid()
plt.xlabel("PC1")
plt.ylabel("PC2")
plt.show()
脚注
1説明変数とは、求める結果である目的変数を説明する変数のことである。言い換えれば、結果を示す目的変数に対して、その原因となっている変数が説明変数である。独立変数とも。
2複数の変数をまとめて作る変数のことである。主成分分析では、各主成分が合成変数である。
3プログラムの設計図であるソースコードが、一般に広く開示されているソフトウェアやライブラリのこと。ライセンスで認められた範囲内であれば、万人が無償で使用できる。また、誰でも開発できるため、利用者目線のアップデートや少人数では気づかないバグの修正がなされやすい。
PythonやJupyter Notebook、Anacondaはすべてオープンソースのソフトウェアである。
4いずれも、機械学習で利用される手法である。分類は画像分類などで、回帰は株価などの将来予測に利用される。クラスタリングは、データの傾向から類似するデータのグループ分けをする際に用いられる。
5Pandasで表形式のデータを格納するデータ型のこと。ExcelでCSV、TSV形式のファイルを読み込んだのような形で、セルや行・列を指定してデータを操作する。
6pythonでは「0」から数字が始まる。
7計測値から平均値を引き、標準偏差で割る操作のこと。