データが少ない医療現場でも大丈夫？ランダムフォレストと回帰分析の違いをMSEとSHAP値で実験比較

import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.ensemble import RandomForestRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_absolute_error
from sklearn.metrics import mean_squared_error
import shap
import matplotlib.pyplot as plt

# データ生成
np.random.seed(42)
n = 2000
fim_in = np.random.randint(20, 90, n)
age = np.random.randint(40, 95, n)
mmse = np.random.randint(5, 30, n)
paralysis = np.random.randint(0, 5, n)

noise = np.random.normal(0, 5, n)
fim_out = (
    fim_in * 3 +
    30 * np.sin(mmse / 10) -
    (age - 35) / 2 +
    (30 - paralysis**2) / 4 +
    mmse * (6 - paralysis) +
    noise
)

fim_out = fim_out / (fim_out.max() / 126)

df = pd.DataFrame({
    'fim_in': fim_in,
    'age': age,
    'mmse': mmse,
    'paralysis': paralysis,
    'fim_out': fim_out
})

# 特徴量と目的変数
features = ["fim_in", "age", "mmse", "paralysis"]
target = "fim_out"

sample_sizes = [50, 100, 150,200,250,300,350,400,450,500,550, 600,650,700,750,800]
num_iterations = 10

results_lr = []
results_rf = []
results_std_lr = []
results_std_rf = []

for size in sample_sizes:
    mae_list_lr = []
    mae_list_rf = []

    for i in range(num_iterations):
        sample_df = df.sample(n=size, random_state=i)
        X = sample_df[features]
        y = sample_df[target]

        X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=i)

        # 線形回帰モデル
        scaler = StandardScaler()
        X_train_scaled = scaler.fit_transform(X_train)
        X_test_scaled = scaler.transform(X_test)

        model_lr = LinearRegression()
        model_lr.fit(X_train_scaled, y_train)
        y_pred_lr = model_lr.predict(X_test_scaled)
        mae_lr = mean_squared_error(y_test, y_pred_lr)
        mae_list_lr.append(mae_lr)

        # ランダムフォレストモデル
        model_rf = RandomForestRegressor(n_estimators=100, random_state=i)
        model_rf.fit(X_train, y_train)
        y_pred_rf = model_rf.predict(X_test)
        mae_rf = mean_squared_error(y_test, y_pred_rf)
        mae_list_rf.append(mae_rf)

        # SHAP解析
        if size == max(sample_sizes) and i == 0:
            print("SHAP解析を実行中...")

            explainer_rf = shap.Explainer(model_rf, X_train)
            shap_values_rf = explainer_rf(X_test)
            shap.summary_plot(shap_values_rf, X_test, plot_type="bar")
            shap.summary_plot(shap_values_rf, X_test)
            plt.title("Random Forest SHAP Feature Importance")
            plt.show()

            explainer_lr = shap.Explainer(model_lr, X_train_scaled)
            shap_values_lr = explainer_lr(X_test_scaled)
            shap.summary_plot(shap_values_lr, X_test_scaled, feature_names=features, plot_type="bar")
            shap.summary_plot(shap_values_lr, X_test_scaled, feature_names=features)
            plt.title("linear regression SHAP Feature Importance")
            plt.show()

    results_lr.append(np.mean(mae_list_lr))
    results_rf.append(np.mean(mae_list_rf))
    results_std_lr.append(np.std(mae_list_lr))
    results_std_rf.append(np.std(mae_list_rf))

# 結果表示
print("サンプルサイズ:", sample_sizes)
print("線形回帰の平均MAE:", results_lr)
print("ランダムフォレストの平均MAE:", results_rf)

# MAE推移の可視化（平均のみ）
plt.plot(sample_sizes, results_lr, label="Linear Regression", marker="o")
plt.plot(sample_sizes, results_rf, label="Random Forest", marker="s")
plt.xlabel("Sample Size")
plt.ylabel("Mean squer Error")
plt.title("モデル性能の比較")
plt.legend()
plt.grid(True)
plt.show()

# MAEと標準偏差の推移（エラーバー付き）
plt.figure(figsize=(10, 6))

results_lr = np.array(results_lr)
results_rf = np.array(results_rf)
results_std_lr = np.array(results_std_lr)
results_std_rf = np.array(results_std_rf)

plt.plot(sample_sizes, results_rf, label="Random Forest", color="blue")
plt.fill_between(sample_sizes, results_rf - results_std_rf, results_rf + results_std_rf, color="blue", alpha=0.2)

plt.plot(sample_sizes, results_lr, label="Linear Regression", color="green")
plt.fill_between(sample_sizes, results_lr - results_std_lr, results_lr + results_std_lr, color="green", alpha=0.2)

plt.xlabel("Sample Size")
plt.ylabel("Mean squared Error (MAE)")
plt.title("Comparison of MSE vs. Sample Size\n(Random Forest vs. Linear Regression)")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()