Part1 基本統計量とその可視化
全体の流れ
Part1では，ChemTHEATRE（https://chem-theatre.com/）に収録されているスナメリのデータを利用して，データの基本的な特性を見る基本統計量を分析し，箱ひげ図を利用して可視化する。
Chap.1 ライブラリの読み込み
%matplotlib inline
import numpy as np # as np, as pdなどはプログラムで呼び出すときの省略名を設定している
import pandas as pd # いちいちnumpyなどと都度入力するのがめんどうなので慣習としてこの様に記載される
import matplotlib.pyplot as plt
import seaborn as sns
１行目は、可視化の際に出力されるグラフがJupyter Notebook内で表示する設定1である。２行目以降で、今回のコードで必要なライブラリを読み込む。いずれもAnacondaに標準でインストールされているライブラリである。
具体的なライブラリの内容は以下の通りである。
| ライブラリ | 概要 | 今回の使用目的 | 公式URL |
|---|---|---|---|
| NumPy | 数値計算ライブラリ | 統計処理上の数値計算に利用 | https://www.numpy.org |
| pandas | データ分析ライブラリ | データ読み込みや整形に利用 | https://pandas.pydata.org |
| Matplotlib | グラフ描画ライブラリ | データの可視化に利用 | https://matplotlib.org |
| seaborn | 統計データ可視化ライブラリ | データの可視化に利用 | https://seaborn.pydata.org |
Chap.2 データの読み込み
先述の通り、今回はスナメリ（Neophocaena phocaenoides）のデータを取り扱う。ChemTHEATREのSample Searchからスナメリに関するデータをTSV形式2でダウンロードする。ChemTHEATREのSample Searchのデータは、測定結果（measureddata）と採集標本（samples）の2種類あるが、今回はその両方を使用する。
ChemTHEATREの操作マニュアルはChemTHEATRE Wikiサイトから入手できる。
簡単な操作ガイドは以下の通り。
- ChemTEHATREのトップページからメニューバーの「Sample Search」をクリックする。
- 「Sample Type」のプルダウンから「Biotic - Mammals - Marine mammals」を選択する。
- 「Scientific Name」のプルダウンからスナメリの学名「Neophocaena phocaenoides」を選択する。
- 「Search」ボタンをクリックする。
- 「Export samples (TSV)」ボタンでサンプル一覧が，「Export measured data (TSV)」ボタンで測定値一覧が，タブ区切りのテキストファイルとしてそれぞれエクスポートされる。
- エクスポートされたファイルを任意のディレクトリに保存し，以降の解析に使用する。
data_file = "measureddata_20190826061813.tsv" #変数に入力する文字列を、各自のmeasureddataのtsvファイル名に変更する
chem = pd.read_csv(data_file, delimiter="\t")
chem = chem.drop(["ProjectID", "ScientificName", "RegisterDate", "UpdateDate"], axis=1) #後でsamplesと結合する際に重複する列の削除
sample_file = "samples_20190826061810.tsv" #変数に入力する文字列を、各自のsamplesのtsvファイル名に変更する
sample = pd.read_csv(sample_file, delimiter="\t")
ダウンロードしたTSVファイルは、pandasの関数（read_csv）で読み込む。この際、各ファイルのパス（上セルでの"data_file"や"sample_file"）は、それぞれの環境に合わせたファイル名に変更する必要がある。サンプルコードでは、ノートブックと同じディレクトリにデータファイルがある環境で記述している。
Chap.3 データの下処理
今回扱うデータは、様々な化学物質の測定結果を一緒くたにしたものである。更に、測定結果と標本のデータがバラバラで存在している。つまり、このままでは扱いづらいので、扱いやすい形に整形する必要がある。具体的には下記の操作が必要である。
- chemにsampleを結合させる。
- 今回扱う化学物質の測定データのみを抽出する。
- 可視化に不要な列を削除する。
Sec.3-1 必要データの抽出
まず、上で示した1.と2.について処理を行う。
結合については、各化学物質を測定した標本のデータフレーム（sample）を各化学物質のデータフレーム（chem）に付け加えるという作業である。ここでは、chemの各データが持つSampleIDに対応するsampleのデータをchemの右側にくっつけるという処理をする。
df = pd.merge(chem, sample, on="SampleID")
今回扱うデータは、各化学物質がかなり細かく分類されている。それぞれの化学物質について計測結果を可視化してもよいが、今回は計測結果の概要をひと目で分かるように可視化したいので、化学物質を大まかに分類したときの各グループの合計値のデータを抽出する。
また、今回扱うデータには、2種類の異なる単位が混在している。異なる単位同士では、比較しても意味がないので、今回は単位ごとにデータを分割する。以下のコードではそれぞれの1行目でUnitを指定してlipid, wet換算のものを抽出した上で、その中でΣから始まる変数だけを抽出している。2行目ではΣから始まるものの、実際は化学物質濃度の合計値ではない変数をデータから削除している。
data_lipid = df[(df["Unit"] == "ng/g lipid") & df["ChemicalName"].str.startswith("Σ")]
data_lipid = data_lipid[(data_lipid["ChemicalName"] != "ΣOH-penta-PCB") & (data_lipid["ChemicalName"] != "ΣOH-hexa-PCB")
& (data_lipid["ChemicalName"] != "ΣOH-hepta-PCB") & (data_lipid["ChemicalName"] != "ΣOH-octa-PCB")]
data_wet = df[(df["Unit"] == "ng/g wet") & df["ChemicalName"].str.startswith("Σ")]
data_wet = data_wet[(data_wet["ChemicalName"] != "ΣOH-penta-PCB") & (data_wet["ChemicalName"] != "ΣOH-hexa-PCB")
& (data_wet["ChemicalName"] != "ΣOH-hepta-PCB") & (data_wet["ChemicalName"] != "ΣOH-octa-PCB")]
データの概形
Sec.3-1で処理した結果（単位がng/g lipidの方のみ）である。不要なデータを削除したことで、6種の化学物質についての131のデータのみが抽出できた。
data_lipid
| MeasuredID | SampleID | ChemicalID | ChemicalName | ExperimentID | MeasuredValue | AlternativeData | Unit | Remarks_x | ProjectID | ... | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 4371 | 24782 | SAA001903 | CH0000033 | ΣDDTs | EXA000001 | 68000.0 | NaN | ng/g lipid | NaN | PRA000030 | ... |
| 4372 | 24768 | SAA001903 | CH0000096 | ΣPCBs | EXA000001 | 5700.0 | NaN | ng/g lipid | NaN | PRA000030 | ... |
| 4381 | 24754 | SAA001903 | CH0000138 | ΣPBDEs | EXA000001 | 170.0 | NaN | ng/g lipid | NaN | PRA000030 | ... |
| 4385 | 24902 | SAA001903 | CH0000142 | ΣHBCDs | EXA000001 | 5.6 | NaN | ng/g lipid | NaN | PRA000030 | ... |
| 4389 | 24810 | SAA001903 | CH0000146 | ΣHCHs | EXA000001 | 1100.0 | NaN | ng/g lipid | NaN | PRA000030 | ... |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 4812 | 25125 | SAA001940 | CH0000152 | ΣCHLs | EXA000001 | 160.0 | NaN | ng/g lipid | NaN | PRA000030 | ... |
| 4816 | 25112 | SAA001941 | CH0000033 | ΣDDTs | EXA000001 | 9400.0 | NaN | ng/g lipid | NaN | PRA000030 | ... |
| 4817 | 25098 | SAA001941 | CH0000096 | ΣPCBs | EXA000001 | 1100.0 | NaN | ng/g lipid | NaN | PRA000030 | ... |
| 4818 | 25140 | SAA001941 | CH0000146 | ΣHCHs | EXA000001 | 41.0 | NaN | ng/g lipid | NaN | PRA000030 | ... |
| 4819 | 25126 | SAA001941 | CH0000152 | ΣCHLs | EXA000001 | 290.0 | NaN | ng/g lipid | NaN | PRA000030 | ... |
131 rows × 74 columns
Sec.3-2 不要な列の削除
続いて、不要な列の削除を行う。一般に、情報量が同じなら、扱うデータは容量が小さい方が良いので、無駄なデータは削ったほうがよい。今回扱うデータフレームには、データのない空欄が多く存在する。これは、python上ではN/A（浮動小数点数）として扱われ、情報はないが容量はそれなりに消費する。したがって、データフレームからすべてN/Aの列（情報のない列）をまるごと削除する。
data_lipid = data_lipid.dropna(how='all', axis=1)
data_wet = data_wet.dropna(how='all', axis=1)
データの概形
Sec.3-2の処理をした結果である。Sec.3-1の結果と比較して列の数（表の下に出るcolumnsの値）がかなり小さくなったことがわかる。
data_lipid
| MeasuredID | SampleID | ChemicalID | ChemicalName | ExperimentID | MeasuredValue | Unit | ProjectID | SampleType | TaxonomyID | ... | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 4371 | 24782 | SAA001903 | CH0000033 | ΣDDTs | EXA000001 | 68000.0 | ng/g lipid | PRA000030 | ST004 | 34892 | ... |
| 4372 | 24768 | SAA001903 | CH0000096 | ΣPCBs | EXA000001 | 5700.0 | ng/g lipid | PRA000030 | ST004 | 34892 | ... |
| 4381 | 24754 | SAA001903 | CH0000138 | ΣPBDEs | EXA000001 | 170.0 | ng/g lipid | PRA000030 | ST004 | 34892 | ... |
| 4385 | 24902 | SAA001903 | CH0000142 | ΣHBCDs | EXA000001 | 5.6 | ng/g lipid | PRA000030 | ST004 | 34892 | ... |
| 4389 | 24810 | SAA001903 | CH0000146 | ΣHCHs | EXA000001 | 1100.0 | ng/g lipid | PRA000030 | ST004 | 34892 | ... |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 4812 | 25125 | SAA001940 | CH0000152 | ΣCHLs | EXA000001 | 160.0 | ng/g lipid | PRA000030 | ST004 | 34892 | ... |
| 4816 | 25112 | SAA001941 | CH0000033 | ΣDDTs | EXA000001 | 9400.0 | ng/g lipid | PRA000030 | ST004 | 34892 | ... |
| 4817 | 25098 | SAA001941 | CH0000096 | ΣPCBs | EXA000001 | 1100.0 | ng/g lipid | PRA000030 | ST004 | 34892 | ... |
| 4818 | 25140 | SAA001941 | CH0000146 | ΣHCHs | EXA000001 | 41.0 | ng/g lipid | PRA000030 | ST004 | 34892 | ... |
| 4819 | 25126 | SAA001941 | CH0000152 | ΣCHLs | EXA000001 | 290.0 | ng/g lipid | PRA000030 | ST004 | 34892 | ... |
131 rows × 37 columns
Chap.4 ΣDDTsの可視化
Chap.3で処理したことで、6種類の化学物質について、2種の単位ごとのデータが準備できた。ここでは、まずΣDDTsのng/g lipidについて、データの散らばり具合を可視化する。
Sec.4-1 データの準備
まず可視化するにあたって、data_lipidからΣDDTsのデータを抽出する。
data_ddt = data_lipid[data_lipid["ChemicalName"] == "ΣDDTs"]
ddt_vals = data_ddt.loc[:, "MeasuredValue"]
data_ddt
| MeasuredID | SampleID | ChemicalID | ChemicalName | ExperimentID | MeasuredValue | Unit | ProjectID | SampleType | TaxonomyID | ... | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 4371 | 24782 | SAA001903 | CH0000033 | ΣDDTs | EXA000001 | 68000.0 | ng/g lipid | PRA000030 | ST004 | 34892 | ... |
| 4401 | 24783 | SAA001904 | CH0000033 | ΣDDTs | EXA000001 | 140000.0 | ng/g lipid | PRA000030 | ST004 | 34892 | ... |
| 4431 | 24784 | SAA001905 | CH0000033 | ΣDDTs | EXA000001 | 140000.0 | ng/g lipid | PRA000030 | ST004 | 34892 | ... |
| 4461 | 24785 | SAA001906 | CH0000033 | ΣDDTs | EXA000001 | 130000.0 | ng/g lipid | PRA000030 | ST004 | 34892 | ... |
| 4491 | 24786 | SAA001907 | CH0000033 | ΣDDTs | EXA000001 | 280000.0 | ng/g lipid | PRA000030 | ST004 | 34892 | ... |
5 rows × 37 columns
Sec.4-2 基本統計量
Sec.4-1でデータが抽出できたので、このデータセットの特徴を掴むために、基本統計量の算出をする。基本統計量は、データセットを要約するいくつかの統計量である。今回は、データ総数、平均、標準偏差、第1四分位数、中央値、第3四分位数、最小値、最大値を計算する。
count = ddt_vals.count()
mean = ddt_vals.mean()
std = ddt_vals.std()
q1 = ddt_vals.quantile(0.25)
med = ddt_vals.median()
q3 = ddt_vals.quantile(0.75)
min = ddt_vals.min()
max = ddt_vals.max()
count, mean, std, q1, med, q3, min, max
(24, 102137.5, 81917.43357743236, 41750.0, 70500.0, 140000.0, 9400.0, 280000.0)
上セルでは、基本統計量を各個に計算したが、pandasには基本統計量をまとめて出す関数がある。
avgs = ddt_vals.describe()
avgs
count 24.000000
mean 102137.500000
std 81917.433577
min 9400.000000
25% 41750.000000
50% 70500.000000
75% 140000.000000
max 280000.000000
Name: MeasuredValue, dtype: float64
Sec.4-3 jitter plot
データの散らばり具合を可視化する方法として、全データをプロットする散布図がある (江口: これは散布図ではなくjitter plotでは？）。Pythonではseabornのstripplot関数を使うことで簡単にjitter plotを描くことができる。jitter plotは点どうしが重なってしまう場合にも、点にゆらぎを加えることでそれぞれの天を見やすく記述できる良い可視化手法の1つである。
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
sns.stripplot(x="ChemicalName", y="MeasuredValue", data=data_ddt, color='black', ax=ax) # jitterplotの作図
ax.set_ylabel("ng/g lipid")
plt.show()
Sec.4-4 箱ひげ図
Sec.4-3とは異なり、基本統計量でデータセットを要約したものをプロットするのが、箱ひげ図である。箱ひげ図は、seabornのboxplot関数で描画できる。ここでは先程のjitter plotに箱ひげ図を重ね書きすることで、データの分布をより把握しやすい図を出力している。
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
sns.stripplot(x="ChemicalName", y="MeasuredValue", data=data_lipid[data_lipid["ChemicalName"] == "ΣDDTs"], color='black', ax=ax) # jitterplotの作図
sns.boxplot(x="ChemicalName", y="MeasuredValue", data=data_lipid[data_lipid["ChemicalName"] == "ΣDDTs"], ax=ax) # 箱ひげ図の作図
ax.set_ylabel("ng/g lipid")
plt.show()
Chap.5 スナメリ（Neophocaena phocaenoides）の各化学物質データの可視化
Chap.4ではΣDDTsの値のみを可視化した。それを応用して、data_lipidとdata_wetのデータを可視化する。
Sec.5-1 箱ひげ図（Sec.4-4の応用）
まず、data_lipidのデータについて、Sec.4-4と同じようにグラフを描画する。stripplot・boxplot関数はパラメータxに代入する列名でデータを自動的に分類する機能があるので、ひとまとめにしたデータをそのままパラメータdataに代入すれば良い。
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
sns.stripplot(x="ChemicalName", y="MeasuredValue", data=data_lipid, color='black', ax=ax) # jitterplotの作図
sns.boxplot(x="ChemicalName", y="MeasuredValue", data=data_lipid, ax=ax) # 箱ひげ図の作図
ax.set_ylabel("ng/g lipid")
plt.show()
Sec.5-2 Sec.5-1の改善
上のグラフでは、ΣDDTsのばらつきが大きいので、y軸の間隔が広くなっている。したがって、残りの5種類のグラフが潰れており、見えにくい。そこで、データの値の対数を取ることで、グラフを見やすくする。
まず、対数を取るにあたって、0の対数は取れないので、データセットに0が含まれるか確認する。in文を使うことで、該当する値の有無がTRUE / FALSEで返される。
0 in data_lipid.loc[:, "MeasuredValue"].values
False
上セルでデータセットに０が含まれないことがわかったので、そのまま対数の値を取れる。
data_lipid.loc[:, "MeasuredValue"] = data_lipid.loc[:, "MeasuredValue"].apply(np.log10)
C:\Users\masah\Anaconda3\lib\site-packages\pandas\core\indexing.py:543: SettingWithCopyWarning:
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead
See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
self.obj[item] = s
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
sns.stripplot(x="ChemicalName", y="MeasuredValue", data=data_lipid, color='black', ax=ax) # jitterplotの作図
sns.boxplot(x="ChemicalName", y="MeasuredValue", data=data_lipid, ax=ax) # 箱ひげ図の作図
ax.set_ylabel("log(ng/g) lipid")
plt.show()
data_wetについても同様に、可視化する。まず、データセットに0が含まれる場合、最小値の1/2を代わりに代入する。その後、対数を取り、可視化する。
if 0 in data_wet.loc[:, "MeasuredValue"].values:
data_wet["MeasuredValue"].replace(0, data_wet[data_wet["MeasuredValue"] != 0].loc[:, "MeasuredValue"].values.min() / 2) #最小値の1/2を代入
else:
pass
data_lipid.loc[:, "MeasuredValue"] = data_lipid["MeasuredValue"].apply(np.log10)
data_lipid.loc[:, "Unit"] = "log(ng/g lipid)"
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
sns.stripplot(x="ChemicalName", y="MeasuredValue", data=data_lipid, color='black', ax=ax)
sns.boxplot(x="ChemicalName", y="MeasuredValue", data=data_lipid, ax=ax)
ax.set_ylabel("log(ng/g lipid)")
plt.show()
脚注
1 「マジックコマンド」と呼ばれるJuyterNotebookの機能であり、一般のpythonコードでは利用できないので注意が必要である。
2 tab-separated valuesの略。TAB文字（キーボードの左端にある）でデータを区切ったファイル形式である。表の形をしたデータを保持するのに向いている。類似のファイル形式に、CSV（comma-separated valuesの略。「,（コンマ）」で区切る）がある。