【python】尤度比検定で統計モデルの比較をしよう|統計的仮説検定

import pandas as pd
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn import metrics
import numpy as np
from scipy.stats import chi2
from sklearn.ensemble import RandomForestClassifier
import matplotlib.pyplot as plt

# データセットのロード
data_breast_cancer = load_breast_cancer()
# Pandasによるデータの表示
df_target = pd.DataFrame(data_breast_cancer["target"], columns=["target"])
df_data = pd.DataFrame(data_breast_cancer["data"], columns=data_breast_cancer["feature_names"])
df = pd.concat([df_target, df_data], axis=1)

y = df["target"]
X = df.loc[:, "mean radius":]

# 教師データとテストデータに分割する
X_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.2, random_state=42)
seed=42

#グリッドサーチでハイパーパラメーターのチューニングを行う。
scores = {}
for max_depth in range(3, 10, 1):
    for min_samples_leaf in range(2, 10, 1):
        model = RandomForestClassifier(max_depth=max_depth, 
                          min_samples_leaf=min_samples_leaf, 
                          n_estimators=100, 
                          random_state=seed
                          )
        model.fit(X_train, y_train)
        scores[(max_depth, min_samples_leaf)] = model.score(X_valid, y_valid)

# 検証結果をscoresに格納 
scores = pd.Series(scores)

# 表示
print('ベストスコア:{:.2f}'.format(scores.max()))
print('その時のパラメータ(max_depth, min_samples_leaf):{}'.format(scores.idxmax()))

array_bst = np.array(scores.idxmax())

max_depth = array_bst[0]
min_samples_leaf = array_bst[1]
model_bst = RandomForestClassifier(max_depth=max_depth,
                          min_samples_leaf=min_samples_leaf,
                          n_estimators=100, 
                          random_state=seed
                          )

model_bst.fit(X_train, y_train)

# 過学習確認
pred_rfc_class = model_bst.predict(X_valid)
pred_rfc_proba = model_bst.predict_proba(X_valid)[:, 1]

# 過学習確認
print("教師データ",model_bst.score(X_train, y_train))
print("テストスコア",model_bst.score(X_valid, y_valid))



#feature_importanceで特徴量重要度を算出する
feature = pd.DataFrame()
feature["columns"] = X.columns
feature["score"] = model_bst.feature_importances_
feature = feature.sort_values(by='score', ascending=False)

plt.figure(figsize=(50,40))
plt.plot(feature["columns"], feature["score"])
plt.title("feature_importance")
plt.show()

# データセットの分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# ロジスティック回帰による学習
clf_full = RandomForestClassifier()
clf_full.fit(X_train, y_train)

# 特徴量を削減していないモデルで尤度を算出
y_prob_full = clf_full.predict_proba(X_test)

# 特徴量を削減したモデルをロジスティック回帰で作る
X_train_reduced = X_train[["worst area","worst concave points","mean concave points","worst radius",
                         "worst perimeter","mean concavity","mean area","mean perimeter","mean radius","worst concavity"]]
X_test_reduced = X_test[["worst area","worst concave points","mean concave points","worst radius",
                         "worst perimeter","mean concavity","mean area","mean perimeter","mean radius","worst concavity"]]
clf_reduced = RandomForestClassifier()
clf_reduced.fit(X_train_reduced, y_train)
y_prob_reduced = clf_reduced.predict_proba(X_test_reduced)

# 対数尤度の計算
ll_full = metrics.log_loss(y_test, y_prob_full, labels=[0,1])
ll_reduced = metrics.log_loss(y_test, y_prob_reduced, labels=[0,1])

# 尤度比検定統計量の算出
lrt = 2 * (ll_reduced - ll_full)

# p値の算出
p_value = 1 - chi2.cdf(lrt, 1)

print("Likelihood Ratio Test statistic: ", lrt)
print("P-value: ", p_value)