実践ケモインフォマティクス

〜主にケモインフォマティクス, 個人的な備忘録〜

技術レポート再現(フィンガープリントによる化合物の予測根拠可視化 | 三井化学)


これまでの記事では、Fingerprintの可視化構造式描画のカスタマイズについてまとめてきた。
今回は、それらでまとめてきた手法を駆使して、機会学習による化合物の物性予測根拠の可視化を実施する。

具体的には、第43回ケモインフォマティクス討論会での講演フィンガープリントによる化合物の予測根拠可視化(三井化学)の内容を再現する。

講演(予稿)の内容

線形の回帰手法(機械学習を含む)を用いた化合物の物性予測においてその予測根拠を可視化する手法を提案している。具体的には"フィンガープリントを記述子として用いた機械学習モデルにおいて、各bitの寄与を分子構造の広がりを考慮して統合し化合物構造上に可視化する"というもの。
予稿中では、水溶解度のデータを検証用として、予測根拠の可視化手法としての有用性を化学的な知見で評価している。

講演予稿集より引用

f:id:chemoinfo-ottantacinque:20210212014155p:plain

解析条件など

データ

1290分子の水溶解度(LogS)データ(sdfファイル形式)
金子研究室HPのデータセットを利用)

変数

radius=2 の Morganフィンガープリントを使用

