{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.14","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":35332,"databundleVersionId":3723648,"sourceType":"competition"},{"sourceId":3739819,"sourceType":"datasetVersion","datasetId":2231132}],"dockerImageVersionId":30804,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"pip install powershap -qq","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-12T07:37:53.384468Z","iopub.execute_input":"2025-06-12T07:37:53.384853Z","iopub.status.idle":"2025-06-12T07:38:13.372433Z","shell.execute_reply.started":"2025-06-12T07:37:53.384806Z","shell.execute_reply":"2025-06-12T07:38:13.370993Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom sklearn.preprocessing import StandardScaler\nfrom scipy.stats import skew, kurtosis\nimport contextlib\nimport io\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.mplot3d import Axes3D\nfrom matplotlib.colors import TwoSlopeNorm\nimport seaborn as sns\nfrom sklearn.decomposition import PCA\nfrom sklearn.model_selection import train_test_split, KFold\nfrom sklearn.preprocessing import PolynomialFeatures\nfrom xgboost import XGBClassifier\nfrom sklearn.metrics import roc_auc_score, roc_curve, auc\nfrom sklearn.linear_model import LogisticRegression\nimport itertools\nimport shap\nfrom powershap import PowerShap\nfrom sklearn.ensemble import RandomForestClassifier\nfrom sklearn.neural_network import MLPClassifier\nfrom sklearn.inspection import PartialDependenceDisplay\nimport warnings\n\n\nwarnings.filterwarnings(\"ignore\")","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-06-12T07:38:13.375089Z","iopub.execute_input":"2025-06-12T07:38:13.375995Z","iopub.status.idle":"2025-06-12T07:38:17.374426Z","shell.execute_reply.started":"2025-06-12T07:38:13.375939Z","shell.execute_reply":"2025-06-12T07:38:17.373337Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"data = pd.read_parquet('../input/amex-data-integer-dtypes-parquet-format/train.parquet')\ntarget = pd.read_csv('../input/amex-default-prediction/train_labels.csv')\ndata = data.groupby(\"customer_ID\").last()\ndata = data.merge(target, on = 'customer_ID', how = 'left')\ndel target","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-12T07:38:17.375868Z","iopub.execute_input":"2025-06-12T07:38:17.376550Z","iopub.status.idle":"2025-06-12T07:38:47.941098Z","shell.execute_reply.started":"2025-06-12T07:38:17.376490Z","shell.execute_reply":"2025-06-12T07:38:47.939749Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"data","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-12T07:38:47.943540Z","iopub.execute_input":"2025-06-12T07:38:47.943917Z","iopub.status.idle":"2025-06-12T07:38:48.142406Z","shell.execute_reply.started":"2025-06-12T07:38:47.943880Z","shell.execute_reply":"2025-06-12T07:38:48.140931Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# preprocess","metadata":{}},{"cell_type":"code","source":"# 獲取X和y\ndef get_X_and_y(df): return df.drop(columns=['target', 'customer_ID', 'S_2']), df['target']\n\n# 計算違義率\ndef caculate_default_rate(df): print(f\"Default Rate: {df.mean():.4%}\")\n\n# 拆分訓練集與測試集\ndef separate_train_and_test(X, y): return train_test_split(X, y, test_size=0.2, random_state=42)\n\n# one-hot編碼\ndef one_hot_encode(df):\n    for feature in df.select_dtypes(include=[\"object\"]):\n        df[feature] = df[feature].fillna(\"miss\")\n        for value in df[feature].unique(): df[f'{feature}_{value}'] = (df[feature] == value).astype(int)\n        df.drop(columns=[feature], inplace=True)\n    return df\n\n# 繪製PCA\ndef plot_pca(X_train, X_test, y, filename):\n    pca = PCA(n_components=2).fit_transform(pd.concat([X_train, X_test], axis=0).fillna(pd.concat([X_train, X_test], axis=0).mean()).reset_index(drop=True))\n    df_pca = pd.DataFrame(pca, columns=['PC1', 'PC2']).assign(TARGET=y.loc[X_train.index.union(X_test.index)].reset_index(drop=True))\n    plt.figure(figsize=(10, 6))\n    for t, c in zip([0, 1], ['blue', 'red']):\n        plt.scatter(df_pca.loc[df_pca['TARGET'] == t, 'PC1'], df_pca.loc[df_pca['TARGET'] == t, 'PC2'], color=c, alpha=0.5, label=f'TARGET={t}')\n    plt.xlabel('PC1'), plt.ylabel('PC2'), plt.grid(), plt.legend(title='TARGET', loc='upper right')\n    plt.savefig(filename, transparent=True, dpi=300), plt.show()\n\n# log轉換\ndef log_transformate(df):\n    for feature in df.select_dtypes(include=['number']).columns[~df.isin([0, 1]).all()]:\n        skewed = skew(df[feature].dropna())\n        if skewed > 1:\n            min_val = df[feature].min(skipna=True)\n            if min_val > 1:\n                df[feature] = np.log(df[feature]).where(df[feature] > 0, np.nan)\n            elif 0 < min_val <= 1:\n                df[feature] = np.log(df[feature] + 1).where(df[feature] + 1 > 0, np.nan)\n            else:\n                df[feature] = np.log(df[feature] + abs(min_val) + 1).where(df[feature] + abs(min_val) + 1 > 0, np.nan)\n    return df\n\ndef standardscaler(df):\n    for feature in df.select_dtypes(include=['number']).columns[~df.isin([0, 1]).all()]:\n        mean = df[feature].mean(skipna=True)\n        std = df[feature].std(skipna=True)\n        if pd.notna(std) and std > 1e-8:\n            df[feature] = (df[feature] - mean) / std\n        else:\n            df[feature] = 0\n    return df\n\n# 填充缺失值\ndef missing_value_impute(df):\n    for feature in df.select_dtypes(include=['number']).columns[~df.isin([0, 1]).all()]:\n        df[f\"{feature}_missing\"] = df[feature].isnull().astype(int)\n        df[feature].fillna(0, inplace=True)\n    return df\n\n# 對齊特徵數\ndef align_features(X_train, X_test):\n    for col in set(X_train.columns) - set(X_test.columns): X_test[col] = 0\n    for col in set(X_test.columns) - set(X_train.columns): X_train[col] = 0\n    return X_train[sorted(X_train.columns)], X_test[sorted(X_test.columns)]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-12T07:38:48.144221Z","iopub.execute_input":"2025-06-12T07:38:48.144550Z","iopub.status.idle":"2025-06-12T07:38:48.163493Z","shell.execute_reply.started":"2025-06-12T07:38:48.144514Z","shell.execute_reply":"2025-06-12T07:38:48.162170Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"X, y = get_X_and_y(data)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-12T07:38:48.165583Z","iopub.execute_input":"2025-06-12T07:38:48.165969Z","iopub.status.idle":"2025-06-12T07:38:48.246016Z","shell.execute_reply.started":"2025-06-12T07:38:48.165931Z","shell.execute_reply":"2025-06-12T07:38:48.244597Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"caculate_default_rate(y)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-12T07:38:48.247598Z","iopub.execute_input":"2025-06-12T07:38:48.247976Z","iopub.status.idle":"2025-06-12T07:38:48.255745Z","shell.execute_reply.started":"2025-06-12T07:38:48.247921Z","shell.execute_reply":"2025-06-12T07:38:48.253903Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"X_train, X_test, y_train, y_test = separate_train_and_test(X, y)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-12T07:38:48.257623Z","iopub.execute_input":"2025-06-12T07:38:48.258023Z","iopub.status.idle":"2025-06-12T07:38:48.663296Z","shell.execute_reply.started":"2025-06-12T07:38:48.257983Z","shell.execute_reply":"2025-06-12T07:38:48.662210Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"X_train = one_hot_encode(X_train)\nX_test = one_hot_encode(X_test)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-12T07:38:48.664515Z","iopub.execute_input":"2025-06-12T07:38:48.664846Z","iopub.status.idle":"2025-06-12T07:38:48.670059Z","shell.execute_reply.started":"2025-06-12T07:38:48.664812Z","shell.execute_reply":"2025-06-12T07:38:48.668872Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plot_pca(X_train, X_test, y, 'original_pca')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-12T07:38:48.674182Z","iopub.execute_input":"2025-06-12T07:38:48.674960Z","iopub.status.idle":"2025-06-12T07:39:03.864503Z","shell.execute_reply.started":"2025-06-12T07:38:48.674910Z","shell.execute_reply":"2025-06-12T07:39:03.859288Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"X_train = log_transformate(X_train)\nX_test = log_transformate(X_test)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-12T07:39:03.870005Z","iopub.execute_input":"2025-06-12T07:39:03.872093Z","iopub.status.idle":"2025-06-12T07:39:12.238909Z","shell.execute_reply.started":"2025-06-12T07:39:03.871843Z","shell.execute_reply":"2025-06-12T07:39:12.237184Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"X_train = standardscaler(X_train)\nX_test = standardscaler(X_test)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-12T07:39:12.240840Z","iopub.execute_input":"2025-06-12T07:39:12.241302Z","iopub.status.idle":"2025-06-12T07:39:23.557772Z","shell.execute_reply.started":"2025-06-12T07:39:12.241256Z","shell.execute_reply":"2025-06-12T07:39:23.556506Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"X_train = missing_value_impute(X_train)\nX_test = missing_value_impute(X_test)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-12T07:39:23.559191Z","iopub.execute_input":"2025-06-12T07:39:23.559542Z","iopub.status.idle":"2025-06-12T07:39:32.332330Z","shell.execute_reply.started":"2025-06-12T07:39:23.559507Z","shell.execute_reply":"2025-06-12T07:39:32.330586Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"X_train, X_test = align_features(X_train, X_test)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-12T07:39:32.333924Z","iopub.execute_input":"2025-06-12T07:39:32.334389Z","iopub.status.idle":"2025-06-12T07:39:35.217778Z","shell.execute_reply.started":"2025-06-12T07:39:32.334341Z","shell.execute_reply":"2025-06-12T07:39:35.216401Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plot_pca(X_train, X_test, y, 'preprocessed_pca')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-12T07:39:35.219365Z","iopub.execute_input":"2025-06-12T07:39:35.219831Z","iopub.status.idle":"2025-06-12T07:39:59.520111Z","shell.execute_reply.started":"2025-06-12T07:39:35.219787Z","shell.execute_reply":"2025-06-12T07:39:59.518582Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# preprocessed auc","metadata":{}},{"cell_type":"code","source":"# 比較不同前處理步驟的LR結果\ndef compare_preprocessing_strategies(X, y):\n    def apply_pipeline(X, steps):\n        X = X.copy()\n        for step in steps:\n            X = step(X)\n        return X\n\n    strategies = [\n        (\"one-hot + miss\",               [one_hot_encode, missing_value_impute]),\n        (\"one-hot + log + miss\",     [one_hot_encode, log_transformate, missing_value_impute]),\n        (\"one-hot + Z-std + miss\",     [one_hot_encode, standardscaler, missing_value_impute]),\n        (\"one-hot + log + Z-std + miss\", [one_hot_encode, log_transformate, standardscaler, missing_value_impute]),\n        (\"one-hot + Z-std + log + miss\", [one_hot_encode, standardscaler, log_transformate, missing_value_impute]),\n    ]\n\n    results = []\n\n    for name, pipeline in strategies:\n        # 切分資料\n        X_train_raw, X_test_raw, y_train, y_test = separate_train_and_test(X, y)\n\n        # 分別應用前處理流程\n        X_train = apply_pipeline(X_train_raw, pipeline)\n        X_test = apply_pipeline(X_test_raw, pipeline)\n\n        # 對齊特徵\n        X_train, X_test = align_features(X_train, X_test)\n\n        # 建立模型並訓練\n        model = LogisticRegression(max_iter=1000)\n        model.fit(X_train, y_train)\n\n        # 計算 AUC\n        y_train_pred = model.predict_proba(X_train)[:, 1]\n        y_test_pred = model.predict_proba(X_test)[:, 1]\n        train_auc = roc_auc_score(y_train, y_train_pred)\n        test_auc = roc_auc_score(y_test, y_test_pred)\n\n        results.append({\n            \"preprocess\": name,\n            \"train auc\": train_auc,\n            \"test auc\": test_auc\n        })\n\n    return pd.DataFrame(results)\n\ndef plot_preprocessed_roc_and_auc_df(X_train, y_train, X_test, y_test):\n    models = {\n        \"xgb\": lambda: XGBClassifier(random_state=42),\n        \"lr\": lambda: LogisticRegression(max_iter=1000),\n        \"rf\": lambda: RandomForestClassifier(n_estimators=100, random_state=42),\n        \"nn\": lambda: MLPClassifier(hidden_layer_sizes=(32,), activation='relu', solver='adam', max_iter=500, random_state=42)\n    }\n\n    kf = KFold(n_splits=10, shuffle=True, random_state=42)\n    mean_fpr = np.linspace(0, 1, 100)\n\n    # 改為存儲四種 AUC\n    auc_scores = {\n        \"cv_train\": {},\n        \"cv_validation\": {},\n        \"full_train\": {},\n        \"test\": {}\n    }\n\n    for model_name, model_fn in models.items():\n        train_aucs = []\n        val_aucs = []\n        plt.figure(figsize=(8, 6))\n\n        for train_idx, val_idx in kf.split(X_train):\n            X_kf_train, X_kf_val = X_train.iloc[train_idx], X_train.iloc[val_idx]\n            y_kf_train, y_kf_val = y_train.iloc[train_idx], y_train.iloc[val_idx]\n            model = model_fn()\n            model.fit(X_kf_train, y_kf_train)\n\n            # 計算 train AUC for current fold\n            y_train_pred = model.predict_proba(X_kf_train)[:, 1]\n            train_auc = roc_auc_score(y_kf_train, y_train_pred)\n            train_aucs.append(train_auc)\n\n            # 計算 validation AUC for current fold\n            y_val_pred = model.predict_proba(X_kf_val)[:, 1]\n            fpr, tpr, _ = roc_curve(y_kf_val, y_val_pred)\n            val_auc = auc(fpr, tpr)\n            val_aucs.append(val_auc)\n\n            tpr_interp = np.interp(mean_fpr, fpr, tpr)\n            tpr_interp[0] = 0.0\n            plt.plot(mean_fpr, tpr_interp, alpha=0.3, label=f'{model_name} Fold {len(val_aucs)}')\n\n        # 儲存交叉驗證結果\n        auc_scores[\"cv_train\"][model_name] = np.mean(train_aucs)\n        auc_scores[\"cv_validation\"][model_name] = np.mean(val_aucs)\n\n        # 繪製平均 ROC 曲線\n        tprs = []\n        for train_idx, val_idx in kf.split(X_train):\n            model = model_fn()\n            model.fit(X_train.iloc[train_idx], y_train.iloc[train_idx])\n            y_val_score = model.predict_proba(X_train.iloc[val_idx])[:, 1]\n            fpr, tpr, _ = roc_curve(y_train.iloc[val_idx], y_val_score)\n            tpr_interp = np.interp(mean_fpr, fpr, tpr)\n            tpr_interp[0] = 0.0\n            tprs.append(tpr_interp)\n        \n        mean_tpr = np.mean(tprs, axis=0)\n        mean_tpr[-1] = 1.0\n        \n        mean_auc = auc(mean_fpr, mean_tpr)\n        plt.plot(mean_fpr, mean_tpr, label=f'Mean ROC (AUC = {mean_auc:.4f})')\n\n        # 整體 train fit（不分折）\n        model = model_fn()\n        model.fit(X_train, y_train)\n        y_full_train_pred = model.predict_proba(X_train)[:, 1]\n        full_train_auc = roc_auc_score(y_train, y_full_train_pred)\n        auc_scores[\"full_train\"][model_name] = full_train_auc\n\n        # Test AUC\n        y_test_pred = model.predict_proba(X_test)[:, 1]\n        test_auc = roc_auc_score(y_test, y_test_pred)\n        auc_scores[\"test\"][model_name] = test_auc\n\n        plt.plot([0, 1], [0, 1], color='gray', linestyle='--')\n        plt.xlabel('False Positive Rate')\n        plt.ylabel('True Positive Rate')\n        plt.legend(loc='lower right')\n        plt.savefig(f'{model_name}_roc_curve.png', transparent=True, dpi=300)\n        plt.show()\n\n    # 結果轉成 DataFrame\n    proprocessed_auc_df = pd.DataFrame(auc_scores).reset_index()\n    proprocessed_auc_df.rename(columns={\"index\": \"model\"}, inplace=True)\n    return proprocessed_auc_df\n\ndef plot_preprocessed_roc_and_auc_df(X_train, y_train, X_test, y_test):\n    models = {\n        \"xgb\": lambda: XGBClassifier(random_state=42),\n        \"lr\": lambda: LogisticRegression(max_iter=1000),\n        \"rf\": lambda: RandomForestClassifier(n_estimators=100, random_state=42),\n        \"nn\": lambda: MLPClassifier(hidden_layer_sizes=(32,), activation='relu', solver='adam', max_iter=500, random_state=42)\n    }\n\n    kf = KFold(n_splits=10, shuffle=True, random_state=42)\n    mean_fpr = np.linspace(0, 1, 100)\n\n    # 改為存儲四種 AUC\n    auc_scores = {\n        \"cv_train\": {},\n        \"cv_validation\": {},\n        \"full_train\": {},\n        \"test\": {}\n    }\n\n    for model_name, model_fn in models.items():\n        train_aucs = []\n        val_aucs = []\n        plt.figure(figsize=(8, 6))\n\n        for train_idx, val_idx in kf.split(X_train):\n            X_kf_train, X_kf_val = X_train.iloc[train_idx], X_train.iloc[val_idx]\n            y_kf_train, y_kf_val = y_train.iloc[train_idx], y_train.iloc[val_idx]\n            model = model_fn()\n            model.fit(X_kf_train, y_kf_train)\n\n            # 計算 train AUC for current fold\n            y_train_pred = model.predict_proba(X_kf_train)[:, 1]\n            train_auc = roc_auc_score(y_kf_train, y_train_pred)\n            train_aucs.append(train_auc)\n\n            # 計算 validation AUC for current fold\n            y_val_pred = model.predict_proba(X_kf_val)[:, 1]\n            fpr, tpr, _ = roc_curve(y_kf_val, y_val_pred)\n            val_auc = auc(fpr, tpr)\n            val_aucs.append(val_auc)\n\n            tpr_interp = np.interp(mean_fpr, fpr, tpr)\n            tpr_interp[0] = 0.0\n            plt.plot(mean_fpr, tpr_interp, alpha=0.3, label=f'{model_name} Fold {len(val_aucs)}')\n\n        # 儲存交叉驗證結果\n        auc_scores[\"cv_train\"][model_name] = np.mean(train_aucs)\n        auc_scores[\"cv_validation\"][model_name] = np.mean(val_aucs)\n\n        # 繪製平均 ROC 曲線\n        tprs = []\n        for train_idx, val_idx in kf.split(X_train):\n            model = model_fn()\n            model.fit(X_train.iloc[train_idx], y_train.iloc[train_idx])\n            y_val_score = model.predict_proba(X_train.iloc[val_idx])[:, 1]\n            fpr, tpr, _ = roc_curve(y_train.iloc[val_idx], y_val_score)\n            tpr_interp = np.interp(mean_fpr, fpr, tpr)\n            tpr_interp[0] = 0.0\n            tprs.append(tpr_interp)\n        \n        mean_tpr = np.mean(tprs, axis=0)\n        mean_tpr[-1] = 1.0\n        \n        mean_auc = auc(mean_fpr, mean_tpr)\n        plt.plot(mean_fpr, mean_tpr, label=f'Mean ROC (AUC = {mean_auc:.4f})')\n\n        # 整體 train fit（不分折）\n        model = model_fn()\n        model.fit(X_train, y_train)\n        y_full_train_pred = model.predict_proba(X_train)[:, 1]\n        full_train_auc = roc_auc_score(y_train, y_full_train_pred)\n        auc_scores[\"full_train\"][model_name] = full_train_auc\n\n        # Test AUC\n        y_test_pred = model.predict_proba(X_test)[:, 1]\n        test_auc = roc_auc_score(y_test, y_test_pred)\n        auc_scores[\"test\"][model_name] = test_auc\n\n        plt.plot([0, 1], [0, 1], color='gray', linestyle='--')\n        plt.xlabel('False Positive Rate')\n        plt.ylabel('True Positive Rate')\n        plt.legend(loc='lower right')\n        plt.savefig(f'{model_name}_roc_curve.png', transparent=True, dpi=300)\n        plt.show()\n\n    # 結果轉成 DataFrame\n    proprocessed_auc_df = pd.DataFrame(auc_scores).reset_index()\n    proprocessed_auc_df.rename(columns={\"index\": \"model\"}, inplace=True)\n    return proprocessed_auc_df","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# compare_preprocessing_df = compare_preprocessing_strategies(X, y)\n# compare_preprocessing_df.to_csv(\"compare_preprocessing_df.csv\", index=False)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# proprocessed_auc_df = plot_preprocessed_roc_and_auc_df(X_train, y_train, X_test, y_test)\n# proprocessed_auc_df.to_csv(\"proprocessed_auc_df.csv\", index=False)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# select-than-interact","metadata":{}},{"cell_type":"code","source":"# 四種方法的特徵選擇\ndef select_features(X_train, y_train):\n    def get_feature_importance(model, X, method='coef'):\n        if method == 'coef':\n            importance = np.abs(model.coef_).flatten()\n        elif method == 'shap':\n            explainer = shap.Explainer(model, X)\n            importance = np.abs(explainer(X).values).mean(axis=0)\n        else:\n            importance = model.feature_importances_\n        return pd.DataFrame({'feature': X.columns, 'importance': importance}).sort_values(by='importance', ascending=False)\n    def generate_feature_groups(importance_df):\n        return {i: importance_df.head(i)['feature'].tolist() for i in range(2, len(importance_df) + 1)}\n    xgb_model = XGBClassifier(random_state=42)\n    xgb_model.fit(X_train, y_train)\n    xgb_importance_df = get_feature_importance(xgb_model, X_train, method='xgb')\n    feature_groups_xgb = generate_feature_groups(xgb_importance_df)\n    lr_model = LogisticRegression(max_iter=1000, random_state=42)\n    lr_model.fit(X_train, y_train)\n    lr_importance_df = get_feature_importance(lr_model, X_train, method='coef')\n    feature_groups_lr = generate_feature_groups(lr_importance_df)\n    shap_xgb_importance_df = get_feature_importance(xgb_model, X_train, method='shap')\n    feature_groups_shap_xgb = generate_feature_groups(shap_xgb_importance_df)\n    shap_lr_importance_df = get_feature_importance(lr_model, X_train, method='shap')\n    feature_groups_shap_lr = generate_feature_groups(shap_lr_importance_df)\n    return {\n        'xgb': feature_groups_xgb,\n        'lr': feature_groups_lr,\n        'shap_xgb': feature_groups_shap_xgb,\n        'shap_lr': feature_groups_shap_lr\n    }\n\n# 創建 only selection dataframe\ndef make_only_selection_df(select_result, X_train, X_test, y_train, y_test):\n    def train_and_evaluate(model, feature_subset):\n        X_train_selected = X_train[feature_subset]\n        X_test_selected = X_test[feature_subset]\n        model.fit(X_train_selected, y_train)\n        y_pred = model.predict_proba(X_test_selected)[:, 1]\n        return roc_auc_score(y_test, y_pred)\n    results_lr = {key: [] for key in select_result.keys()}\n    results_xgb = {key: [] for key in select_result.keys()}\n    feature_counts = [25, 50, 75, 100, 125, 150, 175, 200, 225, 250, 275, 300]\n    for num_features in feature_counts:\n        for key, feature_dict in select_result.items():\n            selected_features = feature_dict[num_features]\n            lr_model = LogisticRegression(max_iter=1000, random_state=42)\n            xgb_model = XGBClassifier(random_state=42)\n            results_lr[key].append(train_and_evaluate(lr_model, selected_features))\n            results_xgb[key].append(train_and_evaluate(xgb_model, selected_features))\n    df_results_lr = pd.DataFrame(results_lr, index=feature_counts)\n    df_results_xgb = pd.DataFrame(results_xgb, index=feature_counts)\n    return df_results_lr, df_results_xgb\n\n# 繪製 only selection auc\ndef plot_only_selection_auc(lr_df, xgb_df):\n    plt.figure(figsize=(18, 18))\n    x = lr_df.index\n    methods = {\n        \"xgb\": {\"color\": \"blue\", \"marker\": \"o\"},\n        \"lr\": {\"color\": \"red\", \"marker\": \"s\"},\n        \"shap_xgb\": {\"color\": \"green\", \"marker\": \"D\"},\n        \"shap_lr\": {\"color\": \"purple\", \"marker\": \"^\"},\n    }\n    for method, style in methods.items():\n        plt.plot(x, lr_df[method], label=f'lr: {method}_selection', color=style[\"color\"], marker=style[\"marker\"], linestyle=\"-\")\n        plt.plot(x, xgb_df[method], label=f'xgb: {method}_selection', color=style[\"color\"], marker=style[\"marker\"], linestyle=\"--\")\n    plt.plot(x, [0.9566] * len(x), label=\"lr: benchmark\", color=\"black\", linewidth=2.5, linestyle=\"-\")\n    plt.plot(x, [0.9554] * len(x), label=\"xgb: benchmark\", color=\"black\", linewidth=2.5, linestyle=\"--\")\n    plt.xlabel(\"Number of features\", fontsize=20)\n    plt.ylabel(\"AUC\", fontsize=20)\n    plt.ylim(0.953, 0.958)\n    plt.grid(True, which='both', linestyle='--', linewidth=0.5)\n    legend = plt.legend(fontsize=18, title=\"Raw features\", title_fontsize=20)\n    plt.xticks(fontsize=20)\n    plt.yticks(fontsize=20)\n    plt.savefig(f'only_selection_auc.png', transparent=True, dpi=300)\n    plt.show()\n\n# 生成多項式特徵\ndef generate_polynomial_dataset(df, features, degree):\n    interactions = pd.DataFrame(index=df.index)\n    for d in range(2, degree + 1):\n        for combination in itertools.combinations_with_replacement(features, d):\n            new_feature_name = \"_x_\".join(combination)\n            interactions[new_feature_name] = df[list(combination)].prod(axis=1)\n    return interactions\n\n# 創建 select than interact dataframe\ndef make_select_than_interact_df(model_type, degree, max_features, select_result, X_train, X_test, y_train, y_test):\n    results = {}\n    for method, feature_groups in select_result.items():\n        for k, features in itertools.islice(feature_groups.items(), max_features):\n            interaction_features_train = generate_polynomial_dataset(X_train, features, degree)\n            interaction_features_test = generate_polynomial_dataset(X_test, features, degree)\n            X_train_augmented = pd.concat([X_train, interaction_features_train], axis=1)\n            X_test_augmented = pd.concat([X_test, interaction_features_test], axis=1)\n            if model_type == 'xgb':\n                model = XGBClassifier(random_state=42)\n            elif model_type == 'lr':\n                model = LogisticRegression(max_iter=1000, random_state=42)\n            model.fit(X_train_augmented, y_train)\n            y_pred_prob = model.predict_proba(X_test_augmented)[:, 1]\n            roc_auc = roc_auc_score(y_test, y_pred_prob)\n            results[f\"{method}_{k}\"] = roc_auc\n            del interaction_features_train, interaction_features_test, X_train_augmented, X_test_augmented\n    data = {method: {int(k.split('_')[-1]): v for k, v in results.items() if k.startswith(method + \"_\")} \n            for method in select_result.keys()}\n    df = pd.DataFrame(data)\n    return df\n\n# 繪製 select than interact auc\ndef plot_select_than_interact_auc(lr_df, xgb_df, degree):\n    plt.figure(figsize=(18, 18))\n    base_x = lr_df.index\n    methods = {\n        \"xgb\": {\"color\": \"blue\", \"marker\": \"o\"},\n        \"lr\": {\"color\": \"red\", \"marker\": \"s\"},\n        \"shap_xgb\": {\"color\": \"green\", \"marker\": \"D\"},\n        \"shap_lr\": {\"color\": \"purple\", \"marker\": \"^\"},\n    }\n    for method, style in methods.items():\n        if degree == 2:    \n            x = 481 + base_x * (base_x + 1) / 2\n        elif degree == 3:\n            x = 481 + (base_x * (base_x + 1) / 2) + (base_x * (base_x + 1) * (base_x + 2) / 6)\n        plt.plot(x, lr_df[method], label=f'lr: {method}_selection', color=style[\"color\"], marker=style[\"marker\"], linestyle=\"-\")\n        plt.plot(x, xgb_df[method], label=f'xgb: {method}_selection', color=style[\"color\"], marker=style[\"marker\"], linestyle=\"--\")\n    plt.plot(x, [0.9566] * len(x), label=\"lr: benchmark\", color=\"black\", linewidth=2.5, linestyle=\"-\")\n    plt.plot(x, [0.9554] * len(x), label=\"xgb: benchmark\", color=\"black\", linewidth=2.5, linestyle=\"--\")\n    plt.xlabel(\"Number of features\", fontsize=20)\n    plt.ylabel(\"AUC\", fontsize=20)\n    plt.ylim(0.953, 0.958)\n    plt.grid(True, which='both', linestyle='--', linewidth=0.5)\n    legend = plt.legend(fontsize=20, title=\"Raw features\", title_fontsize=20)\n    plt.xticks(fontsize=20)\n    plt.yticks(fontsize=20)\n    if degree == 2:    \n        legend = plt.legend(fontsize=20, title=\"Quadratic augment features\", title_fontsize=20)\n        plt.savefig(f'quadratic_select_than_interact_auc.png', transparent=True, dpi=300)\n    elif degree == 3:\n        legend = plt.legend(fontsize=20, title=\"Cubic augment features\", title_fontsize=20)\n        plt.savefig(f'cubic_select_than_interact_auc.png', transparent=True, dpi=300)\n    plt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# select_result = select_features(X_train, y_train)\n# print(\"前 15 個特徵:\", select_result['shap_xgb'][15])\n# print(\"前 30 個特徵:\", select_result['shap_xgb'][30])","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# df_results_lr, df_results_xgb = make_only_selection_df(select_result, X_train, X_test, y_train, y_test)\n# df_results_lr.to_csv(\"lr_only_selection_df.csv\")\n# df_results_xgb.to_csv(\"xgb_only_selection_df.csv\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# plot_only_selection_auc(df_results_lr, df_results_xgb)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# lr_degree2_df = make_select_than_interact_df('lr', 2, 29, select_result, X_train, X_test, y_train, y_test)\n# lr_degree2_df.to_csv(\"lr_degree2_df.csv\")\n# xgb_degree2_df = make_select_than_interact_df('xgb', 2, 29, select_result, X_train, X_test, y_train, y_test)\n# xgb_degree2_df.to_csv(\"xgb_degree2_df.csv\")\n# lr_degree3_df = make_select_than_interact_df('lr', 3, 14, select_result, X_train, X_test, y_train, y_test)\n# lr_degree3_df.to_csv(\"lr_degree3_df.csv\")\n# xgb_degree3_df = make_select_than_interact_df('xgb', 3, 14, select_result, X_train, X_test, y_train, y_test)\n# xgb_degree3_df.to_csv(\"xgb_degree3_df.csv\")","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# plot_select_than_interact_auc(lr_degree2_df, xgb_degree2_df, 2)\n# plot_select_than_interact_auc(lr_degree3_df, xgb_degree3_df, 3)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### LASSO vs PowerSHAP","metadata":{}},{"cell_type":"code","source":"#製作 lasso 與 powershap 的 dataframe\ndef make_lasso_and_powershap_df(X_train, y_train, X_test, y_test, degree, features, lambda_values):\n    auc_result = []\n    if features:\n        interaction_train = generate_polynomial_dataset(X_train, features, degree)\n        interaction_test = generate_polynomial_dataset(X_test, features, degree)\n        X_train = pd.concat([X_train, interaction_train], axis=1)\n        X_test = pd.concat([X_test, interaction_test], axis=1)\n    for lambda_val in lambda_values:\n        model = LogisticRegression(penalty='l1', C=1/lambda_val, solver='liblinear', random_state=42)\n        model.fit(X_train, y_train)\n        selected = np.where(model.coef_[0] != 0)[0]\n        train_auc = roc_auc_score(y_train, model.predict_proba(X_train)[:, 1])\n        test_auc = roc_auc_score(y_test, model.predict_proba(X_test)[:, 1])\n        auc_result.append({'model': f'lasso (λ={lambda_val})', 'p': len(selected), 'train': train_auc, 'test': test_auc})\n    with contextlib.redirect_stdout(io.StringIO()):\n        powershap = PowerShap(model=XGBClassifier(random_state=42), power_iterations=10, alpha=0.05, power_req_iterations=1)\n        powershap.fit(X_train, y_train)\n    selected_cols = X_train.columns[powershap.get_support()]\n    X_train = X_train[selected_cols]\n    X_test = X_test[selected_cols]\n    model = LogisticRegression(random_state=42)\n    model.fit(X_train, y_train)\n    train_auc = roc_auc_score(y_train, model.predict_proba(X_train)[:, 1])\n    test_auc = roc_auc_score(y_test, model.predict_proba(X_test)[:, 1])\n    auc_result.append({'model': f'powershap (α=0.5)', 'p': len(selected_cols), 'train': train_auc, 'test': test_auc})\n    return pd.DataFrame(auc_result)\n    \n#繪製 lasso 和 powershap 的 auc\ndef plot_lasso_vs_powershap_auc(auc_df, degree):\n    lasso_df = auc_df[auc_df['model'].str.contains('lasso')].copy()\n    powershap_df = auc_df[auc_df['model'].str.contains('powershap')].copy()\n    lasso_df = lasso_df.sort_values(by='p')\n    plt.figure(figsize=(10, 6))\n    plt.plot(lasso_df['p'], lasso_df['train'], marker='o', linestyle='-', label='LR: lasso (train)', color='blue')\n    plt.plot(lasso_df['p'], lasso_df['test'], marker='o', linestyle='--', label='LR: lasso (test)', color='blue')\n    plt.scatter(powershap_df['p'], powershap_df['train'], color='red', marker='D', label='LR: powershap (train)')\n    plt.scatter(powershap_df['p'], powershap_df['test'], color='red', marker='X', label='LR: powershap (test)')\n    plt.xlabel('Number of features')\n    plt.ylabel('AUC')\n    # plt.ylim(0.72, 0.76)\n    plt.grid(True, which='both', linestyle='--', linewidth=0.5)\n    if degree == 1:    \n        legend = plt.legend(title=\"Raw features\")\n        plt.savefig(f'raw_lasso_vs_powershap_auc.png', transparent=True, dpi=300)\n    elif degree == 2:    \n        legend = plt.legend(title=\"Quadratic augment features\")\n        plt.savefig(f'quadratic_lasso_vs_powershap_auc.png', transparent=True, dpi=300)\n    elif degree == 3:\n        legend = plt.legend(title=\"Cubic augment features\")\n        plt.savefig(f'cubic_lasso_vs_powershap_auc.png', transparent=True, dpi=300)\n    plt.show()","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"raw_features = []\nraw_lambda = [100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1100, 1200, 1300, 1400, 1500, 1600, 1700, 1800, 1900, 2000, 2100, 2200, 2300, 2400, 2500, 2600, 2700, 2800, 2900, 3000]\nraw_df = make_lasso_and_powershap_df(X_train, y_train, X_test, y_test, 1, raw_features, raw_lambda)\nraw_df.to_csv(\"raw_lasso_vs_powershap_df.csv\")\nplot_lasso_vs_powershap_auc(raw_df, 1)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"quadratic_features = ['P_2', 'B_1', 'B_3', 'B_4', 'D_43', 'D_45', 'S_3', 'D_39', 'D_48', 'B_7', 'B_5', 'D_46', 'D_112', 'D_47', 'D_42', 'D_129', 'D_62', 'B_2', 'D_121', 'R_1', 'S_8', 'S_5', 'D_41', 'B_16', 'S_26', 'D_49', 'B_28', 'B_17', 'D_52', 'D_59']\nquadratic_lambda = [100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1100, 1200, 1300, 1400, 1500, 1600, 1700, 1800, 1900, 2000, 2100, 2200, 2300, 2400, 2500, 2600, 2700, 2800, 2900, 3000]\nquadratic_df = make_lasso_and_powershap_df(X_train, y_train, X_test, y_test, 2, quadratic_features, quadratic_lambda)\nquadratic_df.to_csv(\"quadratic_lasso_vs_powershap_df.csv\")\nplot_lasso_vs_powershap_auc(quadratic_df, 2)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"cubic_features = ['P_2', 'B_1', 'B_3', 'B_4', 'D_43', 'D_45', 'S_3', 'D_39', 'D_48', 'B_7', 'B_5', 'D_46', 'D_112', 'D_47', 'D_42']\ncubic_lambda = [100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 1100, 1200, 1300, 1400, 1500, 1600, 1700, 1800, 1900, 2000, 2100, 2200, 2300, 2400, 2500, 2600, 2700, 2800, 2900, 3000]\ncubic_df = make_lasso_and_powershap_df(X_train, y_train, X_test, y_test, 3, cubic_features, cubic_lambda)\ncubic_df.to_csv(\"cubic_lasso_vs_powershap_df.csv\")\nplot_lasso_vs_powershap_auc(cubic_df, 3)","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}