重回帰と制約条件
ビジネスの要請
多くのビジネスにおいて「この数字は何をいじれば変わるのか」に対する回答を迫られる。この時、因果の検証を目的とした重回帰や、変数への考察から因子分析を行うことが考えられる。しかし、多くの陥穽に気付かず、単純なプログラムのコピペによる結果を解釈することほど哀しいものもない。
ここで予測ではなく、解釈可能性＝特徴量重要度の観点から幾つか検証と理解の整理を試みた
重回帰と信憑性
予測ではなく、解釈可能性＝特徴量重要度の観点から幾つか検証と理解の整理を試みた。ボストン市の住宅価格データを使用する。
# 必須モジュール import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns from sklearn.datasets import load_boston from sklearn.model_selection import train_test_split # データ読み込み boston = load_boston() df = pd.DataFrame(boston.data, columns=boston.feature_names) df['PRICE'] = boston.target # 説明変数、目的変数を定義 X, y = boston.data, boston.target # データ分割 X_train, X_test, y_train, y_test = train_test_split(X, y, train_size =0.8, random_state=0)
重回帰では、各種特徴量がみられる statsmodels のほうがscikit-learnより好み。定数項を入れるのはやや面倒。
import statsmodels.api as sm # 回帰モデル作成 mod = sm.OLS(y_train, sm.add_constant(X_train)) # 訓練 result = mod.fit() # 統計サマリを表示 print(result.summary()) OLS Regression Results ============================================================================== Dep. Variable: PRICE R-squared: 0.773 Model: OLS Adj. R-squared: 0.765 Method: Least Squares F-statistic: 102.2 Date: Sun, 20 Jun 2021 Prob (F-statistic): 9.64e-117 Time: 21:52:40 Log-Likelihood: -1171.5 No. Observations: 404 AIC: 2371. Df Residuals: 390 BIC: 2427. Df Model: 13 Covariance Type: nonrobust ============================================================================== coef std err t P>|t| [0.025 0.975] ------------------------------------------------------------------------------ const 38.0917 5.522 6.898 0.000 27.234 48.949 CRIM -0.1194 0.037 -3.257 0.001 -0.192 -0.047 ZN 0.0448 0.014 3.102 0.002 0.016 0.073 INDUS 0.0055 0.063 0.087 0.931 -0.119 0.130 CHAS 2.3408 0.902 2.595 0.010 0.567 4.115 NOX -16.1236 4.212 -3.828 0.000 -24.404 -7.843 RM 3.7087 0.458 8.106 0.000 2.809 4.608 AGE -0.0031 0.014 -0.218 0.828 -0.031 0.025 DIS -1.3864 0.214 -6.480 0.000 -1.807 -0.966 RAD 0.2442 0.070 3.481 0.001 0.106 0.382 TAX -0.0110 0.004 -2.819 0.005 -0.019 -0.003 PTRATIO -1.0459 0.137 -7.636 0.000 -1.315 -0.777 B 0.0081 0.003 2.749 0.006 0.002 0.014 LSTAT -0.4928 0.054 -9.086 0.000 -0.599 -0.386 ============================================================================== Omnibus: 141.494 Durbin-Watson: 1.996 Prob(Omnibus): 0.000 Jarque-Bera (JB): 629.882 Skew: 1.470 Prob(JB): 1.67e-137 Kurtosis: 8.365 Cond. No. 1.55e+04 ============================================================================== Warnings: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified. [2] The condition number is large, 1.55e+04. This might indicate that there are strong multicollinearity or other numerical problems.
この結果から、決定係数、係数、t値、p値あたりをパッとみつつ説明変数の重要度をみることができる。そしてVIF計算して多重共線性を疑ったり、交絡因子を考えたり、ステップワイズ法なり主成分分析から変数減らしたり、操作変数法云々と、兎に角こねくり回すのが重回帰におけるプロセスになる。
だがビジネス側の要請として、変数を落としたくない！でも係数はマイナスになる筈がない！特徴量重要度として列挙するんだ！という話があるのではないだろうか（あった）
その意味の在処はさておき、そんな話が振ってきたらどうすればよいのか。結論として、制約条件をかけた最適化処理が必要となるであろう。ここの備忘録を記す。
係数のマイナスを許さない法
重回帰における最小二乗法は行列式で比較的容易に表現できる。そこで係数がマイナスなのが嫌ならば、マイナスにしなければ良いのだ。決定係数の行方は如何に。
# 普通の重回帰 class LinearRegression1: def __init__(self): self.w_ = None def fit(self, X, y): # 切片の計算のため先頭列に１を追加 X = np.insert(X, 0, 1, axis=1) self.w_ = np.linalg.inv(X.T @ X) @ X.T @ y def predict(self, X): X = np.insert(X, 0, 1, axis=1) return X @ self.w_ def score(self, X, t): #残差の分散を求める u = ((t - self.predict(X))**2).sum() #実測値の分散を求める v = ((t - t.mean())**2).sum() #決定係数を返す return 1 - u/v # 結果出力 clf1 = LinearRegression1() clf1.fit(X_train, y_train) print(clf1.score(X_test,y_test)) # 横軸がy(実際値)、縦軸がy-y_pred（誤差） y_pred1 = clf1.predict(X) plt.figure(figsize=(6, 6)) plt.hlines([0], min(y)-5, max(y) , color='blue') plt.scatter(y, y - y_pred1, color='red', marker='o', alpha=0.2) plt.xlabel('y') plt.ylabel('y - y_pred1'); # 係数確認 coef_df = pd.DataFrame(clf1.w_[1:], index=boston.feature_names, columns=['clf1']) coef_df
Scipyは、Non-negative least squaresを有している。sklearnで作っておく。
# NLS 係数非負 重回帰 from scipy.optimize import nnls class LinearRegression2: def __init__(self): self.w_ = None def fit(self, X, y): # 先頭列に１を追加 # 追加の理由は先ほどの説明を参照（切片の計算のため） X = np.insert(X, 0, 1, axis=1) #self.w_ = np.linalg.inv(X.T @ X) @ X.T @ y self.w_ = nnls(X,y)[0] def predict(self, X): X = np.insert(X, 0, 1, axis=1) return X @ self.w_ def score(self, X, t): #残差の分散を求める u = ((t - self.predict(X))**2).sum() #実測値の分散を求める v = ((t - t.mean())**2).sum() #決定係数を返す return 1 - u/v clf2 = LinearRegression2() clf2.fit(X_train, y_train) print(clf1.score(X_test,y_test)) coef_df['clf2'] = clf2.w_[1:] coef_df.round(2) # 横軸がy(実際値)、縦軸がy-y_pred（誤差） y_pred2= clf2.predict(X) plt.figure(figsize=(6, 6)) plt.hlines([0], min(y)-5, max(y) , color='blue') plt.scatter(y, y - y_pred2, color='red', marker='o', alpha=0.2) plt.xlabel('y') plt.ylabel('y - y_pred2');
Scipyは、完全な０も困る場合は、lsq_linearで。
from scipy.sparse import rand from scipy.optimize import lsq_linear class LinearRegression3: def __init__(self): self.w_ = None def fit(self, X, y): # 先頭列に１を追加 # 追加の理由は先ほどの説明を参照（切片の計算のため） X = np.insert(X, 0, 1, axis=1) #self.w_ = np.linalg.inv(X.T @ X) @ X.T @ y # 0.01 ~ 100 にする場合はboundsで指定 res = lsq_linear(X, y, bounds=(0.01, 100), lsmr_tol='auto', verbose=1) self.w_ =res.x def predict(self, X): X = np.insert(X, 0, 1, axis=1) return X @ self.w_
ちな決定係数は準に、0.76, 0.41, 0.17という悲惨な結果。しかし良いのだ。。問題は係数を負にしないことだけなのだから、、、
ところでランダムフォレストでお茶は濁せないのか？
ランダムフォレストつまりは決定木における特徴量重要度（Gain）は、多重共線性の影響を受けないわけでは勿論ない。だが、係数として顕現しないので説明時にはごまかせるのではないか…？
# 自由度調整済みr2を算出 def adjusted_r2(X,Y,model): from sklearn.metrics import r2_score import numpy as np r_squared = r2_score(Y, model.predict(X)) adjusted_r2 = 1 - (1-r_squared)*(len(Y)-1)/(len(Y)-X.shape[1]-1) #yhat = model.predict(X) \ #SS_Residual = sum((Y-yhat)**2) \ #SS_Total = sum((Y-np.mean(Y))**2) #r_squared = 1 - (float(SS_Residual))/ SS_Total return adjusted_r2 # 予測モデルの精度確認の各種指標を算出 def get_model_evaluations(X_train,Y_train,X_test,Y_test,model): from sklearn.metrics import explained_variance_score from sklearn.metrics import mean_absolute_error from sklearn.metrics import mean_squared_error from sklearn.metrics import mean_squared_log_error from sklearn.metrics import median_absolute_error # 評価指標確認 # 参考: https://funatsu-lab.github.io/open-course-ware/basic-theory/accuracy-index/ yhat_test = model.predict(X_test) return "adjusted_r2(train) :" + str(adjusted_r2(X_train,Y_train,model)) \ , "adjusted_r2(test) :" + str(adjusted_r2(X_test,Y_test,model)) \ , "平均誤差率(test) :" + str(np.mean(abs(Y_test / yhat_test - 1))) \ , "MAE(test) :" + str(mean_absolute_error(Y_test, yhat_test)) \ , "MedianAE(test) :" + str(median_absolute_error(Y_test, yhat_test)) \ , "RMSE(test) :" + str(np.sqrt(mean_squared_error(Y_test, yhat_test))) \ , "RMSE(test) / MAE(test) :" + str(np.sqrt(mean_squared_error(Y_test, yhat_test)) / mean_absolute_error(Y_test, yhat_test)) #better if result = 1.253 # 必須モジュール import xgboost as xgb import matplotlib.pylab as plt from sklearn.ensemble import RandomForestRegressor # 説明変数、目的変数を定義 X = df.drop('PRICE', axis=1) y = df['PRICE'] # データ分割 X_train, X_test, y_train, y_test = train_test_split(X, y, train_size =0.8, random_state=1) # sklearnによるランダムフォレスト rf_reg = RandomForestRegressor(n_estimators=10) rf_reg = rf_reg.fit(X_train, y_train) fti = rf_reg.feature_importances_ print('Feature Importances:') for i, feat in enumerate(boston['feature_names']): print('\t{0:10s} : {1:>.6f}'.format(feat, fti[i])) # 結果評価 get_model_evaluations(X_train,y_train,X_test,y_test,rf_reg) ###################################### # xgboost による回帰① # Dmatrix作りつつ、説明変数は入れておく dtrain = xgb.DMatrix(X_train, label=y_train,feature_names=boston['feature_names']) dtest = xgb.DMatrix(X_test, label=y_test,feature_names=boston['feature_names']) # 回帰パラメタ設定 xgb_params = {"objective": "reg:linear", "eta": 0.1, "max_depth": 6, "silent": 1} num_rounds = 100 # 学習 gbdt = xgb.train(xgb_params, dtrain, num_rounds) # 特徴量重要度を、weight と gainで出したい fig = plt.figure(figsize=(13, 5)) ax1 = plt.subplot(121) ax1 = xgb.plot_importance(gbdt, ax=ax1,importance_type='weight') ax2 = plt.subplot(122) ax2 = xgb.plot_importance(gbdt, ax=ax2,importance_type='gain') ###################################### # xgboost による回帰② model_A = xgb.XGBRegressor(feature_names=boston['feature_names']) model_A.fit(X_train, y_train) get_model_evaluations(X_train,y_train,X_test,y_test,model_A) y_train_pred = model_A.predict(X_train) y_test_pred = model_A.predict(X_test) from sklearn.metrics import r2_score from sklearn.metrics import mean_squared_error print('MSE train: %.3f, test: %.3f' % ( mean_squared_error(y_train, y_train_pred), mean_squared_error(y_test, y_test_pred))) print('RMSE train: %.3f, test: %.3f' %( np.sqrt(mean_squared_error(y_train, y_train_pred)), np.sqrt(mean_squared_error(y_test, y_test_pred)))) print('R^2 train: %.3f, test: %.3f' % ( r2_score(y_train, y_train_pred), r2_score(y_test, y_test_pred))) # feature importance のプロット import pandas as pd import matplotlib.pyplot as plt importances = pd.Series(model_A.feature_importances_, index = boston.feature_names) importances = importances.sort_values() importances.plot(kind = 'barh') plt.title('XGBRegressorの重要度') plt.show() fig = plt.figure(figsize=(13, 5)) ax1 = plt.subplot(121) ax1 = xgb.plot_importance(model_A, ax=ax1,importance_type='weight') ax2 = plt.subplot(122) ax2 = xgb.plot_importance(model_A, ax=ax2,importance_type='gain') import shap shap.initjs() explainerA = shap.TreeExplainer(model = model_A) shap_valuesA = explainerA.shap_values(X = X_train)
とくに、単調整制約という概念を学んだ。
決定木で別れるノードの重みを、親兄弟の平均を境界として決める事で木全体の単調整を担保できる模様。Xgboostだと、parameterで指定可能というのもそそる仕様だ。
# 単調性制約を課した XGBoost で回帰モデルを構築します. print(X_train.shape) model_B = xgb.XGBRegressor(monotone_constraints=(1,1,1,1,1,1,1,1,1,1,1,1,1)) #model_B = xgb.XGBRegressor(monotone_constraints=(0,0,0,0,0,0,0,0,0,0,0,0,0)) model_B.fit(X_train, y_train) get_model_evaluations(X_train,y_train,X_test,y_test,model_B) shap.initjs() explainerB = shap.TreeExplainer(model = model_B) shap_valuesB = explainerB.shap_values(X = X_train) shap.summary_plot(shap_valuesB, X_train, plot_type = "bar")
結果に対する解釈については、もう少し時間を取ってゆっくり考えることとしたい。一旦。