前処理

  • 以下の変数を削除
    • 分散が0の変数
    • 互いに相関係数が高い(0.95以上)変数ペアの内、片方の変数
      (金子研モジュール使用)
    • 1090以上の分子で同じ値をとる変数
  • 変数は全て標準化(平均:0, 分散:1)
  • Borutaにより変数選択した(参考:Borutaによる変数選択

解析手法

寄与の計算方法

ある原子 i に対する部分構造の寄与を統合した値 A_i を次の式で算出する。


\displaystyle{
A_i = \sum_{n=1}^N(C_n\times\frac{1}{f_n}\times\frac{1}{x_n})
}


 f:部分構造を構成する原子に一つ以上  i を含んでいるフィンガープリントのリスト
C_n  :各フィンガープリントの寄与
 f_n:分子中に含まれる各部分構造の数(n=1,2,...N
x_n  :各部分構造に含まれる原子数

可視化方法

寄与の大きさに応じて以下の色で原子にハイライトする。

  • 正に寄与している原子:
  • 負に寄与している原子:

前準備

conda install -c conda-forge rdkit
pip install dcekit
pip install borutapy

コード

データセットの準備〜前処理まで

import os
import pandas as pd
import numpy as np

import matplotlib.figure as figure
import matplotlib.pyplot as plt

from rdkit import Chem
from rdkit.Chem import AllChem, Draw
from rdkit.Chem.Draw import rdMolDraw2D
from IPython.display import SVG

from boruta import BorutaPy

from sklearn import svm
from sklearn.preprocessing import StandardScaler
from sklearn import preprocessing
from sklearn.ensemble import RandomForestRegressor as rf
from sklearn.model_selection import train_test_split,GridSearchCV
from sklearn.metrics import r2_score, mean_absolute_error

from dcekit.variable_selection import search_high_rate_of_same_values, search_highly_correlated_variables


suppl = Chem.SDMolSupplier(os.getcwd() + '/../../data/logSdataset1290_2d.sdf')
mols_list = [mol for mol in suppl if mol is not None]

# ECFP4フィンガープリントの作成(Extended Connectivity Fingerprint)
ECFP4 = [AllChem.GetMorganFingerprintAsBitVect(mol, 2,1024) for mol in mols_list]
df = pd.DataFrame(np.array(ECFP4, int))

#溶解度(logS)の値を取得
df['logS'] =  [mol.GetProp('logS') if 'logS' in mol.GetPropNames() else 'None' for mol in mols_list]
df = df.astype(float)

#前処理=====================
#欠損のある変数を削除
df.dropna(axis=1,inplace=True)

#分散が0の変数削除
del_num1 = np.where(df.var() == 0)
df = df.drop(df.columns[del_num1], axis=1)

#変数選択(互いに相関関係にある変数の内一方を削除)   
threshold_of_r = 0.95 #変数選択するときの相関係数の絶対値の閾値
corr_var = search_highly_correlated_variables(df, threshold_of_r)
df.drop(df.columns[corr_var], axis=1, inplace=True)

#同じ値を多くもつ変数の削除
#1090化合物をテストデータとするため、1090化合物以上で同じ値をとる変数を削除
rate_of_same_value = []
threshold_of_rate_of_same_value = 1090
for col in df.columns:
    same_value_number = df[col].value_counts()
    rate_of_same_value.append(float(same_value_number[same_value_number.index[0]] / df.shape[0]))
del_var_num = np.where(np.array(rate_of_same_value) >= threshold_of_rate_of_same_value)
df.drop(df.columns[del_var_num], axis=1, inplace=True)
#前処理終わり=====================


#訓練データと学習データに分割
df_X = df.drop(['logS'],axis=1)
df_y = df.loc[:,['logS']]
X = df_X.values
y = df_y.values

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5, random_state=0)

#標準化
scaler_X = StandardScaler()
scaler_y = StandardScaler()

scaled_X_train = scaler_X.fit_transform(X_train)
scaled_y_train = scaler_y.fit_transform(y_train)
scaled_X_test = scaler_X.transform(X_test)
scaled_y_test = scaler_y.transform(y_test)

borutaによる変数選択

(参考:Borutaによる変数選択

#Borutaによる変数選択====================
#pパーセンタイルの最適化(金子研の方法)
corr_list = []
for n in range(10000):
    shadow_features = np.random.rand(scaled_X_train.shape[0]).T
    corr = np.corrcoef(scaled_X_train, shadow_features,rowvar=False)[-1]
    corr = abs(corr[corr < 0.95])
    corr_list.append(corr.max())
    corr_array = np.array(corr_list)
    perc = 100 * (1-corr_array.max())

#Borutaの実施
feat_selector = BorutaPy(rf(), 
                    n_estimators='auto',  # 特徴量の数に比例して、木の本数を増やす
                    verbose=0,
                    alpha=0.05, # 有意水準
                    max_iter=50, # 試行回数
                    perc = perc, #ランダム生成変数の重要度の何%を基準とするか
                    random_state=0
                    )

feat_selector.fit(scaled_X_train, scaled_y_train)

selected_bit = feat_selector.support_
selected_bit_num = df_X.columns[selected_bit]

scaled_X_train_br = scaled_X_train[:,selected_bit]
scaled_X_test_br = scaled_X_test[:,selected_bit]

print('pパーセンタイル:',round(perc,2))
print('boruta後の変数の数:',scaled_X_train_br.shape[1])
"""output
pパーセンタイル: 81.6
boruta後の変数の数: 136

"""

予測モデル構築と精度検証

#クロスバリデーションによるハイパーパラメータ最適化
params = {'C':2 ** np.arange(-5, 5, 3, dtype=float),
          'epsilon':2 ** np.arange(-10, 0, 3, dtype=float)
         }

model_cv = GridSearchCV(svm.SVR(kernel='linear'), params,cv=5, n_jobs=-1)
model_cv.fit(scaled_X_train_br, scaled_y_train)

model = svm.SVR(kernel='linear', 
                C=model_cv.best_params_['C'], 
                epsilon=model_cv.best_params_['epsilon'])

#学習と予測精度の調査
model.fit(scaled_X_train_br, scaled_y_train)

scaled_pred_y_train = np.ndarray.flatten(model.predict(scaled_X_train_br))
scaled_pred_y_test = np.ndarray.flatten(model.predict(scaled_X_test_br))
pred_y_train = scaler_y.inverse_transform(scaled_pred_y_train)
pred_y_test = scaler_y.inverse_transform(scaled_pred_y_test)

print('r2_train:{}'.format(round(r2_score(y_train, pred_y_train),3)))
print('r2_test:{}'.format(round(r2_score(y_test, pred_y_test),3)))
print('MAE_train:{}'.format(round(mean_absolute_error(y_train, pred_y_train)),3))
print('MAE_test:{}'.format(round(mean_absolute_error(y_test, pred_y_test)),3))
"""output
r2_train:0.729
r2_test:0.615
MAE_train:0.84
MAE_test:0.984
"""

plt.figure(figsize=(6,6))
plt.scatter(y_train, pred_y_train,color='blue',alpha=0.3,label='training data')
plt.scatter(y_test, pred_y_test, color='red',alpha=0.3,label='test data')
y_max_ = max(y_train.max(), pred_y_train.max(),y_test.max(), pred_y_test.max())
y_min_ = min(y_train.min(), pred_y_train.min(),y_test.min(), pred_y_test.min())
y_max = y_max_ + (y_max_ - y_min_) * 0.1
y_min = y_min_ - (y_max_ - y_min_) * 0.1

plt.plot([y_min , y_max],[y_min, y_max], 'k-')

plt.ylim(y_min, y_max)
plt.xlim(y_min, y_max)
plt.xlabel('logS(observed)',fontsize=20)
plt.ylabel('logS(predicted)',fontsize=20)
plt.legend(loc='best',fontsize=15)
plt.title('yyplot(logS)',fontsize=20)

plt.savefig('yyplot.png')

plt.show()

予測精度は以下の通りで、予稿集と同程度になった。

r2(決定係数) MAE(誤差平均)
train data 0.729 0.84
test data 0.615 0.984


f:id:chemoinfo-ottantacinque:20210212014819p:plain

回帰係数の可視化

#全データを使ってモデルを再構築
scaler_X = StandardScaler()
scaler_y = StandardScaler()
scaled_X = scaler_X.fit_transform(X[:,selected_bit])
scaled_y = scaler_y.fit_transform(y)
model.fit(scaled_X, scaled_y)
pred_y = model.predict(scaled_X)

#回帰係数の可視化
selected_bit = feat_selector.support_
selected_bit_num = df_X.columns[selected_bit]

#回帰係数の取得
regression_coef = model.coef_[0]

#棒グラフの作成
x = np.arange(len(selected_bit_num))
width = 0.95

fig, ax = plt.subplots(figsize=(20,5))
bar = ax.bar(x, regression_coef, width)
ax.set_xticks(x)
ax.set_xticklabels(selected_bit_num,fontsize=5)
ax.set_ylabel('Regression Coefficient',fontsize=20)
ax.set_title('Regression Coefficient each bits',fontsize=25)
plt.grid(linestyle='--',axis='y')
plt.savefig('Regression Coefficient.png',bbox_inches='tight')
plt.show()
f:id:chemoinfo-ottantacinque:20210212014902p:plain

寄与の計算と可視化

寄与の計算(おさらい)

ある原子i に対する部分構造の寄与を統合した値A_iを次の式で算出する。


\displaystyle{
A_i = \sum_{n=1}^N(C_n\times\frac{1}{f_n}\times\frac{1}{x_n})
}


 f:部分構造を構成する原子に一つ以上  i を含んでいるフィンガープリントのリスト
C_n  :各フィンガープリントの寄与
 f_n:分子中に含まれる各部分構造の数(n=1,2,...N
x_n  :各部分構造に含まれる原子数

可視化方法(おさらい)

寄与の大きさに応じて以下の色で原子にハイライトする。

  • 正に寄与している原子:
  • 負に寄与している原子:

化合物の描画のカスタマイズについては、以下記事を参照
rdMolDraw2Dモジュールを使って構造式描画をカスタマイズ

可視化

適当に一分子ピックアップして可視化してみる。

#ビットと回帰係数の辞書を作っておく
bit_coef = dict(zip(selected_bit_num, regression_coef))

#適当な分子をピックアップ
mol = mols_list[190]
print(Chem.MolToSmiles(mol))
# >>> c1ccc2ocnc2c1

#フィンガープリントを計算
bitI_morgan = {}
fp_morgan = AllChem.GetMorganFingerprintAsBitVect(mol, 2, 1024, bitInfo=bitI_morgan)

#寄与のあるビットの抜き出し
bit_list = list(set(selected_bit_num)&set(bitI_morgan.keys()))
print('寄与のあるビットの数:',len(bit_list))
# >>> 寄与のあるビットの数: 9

#寄与を格納する配列を生成
Ai_list = np.zeros(mol.GetNumAtoms())

for bit in bit_list:

    #フィンガープリントの寄与
    Cn = bit_coef[bit]

    #分子中に含まれる部分構造の数
    fn = len(bitI_morgan[bit])

    for part in bitI_morgan[bit]:
        
        #fingerprintのradiusが0の時はfingerprintの中心原子を抜き出し、寄与を加算する。
        if part[1]==0:
            i = part[0]
            xn = 1
            Ai_list[i] += Cn / fn / xn
        
        #fingerprintのradiouが0以上の時は、該当する原子のリストを抜き出し、寄与を加算する
        else:
            amap={}
            env = Chem.FindAtomEnvironmentOfRadiusN(mol,
                                                    radius=part[1],
                                                    rootedAtAtom=part[0])
            submol=Chem.PathToSubmol(mol,env,atomMap=amap)

            #各部分構造に含まれる原子数
            xn = len(list(amap.keys()))
            
            for i in amap.keys():
                Ai_list[i] += Cn / fn / xn
                
#正規化っぽいことをする
#絶対値の最大値が0.5になるように調整する
Ai_list = Ai_list / abs(Ai_list).max() * 0.5
print('寄与の値')
print(Ai_list)
"""output
寄与の値
[-0.5        -0.44090249 -0.06607348 -0.04066189 -0.04066189 -0.07123291
  0.14772543 -0.04021274  0.08323835]
"""

# ハイライトする原子と色,円の半径を設定
atoms = [i for i in range(len(Ai_list))]
atom_colors = dict()
for i in atoms:
    if Ai_list[i] > 0:
        atom_colors[i] = (1,1-Ai_list[i],1-Ai_list[i])
    else:
        atom_colors[i] = (1+Ai_list[i],1+Ai_list[i],1)

# コンテナの準備
view = rdMolDraw2D.MolDraw2DSVG(300,350)
tm = rdMolDraw2D.PrepareMolForDrawing(mol)

view.DrawMolecule(tm,
                  highlightAtoms=atoms,
                  highlightAtomColors=atom_colors,
                  highlightBonds=[],
                  highlightBondColors={})

#ファイナライズ、保存、描画
view.FinishDrawing()
svg = view.GetDrawingText()
with open('highlighted_sample.svg', 'w') as f:
    f.write(svg)
SVG(svg)

結果

f:id:chemoinfo-ottantacinque:20210212015233p:plain



他にも何分子か描画してみた。

f:id:chemoinfo-ottantacinque:20210212020508p:plain

分極しやすく水との親和性が高いOやNは赤色でハイライトされていることが多く、逆に分極が少なく水との親和性の引く芳香環や炭素鎖は青にハイライトされていることが多い。
化学的な知見とも合うので可視化手法として有用と言える。

寄与の正負を問わないのであればランダムフォレスト系の手法にも適用できそう。

参考

rdMolDraw2Dモジュールを使って構造式描画をカスタマイズ


前回の記事では、原子ごとに類似度の寄与率を可視化することができるメソッドSimilarityMapsを使って、TPSAなどの分子記述子に対する各原子の寄与を可視化する方法についてまとめた。

今回は類似度マップ(SmililarityMaps)を用いずに、描画方法をカスタマイズ(構造式にハイライトをつけたり、原子をインデックス表示にしたり etc...)についてまとめた。
(本記事は「化学の新しいカタチ」の内容を簡潔にまとめたものです。より詳しい内容は、そちらに載っています)



描画の流れ

描画を色々カスタマイズする場合には、rdMolDraw2Dモジュールを使用する。 このモジュールを使用する場合には、以下の手順をふむ必要がある。

参考HPより

  1. コンテナとなるオブジェクトの作成
  2. コンテナの細かい描画設定
  3. コンテナに描画したい分子を追加
  4. コンテナファイナライズして書き込んだデータを取り出す
  5. データを描画,またはファイルとして保存

コード

モジュールのインストール

conda install -c conda-forge rdkit

分子の準備

from rdkit import rdBase, Chem
from rdkit.Chem import AllChem, Draw
from rdkit.Chem.Draw import rdMolDraw2D
from IPython.display import SVG
from matplotlib.colors import ColorConverter

mol = Chem.MolFromSmiles('Cn1c(=O)c2c(ncn2C)n(C)c1=O')
smiles = Chem.MolToSmiles(mol)

# とりあえずDraw.MolToImage(mol)で描画
img = Draw.MolToImage(mol)
img.save('caffeine_molecule.png',bbox_inches='tight')
img
f:id:chemoinfo-ottantacinque:20210211140758p:plain

基本操作

# 1. コンテナを準備(コンテナサイズは300x300とする)
view = rdMolDraw2D.MolDraw2DSVG(300, 300)

# 2.描画設定
option = view.drawOptions()
option.circleAtoms=False
option.continuousHighlight=False

# 3.コンテナに分子を追加
#単一分子の追加
tm = rdMolDraw2D.PrepareMolForDrawing(mol)
#tm = rdMolDraw2D.DrawMolecules(mols)
#複数分子の追加
#tm = rdMolDraw2D.DrawMolecules(mols)

view.DrawMolecule(tm)

# 4.コンテナファイナライズして書き込んだデータを取り出す
view.FinishDrawing()

# 5.6コンテナの描画と保存
svg = view.GetDrawingText()

#保存
with open('caffein_test.svg', 'w') as f:
    f.write(svg)

#描画
SVG(svg)


f:id:chemoinfo-ottantacinque:20210211140654p:plain


何も設定を変えていないのでDraw.MolToImage(mol)使用時と同じ絵になる。

それぞれのメソッドのオプションなど

新しい化学のカタチより転記)

MolDraw2DSVG.DrawMolecule(mol)などのDrawメソッド


オプション 説明
highlightAtoms ハイライトする原子のリスト
highlightAtomColors ハイライトする原子の色
highlightAtomRadii ハイライトする原子の半径を設定
highlightBonds ハイライトする結合のリスト
highlightBondColors ハイライトする結合の色
legend(s) 凡例


MolDraw2DSVG.drawOptions()


オプション

オプション 説明
atomLabelDeuteriumTritium 重水素をD, トリチウムをTで表記
atomLabels 原子Idで指定した原子のラベルを設定
additionalAtomLabelPadding 元素記号と結合との間にどの程度マージンを取るか
circleAtoms ハイライトするアトムを円で表示するかどうか
デフォルト:True
continuousHighlight ハイライト領域を塗りつぶすかどうか
デフォルト:True
fillHighlights ハイライトするアトムの円を塗りつぶすかどうか
デフォルト:True
padding 描画領域のうちのマージンの割合(0〜1)
legendFontSize 凡例のフォントサイズ(px単位)
multipleBondOffset 多重結合の2つ目以降の線をどの程度離して描画するか
オングストローム単位)
addAtomIndices 原子Idを構造式中に表示
デフォルト:False
addBondIndices 結合Idを構造式中に表示
デフォルト:False
annotationFontScale 原子や結合に付記する注釈のフォントサイズ(割合指定)
デフォルト:0.75
addStereoAnnotation 絶対配置R/Sや幾何異性E/Zの表示有無
デフォルト:False


メソッド

メソッド 説明
useBWAtomPalette 構造式を白黒で描画する
useDefaultAtomPalette 構造式描画に用いる色をデフォルトに戻す
updateAtomPalette 構造式描画に用いる色を,原子番号とRGBのタプルを組とする辞書形式で指定する
setBackgroundColour キャンバス色の変更する
setHighlightColour ハイライト領域のデフォルト色を変更する

   

描画をカスタマイズする

原子のラベルとハイライトする色を変更してみる。その他にも細々設定する(コード内に説明付記)
(ハイライトする色はmatplotlibをつかって生成したRGBカラーコードを使う)

#RGBカラーコードを生成
color1 = ColorConverter().to_rgb('cyan')
color2 = ColorConverter().to_rgb('skyblue')
color3 = ColorConverter().to_rgb('yellow')
print('cyan:',color1)
print('skyblue:',color2)
print('yellow:',color3)
"""output
cyan: (0.0, 1.0, 1.0)
skyblue: (0.5294117647058824, 0.807843137254902, 0.9215686274509803)
yellow: (1.0, 0.0, 0.0)
"""

# ハイライトする原子と色,円の半径を設定
atoms = [6,8]
atom_colors = {6: color1, 8: color2}
radii = {6: 0.5, 8: 0.5}

# ハイライトする結合と色を指定
bonds = [1,2,3,4,5]
bonds_colors = {1:color1,2:color2,3:color3,4:color3,5:color3}

# コンテナの準備
view = rdMolDraw2D.MolDraw2DSVG(300,350)
tm = rdMolDraw2D.PrepareMolForDrawing(mol)

# 描画オプションの設定
#フォントサイズをやや小さくする
view.SetFontSize(0.9 * view.FontSize())

option = view.drawOptions()

#原子ラベルを設定する
option.atomLabels[6] = 'N1'
option.atomLabels[8] = 'N2'

#二重結合の幅
option.multipleBondOffset=0.07

#描画領域のマージン
option.padding=0.11

#レジェンドのフォントサイズ
option.legendFontSize=20

# コンテナに分子を登録
view.DrawMolecule(tm,
                  highlightAtoms=atoms,
                  highlightAtomColors=atom_colors,
                  highlightAtomRadii=radii,
                  highlightBonds=bonds,
                  highlightBondColors=bonds_colors,
                  legend='caffeine')


#ファイナライズ、保存、描画
view.FinishDrawing()
svg = view.GetDrawingText()
with open('caffein_test.svg', 'w') as f:
    f.write(svg)
SVG(svg)
f:id:chemoinfo-ottantacinque:20210211140638p:plain

参考


分子記述子への各原子の寄与率を可視化する


原子ごとに類似度の寄与率を可視化することができるメソッドSimilarityMapsを使って、TPSAなどの分子記述子に対する各原子の寄与を可視化する方法についてまとめた。
(本記事は「化学の新しいカタチ」の内容を簡潔にまとめたものです。より詳しい内容は、そちらに載っています。


準備

モジュールのインストール

conda install -c conda-forge rdkit

コード

分子の準備

練習用分子としてカフェイン分子を使用する。

import numpy as np
import matplotlib.pyplot as plt

from rdkit import rdBase, Chem
from rdkit.Chem import AllChem, Draw, rdMolDescriptors
from rdkit.Chem.Draw import SimilarityMaps

mol = Chem.MolFromSmiles('Cn1c(=O)c2c(ncn2C)n(C)c1=O')
smiles = Chem.MolToSmiles(mol)

img = Draw.MolToImage(mol)
img.save('caffeine_molecule.png',bbox_inches='tight')
img
f:id:chemoinfo-ottantacinque:20210210232623p:plain

MolLogPへの寄与の計算と可視化

#_CalcCrippenContribs:各原子ごとに(MolLogPへの寄与,MolMRへの寄与)を格納したタプルが得られる。
crippen = rdMolDescriptors._CalcCrippenContribs(mol)

mol_log = []
mol_mr = []
for x, y in crippen:
    mol_log.append(x)
    mol_mr.append(y)

fig = SimilarityMaps.GetSimilarityMapFromWeights(mol, 
                                                 mol_log, 
                                                 colorMap='PRGn',
                                                 alpha=0.5)
fig.savefig('MollogP.png',bbox_inches='tight')
f:id:chemoinfo-ottantacinque:20210210232949p:plain

棒グラフでの寄与可視化

定量的に分かるように棒グラフでも可視化してみた。

#ついでに各原子の寄与を棒グラフでも表現してみる
atom_list = ['{}{}'.format(atom.GetSymbol(),atom.GetIdx()) for atom in mol.GetAtoms()]

#棒グラフの作成
x = np.arange(len(atom_list))
width = 0.8

fig, ax = plt.subplots()
bar = ax.bar(x, mol_log, width)
ax.set_xticks(x)
ax.set_xticklabels(atom_list)
ax.set_ylabel('Contribution to MollogP',fontsize=12)
plt.grid(linestyle='--',axis='y')
plt.savefig('bar_graph_mollogp.png')
plt.show()
f:id:chemoinfo-ottantacinque:20210210233027p:plain

TPSAへの寄与の計算と可視化

#_CalcTPSAContribs(mol)::各原子ごとにのTPSAへの寄与を格納したタプルが得られる。
tpsa = rdMolDescriptors._CalcTPSAContribs(mol)
fig = SimilarityMaps.GetSimilarityMapFromWeights(mol, 
                                                 weights=tpsa)
fig.savefig('tpsa.png',bbox_inches='tight')
f:id:chemoinfo-ottantacinque:20210210233108p:plain

棒グラフでの寄与可視化

MolLogPと同様に棒グラフでも可視化してみた。

#ついでに各原子の寄与を棒グラフでも表現してみる
atom_list = ['{}{}'.format(atom.GetSymbol(),atom.GetIdx()) for atom in mol.GetAtoms()]

#棒グラフの作成
x = np.arange(len(atom_list))
width = 0.8

fig, ax = plt.subplots()
bar = ax.bar(x, tpsa, width)
ax.set_xticks(x)
ax.set_xticklabels(atom_list)
ax.set_ylabel('Contribution to TPSA',fontsize=12)
plt.grid(linestyle='--',axis='y')
plt.savefig('bar_graph_tpsa.png')
plt.show()
f:id:chemoinfo-ottantacinque:20210210233207p:plain

参考


Fingerprintの可視化について


フィンガープリント(Fingerprint)をRDkitで可視化する方法についてまとめた。 (本記事は「化学の新しいカタチ」の内容を簡潔にまとめたものです。より詳しい内容はそちらに載っています)


フィンガープリントについて

フィンガープリントとは、化合物中に特定の部分構造が含まれるかを 「0 or 1」のビットで表したもの(参考)。
化合物の物性予測の際の変数としてよく用いられるが、あるビットが重要だと分かったとしても、化学的な知見を得るためにはビットと構造を紐づけて可視化することが必須となる。
(フィンガープリントへの変換方法まとめはこちら

準備

モジュールのインストール

conda install -c conda-forge rdkit

分子の準備

題材としてカフェイン分子を利用する。

from rdkit import rdBase, Chem
from rdkit.Chem import AllChem, Draw

mol = Chem.MolFromSmiles('Cn1c(=O)c2c(ncn2C)n(C)c1=O')
smiles = Chem.MolToSmiles(mol)

# とりあえず描画したい時
img = Draw.MolToImage(mol)
img.save('caffeine_molecule.png',bbox_inches='tight')
img
f:id:chemoinfo-ottantacinque:20210209224827p:plain

可視化

可視化に対応している、 「Morganフィンガープリント」 および、 「RDkitフィンガープリント」両方で可視化してみる。

Morganフィンガープリント

### Morganフィンガープリント
#ビットの立っている位置と,関連する部分構造を記録する辞書(Fingerprintの可視化に必要)
bitI_morgan = {}
fp_morgan = AllChem.GetMorganFingerprintAsBitVect(mol, 2, 1024, bitInfo=bitI_morgan)

#ビットの立っている一と関連する部分構造が辞書形式で格納される。
print('Morgan Fingerprintの数:',fp_morgan.GetNumBits())
print('Onのビットの数:',fp_morgan.GetNumOnBits())
"""output
Morgan Fingerprintの数: 1024
Onのビットの数: 24
"""

#それぞれのビットが立っているか調べる
for key in list(bitI_morgan.keys())[:5]:
    print('ビット:{};, 関連する部分構造:{}'.format(key, bitI_morgan[key]))
"""output
ビット:0;, 関連する部分構造:((10, 2),)
ビット:33;, 関連する部分構造:((0, 0), (9, 0), (11, 0), (5, 2))
ビット:121;, 関連する部分構造:((0, 1), (9, 1), (11, 1))
ビット:179;, 関連する部分構造:((6, 2),)
ビット:234;, 関連する部分構造:((1, 2),)
"""

#単一のビット可視化にはDrawMorganbitを使う。
img = Draw.DrawMorganBit(mol, list(bitI_morgan.keys())[1], bitI_morgan)
img.save('bit_single_molecule_morgan.png',bbox_inches='tight')

#複数のビットの可視化
morgan_turples = ((mol, bit, bitI_morgan) for bit in list(bitI_morgan.keys()))
img = Draw.DrawMorganBits(morgan_turples, molsPerRow=5, 
                    legends=['bit: '+str(x) for x in list(bitI_morgan.keys())])
img.save('bit_several_molecule_morgan.png',bbox_inches='tight')
img


単一ビットの描画結果

f:id:chemoinfo-ottantacinque:20210209224846p:plain



複数ビットの描画結果

f:id:chemoinfo-ottantacinque:20210209224901p:plain

MorganFingerprint可視化方法のルール

(参考HPより)

  • 原子の位置は分子内の位置と同じように描画
  • ビット情報の中心原子は薄い青丸で表現
  • 芳香族原子は黄色い丸で表現
  • 環状の脂肪族原子は灰色で表現
  • 直接にはフィンガープリントに含まれないが,原子の結合タイプ決定に影響する部分を薄い灰色で表現

RDkitフィンガープリント

#ビットの立っている位置と,関連する部分構造を記録する辞書(Fingerprintの可視化に必要)
bitI_rdkit = {}
fp_rdkit = Chem.RDKFingerprint(mol, fpSize=1024, bitInfo=bitI_rdkit)

#ビットの立っている一と関連する部分構造が辞書形式で格納される。
print('RDkit Fingerprintの数:',fp_rdkit.GetNumBits())
print('Onのビットの数:',fp_rdkit.GetNumOnBits())
"""output
RDkit Fingerprintの数: 1024
Onのビットの数: 626
"""

#それぞれのビットが立っているか調べる
for key in list(bitI_rdkit.keys())[:5]:
    print('ビット:{};, 関連する部分構造:{}'.format(key, bitI_rdkit[key]))
"""output
ビット:0;, 関連する部分構造:[[1, 3, 4], [3, 4, 9], [3, 4, 5]]
ビット:1;, 関連する部分構造:[[0, 13, 11, 9, 4, 14, 8], [0, 13, 1, 3, 4, 9, 2], [0, 13, 1, 3, 4, 5, 2], [1, 3, 14, 8, 13, 11, 10], [4, 9, 11, 10, 5, 14, 8]]
ビット:2;, 関連する部分構造:[[3], [4]]
ビット:3;, 関連する部分構造:[[1], [5], [6], [7], [9], [11], [13], [14]]
ビット:4;, 関連する部分構造:[[0, 13, 11, 10, 9, 4], [0, 13, 11, 10, 1, 3]]
"""
        
#単一のビット可視化にはDrawRDkitbitを使う。
img = Draw.DrawRDKitBit(mol, list(bitI_rdkit.keys())[0], bitI_rdkit)
img.save('bit_single_molecule_RDKit.png',bbox_inches='tight')

#複数のビットの可視化
rdkit_turples = ((mol, bit, bitI_rdkit) for bit in list(bitI_rdkit.keys())[:12])
img = Draw.DrawRDKitBits(rdkit_turples, molsPerRow=5, 
                   legends=['bit: '+str(x) for x in list(bitI_rdkit.keys())[:12]])
img.save('bit_several_molecule_rdkit.png',bbox_inches='tight')
img


単一ビットの描画結果

f:id:chemoinfo-ottantacinque:20210209224915p:plain



複数ビットの描画結果

f:id:chemoinfo-ottantacinque:20210209224929p:plain

RDKitフィンガープリント可視化のルール

参考HPより)

  • 原子の位置は分子内の位置と同じように描画
  • 描画された結合は全てビット情報の一部
  • 芳香族原子は黄色い丸で表現
  • 脂肪族原子はデフォルト設定では特に決まった表記なし

参考


主成分分析(概要とpython実装)


データ可視化手法としてよく用いられる主成分分析についてまとめた。

主成分分析について

概要

  • 教師なし手法の一つであり、PCA(Principal Component Analysis)と呼ばれる。
  • 多次元のデータを低次元化する手法であり、データの可視化によく用いられる。 (低次元化し2次元にマッピングする事でデータの可視化が可能)
  • データセットの情報量損失を最小限にしつつ低次元化することが可能。
  • データのノイズを寄与率の小さい主成分として除くことが可能。
  • PCA後の変数(主成分)は互いに相関0なので、これを用いて回帰分析すると回帰係数の値が安定になる。(この手法は主成分回帰(PCRと呼ばれる)

方法(ざっくり)

  1. データのばらつきが最も大きくなる方向に新しく軸を設定する。
  2. 1.で設定した軸と直行し、かつデータのばらつきが最も大きくなる方向に2軸目を設定する。
  3. 1,2を繰り返す。
  4. 1,2の過程では、それぞれの主成分がデータ全体をどの程度表せているかの指標(寄与率)が得られる。

主成分分析で得られる指標

  • 固有値:各主成分に対応した固有値であり、主成分の分散に対応する。
  • 寄与率:ある主成分の固有値(分散)が表す情報が、データの全ての情報(分散)の何割を表しているかの指標。
    (言い換えると、ある主成分がデータの何割を表せられているかの指標のこと)
  • 累積寄与率:各主成分の寄与率を大きい順に足したもの。そこまでの主成分でデータの何割を表現できているかを表す。
    (累積寄与率が70〜80%がカバーできていればOKと判断する場合が多い)

注意点

  • 主成分分析の前に変数を標準化する必要がある。

コード

データの準備

水溶解度のデータセットを利用する(金子研究HP から取得可能)

import os
import pandas as pd
import numpy as np
import math

import matplotlib.figure as figure
import matplotlib.pyplot as plt

from rdkit import Chem
from rdkit.ML.Descriptors import MoleculeDescriptors

#データの読み込み
suppl = Chem.SDMolSupplier('logSdataset1290_2d.sdf')
mols = [mol for mol in suppl if mol is not None]

#RDkit記述子の計算
descriptor_names = [descriptor_name[0] for descriptor_name in Chem.Descriptors._descList]
descriptor_calculation = MoleculeDescriptors.MolecularDescriptorCalculator(descriptor_names)
RDkit = [descriptor_calculation.CalcDescriptors(mol) for mol in mols]

#データフレーム作成
df = pd.DataFrame(RDkit, columns = descriptor_names)

df['logS'] =  [mol.GetProp('logS') if 'logS' in mol.GetPropNames() else 'None' for mol in mols]
df = df.astype(float)

#サンプル数を少なくする。
df.sort_values('logS',ascending=False,inplace=True)
df = df.iloc[::10,:]

#相関係数が高い10変数を選抜する。
df_corr = df.corr()
df_corr.sort_values('logS',ascending=False,inplace=True)
selected_features = df_corr.index[0:11]
print(selected_features)

df = df.loc[:,selected_features]

df_X = df.drop('logS',axis=1)
df_y = df['logS']

主成分分析の実施とデータの可視化

from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA #主成分分析器
import matplotlib.ticker as ticker

#標準化
scaler_X = StandardScaler()
X_scaled = scaler_X.fit_transform(df_X)

#主成分分析の実行
pca = PCA()
feature = pca.fit_transform(X_scaled)
df_pca = pd.DataFrame(feature, index=df_X.index, 
                      columns=["PC{}".format(x+1) for x in range(feature.shape[1])])

# データを主成分空間に写像
plt.figure(figsize=(5,4))
plt.scatter(df_pca.iloc[:, 0],df_pca.iloc[:, 1], alpha=0.8,
           c=df_y, vmin=df_y.min(),vmax=df_y.max())
plt.colorbar()

#図にテキストを表示する場合
# for x, y, name in zip(df_pca.iloc[:, 0], df_pca.iloc[:, 1], df_pca.index.tolist()):
#     plt.text(x, y, name)
plt.grid()
plt.xlabel("PC1")
plt.ylabel("PC2")
plt.show()

# 累積寄与率を図示する
plt.figure(figsize=(4,4))
plt.gca().get_xaxis().set_major_locator(ticker.MaxNLocator(integer=True))
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主成分でのマッピング結果
f:id:chemoinfo-ottantacinque:20210208212520p:plain
累積寄与率の表示
f:id:chemoinfo-ottantacinque:20210208212828p:plain

詳細な計算式など

計算式(ざっくり)

以下数式は金子研究室HP記載内容を参考にしています。
より詳細な計算は金子研究室HP主成分分析についてがおすすめです。
また、数式に間違いがあればコメントくだされば助かります。

主成分の計算

以下の様なデータを使った、2変数のPCAを考える。

変数A 変数B
サンプル1  x_A^{(1)}  x_B^{(1)}
サンプル2  x_A^{(2)}  x_B^{(2)}
サンプルn  x_A^{(n)}  x_B^{(n)}


ここで、


\displaystyle{
\boldsymbol{X_A} = 
\begin{bmatrix}
x_A^{(1)} \\
x_A^{(2)} \\
⋮\\
x_A^{(n)} 
\end{bmatrix}  
\qquad
\boldsymbol{X_B} = 
\begin{bmatrix}
x_B^{(1)} \\
x_B^{(2)} \\
⋮\\
x_B^{(n)} 
\end{bmatrix}  
}

と表記する。



i主成分を t_i、第i 主成分に対応するj番目の変数の重み(ローディング)をp_i^{(j)} とすると、


\displaystyle{
t_1 = \boldsymbol{X}_Ap_1^{(1)} + \boldsymbol{X}_Bp_1^{(2)} \\
t_2 = \boldsymbol{X}_Ap_2^{(1)} + \boldsymbol{X}_Bp_2^{(2)} 
}


これを行列で表すと、


\displaystyle{
\begin{bmatrix}
t_1^{(1)}&t_2^{(1)} \\
t_1^{(2)}&t_2^{(2)} \\
⋮&⋮\\
t_1^{(n)}&t_2^{(n)}
\end{bmatrix}
=
\begin{bmatrix}
x_A^{(1)}&x_B^{(1)} \\
x_A^{(2)}&x_B^{(2)} \\
⋮&⋮\\
x_A^{(n)}&x_B^{(n)}
\end{bmatrix}
\begin{bmatrix}
p_1^{(1)}&p_2^{(1)} \\
p_1^{(2)}&p_2^{(2)} \\
\end{bmatrix}
}


 x_i^{(k)} :サンプル k における、 i 番目の変数の値
 t_i^{(k)} :サンプル k における、第 i 主成分の値
 p_i^{(j)} :第 i 主成分に対応する、 j 番目の変数の重み



主成分分析では、データセットのばらつき(分散)が最大の方向を主成分軸とする。
ここで、元のデータの平均が0であれば、軸をとりなおした後のデータの平均も0となることを利用すると、 軸変更後のデータの分散が最大となるためには、主成分の値の二乗和が最大になれば良い
(第1主成分について考えると)つまり、


\displaystyle{
S = \sum_{i=1}^{n}(t_1^{(i)})^2
}


ここで、ローディング(重み)によって、主成分が変わらないように、以下の制約を設ける。


\displaystyle{
(p_1^{(1)})^2 + (p_1^{(2)})^2 = 1
}



上記の制約と、Lagrangeの未定乗数法を利用してSが最大となる条件を求めると、


\displaystyle{
(\boldsymbol{X}^T\boldsymbol{X} - \lambda\boldsymbol{E})\boldsymbol{p_1} = 0
}



\displaystyle{
\boldsymbol{X} = 
\begin{bmatrix}
x_A^{(1)} \\
x_A^{(2)} \\
⋮\\
x_A^{(n)} 
\end{bmatrix} ,
\quad

\boldsymbol{E} = 
\begin{bmatrix}
1&0\\
0&1\\
\end{bmatrix}  ,
\quad

\boldsymbol{p_1} = 
\begin{bmatrix}
p_1^{(1)}\\
p_1^{(2)}\\
\end{bmatrix}  
\qquad
}


となる。

この式が、 p_1^{(1)}=p_1^{(2)}=0 以外の解を持つためには、(\boldsymbol{X}^T\boldsymbol{X} - \lambda\boldsymbol{E})行列式が0である必要がある。

つまりこの問題は、 \lambda 固有値\boldsymbol{p_1} 及び、\boldsymbol{p_2} 固有ベクトルとする固有値問題と考えることができる。これによって\boldsymbol{p_1} , \boldsymbol{p_2} を求め、対応する主成分を計算する。


(累積)寄与率の計算

 i 主成分に対応する固有値\lambda_i  は、その主成分の二乗和(=分散)に等しいので、固有値\lambda_i  を第 i 主成分の持つ情報量と仮定すると、寄与率 c_i 以下の式で求められる。


\displaystyle{
c_i = \frac{\lambda_i}{\sum_{j=1}^m\lambda_j}
}


 m :全ての主成分の数)

また、 i 主成分までの寄与率の和は、第  i 主成分までの累積寄与率 C_i と呼ばれる。


\displaystyle{
C_i = \frac{\sum_{j=1}^i\lambda_j}{\sum_{j=1}^m\lambda_j}
}

計算コード

2変数に絞って、sklearnを使わずに主成分分析してみる。

データの準備

#データの抜き出しと標準化
df_X2 = df_X.iloc[:,:2]
scaler_X = StandardScaler()
X_scaled2 = scaler_X.fit_transform(df_X2)

# 変換前のデータを可視化する。
plt.figure(figsize=(5,5))
plt.scatter(X_scaled2[:,0],X_scaled2[:,1], alpha=0.8)
plt.xlabel("X1")
plt.ylabel("X2")
plt.grid()
plt.show()

主成分分析前のデータはこんな感じ

f:id:chemoinfo-ottantacinque:20210208213519p:plain

主成分分析実行

\boldsymbol{X}^T\boldsymbol{X}固有ベクトル(Eigenvector)と固有値(eigenvalue)を求める。

import numpy as np
import numpy.linalg as LA

#固有値、固有ベクトルの計算
a= np.dot(X_scaled2.T, X_scaled2)
eigenvalue, eigenvector = LA.eig(a)

print('固有値:\n',eigenvalue)
print('----')
print('固有ベクトル:\n',eigenvector)
"""output
固有値:
 [246.22143135  11.77856865]
----
固有ベクトル:
 [[ 0.70710678 -0.70710678]
 [ 0.70710678  0.70710678]]
"""

#主成分の計算
pca = np.dot(X_scaled2, eigenvector)

#寄与率の計算
eigenvalue_sum = eigenvalue.sum()

eigenvalue1 = eigenvalue[0]/eigenvalue_sum
eigenvalue2 = eigenvalue[1]/eigenvalue_sum
print('----')
print('寄与率(第1主成分):',eigenvalue1)
print('寄与率(第2主成分):',eigenvalue2)
"""output
寄与率(第1主成分): 0.9543466331268278
寄与率(第2主成分): 0.04565336687317221
"""

# 変換のデータ
plt.figure(figsize=(5,5))
plt.scatter(pca[:,0],pca[:,1], alpha=0.8)
plt.xlabel("pc1")
plt.ylabel("pc2")
plt.grid()
plt.savefig('after_pca.png')
plt.show()

主成分分析後のマッピング

f:id:chemoinfo-ottantacinque:20210208214614p:plain

参考

Borutaによる変数選択



変数選択は精度の高い予測モデルの構築において非常に重要といえる。本記事では、変数選択手法の一つであるBorutaについてまとめた。

Borutaについて

  • ランダムフォレスト(RF)の変数重要度に基づく変数選択方法
  • 目的変数と関係のない適当な特徴量(shadow features)と、オリジナル変数の重要度を比較することで説明変数を選択する
  • 詳しい説明は金子研HP参照

詳しい手順

  1. 特徴量をコピーして新しい特徴量(shadow features)を作りだす。(複数)
  2. 新しく作成した特徴量(shadow features)の値をランダムに入れ替えて適当な(目的変数と関係のない)特徴量にする。
  3. 全部の特徴量を合わせてRandom forestや勾配ブースティングなどで分類を行い、各特徴量の変数重要度を取得する。
  4. オリジナルの特徴量の変数重要度が Shadow featuresの変数重要度よりも大きければそれをカウントし(hitとし)、再度学習を行う。
    (正確にはShadow features の変数重要度の上位p%。後述)
  5. 1~4を繰り返し、統計的に有意な特徴量をのみを選択する。
    (両側検定でオリジナルの説明変数がshadow featuresと比較して重要か判定する)

pパーセンタイルについて

pパーセンタイル:shadow features の変数重要度の何%を基準にするかの指標。 例えば、p=100にすると、shadow features のうち最も高い変数重要度を基準にhit判定する(下図参照)

p=100の場合のhit判定イメージ
f:id:chemoinfo-ottantacinque:20210203222921p:plain

pパーセンタイルの設定方法

以下は金子研究室HPより。

p=100 とすると説明変数が削除されすぎてモデルの推定性能が低下することがある。(特に、サンプルが少ない時には目的変数と関係性があるshadow features がたまたま生成することもある

上記を考慮して、以下の手順で p を設定するのがベターとのこと。

変数をランダムに並び替えて、目的変数と相関係数を計算することを 10000回ほど行い、その相関係数の絶対値の最大値を rccmaxとした時に下式からpを算出する。


\displaystyle{
p = 100 \times ( 1 - r_{ccmax})
}

pythonでの実装

1290化合物の水溶解度のデータを使って、Borutaによる変数選択あり変数選択なしで予測精度に差が出るか検証する。
(変数選択なしで水溶解度を予測した記事はこちら

準備

データの取得

金子研究室HPよりデータを取得

モジュールのインストール

boruta、RDkit(記述子作成に使用)、DCEkit(金子研モジュール|前処理に使用)をインストールしておく

conda install -c conda-forge rdkit
pip install dcekit
pip install borutapy

解析諸条件

  • 記述子:RDkit記述子を使用
    (ただし、変数重要度が著しく高いMolLogPは除いた)
  • 前処理:分散が0の変数削除など基本的なもののみ。(詳細はコード参照)
  • テストデータの割合(0.5)

コード

モジュールインポート〜データの準備

import os
import pandas as pd
import numpy as np
import math

import matplotlib.figure as figure
import matplotlib.pyplot as plt

from rdkit import Chem
from rdkit.ML.Descriptors import MoleculeDescriptors
from rdkit.Chem import AllChem

from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestRegressor as rf
from sklearn.model_selection import train_test_split,GridSearchCV
from sklearn.metrics import r2_score, mean_absolute_error,mean_squared_error

from dcekit.variable_selection import search_high_rate_of_same_values, search_highly_correlated_variables
from boruta import BorutaPy

#データの読み込み
suppl = Chem.SDMolSupplier('logSdataset1290_2d.sdf')
mols = [mol for mol in suppl if mol is not None]

#RDkit記述子の計算
descriptor_names = [descriptor_name[0] for descriptor_name in Chem.Descriptors._descList]
descriptor_calculation = MoleculeDescriptors.MolecularDescriptorCalculator(descriptor_names)
RDkit = [descriptor_calculation.CalcDescriptors(mol) for mol in mols]
df = pd.DataFrame(RDkit, columns = descriptor_names)
df.drop('MolLogP',axis=1,inplace=True)

#溶解度(logS)の値を取得
df['logS'] =  [mol.GetProp('logS') if 'logS' in mol.GetPropNames() else 'None' for mol in mols]


#ここから変数選択
df = df.astype(float)

#欠損のある変数を削除
df.dropna(axis=1,inplace=True)

#分散が0の変数削除
del_num1 = np.where(df.var() == 0)
df = df.drop(df.columns[del_num1], axis=1)

#変数選択(互いに相関関係にある変数の内一方を削除)   
threshold_of_r = 0.95 #変数選択するときの相関係数の絶対値の閾値
corr_var = search_highly_correlated_variables(df, threshold_of_r)
df.drop(df.columns[corr_var], axis=1, inplace=True)

#同じ値を多くもつ変数の削除
inner_fold_number = 2 #CVでの分割数(予定)
rate_of_same_value = []
threshold_of_rate_of_same_value = (inner_fold_number-1)/inner_fold_number
for col in df.columns:
    same_value_number = df[col].value_counts()
    rate_of_same_value.append(float(same_value_number[same_value_number.index[0]] / df.shape[0]))
del_var_num = np.where(np.array(rate_of_same_value) >= threshold_of_rate_of_same_value)
df.drop(df.columns[del_var_num], axis=1, inplace=True)

print(df.shape)
# => (1290, 64)

64個の特徴量が残った。次にこれらの特徴量をさらにborutaで絞り込む。


データの分割〜borutaによる変数選択

#訓練・テストデータの分割と標準化
df_X = df.drop(['logS'],axis=1)
df_y = df.loc[:,['logS']]
X = df_X.values
y = df_y.values

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5, random_state=0)

#StandardScalerを使って標準化
scaler_X = StandardScaler()
scaler_y = StandardScaler()

scaled_X_train = scaler_X.fit_transform(X_train)
scaled_y_train = scaler_y.fit_transform(y_train)
scaled_X_test = scaler_X.transform(X_test)
scaled_y_test = scaler_y.transform(y_test)

#Borutaによる変数選択
#pパーセンタイルの最適化(金子研の方法)
corr_list = []
for n in range(10000):
    shadow_features = np.random.rand(scaled_X_train.shape[0]).T
    corr = np.corrcoef(scaled_X_train, shadow_features,rowvar=False)[-1]
    corr = abs(corr[corr < 0.95])
    corr_list.append(corr.max())
    corr_array = np.array(corr_list)
    perc = 100 * (1-corr_array.max())

#Borutaの実施
feat_selector = BorutaPy(rf(), 
                    n_estimators='auto',  # 特徴量の数に比例して、木の本数を増やす
                    verbose=0,
                    alpha=0.05, # 有意水準
                    max_iter=100, # 試行回数
                    perc = perc, #ランダム生成変数の重要度の何%を基準とするか
                    random_state=0
                    )

feat_selector.fit(scaled_X_train, scaled_y_train)

scaled_X_train_br = scaled_X_train[:,feat_selector.support_]
scaled_X_test_br = scaled_X_test[:,feat_selector.support_]

print('pパーセンタイル:',round(perc,2))
print('boruta後の変数の数:',scaled_X_train_br.shape[1])
"""output
pパーセンタイル:80.79
boruta後の変数の数: 43
"""

borutaによる変数選択で変数が64から47個まで絞り込まれた。


予測と精度検証

borutaによる変数選択ありの場合と、なしの場合両方でそれぞれ予測を行い、精度を比較する。

model = rf()
model_br = rf()
model.fit(scaled_X_train, scaled_y_train)
model_br.fit(scaled_X_train_br, scaled_y_train)

scaled_pred_y_br = np.ndarray.flatten(model_br.predict(scaled_X_test_br))
pred_y_br = scaler_y.inverse_transform(scaled_pred_y_br)
scaled_pred_y = np.ndarray.flatten(model.predict(scaled_X_test))
pred_y = scaler_y.inverse_transform(scaled_pred_y)

print('r2_br:{}'.format(round(r2_score(y_test, pred_y_test_br),3)))
print('RMSE_br:{}'.format(round(np.sqrt(mean_squared_error(y_test, pred_y_test_br)),3)))

print('r2:{}'.format(round(r2_score(y_test, pred_y_test),3)))
print('RMSE:{}'.format(round(np.sqrt(mean_squared_error(y_test, pred_y_test)),3)))

"""output
r2_br:0.883
RMSE_br:0.697
r2:0.883
RMSE:0.699
"""

結果まとめ

r2 RMSE
borutaで変数選択 0.883 0.697
変数選択なし 0.883 0. 699

今回のデータセットでは精度はほとんど上がらなかったが、精度を落とさずにより少ない変数で予測が行えているだけでも価値があるかもしれない。

参考

多重検定問題を回避した変数選択(Holm法)


前回は、多重検定問題を回避する方法の一つとして知られるBonferroni法(ボンフェローニ法)についてまとめたが、ついでに今回はHolm法(ホルム法)についてまとめた。

多重の多重性について

検定の多重性とは、仮設検定を何度も繰り返すことで、発生する問題のことをさす。

例えば、回帰や分類問題の解析において、無相関検定やt検定などを用いてp値が有意水準を下回った場合にその変数を説明変数と使用する場合を考える。
ある1つの変数を対象として検定をする場合には問題ないが、候補となる多くの変数に対して検定を実施する場合では、本来有意ではない(寄与や相関のない)変数のp値がたまたま有意水準を下回りることで有意であるみなされてしまう

この様に、複数回繰り返された検定全体において帰無仮説が棄却される可能性は、familywise error rateと呼ばれる。 例えば、有意水準α = 0.05の検定を20回繰り返すと、familywise error rate(1回でも帰無仮説棄却される可能性)は 0.642となる。

多重検定を回避する(Familly wise Erro rate)を調整する方法。

Familywise error rateを調整する方法は大きく二つに分けられる。

  1. F統計量やt統計量等の統計量に基づいた方法
  2. 統計量から算出されたp値のみを操作する方法

1 の方法としては、F統計量を用いた(Fisher's LSD)法、t統計量を用いたTukey's HSD法Dunnet法などが知られている。

2 の方法としては、ボンフェローニ(Bonferroni)法ホルム(Holm)法が提唱されている。これらの方法は、統計量に依存せず、どのような検定に対しても利用できる汎用性が高い方法である。

今回はHolm法を使って回帰分析における変数選択を検証する。

Holm法について

Holm法は非常に保守的なBonferroni法の有意水準の補正方法を、改善したもの。
以下の手順で、p値の大きさに従って有意水準\alphaを設定する。

  1. N個の帰無仮説を、p値の小さい順に並べる。
  2. 最もp値が小さい第1順位の帰無仮説有意水準\alpha/N に設定する。
    p値 < \alpha/N の場合:第1順位の帰無仮説を棄却(対立仮説を採択し)、3に進む。
    p値 > \alpha/N の場合:第1順位以下のすべての帰無仮説の判定を保留する。
  3. (第1順位の帰無仮説が棄却された場合)第2順位の帰無仮説有意水準\alpha/(N+1) に設定する。
    p値 < \alpha/(N+1) の場合:第2順位の帰無仮説を棄却(対立仮説を採択し)、4に進む。
    p値 > \alpha/(N+1) の場合:第2順位以下のすべての帰無仮説の判定を保留する。
  4. 上記を繰り返す。

つまり、p値が小さい順位検定を行い、棄却できた場合には、有意水準をどんどん小さく(\alpha をN, N+1, ...で割る)していき、帰無仮説が棄却されにくくしていく方法。

コード

Bostonの家賃データセットを使った家賃予測において、無相関検定で説明変数を選択する。
元の有意水準は0.05に設定し、変数選択したのちPLS回帰(部分的最小二乗回帰)で予測する。

準備

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.datasets import load_boston
from scipy.stats import pearsonr
from scipy.stats import spearmanr

# データ読み込み
boston = load_boston()
df_X = pd.DataFrame(boston.data, columns=boston.feature_names)
df_y = pd.DataFrame(boston.target,columns=['RENT'])

#サンプル数を絞る
df_X = df_X.iloc[::10,:]
df_y = df_y.iloc[::10,:]

print('変数の数:', df_X.shape[1])
#>>> 変数の数: 13

変数選択

目的変数のRENTとその他の変数で無相関検定を実施し、得られたp値がHolm法で補正した有意水準、および補正なしの有意水準を下回った場合にその変数を使用した。

#全ての変数と、目的変数「RENT」で無相関検定(ピアソン)を実施する。

p_values = []

for feature in df_X.columns:
    
    rent = df_y['RENT'].values
    feature = df_X[feature].values

    #ピアソンの相関係数とp値を計算する。
    result = pearsonr(rent, feature) 

    p_values.append(result[1])
    
p_values = np.array(p_values)

#Holm方による変数選択=======
#変数名とp値の辞書を作る
feature_dict = dict(zip(df_X.columns, p_values))


#p値の小さい順位に並び替える
feature_dict = dict(sorted(feature_dict.items(), key=lambda X:X[1]))

print('それぞれの変数のp値(昇順)')
for f, p in zip(feature_dict.keys(), feature_dict.values()):
    print('{}:{}'.format(f,round(p,5)))
"""output
それぞれの変数のp値(昇順)
LSTAT:0.0
RM:0.0
INDUS:0.00078
PTRATIO:0.00283
ZN:0.00391
NOX:0.00439
TAX:0.00542
AGE:0.00591
CRIM:0.01602
CHAS:0.02801
B:0.04196
RAD:0.05783
DIS:0.08387
"""

a = 0.05  
N = 1
selected_features_holm = []
for i, f in enumerate(feature_dict.keys()):
    
    a_holm = a/N
    
    #棄却できれば変数をリストに追加して、有意水準を下げる(Nに1をたす)
    if feature_dict[f] < a_holm:
        selected_features_holm.append(f)
        N += 1
    else:
        break
#Holm方による変数選択ここまで=======
   
#補正をせずに変数選択(比較用)
selected_features = df_X.columns[np.where(p_values < a)[0]]        

print('補正なしのαで選ばれた特徴量:',selected_features)
print('Holm法で選ばれた特徴量:',selected_features_holm)
"""output
補正なしのαで選ばれた特徴量: Index(['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'TAX', 'PTRATIO',
       'B', 'LSTAT'],
      dtype='object')
Holm法で選ばれた特徴量: ['LSTAT', 'RM', 'INDUS', 'PTRATIO', 'ZN', 'NOX', 'TAX', 'AGE']
"""

補正なしαでは11変数、Holm法使用時は8変数が選ばれた。

予測と精度検証

定石ではないが、学習/訓練データを分割せずに検定を実施したため、クロスバリデーションで精度を検証した。
plsの成分数は2で固定し、変数選択しない場合と、変数選択した場合補正ありなし)とでそれぞれ精度を比較した。

from sklearn.model_selection import cross_val_predict
from sklearn.cross_decomposition import PLSRegression
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import r2_score, mean_absolute_error,mean_squared_error

#選ばれた変数セットのデータを作成
df_X_sel = df_X.loc[:,selected_features]
df_X_holm = df_X.loc[:,selected_features_holm]

#標準化
scaler_X = StandardScaler()
scaler_X_sel = StandardScaler()
scaler_X_holm = StandardScaler()
scaler_y = StandardScaler()

scaled_X = scaler_X.fit_transform(df_X)
scaled_X_sel = scaler_X.fit_transform(df_X_sel)
scaled_X_holm = scaler_X.fit_transform(df_X_holm)
scaled_y = scaler_y.fit_transform(df_y)

#予測モデル作成
pls = PLSRegression(n_components=2)

#全変数使用で予測
scaled_pred_y = cross_val_predict(estimator=pls,
                                  X=scaled_X, 
                                  y=scaled_y,
                                  cv=4)
#通常の検定(補正なしα)で選択した変数
scaled_pred_y_sel = cross_val_predict(pls,
                                      scaled_X_sel, 
                                      scaled_y,
                                      cv=4)
#補正α使用して選択した変数
scaled_pred_y_holm = cross_val_predict(pls,
                                      scaled_X_holm, 
                                      scaled_y,
                                      cv=4)

#標準化された値を元に戻す。
pred_y = scaler_y.inverse_transform(scaled_pred_y)
pred_y_sel = scaler_y.inverse_transform(scaled_pred_y_sel)
pred_y_holm = scaler_y.inverse_transform(scaled_pred_y_holm)

#精度を計算する
r2 = r2_score(df_y.values, pred_y)
r2_sel = r2_score(df_y.values, pred_y_sel)
r2_holm = r2_score(df_y.values, pred_y_holm)
rmse = np.sqrt(mean_squared_error(df_y.values, pred_y))
rmse_sel = np.sqrt(mean_squared_error(df_y.values, pred_y_sel))
rmse_holm = np.sqrt(mean_squared_error(df_y.values, pred_y_holm))

print('r2(変数選択なし):',round(r2,3))
print('r2(補正なしα):',round(r2_sel,3))
print('r2(Holm法):',round(r2_holm,3))
print('RMSE(変数選択なし):',round(rmse,3))
print('RMSE(補正なしα):',round(rmse_sel,3))
print('RMSE(Holm法):',round(rmse_holm,3))
"""output
r2(変数選択なし): 0.417
r2(補正なしα): 0.402
r2(Holm法): 0.618
RMSE(変数選択なし): 6.685
RMSE(補正なしα): 6.768
RMSE(Holm法): 5.412
"""


精度評価指標として決定係数(r2)と二乗平均平方根誤差(RMSE)を算出した。 結果を表にまとめると以下のようになる。

r2 RMSE
変数選択なし 0.417 6.685
補正なしαで変数選択 0.402 6.768
Holm法で変数選択 0.618 5.412
Bonferroni法で変数選択 (前回) 0.678 4.971


学習・訓練データの分割をせずに全データを使って検定を実施という条件の下ではあるが、r2, RMSEともに変数選択した方が良好な結果となった。
ちなみに、前回記事で検証したBonferroni検定と比べると、今回の検討に関してはBonferroni法の方が良好な結果となった。

参考


多重検定問題を回避した変数選択(Bonferroni法)


多重検定問題を回避する方法の一つとして知られるBonferroni法(ボンフェローニ法)についてまとめた。

多重の多重性について

検定の多重性とは、仮説検定を何度も繰り返すことで、発生する問題のことをさす。

例えば、回帰や分類問題の解析において、無相関検定やt検定などを用いてp値が有意水準を下回った場合にその変数を説明変数と使用する場合を考える。
ある1つの変数を対象として検定をする場合には問題ないが、候補となる多くの変数に対して検定を実施する場合では、本来、有意ではない(寄与や相関のない)変数のp値がたまたま有意水準を下回り、有意なものとみなされてしまう。

この様に、複数回繰り返された検定全体において帰無仮説が棄却される可能性は、familywise error rateと呼ばれる。 例えば、有意水準α = 0.05の検定を20回繰り返すと、familywise error rate(1回でも帰無仮説棄却される可能性)は0.642となる。

多重検定問題を回避する(Famillywise Error rateを調整する)方法

Familywise error rateを調整する方法は大きく二つに分けられる。

  1. F統計量やt統計量等の統計量に基づいた方法
  2. 統計量から算出されたp値のみを操作する方法

1 の方法としては、F統計量を用いたFisher's LSD、t統計量を用いたTukey's HSD法Dunnet法などが知られている。

2 の方法としては、ボンフェローニ(Bonferroni)法ホルム(Holm)法が提唱されている。 これらの方法は統計量に依存せず、どのような検定に対しても利用できる汎用性が高い方法である。

今回はBonferroni法を使って回帰分析における変数選択を実施した。

Bonferroni法について

検定の全体を踏まえてp値を評価する方法。検定の回数をM回とすると、はじめに設定した有意水準 \alpha をMで割った \alpha' を新たな有意水準として採用する。


\displaystyle{
\alpha' = \frac{\alpha}{M}
}

Bonferroni法の問題点

検定の回数が多いと有意水準が極端に小さくなってしまい「本当は変数間の関係があるのに有意水準が小さすぎて検出できない」という問題がある)

このように第2種の過誤の可能性が高くなるため、「p値 > 補正α」となった場合は、帰無仮説は採択されるのではなく棄却が保留されると考えるのが妥当。

コード

Bostonの家賃データセットを使った家賃予測において、無相関検定で説明変数を選択する。
元の有意水準は0.05に設定し、変数選択したのちPLS回帰(部分的最小二乗回帰)で予測する。

準備

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.datasets import load_boston
from scipy.stats import pearsonr
from scipy.stats import spearmanr

# データ読み込み
boston = load_boston()
df_X = pd.DataFrame(boston.data, columns=boston.feature_names)
df_y = pd.DataFrame(boston.target,columns=['RENT'])

#サンプル数を絞る
df_X = df_X.iloc[::10,:]
df_y = df_y.iloc[::10,:]

print('変数の数:', df_X.shape[1])
#>>> 変数の数: 13

#Bonferroni法で補正したの有意水準と補正なしの有意水準(比較のため)を設定する
a = 0.05
a_bon = 0.05 / (df_X.shape[1])
print('有意水準:',a)
print('補正有意水準:',round(a,5))
"""output
有意水準: 0.05
補正有意水準: 0.00385
"""

変数選択

目的変数のRENTとその他の変数で無相関検定を実施し、得られたp値が Bonferroni法で補正した有意水準、および補正なしの有意水準を下回った場合にその変数を使用する。

#全ての変数と、目的変数「RENT」で無相関検定(ピアソン)を実施する。

p_values = []

for feature in df_X.columns:
    
    rent = df_y['RENT'].values
    feature = df_X[feature].values

    #ピアソンの相関係数とp値を計算する。
    result = pearsonr(rent, feature) 

    p_values.append(result[1])
    
p_values = np.array(p_values)

print('それぞれの変数のp値')
for f, p in zip(df_X.columns, p_values):
    print('{}:{}'.format(f,round(p,5)))
    
"""output
それぞれの変数のp値
CRIM:0.01602
ZN:0.00391
INDUS:0.00078
CHAS:0.02801
NOX:0.00439
RM:0.0
AGE:0.00591
DIS:0.08387
RAD:0.05783
TAX:0.00542
PTRATIO:0.00283
B:0.04196
LSTAT:0.0
"""

selected_features = df_X.columns[np.where(p_values < a)[0]]
selected_features_bon = df_X.columns[np.where(p_values < a_bon)[0]]

print('補正なしαで選ばれた特徴量:',selected_features)
print('補正αで選ばれた特徴量:',selected_features_bon)
"""output
補正なしαで選ばれた特徴量: Index(['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'TAX', 'PTRATIO', 'B', 'LSTAT'],
      dtype='object')
補正αで選ばれた特徴量: Index(['INDUS', 'RM', 'PTRATIO', 'LSTAT'], dtype='object')
"""

補正なしαでは11変数、補正ありαでは4変数が選ばれた。

予測と精度確認

定石ではないが、学習/訓練データを分割せずに検定を実施したため、クロスバリデーションで精度を検証した。 plsの成分数は2で固定し、変数選択しない場合と、変数選択した場合(補正ありとなし)とでそれぞれ精度を比較した。

from sklearn.model_selection import cross_val_predict
from sklearn.cross_decomposition import PLSRegression
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import r2_score, mean_absolute_error,mean_squared_error

#選ばれた変数セットのデータを作成
df_X_sel = df_X.loc[:,selected_features]
df_X_bon = df_X.loc[:,selected_features_bon]

#標準化
scaler_X = StandardScaler()
scaler_X_sel = StandardScaler()
scaler_X_bon = StandardScaler()
scaler_y = StandardScaler()

scaled_X = scaler_X.fit_transform(df_X)
scaled_X_sel = scaler_X.fit_transform(df_X_sel)
scaled_X_bon = scaler_X.fit_transform(df_X_bon)
scaled_y = scaler_y.fit_transform(df_y)

#予測モデル作成
pls = PLSRegression(n_components=2)

#全変数使用で予測
scaled_pred_y = cross_val_predict(estimator=pls,
                                  X=scaled_X, 
                                  y=scaled_y,
                                  cv=4)
#通常の検定(補正なしα)で選択した変数
scaled_pred_y_sel = cross_val_predict(pls,
                                      scaled_X_sel, 
                                      scaled_y,
                                      cv=4)
#補正α使用して選択した変数
scaled_pred_y_bon = cross_val_predict(pls,
                                      scaled_X_bon, 
                                      scaled_y,
                                      cv=4)

#標準化された値を元に戻す。
pred_y = scaler_y.inverse_transform(scaled_pred_y)
pred_y_sel = scaler_y.inverse_transform(scaled_pred_y_sel)
pred_y_bon = scaler_y.inverse_transform(scaled_pred_y_bon)

#精度を計算する
r2 = r2_score(df_y.values, pred_y)
r2_sel = r2_score(df_y.values, pred_y_sel)
r2_bon = r2_score(df_y.values, pred_y_bon)
rmse = np.sqrt(mean_squared_error(df_y.values, pred_y))
rmse_sel = np.sqrt(mean_squared_error(df_y.values, pred_y_sel))
rmse_bon = np.sqrt(mean_squared_error(df_y.values, pred_y_bon))

print('r2(変数選択なし):',round(r2,3))
print('r2(補正なしα):',round(r2_sel,3))
print('r2(補正ありα):',round(r2_bon,3))
print('RMSE(変数選択なし):',round(rmse,3))
print('RMSE(補正なしα):',round(rmse_sel,3))
print('RMSE(補正ありα):',round(rmse_bon,3))
"""output
r2(変数選択なし): 0.417
r2(補正なしα): 0.402
r2(補正ありα): 0.678
RMSE(変数選択なし): 6.685
RMSE(補正なしα): 6.768
RMSE(補正ありα): 4.971
"""


精度評価指標として決定係数(r2)と二乗平均平方根誤差(RMSE)を算出した。 結果を表にまとめると以下のようになる。


r2 RMSE
変数選択なし 0.417 6.685
補正なしαで変数選択 0.402 6.768
補正ありαで変数選択 0.678 4.971


学習・訓練データの分割をせずに全データを使って検定を実施という条件の下ではあるが、r2, RMSEともにBonferroni法で補正した有意水準αを用いて変数選択した方が良好な結果となった

参考


無相関検定(概要とpython実装)

f:id:chemoinfo-ottantacinque:20210131114724p:plain

ある因子(変数)ペアに互に相関があるか検定できる無相関検定についてまとめた。


無相関の検定について

  • 標本から得られた相関係数から「母集団にも同様の相関がある」と言えるかどうかを検定するもの。
  • 帰無仮説を「母相関係数は0である (同様の相関はない)」とおき、その棄却可否を調べる。

無相関検定での注意点

無相関検定のp値と相関係数rは独立に考える。 (無相関検定で評価するのは「相関が0でないかどうか」だけなので、相関が強いかどうかはp値からは主張できない)

相関係数r とp値のパターンとその解釈

  • rの絶対値が大きく、p値も有意水準を下回っている場合
    => 二つの変数の間になんららかの関係があると解釈するのが妥当
  • rの絶対値が小さいが、p値は有意水準を下回っている場合
    => 相関は存在すると考えられるが、関係の強さは弱いと解釈する。
  • p値が有意水準を上回っている場合
    => 標本から得られたデータからどのようなrの値が得られていても、実際に相関関係があるかどうかはこのデータからは結論付けられない。

pythonコード

scikit-learnのボストンの家賃データを利用する。
有意水準を0.05を定め、選択した二つの変数に相関があるか検定する。

変数の準備と相関のチェック

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.datasets import load_boston
from scipy.stats import pearsonr
from scipy.stats import spearmanr

# データ読み込み
boston = load_boston()
df_X = pd.DataFrame(boston.data, columns=boston.feature_names)
df_y = pd.DataFrame(boston.target,columns=['RENT'])
df = pd.concat([df_X,df_y],axis=1)

#全ての変数ペアの相関係数を求める。
df_pearson = df.corr(method='pearson')

#結果をヒートマップで可視化する
fig, ax = plt.subplots(figsize=(10, 8)) 
sns.heatmap(df_pearson, square=True, annot=True, fmt='.2f',annot_kws={'size':10},vmax=1, vmin=-1, center=0)
plt.title('Correlation coefficient (pearson)',fontsize=18)
plt.ylim(df_pearson.shape[1],0)
plt.savefig('Correlation_coefficient_pearson.png')
plt.show()
f:id:chemoinfo-ottantacinque:20210131113448p:plain

無相関検定を実施

相関係数の高いRENT とLSTATについて無相関検定を実施する。

#相関係数の高いRENT とLSTATについて無相関検定を実施する。
rent = df['RENT'].values
lstat = df['LSTAT'].values

#ピアソンの相関係数とp値を計算する。
result = pearsonr(rent, lstat) 

#スピアマンの相関係数とp値を計算する場合
#result = pearsonr(rent, lstat) 

r_value = result[0]
p_value = result[1]
print('相関係数:', r_value)
print('p値:', p_value)
"""output
相関係数: -0.7376627261740148
p値: 5.08110339438697e-88
"""

p値<有意水準(0.05)であるため帰無仮説は棄却され、二つの変数は無相関ではないだろうと考えられる。

極端にサンプル数を減らしてみる

#サンプル抽出する
rent = df['RENT'].values[:5]
lstat = df['LSTAT'].values[:5]

result = pearsonr(rent, lstat) 

r_value = result[0]
p_value = result[1]
print('相関係数:', r_value)
print('p値:', p_value)
"""output
相関係数: -0.7036716841380781
p値: 0.18479027554023453
"""

極端にサンプル数が少ない場合には、p値>有意水準(0.05)となり、「母相関係数は0である (同様の相関はない)」という仮説は棄却されない
つまり少数サンプルでは何も結論づけられないとうこと。

詳細(計算方法)

  1. 仮説を設定する。
    帰無仮説H0:母相関係数は0である (同様の相関はない)
    対立仮説H1:母相関係数は0ではない
  2. 有意水準\alphaを定める。
  3. 下記の式から、統計量tを求める。
    (統計量tは自由度\nu = n-2に従う)
  4. p値を求め棄却有無を判定する。

\displaystyle{
t = \frac{|r|\sqrt{n-2}}{\sqrt{1-r^2}} 
}


import numpy as np
from scipy import stats

#相関係数の高いRENT とLSTATについて無相関検定を実施する。
rent = df['RENT'].values[:8]
lstat = df['LSTAT'].values[:8]

r = np.corrcoef(rent, lstat)[0][1]
n = len(rent)
t = abs(r) * np.sqrt(n-2) / np.sqrt(1-r**2)
nu = n-2
print('相関係数r:',r)
print('サンプル数n:',n)
print('統計量t:',t)
print('自由度ν:',nu)
"""output
相関係数r: -0.48538333300088693
サンプル数n: 8
統計量t: 1.3598759326365544
自由度ν: 6
"""

# x軸の等差数列を生成
X = np.linspace(start=-6, stop=6, num=100)

# t分布
t_sf = stats.t.pdf(x=X, df=nu)

#累積確率関数(cdf)と生存関数(sf)でp値を求める。
one_sided_pval1 = stats.t.cdf(t, nu)
one_sided_pval2 = stats.t.sf(t, nu)

# 可視化
plt.plot(X, t_sf)

X_lower = np.linspace(start=-6, stop=t, num=100)
X_upper = np.linspace(start=t, stop=6, num=100)
t_lower = stats.t.pdf(x=X_lower, df=nu)
t_upper = stats.t.pdf(x=X_upper, df=nu)

base = np.zeros(100)

plt.fill_between(X_lower, base, t_lower, facecolor='lime', alpha=0.5)
plt.fill_between(X_upper, base, t_upper, facecolor='yellow', alpha=0.5)

#求めたt値をプロット
plt.plot(t, 0.00009, 's',label='t value') 

plt.title('t distribution({})'.format(nu), fontsize=15)
plt.ylim(0,0.005)
plt.xlim(3,6)
plt.ylim(0,0.5)
plt.xlim(-6,6)
plt.legend(fontsize=15)
plt.savefig('無相関検定.png')
plt.show()

print('p-value(lower side):',round(one_sided_pval1,4))
print('p-value(upper side):',round(one_sided_pval2,4))
"""
p-value(lower side): 0.8886
p-value(upper side): 0.1114
"""
f:id:chemoinfo-ottantacinque:20210131113856p:plain

算出したp値(黄緑と黄色の面積に該当)はどちらも有意水準(0.05)より大きいので、,「母相関係数は0である (同様の相関はない)」という仮説は棄却されない

参考